inst/examples/bayesian_swissmetro/plot_b11_cnl.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b11. Bayesian cross-nested logit model.
#
# Train is shared between the existing-modes and public-transport nests. Car
# belongs only to the existing nest, while Swissmetro belongs only to the
# public nest. The full CNL expression is specified here and compiled once to
# native Biogeme.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles --data, --python, and --output so this file runs from any working
# directory while keeping the model specification visible below.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

build_b11_cnl_model <- function(database) {
  # Parameter names, starts, bounds, and fixed ASC match native Bayesian b11.
  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_swissmetro <- biogeme_beta("b_time_swissmetro", start = 0, upper = 0)
  b_time_train <- biogeme_beta("b_time_train", start = 0, upper = 0)
  b_time_car <- biogeme_beta("b_time_car", start = 0, upper = 0)
  b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
  b_headway_swissmetro <- biogeme_beta(
    "b_headway_swissmetro", start = 0, upper = 0
  )
  b_headway_train <- biogeme_beta("b_headway_train", start = 0, upper = 0)
  ga_train <- biogeme_beta("ga_train", start = 0)
  ga_swissmetro <- biogeme_beta("ga_swissmetro", start = 0)

  existing_nest_parameter <- biogeme_beta(
    "existing_nest_parameter", start = 1.05, lower = 1, upper = 3
  )
  public_nest_parameter <- biogeme_beta(
    "public_nest_parameter", start = 1.05, lower = 1, upper = 3
  )
  alpha_existing <- biogeme_beta(
    "alpha_existing", start = 0.5, lower = 0, upper = 1
  )
  alpha_public <- 1 - alpha_existing

  utilities <- list(
    `1` = asc_train + b_time_train * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED") +
      b_headway_train * variable("TRAIN_HE") +
      ga_train * variable("GA"),
    `2` = asc_sm + b_time_swissmetro * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED") +
      b_headway_swissmetro * variable("SM_HE") +
      ga_swissmetro * variable("GA"),
    `3` = asc_car + b_time_car * 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")
  )

  # cross_nested_nest() stores the native allocation expressions. Train's two
  # allocations sum to one through alpha_public = 1 - alpha_existing.
  nests <- cross_nested_nests(
    choice_set = c(1L, 2L, 3L),
    nests = list(
      cross_nested_nest(
        nest_parameter = existing_nest_parameter,
        allocation = list(`1` = alpha_existing, `2` = 0, `3` = 1),
        name = "existing"
      ),
      cross_nested_nest(
        nest_parameter = public_nest_parameter,
        allocation = list(`1` = alpha_public, `2` = 1, `3` = 0),
        name = "public"
      )
    )
  )

  cross_nested_logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability,
    nests = nests,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b11_cnl",
      chains = 4,
      bayesian_draws = 4000,
      warmup = 4000,
      calculate_loo = FALSE,
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

# Remove exact stale Bayesian and iteration artifacts before compiling/running.
unlink(file.path(prepared$output, c(
  "b11_cnl.yaml",
  "b11_cnl.nc",
  "b11_cnl.html",
  "__b11_cnl.iter"
)), force = TRUE)

database <- swissmetro_data(prepared$data)
model <- build_b11_cnl_model(database)
fit <- bayesian_estimate(model, model_name = "b11_cnl", control = model$control)

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