Mixed Logit Model for BWS: differences between Apollo and Stata estimations + baseline option issue in Apollo

[tan421]

Hello,

I am estimating a Mixed Logit model for Best-Worst Scaling (BWS) data and would appreciate some guidance on two issues I am encountering when using the apollo package, especially in comparison with Stata’s cmxtmixlogit.

My BWS Design:
- 6 total alternatives: item_1, item_2, item_3, item_4, item_5, item_6
- Each respondent faced 10 choice questions
- Each question contains 3 options
- Each alternative appears five times in every choice set
I estimated the same model specification in (below code using item_4 as baseline):
- Stata with cmxtmixlogit using 500 Halton draws, as it only allows Halton:

Code: Select all

cmxtmixlogit chosen_pair, ///
>         rand(item_1 item_2 item_3 item_5 item_6, normal) ///
>         basealternative(2) intpoint(500) intmethod(halton) ///
>         ltolerance(1e-6) gtolerance(1e-5) tolerance(1e-5)

- Apollo package using both 500 Halton and 2,000 Sobol draws:

Code: Select all

apollo_setWorkDir()

# Initialize code
apollo_initialise()

# Set core controls
apollo_control = list(
  modelName       = "BWS_Item 4 Baseline",
  modelDescr      = "Best-worst model",
  indivID         = "id", 
  panelData       = TRUE,
  mixing          = TRUE,
  outputDirectory = "BWS Output"
)
# Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(b_item_1_m                      = 0, b_item_1_s                      = 1,
                b_item _2_m                       = 0, b_item _2_s                       = 1,
                b_item _3_m                       = 0, b_item_3_s                       = 1,
                b_item_5_m            	       = 0, b_Item_5_s            = 1,
                b_Item_6_m  	               = 0, b_Item_6_s  = 1,
                mu_worst                      = 1
                )

apollo_fixed = c() # if nothing is fixed

# Draw names must match usage in apollo_randCoeff
drawNames = c(
  "draw_item _1", 
  "draw_item _2",
  "draw_item _3", 
  "draw_item_5"
  "draw_item_6"
)

apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 2000,
  interUnifDraws = c(),
  interNormDraws = drawNames,
  intraDrawsType = "sobol",
  intraNDraws = 0
)

apollo_randCoeff = function(apollo_beta, apollo_inputs){
  
  randcoeff = list()

  randcoeff[["b_item_1"]] = b_item_1_m + b_item_1_s * draw_item_1
  randcoeff[["b_item_2"]] = b_item_2_m + b_item_2_s * draw_item_2
  randcoeff[["b_item_3"]] = b_item_3_m + b_item_3_s * draw_item_3
  randcoeff[["b_item_5"]] = b_item_5_m + b_item_5_s * draw_item_5
  randcoeff[["b_item_6"]] = b_item_6_m + b_item_6_s * draw_item_6
  
  return(randcoeff)
}
apollo_inputs = apollo_validateInputs()
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){

  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))

  ### Create list of probabilities P
  P = list()
  
  ### List of utilities
  V = list()
  
  V[["alt1"]] = (
      b_item_1 * (alt1=="item_1") + 
      b_item_2 * (alt1=="item_2") + 
      b_item_3 * (alt1=="item_3 ") + 
      b_item_5 * (alt1=="item_5 ") +
      b_item_6 * (alt1=="item_6 ")
                      )
  V[["alt2"]] = (
      b_item_1 * (alt2=="item_1") + 
      b_item_2 * (alt2=="item_2") + 
      b_item_3 * (alt2=="item_3 ") + 
      b_item_5 * (alt2=="item_5 ") +
      b_item_6 * (alt2=="item_6 ")
                      )
  V[["alt3"]] = (
      b_item_1 * (alt3=="item_1") + 
      b_item_2 * (alt3=="item_2") + 
      b_item_3 * (alt3=="item_3 ") + 
      b_item_5 * (alt3=="item_5 ") +
      b_item_6 * (alt3=="item_6 ")
                      )

  ### Compute probabilities for "best" choice 
  mnl_settings_best = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = apollo_inputs$database$best,
    utilities     = V
  )
  
  P[["choice_best"]] = apollo_mnl(mnl_settings_best, functionality)
  
  ### Compute probabilities for "worst" choice 
  mnl_settings_worst = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=(best!=1), alt2=(best!=2), alt3=(best!=3)),
    choiceVar     = apollo_inputs$database$worst, 
    utilities     = list(alt1 = -mu_worst * V[["alt1"]],
                         alt2 = -mu_worst * V[["alt2"]],
                         alt3 = -mu_worst * V[["alt3"]])
  )
  
  P[["choice_worst"]] = apollo_mnl(mnl_settings_worst, functionality)

  ### Combined model
  P = apollo_combineModels(P, apollo_inputs, functionality)

  ### Take product across observation for same individual
  P = apollo_panelProd(P, apollo_inputs, functionality)
  
  P = apollo_avgInterDraws(P, apollo_inputs, functionality)

  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}
