inst/examples/mdcev_no_outside_good/plot_translated_estimation.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# Translated MDCEV estimation.
#
# This script is the self-contained R counterpart of
# plot_translated_estimation.py. It contains the complete translated-utility
# model specification. rbiogeme compiles the symbolic expression tree once;
# native Python Biogeme performs the likelihood evaluation, differentiation,
# optimization, and reporting.

library(rbiogeme)

# prepare_mdcev_example() is defined in example_utils.R next to this script.
# It reads data.csv by default, accepts --data=, --python=, and --output=,
# configures the selected Python Biogeme environment, and creates a fresh
# output directory so an old result cannot affect this run.
script_arguments <- commandArgs(trailingOnly = FALSE)
script_file <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_file) != 1L) {
  stop("This example must be run as an R script.", call. = FALSE)
}
script_file <- sub("^--file=", "", script_file)
source(file.path(dirname(normalizePath(script_file)), "example_utils.R"))

prepared <- prepare_mdcev_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "translated_estimation"
)

# The native example uses a fresh estimate. These are the native filenames
# that could otherwise be mistaken for reusable optimization output.
stale_files <- c("translated.yaml", "translated.html", "__translated.iter")
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)

database <- prepare_mdcev_database(prepared$data)

# WESML uses the exact multiplier in the native specification.py file.
weight <- variable("weight") * 1.7718243289995812

# Baseline-utility parameters. Names and order match native Biogeme.
cte_shopping <- biogeme_beta("cte_shopping", start = 0)
cte_socializing <- biogeme_beta("cte_socializing", start = 0)
cte_recreation <- biogeme_beta("cte_recreation", start = 0)
number_members_socializing <- biogeme_beta("number_members_socializing", start = 0)
number_members_recreation <- biogeme_beta("number_members_recreation", start = 0)
metropolitan_shopping <- biogeme_beta("metropolitan_shopping", start = 0)
male_shopping <- biogeme_beta("male_shopping", start = 0)
male_socializing <- biogeme_beta("male_socializing", start = 0)
male_recreation <- biogeme_beta("male_recreation", start = 0)
age_15_40_shopping <- biogeme_beta("age_15_40_shopping", start = 0)
age_15_40_recreation <- biogeme_beta("age_15_40_recreation", start = 0)
age_41_60_socializing <- biogeme_beta("age_41_60_socializing", start = 0)
age_41_60_personal <- biogeme_beta("age_41_60_personal", start = 0)
bachelor_socializing <- biogeme_beta("bachelor_socializing", start = 0)
bachelor_personal <- biogeme_beta("bachelor_personal", start = 0)
white_personal <- biogeme_beta("white_personal", start = 0)
spouse_shopping <- biogeme_beta("spouse_shopping", start = 0)
spouse_recreation <- biogeme_beta("spouse_recreation", start = 0)
employed_shopping <- biogeme_beta("employed_shopping", start = 0)
sunday_socializing <- biogeme_beta("sunday_socializing", start = 0)
sunday_personal <- biogeme_beta("sunday_personal", start = 0)

# Complete baseline utility equations. variable() creates native-bound
# symbolic data variables; it does not evaluate a column in R.
shopping <- cte_shopping +
  metropolitan_shopping * variable("metro") +
  male_shopping * variable("male") +
  age_15_40_shopping * variable("age15_40") +
  spouse_shopping * variable("spousepr") +
  employed_shopping * variable("employed")

socializing <- cte_socializing +
  number_members_socializing * variable("hhsize") +
  male_socializing * variable("male") +
  age_41_60_socializing * variable("age41_60") +
  bachelor_socializing * variable("bachigher") +
  sunday_socializing * variable("Sunday")

recreation <- cte_recreation +
  number_members_recreation * variable("hhsize") +
  male_recreation * variable("male") +
  age_15_40_recreation * variable("age15_40") +
  spouse_recreation * variable("spousepr")

personal <- age_41_60_personal * variable("age41_60") +
  bachelor_personal * variable("bachigher") +
  white_personal * variable("white") +
  sunday_personal * variable("Sunday")

baseline_utilities <- list(
  `1` = shopping,
  `2` = socializing,
  `3` = recreation,
  `4` = personal
)

# Observed quantities use the same symbolic t1:t4 / 60 transformation as the
# native specification. number_chosen is already present in data.csv.
consumed_quantities <- list(
  `1` = variable("t1") / 60,
  `2` = variable("t2") / 60,
  `3` = variable("t3") / 60,
  `4` = variable("t4") / 60
)

# Translated MDCEV shape parameters. Gamma is strictly positive, while alpha
# is constrained to the open interval (0, 1), exactly as in the native file.
lowest_positive_value <- 0.0001
gamma_parameters <- list(
  `1` = biogeme_beta("gamma_shopping", start = 1, lower = lowest_positive_value),
  `2` = biogeme_beta("gamma_socializing", start = 1, lower = lowest_positive_value),
  `3` = biogeme_beta("gamma_recreation", start = 1, lower = lowest_positive_value),
  `4` = biogeme_beta("gamma_personal", start = 1, lower = lowest_positive_value)
)
alpha_parameters <- list(
  `1` = biogeme_beta("alpha_shopping", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `2` = biogeme_beta("alpha_socializing", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `3` = biogeme_beta("alpha_recreation", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `4` = biogeme_beta("alpha_personal", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value)
)
scale_parameter <- biogeme_beta("scale", start = 1, lower = lowest_positive_value)

# This explicit object selects native Translated. Its loglikelihood() method
# constructs the MDCEV likelihood on the Python side; no R likelihood is
# evaluated here.
model <- biogeme_mdcev_model(
  database = database,
  model_type = "translated",
  baseline_utilities = baseline_utilities,
  gamma_parameters = gamma_parameters,
  alpha_parameters = alpha_parameters,
  scale_parameter = scale_parameter,
  weights = weight,
  number_of_chosen_alternatives = variable("number_chosen"),
  consumed_quantities = consumed_quantities
)

fit <- mdcev_estimate(
  model = model,
  model_name = "translated",
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "translated",
    tolerance = 0.0004,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)

# These are the native Biogeme post-estimation operations used by the Python
# example. The returned values are ordinary R text/data frames, not Python
# objects.
cat(mdcev_short_summary(fit), "\n")
print(mdcev_parameter_table(fit)[["Estimated parameters"]])
invisible(fit)