inst/examples/hybrid_choice_models/plot_h08_mode_lv_parameter_overrides.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# h08. Build-only hybrid-choice example with an explicit parameter override.
#
# Native Biogeme uses this example to show how a generated measurement loading
# can be replaced by a fixed Beta before estimation. The R version keeps the
# complete specification in this file, stores the override declaratively on a
# generic Biogeme model, and stops before native compilation or estimation.
# Consequently it demonstrates the override without generating Python source
# files or managing Python objects in R.

library(rbiogeme)

# Structural equation for car_centric_attitude. The named draw is compiled by
# native Biogeme only if this build-only model is later used for estimation.
structural_intercept <- biogeme_beta(
  "struct_car_centric_attitude_intercept", start = 0
)
structural_top_manager <- biogeme_beta(
  "struct_car_centric_attitude_top_manager", start = 0
)
structural_car_oriented_parents <- biogeme_beta(
  "struct_car_centric_attitude_car_oriented_parents", start = 0
)
structural_high_education <- biogeme_beta(
  "struct_car_centric_attitude_high_education", start = 0
)
structural_low_education <- biogeme_beta(
  "struct_car_centric_attitude_low_education", start = 0
)
structural_used_to_go_to_school_by_car <- biogeme_beta(
  "struct_car_centric_attitude_used_to_go_to_school_by_car", start = 0
)
structural_sigma_log <- biogeme_beta(
  "struct_car_centric_attitude_sigma_log", start = log(10)
)
car_centric_attitude <- structural_intercept +
  structural_top_manager * variable("top_manager") +
  structural_car_oriented_parents * variable("car_oriented_parents") +
  structural_high_education * variable("high_education") +
  structural_low_education * variable("low_education") +
  structural_used_to_go_to_school_by_car *
    variable("used_to_go_to_school_by_car") +
  exp(structural_sigma_log) *
    draw("struct_car_centric_attitude_draws", "NORMAL_MLHS_ANTI")

# Gaussian measurement equations, with Envir01 serving as the reference
# indicator. Its intercept and loading are fixed by identification; its
# measurement sigma remains a free log-scale parameter, as in native h04.
indicators <- c(
  "Envir01", "Envir02", "Envir06", "Mobil03", "Mobil05", "Mobil08",
  "Mobil09", "Mobil10", "LifSty07"
)
measurement_terms <- vector("list", length(indicators))
names(measurement_terms) <- indicators
measurement_parameters <- vector("list", length(indicators))
names(measurement_parameters) <- indicators
for (indicator_name in indicators) {
  measurement_intercept <- if (indicator_name == "Envir01") {
    0
  } else {
    biogeme_beta(paste0("measurement_intercept_", indicator_name), start = 0)
  }
  measurement_loading <- if (indicator_name == "Envir01") {
    -1
  } else {
    biogeme_beta(
      paste0("measurement_coefficient_car_centric_attitude_", indicator_name),
      start = 0
    )
  }
  measurement_sigma_log <- biogeme_beta(
    paste0("measurement_", indicator_name, "_sigma_log"),
    start = log(10)
  )
  measurement_parameters[[indicator_name]] <- list(
    intercept = measurement_intercept,
    loading = measurement_loading,
    sigma = measurement_sigma_log
  )
  measurement_sigma <- exp(measurement_sigma_log)
  mean_expression <- measurement_intercept + measurement_loading *
    car_centric_attitude
  gaussian_density <- normal_pdf(
    (variable(indicator_name) - mean_expression) / measurement_sigma
  ) / measurement_sigma
  neutral <- (variable(indicator_name) == 6) |
    (variable(indicator_name) == -1)
  measurement_terms[[indicator_name]] <- Elem(
    list(`0` = gaussian_density, `1` = 1), neutral
  )
}

# Preserve the same public parameter discovery order as the simultaneous
# Gaussian model. These neutral terms do not alter the likelihood.
parameter_registration <- 0 * structural_intercept +
  0 * structural_top_manager +
  0 * structural_car_oriented_parents +
  0 * structural_high_education +
  0 * structural_low_education +
  0 * structural_used_to_go_to_school_by_car +
  0 * structural_sigma_log
for (indicator_name in indicators) {
  for (parameter in measurement_parameters[[indicator_name]]) {
    if (inherits(parameter, "biogeme_expression")) parameter_registration <-
      parameter_registration + 0 * parameter
  }
}
conditional_measurement_likelihood <- measurement_terms[[1L]]
for (term in measurement_terms[-1L]) conditional_measurement_likelihood <-
  conditional_measurement_likelihood * term

