plot_b05revenues.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b05. Calculation of public-transportation revenues.
#
# This is the R counterpart of plot_b05revenues.py. The price scenarios and
# model equations are specified here, while native Biogeme evaluates every
# probability, bootstrap simulation, and confidence interval.

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
# selects optima.dat, configures Python, and creates a fresh output directory.
# read_optima_database() is defined in optima.R and reproduces native Optima
# filtering and derived variables. No revenue cache is read by this R port.
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 = "b05revenues"
)
database <- read_optima_database(
  prepared$data_path,
  name = "b05revenues_optima"
)

# Keep all parameter and utility definitions in one visible function. Calling
# it with a different factor changes only the public-transportation cost term;
# parameter names and bounds remain identical to scenarios.py.
build_indicator_specification <- function(factor = 1.0) {
  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)
  mu_no_car <- biogeme_beta("mu_no_car", start = 1, lower = 1, upper = 2)

  # Scaling and comparisons construct native Biogeme expressions.
  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_scenario <- variable("MarginalCostPT") * factor
  marginal_cost_pt_scaled <- marginal_cost_scenario / 10

  # Alternative coding is 0 = PT, 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)
  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")
    )
  )
  prob_pt <- nested_probability(utilities, NULL, nests, 0)
  prob_car <- nested_probability(utilities, NULL, nests, 1)
  prob_sm <- nested_probability(utilities, NULL, nests, 2)
  log_probability <- nested_log_probability(
    utilities,
    availability = NULL,
    nests = nests,
    alternative = variable("Choice")
  )
  list(
    log_probability = log_probability,
    prob_pt = prob_pt,
    prob_car = prob_car,
    prob_sm = prob_sm,
    marginal_cost_scenario = marginal_cost_scenario
  )
}

# Estimate the base scenario afresh. Bootstrap results are required below for
# the native 90% confidence intervals. Use --run-bootstrap=false only when
# running a reduced smoke test that does not request confidence intervals.
base_specification <- build_indicator_specification(factor = 1.0)
estimation_model <- biogeme_model(
  database = database,
  formula = base_specification$log_probability
)
fit <- estimate(
  estimation_model,
  model_name = "b02estimation",
  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 revenue indicators for the Optima example."
  ),
  run_bootstrap = prepared$run_bootstrap
)

if (!isTRUE(prepared$run_bootstrap)) {
  stop(
    "Revenue confidence intervals require native bootstrap results; use --run-bootstrap=true.",
    call. = FALSE
  )
}
bootstrap_betas <- indicator_bootstrap_parameter_draws(fit)

# The native example investigates factors from 0.0 through 4.9 by 0.1. A
# comma-separated --factors=... option is provided for quick smoke tests while
# preserving that exact default grid.
factor_option <- prepared$options$factors
if (is.null(factor_option)) {
  factors <- seq(0.0, 4.9, by = 0.1)
} else {
  factors <- suppressWarnings(as.numeric(strsplit(factor_option, ",", fixed = TRUE)[[1L]]))
  if (length(factors) == 0L || anyNA(factors) || any(!is.finite(factors)) ||
      any(factors < 0)) {
    stop("factors must be a comma-separated list of non-negative numbers.", call. = FALSE)
  }
}

calculate_revenue <- function(factor) {
  specification <- build_indicator_specification(factor)
  simulation_formulas <- list(
    weight = variable("normalized_weight"),
    `Revenue public transportation` =
      specification$prob_pt * specification$marginal_cost_scenario
  )
  simulation_model <- biogeme_model(
    database = database,
    simulations = simulation_formulas
  )
  simulated <- as.data.frame(simulate(
    simulation_model,
    beta = fit,
    control = biogeme_control(
    output_directory = prepared$output,model_name = "b05revenues")
  ), check.names = FALSE)
  intervals <- biogeme_confidence_intervals(
    simulation_model,
    beta_values = bootstrap_betas,
    interval_size = 0.9,
    control = biogeme_control(
    output_directory = prepared$output,model_name = "b05revenues_ci")
  )
  left <- intervals$left
  right <- intervals$right
  revenue <- sum(simulated$weight * simulated$`Revenue public transportation`)
  revenue_left <- sum(left$weight * left$`Revenue public transportation`)
  revenue_right <- sum(right$weight * right$`Revenue public transportation`)
  c(revenue = revenue, lower = revenue_left, upper = revenue_right)
}

revenue_values <- do.call(rbind, lapply(factors, calculate_revenue))
revenue_values <- as.data.frame(revenue_values, check.names = FALSE)
revenue_values$factor <- factors
revenue_values <- revenue_values[, c("factor", "revenue", "lower", "upper")]

current_index <- which.min(abs(revenue_values$factor - 1.0))
current <- revenue_values[current_index, ]
cat(
  sprintf(
    "Total revenues for public transportation (for the sample): %.1f CHF [%.1f CHF, %.1f CHF]\n",
    current$revenue,
    current$lower,
    current$upper
  )
)
largest_index <- which.max(revenue_values$revenue)
largest <- revenue_values[largest_index, ]
cat(
  sprintf(
    "Largest revenue: %.1f obtained with factor %.1f\n",
    largest$revenue,
    largest$factor
  )
)

# Plotting is optional so the example remains runnable in headless sessions.
plot_option <- prepared$options$plot
if (!is.null(plot_option) && indicator_flag_option(plot_option, "plot")) {
  png(file.path(prepared$output, "b05revenues.png"), width = 900, height = 600)
  plot(
    revenue_values$factor,
    revenue_values$revenue,
    type = "l",
    ylim = range(c(revenue_values$lower, revenue_values$upper)),
    xlab = "Public transportation cost factor",
    ylab = "Revenue (CHF)",
    main = "Public transportation revenues"
  )
  lines(revenue_values$factor, revenue_values$lower, lty = 2)
  lines(revenue_values$factor, revenue_values$upper, lty = 2)
  dev.off()
}

invisible(list(fit = fit, revenue_values = revenue_values))