-
-
Save HyewonNamgung/a84cae8ba6ee3a5c8e631ddadb32ff55 to your computer and use it in GitHub Desktop.
The ICLV code for binary mixed logit model with 4 latent variables.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ### Clear memory | |
| rm(list = ls()) | |
| ### Load Apollo library | |
| library(apollo) | |
| ### Initialise code | |
| apollo_initialise() | |
| ### Set core controls | |
| apollo_control = list( | |
| modelName = "Business model ICLV", | |
| modelDescr = "Business model with ICLV", | |
| indivID = "ID", | |
| mixing = TRUE, | |
| nCores = 14 | |
| ) | |
| # ################################################################# # | |
| #### LOAD DATA AND APPLY ANY TRANSFORMATIONS #### | |
| # ################################################################# # | |
| database = read.csv("~/OneDrive - /04.Data/operation_model_1012.csv",header=TRUE) | |
| database$task <- rep(c(1,2,3,4),times = 245) | |
| ## data prepare | |
| database <- database[database$choice!=-1,] | |
| database <- database[database$Q5_c & database$Q5_d &database$Q5_e &database$Q5_h!=-1,] | |
| #database <- database[database$city_level==6|database$city_level==5,] | |
| database$urban <- ifelse(database$city_level==1,1,0) <- ifelse(database$city_level==1,1,0) | |
| database$middle <- ifelse(database$city_level==2|database$city_level==3,1,0) | |
| #database$kure <- ifelse(database$city_level==3,1,0) | |
| database$tenman <-ifelse(database$city_level==4|database$city_level==5,1,0) | |
| #database$tenman_down <- ifelse(database$city_level==5,1,0) | |
| database$rural <- ifelse(database$city_level==6,1,0) | |
| database$sustain1 <- ifelse(database$operation==5,1,0) | |
| database$sustain2 <- ifelse(database$operation==6,1,0) | |
| database$sustain3 <- ifelse(database$operation==7,1,0) | |
| database$sustain4 <- ifelse(database$operation==8,1,0) | |
| database$tech_infra <-ifelse(database$tech_infra==1,1,0) | |
| #database$license <-ifelse(database$license_1 |database$license_2=="yes",1,0) | |
| colnames(database)[1] <- "ID" #change variable name error #change variable name error | |
| ### Subtract mean of indicator variables to centre them on zero | |
| # LV_1 : operation cost | |
| database$Q5_a=database$Q5_a-mean(database$Q5_a) | |
| database$Q5_b=database$Q5_b-mean(database$Q5_b) | |
| # LV_2 : concern about regulation | |
| database$Q5_g=database$Q5_g-mean(database$Q5_g) | |
| database$Q5_h=database$Q5_h-mean(database$Q5_h) | |
| database$Q5_i=database$Q5_i-mean(database$Q5_i) | |
| # LV_3 : thoughts about ride sharing | |
| database$Q12_9=database$Q12_9-mean(database$Q12_9) | |
| database$Q12_10=database$Q12_10-mean(database$Q12_10) | |
| database$Q12_11=database$Q12_11-mean(database$Q12_11) | |
| database$Q12_12=database$Q12_12-mean(database$Q12_12) | |
| # ################################################################# # | |
| #### DEFINE MODEL PARAMETERS #### | |
| # ################################################################# # | |
| ### Vector of parameters, including any that are kept fixed in estimation | |
| apollo_beta = c(asc_origin = 0, | |
| asc_new = 0, | |
| b_license = 0, | |
| b_license_no = 0, | |
| b_dynamic = 0, | |
| b_fare_now = 0, | |
| b_transport = 0, | |
| b_taxi = 0, | |
| b_private = 0, | |
| b_own = 0, | |
| b_goods = 0, | |
| b_passenger = 0, | |
| b_share_no = 0, | |
| b_half = 0, | |
| b_fullav = 0, | |
| b_noav = 0, | |
| lambda_operation = 0, | |
| lambda_regulation = 0, | |
| lambda_rideshare = 0, | |
| lambda_sustain = 0, | |
| gamma_urban_operation = 0, | |
| gamma_middle_operation = 0, | |
| gamma_rural_operation = 0, | |
| gamma_urban_regulation = 0, | |
| gamma_middle_regulation = 0, | |
| gamma_rural_regulation = 0, | |
| gamma_urban_rideshare = 0, | |
| gamma_middle_rideshare= 0, | |
| gamma_rural_rideshare = 0, | |
| gamma_urban_sustain = 0, | |
| gamma_middle_sustain= 0, | |
| gamma_rural_sustain = 0, | |
| #gamma_size_sustain = 0, | |
| gamma_tech_infra_sustain = 0, | |
| zeta_operation_1 = 1, | |
| zeta_operation_2 = 1, | |
| zeta_regulation_1 = 1, | |
| zeta_regulation_2 = 1, | |
| zeta_regulation_3 = 1, | |
| zeta_rideshare_1 = 1, | |
| zeta_rideshare_2 = 1, | |
| zeta_rideshare_3 = 1, | |
| zeta_rideshare_4 = 1, | |
| zeta_sustain_1 = 1, | |
| zeta_sustain_2 = 1, | |
| zeta_sustain_3 = 1, | |
| zeta_sustain_4 = 1, | |
| sigma_operation_1 = 1, | |
| sigma_operation_2 = 1, | |
| sigma_rideshare_1 = 1, | |
| sigma_rideshare_2 = 1, | |
| sigma_rideshare_3 = 1, | |
| sigma_rideshare_4 = 1, | |
| sigma_regulation_1 = 1, | |
| sigma_regulation_2 = 1, | |
| sigma_regulation_3 = 1, | |
| sigma_sustain_1 = 1, | |
| sigma_sustain_2 = 1, | |
| sigma_sustain_3 = 1, | |
| sigma_sustain_4 = 1, | |
| sigma_panel1 = 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_origin", "b_license_no", "b_fare_now","b_own", "b_share_no", "b_noav", "gamma_rural_regulation", "gamma_rural_rideshare", | |
| "gamma_rural_operation", "gamma_rural_sustain") | |
| # ################################################################# # | |
| #### DEFINE RANDOM COMPONENTS #### | |
| # ################################################################# # | |
| ### Set parameters for generating draws | |
| apollo_draws = list( | |
| interDrawsType="mlhs", | |
| interNDraws=500, | |
| interUnifDraws=c(), | |
| interNormDraws=c("eta_operation", "eta_regulation", "eta_rideshare", "eta_sustain", "draws_panel_1","draws_panel_2"), | |
| intraDrawsType= "halton", | |
| intraNDraws=0, | |
| intraUnifDraws= c(), | |
| intraNormDraws = c () | |
| ) | |
| ### Create random parameters | |
| apollo_randCoeff=function(apollo_beta, apollo_inputs){ | |
| randcoeff = list() | |
| randcoeff[["LV_operation"]] = gamma_urban_operation*hiroshima + gamma_middle_operation*middle + gamma_rural_operation*rural + eta_operation | |
| randcoeff[["LV_regulation"]] = gamma_urban_regulation*hiroshima + gamma_middle_regulation*middle + gamma_rural_regulation*rural + eta_regulation | |
| randcoeff[["LV_rideshare"]] = gamma_urban_rideshare*hiroshima + gamma_middle_rideshare*middle + gamma_rural_rideshare*rural + eta_rideshare | |
| randcoeff[["LV_sustain"]] = gamma_urban_sustain*hiroshima + gamma_middle_sustain*middle + gamma_rural_sustain*rural + gamma_tech_infra_sustain*tech_infra + eta_sustain | |
| randcoeff[["ec1"]] = sigma_panel1 * draws_panel_1 | |
| 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() | |
| ### Likelihood of indicators | |
| normalDensity_settings1 = list(outcomeNormal = Q5_a, | |
| xNormal = zeta_operation_1*LV_operation, | |
| mu = 0, | |
| sigma = sigma_operation_1, | |
| rows = (task==1), | |
| componentName = "indi_operation_1") | |
| normalDensity_settings2 = list(outcomeNormal = Q5_b, | |
| xNormal = zeta_operation_2*LV_operation, | |
| mu = 0, | |
| sigma = sigma_operation_2, | |
| rows = (task==1), | |
| componentName = "indi_operation_2") | |
| normalDensity_settings3 = list(outcomeNormal = Q5_g, | |
| xNormal = zeta_regulation_1*LV_regulation, | |
| mu = 0, | |
| sigma = sigma_regulation_1, | |
| rows = (task==1), | |
| componentName = "indi_regulation_1") | |
| normalDensity_settings4 = list(outcomeNormal = Q5_h, | |
| xNormal = zeta_regulation_2*LV_regulation, | |
| mu = 0, | |
| sigma = sigma_regulation_2, | |
| rows = (task==1), | |
| componentName = "indi_regulation_2") | |
| normalDensity_settings5 = list(outcomeNormal = Q5_i, | |
| xNormal = zeta_regulation_3*LV_regulation, | |
| mu = 0, | |
| sigma = sigma_regulation_3, | |
| rows = (task==1), | |
| componentName = "indi_regulation_3") | |
| normalDensity_settings6 = list(outcomeNormal = Q12_9, | |
| xNormal = zeta_rideshare_1*LV_rideshare, | |
| mu = 0, | |
| sigma = sigma_rideshare_1, | |
| rows = (task==1), | |
| componentName = "indi_rideshare_1") | |
| normalDensity_settings7 = list(outcomeNormal = Q12_10, | |
| xNormal = zeta_rideshare_2*LV_rideshare, | |
| mu = 0, | |
| sigma = sigma_rideshare_2, | |
| rows = (task==1), | |
| componentName = "indi_rideshare_2") | |
| normalDensity_settings8 = list(outcomeNormal = Q12_11, | |
| xNormal = zeta_rideshare_3*LV_rideshare, | |
| mu = 0, | |
| sigma = sigma_rideshare_3, | |
| rows = (task==1), | |
| componentName = "indi_rideshare_3") | |
| normalDensity_settings9 = list(outcomeNormal = Q12_12, | |
| xNormal = zeta_rideshare_4*LV_rideshare, | |
| mu = 0, | |
| sigma = sigma_rideshare_4, | |
| rows = (task==1), | |
| componentName = "indi_rideshare_4") | |
| V = list() | |
| V[['no1']] = 0 | |
| V[['sustain1']] = zeta_sustain_1*LV_sustain | |
| mnl_settings1 = list( | |
| alternatives = c(no1=0, sustain1=1), | |
| avail = list(no1=1,sustain1=1), | |
| choiceVar = sustain1, | |
| V = V, | |
| componentName = "indi_sustain_1") | |
| V = list() | |
| V[['no2']] = 0 | |
| V[['sustain2']] = zeta_sustain_2*LV_sustain | |
| mnl_settings2 = list( | |
| alternatives = c(no2=0, sustain2=1), | |
| avail = list(no2=1,sustain2=1), | |
| choiceVar = sustain2, | |
| V = V, | |
| componentName = "indi_sustain_2") | |
| V = list() | |
| V[['no3']] = 0 | |
| V[['sustain3']] = zeta_sustain_3*LV_sustain | |
| mnl_settings3 = list( | |
| alternatives = c(no3=0, sustain3=1), | |
| avail = list(no3=1,sustain3=1), | |
| choiceVar = sustain3, | |
| V = V, | |
| componentName = "indi_sustain_3") | |
| V = list() | |
| V[['no4']] = 0 | |
| V[['sustain4']] = zeta_sustain_4*LV_sustain | |
| mnl_settings4 = list( | |
| alternatives = c(no4=0, sustain4=1), | |
| avail = list(no4=1,sustain4=1), | |
| choiceVar = sustain4, | |
| V = V, | |
| componentName = "indi_sustain_4") | |
| P[["indi_operation_1"]] = apollo_normalDensity(normalDensity_settings1, functionality) | |
| P[["indi_operation_2"]] = apollo_normalDensity(normalDensity_settings2, functionality) | |
| P[["indi_regulation_1"]] = apollo_normalDensity(normalDensity_settings3, functionality) | |
| P[["indi_regulation_2"]] = apollo_normalDensity(normalDensity_settings4, functionality) | |
| P[["indi_regulation_3"]] = apollo_normalDensity(normalDensity_settings5, functionality) | |
| P[["indi_rideshare_1"]] = apollo_normalDensity(normalDensity_settings6, functionality) | |
| P[["indi_rideshare_2"]] = apollo_normalDensity(normalDensity_settings7, functionality) | |
| P[["indi_rideshare_3"]] = apollo_normalDensity(normalDensity_settings8, functionality) | |
| P[["indi_rideshare_4"]] = apollo_normalDensity(normalDensity_settings9, functionality) | |
| P[["indi_sustain_1"]] = apollo_mnl(mnl_settings1, functionality) | |
| P[["indi_sustain_2"]] = apollo_mnl(mnl_settings2, functionality) | |
| P[["indi_sustain3"]] = apollo_mnl(mnl_settings3, functionality) | |
| P[["indi_sustain4"]] = apollo_mnl(mnl_settings4, functionality) | |
| ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant | |
| V = list() | |
| V[['origin']] = 0 | |
| V[['new']] = (asc_new + b_license*(license_2=="yes") + b_license_no*(license_2=="no")+b_fare_now*(fare_2=="now")+ b_dynamic*(fare_2=="dynamic")+ b_transport*(coop_2=="public") | |
| + b_taxi*(coop_2=="taxi") + b_private*(coop_2=="private") + b_own*(coop_2=="no") + b_goods * (share_2=="goods") | |
| + b_passenger*(share_2 =="people") + b_share_no *(share_2=="no")+ b_half*(AV_2=="lv4") + b_fullav*(AV_2== "lv5") | |
| + b_noav*(AV_2=="no")+lambda_operation*LV_operation+lambda_regulation*LV_regulation+lambda_rideshare*LV_rideshare+lambda_sustain*LV_sustain+ec1) | |
| ### Define settings for MNL model component | |
| mnl_settings = list( | |
| alternatives = c(origin=1,new =2), | |
| avail = list(origin=1, new=1), | |
| choiceVar = choice, | |
| V = V | |
| ) | |
| ### Compute probabilities using MNL model | |
| 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 #### | |
| # ################################################################# # | |
| ### Optional: calculate LL before model estimation | |
| apollo_llCalc(apollo_beta, apollo_probabilities, apollo_inputs) | |
| ## Estimate model | |
| model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs) | |
| # ################################################################# # | |
| #### MODEL OUTPUTS #### | |
| # ################################################################# # | |
| # ----------------------------------------------------------------- # | |
| #---- FORMATTED OUTPUT (TO SCREEN) ---- | |
| # ----------------------------------------------------------------- # | |
| apollo_modelOutput(model, modelOutput_settings = list(printPVal=2)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment