mlogit package fails to recover synthetic mixed logit model

Question

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?

Answer 1

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 
gradient close to zero 

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.

Answer 2

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.