inst/examples/swissmetro/plot_b16_panel_discrete_socio_eco.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b16. Discrete mixture with panel data and socioeconomic class membership
#
# This example mirrors plot_b16_panel_discrete_socio_eco.py. It uses the same
# two-class panel mixture as b15a, but lets native class membership depend on
# the INCOME data variable.

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,
# --output, --draws, and --seed 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_b16_panel_discrete_socio_eco_model <- function(database, number_of_draws, seed) {
  number_of_classes <- 2L
  classes <- seq_len(number_of_classes) - 1L

  # Class-specific random coefficients and constants. These names, starts, and
  # fixed Swissmetro ASC normalizations match the native 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) {
      location + scale * draw(paste0("b_time_rnd_class", class), "NORMAL_ANTI")
    },
    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) {
      location + scale * draw(paste0("asc_car_rnd_class", class), "NORMAL_ANTI")
    },
    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) {
      location + scale * draw(paste0("asc_train_rnd_class", class), "NORMAL_ANTI")
    },
    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) {
      location + scale * draw(paste0("asc_sm_rnd_class", class), "NORMAL_ANTI")
    },
    classes,
    asc_sm,
    asc_sm_s
  )

  # Class 0 has no time coefficient, as in b15a and the native example.
  b_time_rnd[[1L]] <- 0

  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")
  )
  trajectory_probabilities <- lapply(utilities, function(class_utilities) {
    panel_likelihood_trajectory(logit_probability(
      utilities = class_utilities,
      availability = availability,
      alternative = variable("CHOICE")
    ))
  })

  # The class score is a complete native expression. INCOME is not converted
  # or evaluated by R; it is compiled as a Biogeme Variable node.
  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_0 <- logit_probability(
    utilities = list(`0` = score_class_0, `1` = 0),
    alternative = 0
  )
  probability_class_1 <- logit_probability(
    utilities = list(`0` = score_class_0, `1` = 0),
    alternative = 1
  )

  # b16 uses the same native outer Monte Carlo order as b15a.
  conditional_choice_probability <- probability_class_0 * trajectory_probabilities[[1L]] +
    probability_class_1 * trajectory_probabilities[[2L]]
  log_probability <- log(monte_carlo(conditional_choice_probability))

  parameter_groups <- list(
    `Class 0` = c(
      "b_cost_class0",
      "asc_car_class0",
      "asc_car_s_class0",
      "asc_train_class0",
      "asc_train_s_class0",
      "asc_sm_s_class0"
    ),
    `Class 1` = c(
      "b_cost_class1",
      "b_time_class1",
      "b_time_s_class1",
      "asc_car_class1",
      "asc_car_s_class1",
      "asc_train_class1",
      "asc_train_s_class1",
      "asc_sm_s_class1"
    )
  )
  draws <- list(
    biogeme_draws("b_time_rnd_class1", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_car_rnd_class0", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_car_rnd_class1", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_train_rnd_class0", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_train_rnd_class1", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_sm_rnd_class0", "NORMAL_ANTI", number_of_draws, seed),
    biogeme_draws("asc_sm_rnd_class1", "NORMAL_ANTI", number_of_draws, seed)
  )

  biogeme_model(
    database = database,
    formula = log_probability,
    draws = draws,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b16_panel_discrete_socio_eco",
      number_of_draws = number_of_draws,
      seed = seed,
      second_derivatives = "never",
      group_of_parameters = parameter_groups,
      generate_html = TRUE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

number_of_draws <- if (!is.null(prepared$options$draws) && nzchar(prepared$options$draws)) {
  example_integer(prepared$options$draws, "draws")
} else {
  5000L
}
seed <- if (!is.null(prepared$options$seed) && nzchar(prepared$options$seed)) {
  example_integer(prepared$options$seed, "seed")
} else {
  1223L
} 

# Always estimate afresh. Remove only exact b16 artifacts so an old YAML or
# iteration file cannot silently supply the estimates.
stale_files <- c(
  "b16_panel_discrete_socio_eco.yaml",
  "__b16_panel_discrete_socio_eco.iter",
  "b16_panel_discrete_socio_eco.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)

# panel=TRUE declares ID after the native-equivalent Swissmetro filter and
# derived-variable operations. The bridge validates contiguous trajectories.
database <- swissmetro_data(prepared$data, panel = TRUE)
model <- build_b16_panel_discrete_socio_eco_model(database, number_of_draws, seed)

cat(sprintf("Panel identifier: %s\n", database$panel_id))
cat(sprintf("Draws per individual: %d\n", number_of_draws))
cat("Class membership covariate: INCOME\n")

# estimate() delegates the class probabilities, panel products, Monte Carlo
# integration, derivatives, optimization, and reporting to native Biogeme.
fit <- estimate(
  model,
  model_name = "b16_panel_discrete_socio_eco",
  control = model$control
)
print(summary(fit))
print(coef(fit))

invisible(fit)