Important: Read this before posting to this forum

  1. This forum is for questions related to the use of Apollo. We will answer some general choice modelling questions too, where appropriate, and time permitting. We cannot answer questions about how to estimate choice models with other software packages.
  2. There is a very detailed manual for Apollo available at http://www.ApolloChoiceModelling.com/manual.html. This contains detailed descriptions of the various Apollo functions, and numerous examples are available at http://www.ApolloChoiceModelling.com/examples.html. In addition, help files are available for all functions, using e.g. ?apollo_mnl
  3. Before asking a question on the forum, users are kindly requested to follow these steps:
    1. Check that the same issue has not already been addressed in the forum - there is a search tool.
    2. Ensure that the correct syntax has been used. For any function, detailed instructions are available directly in Apollo, e.g. by using ?apollo_mnl for apollo_mnl
    3. Check the frequently asked questions section on the Apollo website, which discusses some common issues/failures. Please see http://www.apollochoicemodelling.com/faq.html
    4. Make sure that R is using the latest official release of Apollo.
  4. If the above steps do not resolve the issue, then users should follow these steps when posting a question:
    1. provide full details on the issue, including the entire code and output, including any error messages
    2. posts will not immediately appear on the forum, but will be checked by a moderator first. This may take a day or two at busy times. There is no need to submit the post multiple times.

Large (M)MNL - Analytical gradients could not be calculated for all components

Ask questions about errors you encouunter. Please make sure to include full details about your model specifications, and ideally your model file.
Post Reply
janikm
Posts: 1
Joined: 24 Aug 2023, 15:16

Large (M)MNL - Analytical gradients could not be calculated for all components

Post by janikm »

Dear Apollo team,

I need to estimate a MXL/MMNL, which I want to preface by estimating the respective MNL, both to compare likelihood and to solve errors from the bottom up, which worked out up until now. Before going to the MXL, I am now at a stage where I am unsure how to continue, as in practice the MNL model can be estimated, but multiple things are weird:

1) I tried to iteratively construct utility functions, but this leads to a warning:
Code:

Code: Select all

database <- fread("database_sample.csv")

apollo_control = list(
  modelName = "MNL_full",
  indivID = "case", # I cannot distinguish individuals, i.e. there is only one case per individual
  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")
default = "bicycle"

# mapping of column names and coefficient names
var_map_alternative_specific = c(
  "SOC" = "start_SOC",
  # "P" = "expected_price_per_min",
  # "V" = "is_closeby",
  "D" = "length",
  "A" = "ascent", 
  "B" = "mean_bikeability"
)

var_map_generic = c(
  "D2C" = "start_d2c",
  "O" = "outbound",
  "W" = "weekend",
  # "H" = "is_holiday", # constant in sample
  "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(
          "SOC", "D", #"V",
          "A", "B", "D2C", "O", "W",
          "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 = "_")
  )
)

# remove bicycle SOC and generic coefficients for the default alternative (instead of fixing them)
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") # don't use price for now
generic_coef = c()

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

# fix default ASC and one level of categorical variable
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
  
  # use below for debugging generated utility functions
  # V_str = list()
  # for (m in apollo_inputs$modes){
  #   V_str[[m]] = paste("ASC", m, sep="_")
  #   for (v in apollo_inputs$alternative_specific_coef){
  #     if (!grepl(paste0("_", m), v, fixed=TRUE)) next # skip other than m coefficients
  #     if (grepl("ASC", v, fixed=TRUE)) next # skip ASC
  #     if (v == "SOC" & m == "bicycle") next # don't include SOC in bicycle utility function
  #     if (v %in% names(apollo_inputs$var_map_generic) & m == "bicycle") next # default
  #     v_raw = sub(paste0("_",m), "", sub("b_", "", v))
  #     as_coef = v_raw %in% names(apollo_inputs$var_map_alternative_specific)
  #     V_str[[m]] = paste(V_str[[m]], "+", v, "*", ifelse(as_coef, paste(apollo_inputs$var_map[v_raw], m, sep="_"), apollo_inputs$var_map[v_raw]))
  #   }
  # 
  #   for (v in apollo_inputs$generic_coef){
  #     if (v %in% names(apollo_inputs$var_map_generic) & m == "bicycle") next # default
  #     v_raw = sub("b_", "", v)
  #     as_coef = v_raw %in% names(apollo_inputs$var_map_alternative_specific)
  #     V_str[[m]] = paste(V_str[[m]], "+", v, "*", ifelse(as_coef, paste(apollo_inputs$var_map[v_raw], m, sep="_"), apollo_inputs$var_map[v_raw]))
  #   }
  # }
  # print(V_str)

  V = list()
  for (m in apollo_inputs$modes){
    V[[m]] = get(paste("ASC", m, sep="_"))
    for (v in apollo_inputs$alternative_specific_coef){
      if (!grepl(paste0("_", m), v, fixed=TRUE)) next # skip other than mode coefficients
      if (grepl("ASC", v, fixed=TRUE)) next # skip ASC, already there
      if (v == "SOC" & m == "bicycle") next # don't include SOC in bicycle utility function
      if (v %in% names(apollo_inputs$var_map_generic) & m == "bicycle") next # default
      v_raw = sub(paste0("_",m), "", sub("b_", "", v))
      as_coef = v_raw %in% names(apollo_inputs$var_map_alternative_specific)
      V[[m]] = V[[m]] + get(v) * get(ifelse(as_coef, paste(apollo_inputs$var_map[v_raw], m, sep="_"), apollo_inputs$var_map[v_raw]))
    }

    for (v in apollo_inputs$generic_coef){
      if (v %in% apollo_inputs$var_map_generic & m == "bicycle") next # default
      v_raw = sub("b_", "", v)
      as_coef = v_raw %in% names(apollo_inputs$var_map_alternative_specific)
      V[[m]] = V[[m]]+ get(v) * get(ifelse(as_coef, paste(apollo_inputs$var_map[v_raw], m, sep="_"), apollo_inputs$var_map[v_raw]))
    }
  }
  
  ### 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(estimationRoutine="BFGS"))
