MNL_full.R

MNL_full.R

database <- fread("database_sample.csv")

apollo_control = list(
  modelName = "MNL_full",
  indivID = "case",
  panelData = FALSE,
  mixing = FALSE,
  nCores = 4
)

apollo_control <- apollo_validateControl(database = database, apollo_control = apollo_control)
database <- apollo_validateData(database, apollo_control, silent = FALSE)

modes = c("car", "bicycle", "ebicycle", "escooter", "emoped")
modes_encoded = c(
  "bicycle" = 1,
  "car" = 2,
  "ebicycle" = 3,
  "escooter" = 4,
  "emoped" = 5
)
default = "bicycle"

# mapping of column names and coefficient names
var_map_alternative_specific = c(
  "D" = "length",
  "A" = "ascent_per_100m", 
  "B" = "mean_bikeability"
)

var_map_generic = c(
  "D2C" = "start_d2c",
  "O" = "outbound",
  "W" = "weekend",
  "E_S" = "ends_in_sharenow_area",
  "E_N" = "ends_in_nextbike_area",
  "E_T" = "ends_in_tier_area",
  "H" = "is_holiday",
  "T_0" = "timebucket0004",
  "T_4" = "timebucket0408",
  "T_8" = "timebucket0812",
  "T_12" = "timebucket1216",
  "T_16" = "timebucket1620",
  "T_20" = "timebucket2024",
  "C" = "current_temp",
  "WS" = "current_wind_speed",
  "R" = "current_precipitation_volume",
  "POI_start_alpha" = "transport_start",
  "POI_start_beta" = "shopping_start",
  "POI_start_gamma" = "food_start",
  "POI_start_delta" = "leisure_and_entertainment_start",
  "POI_start_epsilon" = "education_start",
  "POI_start_zeta" = "services_start",
  "POI_start_eta" = "outdoor_start",
  "POI_end_alpha" = "transport_end",
  "POI_end_beta" = "shopping_end",
  "POI_end_gamma" = "food_end",
  "POI_end_delta" = "leisure_and_entertainment_end",
  "POI_end_epsilon" = "education_end",
  "POI_end_zeta" = "services_end",
  "POI_end_eta" = "outdoor_end"
)

var_map = c(var_map_alternative_specific, var_map_generic)

alternative_specific_coef = c(
  c("ASC_car", "ASC_bicycle", "ASC_ebicycle", "ASC_emoped", "ASC_escooter"), 
  do.call(
    paste, 
    c(
      expand.grid(
        "b",
        c(
          "D", "A", "B", "D2C", "O", "W", "H",
          "E_S", "E_N", "E_T",
          "T_0", "T_4", "T_8",
          "T_12", "T_16", "T_20", "C", "WS", "R",
          "POI_start_alpha", "POI_start_beta","POI_start_gamma",  "POI_start_delta",
          "POI_start_epsilon",  "POI_start_zeta",  "POI_start_eta",
          "POI_end_alpha", "POI_end_beta","POI_end_gamma",
          "POI_end_delta",  "POI_end_epsilon",  "POI_end_zeta",  "POI_end_eta"
        ),
        modes
      ), 
     sep = "_")
  )
)
alternative_specific_coef <- alternative_specific_coef[! alternative_specific_coef %in% c(
  "b_SOC_bicycle", do.call(paste, c(expand.grid("b", names(var_map_generic), default), sep="_")))]

# generic_coef = c("b_P")
generic_coef = c()

apollo_beta = c(
  setNames(rep(0, length(alternative_specific_coef)), 
                         alternative_specific_coef),
  setNames(rep(0, length(generic_coef)), generic_coef)
)

apollo_fixed = c(
  "ASC_bicycle"
  , "b_T_0_car", "b_T_0_ebicycle", "b_T_0_emoped", "b_T_0_escooter"
)
  
apollo_inputs <- apollo_validateInputs()

# retain various variables for use in apollo function scope
apollo_inputs$modes <- modes
apollo_inputs$modes_encoded <- modes_encoded
apollo_inputs$var_map_alternative_specific <- var_map_alternative_specific
apollo_inputs$var_map_generic <- var_map_generic
apollo_inputs$var_map <- var_map
apollo_inputs$alternative_specific_coef <- alternative_specific_coef
apollo_inputs$generic_coef <- generic_coef
  
