error in estimating a latent class model
Posted: 11 Jan 2022, 20:34
Hi, I am estimating a latent class model. In estimation the following error message poped up:
"Testing influence of parameters
Error in apollo_avgInterDraws(P[["model"]], apollo_inputs, functionality) :
No Inter-individuals draws to average over!"
I appreciate your suggestions!
I attached the code here:
##----------------------------------
## estimation using Apollo
##----------------------------------
# /mnt/rdrive/xdu23
# cd /mnt/rdrive/xdu23/Corn\ BE\ Yield\ Impact\
##following the example here: http://apollochoicemodelling.com/files/ ... ce_space.r
rm(list=ls())
setwd("Y:/Corn BE Yield Impact")
#install.packages("apollo")
library(apollo)
load("sample_logit_10000_forced_final.RData")
#head(data_sample)
#here "id" is the choice situation ("unit_id" + "crop_year"); "ids" is unit_id
colnames(data_sample)[which(colnames(data_sample)=="ids")] <- c("unit_id")
temp <- which(data_sample$unit_id == 706193 | data_sample$unit_id == 338891) #for this data points, eva=Inf
data_sample <- data_sample[-temp,]
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "LC_MMNL",
modelDescr = "Latent class with continuous random parameters",
indivID = "unit_id",
mixing = TRUE,
nCores = 1,
outputDirectory = "output",
panelData = TRUE
)
database = data_sample
head(database)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation; asc-alternative specific constant.
apollo_beta=c(asc_pln1 = 0, asc_pln2 = 1.13, asc_pln3 = 4.83, asc_pln4 = 9.93, asc_pln5 = 10.20,
asc_pln6 = 10.18, asc_pln7 = 6.07, asc_pln8 = 1.27, asc_pln9 = 11.25, asc_pln10 = 13.24,
asc_pln11 = 12.74, asc_pln12 = 9.42, asc_pln13 = 4.78, asc_pln14 = -0.46, asc_pln15 = 2.54,
asc_pln16 = 6.89, asc_pln17 = 10.76, asc_pln18 = 12.64, asc_pln19 = 12.86, asc_pln20 = 10.64, asc_pln21 = 5.88,
# b_prm = 0,
# b_sub = 0,
# b_lia = 0,
asc_pln2_shift_acres = -0.00007, asc_pln3_shift_acres = 0.0003, asc_pln4_shift_acres = -0.0015, asc_pln5_shift_acres = -0.00078, asc_pln6_shift_acres = -0.0018,
asc_pln7_shift_acres = -0.00077, asc_pln8_shift_acres = -0.0016, asc_pln9_shift_acres = 0.0016, asc_pln10_shift_acres = 0.0018, asc_pln11_shift_acres = 0.0014,
asc_pln12_shift_acres = 0.0015, asc_pln13_shift_acres = 0.006, asc_pln14_shift_acres = 0.0019, asc_pln15_shift_acres = -0.001, asc_pln16_shift_acres = 0.0016,
asc_pln17_shift_acres = 0.005, asc_pln18_shift_acres = 0.0009, asc_pln19_shift_acres = 0.0009, asc_pln20_shift_acres = 0.0014, asc_pln21_shift_acres = 0.0012,
asc_pln2_shift_myield = -0.01, asc_pln3_shift_myield = -0.015, asc_pln4_shift_myield = -0.002, asc_pln5_shift_myield = 0.002, asc_pln6_shift_myield = 0.015,
asc_pln7_shift_myield = 0.023, asc_pln8_shift_myield = 0.023, asc_pln9_shift_myield = -0.001, asc_pln10_shift_myield = -0.0006,asc_pln11_shift_myield = 0.007,
asc_pln12_shift_myield = 0.013, asc_pln13_shift_myield = 0.016, asc_pln14_shift_myield = -0.002,asc_pln15_shift_myield = -0.011, asc_pln16_shift_myield = -0.012,
asc_pln17_shift_myield = -0.001, asc_pln18_shift_myield = 0.00001,asc_pln19_shift_myield = 0.013, asc_pln20_shift_myield = 0.032, asc_pln21_shift_myield = 0.045,
asc_pln2_shift_vyield = -0.028, asc_pln3_shift_vyield = 0.01, asc_pln4_shift_vyield = 0.0055, asc_pln5_shift_vyield = 0.016, asc_pln6_shift_vyield = 0.007,
asc_pln7_shift_vyield = 0.037, asc_pln8_shift_vyield = 0.021, asc_pln9_shift_vyield = 0.012, asc_pln10_shift_vyield = 0.013, asc_pln11_shift_vyield = 0.019,
asc_pln12_shift_vyield = 0.019, asc_pln13_shift_vyield = 0.03, asc_pln14_shift_vyield = 0.009, asc_pln15_shift_vyield = 0.026, asc_pln16_shift_vyield = 0.018,
asc_pln17_shift_vyield = 0.013, asc_pln18_shift_vyield = 0.017, asc_pln19_shift_vyield = 0.018, asc_pln20_shift_vyield = 0.024, asc_pln21_shift_vyield = 0.03,
delta_a = 0, delta_b = 0, gamma_atten_a = 0, gamma_atten_b = 0, gamma_d03_a = 0, gamma_d03_b = 0, gamma_d09_a = 0, gamma_d09_b = 0, gamma_acres_a = 0, gamma_acres_b = 0,
gamma_myield_a = 0, gamma_myield_b = 0, gamma_vyield_a = 0, gamma_vyield_b = 0, gamma_disas_a = 0, gamma_disas_b = 0,
mu_log_b_prm_a = -2.5, sigma_log_b_prm_a = 0.23, mu_log_b_prm_b = -2.5, sigma_log_b_prm_b = 0.23,
mu_log_b_sub_a = -2.5, sigma_log_b_sub_a = 0.25, mu_log_b_sub_b = -2.5, sigma_log_b_sub_b = 0.25,
mu_log_b_lia_a = -5, sigma_log_b_lia_a = 0.10, mu_log_b_lia_b = -2.5, sigma_log_b_lia_b = 0.10)
### 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_pln1","delta_b", "gamma_atten_b")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("draws_prm","draws_sub","draws_lia"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_prm_a"]] = -exp( mu_log_b_prm_a + sigma_log_b_prm_a * draws_prm )
randcoeff[["b_prm_b"]] = -exp( mu_log_b_prm_b + sigma_log_b_prm_b * draws_prm )
randcoeff[["b_sub_a"]] = exp( mu_log_b_sub_a + sigma_log_b_sub_a * draws_sub )
randcoeff[["b_sub_b"]] = exp( mu_log_b_sub_b + sigma_log_b_sub_b * draws_sub )
randcoeff[["b_lia_a"]] = exp( mu_log_b_lia_a + sigma_log_b_lia_a * draws_lia )
randcoeff[["b_lia_b"]] = exp( mu_log_b_lia_b + sigma_log_b_lia_b * draws_lia )
# randcoeff[["delta_a"]] = delta_a_mu + delta_a_sig * draws_pi
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_prm"]] = list(b_prm_a, b_prm_b)
lcpars[["b_sub"]] = list(b_sub_a, b_sub_b)
lcpars[["b_lia"]] = list(b_lia_a, b_lia_b)
V=list()
V[["class_a"]] = delta_a + gamma_atten_a*atten
V[["class_b"]] = delta_b + gamma_atten_b*atten
#+ gamma_d03_a * D03 + gamma_d09_a * D09 + gamma_disas_a * disas + gamma_acres_a * acres + gamma_myield_a * yield_mean + gamma_vyield_a * yield_var
#+ gamma_d03_b * D03 + gamma_d09_b * D09 + gamma_disas_b * disas + gamma_acres_b * acres + gamma_myield_b * yield_mean + gamma_vyield_b * yield_var
classAlloc_settings = list(
alternatives = c(class_a=1, class_b=2),
avail = 1,
V = 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(pln1=1, pln2=2, pln3=3, pln4=4, pln5=5, pln6=6, pln7=7, pln8=8, pln9=9, pln10=10, pln11=11, pln12=12,pln13=13, pln14=14, pln15=15,pln16=16, pln17=17, pln18=18,pln19=19, pln20=20, pln21=21),
avail = list(pln1=av1, pln2=av2, pln3=av3, pln4=av4, pln5=av5, pln6=av6, pln7=av7, pln8=av8,pln9=av9, pln10=av10, pln11=av11, pln12=av12,pln13=av13,pln14=av14,pln15=av15,pln16=av16,pln17=av17,pln18=av18,pln19=av19,pln20=av20,pln21=av21),
choiceVar = choice
)
### Create alternative specific constants and coefficients using interactions with socio-demographics
asc_pln2_value = asc_pln2 + asc_pln2_shift_acres * acres + asc_pln2_shift_myield * yield_mean + asc_pln2_shift_vyield * yield_var
asc_pln3_value = asc_pln3 + asc_pln3_shift_acres * acres + asc_pln3_shift_myield * yield_mean + asc_pln3_shift_vyield * yield_var
asc_pln4_value = asc_pln4 + asc_pln4_shift_acres * acres + asc_pln4_shift_myield * yield_mean + asc_pln4_shift_vyield * yield_var
asc_pln5_value = asc_pln5 + asc_pln5_shift_acres * acres + asc_pln5_shift_myield * yield_mean + asc_pln5_shift_vyield * yield_var
asc_pln6_value = asc_pln6 + asc_pln6_shift_acres * acres + asc_pln6_shift_myield * yield_mean + asc_pln6_shift_vyield * yield_var
asc_pln7_value = asc_pln7 + asc_pln7_shift_acres * acres + asc_pln7_shift_myield * yield_mean + asc_pln7_shift_vyield * yield_var
asc_pln8_value = asc_pln8 + asc_pln8_shift_acres * acres + asc_pln8_shift_myield * yield_mean + asc_pln8_shift_vyield * yield_var
asc_pln9_value = asc_pln9 + asc_pln9_shift_acres * acres + asc_pln9_shift_myield * yield_mean + asc_pln9_shift_vyield * yield_var
asc_pln10_value = asc_pln10 + asc_pln10_shift_acres * acres + asc_pln10_shift_myield * yield_mean + asc_pln10_shift_vyield * yield_var
asc_pln11_value = asc_pln11 + asc_pln11_shift_acres * acres + asc_pln11_shift_myield * yield_mean + asc_pln11_shift_vyield * yield_var
asc_pln12_value = asc_pln12 + asc_pln12_shift_acres * acres + asc_pln12_shift_myield * yield_mean + asc_pln12_shift_vyield * yield_var
asc_pln13_value = asc_pln13 + asc_pln13_shift_acres * acres + asc_pln13_shift_myield * yield_mean + asc_pln13_shift_vyield * yield_var
asc_pln14_value = asc_pln14 + asc_pln14_shift_acres * acres + asc_pln14_shift_myield * yield_mean + asc_pln14_shift_vyield * yield_var
asc_pln15_value = asc_pln15 + asc_pln15_shift_acres * acres + asc_pln15_shift_myield * yield_mean + asc_pln15_shift_vyield * yield_var
asc_pln16_value = asc_pln16 + asc_pln16_shift_acres * acres + asc_pln16_shift_myield * yield_mean + asc_pln16_shift_vyield * yield_var
asc_pln17_value = asc_pln17 + asc_pln17_shift_acres * acres + asc_pln17_shift_myield * yield_mean + asc_pln17_shift_vyield * yield_var
asc_pln18_value = asc_pln18 + asc_pln18_shift_acres * acres + asc_pln18_shift_myield * yield_mean + asc_pln18_shift_vyield * yield_var
asc_pln19_value = asc_pln19 + asc_pln19_shift_acres * acres + asc_pln19_shift_myield * yield_mean + asc_pln19_shift_vyield * yield_var
asc_pln20_value = asc_pln20 + asc_pln20_shift_acres * acres + asc_pln20_shift_myield * yield_mean + asc_pln20_shift_vyield * yield_var
asc_pln21_value = asc_pln21 + asc_pln21_shift_acres * acres + asc_pln21_shift_myield * yield_mean + asc_pln21_shift_vyield * yield_var
### Loop over classes
for(s in 1:2){
### Compute class-specific utilities
V=list()
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V[["pln1"]] = asc_pln1 + b_prm[[s]] * prm1 + b_sub[[s]] * sub1 + b_lia[[s]] * lia1
V[["pln2"]] = asc_pln2_value + b_prm[[s]] * prm2 + b_sub[[s]] * sub2 + b_lia[[s]] * lia2
V[["pln3"]] = asc_pln3_value + b_prm[[s]] * prm3 + b_sub[[s]] * sub3 + b_lia[[s]] * lia3
V[["pln4"]] = asc_pln4_value + b_prm[[s]] * prm4 + b_sub[[s]] * sub4 + b_lia[[s]] * lia4
V[["pln5"]] = asc_pln5_value + b_prm[[s]] * prm5 + b_sub[[s]] * sub5 + b_lia[[s]] * lia5
V[["pln6"]] = asc_pln6_value + b_prm[[s]] * prm6 + b_sub[[s]] * sub6 + b_lia[[s]] * lia6
V[["pln7"]] = asc_pln7_value + b_prm[[s]] * prm7 + b_sub[[s]] * sub7 + b_lia[[s]] * lia7
V[["pln8"]] = asc_pln8_value + b_prm[[s]] * prm8 + b_sub[[s]] * sub8 + b_lia[[s]] * lia8
V[["pln9"]] = asc_pln9_value + b_prm[[s]] * prm9 + b_sub[[s]] * sub9 + b_lia[[s]] * lia9
V[["pln10"]] = asc_pln10_value + b_prm[[s]] * prm10 + b_sub[[s]] * sub10 + b_lia[[s]] * lia10
V[["pln11"]] = asc_pln11_value + b_prm[[s]] * prm11 + b_sub[[s]] * sub11 + b_lia[[s]] * lia11
V[["pln12"]] = asc_pln12_value + b_prm[[s]] * prm12 + b_sub[[s]] * sub12 + b_lia[[s]] * lia12
V[["pln13"]] = asc_pln13_value + b_prm[[s]] * prm13 + b_sub[[s]] * sub13 + b_lia[[s]] * lia13
V[["pln14"]] = asc_pln14_value + b_prm[[s]] * prm14 + b_sub[[s]] * sub14 + b_lia[[s]] * lia14
V[["pln15"]] = asc_pln15_value + b_prm[[s]] * prm15 + b_sub[[s]] * sub15 + b_lia[[s]] * lia15
V[["pln16"]] = asc_pln16_value + b_prm[[s]] * prm16 + b_sub[[s]] * sub16 + b_lia[[s]] * lia16
V[["pln17"]] = asc_pln17_value + b_prm[[s]] * prm17 + b_sub[[s]] * sub17 + b_lia[[s]] * lia17
V[["pln18"]] = asc_pln18_value + b_prm[[s]] * prm18 + b_sub[[s]] * sub18 + b_lia[[s]] * lia18
V[["pln19"]] = asc_pln19_value + b_prm[[s]] * prm19 + b_sub[[s]] * sub19 + b_lia[[s]] * lia19
V[["pln20"]] = asc_pln20_value + b_prm[[s]] * prm20 + b_sub[[s]] * sub20 + b_lia[[s]] * lia20
V[["pln21"]] = asc_pln21_value + b_prm[[s]] * prm21 + b_sub[[s]] * sub21 + b_lia[[s]] * lia21
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)
### Average across inter-individual draws within classes
P[[paste0("Class_",s)]] = apollo_avgInterDraws(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 ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
apollo_saveOutput(model)
"Testing influence of parameters
Error in apollo_avgInterDraws(P[["model"]], apollo_inputs, functionality) :
No Inter-individuals draws to average over!"
I appreciate your suggestions!
I attached the code here:
##----------------------------------
## estimation using Apollo
##----------------------------------
# /mnt/rdrive/xdu23
# cd /mnt/rdrive/xdu23/Corn\ BE\ Yield\ Impact\
##following the example here: http://apollochoicemodelling.com/files/ ... ce_space.r
rm(list=ls())
setwd("Y:/Corn BE Yield Impact")
#install.packages("apollo")
library(apollo)
load("sample_logit_10000_forced_final.RData")
#head(data_sample)
#here "id" is the choice situation ("unit_id" + "crop_year"); "ids" is unit_id
colnames(data_sample)[which(colnames(data_sample)=="ids")] <- c("unit_id")
temp <- which(data_sample$unit_id == 706193 | data_sample$unit_id == 338891) #for this data points, eva=Inf
data_sample <- data_sample[-temp,]
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "LC_MMNL",
modelDescr = "Latent class with continuous random parameters",
indivID = "unit_id",
mixing = TRUE,
nCores = 1,
outputDirectory = "output",
panelData = TRUE
)
database = data_sample
head(database)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation; asc-alternative specific constant.
apollo_beta=c(asc_pln1 = 0, asc_pln2 = 1.13, asc_pln3 = 4.83, asc_pln4 = 9.93, asc_pln5 = 10.20,
asc_pln6 = 10.18, asc_pln7 = 6.07, asc_pln8 = 1.27, asc_pln9 = 11.25, asc_pln10 = 13.24,
asc_pln11 = 12.74, asc_pln12 = 9.42, asc_pln13 = 4.78, asc_pln14 = -0.46, asc_pln15 = 2.54,
asc_pln16 = 6.89, asc_pln17 = 10.76, asc_pln18 = 12.64, asc_pln19 = 12.86, asc_pln20 = 10.64, asc_pln21 = 5.88,
# b_prm = 0,
# b_sub = 0,
# b_lia = 0,
asc_pln2_shift_acres = -0.00007, asc_pln3_shift_acres = 0.0003, asc_pln4_shift_acres = -0.0015, asc_pln5_shift_acres = -0.00078, asc_pln6_shift_acres = -0.0018,
asc_pln7_shift_acres = -0.00077, asc_pln8_shift_acres = -0.0016, asc_pln9_shift_acres = 0.0016, asc_pln10_shift_acres = 0.0018, asc_pln11_shift_acres = 0.0014,
asc_pln12_shift_acres = 0.0015, asc_pln13_shift_acres = 0.006, asc_pln14_shift_acres = 0.0019, asc_pln15_shift_acres = -0.001, asc_pln16_shift_acres = 0.0016,
asc_pln17_shift_acres = 0.005, asc_pln18_shift_acres = 0.0009, asc_pln19_shift_acres = 0.0009, asc_pln20_shift_acres = 0.0014, asc_pln21_shift_acres = 0.0012,
asc_pln2_shift_myield = -0.01, asc_pln3_shift_myield = -0.015, asc_pln4_shift_myield = -0.002, asc_pln5_shift_myield = 0.002, asc_pln6_shift_myield = 0.015,
asc_pln7_shift_myield = 0.023, asc_pln8_shift_myield = 0.023, asc_pln9_shift_myield = -0.001, asc_pln10_shift_myield = -0.0006,asc_pln11_shift_myield = 0.007,
asc_pln12_shift_myield = 0.013, asc_pln13_shift_myield = 0.016, asc_pln14_shift_myield = -0.002,asc_pln15_shift_myield = -0.011, asc_pln16_shift_myield = -0.012,
asc_pln17_shift_myield = -0.001, asc_pln18_shift_myield = 0.00001,asc_pln19_shift_myield = 0.013, asc_pln20_shift_myield = 0.032, asc_pln21_shift_myield = 0.045,
asc_pln2_shift_vyield = -0.028, asc_pln3_shift_vyield = 0.01, asc_pln4_shift_vyield = 0.0055, asc_pln5_shift_vyield = 0.016, asc_pln6_shift_vyield = 0.007,
asc_pln7_shift_vyield = 0.037, asc_pln8_shift_vyield = 0.021, asc_pln9_shift_vyield = 0.012, asc_pln10_shift_vyield = 0.013, asc_pln11_shift_vyield = 0.019,
asc_pln12_shift_vyield = 0.019, asc_pln13_shift_vyield = 0.03, asc_pln14_shift_vyield = 0.009, asc_pln15_shift_vyield = 0.026, asc_pln16_shift_vyield = 0.018,
asc_pln17_shift_vyield = 0.013, asc_pln18_shift_vyield = 0.017, asc_pln19_shift_vyield = 0.018, asc_pln20_shift_vyield = 0.024, asc_pln21_shift_vyield = 0.03,
delta_a = 0, delta_b = 0, gamma_atten_a = 0, gamma_atten_b = 0, gamma_d03_a = 0, gamma_d03_b = 0, gamma_d09_a = 0, gamma_d09_b = 0, gamma_acres_a = 0, gamma_acres_b = 0,
gamma_myield_a = 0, gamma_myield_b = 0, gamma_vyield_a = 0, gamma_vyield_b = 0, gamma_disas_a = 0, gamma_disas_b = 0,
mu_log_b_prm_a = -2.5, sigma_log_b_prm_a = 0.23, mu_log_b_prm_b = -2.5, sigma_log_b_prm_b = 0.23,
mu_log_b_sub_a = -2.5, sigma_log_b_sub_a = 0.25, mu_log_b_sub_b = -2.5, sigma_log_b_sub_b = 0.25,
mu_log_b_lia_a = -5, sigma_log_b_lia_a = 0.10, mu_log_b_lia_b = -2.5, sigma_log_b_lia_b = 0.10)
### 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_pln1","delta_b", "gamma_atten_b")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("draws_prm","draws_sub","draws_lia"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_prm_a"]] = -exp( mu_log_b_prm_a + sigma_log_b_prm_a * draws_prm )
randcoeff[["b_prm_b"]] = -exp( mu_log_b_prm_b + sigma_log_b_prm_b * draws_prm )
randcoeff[["b_sub_a"]] = exp( mu_log_b_sub_a + sigma_log_b_sub_a * draws_sub )
randcoeff[["b_sub_b"]] = exp( mu_log_b_sub_b + sigma_log_b_sub_b * draws_sub )
randcoeff[["b_lia_a"]] = exp( mu_log_b_lia_a + sigma_log_b_lia_a * draws_lia )
randcoeff[["b_lia_b"]] = exp( mu_log_b_lia_b + sigma_log_b_lia_b * draws_lia )
# randcoeff[["delta_a"]] = delta_a_mu + delta_a_sig * draws_pi
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_prm"]] = list(b_prm_a, b_prm_b)
lcpars[["b_sub"]] = list(b_sub_a, b_sub_b)
lcpars[["b_lia"]] = list(b_lia_a, b_lia_b)
V=list()
V[["class_a"]] = delta_a + gamma_atten_a*atten
V[["class_b"]] = delta_b + gamma_atten_b*atten
#+ gamma_d03_a * D03 + gamma_d09_a * D09 + gamma_disas_a * disas + gamma_acres_a * acres + gamma_myield_a * yield_mean + gamma_vyield_a * yield_var
#+ gamma_d03_b * D03 + gamma_d09_b * D09 + gamma_disas_b * disas + gamma_acres_b * acres + gamma_myield_b * yield_mean + gamma_vyield_b * yield_var
classAlloc_settings = list(
alternatives = c(class_a=1, class_b=2),
avail = 1,
V = 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(pln1=1, pln2=2, pln3=3, pln4=4, pln5=5, pln6=6, pln7=7, pln8=8, pln9=9, pln10=10, pln11=11, pln12=12,pln13=13, pln14=14, pln15=15,pln16=16, pln17=17, pln18=18,pln19=19, pln20=20, pln21=21),
avail = list(pln1=av1, pln2=av2, pln3=av3, pln4=av4, pln5=av5, pln6=av6, pln7=av7, pln8=av8,pln9=av9, pln10=av10, pln11=av11, pln12=av12,pln13=av13,pln14=av14,pln15=av15,pln16=av16,pln17=av17,pln18=av18,pln19=av19,pln20=av20,pln21=av21),
choiceVar = choice
)
### Create alternative specific constants and coefficients using interactions with socio-demographics
asc_pln2_value = asc_pln2 + asc_pln2_shift_acres * acres + asc_pln2_shift_myield * yield_mean + asc_pln2_shift_vyield * yield_var
asc_pln3_value = asc_pln3 + asc_pln3_shift_acres * acres + asc_pln3_shift_myield * yield_mean + asc_pln3_shift_vyield * yield_var
asc_pln4_value = asc_pln4 + asc_pln4_shift_acres * acres + asc_pln4_shift_myield * yield_mean + asc_pln4_shift_vyield * yield_var
asc_pln5_value = asc_pln5 + asc_pln5_shift_acres * acres + asc_pln5_shift_myield * yield_mean + asc_pln5_shift_vyield * yield_var
asc_pln6_value = asc_pln6 + asc_pln6_shift_acres * acres + asc_pln6_shift_myield * yield_mean + asc_pln6_shift_vyield * yield_var
asc_pln7_value = asc_pln7 + asc_pln7_shift_acres * acres + asc_pln7_shift_myield * yield_mean + asc_pln7_shift_vyield * yield_var
asc_pln8_value = asc_pln8 + asc_pln8_shift_acres * acres + asc_pln8_shift_myield * yield_mean + asc_pln8_shift_vyield * yield_var
asc_pln9_value = asc_pln9 + asc_pln9_shift_acres * acres + asc_pln9_shift_myield * yield_mean + asc_pln9_shift_vyield * yield_var
asc_pln10_value = asc_pln10 + asc_pln10_shift_acres * acres + asc_pln10_shift_myield * yield_mean + asc_pln10_shift_vyield * yield_var
asc_pln11_value = asc_pln11 + asc_pln11_shift_acres * acres + asc_pln11_shift_myield * yield_mean + asc_pln11_shift_vyield * yield_var
asc_pln12_value = asc_pln12 + asc_pln12_shift_acres * acres + asc_pln12_shift_myield * yield_mean + asc_pln12_shift_vyield * yield_var
asc_pln13_value = asc_pln13 + asc_pln13_shift_acres * acres + asc_pln13_shift_myield * yield_mean + asc_pln13_shift_vyield * yield_var
asc_pln14_value = asc_pln14 + asc_pln14_shift_acres * acres + asc_pln14_shift_myield * yield_mean + asc_pln14_shift_vyield * yield_var
asc_pln15_value = asc_pln15 + asc_pln15_shift_acres * acres + asc_pln15_shift_myield * yield_mean + asc_pln15_shift_vyield * yield_var
asc_pln16_value = asc_pln16 + asc_pln16_shift_acres * acres + asc_pln16_shift_myield * yield_mean + asc_pln16_shift_vyield * yield_var
asc_pln17_value = asc_pln17 + asc_pln17_shift_acres * acres + asc_pln17_shift_myield * yield_mean + asc_pln17_shift_vyield * yield_var
asc_pln18_value = asc_pln18 + asc_pln18_shift_acres * acres + asc_pln18_shift_myield * yield_mean + asc_pln18_shift_vyield * yield_var
asc_pln19_value = asc_pln19 + asc_pln19_shift_acres * acres + asc_pln19_shift_myield * yield_mean + asc_pln19_shift_vyield * yield_var
asc_pln20_value = asc_pln20 + asc_pln20_shift_acres * acres + asc_pln20_shift_myield * yield_mean + asc_pln20_shift_vyield * yield_var
asc_pln21_value = asc_pln21 + asc_pln21_shift_acres * acres + asc_pln21_shift_myield * yield_mean + asc_pln21_shift_vyield * yield_var
### Loop over classes
for(s in 1:2){
### Compute class-specific utilities
V=list()
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V[["pln1"]] = asc_pln1 + b_prm[[s]] * prm1 + b_sub[[s]] * sub1 + b_lia[[s]] * lia1
V[["pln2"]] = asc_pln2_value + b_prm[[s]] * prm2 + b_sub[[s]] * sub2 + b_lia[[s]] * lia2
V[["pln3"]] = asc_pln3_value + b_prm[[s]] * prm3 + b_sub[[s]] * sub3 + b_lia[[s]] * lia3
V[["pln4"]] = asc_pln4_value + b_prm[[s]] * prm4 + b_sub[[s]] * sub4 + b_lia[[s]] * lia4
V[["pln5"]] = asc_pln5_value + b_prm[[s]] * prm5 + b_sub[[s]] * sub5 + b_lia[[s]] * lia5
V[["pln6"]] = asc_pln6_value + b_prm[[s]] * prm6 + b_sub[[s]] * sub6 + b_lia[[s]] * lia6
V[["pln7"]] = asc_pln7_value + b_prm[[s]] * prm7 + b_sub[[s]] * sub7 + b_lia[[s]] * lia7
V[["pln8"]] = asc_pln8_value + b_prm[[s]] * prm8 + b_sub[[s]] * sub8 + b_lia[[s]] * lia8
V[["pln9"]] = asc_pln9_value + b_prm[[s]] * prm9 + b_sub[[s]] * sub9 + b_lia[[s]] * lia9
V[["pln10"]] = asc_pln10_value + b_prm[[s]] * prm10 + b_sub[[s]] * sub10 + b_lia[[s]] * lia10
V[["pln11"]] = asc_pln11_value + b_prm[[s]] * prm11 + b_sub[[s]] * sub11 + b_lia[[s]] * lia11
V[["pln12"]] = asc_pln12_value + b_prm[[s]] * prm12 + b_sub[[s]] * sub12 + b_lia[[s]] * lia12
V[["pln13"]] = asc_pln13_value + b_prm[[s]] * prm13 + b_sub[[s]] * sub13 + b_lia[[s]] * lia13
V[["pln14"]] = asc_pln14_value + b_prm[[s]] * prm14 + b_sub[[s]] * sub14 + b_lia[[s]] * lia14
V[["pln15"]] = asc_pln15_value + b_prm[[s]] * prm15 + b_sub[[s]] * sub15 + b_lia[[s]] * lia15
V[["pln16"]] = asc_pln16_value + b_prm[[s]] * prm16 + b_sub[[s]] * sub16 + b_lia[[s]] * lia16
V[["pln17"]] = asc_pln17_value + b_prm[[s]] * prm17 + b_sub[[s]] * sub17 + b_lia[[s]] * lia17
V[["pln18"]] = asc_pln18_value + b_prm[[s]] * prm18 + b_sub[[s]] * sub18 + b_lia[[s]] * lia18
V[["pln19"]] = asc_pln19_value + b_prm[[s]] * prm19 + b_sub[[s]] * sub19 + b_lia[[s]] * lia19
V[["pln20"]] = asc_pln20_value + b_prm[[s]] * prm20 + b_sub[[s]] * sub20 + b_lia[[s]] * lia20
V[["pln21"]] = asc_pln21_value + b_prm[[s]] * prm21 + b_sub[[s]] * sub21 + b_lia[[s]] * lia21
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)
### Average across inter-individual draws within classes
P[[paste0("Class_",s)]] = apollo_avgInterDraws(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 ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
apollo_saveOutput(model)