Hybrid Model for Bicycle to Work

Hybrid_model.r

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Clear memory
rm(list = ls())

### Load libraries
library(apollo)
library(tidyverse)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName = "Hybrid Model for Bicycle to Work",
  modelDescr = "ICLV model",
  indivID = "id",
  mixing = TRUE,
  nCores = 4,
  outputDirectory = "output"
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database = read.csv("ils-biketoworkdata.csv", header=TRUE)
#database <- na.omit(database)
#View(database)

### Subtract mean of indicator variables to center them on zero

database$bikelane     = database$bikelane - mean(database$bikelane) 
database$pavement     = database$pavement - mean(database$pavement)
database$safety       = database$safety - mean(database$safety)
database$scenery      = database$scenery - mean(database$scenery)
database$shade        = database$shade - mean(database$shade)
database$shorttime    = database$shorttime - mean(database$shorttime)
database$shortestdist = database$shortestdist - mean(database$shortestdist)
database$lesstraffic  = database$lesstraffic - mean(database$lesstraffic)
database$slowvehspeed = database$slowvehspeed - mean(database$slowvehspeed)
database$resstreet    = database$resstreet - mean(database$resstreet)

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc_yes = 0,
              b_age =0,
              b_gender =0,
              b_income = 0,
              b_hhead =0,
              b_educ =0,
              b_bikeown = 0,
              b_sector = 0,
              lambda_safebiking = 0,
              lambda_environment = 0,
              lambda_easeoftravel = 0,
              lambda_lowtraffic = 0,
              gamma_bikelane = 1,
              gamma_pavement = 1,
              gamma_safety = 1,
              gamma_scenery = 1,
              gamma_shade = 1,
              gamma_shorttime = 1,
              gamma_shortestdist = 1,
              gamma_lesstraffic = 1,
              gamma_slowvehspeed = 1,
              gamma_resstreet = 1,
              sigma_bikelane = 1,
              sigma_pavement = 1,
              sigma_safety = 1,
              sigma_scenery = 1,
              sigma_shade = 1,
              sigma_shorttime = 1,
              sigma_shortestdist = 1,
              sigma_lesstraffic = 1,
              sigma_slowvehspeed = 1,
              sigma_resstreet = 1)

### 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="halton",
  interNDraws=100,
  interUnifDraws=c(),
  interNormDraws=c("eta_safebiking", "eta_environment", "eta_easeoftravel", "eta_lowtraffic"),
  
  intraDrawsType='',
  intraNDraws=0,
  intraUnifDraws=c(),
  intraNormDraws=c()
)

