inst/examples/swissmetro/plot_b23a_binary_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b23a. Binary logit model
#
# This example mirrors plot_b23a_binary_logit.py. The Swissmetro choice is
# removed, leaving Train (1) and Car (3). The complete model specification is
# written below so that the example can be read and run on its own.

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. It does not define the model.
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_b23a_binary_logit_model <- function(database) {
  # swissmetro_data() has already applied the common Swissmetro filter
  # (PURPOSE 1 or 3 and CHOICE not equal to 0) and defined the scaled
  # variables. The binary helper additionally keeps observations for which
  # Train and Car are both available and removes CHOICE == 2 (Swissmetro).
  binary_exclude <- variable("CHOICE") == 2 |
    variable("CAR_AV_SP") == 0 |
    variable("TRAIN_AV_SP") == 0
  database <- biogeme_database_remove(database, binary_exclude)

  # Parameter names, starting values, and fixed-status flags match the native
  # Python example exactly. None of these parameters is fixed.
  asc_car <- biogeme_beta("asc_car", start = 0)
  b_time_car <- biogeme_beta("b_time_car", start = 0)
  b_time_train <- biogeme_beta("b_time_train", start = 0)
  b_cost_car <- biogeme_beta("b_cost_car", start = 0)
  b_cost_train <- biogeme_beta("b_cost_train", start = 0)

  # There are only two utilities. Alternative codes remain the native
  # Swissmetro codes: Train = 1 and Car = 3.
  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")

  log_probability <- 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")
  )

  # biogeme_model() stores the complete symbolic likelihood. The bridge
  # compiles it once, and native Biogeme performs estimation and derivatives.
  # Output files are disabled here because this example always estimates
  # afresh and must not silently reuse a YAML or iteration file.
  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b23a_logit",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

database <- swissmetro_data(prepared$data)
model <- build_b23a_binary_logit_model(database)

# Estimate afresh. The native Python example may load
# saved_results/b23a_logit.yaml, but this R example never implicitly recycles
# an old YAML or iteration file.
fit <- estimate(
  model,
  model_name = "b23a_logit",
  control = model$control
)

print(summary(fit))
print(coef(fit))

invisible(fit)