Page 1 of 1

Mixed logit model with error component variance heterogeneity

Posted: 12 May 2025, 12:42
by zalif
Hello,

I am trying to run a mixed logit model with error component and variance heterogeneity. I followed the example given here to include variance heterogeneity https://www.apollochoicemodelling.com/f ... means#p107 and the model runs, but while all the coefficients are estimated, their standard errors are just NAs. Previously, I ran a mixed logit model with error component (without variance heterogeneity) and a mixed logit model with variance heterogeneity (without error component), and they both worked fine. Could you please help me find the error in my code, or is there some other issue due to including both error components and variance heterogeneity in the same model?

Here is my code:

Code: Select all

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName="MXL_noprotest_error_component",
  modelDescr      = "Mixed logit model, no protest respondents, but with error component ",
  indivID         = "KMG.MID",  
  nCores          = 4,
  outputDirectory = "output",
  seed=12345
)


database = anketa2

apollo_beta = c(mu_b_fallow   = 0,
                sigma_b_fallow = 0.01,
                mu_b_meadow    = 0,
                sigma_b_meadow = 0.01,
                mu_b_lands    = 0,
                sigma_b_lands = 0,
                b_payment     = 0,
                mu_b_asc = 0,
                sigma_b_asc = 0,
                sigma_bsh = 0.01,
                lambda_AE_sig_fallow=0.1,
                lambda_AE_sig_meadow=0.1,
                lambda_AE_sig_lands=0.1,
                lambda_AE_sig_ASC=0.1)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c()

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "halton",
  interNDraws    = 500,
  interUnifDraws = c(),
  interNormDraws = c("draws_fallow","draws_meadow","draws_lands","draws_bsh", "draws_asc"),
  intraDrawsType = "halton",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()

  randcoeff[["b_fallow"]] =( mu_b_fallow + sigma_b_fallow * Q65^lambda_AE_sig_fallow*draws_fallow )
  randcoeff[["b_meadow"]] = ( mu_b_meadow + sigma_b_meadow *Q65^lambda_AE_sig_meadow * draws_meadow )
  randcoeff[["b_lands"]] = ( mu_b_lands + sigma_b_lands *Q65^lambda_AE_sig_lands * draws_lands )
  randcoeff[["b_asc"]] = ( mu_b_asc + sigma_b_asc*Q65^lambda_AE_sig_ASC * draws_asc  )
  randcoeff[["bsh"]] = ( sigma_bsh * draws_bsh )

  return(randcoeff)
}

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){

  ### Function initialisation: do not change the following three commands
  ### 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: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["alt1"]] = b_fallow * fallow_1 + b_meadow * meadow_1 + b_lands * lands_1 + b_payment * payment_1
  V[["alt2"]] = b_fallow * fallow_2 + b_meadow * meadow_2 + b_lands * lands_2 + b_payment * payment_2
  V[["alt3"]] = b_asc+bsh

  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = choice,
    utilities     = V
  )

  ### Compute probabilities using MNL model
  P[["model"]] = apollo_mnl(mnl_settings, functionality)

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

  ### Average across inter-individual draws
  P = apollo_avgInterDraws(P, apollo_inputs, functionality)

  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}

# ################################################################# #
#### MODEL ESTIMATION                                            ####
# ################################################################# #

model5 = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
summary(model5)
# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

output5<-apollo_modelOutput(model5,  modelOutput_settings=list(printCorr=T))
Thank you for your help!

Re: Mixed logit model with error component variance heterogeneity

Posted: 19 May 2025, 18:35
by stephanehess
Hi

happy to help, but really need to see your results to diagnose what the problem is

Stephane

Re: Mixed logit model with error component variance heterogeneity

Posted: 20 May 2025, 07:53
by zalif
Hi,

Thanks a lot, below you can find the code and results.

CODE

Code: Select all


### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName="MXL_noprotest_error_component",
  modelDescr      = "Mixed logit model, no protest respondents, but with error component ",
  indivID         = "KMG.MID",  
  nCores          = 4,
  outputDirectory = "output",
  seed=12345
)


database = anketa2