model <- apollo_estimate(apollo_beta, 
                         apollo_fixed, 
                         apollo_probabilities, 
                         apollo_inputs, 
                         estimate_settings = list(
                         estimationRoutine = "bfgs",
                         maxIterations = 2000,
                         maxFunctionEvaluations = 10000,
                         writeIter = TRUE
                         )
)                        
apollo_modelOutput(model, modelOutput_settings = list(printPVal=2))
apollo_saveOutput(model, saveOutput_settings = list(printPVal=2))

Question 1: Differences in coefficient estimates between Apollo and Stata
Stata Output:

Code: Select all

Fitting fixed parameter model:

Fitting full model:

Iteration 0:  Log simulated-likelihood = -15857.711  (not concave)
Iteration 1:  Log simulated-likelihood = -15497.441  
Iteration 2:  Log simulated-likelihood = -15404.274  (not concave)
Iteration 3:  Log simulated-likelihood = -15177.231  
Iteration 4:  Log simulated-likelihood = -14777.356  
Iteration 5:  Log simulated-likelihood = -14749.626  
Iteration 6:  Log simulated-likelihood = -14749.258  
Iteration 7:  Log simulated-likelihood = -14749.257  
Iteration 8:  Log simulated-likelihood = -14749.257  

Mixed logit choice model                     Number of obs        =     60,180
                                             Number of cases      =     10,030
Panel variable: id                           Number of panels     =      1,003

Time variable: questio~d                     Cases per panel: min =         10
                                                              avg =       10.0
                                                              max =         10

Alternatives variable: alt_id                Alts per case:   min =          6
                                                              avg =        6.0
                                                              max =          6
Integration sequence:          Halton
Integration points:               500             Wald chi2(5)    =    1050.22
Log simulated-likelihood = -14749.257             Prob > chi2     =     0.0000

------------------------------------------------------------------------------
 chosen_pair | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
alt_id       |
      item_1 |   1.710708    .069516    24.61   0.000     1.574459    1.846957
      item_2 |   1.966666   .0852595    23.07   0.000      1.79956    2.133771
      item_3 |   .4849733   .0487527     9.95   0.000     .3894197    .5805269
      item_5 |   .0383341   .0368168     1.04   0.298    -.0338255    .1104938
      item_6 |   .2811412    .039947     7.04   0.000     .2028465    .3594358
-------------+----------------------------------------------------------------
/Normal      |
   sd(item_1)|   1.781794    .069258                      1.651093    1.922842
   sd(item_2)|   2.162087   .0831851                      2.005043    2.331431
   sd(item_3)|   1.193251   .0513884                      1.096666    1.298344
   sd(item_5)|   .6628692   .0477565                      .5755763    .7634011
   sd(item_6)|   .8275077   .0455975                       .742795    .9218814
-------------+----------------------------------------------------------------
1            |
       _cons |   .0383401   .0381446     1.01   0.315    -.0364218    .1131021
