inst/examples/swissmetro/plot_b14_nested_endogenous_sampling.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b14. Nested logit with corrections for endogenous sampling
#
# This example mirrors plot_b14_nested_endogenous_sampling.py. The alternative
# sampling corrections and the nested-logit MEV derivatives are compiled into
# one native Biogeme expression; R does not reimplement the correction formula.

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_b14_nested_endogenous_sampling_model <- function(database) {
  # Parameter names, starts, bounds, and fixed ASC match native b14 exactly.
  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)
  nest_parameter <- biogeme_beta(
    "nest_parameter",
    start = 1,
    lower = 1,
    upper = 10
  )

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

  # Train and Car form the non-trivial existing-modes nest. Swissmetro is a
  # native trivial nest because it is not listed in this nest.
  existing <- nested_nest(
    nest_parameter = nest_parameter,
    alternatives = c(1L, 3L),
    name = "existing"
  )
  nests <- nested_nests(
    choice_set = c(1L, 2L, 3L),
    nests = list(existing)
  )

  # These are the log probabilities that a respondent with each choice is
  # included in the sample. They are the exact constants used by Python b14.
  correction <- list(
    `1` = log(4.42e-2),
    `2` = log(3.36e-3),
    `3` = log(7.5e-3)
  )

  # The bridge calls native get_mev_for_nested() and
  # logmev_endogenous_sampling() after compiling this complete expression.
  log_probability <- nested_endogenous_sampling_log_probability(
    utilities = utilities,
    availability = availability,
    nests = nests,
    correction = correction,
    alternative = variable("CHOICE")
  )

  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b14_nested_endogenous_sampling",
      generate_html = TRUE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

# Always estimate from the expression tree. Remove only exact b14 artifacts so
# an old YAML or iteration file cannot silently be recycled.
stale_files <- c(
  "b14_nested_endogenous_sampling.yaml",
  "__b14_nested_endogenous_sampling.iter",
  "b14_nested_endogenous_sampling.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_b14_nested_endogenous_sampling_model(database)

# estimate() delegates the corrected nested likelihood, derivatives,
# optimization, and reporting to native Biogeme.
fit <- estimate(
  model,
  model_name = "b14_nested_endogenous_sampling",
  control = model$control
)
print(summary(fit))
print(coef(fit))

invisible(fit)