inst/examples/swissmetro/plot_b19_individual_level_parameters.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b19. Calculation of individual-level parameters
#
# This example mirrors plot_b19_individual_level_parameters.py. It first
# estimates the b05a normal-mixture model, then simulates native Monte Carlo
# Numerator and Denominator expressions and forms their ratio for each row.

library(rbiogeme)

# prepare_swissmetro_example() is defined in example_utils.R. It parses the
# command line, validates the data/Python paths, configures the bridge, reads
# the data, and creates a fresh output directory. The --data, --python,
# --output, --draws, and --seed options work from any current working
# directory.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "example_utils.R"))

build_b19_components <- function(database, number_of_draws = 10000L, seed = 1223L) {
  # These parameter names and starting values match the b05a native model.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  b_cost <- biogeme_beta("b_cost", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1)

  # draw() creates the named native normal draw used by both estimation and
  # simulation. The random coefficient remains a symbolic Biogeme expression.
  b_time_rnd <- b_time + b_time_s * draw("b_time_rnd", "NORMAL")
  utilities <- list(
    `1` = asc_train + b_time_rnd * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time_rnd * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time_rnd * variable("CAR_TT_SCALED") +
      b_cost * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )
  prob_chosen <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )

  # MonteCarlo() is compiled to native Biogeme. R only combines the returned
  # columns after simulation; it never evaluates the likelihood or draws.
  numerator <- monte_carlo(b_time_rnd * prob_chosen)
  denominator <- monte_carlo(prob_chosen)
  simulations <- list(
    Numerator = numerator,
    Denominator = denominator,
    Choice = variable("CHOICE")
  )
  draws <- biogeme_draws(
    name = "b_time_rnd",
    draw_type = "NORMAL",
    number_of_draws = number_of_draws,
    seed = seed
  )
  list(
    model = biogeme_model(
      database = database,
      formula = log(denominator),
      draws = draws
    ),
    simulations = simulations,
    draws = draws
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b19_individual_level_parameters"
)

number_of_draws <- if (!is.null(prepared$options$draws) && nzchar(prepared$options$draws)) {
  example_integer(prepared$options$draws, "draws")
} else {
  10000L
}
seed <- if (!is.null(prepared$options$seed) && nzchar(prepared$options$seed)) {
  example_integer(prepared$options$seed, "seed")
} else {
  1223L
}

# Estimate the prerequisite afresh. The native b19 example reads a saved b05a
# YAML file; this R counterpart removes exact artifacts and does not recycle it.
stale_files <- c(
  "b05a_normal_mixture.yaml",
  "__b05a_normal_mixture.iter",
  "b05a_normal_mixture.html",
  "b19_individual_level_parameters.yaml",
  "__b19_individual_level_parameters.iter",
  "b19_individual_level_parameters.html"
)
stale_files <- file.path(prepared$output, stale_files)
stale_files <- stale_files[file.exists(stale_files)]
if (length(stale_files) > 0L) unlink(stale_files, force = TRUE)

database <- swissmetro_data(prepared$data)
components <- build_b19_components(database, number_of_draws, seed)
estimation_control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b05a_normal_mixture",
  number_of_draws = number_of_draws,
  seed = seed,
  analytical_hessian_mode = "automatic",
  generate_html = TRUE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)
cat(sprintf("Number of draws: %s\n", format(number_of_draws, big.mark = "_")))

fit <- estimate(
  components$model,
  model_name = "b05a_normal_mixture",
  control = estimation_control
)
print(summary(fit))
print(coef(fit))

# Simulate the named native expressions at the estimated beta values.
components$model$simulations <- components$simulations
simulation <- simulate(
  components$model,
  beta = fit,
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b19_individual_level_parameters",
    number_of_draws = number_of_draws,
    seed = seed,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)
simulation_values <- as.data.frame(simulation, check.names = FALSE)
simulation_values[["Individual-level parameters"]] <-
  simulation_values[["Numerator"]] / simulation_values[["Denominator"]]

print(utils::head(simulation_values))
invisible(simulation_values)