inst/examples/indicators/plot_b02estimation.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b02. Estimation and direct simulation of a nested-logit model.
#
# This is the R counterpart of plot_b02estimation.py. The complete Optima
# model specification is kept in this file. The helper files next to it only
# handle data preparation, command-line arguments, and output directories.

library(rbiogeme)

script_arguments <- commandArgs(trailingOnly = FALSE)
script_argument <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_argument) != 1L) {
  stop("Run this example as an R script with Rscript.", call. = FALSE)
}
example_directory <- dirname(normalizePath(sub("^--file=", "", script_argument)))

# prepare_indicator_estimation_example() is defined in indicator_utils.R. It
# configures the requested Python interpreter, selects optima.dat by default,
# and creates a fresh output directory so old YAML or iteration files cannot
# be reused. read_optima_database() is defined in optima.R and applies the
# exact native Choice/CarAvail filters and native derived-column definitions.
source(file.path(example_directory, "indicator_utils.R"))
source(file.path(example_directory, "optima.R"))

prepared <- prepare_indicator_estimation_example(
  commandArgs(trailingOnly = TRUE),
  example_directory = example_directory,
  default_model = "b02estimation"
)

database <- read_optima_database(
  prepared$data_path,
  name = "b02estimation_optima"
)

# Parameter definitions match scenarios.py exactly. Names, starting values,
# bounds, and the fixed public-transport ASC are part of the equivalence
# contract and are therefore written explicitly here.
asc_car <- biogeme_beta("asc_car", start = 0)
asc_pt <- biogeme_beta("asc_pt", start = 0, fixed = TRUE)
asc_sm <- biogeme_beta("asc_sm", start = 0)
beta_time_fulltime <- biogeme_beta("beta_time_fulltime", start = 0)
beta_time_other <- biogeme_beta("beta_time_other", start = 0)
beta_dist_male <- biogeme_beta("beta_dist_male", start = 0)
beta_dist_female <- biogeme_beta("beta_dist_female", start = 0)
beta_dist_unreported <- biogeme_beta("beta_dist_unreported", start = 0)
beta_cost <- biogeme_beta("beta_cost", start = 0)

# Scaling keeps the native parameter magnitudes near one. Comparisons create
# native 0/1 indicator expressions; they are not evaluated in R.
time_pt_scaled <- variable("TimePT") / 200
time_car_scaled <- variable("TimeCar") / 200
cost_car_scaled <- variable("CostCarCHF") / 10
distance_scaled <- variable("distance_km") / 5
male <- variable("Gender") == 1
female <- variable("Gender") == 2
unreported_gender <- variable("Gender") == -1
fulltime <- variable("OccupStat") == 1
not_fulltime <- variable("OccupStat") != 1
marginal_cost_pt_scaled <- variable("MarginalCostPT") / 10

# Utility expressions and the alternative mapping are the same as in
# scenarios.py: 0 = public transportation, 1 = car, and 2 = slow modes.
v_pt <- asc_pt + beta_time_fulltime * time_pt_scaled * fulltime +
  beta_time_other * time_pt_scaled * not_fulltime +
  beta_cost * marginal_cost_pt_scaled
v_car <- asc_car + beta_time_fulltime * time_car_scaled * fulltime +
  beta_time_other * time_car_scaled * not_fulltime +
  beta_cost * cost_car_scaled
v_sm <- asc_sm + beta_dist_male * distance_scaled * male +
  beta_dist_female * distance_scaled * female +
  beta_dist_unreported * distance_scaled * unreported_gender
utilities <- list(`0` = v_pt, `1` = v_car, `2` = v_sm)

# The no-car nest contains PT and slow modes. The car nest is trivial and is
# represented explicitly to preserve the native nest name and choice set.
mu_no_car <- biogeme_beta("mu_no_car", start = 1, lower = 1, upper = 2)
nests <- nested_nests(
  choice_set = c(0, 1, 2),
  nests = list(
    nested_nest(mu_no_car, alternatives = c(0, 2), name = "no_car"),
    nested_nest(1, alternatives = 1, name = "car")
  )
)

# lognested() is exposed as nested_log_probability(). Native Biogeme estimates
# this log probability; no R callback runs inside the likelihood, gradient,
# Hessian, optimizer, or bootstrap evaluations.
log_probability <- nested_log_probability(
  utilities = utilities,
  availability = NULL,
  nests = nests,
  alternative = variable("Choice")
)
model <- biogeme_model(database = database, formula = log_probability)

control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b02estimation",
  bootstrap_samples = prepared$bootstrap_samples,
  generate_html = TRUE,
  generate_yaml = TRUE,
  save_iterations = FALSE,
  user_notes = "Nested-logit estimation and direct simulation for the Optima indicators example."
)
fit <- estimate(
  model,
  model_name = "b02estimation",
  control = control,
  run_bootstrap = prepared$run_bootstrap
)

# This is the R presentation equivalent of Biogeme's estimated-parameter
# table. The result itself remains native in its numerical calculations.
print(summary(fit))

# Native get_value_c(..., aggregation = FALSE) returns one value per row. The
# public R wrapper compiles the expression first and returns a regular numeric
# vector only after native evaluation has finished.
simulated_log_probability <- evaluate_biogeme_expression_c(
  model = model,
  expression = log_probability,
  beta = fit,
  aggregation = FALSE,
  number_of_draws = 1000L,
  numerically_safe = FALSE,
  use_jit = TRUE
)
print(utils::head(data.frame(
  log_probability = simulated_log_probability,
  check.names = FALSE
)))

# Native get_value_c(..., aggregation = TRUE) is exposed by
# evaluate_biogeme_expression_c(), which returns the single aggregated scalar.
aggregated_log_likelihood <- evaluate_biogeme_expression_c(
  model = model,
  expression = log_probability,
  beta = fit,
  aggregation = TRUE,
  number_of_draws = 1000L,
  numerically_safe = FALSE,
  use_jit = TRUE
)
cat("Final log likelihood:     ", fit$final_log_likelihood, "\n", sep = "")
cat("Simulated log likelihood: ", aggregated_log_likelihood, "\n", sep = "")

invisible(list(
  fit = fit,
  simulated_log_probability = simulated_log_probability,
  aggregated_log_likelihood = aggregated_log_likelihood
))