I am trying to estimate a LC model with random parameters.
Without the random parameters, the exact same model specification works without a problem.
Unfortunately, I am getting the error:
Code: Select all
Error in inClassProb[[i]]$like * classProb[[i]]$like :
non-conformable arraysWhen using multiple cores, I am getting a different error:
Code: Select all
Error in checkForRemoteErrors(lapply(cl, recvResult)) :
16 nodes produced errors; first error: non-conformable arraysCode: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
apollo_control = list(
modelName = "LC_with_covariates",
modelDescr = "LC model with covariates",
indivID = "RID",
nCores = 1,
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("Apollo Modell/merged_csv.csv",header=TRUE)
### for data dictionary, use ?apollo_swissRouteChoiceData
database <- database[order(database$RID), ]
database$age_group <- cut(database$age,
breaks = c(18, 29, 49, 65, Inf),
labels = c("18-29", "30-49", "50-65", "65+"),
right = FALSE)
# Create dummy variables for each age category
database <- cbind(database,
model.matrix(~ age_group - 1, data = database))
# Rename the dummy variables
names(database)[(ncol(database)-3):ncol(database)] <- c("age_18_29", "age_30_49", "age_50_65", "age_65_plus")
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(log_asc_alt1_a_mu = 2,
log_asc_alt1_b_mu = -5,
log_asc_alt1_c_mu = 1.3,
log_asc_alt1_d_mu = 3.4,
log_asc_alt2_a_mu = 1.5,
log_asc_alt2_b_mu = -5.7,
log_asc_alt2_c_mu = 0.4,
log_asc_alt2_d_mu = 2.36,
log_asc_alt3_a_mu = 1.1,
log_asc_alt3_b_mu = -4.6,
log_asc_alt3_c_mu = 0.88,
log_asc_alt3_d_mu = 3.78,
log_asc_alt1_a_sigma = 0,
log_asc_alt1_b_sigma = 0,
log_asc_alt1_c_sigma = 0,
log_asc_alt1_d_sigma = 0,
log_asc_alt2_a_sigma = 0,
log_asc_alt2_b_sigma = 0,
log_asc_alt2_c_sigma = 0,
log_asc_alt2_d_sigma = 0,
log_asc_alt3_a_sigma = 0,
log_asc_alt3_b_sigma = 0,
log_asc_alt3_c_sigma = 0,
log_asc_alt3_d_sigma = 0,
asc_optout = 0,
beta_cost_a = -0.03,
beta_cost_b = -0.01,
beta_cost_c = -0.18,
beta_cost_d = -0.11,
beta_50m_a = -0.18,
beta_50m_b = 0.51,
beta_50m_c = -1.06,
beta_50m_d = -0.83,
beta_100m_a =-0.04,
beta_100m_b =-0.65,
beta_100m_c =0.06,
beta_100m_d =-0.41,
delta_a = -0.5,
gamma_gender_a = 0.99,
gamma_student_a = -1.31,
gamma_rural_a = 0.39,
gamma_ptsub_a = -0.69,
gamma_18_29_a = 0.47,
gamma_30_49_a = 0.25,
gamma_50_65_a = 0.87,
gamma_65plus_a = 0,
delta_b = 2.51,
gamma_gender_b = -0.14,
gamma_student_b = 1.03,
gamma_rural_b = 0.12,
gamma_ptsub_b = -0.55,
gamma_18_29_b = -3.72,
gamma_30_49_b = -2.53,
gamma_50_65_b = -0.87,
gamma_65plus_b = 0,
delta_c = -1.1,
gamma_gender_c = 0.36,
gamma_student_c = 0.12,
gamma_rural_c = 0.05,
gamma_ptsub_c = 0.13,
gamma_18_29_c = 0.78,
gamma_30_49_c = 0.43,
gamma_50_65_c = 1.04,
gamma_65plus_c = 0,
delta_d = 0,
gamma_gender_d = 0,
gamma_student_d = 0,
gamma_rural_d = 0,
gamma_ptsub_d = 0,
gamma_18_29_d = 0,
gamma_30_49_d = 0,
gamma_50_65_d = 0,
gamma_65plus_d = 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("asc_optout","delta_d","gamma_gender_d", "gamma_student_d", "gamma_rural_d", "gamma_ptsub_d", "gamma_18_29_d", "gamma_30_49_d", "gamma_50_65_d", "gamma_65plus_c", "gamma_65plus_a", "gamma_65plus_b", "gamma_65plus_d")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="halton",
interNDraws=800,
interUnifDraws=c(),
interNormDraws=c("draws_asc_alt1","draws_asc_alt2","draws_asc_alt3"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["asc_alt1_a"]] = -exp(log_asc_alt1_a_mu + log_asc_alt1_a_sigma*draws_asc_alt1)
randcoeff[["asc_alt1_b"]] = -exp(log_asc_alt1_b_mu + log_asc_alt1_b_sigma*draws_asc_alt1)
randcoeff[["asc_alt1_c"]] = -exp(log_asc_alt1_c_mu + log_asc_alt1_c_sigma*draws_asc_alt1)
randcoeff[["asc_alt1_d"]] = -exp(log_asc_alt1_d_mu + log_asc_alt1_d_sigma*draws_asc_alt1)
randcoeff[["asc_alt2_a"]] = -exp(log_asc_alt2_a_mu + log_asc_alt2_a_sigma*draws_asc_alt2)
randcoeff[["asc_alt2_b"]] = -exp(log_asc_alt2_b_mu + log_asc_alt2_b_sigma*draws_asc_alt2)
randcoeff[["asc_alt2_c"]] = -exp(log_asc_alt2_c_mu + log_asc_alt2_c_sigma*draws_asc_alt2)
randcoeff[["asc_alt2_d"]] = -exp(log_asc_alt2_d_mu + log_asc_alt2_d_sigma*draws_asc_alt2)
randcoeff[["asc_alt3_a"]] = -exp(log_asc_alt3_a_mu + log_asc_alt3_a_sigma*draws_asc_alt3)
randcoeff[["asc_alt3_b"]] = -exp(log_asc_alt3_b_mu + log_asc_alt3_b_sigma*draws_asc_alt3)
randcoeff[["asc_alt3_c"]] = -exp(log_asc_alt3_c_mu + log_asc_alt3_c_sigma*draws_asc_alt3)
randcoeff[["asc_alt3_d"]] = -exp(log_asc_alt3_d_mu + log_asc_alt3_d_sigma*draws_asc_alt3)
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["beta_cost"]] = list(beta_cost_a, beta_cost_b, beta_cost_c, beta_cost_d)
lcpars[["asc_alt1"]] = list(asc_alt1_a, asc_alt1_b, asc_alt1_c, asc_alt1_d)
lcpars[["asc_alt2"]] = list(asc_alt2_a, asc_alt2_b, asc_alt2_c, asc_alt2_d)
lcpars[["asc_alt3"]] = list(asc_alt3_a, asc_alt3_b, asc_alt3_c, asc_alt3_d)
lcpars[["beta_50m"]] = list(beta_50m_a, beta_50m_b, beta_50m_c, beta_50m_d)
lcpars[["beta_100m"]] = list(beta_100m_a, beta_100m_b, beta_100m_c, beta_100m_d)
### Utilities of class allocation model
V=list()
V[["class_a"]] = delta_a + gamma_gender_a * gender + gamma_student_a * studentstatus + gamma_rural_a * rural + gamma_ptsub_a * ptsub + gamma_18_29_a * age_18_29 + gamma_30_49_a * age_30_49 + gamma_50_65_a * age_50_65 + gamma_65plus_a * age_65_plus
V[["class_b"]] = delta_b + gamma_gender_b * gender + gamma_student_b * studentstatus + gamma_rural_b * rural + gamma_ptsub_b * ptsub + gamma_18_29_b * age_18_29 + gamma_30_49_b * age_30_49 + gamma_50_65_b * age_50_65 + gamma_65plus_b * age_65_plus
V[["class_c"]] = delta_c + gamma_gender_c * gender + gamma_student_c * studentstatus + gamma_rural_c * rural + gamma_ptsub_c * ptsub + gamma_18_29_c * age_18_29 + gamma_30_49_c * age_30_49 + gamma_50_65_c * age_50_65 + gamma_65plus_c * age_65_plus
V[["class_d"]] = delta_d + gamma_gender_d * gender + gamma_student_d * studentstatus + gamma_rural_d * rural + gamma_ptsub_d * ptsub + gamma_18_29_d * age_18_29 + gamma_30_49_d * age_30_49 + gamma_50_65_d * age_50_65 + gamma_65plus_d * age_65_plus
classAlloc_settings = list(
classes = c(class_a=1, class_b=2, class_c=3, class_d=4),
utilities = V
)
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
return(lcpars)
}
# ################################################################# #
#### 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()
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3, optout=4),
avail = list(alt1=avail_alt1, alt2=avail_alt2, alt3=avail_alt3, optout=1),
choiceVar = pref1
)
### Loop over classes
for(s in 1:4){
### Compute class-specific utilities
V=list()
V[["alt1"]] = asc_alt1[[s]] + beta_cost[[s]]*price_alt1
V[["alt2"]] = asc_alt2[[s]] + beta_cost[[s]]*price_alt2
V[["alt3"]] = asc_alt3[[s]] + beta_cost[[s]]*price_alt3 + beta_50m[[s]]*min_50.0 + beta_100m[[s]]*min_100.0
V[["optout"]] = asc_optout
mnl_settings$utilities = V
mnl_settings$componentName = paste0("Class_",s)
### Compute within-class choice probabilities using MNL model
P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs, functionality)
}
### Compute latent class model probabilities
lc_settings = list(inClassProb = P, classProb=pi_values)
P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
### Average across inter-individual draws in class allocation probabilities
P[["model"]] = apollo_avgInterDraws(P[["model"]], apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
### Estimate model
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
### Show output in screen
apollo_modelOutput(model)
### Save output to file(s)
apollo_saveOutput(model)
Code: Select all
apollo_draws and apollo_randCoeff were found, so apollo_control$mixing was set to TRUE
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.
WARNING: Your database contains some entries that are NA. This may well be intentional, but be advised that if these entries are used in your model, the behaviour
may be unexpected.
WARNING: Your database contains variable(s) "age_group" codified as factors. Apollo does not support factors, and using them inside apollo_probabilities
may lead to NA values in the loglikelihood. If you want to use these variables, we recommend manually transforming them into numeric variables.
All checks on database completed.
Generating inter-individual draws ... Done
>
> # ################################################################# #
> #### 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()
+
+ ### Define settings for MNL model component that are generic across classes
+ mnl_settings = list(
+ alternatives = c(alt1=1, alt2=2, alt3=3, optout=4),
+ avail = list(alt1=avail_alt1, alt2=avail_alt2, alt3=avail_alt3, optout=1),
+ choiceVar = pref1
+ )
+
+ ### Loop over classes
+ for(s in 1:4){
+
+ ### Compute class-specific utilities
+ V=list()
+ V[["alt1"]] = asc_alt1[[s]] + beta_cost[[s]]*price_alt1
+ V[["alt2"]] = asc_alt2[[s]] + beta_cost[[s]]*price_alt2
+ V[["alt3"]] = asc_alt3[[s]] + beta_cost[[s]]*price_alt3 + beta_50m[[s]]*min_50.0 + beta_100m[[s]]*min_100.0
+ V[["optout"]] = asc_optout
+
+ mnl_settings$utilities = V
+ mnl_settings$componentName = paste0("Class_",s)
+
+ ### Compute within-class choice probabilities using MNL model
+ P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
+
+ ### Take product across observation for same individual
+ P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs, functionality)
+
+ }
+
+ ### Compute latent class model probabilities
+ lc_settings = list(inClassProb = P, classProb=pi_values)
+ P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
+
+ ### Average across inter-individual draws in class allocation probabilities
+ P[["model"]] = apollo_avgInterDraws(P[["model"]], apollo_inputs, functionality)
+
+ ### Prepare and return outputs of function
+ P = apollo_prepareProb(P, apollo_inputs, functionality)
+ return(P)
+ }
>
>
> # ################################################################# #
> #### MODEL ESTIMATION ####
> # ################################################################# #
>
> ### Estimate model
> model = apollo_estimate(apollo_beta, apollo_fixed,
+ apollo_probabilities, apollo_inputs)
INFORMATION: You can use multiple processor cores to speed up estimation. To do so, specify the desired number of cores using the setting nCores inside
"apollo_control", for example "nCores=2". This computer has 20 available cores. We recommend using no more than 19 cores.
Preparing user-defined functions.
Testing likelihood function...
Apollo found a model component of type classAlloc without a componentName. The name was set to "classAlloc" by default.
INFORMATION: Setting "avail" is missing, so full availability is assumed.
Overview of choices for MNL model component Class_1:
alt1 alt2 alt3 optout
Times available 3801.00 3814.00 3801.00 5708.00
Times chosen 625.00 927.00 825.00 3331.00
Percentage chosen overall 10.95 16.24 14.45 58.36
Percentage chosen when available 16.44 24.31 21.70 58.36
Overview of choices for MNL model component Class_2:
alt1 alt2 alt3 optout
Times available 3801.00 3814.00 3801.00 5708.00
Times chosen 625.00 927.00 825.00 3331.00
Percentage chosen overall 10.95 16.24 14.45 58.36
Percentage chosen when available 16.44 24.31 21.70 58.36
Overview of choices for MNL model component Class_3:
alt1 alt2 alt3 optout
Times available 3801.00 3814.00 3801.00 5708.00
Times chosen 625.00 927.00 825.00 3331.00
Percentage chosen overall 10.95 16.24 14.45 58.36
Percentage chosen when available 16.44 24.31 21.70 58.36
Overview of choices for MNL model component Class_4:
alt1 alt2 alt3 optout
Times available 3801.00 3814.00 3801.00 5708.00
Times chosen 625.00 927.00 825.00 3331.00
Percentage chosen overall 10.95 16.24 14.45 58.36
Percentage chosen when available 16.44 24.31 21.70 58.36
Summary of class allocation for model component :
Mean prob.
Class_1 0.2222
Class_2 0.4117
Class_3 0.1681
Class_4 0.1980
The class allocation probabilities for model component "model" are calculated at the observation level in 'apollo_lcPars', but are used in
'apollo_probabilities' to multiply within class probabilities that are at the individual level. Apollo will average the class allocation probabilities
across observations for the same individual level before using them to multiply the within-class probabilities. If your class allocation probabilities
are constant across choice situations for the same individual, then this is of no concern. If your class allocation probabilities however vary across
choice tasks, then you should change your model specification in 'apollo_probabilities' to only call 'apollo_panelProd' after calling 'apollo_lc'.
Pre-processing likelihood function...
Testing influence of parameters
Error in inClassProb[[i]]$like * classProb[[i]]$like :
non-conformable arraysI have specified this availability and in the overview the correct times available and times chosen are displayed.
Thank you very much.