I can confirm, using the debug functionality commented out in the code above and a manual call to apollo_probabilities, that the utility functions are generated as wanted, but something weird is happening when observing the same debug print throughout the apollo_estimate call:

Code: Select all

Preparing user-defined functions.
$car
[1] "ASC_car + b_SOC_car * start_SOC_car + b_D_car * length_car + b_A_car * ascent_car + b_B_car * mean_bikeability_car + b_D2C_car * start_d2c + b_O_car * outbound + b_W_car * weekend + 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"

$bicycle
[1] "ASC_bicycle + b_D_bicycle * length_bicycle + b_A_bicycle * ascent_bicycle + b_B_bicycle * mean_bikeability_bicycle"

$ebicycle
[1] "ASC_ebicycle + b_SOC_ebicycle * start_SOC_ebicycle + b_D_ebicycle * length_ebicycle + b_A_ebicycle * ascent_ebicycle + b_B_ebicycle * mean_bikeability_ebicycle + b_D2C_ebicycle * start_d2c + b_O_ebicycle * outbound + b_W_ebicycle * weekend + 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"

$escooter
[1] "ASC_escooter + b_SOC_escooter * start_SOC_escooter + b_D_escooter * length_escooter + b_A_escooter * ascent_escooter + b_B_escooter * mean_bikeability_escooter + b_D2C_escooter * start_d2c + b_O_escooter * outbound + b_W_escooter * weekend + 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"

$emoped
[1] "ASC_emoped + b_SOC_emoped * start_SOC_emoped + b_D_emoped * length_emoped + b_A_emoped * ascent_emoped + b_B_emoped * mean_bikeability_emoped + b_D2C_emoped * start_d2c + b_O_emoped * outbound + b_W_emoped * weekend + 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"

...

$car
[1] "ASC_car + b_POI_end_eta_car * outdoor_end + b_POI_end_eta_car * NA + b_POI_end_eta_car * NA"

$bicycle
[1] "ASC_bicycle + b_B_bicycle * mean_bikeability_bicycle + b_B_bicycle * NA + b_B_bicycle * NA"

$ebicycle
[1] "ASC_ebicycle + b_POI_end_eta_ebicycle * outdoor_end + b_POI_end_eta_ebicycle * NA + b_POI_end_eta_ebicycle * NA"

$escooter
[1] "ASC_escooter + b_POI_end_eta_escooter * outdoor_end + b_POI_end_eta_escooter * NA + b_POI_end_eta_escooter * NA"

$emoped
[1] "ASC_emoped + b_POI_end_eta_emoped * outdoor_end + b_POI_end_eta_emoped * NA + b_POI_end_eta_emoped * NA"

WARNING: The pre-processing of 'apollo_probabilities' failed after
  inserting functions. Your model may still run, but this
  indicates a potential problem. Please contact the developers
  for assistance! 
Right before the warning there are NAs in what seems to be parts of the utility functions as generated by Apollo.
This is where I canceled the process, and went ahead to hardcode utility functions using the outputs generated above:
This replaces the utility function part in apollo_probabilities:

