inst/examples/swissmetro/plot_b15a_panel_discrete.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b15a. Discrete mixture with panel data
#
# This example mirrors plot_b15a_panel_discrete.py. It is a two-class latent
# class model. Each class has normally distributed individual parameters, the
# conditional choice probabilities are multiplied over each panel trajectory,
# and the class mixture is integrated by native Monte Carlo.

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_b15a_panel_discrete_model <- function(database, number_of_draws, seed) {
  number_of_classes <- 2L

  # Each class has its own parameter names, starts, and identification
  # constraints, exactly as in the native Python example.
  b_cost <- lapply(seq_len(number_of_classes) - 1L, function(class) {
    biogeme_beta(paste0("b_cost_class", class), start = 0)
  })
  b_time <- lapply(seq_len(number_of_classes) - 1L, function(class) {
    biogeme_beta(paste0("b_time_class", class), start = 0)
  })
  b_time_s <- lapply(seq_len(number_of_classes) - 1L, 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")
    },
    seq_len(number_of_classes) - 1L,
    b_time,
    b_time_s
  )

  asc_car <- lapply(seq_len(number_of_classes) - 1L, function(class) {
    biogeme_beta(paste0("asc_car_class", class), start = 0)
  })
  asc_car_s <- lapply(seq_len(number_of_classes) - 1L, 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")
    },
    seq_len(number_of_classes) - 1L,
    asc_car,
    asc_car_s
  )

  asc_train <- lapply(seq_len(number_of_classes) - 1L, function(class) {
    biogeme_beta(paste0("asc_train_class", class), start = 0)
  })
  asc_train_s <- lapply(seq_len(number_of_classes) - 1L, 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")
    },
    seq_len(number_of_classes) - 1L,
    asc_train,
    asc_train_s
  )

  asc_sm <- lapply(seq_len(number_of_classes) - 1L, function(class) {
    biogeme_beta(paste0("asc_sm_class", class), start = 0, fixed = TRUE)
  })
  asc_sm_s <- lapply(seq_len(number_of_classes) - 1L, 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")
    },
    seq_len(number_of_classes) - 1L,
    asc_sm,
    asc_sm_s
  )

  # Class 0 has no time coefficient. Replacing its random expression by the
  # literal zero is the native identification restriction.
  b_time_rnd[[1L]] <- 0

  score_class_0 <- biogeme_beta("score_class_0", start = -1.7)
  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
  )

  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 mixture is formed conditional on the draws, then integrated by
  # native Monte Carlo. No trajectory product or integration is evaluated in R.
  choice_probability <- probability_class_0 * trajectory_probabilities[[1L]] +
    probability_class_1 * trajectory_probabilities[[2L]]
  log_probability <- log(monte_carlo(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 = "b15a_panel_discrete",
      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 = "b15a_panel_discrete"
)

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 b15a artifacts so an old YAML or
# iteration file cannot silently supply the estimates.
stale_files <- c(
  "b15a_panel_discrete.yaml",
  "__b15a_panel_discrete.iter",
  "b15a_panel_discrete.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_b15a_panel_discrete_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))

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

invisible(fit)