Advanced models with rbiogeme

rbiogeme contributors

The advanced interfaces use the same rule as the core interface: R stores a neutral, data-only specification and the bridge compiles it into native Biogeme once. Likelihoods, priors, integrations, optimizers, posterior sampling, and reports remain native operations.

The chunks are intentionally not evaluated during package documentation builds: they are complete, runnable examples, but advanced estimation and simulation require a configured Python environment and may create native output files.

Monte Carlo integration and named draws

draw() creates a named draw node inside an expression. biogeme_draws() adds draw metadata such as the draw count, type, seed, or a supplied matrix. A model with a simulation formula can be used to inspect an integral without writing an R callback.

library(rbiogeme)

database <- biogeme_database(
  "integral_demo",
  data.frame(x = c(-1, 0, 1), id = c(1, 2, 3))
)

z <- random_variable("z")
integrand <- exp(-0.5 * z * z) / sqrt(2 * pi)
quadrature_integral <- integrate_normal(
  expression = integrand,
  name = "z",
  number_of_quadrature_points = 30L
)

u <- draw("u", "UNIFORM_HALTON2")
draw_integral <- monte_carlo(u * variable("x") + 1)

model <- biogeme_model(
  database = database,
  simulations = list(
    quadrature = quadrature_integral,
    monte_carlo = draw_integral
  ),
  draws = biogeme_draws(
    name = "u",
    draw_type = "UNIFORM_HALTON2",
    number_of_draws = 512L,
    seed = 1234L
  )
)

