inst/examples/bayesian_swissmetro/plot_b18a_ordinal_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b18a. Bayesian ordinal logit model.
#
# The Swissmetro alternatives are not intrinsically ordered; this example
# nevertheless treats the codes 1, 2, and 3 as ordered categories to mirror
# the native Biogeme illustration.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles --data, --python, and --output so this script is self-contained.
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_b18a_ordinal_logit_model <- function(database) {
  positive_lower_bound <- 1e-5
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  tau1 <- biogeme_beta("tau1", start = -1, upper = 0)
  delta2 <- biogeme_beta("delta2", start = 2, lower = positive_lower_bound)
  tau2 <- tau1 + delta2
  utility <- b_time * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED")

  # ordered_logit_log_probability() compiles to native OrderedLogLogit,
  # including the cutpoint ordering and neutral-label semantics.
  biogeme_model(
    database = database,
    formula = ordered_logit_log_probability(
      eta = utility,
      cutpoints = list(tau1, tau2),
      alternative = variable("CHOICE"),
      categories = c(1, 2, 3),
      neutral_labels = numeric()
    ),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b18a_ordinal_logit",
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

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

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