inst/examples/bayesian_swissmetro/plot_b09_nested.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b09. Bayesian nested logit model.
#
# Train and Car share the non-trivial nest "existing". Swissmetro is the
# remaining trivial nest generated by the native NestsForNestedLogit object.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It prepares data and native Python configuration from --data/--python/--output.
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_b09_nested_model <- function(database) {
  # The parameter names and native bounds are preserved 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, upper = 0)
  b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
  nest_parameter <- biogeme_beta(
    "nest_parameter",
    start = 1,
    lower = 1,
    upper = 3
  )

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

  # nested_nest() and nested_nests() are neutral representations of native
  # OneNestForNestedLogit and NestsForNestedLogit objects.
  nests <- nested_nests(
    choice_set = c(1L, 2L, 3L),
    nests = list(nested_nest(nest_parameter, c(1L, 3L), name = "existing"))
  )
  nested_logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability,
    nests = nests,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b09_nested",
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

# Force a new native Bayesian estimation in this directory.
unlink(file.path(prepared$output, c(
  "b09_nested.yaml",
  "b09_nested.nc",
  "b09_nested.html",
  "__b09_nested.iter"
)), force = TRUE)

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

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

# This call delegates correlation construction to native NestsForNestedLogit.
correlation <- nested_logit_correlation(
  model,
  beta_values = coef(fit),
  alternatives_names = c(`1` = "Train", `2` = "Swissmetro", `3` = "Car")
)
print(correlation)
invisible(fit)