-------------+----------------------------------------------------------------
2            |  (base alternative)
-------------+----------------------------------------------------------------
3            |
       _cons |  -.5118027   .0512539    -9.99   0.000    -.6122585    -.411347
-------------+----------------------------------------------------------------
4            |
       _cons |  -.5782508   .0448611   -12.89   0.000    -.6661769   -.4903247
-------------+----------------------------------------------------------------
5            |
       _cons |  -.6193203   .0516247   -12.00   0.000    -.7205029   -.5181377
-------------+----------------------------------------------------------------
6            |
       _cons |  -1.123095   .0539245   -20.83   0.000    -1.228785   -1.017405
------------------------------------------------------------------------------

Apollo Output:

Code: Select all

## Testing likelihood function...
## Apollo found a model component of type MNL without a componentName. The
## name was set to "choice_best" by default.
##
## Overview of choices for MNL model component choice_best:
## alt1 alt2 alt3
## Times available 10030.00 10030.00 10030.00
## Times chosen 3443.00 3423.00 3164.00
## Percentage chosen overall 34.33 34.13 31.55
## Percentage chosen when available 34.33 34.13 31.55
##
## Apollo found a model component of type MNL without a componentName. The
## name was set to "choice_worst" by default.
##
## Overview of choices for MNL model component choice_worst:
## alt1 alt2 alt3
## Times available 6587.00 6607.00 6866.00
## Times chosen 3962.00 2942.00 3126.00
## Percentage chosen overall 39.50 29.33 31.17
## Percentage chosen when available 60.15 44.53 45.53
##
##
## Pre-processing likelihood function...
##
## Testing influence of parameters
## Starting main estimation
## Initial function value: -16356.51
## Initial gradient value:
## b_item_1_m b_item_1_s
## 616.77024 467.78378
## b_item_2_m b_item_2_s
## 629.12382 576.24770
## b_item_3_m b_item_3_s
## 131.82731 67.10042
## b_item_5_m b_item_5_s
## -307.71328 -61.30739
## b_item_6_m b_item_6_s
## -73.20799 -143.91249
## initial value 16356.514309
## iter 2 value 15285.523654
## iter 3 value 15129.775765
## iter 4 value 15107.136178
## iter 5 value 15100.422868
## iter 6 value 15099.434879
## iter 7 value 15098.221534
## iter 8 value 15059.103342
## iter 9 value 15031.129941
## iter 10 value 15029.095172
## iter 11 value 15028.207653
## iter 12 value 14999.064818
## iter 13 value 14998.527119
## iter 14 value 14981.145767
## iter 15 value 14977.605412
## iter 16 value 14976.562408
## iter 17 value 14976.287103
5
## iter 18 value 14976.278394
## iter 18 value 14976.278226
## iter 18 value 14976.278224
## final value 14976.278224
## converged
##
## Estimated parameters with approximate standard errors from BHHH matrix:
## Estimate BHHH se BHH t-ratio (0)
## b_item_1_m 2.16789 0.09748 22.2390
## b_item_1_s 2.25115 0.08792 25.6037
## b_item_2_m 2.24786 0.10875 20.6709
## b_item_2_s 2.48235 0.09390 26.4371
## b_item_3_m 0.90960 0.06671 13.6357
## b_item_3_s 1.64694 0.05501 29.9383
## b_item_5_m 0.03397 0.04195 0.8099
## b_item_5_s 0.73056 0.04665 15.6601
## b_item_6_m 0.51937 0.05022 10.3413
## b_item_6_s 1.12248 0.04385 25.5998
##
## Final LL: -14976.2782
##
## Calculating log-likelihood at equal shares (LL(0)) for applicable
## models...
## Calculating log-likelihood at observed shares from estimation data
## (LL(c)) for applicable models...
## Calculating LL of each model component...
## Calculating other model fit measures
## Computing covariance matrix using numerical jacobian of analytical
## gradient.
## 0%....25%....50%....75%.100%
## Negative definite Hessian with maximum eigenvalue: -61.223143
##
## Please acknowledge the use of Apollo by citing Hess & Palma (2019) -
## doi.org/10.1016/j.jocm.2019.100170
apollo_modelOutput(model, modelOutput_settings = list(printPVal=2))
## Model run by myname using Apollo 0.3.5 on R 4.4.3 for Darwin.
## Please acknowledge the use of Apollo by citing Hess & Palma (2019)
## DOI 10.1016/j.jocm.2019.100170
## www.ApolloChoiceModelling.com
##
## Model name : BWS_Item 4 Baseline
## Model description : Best-worst model
## Model run at : 2025-09-12 14:37:44.090866
## Estimation method : bfgs
## Model diagnosis : successful convergence
## Optimisation diagnosis : Maximum found
## hessian properties : Negative definite
## maximum eigenvalue : -61.223143
## reciprocal of condition number : 0.0840966
## Number of individuals : 1003
## Number of rows in database : 10030
## Number of modelled outcomes : 20060
## choice_best : 10030
6
## choice_worst : 10030
##
## Number of cores used : 1
## Number of inter-individual draws : 2000 (sobol)
##
## LL(start) : -16356.51
## LL (whole model) at equal shares, LL(0) : -17971.35
## LL (whole model) at observed shares, LL(C) : -21944.51
## LL(final, whole model) : -14976.28
## Rho-squared vs equal shares : 0.1667
## Adj.Rho-squared vs equal shares : 0.1661
## Rho-squared vs observed shares : 0.3175
## Adj.Rho-squared vs observed shares : 0.3173
## AIC : 29972.56
## BIC : 30044.69
##
## LL(0,choice_best) : -11019.08
## LL(final,choice_best) : -9157.77
## LL(0,choice_worst) : -6952.27
## LL(final,choice_worst) : -6283.74
##
## Estimated parameters : 10
## Time taken (hh:mm:ss) : 00:19:37.54
## pre-estimation : 00:00:36.66
## estimation : 00:05:58.73
## post-estimation : 00:13:2.15
## Iterations : 21
##
## Unconstrained optimisation.
##
## Estimates:
## Estimate s.e. t.rat.(0) p(2-sided)
## b_item_1_m 2.16789 0.08913 24.3238 0.0000
## b_item_1_s 2.25115 0.08427 26.7135 0.0000
## b_item_2_m 2.24786 0.09491 23.6831 0.0000
## b_item_2_s 2.48235 0.09207 26.9624 0.0000
## b_item_3_m 0.90960 0.06552 13.8825 0.0000
## b_item_3_s 1.64694 0.06433 25.6006 0.0000
## b_item_5_m 0.03397 0.04159 0.8170 0.4139
## b_item_5_s 0.73056 0.05549 13.1654 0.0000
## b_item_6_m 0.51937 0.05043 10.2999 0.0000
## b_item_6_s 1.12248 0.05480 20.4824 0.0000
## Rob.s.e. Rob.t.rat.(0) p(2-sided)
## b_item_1_m 0.13530 16.0228 0.0000
## b_item_1_s 0.10459 21.5238 0.0000
## b_item_2_m 0.14154 15.8819 0.0000
## b_item_2_s 0.11675 21.2625 0.0000
## b_item_3_m 0.08006 11.3610 0.0000
## b_item_3_s 0.08392 19.6247 0.0000
## b_item_5_m 0.04222 0.8047 0.4210
## b_item_5_s 0.06712 10.8837 0.0000
## b_item_6_m 0.06131 8.4719 0.0000
## b_item_6_s 0.07366 15.2385 0.0000

