# mlogit + logitr packages fail to recover true estimates of mixed logit random coefficient model

I am running Monte-Carlo simulations on a simple DGP of a mixed logit random coefficient model to check if the mlogit and logitr packages are able to recover the true parameters which are determined by the DGP. The DGP is set as follows $$u_{ntj} = \beta_{x1} * x1_{njt} + \beta_{x2_{n}} * x2_{njt} + \epsilon_{njt}$$ where $$u_{ijt}$$ describes the latent utility of individual $$n$$ at time $$t$$ of alternative $$j$$, $$x1_{njt}$$ is an explanatory variable with fixed coefficient, $$x2_{njt}$$ is an explanatory variable with random coefficient / slope and $$\epsilon_{njt} \sim Gumbel(0, {\pi}^2/6)$$ is the error term. $$\beta_{x1}$$ is the fixed coefficient belonging to $$x1_{njt}$$ and $$\beta_{x2_{n}}$$ is the random coefficient belonging to $$x2$$.

For ease of simplicity, I set $$(x1,x2) \sim \mathbf{N}(\mathbf{0}, \mathbf{\Sigma})$$ where $$\mathbf{\Sigma} = \mathbf{I_2}$$ , $$\beta_{x2_{n}} \sim \mathbf{N}(\mu_{x2}, \sigma_{x2})$$ and $$\beta_{x1} = \mu_{x2} = \sigma_{x2} = 1$$.

I set the number of Monte-Carlo repetitions / iterations $$I = 100$$ and varied the number of individuals $$N$$ between $$1000, 5000, 10000$$.

The R code I used is as follows:

# Load necessary packages
library(dplyr)
library(mlogit)
library(logitr)
library(EnvStats)
library(parallel)

