inst/examples/bayesian_swissmetro/plot_b16_panel_discrete_socio_eco.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b16. Bayesian latent-class logit with panel data and socioeconomic class
# membership. The class probabilities depend on the native INCOME variable.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles --data, --python, and --output and leaves the model syntax below
# fully visible and runnable from any 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)), "..", "swissmetro", "example_utils.R"))

build_b16_panel_discrete_socio_eco_model <- function(database) {
  number_of_classes <- 2L
  classes <- seq_len(number_of_classes) - 1L

  # Class-specific random coefficients and constants use the same native
  # names, starts, and fixed Swissmetro ASCs as the Python example.
  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
  )

  # Class 0 has no time coefficient, matching the native identification
  # restriction in b16_panel_discrete_socio_eco.py.
  b_time_rnd[[1L]] <- 0

  class_cte <- biogeme_beta("class_cte", start = 0)
  class_inc <- biogeme_beta("class_inc", start = 0)
  score_class_0 <- class_cte + class_inc * variable("INCOME")
  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]]

  biogeme_model(
    database = database,
    formula = log(conditional_choice_probability),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b16_panel_discrete_socio_eco",
      warmup = 40,
      bayesian_draws = 40,
      chains = 1,
      generate_html = TRUE,
      generate_yaml = TRUE,
      generate_netcdf = TRUE
    )
  )
}

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

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

# panel = TRUE declares ID and validates contiguous individual trajectories.
database <- swissmetro_data(prepared$data, panel = TRUE)
model <- build_b16_panel_discrete_socio_eco_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b16_panel_discrete_socio_eco",
  control = model$control
)

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