inst/examples/hybrid_choice_models/plot_h03_mode_lv_gauss_seq.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# h03. Sequential estimation of a mode-choice model with a Gaussian latent
# variable.
#
# The first stage estimates the h02 Gaussian MIMIC model. Its estimated
# structural parameters are then read from a native YAML result and inserted
# as fixed numeric values in the latent expression below. Only the choice
# parameters are estimated in the second stage, exactly as in the native
# example. With no --mimic-results option, the h02 stage is run automatically
# in a fresh subdirectory so this script is runnable from a clean directory.

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_h03_mode_lv_gauss_seq",
  default_data = file.path(example_directory, "optima.dat")
)
database <- optima_database(prepared$data)

# h02 is repeated here intentionally. An example file can therefore be read
# on its own: the latent structural equation, Gaussian measurement equations,
# identification normalization, and native expression syntax are all visible.
build_h02_mimic_model <- function(database) {
  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)
  )
  latent <- 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")

  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,
      sigma = measurement_sigma_log,
      loading = measurement_loading
    )
    sigma <- exp(measurement_sigma_log)
    mean_expression <- measurement_intercept + measurement_loading * latent
    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
    )
  }
  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_likelihood <- measurement_terms[[1L]]
  for (term in measurement_terms[-1L]) conditional_likelihood <-
    conditional_likelihood * term
  biogeme_model(
    database,
    formula = parameter_registration + log(monte_carlo(conditional_likelihood))
  )
}

requested_mimic_results <- prepared$options$mimic_results
if (!is.null(requested_mimic_results) && nzchar(requested_mimic_results)) {
  mimic_yaml <- normalizePath(requested_mimic_results, mustWork = TRUE)
} else {
  # The nested stage is created only after the final output directory has been
  # checked empty by prepare_optima_example(). It cannot accidentally pick up
  # a previous h02 YAML file.
  mimic_directory <- file.path(prepared$output, "mimic_stage")
  dir.create(mimic_directory, recursive = TRUE, showWarnings = FALSE)
  mimic_model <- build_h02_mimic_model(database)
  mimic_yaml <- file.path(mimic_directory, "plot_h02_lv_mimic_gauss.yaml")
  mimic_control <- optima_estimation_control(
    model_name = "plot_h02_lv_mimic_gauss",
    prepared = prepared,
    numerically_safe = FALSE,
    generate_html = FALSE,
    number_of_draws = prepared$number_of_draws
  )
  invisible(estimate(
    mimic_model,
    model_name = "plot_h02_lv_mimic_gauss",
    control = mimic_control,
    yaml_file_name = mimic_yaml
  ))
  if (!file.exists(mimic_yaml)) {
    stop("The h02 prerequisite did not produce its YAML result: ", mimic_yaml, call. = FALSE)
  }
}

# read_results() converts the standard native YAML into ordinary R values.
# The resulting coefficients are fixed estimates in this sequential stage,
# not Python objects and not re-estimated h02 parameters.
estimated_parameters <- read_results(mimic_yaml)$beta_values

# Reconstruct the sequential latent variable. The native h03 example uses the
# five structural explanatory variables but omits the MIMIC intercept in this
# second-stage expression. Its structural disturbance uses a new named draw.
structural_variables <- c(
  "top_manager", "car_oriented_parents", "high_education", "low_education",
  "used_to_go_to_school_by_car"
)
latent_terms <- lapply(structural_variables, function(name) {
  coefficient_name <- paste0("struct_car_centric_attitude_", name)
  variable(name) * as.numeric(estimated_parameters[[coefficient_name]])
})
latent_deterministic <- latent_terms[[1L]]
for (term in latent_terms[-1L]) latent_deterministic <- latent_deterministic + term
structural_sigma <- exp(
  as.numeric(estimated_parameters[["struct_car_centric_attitude_sigma_log"]])
)
car_centric_attitude <- latent_deterministic +
  structural_sigma * draw("car_centric_attitude_draw", "NORMAL_MLHS_ANTI")

# Choice model: this repeats h01's utility specification and adds the native
# car/attitude interaction to alternative 1. The cost coefficient remains
# fixed at -1 and the utility scale remains positive.
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
)

# h03 uses logit() (probability), then MonteCarlo(), then log(). This differs
# from h01's direct loglogit likelihood because the latent variable is sampled.
conditional_likelihood <- logit_probability(
  utilities = utilities,
  availability = availability,
  alternative = variable("Choice")
)
log_likelihood <- log(monte_carlo(conditional_likelihood))
model <- biogeme_model(database, formula = log_likelihood)
control <- optima_estimation_control(
  model_name = "plot_h03_mode_lv_gauss_seq",
  prepared = prepared,
  numerically_safe = TRUE,
  number_of_draws = prepared$number_of_draws
)
fit <- estimate(model, model_name = "plot_h03_mode_lv_gauss_seq", control = control)
print(summary(fit))
print(coef(fit))
invisible(fit)