inst/examples/indicators/plot_b08arc_elasticities.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b08. Arc elasticity for public transportation cost.
#
# This is the R counterpart of plot_b08arc_elasticities.py. Two complete
# native probability expressions are compared: the base scenario and a 20%
# public-transportation cost increase.

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 = "b08arc_elasticities",
  default_bootstrap_samples = 1L
)
database <- read_optima_database(
  prepared$data_path,
  name = "b08arc_elasticities_optima"
)

# Keep the native scenario construction visible. The factor affects only the
# marginal public-transportation cost; all parameter names and utility terms
# remain unchanged between scenarios.
build_indicator_scenario <- function(factor = 1.0, nests = NULL) {
  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)

  time_pt <- variable("TimePT")
  time_car <- variable("TimeCar")
  marginal_cost_pt <- variable("MarginalCostPT")
  cost_car <- variable("CostCarCHF")
  distance_km <- variable("distance_km")
  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 <- marginal_cost_pt * factor

  v_pt <- asc_pt + beta_time_fulltime * (time_pt / 200) * fulltime +
    beta_time_other * (time_pt / 200) * not_fulltime +
    beta_cost * (marginal_cost_scenario / 10)
  v_car <- asc_car + beta_time_fulltime * (time_car / 200) * fulltime +
    beta_time_other * (time_car / 200) * not_fulltime +
    beta_cost * (cost_car / 10)
  v_sm <- asc_sm + beta_dist_male * (distance_km / 5) * male +
    beta_dist_female * (distance_km / 5) * female +
    beta_dist_unreported * (distance_km / 5) * unreported_gender
  utilities <- list(`0` = v_pt, `1` = v_car, `2` = v_sm)
  if (is.null(nests)) {
    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")
      )
    )
  }
  list(
    utilities = utilities,
    nests = nests,
    prob_pt = nested_probability(utilities, NULL, nests, 0),
    log_probability = nested_log_probability(
      utilities,
      availability = NULL,
      nests = nests,
      alternative = variable("Choice")
    ),
    marginal_cost_scenario = marginal_cost_scenario
  )
}

base_scenario <- build_indicator_scenario(factor = 1.0)
after_scenario <- build_indicator_scenario(
  factor = 1.2,
  nests = base_scenario$nests
)

# Arc elasticity is the finite-change elasticity from the native example.
direct_elas_pt <- (
  after_scenario$prob_pt - base_scenario$prob_pt
) * base_scenario$marginal_cost_scenario /
  (
    base_scenario$prob_pt *
      (after_scenario$marginal_cost_scenario - base_scenario$marginal_cost_scenario)
  )

# Estimate the base model afresh rather than reading b02estimation.yaml.
estimation_model <- biogeme_model(
  database = database,
  formula = base_scenario$log_probability
)
fit <- estimate(
  estimation_model,
  model_name = "b02estimation",
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b02estimation",
    generate_html = TRUE,
    generate_yaml = TRUE,
    save_iterations = FALSE,
    user_notes = "Arc elasticities for the Optima indicators example."
  ),
  run_bootstrap = FALSE
)

simulation_formulas <- list(
  weight = variable("normalized_weight"),
  `Prob. PT` = base_scenario$prob_pt,
  direct_elas_pt = direct_elas_pt
)
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 = "b08arc_elasticities")
), check.names = FALSE)
print(utils::head(simulated_values))

simulated_values$`Weighted prob. PT` <-
  simulated_values$weight * simulated_values$`Prob. PT`
denominator_pt <- sum(simulated_values$`Weighted prob. PT`)
aggregate_direct_elas_pt <- sum(
  simulated_values$`Weighted prob. PT` *
    simulated_values$direct_elas_pt / denominator_pt,
  na.rm = TRUE
)
cat(sprintf(
  "Aggregate direct arc elasticity of public transportation wrt cost: %.3g\n",
  aggregate_direct_elas_pt
))

invisible(list(
  fit = fit,
  simulated_values = simulated_values,
  aggregate_direct_elas_pt = aggregate_direct_elas_pt
))