### Create random parameters
apollo_randCoeff=function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["LV_safebiking"]]   = eta_safebiking + gamma_bikelane * bikelane + gamma_pavement * pavement + gamma_safety * safety + b_age * age + b_bikeown * bikeown + b_gender * gender
  randcoeff[["LV_environment"]]  = eta_environment + gamma_scenery * scenery + gamma_shade * shade + b_age * age + b_bikeown * bikeown + b_gender * gender
  randcoeff[["LV_easeoftravel"]] = eta_easeoftravel + gamma_shorttime * shorttime + gamma_shortestdist * shortestdist + b_age * age + b_bikeown * bikeown + b_gender * gender
  randcoeff[["LV_lowtraffic"]]   = eta_lowtraffic + gamma_lesstraffic * lesstraffic + gamma_slowvehspeed * slowvehspeed + gamma_resstreet * resstreet + b_age * age + b_bikeown * bikeown + b_gender * gender
  
  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=bikelane,
                                 xNormal=gamma_bikelane*LV_safebiking,
                                 mu=0,
                                 sigma=sigma_bikelane,
                                 componentName='indic_safebiking_1')
  
  normalDensity_settings2 = list(outcomeNormal=pavement,
                                 xNormal=gamma_pavement*LV_safebiking,
                                 mu=0,
                                 sigma=sigma_pavement,
                                 componentName='indic_safebiking_2')
  
  normalDensity_settings3 = list(outcomeNormal=safety,
                                 xNormal=gamma_safety*LV_safebiking,
                                 mu=0,
                                 sigma=sigma_safety,
                                 componentName='indic_safebiking_3')
  
  normalDensity_settings4 = list(outcomeNormal=scenery,
                                 xNormal=gamma_scenery*LV_environment,
                                 mu=0,
                                 sigma=sigma_scenery,
                                 componentName='indic_environment_1')
  
  normalDensity_settings5 = list(outcomeNormal=shade,
                                 xNormal=gamma_shade*LV_environment,
                                 mu=0,
                                 sigma=sigma_shade,
                                 componentName='indic_environment_2')
  
  normalDensity_settings6 = list(outcomeNormal=shorttime,
                                 xNormal=gamma_shorttime*LV_easeoftravel,
                                 mu=0,
                                 sigma=sigma_shorttime,
                                 componentName='indic_easoftravel_1')
  
  normalDensity_settings7 = list(outcomeNormal=shortestdist,
                                 xNormal=gamma_shortestdist*LV_easeoftravel,
                                 mu=0,
                                 sigma=sigma_shortestdist,
                                 componentName='indic_easoftravel_2')
  
  
  normalDensity_settings8 = list(outcomeNormal=lesstraffic,
                                 xNormal=gamma_lesstraffic*LV_lowtraffic,
                                 mu=0,
                                 sigma=sigma_lesstraffic,
                                 componentName='indic_lowtraffic_1')
  
  normalDensity_settings9 = list(outcomeNormal=slowvehspeed,
                                 xNormal=gamma_slowvehspeed*LV_lowtraffic,
                                 mu=0,
                                 sigma=sigma_slowvehspeed,
                                 componentName='indic_lowtraffic_2')
  
  normalDensity_settings10= list(outcomeNormal=resstreet,
                                 xNormal=gamma_resstreet*LV_lowtraffic,
                                 mu=0,
                                 sigma=sigma_resstreet,
                                 componentName='indic_lowtraffic_3')
  
  P[["indic_safebiking_1"]]   = apollo_normalDensity(normalDensity_settings1, functionality)
  P[["indic_safebiking_2"]]   = apollo_normalDensity(normalDensity_settings2, functionality)
  P[["indic_safebiking_3"]]   = apollo_normalDensity(normalDensity_settings3, functionality)
  P[["indic_environment_1"]]  = apollo_normalDensity(normalDensity_settings4, functionality)
  P[["indic_environment_2"]]  = apollo_normalDensity(normalDensity_settings5, functionality)
  P[["indic_easeoftravel_1"]] = apollo_normalDensity(normalDensity_settings6, functionality)
  P[["indic_easeoftravel_2"]] = apollo_normalDensity(normalDensity_settings7, functionality)
  P[["indic_lowtraffic_1"]]   = apollo_normalDensity(normalDensity_settings8, functionality)
  P[["indic_lowtraffic_2"]]   = apollo_normalDensity(normalDensity_settings9, functionality)
  P[["indic_lowtraffic_3"]]   = apollo_normalDensity(normalDensity_settings10, functionality)
  
  ### Likelihood of choices
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  
  V = list()
  V[['No']] = 0
  V[['Yes']] = asc_yes + b_age * age + b_gender * gender + b_income * income + + b_hhead * hhead + b_educ * educ + b_sector * sector + b_bikeown * bikeown + 
    lambda_safebiking * LV_safebiking + lambda_environment * LV_environment + lambda_easeoftravel * LV_easeoftravel + lambda_lowtraffic * LV_lowtraffic
  
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives = c(No=0,Yes=1),
    avail = 1,
    choiceVar = biketowork,
    V = V
  )
  
  ### Compute probabilities for MNL model component
  P[["biketowork"]] = apollo_mnl(mnl_settings, functionality)
  
  ### Likelihood of the whole model
  P = apollo_combineModels(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)

# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model)