Code: Select all

  V = list()
  V[["car"]] = ASC_car +
    b_SOC_car * start_SOC_car +
    b_D_car * length_car +
    b_A_car * ascent_car +
    b_B_car * mean_bikeability_car +
    b_D2C_car * start_d2c +
    b_O_car * outbound +
    b_W_car * weekend +
    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_SOC_ebicycle * start_SOC_ebicycle +
    b_D_ebicycle * length_ebicycle +
    b_A_ebicycle * ascent_ebicycle +
    b_B_ebicycle * mean_bikeability_ebicycle +
    b_D2C_ebicycle * start_d2c +
    b_O_ebicycle * outbound +
    b_W_ebicycle * weekend +
    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_SOC_escooter * start_SOC_escooter +
    b_D_escooter * length_escooter +
    b_A_escooter * ascent_escooter +
    b_B_escooter * mean_bikeability_escooter +
    b_D2C_escooter * start_d2c +
    b_O_escooter * outbound +
    b_W_escooter * weekend +
    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_SOC_emoped * start_SOC_emoped +
    b_D_emoped * length_emoped +
    b_A_emoped * ascent_emoped +
    b_B_emoped * mean_bikeability_emoped +
    b_D2C_emoped * start_d2c +
    b_O_emoped * outbound +
    b_W_emoped * weekend +
    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
2) Now, I don't get this warning anymore, but what remains is the following - a warning about not being able to compute analytical gradients for some components, and a singular hessian:

Code: Select all

Preparing user-defined functions.

Testing likelihood function...

Overview of choices for MNL model component :
                                 bicycle      car ebicycle emoped
Times available                  10000.0 10000.00  1.0e+04  1e+04
Times chosen                      3460.0  6347.00  2.1e+01  5e+00
Percentage chosen overall           34.6    63.47  2.1e-01  5e-02
Percentage chosen when available    34.6    63.47  2.1e-01  5e-02
                                 escooter
Times available                  10000.00
Times chosen                       167.00
Percentage chosen overall            1.67
Percentage chosen when available     1.67


Pre-processing likelihood function...
Creating cluster...
Preparing workers for multithreading...
Analytical gradient is different to numerical one. Numerical
  gradients will be used.
Analytical gradients could not be calculated for all components,
  numerical gradients will be used.

Testing influence of parameters...........................................................................................................................
Starting main estimation

BGW is using FD derivatives for model Jacobian. (Caller did not provide derivatives.)


Iterates will be written to: 
 /home/jmuires/source/modechoice/MNL_full_iterations.csv    it    nf     F            RELDF    PRELDF    RELDX    MODEL stppar
     0     1 1.609437912e+04
     1     4 1.089770578e+04 3.229e-01 3.044e-01 1.00e+00   G   7.78e-01
     2     5 7.669113839e+03 2.963e-01 3.799e-01 4.39e-01   G   4.93e+00
     3     6 6.786696142e+03 1.151e-01 1.343e-01 5.69e-01   G   2.84e+00
     4     7 5.961581698e+03 1.216e-01 1.457e-01 5.42e-01   G   5.83e-01
     5     8 5.642541748e+03 5.352e-02 5.860e-02 2.49e-01   G   2.23e-01
     6    10 5.400676085e+03 4.286e-02 4.950e-02 4.22e-01   G   2.99e-02
     7    12 5.359884706e+03 7.553e-03 6.111e-03 3.14e-02   G   2.51e-01
     8    14 5.325231630e+03 6.465e-03 7.497e-03 1.41e-01   G   4.64e-02
     9    16 5.318041150e+03 1.350e-03 1.798e-03 8.60e-03   G   1.24e+00
    10    18 5.313333668e+03 8.852e-04 8.546e-04 1.37e-02   S   1.26e-01
...
   117   146 5.276564207e+03 6.204e-12 6.689e-11 5.12e-06   S   6.49e-02
   118   147 5.276564207e+03 3.040e-11 4.788e-11 1.90e-06   S   1.73e-01
   119   148 5.276564207e+03 1.430e-11 1.693e-11 4.32e-06   S   7.04e-02

***** Relative function convergence *****
Additional convergence test using scaled estimation. Parameters
  will be scaled by their current estimates and additional
  iterations will be performed.

BGW is using FD derivatives for model Jacobian. (Caller did not provide derivatives.)


Iterates will be appended to: 
 /home/jmuires/source/modechoice/MNL_full_iterations.csv    it    nf     F            RELDF    PRELDF    RELDX    MODEL stppar
     0     1 5.276564207e+03
     1     7 5.276564207e+03 9.808e-13 2.158e-12 1.21e-05   G   3.95e+00
     2    15 5.276564207e+03 9.790e-14 5.785e-13 3.39e-06  G-S  1.18e+01
     3    20 5.276564207e+03 8.101e-15 1.669e-16 1.51e-09   G   1.76e+04
     4    26 5.276564207e+03 2.068e-15 3.276e-21 3.00e-14   G   8.89e+08

