inst/examples/bayesian_swissmetro/plot_b15_panel_discrete.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b15. Bayesian discrete mixture (latent-class logit) with panel data.
#
# Each class has its own individual-level normal random parameters. The class
# membership probabilities are written explicitly because the native Bayesian
# example does not use the maximum-likelihood logit-model helper here.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It reads --data, --python, and --output, configures native Biogeme, and
# prepares the data without hiding the model specification 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_b15_panel_discrete_model <- function(database) {
  number_of_classes <- 2L
  classes <- seq_len(number_of_classes) - 1L

  # Every class-specific name is retained exactly from native Python.
  b_cost <- lapply(classes, function(class) {
    biogeme_beta(paste0("b_cost_class", class), start = 0)
  })
  b_time <- lapply(classes, function(class) {
    biogeme_beta(paste0("b_time_class", class), start = 0)
  })
  b_time_s <- lapply(classes, function(class) {
    biogeme_beta(paste0("b_time_s_class", class), start = 1)
  })
  b_time_rnd <- Map(
    function(class, location, scale) {
      distributed_parameter(
        paste0("b_time_rnd_class", class),
        location + scale * draw(paste0("b_time_eps_class", class), "NORMAL")
      )
    },
    classes,
    b_time,
    b_time_s
  )

  asc_car <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_car_class", class), start = 0)
  })
  asc_car_s <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_car_s_class", class), start = 1)
  })
  asc_car_rnd <- Map(
    function(class, location, scale) {
      distributed_parameter(
        paste0("asc_car_rnd_class", class),
        location + scale * draw(paste0("asc_car_eps_class", class), "NORMAL")
      )
    },
    classes,
    asc_car,
    asc_car_s
  )

  asc_train <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_train_class", class), start = 0)
  })
  asc_train_s <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_train_s_class", class), start = 1)
  })
  asc_train_rnd <- Map(
    function(class, location, scale) {
      distributed_parameter(
        paste0("asc_train_rnd_class", class),
        location + scale * draw(paste0("asc_train_eps_class", class), "NORMAL")
      )
    },
    classes,
    asc_train,
    asc_train_s
  )

  asc_sm <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_sm_class", class), start = 0, fixed = TRUE)
  })
  asc_sm_s <- lapply(classes, function(class) {
    biogeme_beta(paste0("asc_sm_s_class", class), start = 1)
  })
  asc_sm_rnd <- Map(
    function(class, location, scale) {
      distributed_parameter(
        paste0("asc_sm_rnd_class", class),
        location + scale * draw(paste0("asc_sm_eps_class", class), "NORMAL")
      )
    },
    classes,
    asc_sm,
    asc_sm_s
  )

  # The first class has no time coefficient. This is the native identification
  # restriction; the corresponding random expression is replaced by zero.
  b_time_rnd[[1L]] <- 0

  score_class_0 <- biogeme_beta("score_class_0", start = -1.7)
  probability_class_1 <- 1 / (1 + exp(score_class_0))
  probability_class_0 <- 1 - probability_class_1

  utility_for_class <- function(class_index) {
    list(
      `1` = asc_train_rnd[[class_index]] +
        b_time_rnd[[class_index]] * variable("TRAIN_TT_SCALED") +
        b_cost[[class_index]] * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm_rnd[[class_index]] +
        b_time_rnd[[class_index]] * variable("SM_TT_SCALED") +
        b_cost[[class_index]] * variable("SM_COST_SCALED"),
      `3` = asc_car_rnd[[class_index]] +
        b_time_rnd[[class_index]] * variable("CAR_TT_SCALED") +
        b_cost[[class_index]] * variable("CAR_CO_SCALED")
    )
  }
  utilities <- lapply(seq_len(number_of_classes), utility_for_class)
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )
  conditional_probability_per_class <- lapply(utilities, function(class_utilities) {
    logit_probability(
      utilities = class_utilities,
      availability = availability,
      alternative = variable("CHOICE")
    )
  })
  conditional_choice_probability <- probability_class_0 * conditional_probability_per_class[[1L]] +
    probability_class_1 * conditional_probability_per_class[[2L]]

  # Native panel Bayesian estimation aggregates this per-observation log
  # mixture over each ID; no R-side product or integration is performed.
  biogeme_model(
    database = database,
    formula = log(conditional_choice_probability),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b15_panel_discrete",
      warmup = 4000,
      bayesian_draws = 4000,
      chains = 4,
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

unlink(file.path(prepared$output, c(
  "b15_panel_discrete.yaml",
  "b15_panel_discrete.nc",
  "b15_panel_discrete.html",
  "__b15_panel_discrete.iter"
)), force = TRUE)

# panel = TRUE declares ID and checks that each individual's observations are
# contiguous before native Biogeme performs panel Bayesian estimation.
database <- swissmetro_data(prepared$data, panel = TRUE)
model <- build_b15_panel_discrete_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b15_panel_discrete",
  control = model$control
)

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