I am estimating a mixed multinomial logistic regression on a stated preference discrete choice experiment with two options and an opt-out option with the latest version of Apollo (0.2.6). I am able to estimate the code without error, until I include intra- draws. When I include intra-individual draws I encounter the following error:
------------------------------------------------------------------------------------------------------------------
Computing covariance matrix using numerical methods (maxLik). This may take a while,
no progress bar displayed.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a
saddle point!
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
Calculating LL of each model component...
Warning message:
In sqrt(diag(varcov)) : NaNs produced
------------------------------------------------------------------------------------------------------------------
I have edited this original post to include the code and output from models with only inter-individual draws first, and the code adding intra-individual draws second:
Inter-individual draws only:
Code: Select all
> #### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
>
> #### Check cores
> parallel::detectCores()
[1] 4
>
> #### Set core controls
> apollo_control = list(
+ modelName ="MXL_1",
+ modelDescr ="Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
+ reporting time continous and quadratic specification,
+ main effects only, add one random attribute parameter at a time",
+ indivID ="RID",
+ mixing = TRUE,
+ nCores = 2,
+ analyticGrad = TRUE # Needs to be explicit for models with inter-intra draws (requires extra RAM).
+ )
>
> # inter-individual draws are used for variation across individuals
> # intra-individual draws for variation across observations for the same individual
>
> #### Set database
> database = df1
>
> ###################################################################
> #### DEFINE MODEL PARAMETERS ####
> ###################################################################
>
> # Utility U_ij = sum[1 to K](beta_jk * X_jk) + epsilon_ij
> #
> # where U_ij is the ultity of profile j for individual i
> # X_jk where k = 1...K are attributes in the dce
>
> ### Vector of parameters, including any that are kept fixed in estimation
> apollo_beta = c(asc_a = 0,
+ asc_b = 0,
+ b_time = 0,
+ b_time_sq = 0,
+ b_anonymous = 0,
+ mu_b_investigate = -3,
+ sigma_b_investigate_inter = -0.01,
+ b_amnesty = 0,
+ b_violate_low = 0,
+ b_violate_med = 0,
+ b_violate_high = 0 )
>
> ### 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("b_violate_low")
>
> #####################################################################
> #### DEFINE RANDOM COMPONENTS ####
> #####################################################################
>
> ### Set parameters for generating draws
> apollo_draws = list(
+ interDrawsType = "halton",
+ interNDraws = 500,
+ interNormDraws = c("draws_investigate_inter"),
+ interUnifDraws = c(),
+ intraDrawsType = "halton",
+ intraNDraws = 500,
+ intraNormDraws = c(), #"draws_investigate_intra"
+ intraUnifDraws = c()
+ )
>
> ### Create random parameters
> apollo_randCoeff = function(apollo_beta, apollo_inputs){
+ randcoeff = list()
+
+ randcoeff[["b_investigate"]] = exp( mu_b_investigate
+ + sigma_b_investigate_inter * draws_investigate_inter)
+
+ return(randcoeff)
+ }
>
> #####################################################################
> #### GROUP AND VALIDATE INPUTS ####
> #####################################################################
>
> apollo_inputs = apollo_validateInputs()
Several observations per individual detected based on the value of RID. Setting panelData in
apollo_control set to TRUE.
All checks on apollo_control completed.
All checks on database completed.
Generating inter-individual draws . Done
>
> #####################################################################
> #### 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[['a']] = asc_a + b_anonymous * anonymous_a + b_time * time_a + b_time_sq * time_a * time_a + b_investigate * investigate_a + b_violate_low * (violate_a == 0) + b_violate_med * (violate_a == 1) + b_violate_high * (violate_a == 2) + b_amnesty * amnesty_a
+ V[['b']] = asc_b + b_anonymous * anonymous_b + b_time * time_b + b_time_sq * time_b * time_b + b_investigate * investigate_b + b_violate_low * (violate_b == 0) + b_violate_med * (violate_b == 1) + b_violate_high * (violate_b == 2) + b_amnesty * amnesty_b
+ V[['c']] = 0
+
+ ### Define settings for MNL model component
+ mnl_settings = list(
+ alternatives = c(a=1, b=2, c=3),
+ avail = list(a=1, b=1, c=1),
+ choiceVar = choice,
+ V = 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 ####
> #####################################################################
>
> # ### Optional speedTest
> # speedTest_settings=list(
> # nDrawsTry = c(100, 200, 500),
> # nCoresTry = 1:10,
> # nRep = 10
> # )
> #
> # apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
>
> start = Sys.time()
> MXL_1 = apollo_estimate(apollo_beta, apollo_fixed,
+ apollo_probabilities, apollo_inputs,
+ estimate_settings=list(hessianRoutine="maxLik"))
Testing likelihood function...
Overview of choices for MNL model component :
a b c
Times available 8032.00 8032.00 8032.00
Times chosen 3412.00 3314.00 1306.00
Percentage chosen overall 42.48 41.26 16.26
Percentage chosen when available 42.48 41.26 16.26
Pre-processing likelihood function...
Preparing workers for multithreading...
Testing influence of parameters..........
Starting main estimation
Initial function value: -8765.926
Initial gradient value:
asc_a asc_b b_time b_time_sq
712.5388242 614.3715108 3207.9337189 9846.3914569
b_anonymous mu_b_investigate sigma_b_investigate_inter b_amnesty
1255.7097360 56.4212148 -0.6178251 1333.0880225
b_violate_med b_violate_high
513.4127116 1110.8300423
initial value 8765.926317
iter 2 value 8582.288541
iter 3 value 7770.770249
iter 4 value 7648.697753
iter 5 value 7616.882131
iter 6 value 7607.037904
iter 7 value 7604.399293
iter 8 value 7570.647068
iter 9 value 7552.916107
iter 10 value 7460.863274
iter 11 value 7445.479063
iter 12 value 7441.466730
iter 13 value 7428.754494
iter 14 value 7399.990025
iter 15 value 7383.376540
iter 16 value 7341.186160
iter 17 value 7340.229107
iter 18 value 7317.309733
iter 19 value 7274.567382
iter 20 value 7259.182689
iter 21 value 7238.101461
iter 22 value 7230.494075
iter 23 value 7228.760639
iter 24 value 7227.822498
iter 25 value 7227.664340
iter 26 value 7227.642890
iter 27 value 7227.618981
iter 28 value 7227.550642
iter 29 value 7227.391476
iter 30 value 7227.181059
iter 31 value 7222.266657
iter 32 value 7220.683561
iter 33 value 7219.954010
iter 34 value 7219.106938
iter 35 value 7218.046924
iter 36 value 7214.099235
iter 37 value 7213.265947
iter 38 value 7213.190276
iter 39 value 7212.664686
iter 40 value 7211.618149
iter 41 value 7211.324080
iter 42 value 7211.320379
iter 42 value 7211.320371
iter 42 value 7211.320371
final value 7211.320371
converged
Additional convergence test using scaled estimation. Parameters will be scaled by their current
estimates and additional iterations will be performed.
initial value 7211.320371
iter 1 value 7211.320370
final value 7211.320370
converged
Estimated parameters:
Estimate
asc_a -1.12378
asc_b -1.15800
b_time 0.37853
b_time_sq -0.08825
b_anonymous 0.71798
mu_b_investigate -0.88344
sigma_b_investigate_inter -1.25881
b_amnesty 0.81616
b_violate_low 0.00000
b_violate_med 0.94421
b_violate_high 1.13150
Computing covariance matrix using numerical methods (maxLik). This may take a while, no progress bar
displayed.
Negative definite Hessian with maximum eigenvalue: -8.168464
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
Calculating LL of each model component...
>
> end = Sys.time()
> print(paste0("Run time in seconds: ",end-start))
[1] "Run time in seconds: 11.4913820823034"
> # ----------------------------------------------------------------- #
> #---- FORMATTED OUTPUT (TO SCREEN) ----
> # ----------------------------------------------------------------- #
>
> apollo_modelOutput(MXL_1)
Model run using Apollo for R, version 0.2.6 on Darwin by aliceellyson
www.ApolloChoiceModelling.com
Model name : MXL_1
Model description : Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
reporting time continous and quadratic specification,
main effects only, add one random attribute parameter at a time
Model run at : 2022-01-04 13:23:50
Estimation method : bfgs
Model diagnosis : successful convergence
Number of individuals : 1004
Number of rows in database : 8032
Number of modelled outcomes : 8032
Number of cores used : 2
Number of inter-individual draws : 500 (halton)
LL(start) : -8765.926
LL(0) : -8824.054
LL(C) : -8227.245
LL(final) : -7211.32
Rho-square (0) : 0.1828
Adj.Rho-square (0) : 0.1816
Rho-square (C) : 0.1235
Adj.Rho-square (C) : 0.1223
AIC : 14442.64
BIC : 14512.55
Estimated parameters : 10
Time taken (hh:mm:ss) : 00:11:29.48
pre-estimation : 00:01:26.79
estimation : 00:04:17.55
post-estimation : 00:05:45.14
Iterations : 45
Min abs eigenvalue of Hessian : 8.168464
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_a -1.12378 0.19922 -5.641 0.21454 -5.238
asc_b -1.15800 0.20010 -5.787 0.21406 -5.410
b_time 0.37853 0.20181 1.876 0.20985 1.804
b_time_sq -0.08825 0.04059 -2.174 0.04228 -2.087
b_anonymous 0.71798 0.03621 19.827 0.04420 16.245
mu_b_investigate -0.88344 0.11710 -7.544 0.10801 -8.179
sigma_b_investigate_inter -1.25881 0.11075 -11.367 0.09834 -12.800
b_amnesty 0.81616 0.03637 22.443 0.04541 17.973
b_violate_low 0.00000 NA NA NA NA
b_violate_med 0.94421 0.07390 12.776 0.08257 11.436
b_violate_high 1.13150 0.05179 21.849 0.05685 19.904
Code: Select all
#### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
>
> #### Check cores
> parallel::detectCores()
[1] 40
>
> #### Set core controls
> apollo_control = list(
+ modelName ="MXL_1",
+ modelDescr ="Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
+ reporting time continous and quadratic specification,
+ main effects only, add one random attribute parameter at a time",
+ indivID ="RID",
+ mixing = TRUE,
+ nCores = 4,
+ analyticGrad = TRUE # Needs to be explicit for models with inter-intra draws (requires extra RAM).
+ )
>
> # inter-individual draws are used for variation across individuals
> # intra-individual draws for variation across observations for the same individual
>
> #### Set database
> database = df1
>
> ###################################################################
> #### DEFINE MODEL PARAMETERS ####
> ###################################################################
>
> # Utility U_ij = sum[1 to K](beta_jk * X_jk) + epsilon_ij
> #
> # where U_ij is the ultity of profile j for individual i
> # X_jk where k = 1...K are attributes in the dce
>
> ### Vector of parameters, including any that are kept fixed in estimation
> apollo_beta = c(asc_a = 0,
+ asc_b = 0,
+ b_time = 0,
+ b_time_sq = 0,
+ b_anonymous = 0,
+ mu_b_investigate = -3,
+ sigma_b_investigate_inter = -0.01,
+ sigma_b_investigate_intra = -0.01,
+ #b_investigate = 0,
+ b_amnesty = 0,
+ b_violate_low = 0,
+ b_violate_med = 0,
+ b_violate_high = 0 )
>
>
> ### 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("b_violate_low")
>
> #####################################################################
> #### DEFINE RANDOM COMPONENTS ####
> #####################################################################
>
> ### Set parameters for generating draws
> apollo_draws = list(
+ interDrawsType = "halton",
+ interNDraws = 100,
+ interNormDraws = c("draws_investigate_inter"),
+ interUnifDraws = c(),
+ intraDrawsType = "halton",
+ intraNDraws = 100,
+ intraNormDraws = c("draws_investigate_intra"),
+ intraUnifDraws = c()
+ )
>
> ### Create random parameters
> apollo_randCoeff = function(apollo_beta, apollo_inputs){
+ randcoeff = list()
+
+ # randcoeff[["asc_a"]] = -exp(mu_asc_a + sigma_asc_a * draws_asc_a)
+ # randcoeff[["asc_b"]] = -exp(mu_asc_b + sigma_asc_b * draws_asc_b)
+ # randcoeff[["b_anonymous"]] = -exp(mu_b_anonymous + sigma_b_anonymous_inter * draws_anonymous_inter )
+ # randcoeff[["b_investigate"]] = -exp( mu_b_investigate + sigma_b_investigate * draws_investigate )
+ randcoeff[["b_investigate"]] = exp( mu_b_investigate
+ + sigma_b_investigate_inter * draws_investigate_inter
+ + sigma_b_investigate_intra * draws_investigate_intra )
+ # randcoeff[["b_amnesty"]] = -exp( mu_b_amnesty + sigma_b_amnesty * draws_amnesty )
+ # randcoeff[["b_violate_med"]] = -exp( mu_b_violate_med + sigma_b_violate_med * draws_violate_med )
+ # randcoeff[["b_violate_high"]] = -exp( mu_b_violate_high + sigma_b_violate_high * draws_violate_high )
+
+ return(randcoeff)
+ }
>
> #####################################################################
> #### GROUP AND VALIDATE INPUTS ####
> #####################################################################
>
> apollo_inputs = apollo_validateInputs()
Several observations per individual detected based on the value of RID. Setting
panelData in apollo_control set to TRUE.
All checks on apollo_control completed.
All checks on database completed.
Generating inter-individual draws . Done
Generating intra-individual draws . Done
>
> #####################################################################
> #### 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[['a']] = asc_a + b_anonymous * anonymous_a + b_time * time_a + b_time_sq * time_a * time_a + b_investigate * investigate_a + b_violate_low * (violate_a == 0) + b_violate_med * (violate_a == 1) + b_violate_high * (violate_a == 2) + b_amnesty * amnesty_a
+ V[['b']] = asc_b + b_anonymous * anonymous_b + b_time * time_b + b_time_sq * time_b * time_b + b_investigate * investigate_b + b_violate_low * (violate_b == 0) + b_violate_med * (violate_b == 1) + b_violate_high * (violate_b == 2) + b_amnesty * amnesty_b
+ V[['c']] = 0
+
+ ### Define settings for MNL model component
+ mnl_settings = list(
+ alternatives = c(a=1, b=2, c=3),
+ avail = list(a=1, b=1, c=1),
+ choiceVar = choice,
+ V = V
+ )
+
+ ### Compute probabilities using MNL model
+ P[['model']] = apollo_mnl(mnl_settings, functionality)
+
+ ### Average across intra-individual draws
+ P = apollo_avgIntraDraws(P, apollo_inputs, 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 ####
> #####################################################################
>
> # ### Optional speedTest
> # speedTest_settings=list(
> # nDrawsTry = c(100, 200, 500),
> # nCoresTry = 1:10,
> # nRep = 10
> # )
> #
> # apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
>
> start = Sys.time()
> MXL_1 = apollo_estimate(apollo_beta, apollo_fixed,
+ apollo_probabilities, apollo_inputs,
+ estimate_settings=list(hessianRoutine="maxLik"))
Testing likelihood function...
Overview of choices for MNL model component :
a b c
Times available 8032.00 8032.00 8032.00
Times chosen 3412.00 3314.00 1306.00
Percentage chosen overall 42.48 41.26 16.26
Percentage chosen when available 42.48 41.26 16.26
Pre-processing likelihood function...
Preparing workers for multithreading...
Testing influence of parameters...........
Starting main estimation
Initial function value: -8765.923
Initial gradient value:
asc_a asc_b b_time
712.5375479 614.3704347 3207.9274789
b_time_sq b_anonymous mu_b_investigate
9846.3712588 1255.7085376 56.4241527
sigma_b_investigate_inter sigma_b_investigate_intra b_amnesty
-0.6257522 -0.5877284 1333.0868410
b_violate_med b_violate_high
513.4122345 1110.8289942
initial value 8765.923215
iter 2 value 8582.286139
iter 3 value 7770.768788
iter 4 value 7648.695940
iter 5 value 7616.880165
iter 6 value 7607.035810
iter 7 value 7604.396709
iter 8 value 7570.621983
iter 9 value 7552.892617
iter 10 value 7460.841522
iter 11 value 7445.452187
iter 12 value 7441.430071
iter 13 value 7428.692010
iter 14 value 7399.774756
iter 15 value 7398.391329
iter 16 value 7368.342956
iter 17 value 7356.778666
iter 18 value 7320.951826
iter 19 value 7311.605807
iter 20 value 7295.858801
iter 21 value 7295.324351
iter 22 value 7294.606567
iter 23 value 7290.550494
iter 24 value 7287.732732
iter 25 value 7286.091869
iter 26 value 7263.040792
iter 27 value 7240.433843
iter 28 value 7226.679718
iter 29 value 7217.269867
iter 30 value 7212.553853
iter 31 value 7211.230939
iter 32 value 7211.039130
iter 33 value 7211.025968
iter 34 value 7211.024206
iter 35 value 7211.019441
iter 36 value 7211.009142
iter 37 value 7210.998412
iter 37 value 7210.998339
final value 7210.998339
converged
Additional convergence test using scaled estimation. Parameters will be scaled by
their current estimates and additional iterations will be performed.
initial value 7210.998339
iter 2 value 7210.997260
iter 3 value 7210.996914
iter 4 value 7210.996713
iter 5 value 7210.996415
iter 6 value 7210.995609
iter 7 value 7210.995316
iter 7 value 7210.995283
iter 8 value 7210.994858
iter 9 value 7210.994734
iter 9 value 7210.994690
iter 9 value 7210.994682
final value 7210.994682
converged
Estimated parameters:
Estimate
asc_a -1.12508
asc_b -1.15906
b_time 0.37878
b_time_sq -0.08832
b_anonymous 0.71892
mu_b_investigate -0.90300
sigma_b_investigate_inter -1.27222
sigma_b_investigate_intra -0.17337
b_amnesty 0.81700
b_violate_low 0.00000
b_violate_med 0.94518
b_violate_high 1.13297
Computing covariance matrix using numerical methods (maxLik). This may take a while,
no progress bar displayed.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a
saddle point!
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
Calculating LL of each model component...
Warning message:
In sqrt(diag(varcov)) : NaNs produced
>
> end = Sys.time()
> print(paste0("Run time in seconds: ",end-start))
[1] "Run time in seconds: 4.84306277420786"
> # ----------------------------------------------------------------- #
> #---- FORMATTED OUTPUT (TO SCREEN) ----
> # ----------------------------------------------------------------- #
>
> apollo_modelOutput(MXL_1)
Model run using Apollo for R, version 0.2.6 on Windows by aellyson
www.ApolloChoiceModelling.com
Model name : MXL_1
Model description : Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
reporting time continous and quadratic specification,
main effects only, add one random attribute parameter at a time
Model run at : 2022-01-04 12:21:15
Estimation method : bfgs
Model diagnosis : successful convergence
Number of individuals : 1004
Number of rows in database : 8032
Number of modelled outcomes : 8032
Number of cores used : 4
Number of inter-individual draws : 100 (halton)
Number of intra-individual draws : 100 (halton)
LL(start) : -8765.923
LL(0) : -8824.054
LL(C) : -8227.245
LL(final) : -7210.995
Rho-square (0) : 0.1828
Adj.Rho-square (0) : 0.1816
Rho-square (C) : 0.1235
Adj.Rho-square (C) : 0.1222
AIC : 14443.99
BIC : 14520.89
Estimated parameters : 11
Time taken (hh:mm:ss) : 04:50:35.02
pre-estimation : 00:31:48.92
estimation : 01:49:42.13
post-estimation : 02:29:3.97
Iterations : 50
Min abs eigenvalue of Hessian : 6.733263
Some eigenvalues of Hessian are positive, indicating potential problems!
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_a -1.12508 0.17426 -6.456 0.17795 -6.3224
asc_b -1.15906 0.17542 -6.607 0.17802 -6.5108
b_time 0.37878 0.17466 2.169 0.16470 2.2998
b_time_sq -0.08832 0.03504 -2.521 0.03309 -2.6692
b_anonymous 0.71892 0.03594 20.003 0.04431 16.2265
mu_b_investigate -0.90300 0.07715 -11.705 0.11395 -7.9243
sigma_b_investigate_inter -1.27222 0.06758 -18.825 0.10503 -12.1125
sigma_b_investigate_intra -0.17337 NaN NaN 0.25349 -0.6839
b_amnesty 0.81700 0.03617 22.588 0.04578 17.8458
b_violate_low 0.00000 NA NA NA NA
b_violate_med 0.94518 0.07058 13.392 0.07433 12.7152
b_violate_high 1.13297 0.04789 23.655 0.05461 20.7454