plot_b03cnl.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# Cross-nested-logit estimation using native sampling of alternatives.
# The full utility and CNL allocation specification is intentionally visible.

library(rbiogeme)

# prepare_sampling_example() is defined in example_utils.R. It provides only
# command-line data/output handling; no model specification is hidden there.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "example_utils.R"))

prepared <- prepare_sampling_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "cnl_10_63"
)
alternatives <- prepared$alternatives
observations <- prepared$observations

# Main sampling uses the downtown partition; the MEV sample contains the union
# of Asian and downtown alternatives (63 alternatives in the native example).
all_alternatives <- sort(unique(as.integer(alternatives$ID)))
asian <- all_alternatives[alternatives$Asian[match(all_alternatives, alternatives$ID)] == 1]
downtown <- all_alternatives[alternatives$downtown[match(all_alternatives, alternatives$ID)] == 1]
asian_and_downtown <- intersect(asian, downtown)
only_asian <- setdiff(asian, downtown)
only_downtown <- setdiff(downtown, asian)
asian_or_downtown <- union(asian, downtown)
partition <- biogeme_sampling_partition(
  segments = list(downtown, setdiff(all_alternatives, downtown)),
  sample_sizes = sampling_segment_sizes(10, 2),
  full_set = all_alternatives
)
mev_partition <- biogeme_sampling_partition(
  segments = list(asian_or_downtown),
  sample_sizes = 63,
  full_set = asian_or_downtown
)

log_dist <- cross_variable(
  "log_dist",
  log(((variable("user_lat") - variable("rest_lat"))^2 +
    (variable("user_lon") - variable("rest_lon"))^2)^0.5)
)

# Complete utility specification, with exact native parameter names.
beta_rating <- biogeme_beta("beta_rating", start = 0)
beta_price <- biogeme_beta("beta_price", start = 0)
beta_chinese <- biogeme_beta("beta_chinese", start = 0)
beta_japanese <- biogeme_beta("beta_japanese", start = 0)
beta_korean <- biogeme_beta("beta_korean", start = 0)
beta_indian <- biogeme_beta("beta_indian", start = 0)
beta_french <- biogeme_beta("beta_french", start = 0)
beta_mexican <- biogeme_beta("beta_mexican", start = 0)
beta_lebanese <- biogeme_beta("beta_lebanese", start = 0)
beta_ethiopian <- biogeme_beta("beta_ethiopian", start = 0)
beta_log_dist <- biogeme_beta("beta_log_dist", start = 0)
utility <- beta_rating * variable("rating") +
  beta_price * variable("price") +
  beta_chinese * variable("category_Chinese") +
  beta_japanese * variable("category_Japanese") +
  beta_korean * variable("category_Korean") +
  beta_indian * variable("category_Indian") +
  beta_french * variable("category_French") +
  beta_mexican * variable("category_Mexican") +
  beta_lebanese * variable("category_Lebanese") +
  beta_ethiopian * variable("category_Ethiopian") +
  beta_log_dist * variable("log_dist")

# In the CNL nests, overlapping alternatives receive alpha = 0.5 and
# alternatives belonging to only one nest receive alpha = 1. Sparse storage
# mirrors native dict_of_alpha mappings.
mu_downtown <- biogeme_beta("mu_downtown", start = 1, lower = 1)
mu_asian <- biogeme_beta("mu_asian", start = 1, lower = 1)
downtown_allocation <- c(
  setNames(rep(0.5, length(asian_and_downtown)), as.character(asian_and_downtown)),
  setNames(rep(1, length(only_downtown)), as.character(only_downtown))
)
asian_allocation <- c(
  setNames(rep(0.5, length(asian_and_downtown)), as.character(asian_and_downtown)),
  setNames(rep(1, length(only_asian)), as.character(only_asian))
)
cnl_nests <- cross_nested_nests(
  choice_set = all_alternatives,
  nests = list(
    cross_nested_nest(mu_downtown, as.list(downtown_allocation), name = "downtown"),
    cross_nested_nest(mu_asian, as.list(asian_allocation), name = "asian")
  ),
  sparse = TRUE
)

control <- biogeme_control(
    output_directory = prepared$output,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)
model <- sampled_alternatives_model(
  alternatives = alternatives,
  individuals = observations,
  choice_column = "cnl_3",
  id_column = "ID",
  utility = utility,
  combined_variables = list(log_dist),
  partition = partition,
  mev_partition = mev_partition,
  mev_sample_sizes = 63,
  nests = cnl_nests,
  biogeme_file_name = file.path(prepared$output, "cnl_10_63.dat"),
  model_type = "cnl",
  control = control
)

fit <- estimate_sampled_alternatives(
  model,
  model_name = "cnl_10_63",
  control = control
)
print(fit)
comparison <- compare_sampling_parameters(fit, sampling_true_parameters())
print(comparison$data)
print(comparison$message)

invisible(fit)