inst/examples/assisted/plot_b09post_processing.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b09post_processing. Re-estimate the Pareto-optimal models from b07.
#
# The Pareto checkpoint is an intentional input to this example. Pass it
# explicitly with --pareto=/path/to/b07everything_assisted.pareto; the script
# never searches for or silently reuses an old checkpoint.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles paths and output setup only. The complete catalog model is written
# below so this example remains self-contained for R users.
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_b09everything_model <- function(database) {
  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 parameter definitions preserve the native names and bounds.
  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)
  lambda_travel_time <- biogeme_beta(
    "lambda_travel_time",
    start = 1,
    lower = -10,
    upper = 10
  )
  square_tt_coef <- biogeme_beta("square_tt_coef", start = 0)
  cube_tt_coef <- biogeme_beta("cube_tt_coef", start = 0)
  power_series <- function(the_variable) {
    the_variable + square_tt_coef * the_variable^2 +
      cube_tt_coef * the_variable * the_variable^3
  }

  # One controller synchronizes the functional form in all three alternatives.
  time_controller <- catalog_controller(
    "train_tt_catalog",
    c("linear", "boxcox", "power")
  )
  train_tt_catalog <- catalog(
    "train_tt_catalog",
    list(
      linear = variable("TRAIN_TT_SCALED"),
      boxcox = boxcox(variable("TRAIN_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("TRAIN_TT_SCALED"))
    ),
    controller = time_controller
  )
  sm_tt_catalog <- catalog(
    "sm_tt_catalog",
    list(
      linear = variable("SM_TT_SCALED"),
      boxcox = boxcox(variable("SM_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("SM_TT_SCALED"))
    ),
    controller = time_controller
  )
  car_tt_catalog <- catalog(
    "car_tt_catalog",
    list(
      linear = variable("CAR_TT_SCALED"),
      boxcox = boxcox(variable("CAR_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("CAR_TT_SCALED"))
    ),
    controller = time_controller
  )
  asc_catalogs <- segmentation_catalogs(
    generic_name = "asc",
    beta_parameters = list(asc_train, asc_car),
    potential_segmentations = list(segmentation_ga, segmentation_luggage),
    maximum_number = 2
  )
  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 * train_tt_catalog +
      b_cost_catalogs$train * variable("TRAIN_COST_SCALED"),
    `2` = b_time_catalogs$swissmetro * sm_tt_catalog +
      b_cost_catalogs$swissmetro * variable("SM_COST_SCALED"),
    `3` = asc_catalogs[[2L]] +
      b_time_catalogs$car * car_tt_catalog +
      b_cost_catalogs$car * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )
  choice <- variable("CHOICE")
  logit <- logit_log_probability(utilities, availability, choice)

  mu_existing <- biogeme_beta("mu_existing", start = 1, lower = 1, upper = 10)
  nested_existing <- nested_log_probability(
    utilities,
    availability,
    nested_nests(
      c(1, 2, 3),
      list(nested_nest(mu_existing, c(1, 3), name = "Existing"))
    ),
    choice
  )
  mu_public <- biogeme_beta("mu_public", start = 1, lower = 1, upper = 10)
  nested_public <- nested_log_probability(
    utilities,
    availability,
    nested_nests(
      c(1, 2, 3),
      list(nested_nest(mu_public, c(1, 2), name = "Public"))
    ),
    choice
  )
  model_catalog <- catalog(
    "model_catalog",
    list(
      logit = logit,
      `nested existing` = nested_existing,
      `nested public` = nested_public
    )
  )
  biogeme_model(
    database = database,
    formula = model_catalog,
    availability = availability,
    control = biogeme_control(
    output_directory = prepared$output,
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b09post_processing"
)
pareto_argument <- prepared$options$pareto
if (is.null(pareto_argument) || !nzchar(pareto_argument)) {
  stop(
    "Provide --pareto=/path/to/b07everything_assisted.pareto; " ,
    "the Pareto checkpoint is an explicit input.",
    call. = FALSE
  )
}
pareto_file <- normalizePath(pareto_argument, mustWork = TRUE)

# Match native read_data(): remove only CHOICE == 0.
database <- swissmetro_data(prepared$data, filter_purpose = FALSE)
model <- build_b09everything_model(database)

# recycle = TRUE is intentional here: native b09 demonstrates re-estimation
# from a Pareto checkpoint and permits native result-file recycling.
plot_file <- if (example_flag(prepared$options$plot, default = TRUE)) {
  file.path(prepared$output, "b09_process_pareto.png")
} else {
  NULL
}
fit <- pareto_post_processing(
  model,
  pareto_file_name = pareto_file,
  model_name = "b09post_processing",
  control = model$control,
  recycle = TRUE,
  plot_file_name = plot_file,
  objective_x = 1L,
  objective_y = 0L,
  label_x = "Nbr of parameters",
  label_y = "Negative log likelihood"
)

cat("Pareto models re-estimated: ", length(fit$results), "\n", sep = "")
print(fit$summary)
if (length(fit$results) > 0L) {
  specification_id <- names(fit$results)[[1L]]
  result <- fit$results[[1L]]
  cat("First specification: ", specification_id, "\n", sep = "")
  print(result)
  standard_errors <- result$standard_errors
  t_statistics <- result$t_statistics
  p_values <- result$p_values
  if (is.null(standard_errors)) {
    standard_errors <- rep(NA_real_, length(result$beta_names))
    t_statistics <- standard_errors
    p_values <- standard_errors
  }
  estimated_parameters <- data.frame(
    Value = unname(result$beta_values),
    `Std err` = unname(standard_errors),
    `t-test` = unname(t_statistics),
    `p-value` = unname(p_values),
    row.names = result$beta_names,
    check.names = FALSE
  )
  print(estimated_parameters)
}
if (!is.null(plot_file)) cat("Pareto plot: ", plot_file, "\n", sep = "")
for (message in fit$pareto_statistics) cat(message, "\n", sep = "")

invisible(fit)