apollo_probabilities = function(apollo_beta, apollo_inputs, functionality="estimate"){
  ### Function initialisation: do not change the following three commands
  ### 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
  Prob = list()

  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  
  V = list()
  V[["car"]] = ASC_car +
    b_D_car * length_car +
    b_A_car * ascent_per_100m_car +
    b_B_car * mean_bikeability_car +
    b_D2C_car * start_d2c +
    b_O_car * outbound +
    b_W_car * weekend +
    b_H_car * is_holiday +
    b_E_S_car * ends_in_sharenow_area +
    b_E_N_car * ends_in_nextbike_area +
    b_E_T_car * ends_in_tier_area +
    b_T_0_car * timebucket0004 +
    b_T_4_car * timebucket0408 +
    b_T_8_car * timebucket0812 +
    b_T_12_car * timebucket1216 +
    b_T_16_car * timebucket1620 +
    b_T_20_car * timebucket2024 +
    b_C_car * current_temp +
    b_WS_car * current_wind_speed +
    b_R_car * current_precipitation_volume +
    b_POI_start_alpha_car * transport_start +
    b_POI_start_beta_car * shopping_start +
    b_POI_start_gamma_car * food_start +
    b_POI_start_delta_car * leisure_and_entertainment_start +
    b_POI_start_epsilon_car * education_start +
    b_POI_start_zeta_car * services_start +
    b_POI_start_eta_car * outdoor_start +
    b_POI_end_alpha_car * transport_end +
    b_POI_end_beta_car * shopping_end +
    b_POI_end_gamma_car * food_end +
    b_POI_end_delta_car * leisure_and_entertainment_end +
    b_POI_end_epsilon_car * education_end +
    b_POI_end_zeta_car * services_end +
    b_POI_end_eta_car * outdoor_end

  V[["bicycle"]] = ASC_bicycle +
    b_D_bicycle * length_bicycle +
    b_A_bicycle * ascent_bicycle +
    b_B_bicycle * mean_bikeability_bicycle

  V[["ebicycle"]] = ASC_ebicycle +
    b_D_ebicycle * length_ebicycle +
    b_A_ebicycle * ascent_per_100m_ebicycle +
    b_B_ebicycle * mean_bikeability_ebicycle +
    b_D2C_ebicycle * start_d2c +
    b_O_ebicycle * outbound +
    b_W_ebicycle * weekend +
    b_H_ebicycle * is_holiday +
    b_E_S_ebicycle * ends_in_sharenow_area +
    b_E_N_ebicycle * ends_in_nextbike_area +
    b_E_T_ebicycle * ends_in_tier_area +
    b_T_0_ebicycle * timebucket0004 +
    b_T_4_ebicycle * timebucket0408 +
    b_T_8_ebicycle * timebucket0812 +
    b_T_12_ebicycle * timebucket1216 +
    b_T_16_ebicycle * timebucket1620 +
    b_T_20_ebicycle * timebucket2024 +
    b_C_ebicycle * current_temp +
    b_WS_ebicycle * current_wind_speed +
    b_R_ebicycle * current_precipitation_volume +
    b_POI_start_alpha_ebicycle * transport_start +
    b_POI_start_beta_ebicycle * shopping_start +
    b_POI_start_gamma_ebicycle * food_start +
    b_POI_start_delta_ebicycle * leisure_and_entertainment_start +
    b_POI_start_epsilon_ebicycle * education_start +
    b_POI_start_zeta_ebicycle * services_start +
    b_POI_start_eta_ebicycle * outdoor_start +
    b_POI_end_alpha_ebicycle * transport_end +
    b_POI_end_beta_ebicycle * shopping_end +
    b_POI_end_gamma_ebicycle * food_end +
    b_POI_end_delta_ebicycle * leisure_and_entertainment_end +
    b_POI_end_epsilon_ebicycle * education_end +
    b_POI_end_zeta_ebicycle * services_end +
    b_POI_end_eta_ebicycle * outdoor_end

  V[["escooter"]] = ASC_escooter +
    b_D_escooter * length_escooter +
    b_A_escooter * ascent_per_100m_escooter +
    b_B_escooter * mean_bikeability_escooter +
    b_D2C_escooter * start_d2c +
    b_O_escooter * outbound +
    b_W_escooter * weekend +
    b_H_escooter * is_holiday +
    b_E_S_escooter * ends_in_sharenow_area +
    b_E_N_escooter * ends_in_nextbike_area +
    b_E_T_escooter * ends_in_tier_area +
    b_T_0_escooter * timebucket0004 +
    b_T_4_escooter * timebucket0408 +
    b_T_8_escooter * timebucket0812 +
    b_T_12_escooter * timebucket1216 +
    b_T_16_escooter * timebucket1620 +
    b_T_20_escooter * timebucket2024 +
    b_C_escooter * current_temp +
    b_WS_escooter * current_wind_speed +
    b_R_escooter * current_precipitation_volume +
    b_POI_start_alpha_escooter * transport_start +
    b_POI_start_beta_escooter * shopping_start +
    b_POI_start_gamma_escooter * food_start +
    b_POI_start_delta_escooter * leisure_and_entertainment_start +
    b_POI_start_epsilon_escooter * education_start +
    b_POI_start_zeta_escooter * services_start +
    b_POI_start_eta_escooter * outdoor_start +
    b_POI_end_alpha_escooter * transport_end +
    b_POI_end_beta_escooter * shopping_end +
    b_POI_end_gamma_escooter * food_end +
    b_POI_end_delta_escooter * leisure_and_entertainment_end +
    b_POI_end_epsilon_escooter * education_end +
    b_POI_end_zeta_escooter * services_end +
    b_POI_end_eta_escooter * outdoor_end

  V[["emoped"]] = ASC_emoped +
    b_D_emoped * length_emoped +
    b_A_emoped * ascent_per_100m_emoped +
    b_B_emoped * mean_bikeability_emoped +
    b_D2C_emoped * start_d2c +
    b_O_emoped * outbound +
    b_W_emoped * weekend +
    b_H_emoped * is_holiday +
    b_E_S_emoped * ends_in_sharenow_area +
    b_E_N_emoped * ends_in_nextbike_area +
    b_E_T_emoped * ends_in_tier_area +
    b_T_0_emoped * timebucket0004 +
    b_T_4_emoped * timebucket0408 +
    b_T_8_emoped * timebucket0812 +
    b_T_12_emoped * timebucket1216 +
    b_T_16_emoped * timebucket1620 +
    b_T_20_emoped * timebucket2024 +
    b_C_emoped * current_temp +
    b_WS_emoped * current_wind_speed +
    b_R_emoped * current_precipitation_volume +
    b_POI_start_alpha_emoped * transport_start +
    b_POI_start_beta_emoped * shopping_start +
    b_POI_start_gamma_emoped * food_start +
    b_POI_start_delta_emoped * leisure_and_entertainment_start +
    b_POI_start_epsilon_emoped * education_start +
    b_POI_start_zeta_emoped * services_start +
    b_POI_start_eta_emoped * outdoor_start +
    b_POI_end_alpha_emoped * transport_end +
    b_POI_end_beta_emoped * shopping_end +
    b_POI_end_gamma_emoped * food_end +
    b_POI_end_delta_emoped * leisure_and_entertainment_end +
    b_POI_end_epsilon_emoped * education_end +
    b_POI_end_zeta_emoped * services_end +
    b_POI_end_eta_emoped * outdoor_end
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = apollo_inputs$modes_encoded, 
    avail         = list(
      car = car_available,
      bicycle = bicycle_available,
      ebicycle = ebicycle_available,
      escooter = escooter_available,
      emoped = emoped_available
    ), 
    choiceVar     = vehicle_type,
    utilities     = V
  )
  
  ### Compute probabilities using MNL model
  Prob[["model"]] = apollo_mnl(mnl_settings, functionality)
  
  ### Prepare and return outputs of function
  Prob = apollo_prepareProb(Prob, apollo_inputs, functionality)
  return(Prob)
}

# ################################################################# #
#### MODEL ESTIMATION                                            ####
# ################################################################# #

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, 
                        apollo_inputs, list(validateGrad = FALSE, estimationRoutine = "BGW"))

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

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

apollo_modelOutput(model)