# Define functions
estimate_mixed_logit <- function(seed, N, T, J, beta_x1, mean_x2, sd_x2, nummultistarts){

# Setting parameters
set.seed(seed)

id <- seq(1:N)
t <- seq(1:T)
alt <- seq(1:J)

# Creation of the data frame
df <- expand.grid(id, t, alt)
names(df) <- c('id', 't', 'alt')

df <- df %>%
group_by(id) %>%
mutate(ran_coef = rnorm(1,mean_x2,sd_x2)) %>%
group_by(id,t,alt) %>%
mutate(x1 = rnorm(1),
x2 = rnorm(1),
e = rgevd(1),
y_star = case_when(alt == 1 ~ 0 + beta_x1 * x1 + ran_coef * x2 + e,
alt == 2 ~ 0 + beta_x1 * x1 + ran_coef * x2 + e,
alt == 3 ~ 0 + beta_x1 * x1 + ran_coef * x2 + e)
) %>%
group_by(id,t) %>%
mutate(max_y_star = max(y_star),
choice = ifelse(max_y_star == y_star,1,0),
id_t = cur_group_id()) %>%
arrange(id,t,alt) %>%
ungroup()

# Mlogit estimation

## Reshaping the data in mlogit format
df_cm <- dfidx(df,
idx = list(c('id_t','id')),
chid.var = 'id',
alt.var = 'alt',
group.var = 't',
choice = 'choice',
drop.index = TRUE)

## Estimation
rpar <- c('n')
names(rpar) <- c('x2')

mlm <- mlogit(choice ~ 1 + x1 + x2, data = df_cm, panel = TRUE, correlation = FALSE, rpar = rpar)

# Logitr estimation
logitr <- logitr(
data     = df,
outcome  = 'choice',
obsID    = 'id_t',
panelID  = 'id',
pars     = c('x1', 'x2'),
randPars = c(x2 = 'n'),
numMultiStarts = nummultistarts
)

return(c(mlm$$coefficients[['x1']], mlm$$coefficients[['x2']], mlm$$coefficients[['sd.x2']], logitr$$coefficients[['x1']], logitr$$coefficients[['x2']], logitr$$coefficients[['sd_x2']], logitr$status)) } # Set global parameters R <- 100 N <- 1000 T <- 20 J <- 3 beta_x1 <- 1 mean_x2 <- 1 sd_x2 <- 1 nummultistarts <- 30 # Run Monte Carlo simulations res <- mclapply(1:R, mc.cores = 39, N = N, T = T, J = J, beta_x1 = beta_x1, mean_x2 = mean_x2, sd_x2 = sd_x2, nummultistarts = nummultistarts, FUN <- function(r, N, T, J, beta_x1, mean_x2, sd_x2, nummultistarts){ return(estimate_mixed_logit(seed = r, N = N, T = T, J = J, beta_x1 = beta_x1, mean_x2 = mean_x2, sd_x2 = sd_x2, nummultistarts = nummultistarts )) } ) # Summarize the results res <- as.data.frame(do.call(rbind, res)) names(res) <- c('mlogit_beta_x1', 'mlogit_mean_x2', 'mlogit_sd_x2', 'logitr_beta_x1', 'logitr_mean_x2', 'logitr_sd_x2', 'logitr_status')  which returns the following results: • N = 1000 summary(res1000) mlogit_beta_x1 mlogit_mean_x2 mlogit_sd_x2 logitr_beta_x1 logitr_mean_x2 logitr_sd_x2 logitr_status Min. :0.9701 Min. :0.8299 Min. :0.9425 Min. :0.9720 Min. :0.9312 Min. :-1.154 Min. :3 1st Qu.:0.9886 1st Qu.:0.9569 1st Qu.:1.0069 1st Qu.:0.9902 1st Qu.:1.0215 1st Qu.: 1.042 1st Qu.:3 Median :0.9986 Median :0.9757 Median :1.0370 Median :1.0005 Median :1.0499 Median : 1.065 Median :3 Mean :0.9984 Mean :0.9724 Mean :1.0352 Mean :1.0004 Mean :1.0476 Mean : 1.021 Mean :3 3rd Qu.:1.0085 3rd Qu.:0.9903 3rd Qu.:1.0568 3rd Qu.:1.0096 3rd Qu.:1.0738 3rd Qu.: 1.090 3rd Qu.:3 Max. :1.0326 Max. :1.0692 Max. :1.1166 Max. :1.0335 Max. :1.1218 Max. : 1.155 Max. :3  • N = 5000 summary(res5000) mlogit_beta_x1 mlogit_mean_x2 mlogit_sd_x2 logitr_beta_x1 logitr_mean_x2 logitr_sd_x2 logitr_status Min. :0.9784 Min. :0.9314 Min. :0.9788 Min. :0.9805 Min. :0.9504 Min. :-1.115 Min. :3 1st Qu.:0.9929 1st Qu.:0.9592 1st Qu.:1.0191 1st Qu.:0.9949 1st Qu.:1.0424 1st Qu.: 1.052 1st Qu.:3 Median :0.9986 Median :0.9720 Median :1.0308 Median :1.0003 Median :1.0534 Median : 1.059 Median :3 Mean :0.9979 Mean :0.9715 Mean :1.0315 Mean :0.9998 Mean :1.0525 Mean : 1.038 Mean :3 3rd Qu.:1.0025 3rd Qu.:0.9836 3rd Qu.:1.0434 3rd Qu.:1.0042 3rd Qu.:1.0614 3rd Qu.: 1.071 3rd Qu.:3 Max. :1.0119 Max. :1.0197 Max. :1.0691 Max. :1.0135 Max. :1.1003 Max. : 1.092 Max. :3  • N = 10000 summary(res10000) mlogit_beta_x1 mlogit_mean_x2 mlogit_sd_x2 logitr_beta_x1 logitr_mean_x2 logitr_sd_x2 logitr_status Min. :0.9839 Min. :0.9412 Min. :1.002 Min. :0.9858 Min. :1.027 Min. :1.031 Min. :3 1st Qu.:0.9930 1st Qu.:0.9596 1st Qu.:1.021 1st Qu.:0.9952 1st Qu.:1.044 1st Qu.:1.051 1st Qu.:3 Median :0.9960 Median :0.9664 Median :1.028 Median :0.9978 Median :1.051 Median :1.060 Median :3 Mean :0.9962 Mean :0.9666 Mean :1.029 Mean :0.9981 Mean :1.050 Mean :1.059 Mean :3 3rd Qu.:0.9997 3rd Qu.:0.9732 3rd Qu.:1.037 3rd Qu.:1.0018 3rd Qu.:1.056 3rd Qu.:1.067 3rd Qu.:3 Max. :1.0050 Max. :0.9916 Max. :1.055 Max. :1.0069 Max. :1.082 Max. :1.082 Max. :3  or shown graphically: whereby I ensured that all maximum likelihood estimations achieved convergence, see logitr_status for instance. My results suggest that both mlogit and logitr estimates for $$x2$$ are biased and even inconsistent since they do not converge to the true parameters for $$x2$$, i.e. $$\mu_{x2}$$ and $$\sigma_{x2}$$ when increasing $$N$$. When running the Monte-Carlo simulation in Stata using the exact same data ($$N = 1000$$) and the command cmmixlogit choice x1, random(x2), I obtain unbiased results, see the following: summary(stata1000) stata_beta_x1 stata_mean_x2 stata_sd_x2 stata_converged Min. :0.9716 Min. :0.8698 Min. :0.8857 Min. :1 1st Qu.:0.9918 1st Qu.:0.9731 1st Qu.:0.9728 1st Qu.:1 Median :0.9997 Median :0.9990 Median :1.0033 Median :1 Mean :1.0015 Mean :1.0007 Mean :1.0063 Mean :1 3rd Qu.:1.0107 3rd Qu.:1.0292 3rd Qu.:1.0345 3rd Qu.:1 Max. :1.0407 Max. :1.1016 Max. :1.1531 Max. :1  Does anybody have an explanation for this outcome or is there a conceptual misunderstanding in my set-up? My best guess is that Stata’s numeric optimization outperforms the ones of both R packages. However, I am surprised that this seems to lead to biased + inconsistent estimates (in finite samples up to $$N = 10000$$). P.S.: I am aware of this post (mlogit package fails to recover synthetic mixed logit model) but changing the optimization method / algorithm to the BHHH algorithm (for the mlogit package) or the number of draws mlogit R argument set to 100) didn’t affect the outcome of bias and inconsistency in my case. Update August 24th 2023: Thanks @jhelvy for your reply and support! According to https://cran.r-project.org/web/packages/mlogit/mlogit.pdf mlogit's default value of draws is 40. Stata’s cmmixlogit is by default using the Hammersley instead of the Halton sequence for the generation of point sets used in the Monte Carlo integration. Moreover, it uses a flexible number of integration points / draws: $$500 + \lfloor 2.5 \sqrt{N_c (\ln(r + 5) + v)}\rfloor$$, where $$N_c$$ is the number of cases, $$r$$ is the number of random coefficients in the model, and $$v$$ is the number of variance parameters, see https://www.stata.com/manuals/cmcmmixlogit.pdf. As you have suggested I replicated my analysis ($$N = 1000$$, $$R = 100$$, same DGP as described above) for varying number of draws. The following chart shows the mean absolute deviation across the $$R$$ number of iterations to the true parameter value for different levels of draws. The significance test refers to an one sample t test testing the null hypothesis that the respective mean of the absolute bias is statistically significant from zero with $$\alpha =0.05$$. TRUE / FALSE are representing that the respective p-value of this test below / above 10%. My conclusion from this additional analysis is that 200 draws seems to be an adequate number to offset the tradeoff between bias and computation time. However, I am wondering that even for 500 draws there is still some statistically significant bias and that logitr seems to be performing less than mlogit. @jhelvy Do you have an explanation for this? My understanding is that mlogit uses pseudo-random numbers in contrast to Halton sequences. Is this maybe explaining the difference although I would expect it to be the opposite, see Chapter 9 of the textbook by Kenneth Train. Finally, you’re absolutely right regarding the necessity to take the absolute value for the sd_x2 parameter. I changed it in my code and ensured that this does not affect my previous results. Update August 28th 2023: I run some additional Monte-Carlo simulations where I (i) altered the distribution of x2 from normal to log-normal and (ii) switched from preference to willingness-to-pay (WTP) space estimation mererly using the logitr package. For the latter, I altered the DGP by introducing a new variable $$x3 \sim \mathbf{N}(0, 1)$$ with $$\beta_{x3_{n}} = -1$$ which can be viewed as the price which enters negatively into the utility function: $$u_{ntj} = \beta_{x1} * x1_{njt} + \beta_{x2_{n}} * x2_{njt} + \beta_{x3_{n}} * x3_{njt} + \epsilon_{njt}$$ Rewriting to the WTP space expression gives: $$u_{ntj} = \lambda_n ( \omega_{x1} * x1_{njt} + \omega_{x2_{n}} * x2_{njt} - x3_{njt}) + \epsilon_{njt}$$ where $$\lambda_n = \beta_{x3_{n}}$$ de facto since I set $$\sigma_{\epsilon} = 1$$, see above. The following table shows the results by specification where for all simulations $$N = 1000$$ and $$R = 100$$ and only the estimates of converged iterations $$r$$ have been taken into account. Space Distribution $$\beta_{x2}$$ Distribution $$\lambda$$ Number of Draws Number of logitr multistarts mlogit bias $$\beta_{x1}$$ mlogit bias $$\mu_{\beta_{x2}}$$ mlogit bias $$\sigma_{\beta_{x2}}$$ mlogit convergence rate logitr bias $$\beta_{x1}$$ logitr bias $$\mu_{\beta_{x2}}$$ logitr bias $$\sigma_{\beta_{x2}}$$ logitr bias $$\lambda$$ / $$\mu_{\lambda}$$ logitr bias $$\sigma_{\lambda}$$ logitr convergence rate Preference $$\mathbf{Lognormal}(1,1)$$ 200 30 0% 5% 1% 100% 0% 1% 1% 52% Preference $$\mathbf{Lognormal}(1,1)$$ 500 200 0% 8% 0% 100% 0% 0% 0% 64% WTP $$\mathbf{N}(1,1)$$ fixed to $$\lambda$$ = 1 200 30 0% 4% 1% 0% 100% WTP $$\mathbf{N}(1,1)$$ $$\mathbf{N}(1,1)$$ 200 30 0% -7% 6% -5% 819% 100% WTP $$\mathbf{N}(1,1)$$ $$\mathbf{Lognormal}(1,1)$$ 200 30 0% -4% 4% -2% 2% 100% WTP $$\mathbf{Lognormal}(1,1)$$ fixed to $$\lambda$$ = 1 200 30 0% -2% 1% 0% 95% WTP $$\mathbf{Lognormal}(1,1)$$ $$\mathbf{N}(1,1)$$ 200 30 0% 2% 1% -9% 831% 68% WTP $$\mathbf{Lognormal}(1,1)$$ $$\mathbf{Lognormal}(1,1)$$ 500 200 0% -2% 1% -1% 6% 22% From my point of view, there are three key findings: • mlogit estimates for $$\mu_{\beta_{x2}}$$ in the preference space when $$x2 \sim \mathbf{Lognormal}(1,1)$$ seem to considerably deviate. • logitr estimates for $$\mu_{\beta_{x2}}$$ and $$\sigma_{\beta_{x2}}$$ in the WTP space when $$x2 \sim \mathbf{N}(1,1)$$ seem to considerably deviate independent from the distribution of $$\lambda$$. • logitr estimates for the nuisance parameter of $$\lambda$$ can be heavily biased which are, however, of less interest. Any thoughts on these additional insights? • Just to add I am using R version 4.2.3 and the following package versions: mlogit (1.1.1), logitr (1.1.0) and EnvStats (2.8.0). Commented Aug 21, 2023 at 11:12 • Thanks for extending the analysis. I ran my code and also still got some bias at 200 draws, so clearly something is off. I still suspect this may have to do with the draws. I also expected that Halton would perform better (not worse) than pseudo-random numbers, but perhaps either my implementation of Halton is off or they just don't do as good a job as people thought. There is an option to use sobol draws in logitr - could you try changing drawType = 'sobol' and replicate your results for logitr? Commented Aug 30, 2023 at 10:33 • On the WTP estimates, I suspect this still may be the draws for the most part. For the ones where you're getting crazy bias levels on the SD parameters (~800%), I suspect that it is because you're modeling the scale parameter with a normal distribution, which really should never be done. The scale parameter must be strictly positive, and if some of the draws are negative then the model won't converge. I should probably not even allow that - I could just insert a check where the user cannot use 'n' for the randScale argument. Commented Aug 30, 2023 at 10:40 • Thanks for your reply @jhelvy! As requested, I additionally run simulations with number of draws equal to 100 and 200 and compared (i) mlogit with pseudo random numbers draws (the default) to mlogit with halton draws and (ii) logitr with halton draws to logitr with sobol draws. For number of draws equal to 200 there is no real difference detectable in the estimates among the two groups. For number of draws equal to 100, there is again no difference for the mlogit estimators but the sobol draw variant of logitr performs slight less than the halton variant. Commented Sep 3, 2023 at 11:58 • I also had a look how Stata’s cmmixlogit performs if I force it to just use 100 draws of either the Hammersley sequence (the default) or the Halton sequence. Both are unbiased than mlogit or logitr variants even for this number of draws. Commented Sep 3, 2023 at 11:58 ## 1 Answer logitr author here. Thanks for bringing this up - it's super important. This issue was also posted as an issue on the package repo here. My immediate thought was that the number of Halton draws being used simply isn't enough to do a good enough job approximating the normal distribution. By default, logitr uses just 50 draws (same as mlogit if I'm not mistaken), and that is a trade off made between resolution and speed. It looks like that choice though may be resulting in a bias of ~5%. I just replicated the same simulation for just logitr setting numDraws = 200 and got much better results. Here are my results for N = 1000 and R = 100:  logitr_beta_x1 logitr_mean_x2 logitr_sd_x2 logitr_status Min. 0.9620884 0.8614243 0.9376236 3 1st Qu. 0.9894498 0.9764979 0.9912670 3 Median 0.9980751 1.0080507 1.0130822 3 Mean 0.9984994 1.0000567 1.0099384 3 3rd Qu. 1.0075022 1.0288857 1.0280742 3 Max. 1.0353412 1.1213550 1.0746411 3  A helpful exercise might be to replicate this using numDraws as a variable and scale it from 50 to maybe 500 and see at what point we start seeing consistent results. That will probably also change with more random parameters, etc, but it might at least provide some justification for setting the default number of draws. I'm happy to increase that number if it makes sense. Btw, below is the code I used for my simulation. I re-wrote it myself as a totally separate approach to try to rule out other issues, like running things in parallel, etc. I thought there might be some chance that the seed setting was causing problems with the parallel processing (which can sometimes happen), but my code runs in serial and I was still getting the bias result with numDraws = 50. Also, note that for the sd_x2 parameter, I keep the absolute value rather than the raw estimate because it can sometimes be negative (same for mlogit). If that happens, they should be interpreted as positive. This is done because by allowing a negative SD parameter the optimization doesn't have to be constrained, which makes it much easier for the solver. Keeping negative estimates might skew the outcome of the simulation exercise. library(logitr) library(dplyr) estimate_mixed_logit <- function( n, t, j, beta_x1, mean_x2, sd_x2, nummultistarts ) { # Same betas for x2 across population beta_x2 <- rnorm(n, mean_x2, sd_x2) # Generate the X data matrix X <- as.data.frame(matrix(rnorm(n*t*j*2), ncol = 2)) names(X) <- c('x1', 'x2') X$$id <- rep(seq(n), each = t*j) X$$alt_id <- rep(seq(j), t*n) X$$q_id <- rep(rep(seq(t), each = j), n) X$$obs_id <- rep(seq(n*t), each = j) # Compute utility by respondent X_resp <- split(X[c('x1', 'x2')], X$id)
utility <- list()
for (i in 1:n) {
beta <- c(beta_x1, beta_x2[i])
utility[[i]] <- as.matrix(X_resp[[i]]) %*% beta
}
V <- unlist(utility)

