advanced-models.R
In rbiogeme: Interface to 'Biogeme'

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE,
  eval = FALSE,
  collapse = TRUE,
  comment = "#>"
)

## ----monte-carlo--------------------------------------------------------------
# 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
#   )
# )

## ----simulated-likelihood-----------------------------------------------------
# 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
#   )
# )

## ----bayesian-----------------------------------------------------------------
# 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)

## ----bayesian-simulation------------------------------------------------------
# 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)

## ----transformations----------------------------------------------------------
# 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")

## ----indexed-selection--------------------------------------------------------
# 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--------------------------------------------------------------------
# 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)

## ----mdcev-forecast-----------------------------------------------------------
# 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-----------------------------------------------------------------
# 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

## ----segmentation-------------------------------------------------------------
# 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-----------------------------------------------------------------
# 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-----------------------------------------------------------------
# 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)

## ----hybrid-choice------------------------------------------------------------
# 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
# )

## ----equivalence-checklist----------------------------------------------------
# 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
# )