Code: Select all
parallel::detectCores()
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
library(foreach)
library(iterators)
library(doParallel)
### Initialise code
apollo_initialise()
# LN-OC,CT,R_b,CF_b&T-EM BFGS_1500
### Set core controls
apollo_control = list(
modelName = "MODEL_X_2000",
modelDescr = "MODEL_X_2000",
indivID = "Person_ID",
nCores = 10,
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("D:\\IST\\Project\\Data\\fourWheelerFaceToFace_V3(cleaning).csv", header = TRUE)
database1= read.csv("D:\\IST\\Project\\Data\\fourWheelerFaceToFace_V3(cleaning).csv", header = TRUE)
### for data dictionary, use ?apollo_swissRouteChoiceData
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1_a = 0.329,
asc_1_b = -0.682,
asc_2_a = 0,
asc_2_b = 0,
beta_pp_a = -0.005 ,
beta_pp_b = -0.015 ,
# beta_oc_a = 0,
# beta_oc_b = 0,
#
betaa_oc_a = -0.054,
# log_oc_a_mu = 0,
#
# log_oc_a_sig = 0,
log_oc_b_mu = -5.566,
log_oc_b_sig = 3.756,
#betaa_oc_b = 0,
# Jn_oc_a_mu = 0,
# Jn_oc_b_mu = 0,
# Jn_oc_a_sig = 0,
# Jn_oc_b_sig = 0,
# tr_oc_p_a = 0,
#tr_oc_p_b = 0,
# tr_oc_q_a = 0,
# tr_oc_q_b = 0,
# beta_ct_a = 0,
# beta_ct_b = 0,
betaa_ct_a = -0.020,
# log_ct_a_mu = 0,
#
# log_ct_a_sig = 0,
log_ct_b_mu = 5.488,
log_ct_b_sig =-16.790,
# betaa_ct_b = 0,
# tr_ct_p_a = 0,
#tr_ct_p_b = 0,
# tr_ct_q_a = 0,
# tr_ct_q_b = 0,
# betaa_r_a = 0,
# beta_r_b = 0,
# log_r_a_mu = 0,
# log_r_a_sig = 0,
# log_r_b_mu = 0,
#
# log_r_b_sig = 0,
# tr_r_p_a = 0,
# tr_r_p_b = 0,
tr_r_q_a = 0.129,
tr_r_q_b = 0.134,
# betaa_r_b = 0,
beta_cf_a = 0.112,
beta_cf_b = 0.337,
# log_cf_a_mu = 0,
# log_cf_a_sig = 0,
# log_cf_b_mu = 0,
#
# log_cf_b_sig = 0,
# tr_cf_p_a = 0,
# tr_cf_p_b = 0,
# betaa_cf_a = 0,
# tr_cf_q_a = 0,
# tr_cf_q_b = 0,
# betaa_cf_b = 0,
# beta_em_a = 0,
# beta_em_b = 0,
# log_em_a_mu = -1.9987,
# log_em_a_sig = -0.08649,
# log_em_b_mu = -10.6245,
#
# log_em_b_sig = -0.1149,
# tr_em_p_a = 0,
#tr_em_p_b = 0,
# tr_em_q_a = 0,
# betaa_em_a = 0,
# tr_em_q_b = 0,
#
beta_priorityLanes_a = 0.215,
beta_priorityLanes_b = -0.099,
beta_freeParking_a = 0.077,
beta_freeParking_b = 0.240,
beta_tollExemption_a = 0,
beta_tollExemption_b = 0,
#
delta_a = -0.841,
# gamma_knowledgeLow_a = 0,
# gamma_knowledgeMid_a = 0,
# gamma_KnowledgeHigh_a = 0,
gamma_voCoded_a = 0.347,
gamma_dailyDistanceLow_a =0 ,
gamma_dailyDistanceMid_a = -0.669,
gamma_dailyDistanceHigh_a = -1.397,
gamma_longDistanceLow_a = 0,
gamma_longDistanceMid_a = 1.013,
gamma_longDistanceHigh_a = 1.732,
# gamma_techEnth_a = 0,
# gamma_enviEnth_a = 0,
# gamma_socIma_a = 0,
# gamma_monBen_a = 0,
gamma_perFee_a = -0.339,
gamma_perRisk_a = -0.608,
# gamma_insUti_a = 0,
delta_b = 0,
# gamma_knowledgeLow_b = 0,
# gamma_knowledgeMid_b = 0,
# gamma_KnowledgeHigh_b = 0,
gamma_voCoded_b = 0,
gamma_dailyDistanceLow_b = 0,
gamma_dailyDistanceMid_b = 0,
gamma_dailyDistanceHigh_b = 0,
gamma_longDistanceLow_b = 0,
gamma_longDistanceMid_b = 0,
gamma_longDistanceHigh_b = 0,
# gamma_techEnth_b = 0,
# gamma_enviEnth_b = 0,
# gamma_socIma_b = 0,
# gamma_monBen_b = 0,
gamma_perFee_b = 0,
gamma_perRisk_b = 0
# gamma_insUti_b = 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_2_a","asc_2_b", "delta_b",
"gamma_dailyDistanceLow_a",
"gamma_longDistanceLow_a",
# "gamma_knowledgeLow_a",
# "gamma_knowledgeLow_b",
# "gamma_knowledgeMid_b",
# "gamma_KnowledgeHigh_b",
"gamma_voCoded_b" ,
"gamma_dailyDistanceLow_b" ,
"gamma_dailyDistanceMid_b",
"gamma_dailyDistanceHigh_b",
"gamma_longDistanceLow_b",
"gamma_longDistanceMid_b" ,
"gamma_longDistanceHigh_b",
"beta_tollExemption_a",
"beta_tollExemption_b",
# # "gamma_techEnth_b",
# # "gamma_enviEnth_b",
# # "gamma_socIma_b",
# # "gamma_monBen_b",
"gamma_perFee_b",
"gamma_perRisk_b"
# "gamma_insUti_b"
)
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="mlhs",
interNDraws= 2000,
interUnifDraws=c("draws_ct_a","draws_r_a", "draws_r_b"),
interNormDraws=c("draws_oc_a"),
# "draws_oc_a", "draws_oc_b", "draws_ct_a", "draws_ct_b", "draws_r_a", "draws_r_b", "draws_cf_a", "draws_cf_b" , "draws_em_a", "draws_em_b"
# "draws_oc_a", "draws_ct_a", "draws_r_a", "draws_cf_a", "draws_em_a"
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
# one class lognormal
randcoeff[["beta_oc_a"]] = betaa_oc_a + 0*(draws_oc_a)
randcoeff[["beta_oc_b"]] = -exp(log_oc_b_mu + log_oc_b_sig*draws_oc_a)
randcoeff[["beta_ct_a"]] = betaa_ct_a + 0*(draws_ct_a)
randcoeff[["beta_ct_b"]] = -exp(log_ct_b_mu + log_ct_b_sig*draws_ct_a)
# randcoeff[["beta_r_a"]] = (betaa_r_a + 0*draws_r_a)
# randcoeff[["beta_r_b"]] = exp(log_r_b_mu + log_r_b_sig*draws_r_a)
# randcoeff[["beta_cf_a"]] = (betaa_cf_a + 0*draws_cf_a)
# randcoeff[["beta_cf_b"]] = exp(log_cf_b_mu + log_cf_b_sig*draws_cf_a)
# randcoeff[["beta_em_a"]] = (log_em_a_mu + 0*draws_em_a)
# randcoeff[["beta_em_b"]] = -exp(log_em_b_mu + log_em_b_sig*draws_em_a)
# lognormal distribution
# randcoeff[["beta_oc_a"]] = -exp(log_oc_a_mu + log_oc_a_sig*draws_oc_a)
# randcoeff[["beta_oc_b"]] = -exp(log_oc_b_mu + log_oc_b_sig*draws_oc_a)
# randcoeff[["beta_ct_a"]] = -exp(log_ct_a_mu + log_ct_a_sig*draws_ct_a)
# randcoeff[["beta_ct_b"]] = -exp(log_ct_b_mu + log_ct_b_sig*draws_ct_a)
# randcoeff[["beta_r_a"]] = exp(log_r_a_mu + log_r_a_sig*draws_r_a)
# randcoeff[["beta_r_b"]] = exp(log_r_b_mu + log_r_b_sig*draws_r_a)
# randcoeff[["beta_cf_a"]] = exp(log_cf_a_mu + log_cf_a_sig*draws_cf_a)
# randcoeff[["beta_cf_b"]] = exp(log_cf_b_mu + log_cf_b_sig*draws_cf_a)
# randcoeff[["beta_em_a"]] = -exp(log_em_a_mu + log_em_a_sig*draws_em_a)
# randcoeff[["beta_em_b"]] = -exp(log_em_b_mu + log_em_b_sig*draws_em_a)
#triangular distribution
# randcoeff[["beta_oc_a"]] = tr_oc_q_a * (draws_oc_a + draws_oc_b)/2
# randcoeff[["beta_oc_b"]] = tr_oc_q_b *(draws_oc_a + draws_oc_b)/2
# randcoeff[["beta_ct_a"]] = tr_ct_q_a * (draws_ct_a + draws_ct_b)/2
# randcoeff[["beta_ct_b"]] = tr_ct_q_b *(draws_ct_a + draws_ct_b)/2
randcoeff[["beta_r_a"]] = tr_r_q_a * (draws_r_a + draws_r_b)/2
randcoeff[["beta_r_b"]] = tr_r_q_b*(draws_r_a + draws_r_b)/2
# randcoeff[["beta_cf_a"]] = tr_cf_q_a *(draws_cf_a + draws_cf_b)/2
# randcoeff[["beta_cf_b"]] = tr_cf_q_b * (draws_cf_a + draws_cf_b)/2
# randcoeff[["beta_em_a"]] = tr_em_q_a * (draws_em_a + draws_em_b)/2
# randcoeff[["beta_em_a"]] = tr_em_q_a* (draws_em_a + draws_em_b)/2
# randcoeff[["beta_em_b"]] = tr_em_q_b * (draws_em_a + draws_em_b)/2
#Johnson's distribution
# randcoeff[['beta_oc_a']] = exp(Jn_oc_a_mu + Jn_oc_a_sig*draws_oc_a)/(1 + exp(Jn_oc_a_mu + Jn_oc_a_sig*draws_oc_a))
# randcoeff[['beta_oc_b']] = exp(Jn_oc_b_mu + Jn_oc_b_sig*draws_oc_a)/(1 + exp(Jn_oc_b_mu + Jn_oc_b_sig*draws_oc_a))
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars = function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["asc_1"]] = list(asc_1_a, asc_1_b)
lcpars[["asc_2"]] = list(asc_2_a, asc_2_b)
lcpars[["beta_pp"]] = list(beta_pp_a, beta_pp_b)
lcpars[["beta_oc"]] = list(beta_oc_a, beta_oc_b)
lcpars[["beta_ct"]] = list(beta_ct_a, beta_ct_b)
lcpars[["beta_r"]] = list(beta_r_a, beta_r_b)
lcpars[["beta_cf"]] = list(beta_cf_a, beta_cf_b)
# lcpars[["beta_em"]] = list(beta_em_a, beta_em_b)
lcpars[["beta_priorityLanes"]] = list(beta_priorityLanes_a, beta_priorityLanes_b)
lcpars[["beta_freeParking"]] = list(beta_freeParking_a, beta_freeParking_b)
lcpars[["beta_tollExemption"]] = list(beta_tollExemption_a, beta_tollExemption_b)
V=list()
V[["class_a"]] = delta_a + gamma_voCoded_a*evoCoded + gamma_dailyDistanceLow_a*dailyDistanceLow + gamma_dailyDistanceMid_a*dailyDistanceMiddle+ gamma_dailyDistanceHigh_a*dailyDistanceHigh + gamma_longDistanceLow_a*longDistanceLow + gamma_longDistanceMid_a*longDistanceMiddle+ gamma_longDistanceHigh_a*longDistanceHigh + gamma_perFee_a*perFeeML + gamma_perRisk_a*perRiskML
#+ gamma_perFee_a*perFeeML + gamma_perRisk_a*perRiskML
#+ +gamma_voCoded_a*evoCoded
#+ gamma_knowledgeLow_a*knowledgeLow + gamma_knowledgeMid_a*knowledgeMiddle + gamma_KnowledgeHigh_a*knowledgeHigh
#+ gamma_monBen_a*monBenML + gamma_perFee_a*perFeeML + gamma_perRisk_a*perRiskML + gamma_insUti_a*insUtiML
V[["class_b"]] = delta_b +gamma_voCoded_b*evoCoded + gamma_dailyDistanceLow_b*dailyDistanceLow + gamma_dailyDistanceMid_b*dailyDistanceMiddle+ gamma_dailyDistanceHigh_b*dailyDistanceHigh + gamma_longDistanceLow_b*longDistanceLow + gamma_longDistanceMid_b*longDistanceMiddle+ gamma_longDistanceHigh_b*longDistanceHigh + gamma_perFee_b*perFeeML + gamma_perRisk_b*perRiskML
#+ gamma_perFee_b*perFeeML + gamma_perRisk_b*perRiskML
#+ + gamma_voCoded_b*evoCoded
#+ gamma_monBen_b*monBenML + gamma_perFee_b*perFeeML + gamma_perRisk_b*perRiskML + gamma_insUti_b*insUtiML
#+ gamma_knowledgeLow_b*knowledgeLow + gamma_knowledgeMid_b*knowledgeMiddle + gamma_KnowledgeHigh_b*knowledgeHigh + gamma_voCoded_b*evoCoded
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
P = list()
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2),
avail = list(alt1=1, alt2=1),
choiceVar = choiceVehicle
)
### Loop over classes
for(s in 1:2){
### Compute class-specific utilities
V=list()
V[["alt1"]] = asc_1[[s]] + beta_pp[[s]]*ppEVI + beta_oc[[s]]*ocEV+ beta_ct[[s]]*ctEV/60 + beta_r[[s]]*rEV/100 + beta_cf[[s]]*cfEV + beta_priorityLanes[[s]]*priorityLanes + beta_freeParking[[s]]*freeParking + beta_tollExemption[[s]]*tollExemption
V[["alt2"]] = asc_2[[s]] + beta_pp[[s]]*ppCVI + beta_oc[[s]]*ocCV + beta_r[[s]]*rCV/100
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 AND OUTPUT ####
# ################################################################# #
### Estimate model
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings = list(estimationRoutine="BFGS",maxIterations= 5000))
### Show output in screen
apollo_modelOutput(model)
### Save output to file(s)
apollo_saveOutput(model)