beezdemand (version 0.2.0)

fit_demand_hurdle: 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"),
  ...
)

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

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).

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).$$

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.

See Also

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

Examples

# \donttest{
# 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)
# }