I wonder whether there are known reasons why Apollo and Stata cmxtmixlogit may produce substantially different estimates for the same BWS mixed logit specification.

Question 2: Baseline-specific identification error in Apollo (but not in Stata)
When estimating the model in Apollo using 2000 Sobol draws, I encounter an identification error unless I use item_4 as the baseline.
If I try to use any of the other 5 options (e.g., item_1) as the baseline,

Code: Select all

# Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
                b_item _2_m                       = 0, b_item _2_s                       = 1,
                b_item _3_m                       = 0, b_item_3_s                       = 1,
b_item_4_m                      = 0, b_item_4_s                      = 1,
                b_item_5_m            	       = 0, b_Item_5_s            = 1,
                b_Item_6_m  	               = 0, b_Item_6_s  = 1,
                mu_worst                      = 1
                )

apollo_fixed = c() # if nothing is fixed

# Draw names must match usage in apollo_randCoeff
drawNames = c(
  "draw_item _2",
  "draw_item _3", 
"draw_item _4",
  "draw_item_5"
  "draw_item_6"
)

apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 2000,
  interUnifDraws = c(),
  interNormDraws = drawNames,
  intraDrawsType = "sobol",
  intraNDraws = 0
)

apollo_randCoeff = function(apollo_beta, apollo_inputs){
  
  randcoeff = list()

  randcoeff[["b_item_2"]] = b_item_2_m + b_item_2_s * draw_item_2
  randcoeff[["b_item_3"]] = b_item_3_m + b_item_3_s * draw_item_3
randcoeff[["b_item_1"]] = b_item_4_m + b_item_4_s * draw_item_4
  randcoeff[["b_item_5"]] = b_item_5_m + b_item_5_s * draw_item_5
  randcoeff[["b_item_6"]] = b_item_6_m + b_item_6_s * draw_item_6
  
  return(randcoeff)
}
apollo_inputs = apollo_validateInputs()
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){

  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))

  ### Create list of probabilities P
  P = list()
  
  ### List of utilities
  V = list()
  
  V[["alt1"]] = (
      b_item_2 * (alt1=="item_2") + 
      b_item_3 * (alt1=="item_3 ") + 
      b_item_4 * (alt1=="item_4") + 
      b_item_5 * (alt1=="item_5 ") +
      b_item_6 * (alt1=="item_6 ")
                      )
  V[["alt2"]] = (
      b_item_2 * (alt2=="item_2") + 
      b_item_3 * (alt2=="item_3 ") + 
       b_item_4 * (alt2=="item_4") + 
      b_item_5 * (alt2=="item_5 ") +
      b_item_6 * (alt2=="item_6 ")
                      )
  V[["alt3"]] = (
      b_item_2 * (alt3=="item_2") + 
      b_item_3 * (alt3=="item_3 ") + 
      b_item_4 * (alt3=="item_4") + 
      b_item_5 * (alt3=="item_5 ") +
      b_item_6 * (alt3=="item_6 ")
                      )

