inst/examples/swissmetro/plot_b22a_multiple_models.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b22a. Assisted specification with a large catalog
#
# This example mirrors plot_b22a_multiple_models.py. The complete b22b catalog
# definition is included below so this script is self-contained. There are
# 504 possible combinations, so native Biogeme uses its assisted-specification
# heuristic rather than enumerating every model.

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. It does not define the model.
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_b22a_multiple_models_model <- function(database) {
  # These parameters and names match native b22b.
  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)
  b_headway <- biogeme_beta("b_headway", start = 0)

  # A segmentation maps data values to readable native parameter suffixes.
  gender_segmentation <- biogeme_database_segmentation(
    database,
    "MALE",
    c(`0` = "female", `1` = "male")
  )
  ga_segmentation <- biogeme_database_segmentation(
    database,
    "GA",
    c(`0` = "without_ga", `1` = "with_ga")
  )
  luggage_segmentation <- biogeme_database_segmentation(
    database,
    "LUGGAGE",
    c(`0` = "no_lugg", `1` = "one_lugg", `3` = "several_lugg")
  )

  # The ASC controller permits zero, one, or two of these segmentations for
  # both ASC catalogs. Native segmentation_catalogs() creates the same
  # synchronized choices and parameter names at compilation time.
  asc_catalogs <- segmentation_catalogs(
    "asc",
    list(asc_car, asc_train),
    list(gender_segmentation, luggage_segmentation, ga_segmentation),
    maximum_number = 2
  )
  asc_car_catalog <- asc_catalogs[[1L]]
  asc_train_catalog <- asc_catalogs[[2L]]

  # Headway is either omitted or included for Train and Swissmetro together.
  headway_controller <- catalog_controller(
    "train_headway_catalog",
    c("without_headway", "with_headway")
  )
  train_headway_catalog <- catalog(
    "train_headway_catalog",
    list(
      without_headway = 0,
      with_headway = b_headway * variable("TRAIN_HE")
    ),
    headway_controller
  )
  sm_headway_catalog <- catalog(
    "sm_headway_catalog",
    list(
      without_headway = 0,
      with_headway = b_headway * variable("SM_HE")
    ),
    headway_controller
  )

  # piecewise() is a symbolic native piecewise_formula node. NULL marks an
  # open endpoint; the finite breakpoints define the interval slopes.
  ell_tt <- biogeme_beta("lambda_tt", start = 1, lower = -10, upper = 10)
  train_tt_options <- list(
    linear = variable("TRAIN_TT_SCALED"),
    log = logzero(variable("TRAIN_TT_SCALED")),
    sqrt = variable("TRAIN_TT_SCALED") ^ 0.5,
    piecewise_1 = piecewise(variable("TRAIN_TT_SCALED"), list(0, 0.1, NULL)),
    piecewise_2 = piecewise(variable("TRAIN_TT_SCALED"), list(0, 0.25, NULL)),
    boxcox = boxcox(variable("TRAIN_TT_SCALED"), ell_tt)
  )
  time_controller <- catalog_controller(
    "train_tt_catalog",
    names(train_tt_options)
  )
  train_tt_catalog <- catalog("train_tt_catalog", train_tt_options, time_controller)
  sm_tt_catalog <- catalog(
    "sm_tt_catalog",
    list(
      linear = variable("SM_TT_SCALED"),
      log = logzero(variable("SM_TT_SCALED")),
      sqrt = variable("SM_TT_SCALED") ^ 0.5,
      piecewise_1 = piecewise(variable("SM_TT_SCALED"), list(0, 0.1, NULL)),
      piecewise_2 = piecewise(variable("SM_TT_SCALED"), list(0, 0.25, NULL)),
      boxcox = boxcox(variable("SM_TT_SCALED"), ell_tt)
    ),
    time_controller
  )
  car_tt_catalog <- catalog(
    "car_tt_catalog",
    list(
      linear = variable("CAR_TT_SCALED"),
      log = logzero(variable("CAR_TT_SCALED")),
      sqrt = variable("CAR_TT_SCALED") ^ 0.5,
      piecewise_1 = piecewise(variable("CAR_TT_SCALED"), list(0, 0.1, NULL)),
      piecewise_2 = piecewise(variable("CAR_TT_SCALED"), list(0, 0.25, NULL)),
      boxcox = boxcox(variable("CAR_TT_SCALED"), ell_tt)
    ),
    time_controller
  )

  # Travel cost has its own shared controller and its own Box-Cox parameter.
  ell_cost <- biogeme_beta("lambda_cost", start = 1, lower = -10, upper = 10)
  train_cost_options <- list(
    linear = variable("TRAIN_COST_SCALED"),
    log = logzero(variable("TRAIN_COST_SCALED")),
    sqrt = variable("TRAIN_COST_SCALED") ^ 0.5,
    piecewise_1 = piecewise(variable("TRAIN_COST_SCALED"), list(0, 0.1, NULL)),
    piecewise_2 = piecewise(variable("TRAIN_COST_SCALED"), list(0, 0.25, NULL)),
    boxcox = boxcox(variable("TRAIN_COST_SCALED"), ell_cost)
  )
  cost_controller <- catalog_controller(
    "train_cost_catalog",
    names(train_cost_options)
  )
  train_cost_catalog <- catalog(
    "train_cost_catalog",
    train_cost_options,
    cost_controller
  )
  sm_cost_catalog <- catalog(
    "sm_cost_catalog",
    list(
      linear = variable("SM_COST_SCALED"),
      log = logzero(variable("SM_COST_SCALED")),
      sqrt = variable("SM_COST_SCALED") ^ 0.5,
      piecewise_1 = piecewise(variable("SM_COST_SCALED"), list(0, 0.1, NULL)),
      piecewise_2 = piecewise(variable("SM_COST_SCALED"), list(0, 0.25, NULL)),
      boxcox = boxcox(variable("SM_COST_SCALED"), ell_cost)
    ),
    cost_controller
  )
  car_cost_catalog <- catalog(
    "car_cost_catalog",
    list(
      linear = variable("CAR_CO_SCALED"),
      log = logzero(variable("CAR_CO_SCALED")),
      sqrt = variable("CAR_CO_SCALED") ^ 0.5,
      piecewise_1 = piecewise(variable("CAR_CO_SCALED"), list(0, 0.1, NULL)),
      piecewise_2 = piecewise(variable("CAR_CO_SCALED"), list(0, 0.25, NULL)),
      boxcox = boxcox(variable("CAR_CO_SCALED"), ell_cost)
    ),
    cost_controller
  )

  utilities <- list(
    `1` = asc_train_catalog + b_time * train_tt_catalog +
      b_cost * train_cost_catalog + train_headway_catalog,
    `2` = b_time * sm_tt_catalog + b_cost * sm_cost_catalog + sm_headway_catalog,
    `3` = asc_car_catalog + b_time * car_tt_catalog + b_cost * car_cost_catalog
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  biogeme_model(
    database = database,
    formula = logit_log_probability(
      utilities = utilities,
      availability = availability,
      alternative = variable("CHOICE")
    ),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b22_multiple_models",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

# count_number_of_specifications() asks native controllers for the complete
# configuration count; it does not enumerate or estimate models in R.
number_of_specifications <- count_number_of_specifications(
  model,
  model_name = "b22_multiple_models",
  control = model$control
)
if (is.null(number_of_specifications)) {
  cat("There are too many possible specifications to be enumerated\n")
} else {
  cat(sprintf("There are %d possible specifications\n", number_of_specifications))
}

# force=TRUE removes the exact Pareto checkpoint and clears native assisted
# specification caches. The heuristic and all quick/final estimations remain
# native Biogeme operations.
pareto_file <- file.path(prepared$output, "b22_multiple_models.pareto")
fit <- assisted_specification(
  model,
  objectives = "aic_bic_dimension",
  pareto_file_name = pareto_file,
  model_name = "b22_multiple_models",
  control = model$control,
  force = TRUE
)

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]]))
  }
}

invisible(fit)