inst/examples/swissmetro/plot_b20_multiple_models.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b20. Estimation of several specifications from synchronized catalogs
#
# This example mirrors plot_b20_multiple_models.py. Catalogs are neutral R
# expression nodes until the complete likelihood is compiled. Native Biogeme
# then enumerates the synchronized catalog choices and performs every model
# estimation and result-processing operation.

library(rbiogeme)

# 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 the output directory. This shared helper is only data
# and run setup; the full b20 model specification is below.
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_b20_multiple_models_model <- function(database) {
  # These four names and starting values are exactly those in native b20.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)

  # A catalog_controller() synchronizes selections across catalogs. Here the
  # same ASC choice (unsegmented or MALE-segmented) is used for Train and Car.
  gender <- biogeme_database_segmentation(
    database,
    variable = "MALE",
    mapping = c(`0` = "female", `1` = "male")
  )
  asc_controller <- catalog_controller("asc", c("no_seg", "MALE"))
  asc_train_catalog <- catalog(
    "segmented_asc_train",
    list(
      no_seg = asc_train,
      MALE = segment_beta(asc_train, list(gender))
    ),
    controller = asc_controller
  )
  asc_car_catalog <- catalog(
    "segmented_asc_car",
    list(
      no_seg = asc_car,
      MALE = segment_beta(asc_car, list(gender))
    ),
    controller = asc_controller
  )

  # A second shared controller forces Train and Swissmetro to use either a
  # linear travel-time variable or its native log transformation together.
  time_controller <- catalog_controller(
    "train_tt_catalog",
    c("linear", "log")
  )
  train_tt_catalog <- catalog(
    "train_tt_catalog",
    list(
      linear = variable("TRAIN_TT_SCALED"),
      log = log(variable("TRAIN_TT_SCALED"))
    ),
    controller = time_controller
  )
  sm_tt_catalog <- catalog(
    "sm_tt_catalog",
    list(
      linear = variable("SM_TT_SCALED"),
      log = log(variable("SM_TT_SCALED"))
    ),
    controller = time_controller
  )

  # Catalog objects participate in ordinary symbolic arithmetic. No R
  # callback is evaluated when these utilities are estimated in Python.
  utilities <- list(
    `1` = asc_train_catalog + b_time * train_tt_catalog +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = b_time * sm_tt_catalog + b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car_catalog + b_time * variable("CAR_TT_SCALED") +
      b_cost * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  # logit_log_probability() maps directly to native models.loglogit, the
  # numerically stable log-likelihood constructor used by Python b20.
  log_probability <- logit_log_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )

  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b20multiple_models",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b20multiple_models"
)
database <- swissmetro_data(prepared$data)
model <- build_b20_multiple_models_model(database)

# force = TRUE makes every catalog configuration call native estimate(), so
# an old YAML or iteration file in this directory cannot be reused.
fit <- estimate_catalog(
  model,
  model_name = "b20multiple_models",
  control = model$control,
  force = TRUE
)

cat(sprintf("A total of %d models have been estimated:\n", length(fit$results)))
for (configuration in names(fit$results)) {
  result <- fit$results[[configuration]]
  cat(sprintf(
    "%s: LL=%.2f K=%d\n",
    configuration,
    result$final_log_likelihood,
    length(result$beta_names)
  ))
}

# These tables and the Pareto set are produced by native Biogeme result
# processing and returned as ordinary R tables/vectors for inspection.
print(fit$summary)
for (name in names(fit$description)) {
  if (!identical(name, unname(fit$description[[name]]))) {
    cat(sprintf("%s: %s\n", name, fit$description[[name]]))
  }
}
cat(sprintf("Out of them, %d are non dominated.\n", length(fit$non_dominated)))
for (configuration in fit$non_dominated) cat(configuration, "\n", sep = "")
print(fit$non_dominated_summary)
cat(fit$latex, "\n")

invisible(fit)