Apollo returns the following error:

Code: Select all

Error in apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, ...) : 
SPECIFICATION ISSUE - Parameters b_item_4_m, b_item_4_s do not influence the log-likelihood of your model!

This error occurs in both restricted and unrestricted models. The same issue arises regardless of whether I use Halton or Sobol draws in Apollo. However, when I estimate the same model in Stata (cmxtmixlogit), I can change the baseline option freely and do not encounter any analogous error.
The paper I followed used the least popularly chosen option as its baseline. However, in my case, the least-chosen option is not item_4, choice shares vary across items, and item_4 does not appear to be special based on observed choice frequencies.

Thank you very much. I look forward to your insights. Let me know if more details would be helpful.

[stephanehess]

Hi

as a first step, can you confirm that you get exactly the same results with MNL in both packages?

Stephane

[tan421]

Hello Dr. Hess,

Thanks for your reply. The MNL outputs from both Stata mlogit and Apollo are relatively comparable only for the coefficients of items 1-3. Please see the outputs below:

stata

Code: Select all

********************************************************************************
. //Use item_4 as baseline
. *****************
. mlogit chosen_pair item_1 item_2 item_3 item_5 item_6

Iteration 0:  Log likelihood = -27114.774  
Iteration 1:  Log likelihood =  -25524.13  
Iteration 2:  Log likelihood = -25454.221  
Iteration 3:  Log likelihood = -25454.001  
Iteration 4:  Log likelihood = -25454.001  

Multinomial logistic regression                        Number of obs =  60,180
                                                       LR chi2(5)    = 3321.54
                                                       Prob > chi2   =  0.0000