***** False convergence *****
The estimation of the scaled model failed, and the unscaled
  version will be returned instead.

Estimated parameters:
                                Estimate
ASC_car                         2.239579
ASC_bicycle                     0.000000
ASC_ebicycle                   -2.782581
ASC_emoped                     -2.110107
ASC_escooter                    6.039018
b_SOC_car                      -0.017630
b_D_car                       -4.384e-05
b_A_car                         0.004233
b_B_car                        -1.206232
b_D2C_car                     1.4587e-04
b_O_car                         0.147830
b_W_car                        -0.129533
b_T_0_car                       0.000000
b_T_4_car                       0.678080
b_T_8_car                       0.918854
b_T_12_car                      0.824404
b_T_16_car                      0.902038
b_T_20_car                      0.555138
b_C_car                        -0.044104
b_WS_car                        0.039057
b_R_car                         0.080331
b_POI_start_alpha_car          -0.543449
b_POI_start_beta_car            1.362755
b_POI_start_gamma_car          -1.147940
b_POI_start_delta_car          -0.052960
b_POI_start_epsilon_car        -1.819648
b_POI_start_zeta_car            0.571766
b_POI_start_eta_car            -1.187990
b_POI_end_alpha_car            -0.569859
b_POI_end_beta_car              1.066813
b_POI_end_gamma_car            -1.262325
b_POI_end_delta_car            -0.504386
b_POI_end_epsilon_car          -1.617586
b_POI_end_zeta_car              0.980248
b_POI_end_eta_car              -0.898540
b_D_bicycle                  -3.9293e-04
b_A_bicycle                     0.004089
b_B_bicycle                    -0.502259
b_SOC_ebicycle                 -0.018573
b_D_ebicycle                 -2.7892e-04
b_A_ebicycle                    0.003795
b_B_ebicycle                    0.793691
b_D2C_ebicycle                5.2105e-04
b_O_ebicycle                    0.511662
b_W_ebicycle                   -0.905533
b_T_0_ebicycle                  0.000000
b_T_4_ebicycle                 -2.019232
b_T_8_ebicycle                 -1.654777
b_T_12_ebicycle                -1.839972
b_T_16_ebicycle                -2.070412
b_T_20_ebicycle                -0.845246
b_C_ebicycle                   -0.028911
b_WS_ebicycle                   0.197240
b_R_ebicycle                   -1.577646
b_POI_start_alpha_ebicycle     -1.264917
b_POI_start_beta_ebicycle      -3.426741
b_POI_start_gamma_ebicycle     -1.354868
b_POI_start_delta_ebicycle      3.031877
b_POI_start_epsilon_ebicycle   -3.675478
b_POI_start_zeta_ebicycle      -2.896654
b_POI_start_eta_ebicycle        4.386000
b_POI_end_alpha_ebicycle       -1.105115
b_POI_end_beta_ebicycle        -0.943038
b_POI_end_gamma_ebicycle       -0.608778
b_POI_end_delta_ebicycle        0.681714
b_POI_end_epsilon_ebicycle     -1.759777
b_POI_end_zeta_ebicycle        -1.216601
b_POI_end_eta_ebicycle          2.233369
b_SOC_escooter                 -0.063476
b_D_escooter                 -5.3448e-04
b_A_escooter                   -0.012614
b_B_escooter                    0.505693
b_D2C_escooter                1.9038e-04
b_O_escooter                   -0.034644
b_W_escooter                   -0.507325
b_T_0_escooter                  0.000000
b_T_4_escooter                 -1.546724
b_T_8_escooter                 -1.430709
b_T_12_escooter                -2.008773
b_T_16_escooter                -2.823806
b_T_20_escooter                -1.479479
b_C_escooter                  -8.498e-05
b_WS_escooter                  -0.018488
b_R_escooter                 -138.492939
b_POI_start_alpha_escooter     -1.146293
b_POI_start_beta_escooter      -1.257266
b_POI_start_gamma_escooter     -0.142349
b_POI_start_delta_escooter      1.699514
b_POI_start_epsilon_escooter   -1.952530
b_POI_start_zeta_escooter      -1.570653
b_POI_start_eta_escooter        2.739583
b_POI_end_alpha_escooter        0.499049
b_POI_end_beta_escooter         0.006354
b_POI_end_gamma_escooter       -0.982700
b_POI_end_delta_escooter       -0.578088
b_POI_end_epsilon_escooter     -0.907645
b_POI_end_zeta_escooter         0.007862
b_POI_end_eta_escooter         -1.154322
b_SOC_emoped                   -0.013492
b_D_emoped                    -8.282e-05
b_A_emoped                      0.007010
b_B_emoped                     -0.063563
b_D2C_emoped                 -1.5095e-04
b_O_emoped                     -0.485589
b_W_emoped                    -15.118535
b_T_0_emoped                    0.000000
b_T_4_emoped                  -16.916359
b_T_8_emoped                   -2.978931
b_T_12_emoped                  -2.943153
b_T_16_emoped                  -3.318909
b_T_20_emoped                 -17.055849
b_C_emoped                      0.109977
b_WS_emoped                     0.330498
b_R_emoped                      0.009859
b_POI_start_alpha_emoped        3.939417
b_POI_start_beta_emoped        -1.097240
b_POI_start_gamma_emoped       -3.233876
b_POI_start_delta_emoped       -9.084779
b_POI_start_epsilon_emoped     -2.095627
b_POI_start_zeta_emoped        -1.036817
b_POI_start_eta_emoped          3.084146
b_POI_end_alpha_emoped         -5.515382
b_POI_end_beta_emoped           5.072686
b_POI_end_gamma_emoped         -8.525597
b_POI_end_delta_emoped          0.754469
b_POI_end_epsilon_emoped        0.404757
b_POI_end_zeta_emoped           0.987306
b_POI_end_eta_emoped           -1.418269

