inst/examples/assisted/plot_b10_parameter_overrides.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b10_parameter_overrides. Fix generated coefficients for a missing category.
#
# The didactic variable has_pt_subscr is coded 1 (subscription), 2 (no
# subscription), and -99 (missing). Native segmentation therefore generates a
# minus_99 coefficient. The override is applied to the complete native
# expression graph, so every catalog branch is handled consistently.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It supplies data paths and output setup only; the model specification below
# is complete and can be read without consulting another model file.
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 = "b10_parameter_overrides"
)

# Match native read_data(): remove only CHOICE == 0 before deriving the
# didactic subscription variable and changing the first retained observation.
database <- swissmetro_data(prepared$data, filter_purpose = FALSE)
database$data$has_pt_subscr <- ifelse(database$data$GA == 1, 1, 2)
database$data$has_pt_subscr[[1L]] <- -99

segmentation_pt_subscription <- biogeme_database_segmentation(
  database,
  "has_pt_subscr",
  c(`2` = "no_pt_subscr", `1` = "pt_subscr", `-99` = "minus_99"),
  reference = "no_pt_subscr"
)

# Native generated-segmentation names and parameter names are preserved.
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_pt_subscription),
  maximum_number = 1,
  selected_name = "has_pt_subscr"
)

v_train <- asc_catalogs[[1L]] +
  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_catalogs[[2L]] +
  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")
)
log_probability <- logit_log_probability(
  utilities,
  availability,
  variable("CHOICE")
)

control <- biogeme_control(
    output_directory = prepared$output,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)
base_model <- biogeme_model(
  database = database,
  formula = log_probability,
  availability = availability,
  control = control
)

# Discover the actual generated Beta names through native expression
# inspection. No names are guessed from the catalog labels.
all_parameter_names <- biogeme_native_parameter_names(
  base_model,
  controls = control
)$all
missing_parameter_names <- sort(
  all_parameter_names[endsWith(all_parameter_names, "_minus_99")]
)
if (length(missing_parameter_names) == 0L) {
  stop("The didactic minus_99 segmentation did not generate any parameters.", call. = FALSE)
}
overrides <- setNames(
  lapply(missing_parameter_names, function(name) 0),
  missing_parameter_names
)

# parameter_overrides is declarative: the bridge compiles these replacements
# and native Biogeme applies them before estimation. No R callback is run
# inside the likelihood or optimizer.
model <- biogeme_model(
  database = database,
  formula = log_probability,
  availability = availability,
  parameter_overrides = overrides,
  control = control
)
remaining_parameter_names <- biogeme_native_parameter_names(
  model,
  controls = control
)$all

cat("Generated minus_99 parameters: ", paste(missing_parameter_names, collapse = ", "), "\n", sep = "")
cat(
  "Parameters remaining after overrides: ",
  paste(intersect(missing_parameter_names, remaining_parameter_names), collapse = ", "),
  "\n",
  sep = ""
)

fit <- estimate(
  model,
  model_name = "b10_parameter_overrides",
  control = control
)
print(fit)

invisible(fit)