inst/examples/hybrid_choice_models/plot_h04_mode_lv_gauss_simult.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# h04. Simultaneous Gaussian hybrid mode-choice model.
#
# This is the R counterpart of plot_h04_mode_lv_gauss_simult.py. The complete
# latent-variable, Gaussian measurement, normalization, and mode-choice
# specification is intentionally visible in this file. All expressions remain
# symbolic in R and are compiled once; native Biogeme performs the joint
# likelihood evaluation, Monte-Carlo integration, derivatives, and estimation.

library(rbiogeme)

script_arguments <- commandArgs(trailingOnly = FALSE)
script_argument <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_argument) != 1L) {
  stop("Run this example as an R script with Rscript.", call. = FALSE)
}
example_directory <- dirname(normalizePath(sub("^--file=", "", script_argument)))
source(file.path(example_directory, "example_utils.R"))
source(file.path(example_directory, "optima.R"))

prepared <- prepare_optima_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "plot_h04_mode_lv_gauss_simult",
  default_data = file.path(example_directory, "optima.dat")
)
database <- optima_database(prepared$data)

# Structural equation for car_centric_attitude. The names are the native
# Biogeme names and are preserved in the result YAML.
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")

# The latent variable is measured by the nine indicators used in the native
# semantic specification. Envir01 is the reference indicator: its intercept is
# fixed at 0 and its loading at -1 for identification. Every Gaussian sigma is
# estimated on the log scale, so exp(sigma_log) is strictly positive.
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_sigma_log <- biogeme_beta(
    paste0("measurement_", indicator_name, "_sigma_log"),
    start = log(10)
  )
  measurement_loading <- if (indicator_name == "Envir01") {
    -1
  } else {
    biogeme_beta(
      paste0("measurement_coefficient_car_centric_attitude_", indicator_name),
      start = 0
    )
  }
  measurement_parameters[[indicator_name]] <- list(
    intercept = measurement_intercept,
    loading = measurement_loading,
    sigma = measurement_sigma_log
  )

  sigma <- exp(measurement_sigma_log)
  mean_expression <- measurement_intercept + measurement_loading *
    car_centric_attitude
  gaussian_density <- normal_pdf(
    (variable(indicator_name) - mean_expression) / sigma
  ) / sigma
  neutral <- (variable(indicator_name) == 6) |
    (variable(indicator_name) == -1)
  measurement_terms[[indicator_name]] <- Elem(
    list(`0` = gaussian_density, `1` = 1), neutral
  )
}

# The native semantic builder declares structural parameters first and then
# measurement parameters in intercept -> loading -> sigma order. The neutral
# zero terms preserve that public ordering without changing 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

# The choice model is the same three-alternative specification as h01, with
# the latent attitude entering the car utility. The common cost coefficient is
# fixed at -1 and the positive scale is estimated.
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
)

# logit_probability() creates the native conditional logit probability. The
# measurement and choice probabilities are multiplied before Monte Carlo
# integration: this is simultaneous, joint estimation.
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))

model <- biogeme_model(database, formula = log_likelihood)
control <- optima_estimation_control(
  model_name = "plot_h04_mode_lv_gauss_simult",
  prepared = prepared,
  numerically_safe = NULL,
  number_of_draws = prepared$number_of_draws
)

# estimate() always runs a fresh native estimation. Existing YAML and
# iteration files are never silently recycled.
fit <- estimate(model, model_name = "plot_h04_mode_lv_gauss_simult", control = control)
print(summary(fit))
print(coef(fit))
invisible(fit)