Final LL: -5276.5642

Calculating log-likelihood at equal shares (LL(0)) for applicable
  models...
Calculating log-likelihood at observed shares from estimation
  data (LL(c)) for applicable models...
Calculating LL of each model component...
Calculating other model fit measures
INFORMATION: Your model took more than 10 minutes to estimate, so it was saved
  to file /home/jmuires/source/modechoice/MNL_full_model.rds
  before calculating its covariance matrix. If calculation of the
  covariance matrix fails or is stopped before finishing, you can
  load the model up to this point using apollo_loadModel. You may
  also want to inspect the approximate BHHH standard errors shown
  above to determine whether you wish to continue this process. 
Computing covariance matrix using numerical methods (numDeriv).
 0%....25%....50%....75%....100%
WARNING: Singular Hessian, cannot calculate s.e. 
Could not write hessian to a file.
Negative definite Hessian with maximum eigenvalue: 0
Computing score matrix...

Your model was estimated using the BGW algorithm. Please
  acknowledge this by citing Bunch et al. (1993) - DOI
  10.1145/151271.151279
  
 > apollo_modelOutput(model)
Model run by jmuires using Apollo 0.3.0 on R 4.3.1 for Linux.
Please acknowledge the use of Apollo by citing Hess & Palma (2019)
  DOI 10.1016/j.jocm.2019.100170
  www.ApolloChoiceModelling.com

Model name                                  : MNL_full
Model description                           : No model description provided in apollo_control
Model run at                                : 2023-08-30 17:41:51.490785
Estimation method                           : bgw
Model diagnosis                             : Relative function convergence
Optimisation diagnosis                      : Maximum found
     hessian properties                     : Negative definite
     maximum eigenvalue                     : 0
     reciprocal of condition number         : 1.93177e-20
Number of individuals                       : 10000
Number of rows in database                  : 10000
Number of modelled outcomes                 : 10000

Number of cores used                        :  4 
Model without mixing

LL(start)                                   : -16094.38
LL at equal shares, LL(0)                   : -16094.38
LL at observed shares, LL(C)                : -7408.43
LL(final)                                   : -5276.56
Rho-squared vs equal shares                  :  0.6721 
Adj.Rho-squared vs equal shares              :  0.6645 
Rho-squared vs observed shares               :  0.2878 
Adj.Rho-squared vs observed shares           :  0.2717 
AIC                                         :  10799.13 
BIC                                         :  11686 

Estimated parameters                        : 123
Time taken (hh:mm:ss)                       :  02:20:10.36 
     pre-estimation                         :  00:02:11.69 
     estimation                             :  00:30:38.79 
     post-estimation                        :  01:47:19.88 
Iterations                                  :  119  

Unconstrained optimisation.

Estimates:
                                Estimate        s.e.   t.rat.(0)
