inst/examples/hybrid_choice_models/plot_h01_mode_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# h01. Baseline mode-choice logit model.
#
# This is the R counterpart of plot_h01_mode_logit.py. The choice model is
# written here so the example is self-contained: the shared files only load
# and prepare Optima data. variable() and biogeme_beta() create neutral
# symbolic nodes; native Biogeme receives the complete tree at estimate().

library(rbiogeme)

script_arguments <- commandArgs(trailingOnly = FALSE)
script_argument <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_argument) != 1L) {
  stop("Run this example as an R script with Rscript.", call. = FALSE)
}
example_directory <- dirname(normalizePath(sub("^--file=", "", script_argument)))
source(file.path(example_directory, "example_utils.R"))
source(file.path(example_directory, "optima.R"))

prepared <- prepare_optima_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "plot_h01_mode_logit",
  default_data = file.path(example_directory, "optima.dat")
)
database <- optima_database(prepared$data)

# The native specification fixes the common cost coefficient at -1 and
# estimates a positive scale parameter. A fixed R parameter is compiled as a
# native fixed Biogeme Beta, not as an R-side constant substitution.
choice_beta_cost <- biogeme_beta(
  "choice_beta_cost", start = -1, upper = 0, fixed = TRUE
)
choice_asc_car <- biogeme_beta("choice_asc_car", start = 0)
choice_asc_pt <- biogeme_beta("choice_asc_pt", start = 0)
choice_beta_dist_work <- biogeme_beta(
  "choice_beta_dist_work", start = 0, upper = 0
)
choice_beta_dist_other_purposes <- biogeme_beta(
  "choice_beta_dist_other_purposes", start = 0, upper = 0
)
choice_scale_parameter <- biogeme_beta(
  "choice_scale_parameter", start = 1, lower = 0.0001
)
choice_beta_time_car <- biogeme_beta(
  "choice_beta_time_car", start = 0, upper = 0
)
choice_beta_time_pt <- biogeme_beta(
  "choice_beta_time_pt", start = 0, upper = 0
)

work_trip <- variable("PurpHWH") == 1
other_trip_purposes <- variable("PurpHWH") != 1
choice_beta_dist <- choice_beta_dist_work * work_trip +
  choice_beta_dist_other_purposes * other_trip_purposes

# These are the native utilities for alternatives 0 (public transport),
# 1 (car), and 2 (slow modes). R arithmetic is overloaded, so this builds a
# symbolic expression rather than calculating a vector in R.
v_public_transport <- choice_asc_pt +
  choice_beta_time_pt * variable("TimePT_hour") +
  choice_beta_cost * variable("MarginalCostPT")
v_car <- choice_asc_car +
  choice_beta_time_car * variable("TimeCar_hour") +
  choice_beta_cost * variable("CostCarCHF")
v_slow_modes <- choice_beta_dist * variable("distance_km")

utilities <- list(
  `0` = choice_scale_parameter * v_public_transport,
  `1` = choice_scale_parameter * v_car,
  `2` = choice_scale_parameter * v_slow_modes
)
availability <- list(
  `0` = 1,
  `1` = variable("car_is_available"),
  `2` = 1
)

model <- logit_model(
  database = database,
  choice = "Choice",
  utilities = utilities,
  availability = availability
)
control <- optima_estimation_control(
  model_name = "plot_h01_mode_logit",
  prepared = prepared,
  # Native h01 leaves these numerical controls at their Biogeme defaults.
  numerically_safe = NULL,
  second_derivatives = NULL,
  max_iterations = NULL
)

# estimate() always starts a native estimation. This example therefore does
# not silently recycle a YAML or iteration file from an earlier run.
fit <- estimate(model, model_name = "plot_h01_mode_logit", control = control)
print(summary(fit))
print(coef(fit))
invisible(fit)