# Compute probabilities of choosing each alternative

expV <- exp(V)
sumExpV <- rowsum(expV, group = X$$obs_id, reorder = FALSE) reps <- table(X$$obs_id)
p <- expV / sumExpV[rep(seq_along(reps), reps),]

# Simulate choices according to probabilities

p_obs <- split(p, X$$obs_id) choices <- list() for (i in seq_len(length(p_obs))) { choice <- rep(0, reps[i]) choice[sample(seq(reps[i]), 1, prob = p_obs[[i]])] <- 1 choices[[i]] <- choice } X$$choice <- unlist(choices)

# Logitr estimation

suppressMessages(
model <- logitr(
data     = as.data.frame(X),
outcome  = 'choice',
obsID    = 'obs_id',
panelID  = 'id',
pars     = c('x1', 'x2'),
randPars = c(x2 = 'n'),
numDraws = 200,
numMultiStarts = nummultistarts
))

return(c(
model$$coefficients[['x1']], model$$coefficients[['x2']],
abs(model$$coefficients[['sd_x2']]), model$$status
))

}

R <- 100
results <- list()
for (i in 1:R) {
results[[i]] <- estimate_mixed_logit(
n = 1000,
t = 20,
j = 3,
beta_x1 = 1,
mean_x2 = 1,
sd_x2 = 1,
nummultistarts = 30
)
}

# Summarize the results

res <- as.data.frame(do.call(rbind, results))
names(res) <- c('logitr_beta_x1', 'logitr_mean_x2', 'logitr_sd_x2', 'logitr_status')

apply(res, 2, summary)

• Thanks @jhelvy for your reply. I have posted two updates. Would be fantastic to get your view on that! Commented Aug 28, 2023 at 11:45