Skip to contents

Predict Method for Hurdle Demand Models

Source: R/hurdle-methods.R
predict.beezdemand_hurdle.Rd

Returns predictions from a fitted hurdle demand model.

Usage

# S3 method for class 'beezdemand_hurdle'
predict(
  object,
  newdata = NULL,
  type = c("demand", "response", "link", "parameters", "probability"),
  prices = NULL,
  marginal = FALSE,
  marginal_method = c("normal", "kde", "empirical"),
  correction = TRUE,
  seed = 42L,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

object

An object of class beezdemand_hurdle.

newdata

Optional data frame containing a price column matching the fitted object's x_var. If newdata includes the id column, subject-specific predictions are returned; otherwise population predictions are returned. If newdata is NULL, returns predictions for all subjects across a price grid.

type

One of:

"demand"

(default) Predicted expected consumption = (1 - P0) * response (the marginal expectation; see Scoring predictions)

"response"

Predicted consumption conditional on consuming (part II), \(E[Y \mid Y > 0]\)

"link"

Predicted log-consumption (linear predictor of part II)

"probability"

Predicted probability of zero consumption (part I)

"parameters"

Subject-specific parameters (no .fitted column)

The default changed from "response" to "demand" in beezdemand 0.3.0; omitting type emits a one-time-per-session message naming the change.

prices

Optional numeric vector of prices used only when newdata = NULL.

marginal

Logical; if TRUE, computes population-averaged (marginal) predictions by integrating over the random effects distribution. For type = "probability", uses KDE/Normal/Empirical integration of the binary component. For type = "response" and type = "demand", uses Monte Carlo integration over all random effects, producing the population-average demand curve (accounting for Jensen's inequality). Default is FALSE, which gives conditional (RE = 0) predictions representing a "typical" subject at the center of the RE distribution.

marginal_method

Character. Method for marginal integration; one of "normal" (default; integrate over the model-assumed N(0, sigma_a) distribution of the zero-component intercept, over the whole real line), "kde" (kernel density estimate of the shrunken BLUPs), or "empirical" (simple average over the BLUPs). "normal" is the model-consistent choice: it integrates over the same distribution the fitted likelihood integrates over. "kde" and "empirical" are descriptive summaries of the shrunken BLUPs, which understate the random-effect spread (see Details). The default was "kde" in the development versions before 0.3.0. Ignored when marginal = FALSE.

correction

Logical; if TRUE (default), applies the lognormal retransformation correction exp(sigma_e^2 / 2) when back-transforming from the log scale to the natural consumption scale. This produces the arithmetic mean (conditional on Q > 0). Set to FALSE to obtain the median (geometric mean), which is useful for individual-level "most likely" predictions. Only applies to type = "response" and type = "demand".

seed

Integer or NULL. Random seed for Monte Carlo marginal predictions (default 42L). Set to NULL to use current RNG state. The global RNG state is preserved and restored after the call.

se.fit

Logical; if TRUE, includes a .se.fit column (delta-method via sdreport when available).

interval

One of "none" (default) or "confidence".

level

Confidence level when interval = "confidence".

...

Unused.

Value

For type = "parameters", a tibble of subject-level parameters. For type = "probability" with marginal = TRUE, a tibble with columns for price, prob_zero, and .fitted (no subject column). Otherwise, a tibble containing the newdata columns plus .fitted and helper columns predicted_log_consumption, predicted_consumption, prob_zero, and expected_consumption. When requested, .se.fit and .lower/.upper are included. Marginal results carry a marginal_method attribute; Monte Carlo marginal results (type = "response" / "demand") also carry a logical re_cov_fallback attribute that is TRUE when the random-effect draws had to use an uncorrelated (diagonal) covariance because the estimated covariance was not positive definite.

Details

Retransformation correction

The hurdle model specifies Gaussian errors on log-consumption (Part II): log(Q) ~ N(mu, sigma_e^2). The conditional distribution of Q given Q > 0 is therefore lognormal. The arithmetic mean of a lognormal is exp(mu + sigma_e^2/2), not exp(mu). Using exp(mu) returns the median (geometric mean), which systematically underestimates the arithmetic mean by a factor of exp(sigma_e^2/2). This correction is applied by default when type = "response" or type = "demand". Set correction = FALSE to obtain the median instead.

This is a parametric correction assuming normality of log-scale residuals (Duan, 1983). Under the model's normality assumption, this is equivalent to Duan's nonparametric smearing estimator.

Marginal P(zero)

The conditional P(zero) curve (when marginal = FALSE) sets the random intercept to zero, which produces a near step-function that misrepresents the observed fraction of zero responses. The marginal curve integrates over the random effect distribution, answering "what fraction of the population has stopped buying at this price?"

The "kde" and "empirical" methods integrate over empirical Bayes estimates (BLUPs) of the random intercepts. They are descriptive rather than model-consistent: BLUPs are shrunk toward zero compared to the true random effects (more so for subjects with few observations), so these methods understate the RE spread, and the marginal curve they produce summarises the fitted subjects rather than the population-level quantity the model defines. The "normal" method integrates over the model-assumed N(0, sigma_a) distribution, which is the model-consistent marginal (the same one the fitted likelihood integrates over) and is the choice to use when the marginal curve is reported as an estimate. It can be wrong only in the way the model itself is wrong, that is, if the normality assumption fails. The default remains "kde" for continuity with earlier versions. Use plot_qq() to assess RE normality.

Conditional vs. marginal demand predictions

Population-level demand predictions (when no subject ID is provided) can be computed in two ways:

  • Conditional (default, marginal = FALSE): Sets all random effects to zero and evaluates the demand function at the fixed-effect (population) parameters. For nonlinear models, this corresponds to the conditional mode rather than the population-average mean.

  • Marginal (marginal = TRUE): Integrates the prediction over the estimated random-effects distribution via Monte Carlo sampling. This gives the population-average demand curve. Due to Jensen's inequality, this curve lies above the conditional curve when the demand function is convex in the random effects (which it is for exponential demand with log-normal Q0 and alpha).

The conditional prediction is appropriate for characterizing a "typical" subject. The marginal prediction is appropriate for predicting aggregate consumption in a population.

Scoring predictions

A hurdle model has two natural "predicted consumption" quantities and they answer different questions. type = "response" is the conditional positive mean \(E[Y \mid Y > 0]\): consumption given that any is purchased. type = "demand" is the marginal expectation \((1 - p_0) E[Y \mid Y > 0]\), which weights in the probability of buying nothing. Observed consumption includes zeros, so predictions scored against raw data (cross-validation error, calibration, model comparison on predictions) must use type = "demand"; scoring with "response" systematically overstates error wherever \(p_0\) is large (typically at high prices). This is why "demand" is the default as of beezdemand 0.3.0.

Examples

# \donttest{
data(apt)
fit <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
#> Sample size may be too small for reliable estimation.
#>   Subjects: 10, Parameters: 12, Recommended minimum: 60 subjects.
#>   Consider using more subjects or the simpler 2-RE model.
#> Fitting HurdleDemand3RE model...
#>   Part II: zhao_exponential
#>   Subjects: 10, Observations: 160
#>   Fixed parameters: 12, Random effects per subject: 3
#>   Optimizing...
#>   Converged in 81 iterations
#>   Computing standard errors...
#> Done. Log-likelihood: 32.81

# Get subject-specific parameters
pars <- predict(fit, type = "parameters")

# Predict demand at specific prices
demand <- predict(fit, type = "demand", prices = c(0, 0.5, 1, 2, 5))
# }