inst/examples/assisted/plot_simple_example.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# simple_example. Catalogs and assisted specification building blocks.
#
# This is a didactic counterpart of native plot_simple_example.py. It shows
# how several symbolic model alternatives can be kept in one catalog. The R
# side only creates the expression tree and asks native Biogeme for counts and
# exact configuration identifiers; it does not implement a catalog engine.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It supplies data and command-line/output configuration. All model equations
# and catalog definitions remain in this self-contained script.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "assisted_simple_example"
)

# Use filter_purpose = FALSE because native biogeme.data.swissmetro.read_data()
# removes CHOICE == 0 but retains every PURPOSE value for this tutorial.
database <- swissmetro_data(prepared$data, filter_purpose = FALSE)

# This helper mirrors the native tutorial's CentralController listing. The
# count and identifiers come from native Biogeme after the complete expression
# is compiled; the R package does not guess or recreate catalog combinations.
print_all_configurations <- function(expression, label) {
  model <- biogeme_model(database, formula = expression)
  total <- count_number_of_specifications(model)
  ids <- catalog_configuration_ids(model)
  cat("\n", label, "\n", sep = "")
  cat("Total: ", total, " configurations\n", sep = "")
  if (length(ids) > 0L) cat(paste(ids, collapse = "\n"), "\n", sep = "")
  invisible(model)
}

# Parameters used by the baseline logit and nested-logit specifications.
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)

# Utilities are ordinary R syntax overloaded for symbolic Biogeme nodes.
v_train <- asc_train + b_time * variable("TRAIN_TT_SCALED") +
  b_cost * variable("TRAIN_COST_SCALED")
v_swissmetro <- b_time * variable("SM_TT_SCALED") +
  b_cost * variable("SM_COST_SCALED")
v_car <- asc_car + b_time * variable("CAR_TT_SCALED") +
  b_cost * variable("CAR_CO_SCALED")
utilities <- list(`1` = v_train, `2` = v_swissmetro, `3` = v_car)
availability <- list(
  `1` = variable("TRAIN_AV_SP"),
  `2` = variable("SM_AV"),
  `3` = variable("CAR_AV_SP")
)
choice <- variable("CHOICE")

# Each branch below is a complete native log-likelihood expression.
log_probability_logit <- logit_log_probability(
  utilities, availability, choice
)
mu_existing <- biogeme_beta("mu_existing", start = 1, lower = 1, upper = 10)
existing <- nested_nest(mu_existing, alternatives = c(1, 3))
nests <- nested_nests(choice_set = c(1, 2, 3), nests = list(existing))
log_probability_nested <- nested_log_probability(
  utilities, availability, nests, choice
)

# A catalog is a symbolic expression with one named native specification per
# branch. Its default branch is the first one, but native configuration IDs can
# select any branch without exposing Python objects to the R user.
model_catalog <- catalog(
  "model_catalog",
  list(logit = log_probability_logit, nested = log_probability_nested)
)
print(model_catalog)
model_catalog_model <- biogeme_model(database, formula = model_catalog)
selected_configuration <- "model_catalog:nested"
cat(
  "Selected configuration is resolved natively: ",
  selected_configuration,
  "\n",
  sep = ""
)
print(biogeme_native_parameter_names(
  model_catalog_model,
  configuration_id = selected_configuration
))
for (specification_name in names(model_catalog$expressions)) {
  cat("Specification ", specification_name, ": ", sep = "")
  print(model_catalog$expressions[[specification_name]])
}
print_all_configurations(model_catalog, "Logit and nested-logit catalog")

# Nonlinear specifications: the same parameter can multiply a linear,
# Box--Cox, or squared travel-time expression. boxcox() is compiled to native
# Biogeme; R does not evaluate the transformation on the data.
lambda_travel_time <- biogeme_beta(
  "lambda_travel_time", start = 1, lower = -10, upper = 10
)
train_tt_catalog <- catalog(
  "train_tt_catalog",
  list(
    linear = variable("TRAIN_TT"),
    boxcox = boxcox(variable("TRAIN_TT"), lambda_travel_time),
    squared = variable("TRAIN_TT") * variable("TRAIN_TT")
  )
)
asc_train_nonlinear <- biogeme_beta("ASC_TRAIN", start = 0)
b_time_nonlinear <- biogeme_beta("B_TIME", start = 0, upper = 0)
v_train_catalog <- asc_train_nonlinear + b_time_nonlinear * train_tt_catalog
print_all_configurations(v_train_catalog, "Nonlinear train-time catalog")

# Two independent catalogs create all combinations of their choices.
car_tt_catalog <- catalog(
  "car_tt_catalog",
  list(
    linear = variable("CAR_TT"),
    boxcox = boxcox(variable("CAR_TT"), lambda_travel_time),
    squared = variable("CAR_TT") * variable("CAR_TT")
  )
)
print_all_configurations(
  train_tt_catalog + car_tt_catalog,
  "Unsynchronized train/car time catalogs"
)

