plot_b03_scale.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b03. Moneymetric and heteroscedastic specification
#
# This example estimates a group-specific scale. Fixing the cost coefficient
# at -1 expresses utility in CHF (a moneymetric specification), while the two
# scale parameters allow different error variances for GROUP 3 and all other
# observations. Estimation remains entirely in 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_b03_scale_model <- function(database) {
  # Fixing b_cost at -1 normalizes utility in monetary units. The two scale
  # parameters are estimated above the positive lower bound used by Python.
  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 = -1, fixed = TRUE)
  scale_not_group3 <- biogeme_beta(
    "scale_not_group3", start = 1, lower = 0.001
  )
  scale_group3 <- biogeme_beta("scale_group3", start = 1, lower = 0.001)

  # These utility expressions are symbolic trees. variable() does not read
  # the data in R; the full tree is compiled once by the Python bridge.
  v_train <- asc_train + b_time * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED")
  v_swissmetro <- asc_sm + b_time * variable("SM_TT_SCALED") +
    b_cost * variable("SM_COST_SCALED")
  v_car <- asc_car + b_time * variable("CAR_TT_SCALED") +
    b_cost * variable("CAR_CO_SCALED")

  # GROUP is retained as a raw database column by swissmetro_data(). The
  # comparisons are native 0/1 indicator expressions, so scale is selected
  # per observation inside the native likelihood evaluation.
  scale <- (variable("GROUP") != 3) * scale_not_group3 +
    (variable("GROUP") == 3) * scale_group3

  # Multiplying every utility by scale is the heteroscedastic specification;
  # it is not a separate R-side likelihood calculation.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = scale * v_train,
      `2` = scale * v_swissmetro,
      `3` = scale * v_car
    ),
    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 work from any current working directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b03_scale"
)

# estimate() always performs fresh native estimation. Remove only exact b03
# artifacts so a reused output directory cannot silently recycle an old YAML
# or iteration file.
stale_files <- c("b03_scale.yaml", "__b03_scale.iter", "b03_scale.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_b03_scale_model(database)
control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b03_scale",
  user_notes = paste0(
    "Illustrates a moneymetric and heteroscedastic specification. A different ",
    "scale is associated with different segments of the sample. The utility ",
    "function is expressed in CHF."
  ),
  generate_html = TRUE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

# The complete expression graph is compiled once. Native Python Biogeme then
# performs estimation and produces the reported diagnostics.
fit <- estimate(model, model_name = "b03_scale", control = control)

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