apollo_beta = c(mu_b_fallow   = 0,
                sigma_b_fallow = 0.01,
                mu_b_meadow    = 0,
                sigma_b_meadow = 0.01,
                mu_b_lands    = 0,
                sigma_b_lands = 0,
                b_payment     = 0,
                mu_b_asc = 0,
                sigma_b_asc = 0,
                sigma_bsh = 0.01,
                lambda_AE_sig_fallow=0.1,
                lambda_AE_sig_meadow=0.1,
                lambda_AE_sig_lands=0.1,
                lambda_AE_sig_ASC=0.1)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c()

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "halton",
  interNDraws    = 500,
  interUnifDraws = c(),
  interNormDraws = c("draws_fallow","draws_meadow","draws_lands","draws_bsh", "draws_asc"),
  intraDrawsType = "halton",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()

  randcoeff[["b_fallow"]] =( mu_b_fallow + sigma_b_fallow * Q65^lambda_AE_sig_fallow*draws_fallow )
  randcoeff[["b_meadow"]] = ( mu_b_meadow + sigma_b_meadow *Q65^lambda_AE_sig_meadow * draws_meadow )
  randcoeff[["b_lands"]] = ( mu_b_lands + sigma_b_lands *Q65^lambda_AE_sig_lands * draws_lands )
  randcoeff[["b_asc"]] = ( mu_b_asc + sigma_b_asc*Q65^lambda_AE_sig_ASC * draws_asc  )
  randcoeff[["bsh"]] = ( sigma_bsh * draws_bsh )

  return(randcoeff)
}

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){

  ### Function initialisation: do not change the following three commands
  ### 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: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["alt1"]] = b_fallow * fallow_1 + b_meadow * meadow_1 + b_lands * lands_1 + b_payment * payment_1+bsh
  V[["alt2"]] = b_fallow * fallow_2 + b_meadow * meadow_2 + b_lands * lands_2 + b_payment * payment_2+bsh
  V[["alt3"]] = b_asc

  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = choice,
    utilities     = V
  )

  ### Compute probabilities using MNL model
  P[["model"]] = apollo_mnl(mnl_settings, functionality)

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

  ### Average across inter-individual draws
  P = apollo_avgInterDraws(P, apollo_inputs, functionality)

  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}

# ################################################################# #
#### MODEL ESTIMATION                                            ####
# ################################################################# #

model5 = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
m5<-summary(model5)
# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

output5<-apollo_modelOutput(model5)

RESULTS

Code: Select all

> model5 = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
Preparing user-defined functions.

Testing likelihood function...

Overview of choices for MNL model component :
                                    alt1    alt2    alt3
Times available                  3160.00 3160.00 3160.00
Times chosen                     1029.00  928.00 1203.00
Percentage chosen overall          32.56   29.37   38.07
Percentage chosen when available   32.56   29.37   38.07


Pre-processing likelihood function...
Creating cluster...
Preparing workers for multithreading...

Testing influence of parameters
Starting main estimation

BGW using analytic model derivatives supplied by caller...


Iterates will be written to: 
 output/MXL_noprotest_error_component_iterations.csv
    it    nf     F            RELDF    PRELDF    RELDX    MODEL stppar
     0     1 3.463607086e+03
     1     2 3.372046602e+03 2.644e-02 6.220e-03 2.25e-01   G   3.49e+00
     2     3 3.038470325e+03 9.892e-02 4.697e-02 6.83e-01   G   1.70e-01
     3     4 2.775230068e+03 8.664e-02 5.601e-02 6.78e-01   G   0.00e+00
     4     5 2.536363359e+03 8.607e-02 5.698e-02 4.58e-01   G   0.00e+00
     5     6 2.389848325e+03 5.777e-02 4.577e-02 1.56e-01   G   -2.68e-03
     6     7 2.300778569e+03 3.727e-02 3.011e-02 2.46e-01   G   -6.71e-04
     7     8 2.273480860e+03 1.186e-02 2.425e-02 2.50e-01   S   0.00e+00
     8    10 2.270554439e+03 1.287e-03 3.589e-03 4.43e-02  S-G  -8.80e-05
     9    11 2.265368559e+03 2.284e-03 1.772e-03 2.41e-02   G   -1.10e-05
    10    12 2.263056786e+03 1.020e-03 9.104e-04 2.58e-02   G   -2.75e-06
    11    13 2.262794468e+03 1.159e-04 4.100e-04 1.49e-02   G   -6.88e-07
    12    14 2.262385390e+03 1.808e-04 2.919e-04 1.08e-02   S   0.00e+00
    13    15 2.262210328e+03 7.738e-05 1.301e-04 6.11e-03   S   0.00e+00
    14    16 2.262149838e+03 2.674e-05 2.556e-05 2.25e-03   S   0.00e+00
    15    17 2.262148414e+03 6.295e-07 5.379e-07 9.04e-04   S   0.00e+00
    16    18 2.262148268e+03 6.476e-08 6.165e-08 3.85e-04   S   0.00e+00
    17    19 2.262148264e+03 1.624e-09 1.450e-09 3.56e-05   S   0.00e+00
    18    20 2.262148264e+03 1.119e-10 1.104e-10 1.23e-05   S   0.00e+00
    19    21 2.262148264e+03 8.803e-12 9.342e-12 2.41e-06   S   0.00e+00

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

Estimated parameters:
                        Estimate
mu_b_fallow             -0.32190
sigma_b_fallow           0.32187
mu_b_meadow             -0.22844
sigma_b_meadow           0.23559
mu_b_lands              -0.36102
sigma_b_lands            0.31957
b_payment                0.01103
mu_b_asc                -1.90116
sigma_b_asc             -0.46255
sigma_bsh                2.71115
lambda_AE_sig_fallow    -0.06614
lambda_AE_sig_meadow    -0.05015
lambda_AE_sig_lands     -0.05175
lambda_AE_sig_ASC     -349.22272