# Supplying the same controller synchronizes the two catalogs: native Biogeme
# selects the same named transformation in both places.
car_tt_catalog_synchronized <- catalog(
  "car_tt_catalog",
  list(
    linear = variable("CAR_TT"),
    boxcox = boxcox(variable("CAR_TT"), lambda_travel_time),
    squared = variable("CAR_TT") * variable("CAR_TT")
  ),
  controller = train_tt_catalog$controller
)
print_all_configurations(
  train_tt_catalog + car_tt_catalog_synchronized,
  "Synchronized train/car time catalogs"
)

# A generic coefficient catalog has a generic and an alternative-specific
# branch. The helper preserves native parameter names and controller order.
generic_catalogs <- generic_alt_specific_catalogs(
  generic_name = "coefficients",
  beta_parameters = list(b_time_nonlinear, b_cost),
  alternatives = c("train", "car")
)
b_time_catalogs <- generic_catalogs[[1L]]
b_cost_catalogs <- generic_catalogs[[2L]]
print_all_configurations(
  b_time_catalogs$train * variable("TRAIN_TT") +
    b_cost_catalogs$train * variable("TRAIN_COST") +
    b_time_catalogs$car * variable("CAR_TT") +
    b_cost_catalogs$car * variable("CAR_CO"),
  "Synchronized generic/alternative-specific catalogs"
)

# Separate calls create independent controllers, so time and cost choices can
# vary independently even though each parameter's train/car branches remain
# synchronized with one another.
time_catalogs <- generic_alt_specific_catalogs(
  generic_name = "time_coefficient",
  beta_parameters = list(b_time_nonlinear),
  alternatives = c("train", "car")
)[[1L]]
cost_catalogs <- generic_alt_specific_catalogs(
  generic_name = "cost_coefficient",
  beta_parameters = list(b_cost),
  alternatives = c("train", "car")
)[[1L]]
print_all_configurations(
  time_catalogs$train * variable("TRAIN_TT") +
    cost_catalogs$train * variable("TRAIN_COST") +
    time_catalogs$car * variable("CAR_TT") +
    cost_catalogs$car * variable("CAR_CO"),
  "Independent time and cost catalogs"
)

# Native generate_segmentation() uses data values, but its segmentation
# definition is symbolic. Define COMMUTERS as a native derived variable so no
# R-side precomputed segment column is needed.
database <- biogeme_database_define_variable(
  database,
  "COMMUTERS",
  variable("PURPOSE") == 1
)
segmentation_purpose <- biogeme_database_segmentation(
  database,
  "COMMUTERS",
  c(`0` = "non_commuters", `1` = "commuters"),
  reference = "non_commuters"
)
segmentation_luggage <- biogeme_database_segmentation(
  database,
  "LUGGAGE",
  c(`0` = "no_lugg", `1` = "one_lugg", `3` = "several_lugg"),
  reference = "no_lugg"
)

# Segmentations are catalog choices. maximum_number controls whether zero,
# one, or two potential segmentations may be selected by native Biogeme.
asc_segmented <- segmentation_catalogs(
  generic_name = "asc",
  beta_parameters = list(asc_train_nonlinear, asc_car),
  potential_segmentations = list(segmentation_purpose, segmentation_luggage),
  maximum_number = 2
)
print_all_configurations(
  asc_segmented[[1L]] + asc_segmented[[2L]],
  "ASC catalogs with up to two segmentations"
)
asc_segmented_one <- segmentation_catalogs(
  generic_name = "asc",
  beta_parameters = list(asc_train_nonlinear, asc_car),
  potential_segmentations = list(segmentation_purpose, segmentation_luggage),
  maximum_number = 1
)
print_all_configurations(
  asc_segmented_one[[1L]] + asc_segmented_one[[2L]],
  "ASC catalogs with at most one segmentation"
)

# The same segmentation choices can be nested below an alternative-specific
# coefficient catalog. This mirrors the final two native tutorial blocks.
b_time_segmented_one <- generic_alt_specific_catalogs(
  generic_name = "b_time",
  beta_parameters = list(b_time_nonlinear),
  alternatives = c("train", "car"),
  potential_segmentations = list(segmentation_purpose, segmentation_luggage),
  maximum_number = 1
)[[1L]]
print_all_configurations(
  b_time_segmented_one$train,
  "Alternative-specific time coefficient with one segmentation"
)
b_time_segmented_two <- generic_alt_specific_catalogs(
  generic_name = "b_time",
  beta_parameters = list(b_time_nonlinear),
  alternatives = c("train", "car"),
  potential_segmentations = list(segmentation_purpose, segmentation_luggage),
  maximum_number = 2
)[[1L]]
print_all_configurations(
  b_time_segmented_two$train,
  "Alternative-specific time coefficient with two segmentations"
)

invisible(list(model = model_catalog_model, database = database))