inst/examples/swissmetro/plot_b18b_ordinal_probit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b18b. Ordinal probit model
#
# This example mirrors plot_b18b_ordinal_probit.py. It is intentionally an
# ordinal-probit illustration: the Swissmetro choice codes are not intrinsically
# ordered, but the native example treats 1 -> 2 -> 3 as ordered categories.

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, and
# --output 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_b18b_ordinal_probit_model <- function(database) {
  # Parameter names, starts, and bounds match the native Python example.
  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 = 0)

  # The second cutpoint is symbolic. OrderedProbit uses the native normal CDF;
  # the complete response graph is compiled before native estimation begins.
  tau2 <- tau1 + delta2
  utility <- b_time * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED")
  log_probability <- ordered_probit_log_probability(
    eta = utility,
    cutpoints = list(tau1, tau2),
    alternative = variable("CHOICE"),
    categories = c(1, 2, 3),
    neutral_labels = numeric()
  )

  # biogeme_model() is the generic formula interface. The bridge delegates
  # probability evaluation, derivatives, optimization, and reporting to
  # native Biogeme without running an R callback.
  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b18b_ordinal_probit",
      generate_html = TRUE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

# Always estimate afresh. Remove only exact b18b artifacts so an old YAML or
# iteration file cannot silently supply the estimates.
stale_files <- c(
  "b18b_ordinal_probit.yaml",
  "__b18b_ordinal_probit.iter",
  "b18b_ordinal_probit.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)
model <- build_b18b_ordinal_probit_model(database)
fit <- estimate(
  model,
  model_name = "b18b_ordinal_probit",
  control = model$control
)
print(summary(fit))
print(coef(fit))

invisible(fit)