Code: Select all
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc_cin_1 = 0, asc_cin_2 = 0,
asc_cout_1 = 0, asc_cout_2 = 0,
asc_tp_1 = 0, asc_tp_2 = 0,
asc_pr_1 = 0, asc_pr_2 = 0,
asc_bike_1 = 0, asc_bike_2 = 0,
asc_walk_1 = 0, asc_walk_2 = 0,
asc_notravel_1 = 0, asc_notravel_2 = 0,
b_cost_cin_1 = 0, b_cost_cin_2 = 0,
b_cost_cout_1 = 0, b_cost_cout_2 = 0,
b_cost_tp_1 = 0, b_cost_tp_2 = 0,
b_cost_pr_1 = 0, b_cost_pr_2 = 0,
b_tvia_c_1 = 0, b_tvia_c_2 = 0,
b_tvia_tp_1 = 0, b_tvia_tp_2 = 0,
b_tvia_pr_1 = 0, b_tvia_pr_2 = 0,
b_tvia_bike_1 = 0, b_tvia_bike_2 = 0,
b_tvia_walk_1 = 0, b_tvia_walk_2 = 0,
b_tbus_cin_1 = 0, b_tbus_cin_2 = 0,
b_tbus_cout_1 = 0, b_tbus_cout_2 = 0,
b_tesp_1 = 0, b_tesp_2 = 0,
ETA_Acept_1 = 1, ETA_Acept_2 = 1,
#gamma_MALE = 0,
#gamma_OPINION_ZBE_GOOD = 0,
#gamma_LICENSE = 0,
#gamma_CocheSost = 0,
#gamma_FREQ_ZBE_RESIDENTE = 0,
#gamma_DISPOSICION_CAMBIO_ZBE_SI = 0,
#gamma_MOTIVO_TRABAJO = 0,
# CLASES PARA EL LC
delta_1 = 0, delta_2 = 0,
# Parámetros HDC
Zeta_Acept1 = 0,
Zeta_Acept2 = 0,
Zeta_Acept3 = 0,
Zeta_Acept4 = 0,
Zeta_Acept5 = 0,
Zeta_Acept10 = 0,
tau_Acept1_1 = 1,
tau_Acept1_2 = 2,
tau_Acept1_3 = 3,
tau_Acept1_4 = 4,
tau_Acept2_1 = 1,
tau_Acept2_2 = 2,
tau_Acept2_3 = 3,
tau_Acept2_4 = 4,
tau_Acept3_1 = 1,
tau_Acept3_2 = 2,
tau_Acept3_3 = 3,
tau_Acept3_4 = 4,
tau_Acept4_1 = 1,
tau_Acept4_2 = 2,
tau_Acept4_3 = 3,
tau_Acept4_4 = 4,
tau_Acept5_1 = 1,
tau_Acept5_2 = 2,
tau_Acept5_3 = 3,
tau_Acept5_4 = 4,
tau_Acept10_1 = 1,
tau_Acept10_2 = 2,
tau_Acept10_3 = 3,
tau_Acept10_4 = 4,
sigma_panel_1 = 1.2, sigma_panel_2 = 1.2
#sigma_public = 1,
#sigma_coche = 1
#b_tespe_tp = 0,
#b_tesp_pr = 0
#b_pend_bike = 0,
#b_pend_walk = 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_notravel_1","asc_notravel_2")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("eta1","eta2", "eta3","eta4", "eta5","eta6", "eta7", "eta11"), # "eta8", "eta9", "eta10",
intraDrawsType = "mlhs",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff=function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["ec1_1"]] = sigma_panel_1 * eta1
randcoeff[["ec1_2"]] = sigma_panel_2 * eta1
randcoeff[["ec2_1"]] = sigma_panel_1 * eta2
randcoeff[["ec2_2"]] = sigma_panel_2 * eta2
randcoeff[["ec3_1"]] = sigma_panel_1 * eta3
randcoeff[["ec3_2"]] = sigma_panel_2 * eta3
randcoeff[["ec4_1"]] = sigma_panel_1 * eta4
randcoeff[["ec4_2"]] = sigma_panel_2 * eta4
randcoeff[["ec5_1"]] = sigma_panel_1 * eta5
randcoeff[["ec5_2"]] = sigma_panel_2 * eta5
randcoeff[["ec6_1"]] = sigma_panel_1 * eta6
randcoeff[["ec6_2"]] = sigma_panel_2 * eta6
randcoeff[["ec7_1"]] = sigma_panel_1 * eta7
randcoeff[["ec7_2"]] = sigma_panel_2 * eta7
# randcoeff[["ec8"]] = sigma_coche * eta8
# randcoeff[["ec9"]] = sigma_public * eta9
# randcoeff[["ec10"]] = sigma_active * eta10
randcoeff[["Acept"]] = eta11
# gamma_MALE * MALE + gamma_OPINION_ZBE_GOOD * OPINION_ZBE_GOOD + gamma_LICENSE * LICENSE +
#gamma_CocheSost * CocheSost + gamma_FREQ_ZBE_RESIDENTE * FREQ_ZBE_RESIDENTE + gamma_DISPOSICION_CAMBIO_ZBE_SI * DISPOSICION_CAMBIO_ZBE_SI +
#gamma_MOTIVO_TRABAJO * MOTIVO_TRABAJO
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars = function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["asc_cin"]] = list(asc_cin_1, asc_cin_2)
lcpars[["asc_cout"]] = list(asc_cout_1, asc_cout_2)
lcpars[["asc_tp"]] = list(asc_tp_1, asc_tp_2)
lcpars[["asc_pr"]] = list(asc_pr_1, asc_pr_2)
lcpars[["asc_bike"]] = list(asc_bike_1, asc_bike_2)
lcpars[["asc_walk"]] = list(asc_walk_1, asc_walk_2)
lcpars[["asc_notravel"]] = list(asc_notravel_1, asc_notravel_2)
lcpars[["b_cost_cin"]] = list(b_cost_cin_1, b_cost_cin_2)
lcpars[["b_cost_cout"]] = list(b_cost_cout_1, b_cost_cout_2)
lcpars[["b_cost_tp"]] = list(b_cost_tp_1, b_cost_tp_2)
lcpars[["b_cost_pr"]] = list(b_cost_pr_1, b_cost_pr_2)
lcpars[["b_tvia_c"]] = list(b_tvia_c_1, b_tvia_c_2)
lcpars[["b_tvia_tp"]] = list(b_tvia_tp_1, b_tvia_tp_2)
lcpars[["b_tvia_pr"]] = list(b_tvia_pr_1, b_tvia_pr_2)
lcpars[["b_tvia_bike"]] = list(b_tvia_bike_1, b_tvia_bike_2)
lcpars[["b_tvia_walk"]] = list(b_tvia_walk_1, b_tvia_walk_2)
lcpars[["b_tbus_cin"]] = list(b_tbus_cin_1, b_tbus_cin_2)
lcpars[["b_tbus_cout"]] = list(b_tbus_cout_1, b_tbus_cout_2)
lcpars[["b_tesp"]] = list(b_tesp_1, b_tesp_2)
lcpars[["ETA_Acept"]] = list(ETA_Acept_1, ETA_Acept_2)
#lcpars[["Acept"]] = list(Acept_1, Acept_2)
lcpars[["ec1"]] = list(ec1_1, ec1_2)
lcpars[["ec2"]] = list(ec2_1, ec2_2)
lcpars[["ec3"]] = list(ec3_1, ec3_2)
lcpars[["ec4"]] = list(ec4_1, ec4_2)
lcpars[["ec5"]] = list(ec5_1, ec5_2)
lcpars[["ec6"]] = list(ec6_1, ec6_2)
lcpars[["ec7"]] = list(ec7_1, ec7_2)
V=list()
V[["class_a"]] = delta_1
V[["class_b"]] = delta_2
classAlloc_settings = list(
classes = c(class_a=1, class_b=2),
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 para cada individuo
P = list()
### Likelihood of indicators - Modelos de medicion
{ ol_settings1 = list(outcomeOrdered=latente1,
V=Zeta_Acept1 * Acept,
tau=c(tau_Acept1_1, tau_Acept1_2, tau_Acept1_3, tau_Acept1_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente1")}
{ ol_settings2 = list(outcomeOrdered=latente2,
V=Zeta_Acept2 * Acept,
tau=c(tau_Acept2_1, tau_Acept2_2, tau_Acept2_3, tau_Acept2_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente2")}
{ ol_settings3 = list(outcomeOrdered=latente3,
V=Zeta_Acept3 * Acept,
tau=c(tau_Acept3_1, tau_Acept3_2, tau_Acept3_3, tau_Acept3_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente3")}
{ ol_settings4 = list(outcomeOrdered=latente4,
V=Zeta_Acept4 * Acept,
tau=c(tau_Acept4_1, tau_Acept4_2, tau_Acept4_3, tau_Acept4_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente4")}
{ ol_settings5 = list(outcomeOrdered=latente5,
V=Zeta_Acept5 * Acept,
tau=c(tau_Acept5_1, tau_Acept5_2, tau_Acept5_3, tau_Acept5_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente5")}
{ ol_settings10 = list(outcomeOrdered=latente10,
V=Zeta_Acept10 * Acept,
tau=c(tau_Acept10_1, tau_Acept10_2, tau_Acept10_3, tau_Acept10_4 ),
coding = c(1,2,3,4,5),
rows=(OBSERVACION==1),
componentName = "latente10")}
{
P[["latente1"]] = apollo_ol(ol_settings1, functionality)
P[["latente2"]] = apollo_ol(ol_settings2, functionality)
P[["latente3"]] = apollo_ol(ol_settings3, functionality)
P[["latente4"]] = apollo_ol(ol_settings4, functionality)
P[["latente5"]] = apollo_ol(ol_settings5, functionality)
P[["latente10"]] = apollo_ol(ol_settings10, functionality)
}
S <- 2
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(cin=1,cout=2, tp=3, pr=4, bike=5, walk=6, notravel=99),
avail = list(cin= C_IN_AVAIL,cout= C_OUT_AVAIL, tp=TP_AVAIL, pr=PR_AVAIL, bike=BIKE_AVAIL, walk=WALK_AVAIL, notravel=1),
choiceVar = choice
)
for(s in 1:S){
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[["cin"]] = asc_cin[[s]] + b_cost_cin[[s]] * C_in_cost + b_tvia_c[[s]] * C_in_tvia + b_tbus_cin[[s]] * C_in_tbusq + ec1[[s]] + ETA_Acept[[s]] * Acept #+ ec8
V[["cout"]] = asc_cout[[s]] + b_cost_cout[[s]] * C_out_cost + b_tvia_c[[s]] * C_out_tvia + b_tbus_cout[[s]] * C_out_tbusq + ec2[[s]] #+ ec8
V[["tp"]] = asc_tp[[s]] + b_cost_tp[[s]] * TP_cost + b_tvia_tp[[s]] * TP_tvia + b_tesp[[s]] * TP_tespe + ec3[[s]] #+ec9
V[["pr"]] = asc_pr[[s]] + b_cost_pr[[s]] * PR_cost + b_tvia_pr[[s]] * PR_tvia + b_tesp[[s]] * PR_tespe + ec4[[s]] #+ec9
V[["bike"]] = asc_bike[[s]] + b_tvia_bike[[s]] * Bici_tvia + ec5[[s]] #+ Beta_pend_bike * LLANO_BIKE
V[["walk"]] = asc_walk[[s]] + b_tvia_walk[[s]] * Andar_tvia + ec6[[s]] #+ Beta_pend_walk * LLANO_WALK
V[["notravel"]] = asc_notravel[[s]] + ec7[[s]]
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[paste0("Class_", 1:S)], classProb=pi_values)
P[["choice"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
### Comment out as necessary
P = apollo_combineModels(P, apollo_inputs, functionality)
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
#### Compute latent class model probabilities
lc_settings = list(inClassProb = P[paste0("Class_", 1:S)], classProb=pi_values)
P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
### Average across inter-individual draws in class allocation probabilities
# P[["model"]] = apollo_combineModels(P[["model"]], apollo_inputs, functionality)
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)