Final LL: -2262.1483

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%
WARNING: Singular Hessian, cannot calculate s.e. 
Could not write hessian to a file.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a
  saddle point! 
Computing score matrix...

> output5<-apollo_modelOutput(model5)
Model run by ZAlif using Apollo 0.3.5 on R 4.3.3 for Windows.
Please acknowledge the use of Apollo by citing Hess & Palma (2019)
  DOI 10.1016/j.jocm.2019.100170
  www.ApolloChoiceModelling.com

Model name                                  : MXL_noprotest_error_component
Model description                           : Mixed logit model, no protest respondents, but with error component 
Model run at                                : 2025-05-20 08:21:52.477646
Estimation method                           : bgw
Model diagnosis                             : Relative function convergence
Optimisation diagnosis                      : Saddle point found
     hessian properties                     : Some eigenvalues are positive and others negative
     maximum eigenvalue                     : 0
     reciprocal of condition number         : not calculated (Hessian is not negative definite)
Number of individuals                       : 395
Number of rows in database                  : 3160
Number of modelled outcomes                 : 3160

Number of cores used                        :  4 
Number of inter-individual draws            : 500 (halton)

LL(start)                                   : -3463.61
LL at equal shares, LL(0)                   : -3471.61
LL at observed shares, LL(C)                : -3453.4
LL(final)                                   : -2262.15
Rho-squared vs equal shares                  :  0.3484 
Adj.Rho-squared vs equal shares              :  0.3444 
Rho-squared vs observed shares               :  0.345 
Adj.Rho-squared vs observed shares           :  0.3415 
AIC                                         :  4552.3 
BIC                                         :  4637.11 

Estimated parameters                        : 14
Time taken (hh:mm:ss)                       :  00:03:41.1 
     pre-estimation                         :  00:00:27.59 
     estimation                             :  00:00:28.74 
     post-estimation                        :  00:02:44.77 
Iterations                                  :  19  

Unconstrained optimisation.

Estimates:
                        Estimate        s.e.   t.rat.(0)    Rob.s.e. Rob.t.rat.(0)
mu_b_fallow             -0.32190          NA          NA          NA            NA
sigma_b_fallow           0.32187          NA          NA          NA            NA
mu_b_meadow             -0.22844          NA          NA          NA            NA
sigma_b_meadow           0.23559          NA          NA          NA            NA
mu_b_lands              -0.36102          NA          NA          NA            NA
sigma_b_lands            0.31957          NA          NA          NA            NA
b_payment                0.01103          NA          NA          NA            NA
mu_b_asc                -1.90116          NA          NA          NA            NA
sigma_b_asc             -0.46255          NA          NA          NA            NA
sigma_bsh                2.71115          NA          NA          NA            NA
lambda_AE_sig_fallow    -0.06614          NA          NA          NA            NA
lambda_AE_sig_meadow    -0.05015          NA          NA          NA            NA
lambda_AE_sig_lands     -0.05175          NA          NA          NA            NA
lambda_AE_sig_ASC     -349.22272          NA          NA          NA            NA
In addition, I would like to check whether I am specifying error components correctly. I am using them to incorporate status quo effects, as per Scarpa et al., 2005. To do so, I have specified a random coefficient "bsh" and added it to the first two alternatives to induce correlation between them, but leaving it out of the third, SQ alternative. Is this the correct way to specify this model? (In my original post I just added an error component to the SQ option, but upon further reading I now assume this was incorrect.)

Thanks a lot for your help!

Re: Mixed logit model with error component variance heterogeneity

Posted: 28 May 2025, 13:10
by stephanehess
Hi

what is the range of Q65?

In relation to the error component term, what you are doing would create correlation between the first two alternatives, but also heteroskedasticity. If you want to avoid the latter, then you could also add an additional error component to alt3, defined as

Code: Select all

randcoeff[["sig3"]] = ( sigma_bsh * draws_alt3 )
where draws_alt3 needs to be added to apollo_draws. Then maybe switch from halton to mlhs given the high dimensionality

Stephane

Re: Mixed logit model with error component variance heterogeneity

Posted: 04 Jun 2025, 07:49
by zalif
Hello,

Q65 is a binary variable with values 0/1.

Is it preferred to avoid heteroskedasticity or is this more up to personal preference?

Thank you!

Živa

Re: Mixed logit model with error component variance heterogeneity

Posted: 13 Jun 2025, 16:52
by stephanehess
Hi

in relation to heteroskedasticity, this is an empirical question.

But in relation to the other point, you seem to be estimating an exponent on a 0/1 variable, and this then of course goes wrong when the exponent goes negative

Stephane