Log likelihood = -25454.001                            Pseudo R2     =  0.0612

------------------------------------------------------------------------------
 chosen_pair | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
0            |  (base outcome)
-------------+----------------------------------------------------------------
1            |
      item_1 |   1.027262   .0261566    39.27   0.000     .9759963    1.078528
      item_2 |   1.039888   .0261739    39.73   0.000     .9885883    1.091188
      item_3 |   .5241555   .0255685    20.50   0.000     .4740421    .5742689
      item_5 |    .012725   .0252241     0.50   0.614    -.0367133    .0621634
      item_6 |   .2767956   .0253679    10.91   0.000     .2270755    .3265157
       _cons |  -1.752187   .0123412  -141.98   0.000    -1.776376   -1.727999
------------------------------------------------------------------------------

Apollo:

Code: Select all

Overview of choices for MNL model component choice_best:
                                     alt1     alt2     alt3
Times available                  10030.00 10030.00 10030.00
Times chosen                      3443.00  3423.00  3164.00
Percentage chosen overall           34.33    34.13    31.55
Percentage chosen when available    34.33    34.13    31.55

Apollo found a model component of type MNL without a componentName. The
  name was set to "choice_worst" by default.

Overview of choices for MNL model component choice_worst:
                                    alt1    alt2    alt3
Times available                  6587.00 6607.00 6866.00
Times chosen                     3962.00 2942.00 3126.00
Percentage chosen overall          39.50   29.33   31.17
Percentage chosen when available   60.15   44.53   45.53


Pre-processing likelihood function...

Testing influence of parameters
Starting main estimation

BGW using analytic model derivatives supplied by caller...


Iterates will be written to: 
 BWS MNL Output/TIS_Apollo Code_iterations.csv
    it    nf     F            RELDF    PRELDF    RELDX    MODEL stppar
     0     1 1.797134748e+04
     1     3 1.786391165e+04 5.978e-03 5.953e-03 1.00e+00   G   1.93e+01
     2     5 1.701213519e+04 4.768e-02 4.475e-02 8.58e-01   G   7.62e-01
     3     6 1.654714421e+04 2.733e-02 2.687e-02 4.81e-01   S   0.00e+00
     4     7 1.654321390e+04 2.375e-04 2.311e-04 3.21e-02   S   0.00e+00
     5     8 1.654320078e+04 7.929e-07 7.897e-07 1.24e-03   S   0.00e+00
     6     9 1.654320077e+04 3.966e-10 4.107e-10 2.73e-05   G   0.00e+00
     7    10 1.654320077e+04 2.460e-12 2.588e-12 1.71e-06   G   0.00e+00

***** Relative function convergence *****

Estimated parameters with approximate standard errors from BHHH matrix:
                              Estimate     BHHH se BHH t-ratio (0)
b_item_1                       1.20606     0.02387          50.526
b_item_2                       1.21631     0.02245          54.169
b_item_3                       0.61623     0.02385          25.842
b_item_5		      	      0.05610     0.03154           1.779
b_item_6			      0.41082     0.02674          15.362

Final LL: -16543.2008

Calculating log-likelihood at equal shares (LL(0)) for applicable
  models...
Calculating log-likelihood at observed shares from estimation data (LL(c))
  for applicable models...
Calculating LL of each model component...
Calculating other model fit measures
Computing covariance matrix using numerical jacobian of analytical
  gradient.
 0%....25%....50%...100%
Negative definite Hessian with maximum eigenvalue: -327.934805

I am only able to compare the MNL outputs using item_4 as a baseline since Apollo shows the same error message if I use others:

Error in apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, : SPECIFICATION ISSUE - Parameter b_item_4 does not influence the log-likelihood of your model!

[stephanehess]

Hi

they are clearly not the same model. You need to get exactly the same results for MNL. Can you show us the specification for both packages? I note the Stata one has a constant for example

Regarding item_4, this will work, but you need to include it in the utility functions of course

Stephane