inst/examples/assisted/plot_b05alt_spec_segmentation.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b05alt_spec_segmentation. Combine segmentation and alternative-specific
# catalogs.
#
# Native Biogeme estimates 4 constant specifications, 6 time-coefficient
# specifications, and 2 cost-coefficient specifications: 48 combinations.
# R keeps the full symbolic specification visible; native Biogeme performs
# catalog resolution, estimation, and reporting.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It supplies data and run configuration only. All model syntax is specified
# below so this example can be read and reused on its own.
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"))

build_b05alt_spec_segmentation_model <- function(database) {
  # Match native read_data() and create COMMUTERS as a native derived column.
  database <- biogeme_database_define_variable(
    database,
    "COMMUTERS",
    variable("PURPOSE") == 1
  )
  segmentation_ga <- biogeme_database_segmentation(
    database,
    "GA",
    c(`0` = "noGA", `1` = "GA"),
    reference = "noGA"
  )
  segmentation_luggage <- biogeme_database_segmentation(
    database,
    "LUGGAGE",
    c(`0` = "no_lugg", `1` = "one_lugg", `3` = "several_lugg"),
    reference = "no_lugg"
  )
  segmentation_first <- biogeme_database_segmentation(
    database,
    "FIRST",
    c(`0` = "2nd_class", `1` = "1st_class"),
    reference = "2nd_class"
  )
  segmentation_purpose <- biogeme_database_segmentation(
    database,
    "COMMUTERS",
    c(`0` = "non_commuters", `1` = "commuters"),
    reference = "non_commuters"
  )

  # These four base Betas have the same names and starts as native b05.
  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)

  asc_catalogs <- segmentation_catalogs(
    generic_name = "asc",
    beta_parameters = list(asc_train, asc_car),
    potential_segmentations = list(segmentation_ga, segmentation_luggage),
    maximum_number = 2
  )

  # This helper nests segmentation choices below the generic/altspec choice.
  # One segmentation at most is allowed, matching native b05 exactly.
  b_time_catalogs <- generic_alt_specific_catalogs(
    generic_name = "b_time",
    beta_parameters = list(b_time),
    alternatives = c("train", "swissmetro", "car"),
    potential_segmentations = list(segmentation_first, segmentation_purpose),
    maximum_number = 1
  )[[1L]]
  b_cost_catalogs <- generic_alt_specific_catalogs(
    generic_name = "b_cost",
    beta_parameters = list(b_cost),
    alternatives = c("train", "swissmetro", "car")
  )[[1L]]

  utilities <- list(
    `1` = asc_catalogs[[1L]] +
      b_time_catalogs$train * variable("TRAIN_TT_SCALED") +
      b_cost_catalogs$train * variable("TRAIN_COST_SCALED"),
    `2` = b_time_catalogs$swissmetro * variable("SM_TT_SCALED") +
      b_cost_catalogs$swissmetro * variable("SM_COST_SCALED"),
    `3` = asc_catalogs[[2L]] +
      b_time_catalogs$car * variable("CAR_TT_SCALED") +
      b_cost_catalogs$car * variable("CAR_CO_SCALED")
  )
  log_probability <- logit_log_probability(
    utilities = utilities,
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    ),
    alternative = variable("CHOICE")
  )
  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b05alt_spec_segmentation",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

# Native biogeme.data.swissmetro.read_data() removes only CHOICE == 0.
database <- swissmetro_data(prepared$data, filter_purpose = FALSE)
model <- build_b05alt_spec_segmentation_model(database)

# Estimate all 48 combinations afresh; no prior YAML or iteration file is read.
fit <- estimate_catalog(
  model,
  model_name = "b05alt_spec_segmentation",
  control = model$control,
  force = TRUE
)

cat("A total of ", length(fit$results), " models have been estimated.\n", sep = "")
for (configuration in names(fit$results)) {
  result <- fit$results[[configuration]]
  cat(
    configuration,
    ": LL=",
    formatC(result$final_log_likelihood, digits = 2, format = "f"),
    " K=",
    length(result$beta_names),
    "\n",
    sep = ""
  )
}
print(fit$summary)
for (name in names(fit$description)) {
  if (!identical(name, unname(fit$description[[name]]))) {
    cat(name, "\t", fit$description[[name]], "\n", sep = "")
  }
}
cat("Non dominated models:\n")
for (configuration in fit$non_dominated) cat(configuration, "\n", sep = "")
print(fit$non_dominated_summary)
cat(fit$latex, "\n")

invisible(fit)