Nothing
inst/doc/apollo-benchmark.R
In rmdcev: Kuhn-Tucker and Multiple Discrete-Continuous Extreme Value Models
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = FALSE # heavy computation — run manually
)
## ----packages-----------------------------------------------------------------
# library(rmdcev)
# library(apollo)
# library(dplyr)
# library(tidyr)
# library(bench) # for timing; install.packages("bench") if needed
## ----fn-to-wide---------------------------------------------------------------
# to_wide <- function(data_long, nalts) {
# data_long <- as.data.frame(data_long)
#
# # Individual-level covariates (b-cols that don't vary within id)
# # b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific
# all_b <- grep("^b[0-9]", names(data_long), value = TRUE)
# # detect individual-specific cols: same value for all rows of an id
# ind_b <- all_b[sapply(all_b, function(v) {
# all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1))
# })]
#
# ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b)))
#
# alt_vars <- data_long %>%
# select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>%
# pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b)))
#
# db <- left_join(ind_vars, alt_vars, by = "id")
#
# # Numeraire quantity
# pj_xj <- Reduce("+", lapply(1:nalts, function(j) {
# db[[paste0("price_", j)]] * db[[paste0("quant_", j)]]
# }))
# db$quant_num <- db$income - pj_xj
# db
# }
## ----fn-apollo-prob-----------------------------------------------------------
# make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"),
# npsi_j = 1, ind_col = "b2") {
# model <- match.arg(model)
# alt_names <- paste0("alt", 1:nalts)
# all_names <- c("outside", alt_names)
#
# function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
#
# # Utility for each non-numeraire alternative
# V <- setNames(vector("list", nalts + 1), all_names)
# V[["outside"]] <- 0
# for (j in seq_len(nalts)) {
# vj <- get(ind_col) * psi_ind # individual-specific component
# for (k in seq_len(npsi_j)) {
# vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k))
# }
# V[[paste0("alt", j)]] <- vj
# }
#
# # Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others
# alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0
# alpha <- setNames(
# c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)),
# all_names
# )
#
# gamma <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))),
# all_names
# )
#
# sigma <- exp(ln_sigma)
#
# cChoice <- setNames(
# c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))),
# all_names
# )
# cost <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))),
# all_names
# )
# avail <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) 1)),
# all_names
# )
#
# mdcev_settings <- list(
# alternatives = all_names,
# budget = income,
# V = V,
# alpha = alpha,
# gamma = gamma,
# sigma = sigma,
# continuousChoice = cChoice,
# cost = cost,
# avail = avail
# )
#
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }
# }
## ----fn-apollo-beta-----------------------------------------------------------
# make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) {
# model <- match.arg(model)
# betas <- c(psi_ind = 0)
# for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0
# for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0
# if (model == "gamma") betas["ln_alpha_out"] <- 0
# betas["ln_sigma"] <- 0
# betas
# }
## ----fn-fit-both--------------------------------------------------------------
# fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") {
# nalts <- sim$data %>% distinct(alt) %>% nrow()
#
# # ---- rmdcev ----
# rmdcev_formula <- as.formula(paste0(
# "~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ")
# ))
# t_rmdcev <- system.time(
# fit_r <- mdcev(
# formula = rmdcev_formula,
# data = sim$data,
# psi_ascs = 0,
# model = model,
# algorithm = "MLE",
# print_iterations = FALSE,
# backend = "rstan"
# )
# )
#
# # ---- Apollo ----
# database <- to_wide(sim$data, nalts)
# db_col <- paste0("b", npsi_j + 1) # individual-specific column
#
# apollo_control <- list(
# modelName = paste0("apollo_", model),
# modelDescr = paste0("Apollo ", model),
# indivID = "id",
# outputDirectory = tempdir(),
# nCores = 1
# )
# apollo_beta <- make_apollo_beta(nalts, model, npsi_j)
# apollo_fixed <- character(0)
#
# apollo_inputs <- apollo_validateInputs(
# apollo_beta = apollo_beta,
# apollo_fixed = apollo_fixed,
# database = database,
# apollo_control = apollo_control
# )
#
# apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col)
#
# t_apollo <- system.time(
# fit_a <- apollo_estimate(
# apollo_beta, apollo_fixed, apollo_prob, apollo_inputs,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L)
# )
# )
#
# # ---- Extract estimates ----
# r_psi <- as.numeric(fit_r$stan_fit$par$psi)
# r_gamma <- as.numeric(fit_r$stan_fit$par$gamma)
# r_scale <- as.numeric(fit_r$stan_fit$par$scale)
#
# a_est <- fit_a$estimate
# a_gamma <- exp(a_est[grep("^ln_g", names(a_est))])
# a_scale <- exp(a_est["ln_sigma"])
#
# # Map column names: psi_j1 = alt-varying, psi_ind = individual-varying
# a_psi_j <- a_est[grep("^psi_j", names(a_est))]
# a_psi_ind <- a_est["psi_ind"]
# a_psi <- c(a_psi_j, a_psi_ind)
#
# list(
# ll = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout),
# psi = rbind(rmdcev = r_psi, apollo = a_psi),
# gamma = rbind(rmdcev = r_gamma[1:min(5, nalts)],
# apollo = a_gamma[1:min(5, nalts)]),
# scale = c(rmdcev = r_scale, apollo = as.numeric(a_scale)),
# time_s = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"])
# )
# }
## ----param-recovery-----------------------------------------------------------
# set.seed(42)
# true_psi_j <- 0.5
# true_psi_ind <- -1.5
# true_gamma <- rep(5, 5)
# true_scale <- 1.0
#
# sim_h0 <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = 500,
# nalts = 5,
# psi_j_parms = true_psi_j,
# psi_i_parms = true_psi_ind,
# gamma_parms = true_gamma,
# scale_parms = true_scale
# )
#
# cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1)
#
# cat("=== hybrid0 parameter comparison ===\n")
# cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
# cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"]))
# cat(sprintf(" LL difference = %.6f\n",
# abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"])))
#
# cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n")
# print(round(cmp_h0$psi, 4))
#
# cat("\nGamma (true: all 5)\n")
# print(round(cmp_h0$gamma, 4))
#
# cat(sprintf("\nScale (true: 1.0): rmdcev=%.4f apollo=%.4f\n",
# cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"]))
## ----param-recovery-gamma-----------------------------------------------------
# set.seed(42)
# true_alpha <- 0.5
#
# sim_g <- GenerateMDCEVData(
# model = "gamma",
# nobs = 500,
# nalts = 5,
# psi_j_parms = true_psi_j,
# psi_i_parms = true_psi_ind,
# gamma_parms = true_gamma,
# alpha_parms = true_alpha,
# scale_parms = true_scale
# )
#
# cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1)
#
# cat("=== gamma parameter comparison ===\n")
# cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
# cmp_g$ll["rmdcev"], cmp_g$ll["apollo"]))
# cat(sprintf(" LL difference = %.6f\n",
# abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"])))
#
# cat("\nPsi:\n"); print(round(cmp_g$psi, 4))
# cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4))
# cat(sprintf("\nScale: rmdcev=%.4f apollo=%.4f\n",
# cmp_g$scale["rmdcev"], cmp_g$scale["apollo"]))
## ----timing-setup-------------------------------------------------------------
# grid <- expand.grid(
# n = c(500, 1000, 5000),
# J = c(50, 100),
# rep = 1:3,
# stringsAsFactors = FALSE
# ) %>%
# filter(!(n == 1000 & J == 100)) # keep n=1000 for J=50 only
## ----timing-mle---------------------------------------------------------------
# timing_mle <- lapply(seq_len(nrow(grid)), function(i) {
# n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i]
# set.seed(r * 100 + i)
#
# sim <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = n,
# nalts = J,
# psi_j_parms = 0.5,
# psi_i_parms = -1.5,
# gamma_parms = rep(5, J),
# scale_parms = 1.0
# )
#
# res <- fit_both(sim, model = "hybrid0", npsi_j = 1)
# data.frame(n=n, J=J, rep=r,
# rmdcev_s = res$time_s["rmdcev"],
# apollo_s = res$time_s["apollo"],
# ll_diff = abs(res$ll["rmdcev"] - res$ll["apollo"]))
# })
#
# timing_mle <- do.call(rbind, timing_mle)
## ----timing-mle-summary-------------------------------------------------------
# summary_mle <- timing_mle %>%
# group_by(n, J) %>%
# summarise(
# rmdcev_mean_s = round(mean(rmdcev_s), 1),
# apollo_mean_s = round(mean(apollo_s), 1),
# speedup = round(mean(apollo_s) / mean(rmdcev_s), 1),
# ll_diff_mean = signif(mean(ll_diff), 3),
# .groups = "drop"
# )
#
# cat("=== MLE estimation timing (seconds, 3 replicates) ===\n")
# print(summary_mle, n = Inf)
## ----sim-apollo-prob5---------------------------------------------------------
# # Concrete apollo_probabilities for J=5 (used only in this section)
# apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
#
# V <- list(
# outside = 0,
# alt1 = psi_ind * b2 + psi_j1 * b1_1,
# alt2 = psi_ind * b2 + psi_j1 * b1_2,
# alt3 = psi_ind * b2 + psi_j1 * b1_3,
# alt4 = psi_ind * b2 + psi_j1 * b1_4,
# alt5 = psi_ind * b2 + psi_j1 * b1_5
# )
# mdcev_settings <- list(
# alternatives = c("outside","alt1","alt2","alt3","alt4","alt5"),
# budget = income, V = V,
# alpha = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0),
# gamma = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3),
# alt4=exp(ln_g4), alt5=exp(ln_g5)),
# sigma = exp(ln_sigma),
# continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2,
# alt3=quant_3, alt4=quant_4, alt5=quant_5),
# cost = list(outside=1, alt1=price_1, alt2=price_2,
# alt3=price_3, alt4=price_4, alt5=price_5),
# avail = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1)
# )
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }
## ----sim-convergence----------------------------------------------------------
# set.seed(42)
# J5 <- 5
# n500 <- 500
#
# sim5 <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = n500,
# nalts = J5,
# psi_j_parms = 0.5,
# psi_i_parms = -1.5,
# gamma_parms = rep(5, J5),
# scale_parms = 1.0
# )
#
# # ── Fit both packages ──────────────────────────────────────────────────────
# fit_r5 <- mdcev(
# formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0,
# model = "hybrid0", algorithm = "MLE",
# print_iterations = FALSE, backend = "rstan"
# )
#
# database5 <- to_wide(sim5$data, J5)
# ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id",
# outputDirectory = tempdir(), nCores = 1)
# ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1)
# ai5 <- suppressMessages(
# apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0),
# database = database5, apollo_control = ac5))
# fit_a5 <- suppressMessages(suppressWarnings(
# apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L))
# ))
#
# # ── 1% price-change policy ─────────────────────────────────────────────────
# baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500)
# price_change <- 0.01 * colMeans(baseline_prices) # 1% of mean price per alt
#
# # rmdcev: additive policy vector (0 for numeraire, then per-alt changes)
# policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE)
# policies5$price_p[[2]] <- c(0, price_change)
#
# df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1)
#
# t_rmdcev5 <- system.time(
# demand_r5 <- mdcev.sim(
# df_sim5$df_indiv, df_common = df_sim5$df_common,
# sim_options = df_sim5$sim_options,
# cond_err = 0, nerrs = 100, sim_type = "demand"
# )
# )
#
# # Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index]
# mean_demand <- function(demand_list, policy_idx, alt_idx) {
# mean(sapply(demand_list, function(ind)
# mean(sapply(ind, function(s) s[policy_idx, alt_idx]))))
# }
# d_base_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a))
# d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a))
# delta_r <- d_policy_r - d_base_r
#
# # Apollo: predict at baseline and policy databases
# # Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt)
# database5_pol <- database5
# for (j in 1:J5) {
# database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j]
# }
# # Recompute quant_num at observed quantities under new prices
# pj_xj_pol <- Reduce("+", lapply(1:J5, function(j)
# database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]]))
# database5_pol$quant_num <- database5_pol$income - pj_xj_pol
#
# ai5_pol <- suppressMessages(
# apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0),
# database = database5_pol, apollo_control = ac5))
#
# t_apollo5_base <- system.time(
# pred_base <- suppressMessages(suppressWarnings(
# apollo_prediction(fit_a5, apollo_probabilities_j5, ai5,
# prediction_settings = list(nRep = 100, silent = TRUE))
# ))
# )
# t_apollo5_pol <- system.time(
# pred_pol <- suppressMessages(suppressWarnings(
# apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol,
# prediction_settings = list(nRep = 100, silent = TRUE))
# ))
# )
# t_apollo5 <- t_apollo5_base + t_apollo5_pol
#
# delta_a <- sapply(1:J5, function(j) {
# mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) -
# mean(pred_base[[paste0("alt", j, "_cont_mean")]])
# })
#
# cat("=== Demand change: 1% price increase (n=500, J=5) ===\n")
# cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff"))
# for (j in 1:J5) {
# cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j],
# abs(delta_r[j] - delta_a[j])))
# }
# cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a))))
# cat(sprintf("rmdcev %% change: %s%%\n",
# paste(round(100 * delta_r / d_base_r, 2), collapse = ", ")))
#
# cat(sprintf("\nTiming: rmdcev %.2f s | Apollo (2 predictions) %.2f s\n",
# t_rmdcev5["elapsed"], t_apollo5["elapsed"]))
## ----make-concrete-prob-------------------------------------------------------
# # Dynamically generate a concrete apollo_probabilities function for any J.
# # apollo_prediction() does source-code inspection and cannot use factory closures,
# # so we build and eval() a function body with the right number of alternatives.
# make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") {
# all_str <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")")
# V_items <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col,
# " + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n ")
# alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ")
# gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ")
# choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ")
# cost_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ")
# avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ")
#
# eval(parse(text = sprintf(
# 'function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
# V <- list(outside = 0, %s)
# mdcev_settings <- list(
# alternatives = %s, budget = income, V = V,
# alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma),
# continuousChoice = list(%s), cost = list(%s), avail = list(%s)
# )
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items)))
# }
## ----timing-sim-grid----------------------------------------------------------
# timing_sim <- lapply(c(50, 100), function(J) {
# n <- 500
# set.seed(J)
#
# sim <- GenerateMDCEVData(
# model = "hybrid0", nobs = n, nalts = J,
# psi_j_parms = 0.5, psi_i_parms = -1.5,
# gamma_parms = rep(5, J), scale_parms = 1.0
# )
#
# # ── rmdcev ──────────────────────────────────────────────────────────────
# fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0",
# algorithm = "MLE", print_iterations = FALSE, backend = "rstan")
# bp <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n)
# pc <- 0.01 * colMeans(bp)
# pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE)
# pols$price_p[[2]] <- c(0, pc)
# df_s <- PrepareSimulationData(fit_r, pols, nsims = 1)
#
# rmdcev_times <- sapply(1:3, function(r)
# system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common,
# sim_options = df_s$sim_options,
# cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"])
#
# # ── Apollo ──────────────────────────────────────────────────────────────
# database <- to_wide(sim$data, J)
# ac <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id",
# outputDirectory = tempdir(), nCores = 1)
# ab <- make_apollo_beta(J, npsi_j = 1)
# pfn <- make_apollo_prob_concrete(J)
# ai <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0),
# database = database, apollo_control = ac)
# fit_a <- apollo_estimate(ab, character(0), pfn, ai,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L))
#
# # Policy database: 1% additive price increase
# db_pol <- database
# for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j]
# pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j)
# db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]]))
# db_pol$quant_num <- db_pol$income - pj_xj_pol
# ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0),
# database = db_pol, apollo_control = ac)
#
# apollo_times <- sapply(1:3, function(r) {
# t1 <- system.time(apollo_prediction(fit_a, pfn, ai,
# prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
# t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol,
# prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
# t1 + t2
# })
#
# data.frame(n = n, J = J,
# rmdcev_mean = round(mean(rmdcev_times), 2),
# rmdcev_sd = round(sd(rmdcev_times), 2),
# apollo_mean = round(mean(apollo_times), 2),
# apollo_sd = round(sd(apollo_times), 2),
# speedup = round(mean(apollo_times) / mean(rmdcev_times), 1))
# })
#
# timing_sim <- do.call(rbind, timing_sim)
# print(timing_sim, row.names = FALSE)
## ----session------------------------------------------------------------------
# sessionInfo()
Try the rmdcev package in your browser
Any scripts or data that you put into this service are public.
rmdcev documentation built on Sept. 15, 2026, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = FALSE # heavy computation — run manually
)
## ----packages-----------------------------------------------------------------
# library(rmdcev)
# library(apollo)
# library(dplyr)
# library(tidyr)
# library(bench) # for timing; install.packages("bench") if needed
## ----fn-to-wide---------------------------------------------------------------
# to_wide <- function(data_long, nalts) {
# data_long <- as.data.frame(data_long)
#
# # Individual-level covariates (b-cols that don't vary within id)
# # b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific
# all_b <- grep("^b[0-9]", names(data_long), value = TRUE)
# # detect individual-specific cols: same value for all rows of an id
# ind_b <- all_b[sapply(all_b, function(v) {
# all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1))
# })]
#
# ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b)))
#
# alt_vars <- data_long %>%
# select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>%
# pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b)))
#
# db <- left_join(ind_vars, alt_vars, by = "id")
#
# # Numeraire quantity
# pj_xj <- Reduce("+", lapply(1:nalts, function(j) {
# db[[paste0("price_", j)]] * db[[paste0("quant_", j)]]
# }))
# db$quant_num <- db$income - pj_xj
# db
# }
## ----fn-apollo-prob-----------------------------------------------------------
# make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"),
# npsi_j = 1, ind_col = "b2") {
# model <- match.arg(model)
# alt_names <- paste0("alt", 1:nalts)
# all_names <- c("outside", alt_names)
#
# function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
#
# # Utility for each non-numeraire alternative
# V <- setNames(vector("list", nalts + 1), all_names)
# V[["outside"]] <- 0
# for (j in seq_len(nalts)) {
# vj <- get(ind_col) * psi_ind # individual-specific component
# for (k in seq_len(npsi_j)) {
# vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k))
# }
# V[[paste0("alt", j)]] <- vj
# }
#
# # Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others
# alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0
# alpha <- setNames(
# c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)),
# all_names
# )
#
# gamma <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))),
# all_names
# )
#
# sigma <- exp(ln_sigma)
#
# cChoice <- setNames(
# c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))),
# all_names
# )
# cost <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))),
# all_names
# )
# avail <- setNames(
# c(list(1), lapply(seq_len(nalts), function(j) 1)),
# all_names
# )
#
# mdcev_settings <- list(
# alternatives = all_names,
# budget = income,
# V = V,
# alpha = alpha,
# gamma = gamma,
# sigma = sigma,
# continuousChoice = cChoice,
# cost = cost,
# avail = avail
# )
#
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }
# }
## ----fn-apollo-beta-----------------------------------------------------------
# make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) {
# model <- match.arg(model)
# betas <- c(psi_ind = 0)
# for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0
# for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0
# if (model == "gamma") betas["ln_alpha_out"] <- 0
# betas["ln_sigma"] <- 0
# betas
# }
## ----fn-fit-both--------------------------------------------------------------
# fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") {
# nalts <- sim$data %>% distinct(alt) %>% nrow()
#
# # ---- rmdcev ----
# rmdcev_formula <- as.formula(paste0(
# "~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ")
# ))
# t_rmdcev <- system.time(
# fit_r <- mdcev(
# formula = rmdcev_formula,
# data = sim$data,
# psi_ascs = 0,
# model = model,
# algorithm = "MLE",
# print_iterations = FALSE,
# backend = "rstan"
# )
# )
#
# # ---- Apollo ----
# database <- to_wide(sim$data, nalts)
# db_col <- paste0("b", npsi_j + 1) # individual-specific column
#
# apollo_control <- list(
# modelName = paste0("apollo_", model),
# modelDescr = paste0("Apollo ", model),
# indivID = "id",
# outputDirectory = tempdir(),
# nCores = 1
# )
# apollo_beta <- make_apollo_beta(nalts, model, npsi_j)
# apollo_fixed <- character(0)
#
# apollo_inputs <- apollo_validateInputs(
# apollo_beta = apollo_beta,
# apollo_fixed = apollo_fixed,
# database = database,
# apollo_control = apollo_control
# )
#
# apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col)
#
# t_apollo <- system.time(
# fit_a <- apollo_estimate(
# apollo_beta, apollo_fixed, apollo_prob, apollo_inputs,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L)
# )
# )
#
# # ---- Extract estimates ----
# r_psi <- as.numeric(fit_r$stan_fit$par$psi)
# r_gamma <- as.numeric(fit_r$stan_fit$par$gamma)
# r_scale <- as.numeric(fit_r$stan_fit$par$scale)
#
# a_est <- fit_a$estimate
# a_gamma <- exp(a_est[grep("^ln_g", names(a_est))])
# a_scale <- exp(a_est["ln_sigma"])
#
# # Map column names: psi_j1 = alt-varying, psi_ind = individual-varying
# a_psi_j <- a_est[grep("^psi_j", names(a_est))]
# a_psi_ind <- a_est["psi_ind"]
# a_psi <- c(a_psi_j, a_psi_ind)
#
# list(
# ll = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout),
# psi = rbind(rmdcev = r_psi, apollo = a_psi),
# gamma = rbind(rmdcev = r_gamma[1:min(5, nalts)],
# apollo = a_gamma[1:min(5, nalts)]),
# scale = c(rmdcev = r_scale, apollo = as.numeric(a_scale)),
# time_s = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"])
# )
# }
## ----param-recovery-----------------------------------------------------------
# set.seed(42)
# true_psi_j <- 0.5
# true_psi_ind <- -1.5
# true_gamma <- rep(5, 5)
# true_scale <- 1.0
#
# sim_h0 <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = 500,
# nalts = 5,
# psi_j_parms = true_psi_j,
# psi_i_parms = true_psi_ind,
# gamma_parms = true_gamma,
# scale_parms = true_scale
# )
#
# cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1)
#
# cat("=== hybrid0 parameter comparison ===\n")
# cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
# cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"]))
# cat(sprintf(" LL difference = %.6f\n",
# abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"])))
#
# cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n")
# print(round(cmp_h0$psi, 4))
#
# cat("\nGamma (true: all 5)\n")
# print(round(cmp_h0$gamma, 4))
#
# cat(sprintf("\nScale (true: 1.0): rmdcev=%.4f apollo=%.4f\n",
# cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"]))
## ----param-recovery-gamma-----------------------------------------------------
# set.seed(42)
# true_alpha <- 0.5
#
# sim_g <- GenerateMDCEVData(
# model = "gamma",
# nobs = 500,
# nalts = 5,
# psi_j_parms = true_psi_j,
# psi_i_parms = true_psi_ind,
# gamma_parms = true_gamma,
# alpha_parms = true_alpha,
# scale_parms = true_scale
# )
#
# cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1)
#
# cat("=== gamma parameter comparison ===\n")
# cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
# cmp_g$ll["rmdcev"], cmp_g$ll["apollo"]))
# cat(sprintf(" LL difference = %.6f\n",
# abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"])))
#
# cat("\nPsi:\n"); print(round(cmp_g$psi, 4))
# cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4))
# cat(sprintf("\nScale: rmdcev=%.4f apollo=%.4f\n",
# cmp_g$scale["rmdcev"], cmp_g$scale["apollo"]))
## ----timing-setup-------------------------------------------------------------
# grid <- expand.grid(
# n = c(500, 1000, 5000),
# J = c(50, 100),
# rep = 1:3,
# stringsAsFactors = FALSE
# ) %>%
# filter(!(n == 1000 & J == 100)) # keep n=1000 for J=50 only
## ----timing-mle---------------------------------------------------------------
# timing_mle <- lapply(seq_len(nrow(grid)), function(i) {
# n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i]
# set.seed(r * 100 + i)
#
# sim <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = n,
# nalts = J,
# psi_j_parms = 0.5,
# psi_i_parms = -1.5,
# gamma_parms = rep(5, J),
# scale_parms = 1.0
# )
#
# res <- fit_both(sim, model = "hybrid0", npsi_j = 1)
# data.frame(n=n, J=J, rep=r,
# rmdcev_s = res$time_s["rmdcev"],
# apollo_s = res$time_s["apollo"],
# ll_diff = abs(res$ll["rmdcev"] - res$ll["apollo"]))
# })
#
# timing_mle <- do.call(rbind, timing_mle)
## ----timing-mle-summary-------------------------------------------------------
# summary_mle <- timing_mle %>%
# group_by(n, J) %>%
# summarise(
# rmdcev_mean_s = round(mean(rmdcev_s), 1),
# apollo_mean_s = round(mean(apollo_s), 1),
# speedup = round(mean(apollo_s) / mean(rmdcev_s), 1),
# ll_diff_mean = signif(mean(ll_diff), 3),
# .groups = "drop"
# )
#
# cat("=== MLE estimation timing (seconds, 3 replicates) ===\n")
# print(summary_mle, n = Inf)
## ----sim-apollo-prob5---------------------------------------------------------
# # Concrete apollo_probabilities for J=5 (used only in this section)
# apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
#
# V <- list(
# outside = 0,
# alt1 = psi_ind * b2 + psi_j1 * b1_1,
# alt2 = psi_ind * b2 + psi_j1 * b1_2,
# alt3 = psi_ind * b2 + psi_j1 * b1_3,
# alt4 = psi_ind * b2 + psi_j1 * b1_4,
# alt5 = psi_ind * b2 + psi_j1 * b1_5
# )
# mdcev_settings <- list(
# alternatives = c("outside","alt1","alt2","alt3","alt4","alt5"),
# budget = income, V = V,
# alpha = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0),
# gamma = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3),
# alt4=exp(ln_g4), alt5=exp(ln_g5)),
# sigma = exp(ln_sigma),
# continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2,
# alt3=quant_3, alt4=quant_4, alt5=quant_5),
# cost = list(outside=1, alt1=price_1, alt2=price_2,
# alt3=price_3, alt4=price_4, alt5=price_5),
# avail = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1)
# )
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }
## ----sim-convergence----------------------------------------------------------
# set.seed(42)
# J5 <- 5
# n500 <- 500
#
# sim5 <- GenerateMDCEVData(
# model = "hybrid0",
# nobs = n500,
# nalts = J5,
# psi_j_parms = 0.5,
# psi_i_parms = -1.5,
# gamma_parms = rep(5, J5),
# scale_parms = 1.0
# )
#
# # ── Fit both packages ──────────────────────────────────────────────────────
# fit_r5 <- mdcev(
# formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0,
# model = "hybrid0", algorithm = "MLE",
# print_iterations = FALSE, backend = "rstan"
# )
#
# database5 <- to_wide(sim5$data, J5)
# ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id",
# outputDirectory = tempdir(), nCores = 1)
# ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1)
# ai5 <- suppressMessages(
# apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0),
# database = database5, apollo_control = ac5))
# fit_a5 <- suppressMessages(suppressWarnings(
# apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L))
# ))
#
# # ── 1% price-change policy ─────────────────────────────────────────────────
# baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500)
# price_change <- 0.01 * colMeans(baseline_prices) # 1% of mean price per alt
#
# # rmdcev: additive policy vector (0 for numeraire, then per-alt changes)
# policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE)
# policies5$price_p[[2]] <- c(0, price_change)
#
# df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1)
#
# t_rmdcev5 <- system.time(
# demand_r5 <- mdcev.sim(
# df_sim5$df_indiv, df_common = df_sim5$df_common,
# sim_options = df_sim5$sim_options,
# cond_err = 0, nerrs = 100, sim_type = "demand"
# )
# )
#
# # Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index]
# mean_demand <- function(demand_list, policy_idx, alt_idx) {
# mean(sapply(demand_list, function(ind)
# mean(sapply(ind, function(s) s[policy_idx, alt_idx]))))
# }
# d_base_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a))
# d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a))
# delta_r <- d_policy_r - d_base_r
#
# # Apollo: predict at baseline and policy databases
# # Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt)
# database5_pol <- database5
# for (j in 1:J5) {
# database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j]
# }
# # Recompute quant_num at observed quantities under new prices
# pj_xj_pol <- Reduce("+", lapply(1:J5, function(j)
# database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]]))
# database5_pol$quant_num <- database5_pol$income - pj_xj_pol
#
# ai5_pol <- suppressMessages(
# apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0),
# database = database5_pol, apollo_control = ac5))
#
# t_apollo5_base <- system.time(
# pred_base <- suppressMessages(suppressWarnings(
# apollo_prediction(fit_a5, apollo_probabilities_j5, ai5,
# prediction_settings = list(nRep = 100, silent = TRUE))
# ))
# )
# t_apollo5_pol <- system.time(
# pred_pol <- suppressMessages(suppressWarnings(
# apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol,
# prediction_settings = list(nRep = 100, silent = TRUE))
# ))
# )
# t_apollo5 <- t_apollo5_base + t_apollo5_pol
#
# delta_a <- sapply(1:J5, function(j) {
# mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) -
# mean(pred_base[[paste0("alt", j, "_cont_mean")]])
# })
#
# cat("=== Demand change: 1% price increase (n=500, J=5) ===\n")
# cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff"))
# for (j in 1:J5) {
# cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j],
# abs(delta_r[j] - delta_a[j])))
# }
# cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a))))
# cat(sprintf("rmdcev %% change: %s%%\n",
# paste(round(100 * delta_r / d_base_r, 2), collapse = ", ")))
#
# cat(sprintf("\nTiming: rmdcev %.2f s | Apollo (2 predictions) %.2f s\n",
# t_rmdcev5["elapsed"], t_apollo5["elapsed"]))
## ----make-concrete-prob-------------------------------------------------------
# # Dynamically generate a concrete apollo_probabilities function for any J.
# # apollo_prediction() does source-code inspection and cannot use factory closures,
# # so we build and eval() a function body with the right number of alternatives.
# make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") {
# all_str <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")")
# V_items <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col,
# " + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n ")
# alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ")
# gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ")
# choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ")
# cost_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ")
# avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ")
#
# eval(parse(text = sprintf(
# 'function(apollo_beta, apollo_inputs, functionality = "estimate") {
# apollo_attach(apollo_beta, apollo_inputs)
# on.exit(apollo_detach(apollo_beta, apollo_inputs))
# P <- list()
# V <- list(outside = 0, %s)
# mdcev_settings <- list(
# alternatives = %s, budget = income, V = V,
# alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma),
# continuousChoice = list(%s), cost = list(%s), avail = list(%s)
# )
# P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
# P <- apollo_prepareProb(P, apollo_inputs, functionality)
# return(P)
# }', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items)))
# }
## ----timing-sim-grid----------------------------------------------------------
# timing_sim <- lapply(c(50, 100), function(J) {
# n <- 500
# set.seed(J)
#
# sim <- GenerateMDCEVData(
# model = "hybrid0", nobs = n, nalts = J,
# psi_j_parms = 0.5, psi_i_parms = -1.5,
# gamma_parms = rep(5, J), scale_parms = 1.0
# )
#
# # ── rmdcev ──────────────────────────────────────────────────────────────
# fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0",
# algorithm = "MLE", print_iterations = FALSE, backend = "rstan")
# bp <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n)
# pc <- 0.01 * colMeans(bp)
# pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE)
# pols$price_p[[2]] <- c(0, pc)
# df_s <- PrepareSimulationData(fit_r, pols, nsims = 1)
#
# rmdcev_times <- sapply(1:3, function(r)
# system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common,
# sim_options = df_s$sim_options,
# cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"])
#
# # ── Apollo ──────────────────────────────────────────────────────────────
# database <- to_wide(sim$data, J)
# ac <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id",
# outputDirectory = tempdir(), nCores = 1)
# ab <- make_apollo_beta(J, npsi_j = 1)
# pfn <- make_apollo_prob_concrete(J)
# ai <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0),
# database = database, apollo_control = ac)
# fit_a <- apollo_estimate(ab, character(0), pfn, ai,
# estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
# printLevel = 0L))
#
# # Policy database: 1% additive price increase
# db_pol <- database
# for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j]
# pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j)
# db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]]))
# db_pol$quant_num <- db_pol$income - pj_xj_pol
# ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0),
# database = db_pol, apollo_control = ac)
#
# apollo_times <- sapply(1:3, function(r) {
# t1 <- system.time(apollo_prediction(fit_a, pfn, ai,
# prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
# t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol,
# prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
# t1 + t2
# })
#
# data.frame(n = n, J = J,
# rmdcev_mean = round(mean(rmdcev_times), 2),
# rmdcev_sd = round(sd(rmdcev_times), 2),
# apollo_mean = round(mean(apollo_times), 2),
# apollo_sd = round(sd(apollo_times), 2),
# speedup = round(mean(apollo_times) / mean(rmdcev_times), 1))
# })
#
# timing_sim <- do.call(rbind, timing_sim)
# print(timing_sim, row.names = FALSE)
## ----session------------------------------------------------------------------
# sessionInfo()
Try the rmdcev package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.