ASC_car                         2.239579          NA          NA
ASC_bicycle                     0.000000          NA          NA
ASC_ebicycle                   -2.782581          NA          NA
ASC_emoped                     -2.110107          NA          NA
ASC_escooter                    6.039018          NA          NA
b_SOC_car                      -0.017630          NA          NA
b_D_car                       -4.384e-05          NA          NA
b_A_car                         0.004233          NA          NA
b_B_car                        -1.206232          NA          NA
b_D2C_car                     1.4587e-04          NA          NA
b_O_car                         0.147830          NA          NA
b_W_car                        -0.129533          NA          NA
b_T_0_car                       0.000000          NA          NA
b_T_4_car                       0.678080          NA          NA
b_T_8_car                       0.918854          NA          NA
b_T_12_car                      0.824404          NA          NA
b_T_16_car                      0.902038          NA          NA
b_T_20_car                      0.555138          NA          NA
b_C_car                        -0.044104          NA          NA
b_WS_car                        0.039057          NA          NA
b_R_car                         0.080331          NA          NA
b_POI_start_alpha_car          -0.543449          NA          NA
b_POI_start_beta_car            1.362755          NA          NA
b_POI_start_gamma_car          -1.147940          NA          NA
b_POI_start_delta_car          -0.052960          NA          NA
b_POI_start_epsilon_car        -1.819648          NA          NA
b_POI_start_zeta_car            0.571766          NA          NA
b_POI_start_eta_car            -1.187990          NA          NA
b_POI_end_alpha_car            -0.569859          NA          NA
b_POI_end_beta_car              1.066813          NA          NA
b_POI_end_gamma_car            -1.262325          NA          NA
b_POI_end_delta_car            -0.504386          NA          NA
b_POI_end_epsilon_car          -1.617586          NA          NA
b_POI_end_zeta_car              0.980248          NA          NA
b_POI_end_eta_car              -0.898540          NA          NA
b_D_bicycle                  -3.9293e-04          NA          NA
b_A_bicycle                     0.004089          NA          NA
b_B_bicycle                    -0.502259          NA          NA
b_SOC_ebicycle                 -0.018573          NA          NA
b_D_ebicycle                 -2.7892e-04          NA          NA
b_A_ebicycle                    0.003795          NA          NA
b_B_ebicycle                    0.793691          NA          NA
b_D2C_ebicycle                5.2105e-04          NA          NA
b_O_ebicycle                    0.511662          NA          NA
b_W_ebicycle                   -0.905533          NA          NA
b_T_0_ebicycle                  0.000000          NA          NA
b_T_4_ebicycle                 -2.019232          NA          NA
b_T_8_ebicycle                 -1.654777          NA          NA
b_T_12_ebicycle                -1.839972          NA          NA
b_T_16_ebicycle                -2.070412          NA          NA
b_T_20_ebicycle                -0.845246          NA          NA
b_C_ebicycle                   -0.028911          NA          NA
b_WS_ebicycle                   0.197240          NA          NA
b_R_ebicycle                   -1.577646          NA          NA
b_POI_start_alpha_ebicycle     -1.264917          NA          NA
b_POI_start_beta_ebicycle      -3.426741          NA          NA
b_POI_start_gamma_ebicycle     -1.354868          NA          NA
b_POI_start_delta_ebicycle      3.031877          NA          NA
b_POI_start_epsilon_ebicycle   -3.675478          NA          NA
b_POI_start_zeta_ebicycle      -2.896654          NA          NA
b_POI_start_eta_ebicycle        4.386000          NA          NA
b_POI_end_alpha_ebicycle       -1.105115          NA          NA
b_POI_end_beta_ebicycle        -0.943038          NA          NA
b_POI_end_gamma_ebicycle       -0.608778          NA          NA
b_POI_end_delta_ebicycle        0.681714          NA          NA
b_POI_end_epsilon_ebicycle     -1.759777          NA          NA
b_POI_end_zeta_ebicycle        -1.216601          NA          NA
b_POI_end_eta_ebicycle          2.233369          NA          NA
b_SOC_escooter                 -0.063476          NA          NA
b_D_escooter                 -5.3448e-04          NA          NA
b_A_escooter                   -0.012614          NA          NA
b_B_escooter                    0.505693          NA          NA
b_D2C_escooter                1.9038e-04          NA          NA
b_O_escooter                   -0.034644          NA          NA
b_W_escooter                   -0.507325          NA          NA
b_T_0_escooter                  0.000000          NA          NA
b_T_4_escooter                 -1.546724          NA          NA
b_T_8_escooter                 -1.430709          NA          NA
b_T_12_escooter                -2.008773          NA          NA
b_T_16_escooter                -2.823806          NA          NA
b_T_20_escooter                -1.479479          NA          NA
b_C_escooter                  -8.498e-05          NA          NA
b_WS_escooter                  -0.018488          NA          NA
b_R_escooter                 -138.492939          NA          NA
b_POI_start_alpha_escooter     -1.146293          NA          NA
b_POI_start_beta_escooter      -1.257266          NA          NA
b_POI_start_gamma_escooter     -0.142349          NA          NA
b_POI_start_delta_escooter      1.699514          NA          NA
b_POI_start_epsilon_escooter   -1.952530          NA          NA
b_POI_start_zeta_escooter      -1.570653          NA          NA
b_POI_start_eta_escooter        2.739583          NA          NA
b_POI_end_alpha_escooter        0.499049          NA          NA
b_POI_end_beta_escooter         0.006354          NA          NA
b_POI_end_gamma_escooter       -0.982700          NA          NA
b_POI_end_delta_escooter       -0.578088          NA          NA
b_POI_end_epsilon_escooter     -0.907645          NA          NA
b_POI_end_zeta_escooter         0.007862          NA          NA
b_POI_end_eta_escooter         -1.154322          NA          NA
b_SOC_emoped                   -0.013492          NA          NA
b_D_emoped                    -8.282e-05          NA          NA
b_A_emoped                      0.007010          NA          NA
b_B_emoped                     -0.063563          NA          NA
b_D2C_emoped                 -1.5095e-04          NA          NA
b_O_emoped                     -0.485589          NA          NA
b_W_emoped                    -15.118535          NA          NA
b_T_0_emoped                    0.000000          NA          NA
b_T_4_emoped                  -16.916359          NA          NA
b_T_8_emoped                   -2.978931          NA          NA
b_T_12_emoped                  -2.943153          NA          NA
b_T_16_emoped                  -3.318909          NA          NA
b_T_20_emoped                 -17.055849          NA          NA
b_C_emoped                      0.109977          NA          NA
b_WS_emoped                     0.330498          NA          NA
b_R_emoped                      0.009859          NA          NA
b_POI_start_alpha_emoped        3.939417          NA          NA
b_POI_start_beta_emoped        -1.097240          NA          NA
b_POI_start_gamma_emoped       -3.233876          NA          NA
b_POI_start_delta_emoped       -9.084779          NA          NA
b_POI_start_epsilon_emoped     -2.095627          NA          NA
b_POI_start_zeta_emoped        -1.036817          NA          NA
b_POI_start_eta_emoped          3.084146          NA          NA
b_POI_end_alpha_emoped         -5.515382          NA          NA
b_POI_end_beta_emoped           5.072686          NA          NA
b_POI_end_gamma_emoped         -8.525597          NA          NA
b_POI_end_delta_emoped          0.754469          NA          NA
b_POI_end_epsilon_emoped        0.404757          NA          NA
b_POI_end_zeta_emoped           0.987306          NA          NA
b_POI_end_eta_emoped           -1.418269          NA          NA
Note, that this is the latest version of Apollo (0.3.0, installed through CRAN) which seems to implement some different estimation routine (BGW?) than previously. This is why I additionally tried the same with BFGS, with similar results (no analytical gradient, singular hessian) except for the fact this time it indicates a potential saddle point convergence.

