2. Data management, model description and testing

Data management

mlogit is loaded using:

library("mlogit")

It comes with several data sets that we’ll use to illustrate the features of the library. Data sets used for multinomial logit estimation concern some individuals, that make one or a sequential choice of one alternative among a set of mutually exclusive alternatives. The determinants of these choices are covariates that can depend on the alternative and the choice situation, only on the alternative or only on the choice situation.

To illustrate this typology of the covariates, consider the case of repeated choice of destinations for vacations by families:

• the length of the vacation, the season are choice situation specific variables,
• income, family size are individual specific variables,
• distance to destination, cost are alternative specific variables.

The unit of observation is therefore the choice situation, and it is also the individual if there is only one choice situation per individual observed, which is often the case.

Such data have therefore a specific structure that can be characterized by three indexes: the alternative, the choice situation and the individual. These three indexes will be denoted alt, chid and id. Note that the distinction between chid and id is only relevant if we have repeated observations for the same individual.

Data sets can have two different shapes: a wide shape (one row for each choice situation) or a long shape (one row for each alternative and, therefore, as many rows as there are alternatives for each choice situation).

mlogit deals with both format. It depends on the dfidx package which takes as first argument a data.frame and returns a dfidx object, which is a data.frame in “long” format with a special data frame column which contains the indexes.

data("Train", package = "mlogit")
Train$choiceid <- 1:nrow(Train)
head(Train, 3)

##   id choiceid choice price_A time_A change_A comfort_A price_B
## 1  1        1      A    2400    150        0         1    4000
## 2  1        2      A    2400    150        0         1    3200
## 3  1        3      A    2400    115        0         1    4000
##   time_B change_B comfort_B
## 1    150        0         1
## 2    130        0         1
## 3    115        0         0

Tr <- dfidx(Train, shape = "wide", varying = 4:11, sep = "_",
            idx = list(c("choiceid", "id")), idnames = c(NA, "alt"))

Tr$price <- Tr$price / 100 * 2.20371
Tr$time <- Tr$time / 60

head(Tr, 3)

## ~~~~~~~
##  first 3 observations out of 5858 
## ~~~~~~~
##   choice    price time change comfort idx
## 1      A 52.88904  2.5      0       1 1:A
## 2      A 88.14840  2.5      0       1 1:B
## 3      A 52.88904  2.5      0       1 2:A
## 
## ~~~ indexes ~~~~
##   choiceid id alt
## 1        1  1   A
## 2        1  1   B
## 3        2  1   A
## indexes:  1, 1, 2

head(idx(Tr), 3)

##   choiceid id alt
## 1        1  1   A
## 2        1  1   B
## 3        2  1   A
## indexes:  1, 1, 2

data("ModeCanada", package = "mlogit")
head(ModeCanada)

##   case   alt choice dist  cost ivt ovt freq income urban noalt
## 1    1 train      0   83 28.25  50  66    4     45     0     2
## 2    1   car      1   83 15.77  61   0    0     45     0     2
## 3    2 train      0   83 28.25  50  66    4     25     0     2
## 4    2   car      1   83 15.77  61   0    0     25     0     2
## 5    3 train      0   83 28.25  50  66    4     70     0     2
## 6    3   car      1   83 15.77  61   0    0     70     0     2

MC <- dfidx(ModeCanada, subset = noalt == 4,
            alt.levels = c("train", "air", "bus", "car"))

MC <- dfidx(ModeCanada, subset = noalt == 4, idx = list(NA, "alt"))

MC <- dfidx(ModeCanada, subset = noalt == 4, idx = "case",
            alt.levels = c("train", "air", "bus", "car"))

MC <- dfidx(ModeCanada, subset = noalt == 4, idx = c("case", "alt"))

MC <- dfidx(ModeCanada, subset = noalt == 4)

MC <- dfidx(ModeCanada, subset = noalt == 4, idx = c("case", "alt"),
            drop.index = FALSE)
head(MC)

