inst/examples/bayesian_swissmetro/plot_b12_panel.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b12. Bayesian mixture of logit models with panel data.
#
# The random time coefficient and alternative-specific constants vary by
# individual. Bayesian estimation samples these distributed parameters
# directly, so this example does not add a Monte-Carlo integration node.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It reads --data, --python, and --output, configures native Biogeme, and
# prepares the Swissmetro data while keeping this model specification visible.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

build_b12_panel_model <- function(database) {
  positive_lower_bound <- 1e-5

  # biogeme_beta() is the R representation of native Beta. The fixed flag and
  # bounds preserve the Python constructor arguments exactly.
  b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
  b_time <- biogeme_beta("b_time", start = 0, upper = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1, lower = positive_lower_bound)
  b_time_eps <- draw("b_time_eps", "NORMAL")
  b_time_rnd <- distributed_parameter(
    "b_time_rnd",
    b_time + b_time_s * b_time_eps
  )

  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_car_s <- biogeme_beta("asc_car_s", start = 1, lower = positive_lower_bound)
  asc_car_eps <- draw("asc_car_eps", "NORMAL")
  asc_car_rnd <- distributed_parameter(
    "asc_car_rnd",
    asc_car + asc_car_s * asc_car_eps
  )

  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_train_s <- biogeme_beta("asc_train_s", start = 1, lower = positive_lower_bound)
  asc_train_eps <- draw("asc_train_eps", "NORMAL")
  asc_train_rnd <- distributed_parameter(
    "asc_train_rnd",
    asc_train + asc_train_s * asc_train_eps
  )

  asc_sm <- biogeme_beta("asc_sm", start = 0)
  asc_sm_s <- biogeme_beta("asc_sm_s", start = 1, lower = positive_lower_bound)
  asc_sm_eps <- draw("asc_sm_eps", "NORMAL")
  asc_sm_rnd <- distributed_parameter(
    "asc_sm_rnd",
    asc_sm + asc_sm_s * asc_sm_eps
  )

  utilities <- list(
    `1` = asc_train_rnd + b_time_rnd * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm_rnd + b_time_rnd * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car_rnd + 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")
  )

  # The conditional logit log-likelihood is compiled once and evaluated by
  # native Biogeme after the panel database has prepared individual draws.
  biogeme_model(
    database = database,
    formula = logit_log_probability(
      utilities = utilities,
      availability = availability,
      alternative = variable("CHOICE")
    ),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b12_panel",
      warmup = 10,
      bayesian_draws = 10,
      chains = 4,
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

# Remove only this example's artifacts so a previous YAML, NetCDF, HTML, or
# iteration file cannot silently change the result of a clean run.
unlink(file.path(prepared$output, c(
  "b12_panel.yaml",
  "b12_panel.nc",
  "b12_panel.html",
  "__b12_panel.iter"
)), force = TRUE)

# panel = TRUE declares ID as the panel identifier after the native-equivalent
# filter and derived-variable operations. The bridge validates contiguous IDs.
database <- swissmetro_data(prepared$data, panel = TRUE)
model <- build_b12_panel_model(database)
fit <- bayesian_estimate(model, model_name = "b12_panel", control = model$control)

print(summary(fit))
print(coef(fit))
print(bayesian_stored_variables(fit))
invisible(fit)