fit_demand_hurdle: Fit Two-Part Mixed Effects Hurdle Demand Model

In beezdemand: Behavioral Economic Easy Demand

Fit Two-Part Mixed Effects Hurdle Demand Model

Description

Fits a two-part hurdle model for demand data using TMB (Template Model Builder). Part I models the probability of zero consumption using logistic regression. Part II models log-consumption given positive response using a nonlinear mixed effects model.

Usage

fit_demand_hurdle(
  data,
  y_var,
  x_var,
  id_var,
  random_effects = c("zeros", "q0", "alpha"),
  epsilon = 0.001,
  start_values = NULL,
  tmb_control = list(max_iter = 200, eval_max = 1000, trace = 0),
  verbose = 1,
  part2 = c("zhao_exponential", "exponential", "simplified_exponential"),
  ...
)

Arguments

data	A data frame containing the demand data.
y_var	Character string specifying the column name for consumption values.
x_var	Character string specifying the column name for price.
id_var	Character string specifying the column name for subject IDs.
random_effects	Character vector specifying which random effects to include. Options are "zeros" (a_i for Part I), "q0" (b_i for intensity), and "alpha" (c_i for elasticity). Default is c("zeros", "q0", "alpha") for the full 3-random-effect model. Use c("zeros", "q0") for the simplified 2-random-effect model (fixed alpha across subjects).
epsilon	Small constant added to price before log transformation in Part I. Used to handle zero prices: log(price + epsilon). Default is 0.001.
start_values	Optional named list of starting values for optimization. If NULL (default), sensible defaults are used.
tmb_control	List of control parameters for TMB optimization: max_iter Maximum number of optimization iterations (default 200) eval_max Maximum number of function evaluations (default 1000) trace Print optimization trace: 0 = none, 1 = some (default 0)
verbose	Integer controlling output verbosity: 0 = silent, 1 = progress messages, 2 = detailed optimization trace. Default is 1.
part2	Character string selecting the Part II mean function. Options are "zhao_exponential" (default; no Q0 normalization in the exponent), "exponential" (HS-standardized; Q0 inside the exponent), and "simplified_exponential" (SND/log-linear; no k parameter).
...	Additional arguments (reserved for future use).

Details

The model structure is:

Part I (Binary - probability of zero consumption):

logit(\pi_{ij}) = \beta_0 + \beta_1 \cdot \log(price + \epsilon) + a_i

Part II (Continuous - log consumption given positive):

With 3 random effects (random_effects = c("zeros", "q0", "alpha")):

\log(Q_{ij}) = (\log Q_0 + b_i) + k \cdot (\exp(-\alpha_i \cdot price) - 1) + \epsilon_{ij}

where \alpha_i = \exp(\log(\alpha) + c_i) and k = \exp(\log(k)).

With 2 random effects (random_effects = c("zeros", "q0")):

\log(Q_{ij}) = (\log Q_0 + b_i) + k \cdot (\exp(-\alpha \cdot price) - 1) + \epsilon_{ij}

where \alpha = \exp(\log(\alpha)) and k = \exp(\log(k)).

Random effects follow a multivariate normal distribution with unstructured covariance matrix. Use compare_hurdle_models for likelihood ratio tests comparing nested models.

Value

An object of class beezdemand_hurdle containing:

model

List with coefficients, se, variance_components, correlations

random_effects

Matrix of empirical Bayes random effect estimates

subject_pars

Data frame of subject-specific parameters including Q0, alpha, breakpoint, Pmax, Omax

tmb_obj

TMB objective function object

opt

Optimization result from nlminb

sdr

TMB sdreport object

call

The matched call

data

Original data used for fitting

param_info

List with y_var, x_var, id_var, n_subjects, n_obs, etc.

converged

Logical indicating convergence

loglik

Log-likelihood at convergence

AIC, BIC

Information criteria

error_message

Error message if fitting failed, NULL otherwise

Parameterization and comparability

The TMB backend estimates positive-constrained parameters on the natural-log scale: \log(Q_0), \log(\alpha), and \log(k). Reporting methods (summary(), tidy(), coef()) can back-transform to the natural scale or present parameters on the \log_{10} scale.

To compare \alpha estimates with models fit in \log_{10} space, use:

\log_{10}(\alpha) = \log(\alpha) / \log(10).

See Also

summary.beezdemand_hurdle, predict.beezdemand_hurdle, plot.beezdemand_hurdle, compare_hurdle_models, simulate_hurdle_data

Examples

# Load example data
data(apt)

# Fit full model with 3 random effects
fit3 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0", "alpha"))

# Fit simplified model with 2 random effects (fixed alpha)
fit2 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0"))

# View results
summary(fit3)

# Compare models with likelihood ratio test
compare_hurdle_models(fit3, fit2)

beezdemand documentation built on March 3, 2026, 9:07 a.m.