inst/examples/bayesian_swissmetro/plot_b23a_binary_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b23a. Bayesian binary logit model.
#
# The common Swissmetro preparation is followed by the binary-data filter,
# leaving Train (1) and Car (3). The binary alternatives retain their native
# Swissmetro codes rather than being renumbered to 1 and 2.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It reads --data, --python, and --output and prepares a clean run directory.
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_b23a_binary_logit_model <- function(database) {
  # swissmetro_binary.py removes Swissmetro choices and observations for which
  # Train or Car is unavailable. The operation remains a native database
  # expression and preserves filtered row identifiers.
  binary_exclude <- variable("CHOICE") == 2 |
    variable("CAR_AV_SP") == 0 |
    variable("TRAIN_AV_SP") == 0
  database <- biogeme_database_remove(database, binary_exclude)

  # Bounds and parameter names match the native Bayesian example exactly.
  asc_car <- biogeme_beta("asc_car", start = 0)
  b_time_car <- biogeme_beta("b_time_car", start = 0, upper = 0)
  b_time_train <- biogeme_beta("b_time_train", start = 0, upper = 0)
  b_cost_car <- biogeme_beta("b_cost_car", start = 0, upper = 0)
  b_cost_train <- biogeme_beta("b_cost_train", start = 0, upper = 0)

  v_train <- b_time_train * variable("TRAIN_TT_SCALED") +
    b_cost_train * variable("TRAIN_COST_SCALED")
  v_car <- asc_car + b_time_car * variable("CAR_TT_SCALED") +
    b_cost_car * variable("CAR_CO_SCALED")

  # logit_log_probability() compiles to native loglogit with the two native
  # alternative codes and their availability expressions.
  biogeme_model(
    database = database,
    formula = logit_log_probability(
      utilities = list(`1` = v_train, `3` = v_car),
      availability = list(
        `1` = variable("TRAIN_AV_SP"),
        `3` = variable("CAR_AV_SP")
      ),
      alternative = variable("CHOICE")
    ),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b23a_binary_logit",
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b23a_binary_logit"
)
unlink(file.path(prepared$output, c(
  "b23a_binary_logit.yaml",
  "b23a_binary_logit.nc",
  "b23a_binary_logit.html",
  "__b23a_binary_logit.iter"
)), force = TRUE)

database <- swissmetro_data(prepared$data)
model <- build_b23a_binary_logit_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b23a_binary_logit",
  control = model$control
)

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