inst/examples/indicators/plot_b04market_shares.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b04. Calculation of market shares and confidence intervals.
#
# This is the R counterpart of plot_b04market_shares.py. The nested-logit
# specification is written out here; native Biogeme performs estimation,
# simulation, bootstrap handling, and confidence-interval quantiles.

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 documents the exact native
# data filter and derived variables used by all indicators examples.
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 = "b04market_shares"
)
database <- read_optima_database(
  prepared$data_path,
  name = "b04market_shares_optima"
)

# Parameter definitions match scenarios.py exactly. Parameter names and the
# fixed public-transport ASC are part of the native equivalence contract.
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 indicators are native symbolic expressions. Nothing here is
# evaluated by R before the complete expression tree reaches Biogeme.
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

# 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")
  )
)
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")
)

# Estimate afresh. The simulation-only formulas are kept in a separate model
# so normalized_weight is returned as a column, not used as an estimation
# weight. Bootstrap draws are generated by native Biogeme for the intervals.
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 = "Nested-logit estimation and market-share indicators for the Optima example."
  ),
  run_bootstrap = prepared$run_bootstrap
)

simulation_formulas <- list(
  weight = variable("normalized_weight"),
  `Prob. PT` = prob_pt,
  `Prob. car` = prob_car,
  `Prob. SM` = prob_sm
)
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 = "b04market_shares")
), check.names = FALSE)

# indicator_bootstrap_parameter_draws() only names the serialized native
# bootstrap rows. biogeme_confidence_intervals() delegates all repeated native
# simulation and quantile calculations to Biogeme's public API.
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 = "b04market_shares_ci")
)
left <- intervals$left
right <- intervals$right

# The reported market shares are weighted means of individual probabilities,
# exactly as in the native Pandas example.
simulated_values$`Weighted prob. car` <-
  simulated_values$weight * simulated_values$`Prob. car`
simulated_values$`Weighted prob. PT` <-
  simulated_values$weight * simulated_values$`Prob. PT`
simulated_values$`Weighted prob. SM` <-
  simulated_values$weight * simulated_values$`Prob. SM`
left$`Weighted prob. car` <- left$weight * left$`Prob. car`
left$`Weighted prob. PT` <- left$weight * left$`Prob. PT`
left$`Weighted prob. SM` <- left$weight * left$`Prob. SM`
right$`Weighted prob. car` <- right$weight * right$`Prob. car`
right$`Weighted prob. PT` <- right$weight * right$`Prob. PT`
right$`Weighted prob. SM` <- right$weight * right$`Prob. SM`

market_share_car <- mean(simulated_values$`Weighted prob. car`)
market_share_car_left <- mean(left$`Weighted prob. car`)
market_share_car_right <- mean(right$`Weighted prob. car`)
market_share_pt <- mean(simulated_values$`Weighted prob. PT`)
market_share_pt_left <- mean(left$`Weighted prob. PT`)
market_share_pt_right <- mean(right$`Weighted prob. PT`)
market_share_sm <- mean(simulated_values$`Weighted prob. SM`)
market_share_sm_left <- mean(left$`Weighted prob. SM`)
market_share_sm_right <- mean(right$`Weighted prob. SM`)

cat(
  sprintf(
    "Market share for car: %.1f%% [%.1f%%, %.1f%%]\n",
    100 * market_share_car,
    100 * market_share_car_left,
    100 * market_share_car_right
  )
)
cat(
  sprintf(
    "Market share for PT:  %.1f%% [%.1f%%, %.1f%%]\n",
    100 * market_share_pt,
    100 * market_share_pt_left,
    100 * market_share_pt_right
  )
)
cat(
  sprintf(
    "Market share for SM:   %.1f%% [%.1f%%, %.1f%%]\n",
    100 * market_share_sm,
    100 * market_share_sm_left,
    100 * market_share_sm_right
  )
)

invisible(list(
  fit = fit,
  simulated_values = simulated_values,
  left = left,
  right = right,
  market_shares = c(PT = market_share_pt, car = market_share_car, SM = market_share_sm)
))