plot_b09wtp.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b09. Willingness to pay from utility derivatives.
#
# This is the R counterpart of plot_b09wtp.py. The model specification is
# complete in this file. Derive(), simulation, bootstrap draws, and confidence
# intervals are delegated to native Biogeme; R performs the final subgroup
# weighting and reporting.

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 Python and creates a fresh output directory. read_optima_database()
# is defined in optima.R and applies the native Optima preparation.
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 = "b09wtp"
)
database <- read_optima_database(
  prepared$data_path,
  name = "b09wtp_optima"
)

# Parameter and utility definitions match scenarios.py. Parameter names,
# starting values, bounds, and the fixed PT ASC are preserved exactly.
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 indicator expressions remain symbolic until the Python bridge
# compiles the complete model tree.
time_pt <- variable("TimePT")
time_car <- variable("TimeCar")
marginal_cost_pt <- variable("MarginalCostPT")
cost_car <- variable("CostCarCHF")
distance_km <- variable("distance_km")
time_pt_scaled <- time_pt / 200
time_car_scaled <- time_car / 200
cost_car_scaled <- cost_car / 10
distance_scaled <- 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 <- marginal_cost_pt / 10

# Native 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")
  )
)
log_probability <- nested_log_probability(
  utilities,
  availability = NULL,
  nests = nests,
  alternative = variable("Choice")
)

# WTP is the ratio of the marginal utility of time to the marginal utility of
# cost. Derive() is evaluated by native Biogeme, not by an R callback.
wtp_pt_time <- Derive(v_pt, "TimePT") / Derive(v_pt, "MarginalCostPT")
wtp_car_time <- Derive(v_car, "TimeCar") / Derive(v_car, "CostCarCHF")

# Estimate afresh with native bootstrap replications. The old b02estimation
# YAML file is never read or reused.
estimation_model <- biogeme_model(database = database, formula = 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 = "Willingness-to-pay indicators for the Optima example."
  ),
  run_bootstrap = prepared$run_bootstrap
)
if (!isTRUE(prepared$run_bootstrap)) {
  stop(
    "WTP confidence intervals require native bootstrap results; use --run-bootstrap=true.",
    call. = FALSE
  )
}

simulation_formulas <- list(
  weight = variable("normalized_weight"),
  `WTP PT time` = wtp_pt_time,
  `WTP CAR time` = wtp_car_time
)
simulation_model <- biogeme_model(database = database, simulations = simulation_formulas)
simulated_values <- as.data.frame(simulate(
  simulation_model,
  beta = fit,
  control = biogeme_control(
    output_directory = prepared$output,model_name = "b09wtp")
), check.names = FALSE)
print(utils::head(simulated_values))

# The native example reports time WTP per hour, hence the factor 60.
wtpcar <- mean(60 * simulated_values$`WTP CAR time` * simulated_values$weight)
bootstrap_betas <- indicator_bootstrap_parameter_draws(fit)
intervals <- biogeme_confidence_intervals(
  simulation_model,
  beta_values = bootstrap_betas,
  interval_size = 0.9,
  control = biogeme_control(
    output_directory = prepared$output,model_name = "b09wtp_ci")
)
left <- intervals$left
right <- intervals$right
print(utils::head(left))
print(utils::head(right))
wtpcar_left <- mean(60 * left$`WTP CAR time` * left$weight)
wtpcar_right <- mean(60 * right$`WTP CAR time` * right$weight)
cat(sprintf(
  "Average WTP for car: %.3g CI:[%.3g, %.3g]\n",
  wtpcar,
  wtpcar_left,
  wtpcar_right
))

unique_wtp <- unique(60 * simulated_values$`WTP CAR time`)
cat("Unique values:      ", paste(sprintf("%.3g", unique_wtp), collapse = ", "), "\n", sep = "")

# Normalize weights within each subgroup exactly as in the native example.
wtp_for_subgroup <- function(filter) {
  filter <- as.logical(filter)
  size <- sum(filter)
  if (size == 0L) stop("The requested WTP subgroup is empty.", call. = FALSE)
  subgroup_weight <- simulated_values$weight[filter]
  normalized_weight <- subgroup_weight * size / sum(subgroup_weight)
  value <- mean(60 * simulated_values$`WTP CAR time`[filter] * normalized_weight)
  lower <- mean(60 * left$`WTP CAR time`[filter] * normalized_weight)
  upper <- mean(60 * right$`WTP CAR time`[filter] * normalized_weight)
  c(value = value, lower = lower, upper = upper)
}

workers <- wtp_for_subgroup(database$data$OccupStat == 1)
females <- wtp_for_subgroup(database$data$Gender == 2)
males <- wtp_for_subgroup(database$data$Gender == 1)
cat(sprintf(
  "WTP car for workers: %.3g CI:[%.3g, %.3g]\n",
  workers[["value"]], workers[["lower"]], workers[["upper"]]
))
cat(sprintf(
  "WTP car for females: %.3g CI:[%.3g, %.3g]\n",
  females[["value"]], females[["lower"]], females[["upper"]]
))
cat(sprintf(
  "WTP car for males  : %.3g CI:[%.3g, %.3g]\n",
  males[["value"]], males[["lower"]], males[["upper"]]
))

# Plotting is optional for headless sessions. The histogram is weighted by
# normalized_weight, matching the native example's population display.
plot_option <- prepared$options$plot
if (!is.null(plot_option) && indicator_flag_option(plot_option, "plot")) {
  breaks <- pretty(range(60 * simulated_values$`WTP CAR time`), n = 20)
  groups <- cut(60 * simulated_values$`WTP CAR time`, breaks, include.lowest = TRUE)
  weighted_counts <- tapply(simulated_values$weight, groups, sum)
  mids <- (breaks[-1L] + breaks[-length(breaks)]) / 2
  png(file.path(prepared$output, "b09wtp.png"), width = 900, height = 600)
  plot(
    mids,
    as.numeric(weighted_counts),
    type = "h",
    lwd = 8,
    xlab = "WTP (CHF/hour)",
    ylab = "Individuals",
    main = "Willingness to pay for car travel time"
  )
  dev.off()
}

invisible(list(
  fit = fit,
  simulated_values = simulated_values,
  left = left,
  right = right,
  average_car_wtp = c(value = wtpcar, lower = wtpcar_left, upper = wtpcar_right),
  subgroups = rbind(workers = workers, females = females, males = males)
))