inst/examples/indicators/plot_b06point_elasticities.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b06. Direct point elasticities.
#
# This is the R counterpart of plot_b06point_elasticities.py. Derive() creates
# native symbolic derivatives; Biogeme performs their evaluation and the R
# code only carries out the final reported weighted aggregates.

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

# These parameter definitions and utility equations are identical to
# scenarios.py. The ASC for PT is fixed; all other model parameters are free.
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)

# Comparisons and scaling remain neutral Biogeme expressions until compilation.
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

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

# Derive() is the R spelling of native Biogeme's Derive operator. It is
# compiled symbolically, so no R derivative callback runs during simulation.
direct_elas_pt_time <- Derive(prob_pt, "TimePT") * time_pt / prob_pt
direct_elas_pt_cost <- Derive(prob_pt, "MarginalCostPT") * marginal_cost_pt / prob_pt
direct_elas_car_time <- Derive(prob_car, "TimeCar") * time_car / prob_car
direct_elas_car_cost <- Derive(prob_car, "CostCarCHF") * cost_car / prob_car
direct_elas_sm_dist <- Derive(prob_sm, "distance_km") * distance_km / prob_sm

# Estimate the point values afresh. This is the same likelihood used to create
# b02estimation.yaml in the native documentation, but never loads that file.
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",
    generate_html = TRUE,
    generate_yaml = TRUE,
    save_iterations = FALSE,
    user_notes = "Direct point elasticities for the Optima indicators example."
  ),
  run_bootstrap = FALSE
)

simulation_formulas <- list(
  weight = variable("normalized_weight"),
  `Prob. car` = prob_car,
  `Prob. public transportation` = prob_pt,
  `Prob. slow modes` = prob_sm,
  direct_elas_pt_time = direct_elas_pt_time,
  direct_elas_pt_cost = direct_elas_pt_cost,
  direct_elas_car_time = direct_elas_car_time,
  direct_elas_car_cost = direct_elas_car_cost,
  direct_elas_sm_dist = direct_elas_sm_dist
)
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 = "b06point_elasticities")
), check.names = FALSE)
print(utils::head(simulated_values))

# Compute the same weighted aggregate indicators as the native example.
simulated_values$`Weighted prob. car` <-
  simulated_values$weight * simulated_values$`Prob. car`
simulated_values$`Weighted prob. PT` <-
  simulated_values$weight * simulated_values$`Prob. public transportation`
simulated_values$`Weighted prob. SM` <-
  simulated_values$weight * simulated_values$`Prob. slow modes`
denominator_car <- sum(simulated_values$`Weighted prob. car`)
denominator_pt <- sum(simulated_values$`Weighted prob. PT`)
denominator_sm <- sum(simulated_values$`Weighted prob. SM`)

aggregate_car_time <- sum(
  simulated_values$`Weighted prob. car` *
    simulated_values$direct_elas_car_time / denominator_car
)
aggregate_car_cost <- sum(
  simulated_values$`Weighted prob. car` *
    simulated_values$direct_elas_car_cost / denominator_car
)
aggregate_pt_time <- sum(
  simulated_values$`Weighted prob. PT` *
    simulated_values$direct_elas_pt_time / denominator_pt
)
aggregate_pt_cost <- sum(
  simulated_values$`Weighted prob. PT` *
    simulated_values$direct_elas_pt_cost / denominator_pt
)
aggregate_sm_distance <- sum(
  simulated_values$`Weighted prob. SM` *
    simulated_values$direct_elas_sm_dist / denominator_sm
)

cat(sprintf("Aggregate direct point elasticity of car wrt time: %.3g\n", aggregate_car_time))
cat(sprintf("Aggregate direct point elasticity of car wrt cost: %.3g\n", aggregate_car_cost))
cat(sprintf("Aggregate direct point elasticity of PT wrt time: %.3g\n", aggregate_pt_time))
cat(sprintf("Aggregate direct point elasticity of PT wrt cost: %.3g\n", aggregate_pt_cost))
cat(sprintf("Aggregate direct point elasticity of SM wrt distance: %.3g\n", aggregate_sm_distance))

invisible(list(
  fit = fit,
  simulated_values = simulated_values,
  aggregate = c(
    car_time = aggregate_car_time,
    car_cost = aggregate_car_cost,
    pt_time = aggregate_pt_time,
    pt_cost = aggregate_pt_cost,
    sm_distance = aggregate_sm_distance
  )
))