## ~~~~~~~
##  first 10 observations out of 11116 
## ~~~~~~~
##    case   alt choice dist   cost ivt ovt freq income urban noalt
## 1   109 train      0  377  58.25 215  74    4     45     0     4
## 2   109   air      1  377 142.80  56  85    9     45     0     4
## 3   109   bus      0  377  27.52 301  63    8     45     0     4
## 4   109   car      0  377  71.63 262   0    0     45     0     4
## 5   110 train      0  377  58.25 215  74    4     70     0     4
## 6   110   air      1  377 142.80  56  85    9     70     0     4
## 7   110   bus      0  377  27.52 301  63    8     70     0     4
## 8   110   car      0  377  71.63 262   0    0     70     0     4
## 9   111 train      0  377  58.25 215  74    4     35     0     4
## 10  111   air      1  377 142.80  56  85    9     35     0     4
##         idx
## 1  109:rain
## 2   109:air
## 3   109:bus
## 4   109:car
## 5  110:rain
## 6   110:air
## 7   110:bus
## 8   110:car
## 9  111:rain
## 10  111:air
## 
## ~~~ indexes ~~~~
##    case   alt
## 1   109 train
## 2   109   air
## 3   109   bus
## 4   109   car
## 5   110 train
## 6   110   air
## 7   110   bus
## 8   110   car
## 9   111 train
## 10  111   air
## indexes:  1, 2

Model description

Standard formulas are not very practical to describe random utility models, as these models may use different sets of covariates. Actually, working with random utility models, one has to consider at most four sets of covariates:

1. alternative and choice situation specific covariates \(x_{ij}\) with generic coefficients \(\beta\) and and alternative specific covariates \(t_j\) with a generic coefficient \(\nu\),
2. choice situation specific covariates \(z_i\) with alternative specific coefficients \(\gamma_j\),
3. alternative and choice situation specific covariates \(w_{ij}\) with alternative specific coefficients \(\delta_j\),
4. choice situation specific covariates \(v_i\) that influence the variance of the errors.

The first three sets of covariates enter the observable part of the utility which can be written, alternative \(j\):

\[ V_{ij}=\alpha_j + \beta x_{ij} + \nu t_j + \gamma_j z_i + \delta_j w_{ij} . \]

As the absolute value of utility is irrelevant, only utility differences are useful to modelise the choice for one alternative. For two alternatives \(j\) and \(k\), we obtain:

\[ V_{ij}-V_{ik}=(\alpha_j-\alpha_k) + \beta (x_{ij}-x_{ik}) + (\gamma_j-\gamma_k) z_i + (\delta_j w_{ij} - \delta_k w_{ik}) + \nu(t_j - t_k). \]

It is clear from the previous expression that coefficients of choice situation specific variables (the intercept being one of those) should be alternative specific, otherwise they would disappear in the differentiation. Moreover, only differences of these coefficients are relevant and can be identified. For example, with three alternatives 1, 2 and 3, the three coefficients \(\gamma_1, \gamma_2, \gamma_3\) associated to a choice situation specific variable cannot be identified, but only two linear combinations of them. Therefore, one has to make a choice of normalization and the simplest one is just to set \(\gamma_1 = 0\).

Coefficients for alternative and choice situation specific variables may (or may not) be alternative specific. For example, transport time is alternative specific, but 10 mn in public transport may not have the same impact on utility than 10 mn in a car. In this case, alternative specific coefficients are relevant. Monetary cost is also alternative specific, but in this case, one can consider than 1$ is 1$ whatever it is spent for the use of a car or in public transports. In this case, a generic coefficient is relevant.

The treatment of alternative specific variables don’t differ much from the alternative and choice situation specific variables with a generic coefficient. However, if some of these variables are introduced, the \(\nu\) parameter can only be estimated in a model without intercepts to avoid perfect multicolinearity.

Individual-related heteroscedasticity (see Swait and Louviere 1993) can be addressed by writing the utility of choosing \(j\) for individual \(i\): \(U_{ij}=V_{ij} + \sigma_i \epsilon_{ij}\), where \(\epsilon\) has a variance that doesn’t depend on \(i\) and \(j\) and \(\sigma_i^2 = f(v_i)\) is a parametric function of some individual-specific covariates. Note that this specification induce choice situation heteroscedasticity, also denoted scale heterogeneity.³. As the overall scale of utility is irrelevant, the utility can also be writen as: \(U_{ij}^* = U_{ij} / \sigma_i = V_{ij}/\sigma_i + \epsilon_{ij}\), i.e., with homoscedastic errors. if \(V_{ij}\) is a linear combination of covariates, the associated coefficients are then divided by \(\sigma_i\).