# Native mode-choice utilities with the latent attitude entering the car
# utility. Alternatives are coded 0 (public transport), 1 (car), and 2 (slow).
choice_beta_cost <- biogeme_beta(
  "choice_beta_cost", start = -1, upper = 0, fixed = TRUE
)
choice_asc_car <- biogeme_beta("choice_asc_car", start = 0)
choice_asc_pt <- biogeme_beta("choice_asc_pt", start = 0)
choice_beta_dist_work <- biogeme_beta(
  "choice_beta_dist_work", start = 0, upper = 0
)
choice_beta_dist_other_purposes <- biogeme_beta(
  "choice_beta_dist_other_purposes", start = 0, upper = 0
)
choice_scale_parameter <- biogeme_beta(
  "choice_scale_parameter", start = 1, lower = 0.0001
)
choice_beta_time_car <- biogeme_beta(
  "choice_beta_time_car", start = 0, upper = 0
)
choice_beta_time_pt <- biogeme_beta(
  "choice_beta_time_pt", start = 0, upper = 0
)
choice_beta_car_centric_attitude_car <- biogeme_beta(
  "choice_beta_car_centric_attitude_car", start = 0
)
work_trip <- variable("PurpHWH") == 1
other_trip_purposes <- variable("PurpHWH") != 1
choice_beta_dist <- choice_beta_dist_work * work_trip +
  choice_beta_dist_other_purposes * other_trip_purposes
v_public_transport <- choice_asc_pt +
  choice_beta_time_pt * variable("TimePT_hour") +
  choice_beta_cost * variable("MarginalCostPT")
v_car <- choice_asc_car +
  choice_beta_time_car * variable("TimeCar_hour") +
  choice_beta_cost * variable("CostCarCHF") +
  choice_beta_car_centric_attitude_car * car_centric_attitude
v_slow_modes <- choice_beta_dist * variable("distance_km")
utilities <- list(
  `0` = choice_scale_parameter * v_public_transport,
  `1` = choice_scale_parameter * v_car,
  `2` = choice_scale_parameter * v_slow_modes
)
availability <- list(
  `0` = 1,
  `1` = variable("car_is_available"),
  `2` = 1
)
conditional_choice_likelihood <- logit_probability(
  utilities = utilities,
  availability = availability,
  alternative = variable("Choice")
)
combined_conditional_likelihood <- conditional_measurement_likelihood *
  conditional_choice_likelihood
log_likelihood <- parameter_registration +
  log(monte_carlo(combined_conditional_likelihood))

# This example has no data in the native version. The generic R model API
# requires a database, so use one symbolic row solely to validate the variable
# names and retain the declarative override. No values are evaluated because
# this script does not compile, simulate, or estimate the model.
symbolic_columns <- c(
  "top_manager", "car_oriented_parents", "high_education",
  "low_education", "used_to_go_to_school_by_car", "Envir01", "Envir02",
  "Envir06", "Mobil03", "Mobil05", "Mobil08", "Mobil09", "Mobil10",
  "LifSty07", "PurpHWH", "TimePT_hour", "MarginalCostPT", "Choice",
  "car_is_available", "TimeCar_hour", "CostCarCHF", "distance_km"
)
symbolic_data <- as.data.frame(
  setNames(lapply(symbolic_columns, function(...) 0), symbolic_columns),
  check.names = FALSE
)
symbolic_database <- biogeme_database("h08_symbolic", symbolic_data)

# Replace one generated measurement loading by a fixed Beta with the same
# native name. The generic model stores this override for the bridge to apply
# at the next native compilation; this build-only example intentionally stops
# before that compilation, matching native h08.
loading_name <- "measurement_coefficient_car_centric_attitude_Envir02"
overrides <- setNames(
  list(biogeme_beta(loading_name, start = 0.5, fixed = TRUE)),
  loading_name
)
unoverridden_model <- biogeme_model(
  symbolic_database,
  formula = log_likelihood
)
before <- names(unoverridden_model$parameters)
model <- biogeme_model(
  symbolic_database,
  formula = log_likelihood,
  parameter_overrides = overrides
)
after <- names(model$parameters)
parameter_table <- biogeme_model_parameters(model)
overridden_loading <- parameter_table[parameter_table$name == loading_name, , drop = FALSE]

cat("Overridden parameter: ", loading_name, "\n", sep = "")
cat("Parameter present before override: ", loading_name %in% before, "\n", sep = "")
cat("Parameter present after override: ", loading_name %in% after, "\n", sep = "")
cat("Overridden parameter fixed: ", overridden_loading$fixed[[1L]], "\n", sep = "")
cat("Number of Betas before override: ", length(unique(before)), "\n", sep = "")
cat("Number of Betas after override: ", length(unique(after)), "\n", sep = "")
cat("Final expression class: ", class(log_likelihood)[[1L]], "\n", sep = "")