inst/examples/swissmetro/plot_b01a_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01a. Estimation of a multinomial logit model
#
# This example demonstrates the complete R-side specification of the
# Swissmetro model. Expressions such as variable("TRAIN_TT_SCALED") are
# symbolic Biogeme expressions: they do not evaluate the data in R. The
# complete expression tree is compiled once by the Python bridge, and native
# Python Biogeme performs the numerical estimation.

library(rbiogeme)

# Locate the shared command-line/data helper relative to this script. The
# model specification itself remains in this file so the example is readable
# without opening the package implementation.
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_b01a_model <- function(database) {
  # biogeme_beta() creates a symbolic model parameter. The Swissmetro ASC is
  # fixed at zero to identify the model; fixed = TRUE is the R equivalent of
  # a native Biogeme Beta with a fixed value.
  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_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)

  # Ordinary R arithmetic is overloaded for Biogeme expressions. Each
  # utility below is therefore a symbolic expression, not a numeric vector.
  # The backtick-quoted names are the alternative identifiers 1, 2, and 3.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm + b_time * variable("SM_TT_SCALED") +
        b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * 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")
    )
  )
}

# 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 can be supplied explicitly from any working directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01a_logit"
)
database <- swissmetro_data(prepared$data)

# Build the model and choose report/output controls. estimate() always calls
# native Biogeme estimation; it does not load an existing YAML result.
model <- build_b01a_model(database)
control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b01a_logit",
  generate_html = TRUE,
  generate_yaml = TRUE,
  save_iterations = FALSE
)
fit <- estimate(model, model_name = "b01a_logit", control = control)

# The returned object is an ordinary R result wrapper containing native
# Biogeme estimates and statistics.
print(summary(fit))
print(coef(fit))
invisible(fit)