inst/examples/swissmetro/plot_b01b_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01b. Linear utilities and automatic parameter segmentation
#
# This example uses the same three-alternative Swissmetro model as b01a, but
# illustrates LinearUtility and segmented alternative-specific constants.
# These R calls construct symbolic nodes that are compiled to the native
# Biogeme LinearUtility and Segmentation objects by the bridge.

library(rbiogeme)

# The shared helper contains only command-line parsing and data preparation;
# the complete b01b model specification is kept below in this script.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "example_utils.R"))

build_b01b_model <- function(database) {
  # A segmentation maps the values of a data variable to labels. The labels
  # become part of the generated parameter names, for example
  # asc_train_diff_male and asc_train_diff_with_ga.
  male <- biogeme_database_segmentation(
    database,
    "MALE",
    c(`0` = "female", `1` = "male")
  )
  ga <- biogeme_database_segmentation(
    database,
    "GA",
    c(`0` = "without_ga", `1` = "with_ga")
  )
  segmentations <- list(male, ga)

  # Starting values are part of the model specification and match the native
  # Python example.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = -1.28)
  b_cost <- biogeme_beta("b_cost", start = -1.08)

  # segment_beta() applies both segmentation schemes to an ASC. It leaves
  # the generic time and cost coefficients unsegmented.
  segmented_asc_car <- segment_beta(
    asc_car,
    segmentations
  )
  segmented_asc_train <- segment_beta(
    asc_train,
    segmentations
  )

  # A linear_term pairs one symbolic coefficient with one symbolic variable.
  # linear_utility() combines those pairs into a native LinearUtility node.
  term <- function(beta, name) linear_term(beta, variable(name))

  # As in b01a, the named list maps each utility and availability expression
  # to an alternative identifier. Numeric alternative names use backticks.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = segmented_asc_train + linear_utility(list(
        term(b_time, "TRAIN_TT_SCALED"),
        term(b_cost, "TRAIN_COST_SCALED")
      )),
      `2` = linear_utility(list(
        term(b_time, "SM_TT_SCALED"),
        term(b_cost, "SM_COST_SCALED")
      )),
      `3` = segmented_asc_car + linear_utility(list(
        term(b_time, "CAR_TT_SCALED"),
        term(b_cost, "CAR_CO_SCALED")
      ))
    ),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

# prepare_swissmetro_example() is defined in example_utils.R. It parses the
# command line, validates the data/Python paths, configures the bridge, reads
# the data, and creates a fresh output directory before the shared
# Swissmetro filtering and derived variables are applied below.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01b_logit"
)
arguments <- prepared$options
bootstrap_samples <- example_integer(
  if (is.null(arguments$bootstrap_samples)) "100" else arguments$bootstrap_samples,
  "bootstrap_samples"
)
run_bootstrap <- example_flag(arguments$run_bootstrap, default = TRUE)
database <- swissmetro_data(prepared$data)

# The model is compiled once when estimate() is called. No R callback is
# executed inside the native likelihood, gradient, Hessian, or bootstrap
# evaluations.
model <- build_b01b_model(database)
control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b01b_logit",
  second_derivatives = "never",
  bootstrap_samples = bootstrap_samples,
  variance_covariance_type = "BHHH",
  user_notes = paste(
    "Example of a logit model with LinearUtility and automatic",
    "segmentation of alternative-specific constants."
  ),
  generate_html = TRUE,
  generate_yaml = TRUE,
  save_iterations = FALSE
)
fit <- estimate(
  model,
  model_name = "b01b_logit",
  control = control,
  run_bootstrap = run_bootstrap
)

# The default is the 100 bootstrap replications used by the Python example;
# pass --bootstrap-samples=3 for a quick smoke test.
print(summary(fit))
print(coef(fit))
if (!is.null(fit$variance_covariance)) print(vcov(fit))
invisible(fit)