A logit model with only choice situation specific variables is sometimes called a multinomial logit model, one with only alternative specific variables a conditional logit model and one with both kind of variables a mixed logit model. This is seriously misleading: conditional logit model is also a logit model for longitudinal data in the statistical literature and mixed logit is one of the names of a logit model with random parameters. Therefore, in what follows, we’ll use the name multinomial logit model for the model we’ve just described whatever the nature of the explanatory variables used.

mlogit package provides objects of class mFormula which are built upon Formula objects provided by the Formula package.⁴ The Formula package provides richer formulas, which accept multiple responses (a feature not used here) and multiple set of covariates. It has in particular specific model.frame and model.matrix methods which can be used with one or several sets of covariates.

To illustrate the use of mFormula objects, we use again the ModeCanada data set and consider three sets of covariates that will be indicated in a three-part formula, which refers to the first three items of the four points list at start of this section.

• cost (monetary cost) is an alternative specific covariate with a generic coefficient (part 1),
• income and urban are choice situation specific covariates (part 2),
• ivt (in vehicle travel time) is alternative specific and alternative specific coefficients are expected (part 3).

library("Formula")
f <- Formula(choice ~ cost | income + urban | ivt)

Some parts of the formula may be omitted when there is no ambiguity. For example, the following sets of formulas are identical:

f2 <- Formula(choice ~ cost + ivt | income + urban)
f2 <- Formula(choice ~ cost + ivt | income + urban | 0)

f3 <- Formula(choice ~ 0 | income | 0)
f3 <- Formula(choice ~ 0 | income)

f4 <- Formula(choice ~ cost + ivt)
f4 <- Formula(choice ~ cost + ivt | 1)
f4 <- Formula(choice ~ cost + ivt | 1 | 0)

By default, an intercept is added to the model, it can be removed by using + 0 or - 1 in the second part.

f5 <- Formula(choice ~ cost | income + 0 | ivt)
f5 <- Formula(choice ~ cost | income - 1 | ivt)

model.frame and model.matrix methods are provided for mFormula objects. The latter is of particular interest, as illustrated in the following example:

f <- Formula(choice ~ cost | income  | ivt)
mf <- model.frame(MC, f)
head(model.matrix(mf), 4)

##   (Intercept)   cost
## 1           1  58.25
## 2           1 142.80
## 3           1  27.52
## 4           1  71.63

The model matrix contains \(J-1\) columns for every choice situation specific variables (income and the intercept), which means that the coefficient associated to the first alternative (air) is set to 0. It contains only one column for cost because we want a generic coefficient for this variable. It contains \(J\) columns for ivt, because it is an alternative specific variable for which we want alternative specific coefficients.

Testing

As for all models estimated by maximum likelihood, three testing procedures may be applied to test hypothesis about models fitted using mlogit. The set of hypothesis tested defines two models: the unconstrained model that doesn’t take these hypothesis into account and the constrained model that impose these hypothesis.

This in turns define three principles of tests: the Wald test, based only on the unconstrained model, the Lagrange multiplier test (or score test), based only on the constrained model and the likelihood ratio test, based on the comparison of both models.

Two of these three tests are implemented in the lmtest package (Zeileis and Hothorn 2002): waldtest and lrtest. The Wald test is also implemented as linearHypothesis in package car (???), with a fairly different syntax. We provide special methods of waldtest and lrtest for mlogit objects and we also provide a function for the Lagrange multiplier (or score) test called scoretest.

We’ll see later that the score test is especially useful for mlogit objects when one is interested in extending the basic multinomial logit model because, in this case, the unconstrained model may be difficult to estimate. For the presentation of further tests, we provide a convenient statpval function which extract the statistic and the p-value from the objects returned by the testing function, which can be either of class anova or htest.

statpval <- function(x){
    if (inherits(x, "anova")) 
        result <- as.matrix(x)[2, c("Chisq", "Pr(>Chisq)")]
    if (inherits(x, "htest")) result <- c(x$statistic, x$p.value)
    names(result) <- c("stat", "p-value")
    round(result, 3)
}