# mlogit package fails to recover synthetic mixed logit model

I am generating data from the following synthetic mixed effects model for the utility of agent $h$ choosing transportation mode $j$:

$U_{hj} = \beta_{\text{price}} X_{h,\text{price}_j} + \beta_{\text{time}} X_{h,\text{time}_j} + \gamma_h X_{h,\text{bus}_j}+ \epsilon_{hj} \\ \gamma_h \sim \text{N}(0, 10)$

where $X_{h,\text{price}_j},X_{h,\text{time}_j}$ are agent $h$'s price and time for transporation mode $j$. $X_{h,\text{bus}_j}$ is an indicator function which is equal to 1 if mode $j$ is a bus.

The true model has $\beta_{\text{price}} = \beta_{\text{time}} = -1$.

For each agent, I generate the price and time covariates ($X_{h,\text{price}_j},X_{h,\text{time}_j}$) for each mode of transportation (red bus, train, car) and I sample the random effect $\gamma_h$ for each agent from $N(0,10)$, and then sample the transportation choice from the corresponding logit model.

Here's the R code that generates the choice data:

library(mlogit)
set.seed(1987)
n.samples <- 10000
re.sd     <- 10 # random effect standard deviation

## generate price, time covariates
mode       <- c('red.bus', 'train', 'car')
is.bus     <- c(1, 0, 0)
price.data <- rnorm(3*n.samples) + c(2,5,10)
time.data  <- rnorm(3*n.samples) + c(10, 5, 2)
obs          <- data.frame(mode=mode,
price=price.data,
time=time.data,
bus=is.bus)

## generate random effect + choice data
obs$agent.id <- rep(1:n.samples, each=3) obs$choice   <- NA
obs$util <- NA for( a.id in unique(obs$agent.id) ) {
xx.sub <- obs[obs$agent.id == a.id,] obs$util[obs$agent.id == a.id] <- xx.sub$price*-1 + xx.sub$time*-1 + xx.sub$bus*rnorm(1,sd=re.sd)
uu <- obs$util[obs$agent.id == a.id]
p.vec <- exp(uu)/sum(exp(uu))
obs$choice[obs$agent.id == a.id] <- rmultinom(1, 1, p.vec)==1
}

logit.data <- mlogit.data(obs,
shape   = "long",
choice  = "choice",
varying = which( colnames(obs) %in% c('price', 'time', 'bus') ),
alt.var = 'mode')


Here are the first few rows of the generated dataset:

> head(logit.data)
mode      price     time bus agent.id choice       util
1.red.bus red.bus  0.7113747 9.700200   1        1  FALSE -21.306406
1.train     train  5.2730015 5.908244   0        1   TRUE -11.181246
1.car         car 12.8485015 2.256558   0        1  FALSE -15.105060
2.red.bus red.bus  1.7021472 9.642485   1        2   TRUE  -7.677141
2.train     train  5.2443204 5.671528   0        2  FALSE -10.915848
2.car         car 10.2686018 2.250377   0        2  FALSE -12.518979


I attempt to fit the correctly specified model with the mlogit package, but I find that the estimates are wrong (the standard deviation of the random effect, in particular):

m.mixed <- mlogit(choice ~ price + time + bus | 0,
data=logit.data,
rpar= c(bus = 'n'),
R = 300, halton = NA)

summary(m.mixed)
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data,
rpar = c(bus = "n"), R = 300, halton = NA)

Frequencies of alternatives:
car red.bus   train
0.1317  0.4084  0.4599

bfgs method
8 iterations, 0h:1m:55s
g'(-H)^-1g =   594
last step couldn't find higher value

Coefficients :
Estimate Std. Error  t-value  Pr(>|t|)
price  -1.253720   0.026698 -46.9592 < 2.2e-16 ***
time   -1.329888   0.032389 -41.0603 < 2.2e-16 ***
bus     1.402522   0.318802   4.3993 1.086e-05 ***
sd.bus 19.455180   2.107521   9.2313 < 2.2e-16 ***


The true sd.bus is 10, but mlogit estimates it at almost 20.

Why isn't mlogit recovering the true model?

