inst/examples/hybrid_choice_models/plot_h02_lv_mimic_gauss.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# h02. Gaussian MIMIC model with one latent variable.
#
# This example keeps the complete structural and measurement specification in
# the script. It is the R counterpart of plot_h02_lv_mimic_gauss.py. The
# semantic native resolver is represented here by ordinary neutral nodes:
# native Biogeme still receives the final expression graph and performs the
# Gaussian densities, Monte-Carlo integration, differentiation, and estimate.

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

# Structural equation for car_centric_attitude. The parameter names mirror the
# native latent-variable builder and are part of the equivalence contract.
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 native h02 specification uses these nine indicators. Every sigma is
# positive through a log-exp transformation, while Envir01 supplies the
# identification normalization (intercept 0 and loading -1).
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 * car_centric_attitude
  standardized_residual <-
    (variable(indicator_name) - mean_expression) / sigma
  gaussian_density <- normal_pdf(standardized_residual) / sigma

  # Likert labels 6 and -1 are neutral/missing in native Biogeme and must
  # contribute one to a likelihood product. Elem() is symbolic indexed
  # selection; it is not an R vector lookup.
  neutral <- (variable(indicator_name) == 6) |
    (variable(indicator_name) == -1)
  measurement_terms[[indicator_name]] <- Elem(
    list(`0` = gaussian_density, `1` = 1),
    neutral
  )
}

# Native latent_variables builds the parameter table in structural order and
# then, for each indicator, in intercept -> sigma -> loading order. These
# zero-valued registration terms preserve that public result ordering while
# leaving the likelihood unchanged; they are part of the symbolic graph and
# are compiled by Biogeme together with the actual measurement equations.
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
  }
}

# MultipleProduct in native Biogeme is represented by a product tree built in
# R. The resulting tree is compiled once; no R loop runs during estimation.
conditional_likelihood <- measurement_terms[[1L]]
for (term in measurement_terms[-1L]) conditional_likelihood <-
  conditional_likelihood * term
log_likelihood <- parameter_registration + log(monte_carlo(conditional_likelihood))

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

# estimate() is deliberately used instead of estimate_or_load(): a clean
# output directory is required for every run and old YAML/iteration files are
# never silently reused.
fit <- estimate(model, model_name = "plot_h02_lv_mimic_gauss", control = control)
print(summary(fit))
print(coef(fit))
invisible(fit)