plot_b02_weight.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b02. Estimation with weights: WESML
#
# This example reproduces the native Swissmetro weighted-exogenous-sample
# maximum-likelihood (WESML) model. The utilities are specified in R, while
# the likelihood and weight expressions are compiled once and evaluated by
# native Python Biogeme.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# complete model specification remains in this script.
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_b02_weight_model <- function(database) {
  # biogeme_beta() creates symbolic native parameters. Fixing the Swissmetro
  # ASC at zero identifies the utility scale, exactly as in the Python model.
  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)

  # R arithmetic is overloaded for Biogeme expressions. These are symbolic
  # utility trees, not R vectors evaluated before the native estimation.
  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")
  )

  # GROUP is a raw column retained by swissmetro_data(). The comparisons are
  # native indicator expressions. Group 1 receives zero weight, while groups
  # 2 and 3 receive the exact constants used by plot_b02_weight.py.
  weight_group_2 <- 8.890991e-01
  weight_group_3 <- 1.2
  weight <- weight_group_2 * (variable("GROUP") == 2) +
    weight_group_3 * (variable("GROUP") == 3)

  # weight= attaches the expression as the native Biogeme WESML formula; it
  # does not create a second likelihood or weighting implementation in R.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability,
    weight = weight
  )
}

# 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.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b02_weight"
)

# estimate() always performs fresh native estimation. Remove only exact b02
# artifacts so a reused output directory cannot silently recycle an old YAML
# or iteration file.
stale_files <- c("b02_weight.yaml", "__b02_weight.iter", "b02_weight.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_b02_weight_model(database)
control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b02_weight",
  user_notes = paste0(
    "Example of a logit model with three alternatives: Train, Car and ",
    "Swissmetro. Weighted Exogenous Sample Maximum Likelihood estimator (WESML)"
  ),
  generate_html = TRUE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

# The complete expression graph is compiled once. Native Python Biogeme then
# performs the weighted estimation and all requested post-estimation work.
fit <- estimate(model, model_name = "b02_weight", control = control)

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