hurdle-demand-models.R
In beezdemand: Behavioral Economic Easy Demand

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    dev = "png",
    dpi = 144,
    fig.width = 7,
    fig.height = 5,
    warning = FALSE,
    message = FALSE
)

library(beezdemand)

## ----data-example, echo=TRUE--------------------------------------------------
library(beezdemand)

# Example data structure
knitr::kable(
    head(apt, 10),
    caption = "Example APT data structure (first 10 rows)"
)

## ----basic-fit----------------------------------------------------------------
# Fit 2-RE model (simpler, faster)
fit2 <- fit_demand_hurdle(
    data = apt,
    y_var = "y",
    x_var = "x",
    id_var = "id",
    random_effects = c("zeros", "q0"), # 2 random effects
    verbose = 0
)

## ----fit-3re, eval=FALSE------------------------------------------------------
# # Fit 3-RE model (more flexible)
# fit3 <- fit_demand_hurdle(
#     data = apt,
#     y_var = "y",
#     x_var = "x",
#     id_var = "id",
#     random_effects = c("zeros", "q0", "alpha"), # 3 random effects
#     verbose = 1
# )

## ----summary-example----------------------------------------------------------
# View summary
summary(fit2)

# Extract coefficients
coef(fit2)

# Standardized tidy summaries
tidy(fit2) |> head()
glance(fit2)

# Get subject-specific parameters
head(get_subject_pars(fit2))

## ----diagnostics--------------------------------------------------------------
check_demand_model(fit2)
plot_residuals(fit2)$fitted
plot_qq(fit2)

## ----model-comparison, eval=FALSE---------------------------------------------
# # Compare nested models
# compare_hurdle_models(fit3, fit2)
# 
# # Unified model comparison (AIC/BIC + LRT when appropriate)
# compare_models(fit3, fit2)
# 
# # Output:
# # Likelihood Ratio Test
# # =====================
# #            Model n_RE    LogLik df     AIC     BIC
# #    Full (3 RE)    3 -1234.56 12 2493.12 2543.21
# # Reduced (2 RE)    2 -1245.78  9 2509.56 2547.89
# #
# # LR statistic: 22.44
# # df: 3
# # p-value: 5.2e-05

## ----plotting-demand, fig.cap="Population demand curve from 2-RE hurdle model."----
# Population demand curve
plot(fit2, type = "demand")

## ----plotting-probability, fig.cap="Probability of zero consumption as a function of price."----
# Probability of zero consumption
plot(fit2, type = "probability")

## ----plotting-extra, eval=FALSE-----------------------------------------------
# # Distribution of individual parameters
# plot(fit2, type = "parameters")
# plot(fit2, type = "parameters", parameters = c("Q0", "alpha", "Pmax"))
# 
# # Individual demand curves
# plot(fit2, type = "individual", subjects = c("1", "2", "3", "4"))
# 
# # Single subject with observed data
# plot_subject(fit2, subject_id = "1", show_data = TRUE, show_population = TRUE)

## ----simulate-basic, eval=FALSE-----------------------------------------------
# # Simulate with default parameters
# sim_data <- simulate_hurdle_data(
#     n_subjects = 100,
#     seed = 123
# )
# 
# head(sim_data)
# #   id    x         y delta       a_i        b_i
# # 1  1 0.00 12.345678     0 -0.523456  0.1234567
# # 2  1 0.50 11.234567     0 -0.523456  0.1234567
# # ...
# 
# # Custom parameters
# sim_custom <- simulate_hurdle_data(
#     n_subjects = 100,
#     logQ0 = log(15), # Q0 = 15
#     alpha = 0.1, # Lower elasticity
#     k = 3, # Higher k (ensures Pmax exists)
#     stop_at_zero = FALSE, # All prices for all subjects
#     seed = 456
# )

## ----monte-carlo, eval=FALSE--------------------------------------------------
# # Run Monte Carlo study
# mc_results <- run_hurdle_monte_carlo(
#     n_sim = 100, # Number of simulations
#     n_subjects = 100, # Subjects per simulation
#     n_random_effects = 2, # 2-RE model
#     verbose = TRUE,
#     seed = 123
# )
# 
# # View summary
# print_mc_summary(mc_results)
# 
# # Monte Carlo Simulation Summary
# # ==============================
# #
# # Simulations: 100 attempted, 95 converged (95.0%)
# #
# #    Parameter True Mean_Est   Bias Rel_Bias% Emp_SE Mean_SE SE_Ratio Coverage_95%  N
# #        beta0 -2.0    -2.01  -0.01      0.5   0.12    0.11     0.92         94.7 95
# #        beta1  1.0     1.02   0.02      2.0   0.08    0.08     1.00         95.8 95
# #        logQ0  2.3     2.29  -0.01     -0.4   0.05    0.05     1.00         94.7 95
# #            k  2.0     2.03   0.03      1.5   0.15    0.14     0.93         93.7 95
# #        alpha  0.5     0.51   0.01      2.0   0.04    0.04     1.00         95.8 95
# # ...

## ----integration, eval=FALSE--------------------------------------------------
# # Fit hurdle model
# hurdle_fit <- fit_demand_hurdle(data,
#     y_var = "y", x_var = "x", id_var = "id",
#     random_effects = c("zeros", "q0"), verbose = 0
# )
# 
# # Extract subject parameters
# hurdle_pars <- get_subject_pars(hurdle_fit)
# 
# # These can be merged with other analyses
# # e.g., correlate with individual characteristics
# cor(hurdle_pars$Q0, subject_characteristics$age)
# cor(hurdle_pars$breakpoint, subject_characteristics$dependence_score)

## ----export, eval=FALSE-------------------------------------------------------
# # Subject parameters
# write.csv(get_subject_pars(hurdle_fit), "hurdle_subject_parameters.csv")
# 
# # Summary statistics
# summ <- get_hurdle_param_summary(hurdle_fit)
# print(summ)

## ----control-params, eval=FALSE-----------------------------------------------
# fit <- fit_demand_hurdle(
#     data,
#     y_var = "y",
#     x_var = "x",
#     id_var = "id",
#     epsilon = 0.001, # Constant for log(price + epsilon)
#     tmb_control = list(
#         max_iter = 200, # Maximum iterations
#         eval_max = 1000, # Maximum function evaluations
#         trace = 0 # Optimization trace level
#     ),
#     verbose = 1 # 0 = silent, 1 = progress, 2 = detailed
# )

## ----custom-starts, eval=FALSE------------------------------------------------
# custom_starts <- list(
#     beta0 = -3.0,
#     beta1 = 1.5,
#     logQ0 = log(8),
#     k = 2.5,
#     alpha = 0.1,
#     logsigma_a = 0.5,
#     logsigma_b = -0.5,
#     logsigma_e = -1.0,
#     rho_ab_raw = 0
# )
# 
# fit <- fit_demand_hurdle(data,
#     y_var = "y", x_var = "x", id_var = "id",
#     start_values = custom_starts, verbose = 1
# )