simulate(
  model,
  beta = numeric(),
  control = biogeme_control(
    output_directory = tempfile("rbiogeme-mc-"),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)

For simulation-based likelihoods, place the Monte Carlo average inside the likelihood formula, for example:

b <- biogeme_beta("b", start = 0)
conditional_probability <- exp(b * variable("x") * draw("x_draw", "NORMAL"))

simulated_likelihood <- biogeme_model(
  database = database,
  formula = log(monte_carlo(conditional_probability)),
  draws = biogeme_draws(
    name = "x_draw",
    draw_type = "NORMAL_ANTI",
    number_of_draws = 128L,
    seed = 1223L
  )
)

random_variable() and integrate_normal() describe native quadrature. monte_carlo() describes a native draw average. Random-draw results can vary according to the native draw design and seed; deterministic models should match native results up to floating-point precision.

Bayesian estimation

Bayesian priors are data-only descriptors. They do not contain R functions or callbacks. Attach one to a parameter with prior = biogeme_prior(...) and use bayesian_estimate() with native Bayesian controls.

data <- data.frame(
  choice = c(1, 2, 1, 2, 1, 2),
  x = c(1, 2, 1, 3, 2, 1)
)
database <- biogeme_database("bayesian_demo", data)

asc_2 <- biogeme_beta(
  "asc_2",
  start = 0,
  prior = biogeme_prior("normal", sigma = 2)
)
b_x <- biogeme_beta(
  "b_x",
  start = 0,
  prior = biogeme_prior("student_t", sigma = 3, nu = 5)
)

model <- logit_model(
  database,
  choice = "choice",
  utilities = list(
    `1` = 0,
    `2` = asc_2 + b_x * variable("x")
  )
)

bayesian_fit <- bayesian_estimate(
  model,
  model_name = "bayesian_demo",
  controls = list(
    bayesian_draws = 1000L,
    warmup = 500L,
    chains = 2L,
    target_accept = 0.9,
    calculate_likelihood = TRUE,
    calculate_waic = TRUE,
    calculate_loo = TRUE,
    output_directory = tempfile("rbiogeme-bayesian-")
  )
)

summary(bayesian_fit)
coef(bayesian_fit)
bayesian_stored_variables(bayesian_fit)

The Bayesian result contains native posterior summaries and paths to native output files. Posterior draws remain in the NetCDF file managed by native Biogeme; ordinary R code does not need to create or manage PyMC or ArviZ objects.

Posterior simulation uses the same model and a named simulation list:

simulation_model <- biogeme_model(
  database = database,
  formula = logit_log_probability(
    utilities = list(`1` = 0, `2` = asc_2 + b_x * variable("x")),
    alternative = variable("choice")
  ),
  simulations = list(
    probability = logit_probability(
      utilities = list(`1` = 0, `2` = asc_2 + b_x * variable("x")),
      alternative = variable("choice")
    ),
    marginal_index = b_x * variable("x")
  )
)

posterior_simulation <- simulate_bayesian(
  simulation_model,
  bayesian_results = bayesian_fit,
  percentage_of_draws_to_use = 10,
  lower_quantile = 0.025,
  upper_quantile = 0.975
)
as.data.frame(posterior_simulation)

Piecewise, Box–Cox, and symbolic derivatives

Use expression functions when a native model requires a transformation or a derivative. The result is still a symbolic node and is compiled along with the rest of the model.

x <- variable("x")
lambda <- biogeme_beta("lambda", start = 1)

piecewise_x <- piecewise(
  expression = x,
  thresholds = c(NULL, 0, 10, NULL),
  betas = list(
    biogeme_beta("piece_1", start = 1),
    biogeme_beta("piece_2", start = 1),
    biogeme_beta("piece_3", start = 1)
  ),
  transform = "formula"
)

boxcox_x <- boxcox(x, lambda)
derivative <- derive(boxcox_x, "x")

Elem(mapping, index) selects one expression from a named integer mapping and is useful for category-dependent specifications:

category_beta <- Elem(
  mapping = list(
    `1` = biogeme_beta("b_category_1", start = 0),
    `2` = biogeme_beta("b_category_2", start = 0)
  ),
  index = variable("category")
)

MDCEV estimation and forecasting

MDCEV model constructors describe the native variants explicitly. Baseline utilities, gamma parameters, observed quantities, prices, and the number of chosen alternatives are named by alternative code.

mdcev_data <- data.frame(
  chosen = c(1, 1, 2, 2),
  quantity_1 = c(1, 2, 0, 1),
  quantity_2 = c(0, 1, 2, 1),
  price_1 = c(2, 2, 3, 2),
  price_2 = c(3, 3, 2, 3),
  x = c(1, 2, 1, 3)
)
mdcev_database <- biogeme_database("mdcev_demo", mdcev_data)

gamma_1 <- biogeme_beta("gamma_1", start = 1, lower = 0)
gamma_2 <- biogeme_beta("gamma_2", start = 1, lower = 0)

mdcev_model <- biogeme_mdcev_model(
  database = mdcev_database,
  model_type = "gamma_profile",
  baseline_utilities = list(
    `1` = biogeme_beta("asc_1", start = 0) + variable("x"),
    `2` = biogeme_beta("asc_2", start = 0) + variable("x")
  ),
  gamma_parameters = list(`1` = gamma_1, `2` = gamma_2),
  prices = list(`1` = variable("price_1"), `2` = variable("price_2")),
  number_of_chosen_alternatives = variable("chosen"),
  consumed_quantities = list(
    `1` = variable("quantity_1"),
    `2` = variable("quantity_2")
  )
)

mdcev_fit <- mdcev_estimate(
  mdcev_model,
  model_name = "mdcev_demo",
  control = biogeme_control(
    output_directory = tempfile("rbiogeme-mdcev-"),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)
mdcev_short_summary(mdcev_fit)
mdcev_parameter_table(mdcev_fit)

Forecasting accepts native Gumbel draws supplied by the user or generates them through the native bridge. The two native algorithms can be compared with brute_force = TRUE and brute_force = FALSE:

epsilons <- mdcev_generate_epsilons(
  model = mdcev_model,
  number_of_observations = nrow(mdcev_data),
  number_of_draws = 128L,
  seed = 1234L
)

forecast <- mdcev_forecast(
  model = mdcev_model,
  fit = mdcev_fit,
  total_budget = 10,
  epsilons = epsilons,
  brute_force = FALSE
)
mdcev_forecast_describe(forecast)
mdcev_validate_forecast(
  model = mdcev_model,
  fit = mdcev_fit,
  total_budget = 10,
  epsilons = epsilons
)

Catalogs, segmentation, and assisted specification

Catalogs allow native Biogeme to estimate multiple synchronized specifications without exposing Python controller objects. A catalog is a named mapping of expressions; shared controllers select the same specification in several catalogs.

catalog_controller <- biogeme_catalog_controller(
  name = "time_form",
  specification_names = c("linear", "log")
)

time_catalog <- biogeme_catalog(
  name = "time_catalog",
  expressions = list(
    linear = variable("x"),
    log = log(variable("x") + 1)
  ),
  controller = catalog_controller
)

catalog_database <- biogeme_database(
  "catalog_demo",
  data.frame(
    choice = c(1, 2, 1, 2),
    x = c(1, 2, 3, 4),
    category = c(1, 2, 1, 2)
  )
)

catalog_model <- biogeme_model(
  database = catalog_database,
  formula = logit_log_probability(
    utilities = list(`1` = 0, `2` = time_catalog),
    alternative = variable("choice")
  )
)

catalog_fit <- estimate_catalog(
  catalog_model,
  model_name = "catalog_demo",
  force = TRUE,
  control = biogeme_control(output_directory = tempfile("rbiogeme-catalog-"))
)
catalog_fit$summary
catalog_fit$non_dominated

For discrete segmentation, define a mapping from observed values to segment names and use segment_beta() to generate the segmented parameter. A reference segment is omitted from the parameterization when supplied.

segmentation <- biogeme_segmentation(
  variable = variable("category"),
  mapping = c(`1` = "low", `2` = "high"),
  reference = "low"
)

b_income <- biogeme_beta("b_income", start = 0)
segmented_income <- segment_beta(
  beta = b_income,
  segmentations = list(segmentation)
)

assisted_specification() combines native catalog search, objective values, validity rules, Pareto persistence, and final estimation. The Pareto file is a deliberate checkpoint, so give it a fresh path for a reproducible search or set force = FALSE when resuming is the intended behavior.

assisted_fit <- assisted_specification(
  model = catalog_model,
  objectives = "loglikelihood_dimension",
  validity = NULL,
  pareto_file_name = file.path(tempdir(), "rbiogeme-assisted.pareto"),
  model_name = "assisted_demo",
  force = TRUE,
  control = biogeme_control(output_directory = tempfile("rbiogeme-assisted-"))
)
assisted_fit$summary

Sampling of alternatives

Sampling-of-alternatives models use two ordinary R data frames: one for the alternative universe and one for individuals. A partition describes strata and sample sizes. Native Biogeme regenerates the sampled choice sets and constructs the sampled likelihood; the R interface does not sample in an R callback.

alternatives <- data.frame(
  alternative_id = 1:6,
  alt_time = c(5, 7, 6, 8, 9, 4),
  alt_cost = c(2, 3, 4, 2, 5, 3),
  nest = c(1, 1, 1, 2, 2, 2)
)
individuals <- data.frame(
  choice = c(1, 5, 3, 6),
  income = c(1, 2, 1, 3)
)

partition <- biogeme_sampling_partition(
  segments = list(c(1, 2, 3), c(4, 5, 6)),
  sample_sizes = c(2L, 2L)
)

b_alt_time <- biogeme_beta("b_alt_time", start = 0)
b_alt_cost <- biogeme_beta("b_alt_cost", start = 0)
sampled_model <- sampled_alternatives_model(
  alternatives = alternatives,
  individuals = individuals,
  choice_column = "choice",
  id_column = "alternative_id",
  utility = b_alt_time * variable("alt_time") + b_alt_cost * variable("alt_cost"),
  partition = partition,
  biogeme_file_name = file.path(tempdir(), "rbiogeme-sampled.dat"),
  model_type = "logit",
  control = biogeme_control(output_directory = tempfile("rbiogeme-sampled-"))
)

sampled_fit <- estimate_sampled_alternatives(
  sampled_model,
  model_name = "sampled_demo",
  control = biogeme_control(output_directory = tempfile("rbiogeme-sampled-fit-"))
)
coef(sampled_fit)

cross_variable() describes an individual-by-alternative expression when the utility uses both data frames. sampling_segment_sizes() is a helper for creating nearly equal sample sizes; it does not perform the random sampling.

Hybrid-choice specifications

Hybrid-choice models combine a choice likelihood with measurement equations for latent variables. In R, write the latent-variable expressions and the ordered-logit or ordered-probit measurement likelihood explicitly, then join the components in a generic biogeme_model() formula. The ordered-response constructors use numeric category codes and cutpoint expressions:

hybrid_database <- biogeme_database(
  "hybrid_demo",
  data.frame(
    choice = c(1, 2, 1, 2),
    indicator = c(1, 2, 3, 2),
    x = c(1, 2, 1, 3)
  )
)

latent <- biogeme_beta("latent", start = 0) +
  biogeme_beta("latent_x", start = 0) * variable("x")
cut_1 <- biogeme_beta("cut_1", start = -1)
cut_2 <- biogeme_beta("cut_2", start = 1)

indicator_log_probability <- ordered_probit_log_probability(
  eta = latent,
  cutpoints = list(cut_1, cut_2),
  alternative = variable("indicator"),
  categories = c(1, 2, 3)
)

choice_log_probability <- logit_log_probability(
  utilities = list(`1` = 0, `2` = latent),
  alternative = variable("choice")
)

hybrid_model <- biogeme_model(
  database = hybrid_database,
  formula = choice_log_probability + indicator_log_probability
)

The same pattern supports a simultaneous measurement-and-choice likelihood. When a native example uses sequential estimation or parameter overrides, keep the stages explicit with separate models and use parameter_overrides in the generic constructor rather than changing parameter names.

Equivalence and reproducibility checklist

For a comparison with a native result, keep the following fixed:

equivalence_directory <- tempfile("rbiogeme-equivalence-")
dir.create(equivalence_directory)

equivalence_control <- biogeme_control(
  output_directory = equivalence_directory,
  seed = 1223L,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

Use the same data rows, derived columns, filters, panel ordering, alternative codes, starting values, parameter names, draw type, draw count, and seed. For random-draw or sampling models, document the accepted simulation noise. Never use estimate_or_load() with an implicit or shared path in an equivalence test; use estimate() or set force = TRUE.

The package’s reference pages document each constructor and operation. The guides here explain the R syntax and workflow without requiring users to look up a separate implementation manual.