I have a question about my model estimation. I have already estimated the MNL and MXL models, and now I’m trying to estimate the HCM. Specifically, I would like to estimate:
I’ve attempted to estimate the model, but I keep running into errors. Would you be able to advise me on how to approach this, or point me towards what might be going wrong?
Code: Select all
apollo_control = list(
modelName = "mmnl241a", # or any name you like
modelDescr = "mmnl for trust in green hotels",
indivID = "ID", # replace with your actual respondent ID column
nCores = 7,
outputDirectory = "output") # or more, depending on your machine
database$BP_1=database$BP_1-mean(database$BP_1)
database$BP_2=database$BP_2-mean(database$BP_2)
database$BP_3=database$BP_3-mean(database$BP_3)
database$BP_4=database$BP_4-mean(database$BP_4)
database$BP_5=database$BP_5-mean(database$BP_5)
database$BP_6=database$BP_6-mean(database$BP_6)
database$BP_7=database$BP_7-mean(database$BP_7)
database$SC_1=database$SC_1-mean(database$SC_1)
database$SC_2=database$SC_2-mean(database$SC_2)
database$SC_3=database$SC_3-mean(database$SC_3)
database$SC_4=database$SC_4-mean(database$SC_4)
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(mu_log_b_HT2 = 0.1526, mu_log_b_HT3 = 0.1085,
sigma_log_b_HT2 = 0.08571, sigma_log_b_HT3 = 0.08259, # Use abs(3.1) if you had a valid positive estimate
mu_log_b_GM2 = 0.5232, mu_log_b_GM3 = 0.4581,
sigma_log_b_GM2 = 0.08760, sigma_log_b_GM3 = 0.08471,
mu_log_b_TC2 = 0.8261, mu_log_b_TC3 = 1.0635,
sigma_log_b_TC2 = 0.07089, sigma_log_b_TC3 = 0.09575,
lambda_BP = 0, lambda_BP_SC = 0,
zeta_BP_1 = 0.54592, zeta_BP_2 = 0.46008, zeta_BP_3 = 0.31626, zeta_BP_4 = 0.51507,
zeta_BP_5 = 0.54111, zeta_BP_6 = 0.38038, zeta_BP_7 = 0.47925,
zeta_SC_1 = 0.35201, zeta_SC_2 = 0.31958, zeta_SC_3 = 0.83219, zeta_SC_4 = 0.74095,
sigma_BP1 = 1, sigma_BP2 = 1, sigma_BP3 = 1, sigma_BP4 = 1,sigma_BP5 = 1, sigma_BP6 = 1,sigma_BP7 = 1,
sigma_SC1 = 1, sigma_SC2 = 1, sigma_SC3 = 1, sigma_SC4 = 1,
gamma_SC_univ = 0, gamma_SC_reason = 0, gamma_SC_age = 0, gamma_HT = 0, gamma_TC = 0, gamma_GM = 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()
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "sobol",
interNDraws = 2000,
interUnifDraws = c(),
interNormDraws = c("eta_BP", "eta_SC","draws_HT2", "draws_GM2", "draws_TC2","draws_HT3", "draws_GM3", "draws_TC3"),
intraDrawsType = "sobol",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_HT2"]] = mu_log_b_HT2 + sigma_log_b_HT2 * draws_HT2
randcoeff[["b_GM2"]] = mu_log_b_GM2 + sigma_log_b_GM2 * draws_GM2
randcoeff[["b_TC2"]] = mu_log_b_TC2 + sigma_log_b_TC2 * draws_TC2
randcoeff[["b_HT3"]] = mu_log_b_HT3 + sigma_log_b_HT3 * draws_HT3
randcoeff[["b_GM3"]] = mu_log_b_GM3 + sigma_log_b_GM3 * draws_GM3
randcoeff[["b_TC3"]] = mu_log_b_TC3 + sigma_log_b_TC3 * draws_TC3
randcoeff[["BP"]] = gamma_HT*HT + gamma_TC*TC + gamma_GM*GM + eta_BP
randcoeff[["SC"]] = gamma_SC_univ*Education + gamma_SC_reason*Reason + gamma_SC_age*Age + eta_SC
return(randcoeff)
}
# ################################################################# #
#### 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()
normalDensity_settings1 = list(outcomeNormal = BP_1,
xNormal = zeta_BP_1*BP,
mu = 0,
sigma = sigma_BP1,
rows = (task==1),
componentName = "indic_BP_1")
normalDensity_settings2 = list(outcomeNormal = BP_2,
xNormal = zeta_BP_2*BP,
mu = 0,
sigma = sigma_BP2,
rows = (task==1),
componentName = "indic_BP_2")
normalDensity_settings3 = list(outcomeNormal = BP_3,
xNormal = zeta_BP_3*BP,
mu = 0,
sigma = sigma_BP3,
rows = (task==1),
componentName = "indic_BP_3")
normalDensity_settings4 = list(outcomeNormal = BP_4,
xNormal = zeta_BP_4*BP,
mu = 0,
sigma = sigma_BP4,
rows = (task==1),
componentName = "indic_BP_4")
normalDensity_settings5 = list(outcomeNormal = BP_5,
xNormal = zeta_BP_5*BP,
mu = 0,
sigma = sigma_BP5,
rows = (task==1),
componentName = "indic_BP_5")
normalDensity_settings6 = list(outcomeNormal = BP_6,
xNormal = zeta_BP_6*BP,
mu = 0,
sigma = sigma_BP6,
rows = (task==1),
componentName = "indic_BP_6")
normalDensity_settings7 = list(outcomeNormal = BP_7,
xNormal = zeta_BP_7*BP,
mu = 0,
sigma = sigma_BP7,
rows = (task==1),
componentName = "indic_BP_7")
normalDensity_settings8 = list(outcomeNormal = SC_1,
xNormal = zeta_SC_1*SC,
mu = 0,
sigma = sigma_SC1,
rows = (task==1),
componentName = "indic_SC_1")
normalDensity_settings9 = list(outcomeNormal = SC_2,
xNormal = zeta_SC_2*SC,
mu = 0,
sigma = sigma_SC2,
rows = (task==1),
componentName = "indic_SC_2")
normalDensity_settings10 = list(outcomeNormal = SC_3,
xNormal = zeta_SC_3*SC,
mu = 0,
sigma = sigma_SC3,
rows = (task==1),
componentName = "indic_SC_3")
normalDensity_settings11 = list(outcomeNormal = SC_4,
xNormal = zeta_SC_4*SC,
mu = 0,
sigma = sigma_SC4,
rows = (task==1),
componentName = "indic_SC_4")
P[["indic_BP_1"]] = apollo_normalDensity(normalDensity_settings1, functionality)
P[["indic_BP_1"]] = apollo_panelProd(P[["indic_BP_1"]] , apollo_inputs, functionality)
P[["indic_BP_2"]] = apollo_normalDensity(normalDensity_settings2, functionality)
P[["indic_BP_2"]] = apollo_panelProd(P[["indic_BP_2"]], apollo_inputs, functionality)
P[["indic_BP_3"]] = apollo_normalDensity(normalDensity_settings3, functionality)
P[["indic_BP_3"]] = apollo_panelProd(P[["indic_BP_3"]], apollo_inputs, functionality)
P[["indic_BP_4"]] = apollo_normalDensity(normalDensity_settings4, functionality)
P[["indic_BP_4"]] = apollo_panelProd(P[["indic_BP_4"]] , apollo_inputs, functionality)
P[["indic_BP_5"]] = apollo_normalDensity(normalDensity_settings5, functionality)
P[["indic_BP_5"]] = apollo_panelProd(P[["indic_BP_5"]] , apollo_inputs, functionality)
P[["indic_BP_6"]] = apollo_normalDensity(normalDensity_settings6, functionality)
P[["indic_BP_6"]] = apollo_panelProd(P[["indic_BP_6"]], apollo_inputs, functionality)
P[["indic_BP_7"]] = apollo_normalDensity(normalDensity_settings7, functionality)
P[["indic_BP_7"]] = apollo_panelProd(P[["indic_BP_7"]], apollo_inputs, functionality)
P[["indic_SC_1"]] = apollo_normalDensity(normalDensity_settings8, functionality)
P[["indic_SC_1"]] = apollo_panelProd(P[["indic_SC_1"]] , apollo_inputs, functionality)
P[["indic_SC_2"]] = apollo_normalDensity(normalDensity_settings9, functionality)
P[["indic_SC_2"]] = apollo_panelProd(P[["indic_SC_2"]], apollo_inputs, functionality)
P[["indic_SC_3"]] = apollo_normalDensity(normalDensity_settings10, functionality)
P[["indic_SC_3"]] = apollo_panelProd(P[["indic_SC_3"]], apollo_inputs, functionality)
P[["indic_SC_4"]] = apollo_normalDensity(normalDensity_settings11, functionality)
P[["indic_SC_4"]] = apollo_panelProd(P[["indic_SC_4"]] , apollo_inputs, functionality)
### Likelihood of choices
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[['alt1']] = (b_HT2*(HT_1=="2") + b_HT3*(HT_1=="3") + b_GM2*(GM_1=="2") + b_GM3*(GM_1=="3")+ b_TC2*(TC_1=="2") + b_TC3*(TC_1=="3")+lambda_BP*BP+lambda_BP_SC*BP*SC)
V[['alt2']] =(b_HT2*(HT_2=="2") + b_HT3*(HT_2=="3") + b_GM2*(GM_2=="2") + b_GM3*(GM_2=="3")+ b_TC2*(TC_2=="2") + b_TC3*(TC_2=="3")+lambda_BP*BP+lambda_BP_SC*BP*SC)
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2),
avail = list(alt1=1, alt2=1),
choiceVar = choice,
V = V
)
### Compute probabilities for MNL model component
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 ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs,
estimate_settings=list(estimationRoutine="bfgs",maxIterations=100000))