I am estimating a mixed logit model with the latest version of Apollo (0.2.5) and 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...
Error in apollo_inputs$apollo_control$mixing && is.function(apollo_inputs$apollo_randCoeff) :
invalid 'x' type in 'x && y'
==========
Because of this, I cannot obtain the output files.
What is hard to understand is that this error occurs only when I set "interNDraws" to 1000 and estimate the model with the full dataset (19860 observations, 4700 individuals).
I didn't encounter the error when I reduced the number of draws (like 10, 100, 250) or the sample size.
What do you think possible causes?
Kind regards,
---
The code is as follows:
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
library(tidyverse)
library(magrittr)
### Initialise code
apollo_initialise()
modelname = "test"
### Set core controls
apollo_control = list(
modelName = modelname,
modelDescr ="MXL model",
indivID ="no",
seed = 123,
mixing = TRUE,
nCores = parallel::detectCores() - 1
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("dataset.csv", header=TRUE)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(ASC_1 =0,
ASC_2 =0,
ASC_3_1 =0,
ASC_3_2 =0,
b_1 =0,
b_2 =0,
b_3 =0,
b_4 =0,
b_5 =0,
b_6 =0,
b_7 =0,
b_8 =0,
b_9 =-10,
b_10 =-10,
b_11 =0,
b_12 =0,
b_13 =0,
b_14 =0,
b_15 =0,
b_16 =0,
b_17 =0,
b_18 =0,
b_19 =0,
b_20 =0,
b_21 =0,
b_22 =0,
b_23 =0,
b_24 =0,
b_25 =0,
b_26 =0,
b_27 =0,
b_28 =0,
b_29 =0,
b_30 =0,
b_31 =0,
b_32 =0,
b_33 =0,
b_34 =0,
b_35 =0,
b_36 =0,
b_37 =0,
b_38 =0,
b_39 =0,
b_40 =0,
b_41 =0,
v_time_mean =0,
v_time_sd =0,
v_time_sd2 =0,
v_time_sd3 =0,
v_time2_mean =0,
v_time2_sd =0,
v_time2_sd2 =0,
v_time2_sd3 =0,
xi_1_sd = 0.1,
xi_2_sd = 0.1,
xi_3_sd = 0.1,
xi_3_1_sd = 0.1,
xi_3_2_sd = 0.1,
xi_3_3_sd = 0.1,
xi_2_1_sd = 0.1,
gamma_1 = 0,
gamma_2 = 0
)
apollo_draws = list(
## inter-personal heterogeneity
interDrawsType = "sobolOwenFaureTezuka",
interNDraws = 1000,
interUnifDraws = c("draws_b_cost"),
interNormDraws = c("draws_v_time", "draws_v_time2",
"draws_xi_1", "draws_xi_2", "draws_xi_3",
"draws_xi_3_1", "draws_xi_3_2", "draws_xi_3_3", "draws_xi_2_1"),
## intra-personal heterogeneity
intraDrawsType = "mlhs",
intraNDraws = 100,
intraUnifDraws =c(),
intraNormDraws = c()
)
### Defining a function that returns random parameters
apollo_randCoeff = function(apollo_beta , apollo_inputs){
randcoeff = list ()
randcoeff[["v_time"]] = (v_time_mean + v_time_sd * draws_v_time +
v_time_sd2 * draws_v_time^2 + v_time_sd3 * draws_v_time^3) * (ref_frequency)^gamma_1
randcoeff[["v_time2"]] = v_time2_mean +
v_time2_sd * draws_v_time2 +
v_time2_sd2 * draws_v_time2^2 +
v_time2_sd3 * draws_v_time2^3
randcoeff[["b_cost"]] = -exp(b_9 + (b_10-b_9) * draws_b_cost) * (ref_income)^gamma_2
randcoeff[["xi_1"]] = xi_1_sd * draws_xi_1
randcoeff[["xi_2"]] = xi_2_sd * draws_xi_2
randcoeff[["xi_3"]] = xi_3_sd * draws_xi_3
randcoeff[["xi_3_1"]] = xi_3_1_sd * draws_xi_3_1
randcoeff[["xi_3_2"]] = xi_3_2_sd * draws_xi_3_2
randcoeff[["xi_3_3"]] = xi_3_3_sd * draws_xi_3_3
randcoeff[["xi_2_1"]] = xi_2_1_sd * draws_xi_2_1
return(randcoeff)
}
### 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_17", "b_18")
### Read in starting values for at least some parameters from existing model output file
apollo_beta=apollo_readBeta(apollo_beta,
apollo_fixed,
modelname,
overwriteFixed=FALSE)
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
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()
# ### Create alternative specific constants and coefficients using interactions with socio-demographics
V_3 = ASC_3_2 + xi_3 +
b_cost * v_time2 * time2 +
b_2 * youth +
b_4 * elderly +
b_25 * male +
b_41 * capital +
b_8 * s1 +
b_18 * freq +
b_27 * free +
b_39 * teleworker
V_3_2 = V_3 + b_cost * v_time * time2 + b_cost * (cost_3 + add_cost * holiday2) + (b_19 + xi_2_1) * holiday2
V_3_3 = V_3 + b_cost * v_time * time3 + b_cost * (cost_3 + add_cost * holiday3) + (b_19 + xi_2_1) * holiday3
V_3_4 = V_3 + b_cost * v_time * time4 + b_cost * (cost_3 + add_cost * holiday4) + (b_19 + xi_2_1) * holiday4
V_3_5 = V_3 + b_cost * v_time * time5 + b_cost * (cost_3 + add_cost * holiday5) + (b_19 + xi_2_1) * holiday5
V_3_6 = V_3 + b_cost * v_time * time6 + b_cost * (cost_3 + add_cost * holiday6) + (b_19 + xi_2_1) * holiday6
V_3_7 = V_3 + b_cost * v_time * time7 + b_cost * (cost_3 + add_cost * holiday7) + (b_19 + xi_2_1) * holiday7
V_3_8 = V_3 + b_cost * v_time * time8 + b_cost * (cost_3 + add_cost * holiday8) + (b_19 + xi_2_1) * holiday8
V_m = ASC_3_1 + xi_3_1 +
b_5 * alone +
b_11 * live_child +
b_13 * couple_double +
b_15 * live_family +
b_22 * house +
b_32 * prof +
b_34 * public +
b_36 * self +
b_30 * parttime +
b_28 * nojob
V_n = ASC_3_2 + xi_3_2 +
b_cost * add_cost_n +
b_6 * alone +
b_12 * live_child +
b_14 * couple_double +
b_16 * live_family +
b_23 * house +
b_33 * prof +
b_35 * public +
b_37 * self +
b_31 * parttime +
b_29 * nojob
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[['V_1']] = ASC_1 + xi_1 +
b_cost * v_time * time1 +
b_cost * cost1 +
b_1 * youth +
b_3 * elderly +
b_24 * male +
b_40 * capital +
b_7 * s1 +
b_17 * freq +
b_26 * free +
b_38 * teleworker
V[['V_2']] = b_cost * v_time * time2 + b_cost * cost2 + xi_2
V[['V_3_2_1']]=V_3_2 + V_m + b_20 * holiday2
V[['V_3_2_2']]=V_3_2 + xi_3_3
V[['V_3_2_3']]=V_3_2 + V_n + b_21 * holiday2
V[['V_3_3_1']]=V_3_3 + V_m + b_20 * holiday3
V[['V_3_3_2']]=V_3_3 + xi_3_3
V[['V_3_3_3']]=V_3_3 + V_n + b_21 * holiday3
V[['V_3_4_1']]=V_3_4 + V_m + b_20 * holiday4
V[['V_3_4_2']]=V_3_4 + xi_3_3
V[['V_3_4_3']]=V_3_4 + V_n + b_21 * holiday4
V[['V_3_5_1']]=V_3_5 + V_m + b_20 * holiday5
V[['V_3_5_2']]=V_3_5 + xi_3_3
V[['V_3_5_3']]=V_3_5 + V_n + b_21 * holiday5
V[['V_3_6_1']]=V_3_6 + V_m + b_20 * holiday6
V[['V_3_6_2']]=V_3_6 + xi_3_3
V[['V_3_6_3']]=V_3_6 + V_n + b_21 * holiday6
V[['V_3_7_1']]=V_3_7 + V_m + b_20 * holiday7
V[['V_3_7_2']]=V_3_7 + xi_3_3
V[['V_3_7_3']]=V_3_7 + V_n + b_21 * holiday7
V[['V_3_8_1']]=V_3_8 + V_m + b_20 * holiday8
V[['V_3_8_2']]=V_3_8 + xi_3_3
V[['V_3_8_3']]=V_3_8 + V_n + b_21 * holiday8
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(V_1 = 1000,
V_2 = 2000,
V_3_2_1 = 3021,
V_3_2_2 = 3022,
V_3_2_3 = 3023,
V_3_3_1 = 3031,
V_3_3_2 = 3032,
V_3_3_3 = 3033,
V_3_4_1 = 3041,
V_3_4_2 = 3042,
V_3_4_3 = 3043,
V_3_5_1 = 3051,
V_3_5_2 = 3052,
V_3_5_3 = 3053,
V_3_6_1 = 3061,
V_3_6_2 = 3062,
V_3_6_3 = 3063,
V_3_7_1 = 3071,
V_3_7_2 = 3072,
V_3_7_3 = 3073,
V_3_8_1 = 3081,
V_3_8_2 = 3082,
V_3_8_3 = 3083),
avail = list(V_1 = 1,
V_2 = 1,
V_3_2_1 = 1,
V_3_2_2 = 1,
V_3_2_3 = 1,
V_3_3_1 = 1,
V_3_3_2 = 1,
V_3_3_3 = 1,
V_3_4_1 = 1,
V_3_4_2 = 1,
V_3_4_3 = 1,
V_3_5_1 = 1,
V_3_5_2 = 1,
V_3_5_3 = 1,
V_3_6_1 = 1,
V_3_6_2 = 1,
V_3_6_3 = 1,
V_3_7_1 = 1,
V_3_7_2 = 1,
V_3_7_3 = 1,
V_3_8_1 = 1,
V_3_8_2 = 1,
V_3_8_3 = 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 ####
# ################################################################# #
# apollo_beta = apollo_searchStart(apollo_beta,
# apollo_fixed,
# apollo_probabilities,
# apollo_inputs)#,
# #searchStart_settings)
model = apollo_estimate(apollo_beta,
apollo_fixed,
apollo_probabilities,
apollo_inputs,
estimate_settings=list(hessianRoutine="maxLik",maxIterations=300)
)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
apollo_saveOutput(model)