Code: Select all

WARNING: Singular Hessian, cannot calculate s.e. 
Could not write hessian to a file.
WARNING: Some eigenvalues of the Hessian are positive, indicating
  convergence to a saddle point! 
  
The sample from the dataset used above can be downloaded from here: https://gist.github.com/janikmu/6c5323e ... 1a084dc2d9
Note that choices are heavily unbalanced, but the full dataset has a very large n (revealed preference), so this should be okay in the long run.
Any help is very appreciated! :-)
janikm
Posts: 1
Joined: 24 Aug 2023, 15:16

Re: Large (M)MNL - Analytical gradients could not be calculated for all components

Post by janikm »

Thanks for the quick and helpful answer!
dpalma wrote: 31 Aug 2023, 20:23 The problem seems to be on the definition of mnl_settings$alternatives = apollo_inputs$modes_encoded. In the code you sent modes_encoded was not defined, so I created it as modes_encoded = setNames(1:5, modes). That allowed for the analytical gradient to be constructed.
I'm sorry, I forgot to include this in the code snippet, I recoded choice factors like this to retain the correct levels before sampling and exporting the sample database:
# encode factors to numeric, but retain encoding
modes_encoded <- setNames(1:5, levels(database$vehicle_type))
database$vehicle_type <- as.numeric(database$vehicle_type)
dpalma wrote: 31 Aug 2023, 20:23 However, its value was different to the value of the numerical gradient. It was quite close, but above the threshold defined by Apollo. So I set estimate_settings$validateGrad = FALSE to force Apollo to use the analytical gradient, and then it worked.
Thanks for hinting me at the validateGrad setting, this solved that particular problem for me!
dpalma wrote: 31 Aug 2023, 20:23 Still, the model seems to be ill-defined (at least with the sample database). It seems like some parameters are unbounded (e.g. b_R_emoped), which is consistent with you saying that the data is highly unbalanced. So you might want to look further into the definition of the utilities. On the other hand, maybe this problem will go away when using the full dataset.
I also managed to solve the singular hessian problem by getting rid of one collinearity in the data that went under the radar and employing a more balanced sample of the dataset (uploaded here: https://gist.github.com/janikmu/6c5323e ... 1a084dc2d9). However, now there is an issue with NaNs:
WARNING: Some eigenvalues of the Hessian are complex, indicating that the Hessian is not symmetrical.
Computing score matrix...

Your model was estimated using the BGW algorithm. Please acknowledge this by citing Bunch et al. (1993) - DOI 10.1145/151271.151279
Warning message:
In sqrt(diag(varcov)) : NaNs produced
The code is here: https://gist.github.com/janikmu/7af0251 ... c42d82875d
Last edited by janikm on 06 Sep 2023, 15:23, edited 3 times in total.


Last bumped by janikm on 06 Sep 2023, 12:31.
dpalma
Posts: 190
Joined: 24 Apr 2020, 17:54

Re: Large (M)MNL - Analytical gradients could not be calculated for all components

Post by dpalma »

Hi,

While Apollo can manage to get analytical derivatives for utilities defined via loops, this functionality is mostly for simple loops. When you use conditionals statements (ifs) inside the loop, it is very likely that Apollo won't be able to get the analytical derivative of that. So it's good that you explicitly wrote the utility functions.

The problem seems to be on the definition of mnl_settings$alternatives = apollo_inputs$modes_encoded. In the code you sent modes_encoded was not defined, so I created it as modes_encoded = setNames(1:5, modes). That allowed for the analytical gradient to be constructed. However, its value was different to the value of the numerical gradient. It was quite close, but above the threshold defined by Apollo. So I set estimate_settings$validateGrad = FALSE to force Apollo to use the analytical gradient, and then it worked.

Still, the model seems to be ill-defined (at least with the sample database). It seems like some parameters are unbounded (e.g. b_R_emoped), which is consistent with you saying that the data is highly unbalanced. So you might want to look further into the definition of the utilities. On the other hand, maybe this problem will go away when using the full dataset.

Anyhow, the working code is below.

Best wishes
David

Code: Select all

rm(list=ls())
database <- read.table("database_sample.csv", sep=";", header=TRUE, dec=",")
library(apollo)

apollo_control = list(
  modelName = "MNL_full",
  indivID = "case", # I cannot distinguish individuals, i.e. there is only one case per individual
  panelData = FALSE,
  mixing = FALSE,
  #nCores = 4
  debug = TRUE
)

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")
default = "bicycle"
modes_encoded = setNames(1:5, modes)

# mapping of column names and coefficient names
var_map_alternative_specific = c(
  "SOC" = "start_SOC",
  # "P" = "expected_price_per_min",
  # "V" = "is_closeby",
  "D" = "length",
  "A" = "ascent", 
  "B" = "mean_bikeability"
)

var_map_generic = c(
  "D2C" = "start_d2c",
  "O" = "outbound",
  "W" = "weekend",
  # "H" = "is_holiday", # constant in sample
  "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(
          "SOC", "D", #"V",
          "A", "B", "D2C", "O", "W",
          "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 = "_")
  )
)