• Interesting question. My first thought was that it might just be due to a single data point, so I re-ran your simulation 20 times (using set.seed(i) for the ì-th run), and collecting the sd.bus estimates. The minimum is 9.0, the mean 18.0, the maximum 33.0. Something indeed seems to be off here. Aug 2, 2018 at 19:59
• I think rep(1:n.samples, each=length(price)) should be rep(1:n.samples, each = 3). Aug 2, 2018 at 21:46
• Thanks @BenjaminChristoffersen - sloppy editing on my part. And just to be clear to other readers: the problem I wrote about remains. Aug 3, 2018 at 1:20

You appear to have hit upon an unlucky combination of optimization parameters, specifically, with respect to the Halton pseudo-random sequence, which may possibly have a bug in it. BFGS appears to be stopping prematurely with R = 300, but not with other significantly smaller or larger values. Fortunately, you don't need large values of R (or Halton at all) in this case.

On my initial run through your code, I got the same results as you did, with runtime statistics indicating that 8 iterations were required for BFGS to converge. I then changed R in the function call to equal 30:

> m.mixed <- mlogit(choice ~ price + time + bus | 0,
+                   data=logit.data,
+                   rpar= c(bus = 'n'),
+                   R = 30, halton = NA)
>
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data,
rpar = c(bus = "n"), R = 30, halton = NA)

Frequencies of alternatives:
car red.bus   train
0.1317  0.4084  0.4599

bfgs method
22 iterations, 0h:1m:28s
g'(-H)^-1g = 3.83E-07

Coefficients :
Estimate Std. Error  z-value Pr(>|z|)
price  -0.988457   0.026861 -36.7989   <2e-16 ***
time   -0.990255   0.032661 -30.3195   <2e-16 ***
bus    -0.118121   0.227826  -0.5185   0.6041
sd.bus 10.369252   0.846377  12.2513   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -8743.9


Note that not only are the coefficient estimates correct, in the sense of being reasonably close to the actual values given their standard errors, the runtime is 40% that of the R = 300 case, although requiring 22 BFGS iterations instead of 8.

With R = 100 22 iterations were also needed, the runtime increased by a little over 30% to about half required in the initial case, and the coefficient estimates were essentially the same as in the R = 30 run:

Coefficients :
Estimate Std. Error  z-value Pr(>|z|)
price  -0.989468   0.026877 -36.8144   <2e-16 ***
time   -0.991518   0.032679 -30.3410   <2e-16 ***
bus    -0.108744   0.226858  -0.4793   0.6317
sd.bus 10.316022   0.841978  12.2521   <2e-16 ***


Avoiding the default Halton sequence altogether also fixed the problem:

> m.mixed <- mlogit(choice ~ price + time + bus | 0,
+                   data=logit.data,
+                   rpar= c(bus = 'n'),
+                   R = 300, halton = NULL)
>
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data,
rpar = c(bus = "n"), R = 300, halton = NULL)

*** blah blah blah ***

Coefficients :
Estimate Std. Error  z-value Pr(>|z|)
price  -0.988426   0.026859 -36.8006   <2e-16 ***
time   -0.990478   0.032655 -30.3314   <2e-16 ***
bus    -0.143370   0.229167  -0.6256   0.5316
sd.bus 10.382572   0.847424  12.2519   <2e-16 ***


The runtime in this case was still 10% lower than in the Halton case, and BFGS required 23 iterations to converge. That 8 iterations in the R = 300 case is definitely an outlier.

Setting R = 400 and halton = NA also generated "correct" results, however, R = 299 broke the estimation again, BFGS requiring 12 iterations and the estimate of sd.bus = 15.248....

EDIT:

I also tried a different seed with 4000 samples, but R = 300 and halton = NA still generated bad results, even worse ones than in the original case as it happened. Reparameterizing the call to specify the Halton prime and drop parameters gave erratic results; prime=11 worked, but prime=29 failed miserably with R=300. I then went through the mlogit R code (thanks for finding it, @khoda!), but the Halton sequence code works correctly.

Multiple other tests, combined with the OP's tests noted in comments, leads me to conclude that Halton doesn't work consistently well in the one-dimensional case, at least for this problem. In actual practice, where the true parameters are unavailable for comparison with the estimates, it would be necessary to try several different parameterizations of halton and R in the mlogit call, and check for consistency of the results (and the value of the log likelihood, I suspect). Avoiding Halton altogether and specifying an increasing sequence of R values for the random number generator until stable estimates are achieved is an alternative that would also likely be workable, runtime considerations aside.