# remove bicycle SOC and generic coefficients for the default alternative (instead of fixing them)
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") # don't use price for now
generic_coef = c()

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

# fix default ASC and one level of categorical variable
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_SOC_car * start_SOC_car +
    b_D_car * length_car +
    b_A_car * ascent_car +
    b_B_car * mean_bikeability_car +
    b_D2C_car * start_d2c +
    b_O_car * outbound +
    b_W_car * weekend +
    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_SOC_ebicycle * start_SOC_ebicycle +
    b_D_ebicycle * length_ebicycle +
    b_A_ebicycle * ascent_ebicycle +
    b_B_ebicycle * mean_bikeability_ebicycle +
    b_D2C_ebicycle * start_d2c +
    b_O_ebicycle * outbound +
    b_W_ebicycle * weekend +
    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_SOC_escooter * start_SOC_escooter +
    b_D_escooter * length_escooter +
    b_A_escooter * ascent_escooter +
    b_B_escooter * mean_bikeability_escooter +
    b_D2C_escooter * start_d2c +
    b_O_escooter * outbound +
    b_W_escooter * weekend +
    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_SOC_emoped * start_SOC_emoped +
    b_D_emoped * length_emoped +
    b_A_emoped * ascent_emoped +
    b_B_emoped * mean_bikeability_emoped +
    b_D2C_emoped * start_d2c +
    b_O_emoped * outbound +
    b_W_emoped * weekend +
    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))
Post Reply