• This is great, thank you! I believe the Halton-specific code lives in the function make.random.nb() in the file mlogit/R/mlogit.rpar.R. Let me know if you spot an error! Given the popularity of this package, I am quite surprised by these results. I suppose that on a real/applied problem, this suggests that we should turn Halton sampling off? Aug 6, 2018 at 16:41
• Either that or don't use a large R! When you're using structured pseudo-random numbers, you don't actually need that many for effective integration, in my experience (others may disagree.) For example, for an $n$-dimensional numerical integration using a Sobol' sequence (which I prefer to Halton), the rate of convergence is $(\log N)^n/N$. So $N=64$ gives you as much accuracy as 236 independent random variates in the one-dimensional case. I'll look at the specific code later on. Aug 6, 2018 at 16:49
• If I turn Halton sampling off as you did in your final code snippet, I find that setting R = 10, 50, 100, or 200 yields inaccurate estimates of sd.bus, but R = 300 or 500 is better. With Halton sampling on, R = 10, 50, 100, 200 yields sd.bus estimates of 0.65, 11.02, 11.06, and 25.6, respectively. There doesn't seem to be a clean recipe for obtaining accurate estimates of the true model. The correct prescription for a practitioner is unclear to me.. Aug 6, 2018 at 18:29
• I have checked the Halton code and run it separately, and it looks good, if slow because it's in R. I have also tried re-parameterizing the Halton call (halton=list(prime=3, drop=10), for example), and for some values of prime (e.g., 11) it works, for others (e.g., 2 (the default), 29) it fails. 29, for example, gives a reasonable sd.bus estimate but price and time are around -0.7. ATM, I have to agree with you; the correct prescription for a practitioner is unclear to me as well. Aug 6, 2018 at 19:09
• Is this a byproduct of an incorrect implementation? Or is this expected with the method of maximum simulated log-likelihood? I am really hoping to pin down the reason for the inaccurate estimates. Aug 6, 2018 at 21:42

I believe that the BFGS implementation is the culprit here. My first two clues were:

1. Calling mlogit() with the argument method='bhhh' instead of the default bfgs resulted in much more accurate estimates.
2. When I obtained inaccurate estimates from bfgs, the stop condition for the optimizer was last step couldn't find higher value, suggesting that the BFGS step was not an ascent direction.

I followed the methodology found in L-BFGS-B FORTRAN SUBROUTINES FOR LARGE SCALE BOUND CONSTRAINED OPTIMIZATION, which I quote here:

If the line search is unable to find a point with a sufficiently lower value of the objective after 20 evaluations of the objective function we conclude that the current direction is not useful. In this case all correction vectors are discarded and the iteration is restarted along the steepest descent direction

I updated the mlogit.optim() in the mlogit/R/mlogit.tools.R function so that the BFGS approximation of the inverse Hessian is reset to the identity if an ascent step is not found in the line search. I capped the maximum number of reset to 10 (in these tests, I never hit the max). The update looks like this (the first line in this copy/paste is unchanged):

    # eval the function and compute the gradient and the hessian
x <- eval(f, parent.frame())
if (is.null(x)){
if(method == 'bfgs' && num.bfgs.reset < 10) {
num.bfgs.reset <- num.bfgs.reset + 1
Hm1            <- diag(nrow(Hm1))
x              <- oldx
next # try again
} else {
## x is null if steptol is reached
code = 3
break
}
}


Running the same simulation as before, I get a much more accurate estimate of the standard deviation:

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data,
rpar = c(bus = "n"), R = 300, halton = NA, method = "bfgs")

Frequencies of alternatives:
car red.bus   train
0.1317  0.4084  0.4599

bfgs method
21 iterations, 0h:3m:50s
g'(-H)^-1g = 1.8E-06
successive function values within tolerance limits

Coefficients :
Estimate Std. Error  z-value Pr(>|z|)
price  -0.989864   0.026878 -36.8276   <2e-16 ***
time   -0.991994   0.032682 -30.3530   <2e-16 ***
bus    -0.118717   0.228376  -0.5198   0.6032
sd.bus 10.357554   0.847985  12.2143   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -8743.2

random coefficients
Min.   1st Qu.     Median       Mean  3rd Qu. Max.
bus -Inf -7.104781 -0.1187173 -0.1187173 6.867347  Inf


I get similar results if I change R, generate new datasets, etc.