Package 'beezdemand'

Description

AIC for Hurdle Demand Model

Usage

## S3 method for class 'beezdemand_hurdle'
AIC(object, ..., k = 2)

Arguments

Value

A numeric AIC value.

Description

Creates annotation layer

Usage

annotation_logticks2(
  base = 10,
  sides = "bl",
  scaled = TRUE,
  short = unit(0.1, "cm"),
  mid = unit(0.2, "cm"),
  long = unit(0.3, "cm"),
  colour = "black",
  size = 0.5,
  linetype = 1,
  alpha = 1,
  data = data.frame(x = NA),
  color = NULL,
  ...
)

Arguments

Details

Inherit and extend layer for use in ggplot draw

Value

ggplot2 layer

Author(s)

Shawn Gilroy [email protected]

Description

Compare nested hurdle demand models using likelihood ratio tests.

Usage

## S3 method for class 'beezdemand_hurdle'
anova(object, ...)

Arguments

Details

All models must be fit to the same data. Models are ordered by degrees of freedom, and sequential likelihood ratio tests are performed.

Value

An object of class anova.beezdemand containing:

table

Data frame with model comparison statistics

lrt

Likelihood ratio test results

Examples

data(apt)
fit2 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0"))
fit3 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0", "alpha"))
anova(fit2, fit3)

Description

Compare nested NLME demand models using likelihood ratio tests.

Usage

## S3 method for class 'beezdemand_nlme'
anova(object, ...)

Arguments

Details

For NLME models, this method delegates to nlme::anova.lme() on the underlying model objects when possible.

Value

An object of class anova.beezdemand containing model comparison statistics.

Examples

data(ko)
fit1 <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                         id_var = "monkey", equation_form = "zben",
                         random_effects = Q0 ~ 1)
fit2 <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                         id_var = "monkey", equation_form = "zben",
                         random_effects = Q0 + alpha ~ 1)
anova(fit1, fit2)

Description

For a single fit, computes joint Wald-chi-square tests on grouped fixed-effect coefficients. For multiple fits passed via ..., performs sequential likelihood-ratio tests on nested models.

Usage

## S3 method for class 'beezdemand_tmb'
anova(
  object,
  ...,
  test = c("Wald", "LRT", "AIC"),
  terms = NULL,
  group_by = c("auto", "parameter", "term")
)

Arguments

Details

The Wald statistic for a coefficient block is $W = \beta_g' \Sigma_{gg}^{-1} \beta_g$, asymptotically $\chi^2$ on length(beta_g) df. An exactly rank-deficient (perfectly collinear) block has a singular $\Sigma_{gg}$ and triggers an explicit error. A near-collinear block is not detected: $\Sigma_{gg}$ stays invertible and $W$ becomes large and unstable, so such a value should be interpreted with caution. For multiple fits, the likelihood-ratio test screens for detectable non-nesting (equal or decreasing degrees of freedom, or a larger model with lower log-likelihood) but cannot prove nesting from log-likelihood and df alone. Pass genuinely nested models. Rows of the multiple-fit table are ordered by ascending degrees of freedom, and the Model column labels them Model1, Model2, ... in that order.

Value

For a single fit, a tibble with Group, Chisq, df, p.value. For multiple fits, a tibble with Model, df, AIC, Chisq, `Pr(>Chisq)`.

See Also

anova.beezdemand_nlme(), confint.beezdemand_tmb().

Examples

data(apt_full)
# 40 subjects per gender keep the example fast; use the full data in practice
ids <- unique(apt_full[c("id", "gender")])
ids <- ids[ids$gender %in% c("Male", "Female"), ]
keep <- unlist(lapply(split(ids$id, ids$gender), head, 40))
dat <- apt_full[apt_full$id %in% keep, ]
fit <- fit_demand_tmb(dat, equation = "exponential",
                      factors = "gender", verbose = 0)
anova(fit)
anova(fit, group_by = "parameter")

Description

A dataset containing alcohol purchase task data for a small number of participants

Usage

apt

Format

Long-form data.frame with columns: id, x, y. Participants were asked how many standard sized alcoholic beverages they would buy at various prices.

Description

A larger dataset containing alcohol purchase task data with demographic covariates. Suitable for testing hurdle models and mixed-effects models with covariates.

Usage

apt_full

Format

A data frame with 18,700 rows and 8 columns:

id

Unique participant identifier (1-1100)

gender

Participant gender (Male/Female)

age

Participant age in years

binges

Number of binge drinking episodes

totdrinks

Total number of drinks consumed

tothours

Total hours spent drinking

x

Price point for the purchase task

y

Number of drinks participant would purchase at price x

Examples

data(apt_full)
# Use a subset for quick demonstration
apt_sub <- apt_full[apt_full$id %in% unique(apt_full$id)[1:20], ]
fit <- fit_demand_hurdle(apt_sub, y_var = "y", x_var = "x", id_var = "id")
summary(fit)

Description

Returns the original data with fitted values and residuals from individual demand curve fits. This enables easy model diagnostics and visualization with the tidyverse.

Usage

## S3 method for class 'beezdemand_fixed'
augment(x, newdata = NULL, ...)

Arguments

Details

For "hs" equation models where fitting is done on the log10 scale, fitted values are back-transformed to the natural scale.

Value

A tibble containing the original data plus:

.fitted

Fitted demand values on the response scale

.resid

Residuals (observed - fitted)

Examples

data(apt)
fit <- fit_demand_fixed(apt, y_var = "y", x_var = "x", id_var = "id")
augmented <- augment(fit)

# Plot residuals by subject
library(ggplot2)
ggplot(augmented, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~id) +
  geom_hline(yintercept = 0, linetype = "dashed")

Description

Returns the original data with fitted values, residuals, and predictions from a hurdle demand model, in the format expected by tidyverse plotting and diagnostic tools.

Usage

## S3 method for class 'beezdemand_hurdle'
augment(x, newdata = NULL, component = c("combined", "continuous"), ...)

Arguments

Details

Residual types for hurdle models

The hurdle model has two components, each requiring different diagnostic approaches:

• Continuous residuals (component = "continuous"): Standard log-scale residuals log(y) - mu for observations where y > 0. Zeros are excluded (NA). Assesses whether the lognormal conditional distribution is well- specified.

• Randomized quantile residuals (component = "combined", default): Following Dunn & Smyth (1996), maps each observation through the fitted hurdle CDF and then the standard normal quantile function. If the model is correctly specified, these residuals are exactly N(0,1) regardless of which component generated the observation. For zeros, a uniform random variate within ⁠[0, P(zero)]⁠ breaks ties that would otherwise create a spike in the QQ plot.

Value

A tibble containing the original data plus:

.fitted

Fitted demand values (natural scale)

.fitted_link

Fitted values on log scale (Part II mean)

.fitted_prob

Predicted probability of consumption (1 - P(zero))

.resid

Residuals (type depends on component; see above)

.resid_response

Residuals on response scale (y - .fitted)

Examples

data(apt)
fit <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
augmented <- augment(fit)

# Plot residuals
library(ggplot2)
ggplot(augmented, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed")

Description

Returns the original data with fitted values and residuals from a nonlinear mixed-effects demand model. This enables easy model diagnostics and visualization with the tidyverse.

Usage

## S3 method for class 'beezdemand_nlme'
augment(x, newdata = NULL, ...)

Arguments

Details

The fitted values and residuals are on the same scale as the response variable used in the model. For equation_form = "zben", this is the LL4-transformed scale. For equation_form = "simplified" or "exponentiated", this is the natural consumption scale.

To back-transform predictions to the natural scale for "zben" models, use: ll4_inv(augmented$.fitted)

Value

A tibble containing the original data plus:

.fitted

Fitted values on the model scale (may be transformed, e.g., LL4)

.resid

Residuals on the model scale

.fixed

Fitted values from fixed effects only (population-level)

Examples

data(ko)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", factors = "dose", equation_form = "zben")
augmented <- augment(fit)

# Plot residuals
library(ggplot2)
ggplot(augmented, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed")

Description

Augment a beezdemand_tmb Model

Usage

## S3 method for class 'beezdemand_tmb'
augment(x, newdata = NULL, ...)

Arguments

Value

A tibble with original data plus .fitted, .resid, and .std_resid columns. Residuals are computed on the model's native scale (log scale for "exponential", natural/LL4 scale for others) to match the C++ likelihood.

Examples

data(apt)
fit <- fit_demand_tmb(apt, equation = "exponential", verbose = 0)
head(augment(fit))

Description

Augment a Cross-Price Demand Model (Linear)

Usage

## S3 method for class 'cp_model_lm'
augment(x, ...)

Arguments

Value

A tibble with the original modelling data and added .fitted and .resid columns.

Description

Augment a Cross-Price Demand Model (Mixed-Effects)

Usage

## S3 method for class 'cp_model_lmer'
augment(x, ...)

Arguments

Value

A tibble with the original modelling data and added .fitted, .resid, and .fixed (population-level prediction with random effects set to zero) columns.

Description

Augment a Cross-Price Demand Model (Nonlinear)

Usage

## S3 method for class 'cp_model_nls'
augment(x, ...)

Arguments

Value

A tibble with the original modelling data and added .fitted and .resid columns. Both are on the model's response scale: raw y for equation = "exponentiated" / "additive", and log10(y) for equation = "exponential" (which is fit to log10(y)).

Description

Methods for printing, summarizing, and visualizing objects of class beezdemand_descriptive created by get_descriptive_summary().

Usage

## S3 method for class 'beezdemand_descriptive'
print(x, ...)

## S3 method for class 'beezdemand_descriptive'
summary(object, ...)

## S3 method for class 'beezdemand_descriptive'
plot(x, x_trans = "identity", y_trans = "identity", show_zeros = FALSE, ...)

Arguments

Details

Print Method

Displays a compact summary showing the number of subjects and prices analyzed, plus a preview of the statistics table.

Summary Method

Provides extended information including:

• Data summary (subjects, prices analyzed)

• Distribution of means across prices (min, median, max)

• Proportion of zeros by price (range)

• Missing data summary

Plot Method

Creates a boxplot showing the distribution of consumption at each price point. Features:

• Red cross markers indicate means

• Boxes show median and quartiles

• Whiskers extend to 1.5 * IQR

• Supports axis transformations (log, sqrt, etc.)

• Uses modern beezdemand styling via theme_apa()

Value

• print() - Returns the object invisibly (called for side effects)

• summary() - Returns a list with extended summary information

• plot() - Returns a ggplot2 object

See Also

get_descriptive_summary()

Examples

data(apt, package = "beezdemand")
desc <- get_descriptive_summary(apt)

# Print compact summary
print(desc)

# Extended summary
summary(desc)

# Default boxplot
plot(desc)

# With log-transformed y-axis
plot(desc, y_trans = "log10")

# With pseudo-log y-axis (defined at zero)
plot(desc, y_trans = "pseudo_log")

Description

Methods for printing, summarizing, and visualizing objects of class beezdemand_empirical created by get_empirical_measures().

Usage

## S3 method for class 'beezdemand_empirical'
print(x, ...)

## S3 method for class 'beezdemand_empirical'
summary(object, ...)

## S3 method for class 'beezdemand_empirical'
plot(x, type = "histogram", ...)

Arguments

Details

Print Method

Displays a compact summary showing the number of subjects analyzed and a preview of the empirical measures table.

Summary Method

Provides extended information including:

• Data summary (subjects, zero consumption patterns, completeness)

• Descriptive statistics for each empirical measure (min, median, mean, max, SD)

• Missing data patterns

Plot Method

Creates visualizations of empirical measures across subjects.

Histogram type (default):

• Six-panel faceted plot showing distribution of each measure

• Helps identify central tendencies and outliers

Matrix type:

• Scatterplot matrix (pairs plot) showing relationships between measures

• Useful for identifying correlated demand metrics

• Lower triangle: scatterplots with smoothed trend lines

• Diagonal: density plots

• Upper triangle: correlation coefficients

Value

• print() - Returns the object invisibly (called for side effects)

• summary() - Returns a list with extended summary information

• plot() - Returns a ggplot2 object

See Also

get_empirical_measures()

Examples

data(apt, package = "beezdemand")
emp <- get_empirical_measures(apt)

# Print compact summary
print(emp)

# Extended summary
summary(emp)

# Histogram of measure distributions
plot(emp)

# Scatterplot matrix
plot(emp, type = "matrix")

Description

BIC for Hurdle Demand Model

Usage

## S3 method for class 'beezdemand_hurdle'
BIC(object, ...)

Arguments

Value

A numeric BIC value.

Description

Computes confidence intervals on derived demand metrics (Pmax, Omax, Qmax, EV, elasticity-at-Pmax) for a TMB mixed-effects demand fit, via a parametric bootstrap. Draws of the fixed-effect parameter vector are taken from the asymptotic Gaussian sampling approximation $N(\hat\beta, \hat\Sigma)$ (the distribution Wald intervals assume; no prior is involved), mapped to per-condition $(Q_0, \alpha, k)$ through the model's fixed-effect design, passed through the canonical Pmax/Omax engine, and summarized by empirical quantiles.

For a factor-expanded fit, one CI is returned per factor cell (e.g. one row per gender level); for an intercept-only fit, one row per statistic with condition = NA. The per-cell point estimate reproduces calc_group_metrics(fit, at = cell); the bootstrap supplies only the interval.

Usage

boot_demand(
  fit,
  statistics = c("Pmax", "Omax", "EV"),
  method = c("parametric"),
  R = 1000L,
  ci_level = 0.95,
  at = NULL,
  seed = NULL,
  ...
)

Arguments

Details

The parametric bootstrap is asymptotically equivalent to the delta method but avoids its linearization, so it is preferable for the strongly nonlinear derived metrics (Pmax/Omax via Lambert-W). Draws are fixed-effect-only (population / per-condition metrics); per-subject metric CIs would require random-effect-aware draws and are out of scope for now.

When k is estimated, its uncertainty is propagated (the log_k column is in the draw matrix); when k is fixed, the fixed value is used. The point estimate always uses the point k, matching calc_group_metrics().

Note that percentile intervals of a nonlinear transform are not guaranteed to bracket the point estimate; conf.low <= conf.high always holds, but ⁠conf.low <= estimate <= conf.high⁠ may not at boundary cases.

Some draws can leave a metric's domain (e.g. Pmax is undefined when a drawn k falls below the Lambert-W threshold). Such non-finite draws are excluded from the quantiles; the per-row count of excluded draws is recorded in attr(x, "n_nonfinite") (so the realized draw count is R minus that). If every draw of a requested metric/condition is non-finite, an error is raised because the interval is undefined.

Value

A tibble with one row per ⁠(statistic, condition)⁠:

statistic

Metric name.

condition

Factor-cell label (e.g. "gender=Male"); NA when the fit has no factors.

estimate

Point estimate (from the coefficient vector and point k).

conf.low, conf.high

Empirical-quantile interval bounds.

level

The confidence level used.

The bootstrap settings are attached as attributes "method", "R", and "seed", plus "n_nonfinite" (a per-row count of draws excluded as non-finite).

See Also

calc_group_metrics, fit_demand_tmb, confint.beezdemand_tmb

Examples

data(apt, package = "beezdemand")
fit <- fit_demand_tmb(apt, equation = "exponential", verbose = 0)
boot_demand(fit, statistics = c("Pmax", "Omax", "EV"), R = 500, seed = 1)

Description

Calculates group-level (population) Omax and Pmax from a fitted hurdle demand model.

Usage

calc_group_metrics(object, ...)

Arguments

Value

A named list with group-level Pmax, Omax, and Qmax.

See Also

calc_omax_pmax, fit_demand_hurdle

Examples

data(apt)
fit <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
calc_group_metrics(fit)

Description

Calculate Population-Level Demand Metrics for TMB Model

Usage

## S3 method for class 'beezdemand_tmb'
calc_group_metrics(object, at = NULL, ...)

Arguments

Value

A list with Pmax, Omax, Qmax, elasticity_at_pmax, method, and conditioned_on describing the reference point used. The conditioned_on field reports the actual conditioning applied (covariate values used, factor treatment per factor) so programmatic consumers do not have to re-derive it.

Marginalization order

For derived metrics (Pmax/Omax/Qmax) that depend nonlinearly on Q0 and alpha jointly, this function marginalizes parameters first then derives metrics:

1. Compute log-Q0 and log-alpha EMMs at each cell of the reference grid produced by .tmb_build_emm_ref_grid().

2. Marginalize each parameter across factor cells with equal weights (matches the emmeans default).

3. Derive Pmax/Omax/Qmax from the marginalized log-parameters at the user-supplied (or training-mean default) covariate point.

The result is "metrics evaluated at the average parameter values" rather than "average metrics across cells". The two answers differ for nonlinear transforms. The convention matches the parameter-level marginalization used by get_demand_param_emms().

Marginalisation policy (TMB vs NLME)

This method averages the log-scale linear predictors with equal weight over the factor cells observed in the fitting data and then exponentiates. The NLME method (calc_group_metrics.beezdemand_nlme()) follows the same policy: per-cell EMMs from emmeans, rows for unobserved cells dropped, equal-weight geometric mean over the remaining cells. Cells of the full factorial crossing with no subjects are never averaged over by either backend. The two agree for any design fit in log space and can differ only when the NLME fit uses param_space = "natural". In both backends continuous covariates are held at the training mean unless at supplies a single value; a multi-value continuous at entry warns and uses its first value.

See Also

fit_demand_tmb, get_demand_param_emms.beezdemand_tmb

Examples

data(apt)
fit <- fit_demand_tmb(apt, equation = "exponential", verbose = 0)
calc_group_metrics(fit)
# Conditioned at a specific covariate value:
# calc_group_metrics(fit_with_cov, at = list(FTND_z = 1))

Description

Calculate Observed Pmax/Omax Grouped by ID

Usage

calc_observed_pmax_omax(
  data,
  id_var = "id",
  price_var = "x",
  consumption_var = "y"
)

Arguments

Value

Data frame with observed pmax/omax for each subject

Examples

data(apt, package = "beezdemand")
calc_observed_pmax_omax(apt, id_var = "id", price_var = "x", consumption_var = "y")

Description

Calculates the maximum expenditure (Omax) and the price at maximum expenditure (Pmax) for the exponential demand model used in the two-part hurdle model.

Usage

calc_omax_pmax(Q0, k, alpha, price_range = NULL)

Arguments

Details

For the demand function:

$Q(p) = Q_0 \cdot \exp(k \cdot (\exp(-\alpha \cdot p) - 1))$

Expenditure is E(p) = p * Q(p). Omax is the maximum of E(p) and Pmax is the price at which this maximum occurs. These are found numerically.

The search range is automatically adjusted based on alpha to ensure the maximum is found. For small alpha values, Pmax can be quite large.

Value

A named list with:

Pmax

Price at maximum expenditure

Omax

Maximum expenditure (price * quantity)

Qmax

Quantity at Pmax

See Also

calc_group_metrics, fit_demand_hurdle

Examples

# Calculate for group-level parameters
calc_omax_pmax(Q0 = 10, k = 2, alpha = 0.5)

# With k > e (~2.718), a strict interior maximum exists
calc_omax_pmax(Q0 = 10, k = 3, alpha = 0.5)

Description

Calculates Amplitude and Persistence latent factors from demand metrics for various beezdemand model objects.

Usage

calculate_amplitude_persistence(
  fit,
  amplitude = c("Intensity", "Q0d", "Q0"),
  persistence = c("BP0", "Pmaxe", "Omaxe", "Alpha"),
  use_inv_alpha = TRUE,
  strict = TRUE,
  min_persistence_components = 2L,
  empirical_y_var = NULL,
  basis_means = NULL,
  basis_sds = NULL,
  ...
)

Arguments

Details

This function calculates Amplitude and Persistence by:

1. Extracting the relevant demand metrics from the fit object.

2. Resolving requested metric columns (case-insensitive, with limited synonym support for common beezdemand outputs).

3. Inverting Alpha if requested (1/Alpha).

4. Standardizing (Z-scoring) the metrics. If basis_means and basis_sds are provided, they are used; otherwise, the current sample's statistics are used.

5. Aggregating the Z-scores into the two latent factors.

Amplitude is defined by the variable specified in amplitude (typically Intensity/Q0). Persistence is defined as the mean of the standardized values of the variables specified in persistence (typically Breakpoint, Pmax, Omax, and 1/Alpha).

Value

A data frame with the original ID and calculated Amplitude and Persistence scores, along with the standardized (Z-scored) constituent metrics.

Models

• beezdemand_fixed: Extracts metrics from fit$results.

• beezdemand_hurdle: Extracts metrics from fit$subject_pars.

• beezdemand_nlme: Calculates subject-specific parameters from fixed and random effects. Parameters Q0 and Alpha are assumed to be on log10 scale for zben and simplified equations and are converted to linear scale. Omax and Pmax are calculated empirically from predictions. Breakpoint is calculated empirically from the raw data.

Examples

data(apt, package = "beezdemand")
fit <- FitCurves(apt, "hs", k = "share")
calculate_amplitude_persistence(fit)

Description

Cross-price style data with cannabis and cigarette context.

Usage

cannabisCigarettes

Format

A data frame with columns including: id, x, y, commodity, and auxiliary fields (Q1035, CigPrice, CanPrice, variable, value).

Examples

# data(cannabisCigarettes)

Description

Changes demand data

Usage

ChangeData(
  dat,
  nrepl = 1,
  replnum = 0.01,
  rem0 = FALSE,
  remq0e = FALSE,
  replfree = NULL,
  xcol = "x",
  ycol = "y",
  idcol = "id"
)

Arguments

Details

Change demand data in various ways. Ways include replacing any number of 0 values with a replacement number (or remove them completely), removing price and consumption at free, replacing free with some number. This will soon replace ReplaceZeros and certain arguments in FitCurves.

Value

Long form dataframe resembling the originally provided dataframe

Author(s)

Brent Kaplan [email protected]

Examples

## Change just the first instance of 0 within each unique value of id with .1
ChangeData(apt, nrepl = 1, replnum = .1)

Description

Performs diagnostic checks on fitted demand models, returning information about convergence, boundary conditions, and residual patterns.

Usage

check_demand_model(object, ...)

## S3 method for class 'beezdemand_hurdle'
check_demand_model(object, ...)

## S3 method for class 'beezdemand_nlme'
check_demand_model(object, ...)

## S3 method for class 'beezdemand_fixed'
check_demand_model(object, ...)

## S3 method for class 'beezdemand_tmb'
check_demand_model(object, ...)

Arguments

Details

The function checks for:

• Convergence status and optimization messages

• Parameters at or near boundaries

• Residual patterns (heteroscedasticity, outliers)

• Random effect variance estimates near zero

• For fit_demand_tmb() fits with estimate_k = TRUE, whether the fit shows the signature of an unidentified k: an implausible k, a decay exponent alpha * Q0 * price that never leaves its linear regime over the observed prices, a degenerate alpha, or log_k resting on a user-supplied bound. This is a screen, not a formal test: no flag means nothing was detected, not that k is identified.

• Correlation matrices near singularity

Value

An object of class beezdemand_diagnostics containing:

convergence

List with convergence status and messages

boundary

List with boundary condition warnings. For beezdemand_tmb fits it also carries k_identification: NULL unless k was estimated as a free parameter, otherwise the screen described in Details.

residuals

Summary statistics for residuals

random_effects

Summary of random effects (if applicable). variances holds random-effect variances: on the log10 scale for beezdemand_tmb fits, and as reported by nlme::VarCorr() for beezdemand_nlme fits. beezdemand_tmb fits also carry sd_log10 (the corresponding standard deviations, matching summary(fit)$variance_components) and sd_internal_log (the raw natural-log-scale standard deviations used by the near-zero degeneracy check).

issues

Character vector of identified issues

recommendations

Character vector of recommendations

Note

This function is named check_demand_model() to avoid potential conflicts with performance::check_model() from the performance package.

See Also

plot_residuals(), plot_qq()

Examples

data(apt)
fit <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
diagnostics <- check_demand_model(fit)
print(diagnostics)

Description

Modern interface for screening cross-price data with standardized output vocabulary aligned with check_systematic_demand().

Usage

check_systematic_cp(
  data,
  trend_threshold = 0.025,
  bounce_threshold_down = 0.1,
  bounce_threshold_up = 0.1,
  bounce_threshold_none = 0.1,
  consecutive_zeros = 2,
  consecutive_nonzeros = 2,
  expected_down = FALSE,
  x_var = "x",
  y_var = "y",
  id_var = "id",
  by = NULL
)

Arguments

Details

If the data contains an id column (or column specified by id_var), each unique ID is checked separately. Otherwise, the entire dataset is treated as a single pattern.

For cross-price data, the wrapper preserves the legacy meaning of check_unsystematic_cp():

• trend_direction and bounce_direction are taken directly from the legacy function outputs.

• trend_pass is set to NA because cross-price systematicity does not use a separate trend “pass/fail” criterion in the same way as purchase-task screening; instead, trend classification determines which bounce rule applies.

• bounce_stat is reported as the proportion relevant to the legacy bounce rule for the detected trend_direction (or expected_down case), computed from the legacy bounce counts and the number of price steps.

Value

An object of class beezdemand_systematicity with the same structure as check_systematic_demand(), with type = "cp".

Examples

data(etm)
check <- check_systematic_cp(etm)

Description

Modern interface for screening purchase task data using Stein et al. (2015) criteria. Returns a structured object with standardized output vocabulary that is consistent with check_systematic_cp().

Usage

check_systematic_demand(
  data,
  trend_threshold = 0.025,
  bounce_threshold = 0.1,
  max_reversals = 0,
  consecutive_zeros = 2,
  x_var = "x",
  y_var = "y",
  id_var = "id",
  by = NULL
)

Arguments

Details

The results tibble contains standardized columns for both demand and cross-price systematicity checks:

id

Subject identifier

type

"demand" for this function

trend_stat

DeltaQ statistic (log-log slope)

trend_threshold

Threshold used

trend_direction

"down", "up", or "none"

trend_pass

Logical: passed trend criterion

bounce_stat

Bounce proportion

bounce_threshold

Threshold used

bounce_direction

"significant" or "none"

bounce_pass

Logical: passed bounce criterion

reversals

Count of reversals from zero

reversals_pass

Logical: passed reversals criterion

returns

NA for demand (CP-specific)

n_positive

Count of positive values

systematic

Logical: passed all criteria

Value

An object of class beezdemand_systematicity with components:

results

Tibble with one row per subject containing systematicity metrics

type

"demand"

call

The original function call

n_total

Total number of subjects

n_systematic

Number of subjects passing all criteria

n_unsystematic

Number of subjects failing at least one criterion

Examples

data(apt)
check <- check_systematic_demand(apt)
print(check)
summary(check)
tidy(check)

# Grouped check: results include group column
data(apt_full)
check_g <- check_systematic_demand(apt_full, by = "gender")
check_g$results

Description

Analyzes whether consumption data shows systematic trends or unsystematic patterns ("bounces") with respect to price. Includes detection of zero-value reversal/return sequences and allows flexible output based on the level of detail requested. See Rzeszutek et al. (in press) for more details.

Usage

check_unsystematic_cp(
  data,
  delta_threshold = 0.025,
  bounce_down_threshold = 0.1,
  bounce_up_threshold = 0.1,
  bounce_none_threshold = 0.1,
  rev_zeroes = 2,
  ret_nums = 2,
  expected_down = FALSE,
  verbose = FALSE,
  detailed = FALSE
)

Arguments

Value

A data frame of class cp_unsystematic with core results:

delta_direction

Character: 'down', 'up', or 'none'.

bounce_direction

Character: 'up', 'down', 'significant', or 'none'.

bounce_any

Logical. TRUE if any bounce pattern detected.

bounce_above

Integer. Number of upward changes meeting threshold.

bounce_below

Integer. Number of downward changes meeting threshold.

reversals

Integer. Detected reversals from 0 to non-0.

returns

Integer. Detected returns from non-0 to 0.

If detailed = TRUE, returns additional columns:

delta_down

Logical. Significant downward trend.

delta_up

Logical. Significant upward trend.

delta_none

Logical. No significant trend.

bounce_up

Logical. Significant bounce up in a downward trend.

bounce_down

Logical. Significant bounce down in an upward trend.

bounce_none

Logical. Significant bounces in no-trend data.

Examples

x_seq <- 10^(seq(-2, 2, length.out = 10))
pattern <- data.frame(x = x_seq, y = c(10, 5, 10, 9, 10, 13, 10, 10, 7, 9))
check_unsystematic_cp(pattern)

Description

Checks to ensure column names are specified

Usage

CheckCols(dat, xcol, ycol, idcol, groupcol = NULL)

Arguments

Details

Check column names

Value

Dataframe

Author(s)

Brent Kaplan [email protected]

Examples

dat <- data.frame(price = 1:5, quantity = c(10, 8, 5, 2, 0), subj = rep(1, 5))
CheckCols(dat, xcol = "price", ycol = "quantity", idcol = "subj")

Description

Applies Stein, Koffarnus, Snider, Quisenberry, & Bickel's (2015) criteria for identification of nonsystematic purchase task data.

Usage

CheckUnsystematic(dat, deltaq = 0.025, bounce = 0.1, reversals = 0, ncons0 = 2)

Arguments

Details

This function applies the 3 criteria proposed by Stein et al., (2015) for identification of nonsystematic purchase task data. The three criteria include trend (deltaq), bounce, and reversals from 0. Also reports number of positive consumption values.

Value

Dataframe

Author(s)

Brent Kaplan [email protected]

Examples

## Using all default values
CheckUnsystematic(apt, deltaq = 0.025, bounce = 0.10, reversals = 0, ncons0 = 2)
## Specifying just 1 zero to flag as reversal
CheckUnsystematic(apt, deltaq = 0.025, bounce = 0.10, reversals = 0, ncons0 = 1)

Description

Methods to extract coefficients from various cross-price demand model objects.

Usage

## S3 method for class 'cp_model_nls'
coef(object, ...)

## S3 method for class 'cp_model_lm'
coef(object, ...)

## S3 method for class 'cp_model_lmer'
coef(object, fixed_only = FALSE, combine = TRUE, ...)

Arguments

Value

Named vector of coefficients

A named numeric vector of model coefficients.

Functions

• coef(cp_model_nls): Extract coefficients from a nonlinear cross-price model

• coef(cp_model_lm): Extract coefficients from a linear cross-price model

• coef(cp_model_lmer): Extract coefficients from a mixed-effects cross-price model

Description

Extract Coefficients from Fixed-Effect Demand Fit

Usage

## S3 method for class 'beezdemand_fixed'
coef(object, report_space = c("internal", "natural", "log10"), ...)

Arguments

Value

A tibble with columns id, term, estimate, estimate_scale, term_display.

Description

Extract Coefficients from Grouped Fixed-Effect Demand Fit

Usage

## S3 method for class 'beezdemand_fixed_grouped'
coef(object, ...)

Arguments

Value

A tibble with group columns prepended

Description

Extract Coefficients from Hurdle Demand Model

Usage

## S3 method for class 'beezdemand_hurdle'
coef(object, report_space = c("natural", "log10", "internal"), ...)

Arguments

Value

Named numeric vector of fixed effect coefficients.

Description

Provides methods to extract fixed effects, random effects, or subject-specific (combined fixed + random) coefficients from a beezdemand_nlme object. This is an S3 method for the generic coef function.

Usage

## S3 method for class 'beezdemand_nlme'
coef(
  object,
  type = "combined",
  report_space = c("internal", "natural", "log10"),
  ...
)

Arguments

Value

Depending on type:

• type="fixed": A named numeric vector of fixed-effect coefficients.

• type="random": A data frame (or list of data frames if multiple levels of grouping) of random effects, as returned by ranef.nlme().

• type="combined": A data frame where rows are subjects (from id_var) and columns are the Q0 and alpha parameters, representing subject-specific estimates (on the log10 scale).

See Also

fixef.beezdemand_nlme, ranef.beezdemand_nlme

Examples

data(ko)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", equation_form = "zben")
coef(fit, type = "fixed")
coef(fit, type = "random")
coef(fit, type = "combined")

Description

Extract coefficients from a fitted beezdemand_tmb model. The type argument selects the return shape. The default, "internal", is unchanged: a named numeric vector of the optimizer's flat parameterization (entries include beta_q0, beta_alpha, logsigma*, and any covariance hyperparameters; intercepts are on the log scale because the optimizer works in unconstrained space). This is the numeric-vector escape hatch consumed by tooling such as car::deltaMethod and multcomp::glht.

Usage

## S3 method for class 'beezdemand_tmb'
coef(object, type = c("internal", "subject", "combined", "fixed"), ...)

Arguments

Details

type = "subject" (alias "combined") returns the per-subject parameter tibble from get_subject_pars (with expanded = NULL, so within-id factor expansion is auto-detected). This is concept-parity with coef.beezdemand_nlme(type = "combined") but not column-identical: it returns resolved per-subject parameters (Q0, alpha, ...), not a per-design-term coefficient matrix. type = "fixed" returns a one-row tibble of the fixed-effect coefficients only (the beta_q0 / beta_alpha block on the internal parameterization), excluding log_k, logsigma*, and rho*.

Scale conversion is not performed here: supplying report_space through ... is an error. Use get_subject_pars or predict.beezdemand_tmb for natural-scale parameters.

Value

For type = "internal", a named numeric vector. For type = "subject"/"combined", a tibble with one row per subject (or one row per subject-by-within-id-factor-level cell when the fit has within-id factor variation). For type = "fixed", a one-row tibble of fixed-effect coefficients.

Examples

data(apt)
fit <- fit_demand_tmb(apt, equation = "exponential", verbose = 0)
coef(fit)                    # raw optimizer vector (default, "internal")
coef(fit, type = "subject")  # per-subject parameter tibble
coef(fit, type = "fixed")    # fixed-effect coefficients

Description

Performs a likelihood ratio test comparing two nested hurdle demand models. Typically used to test whether adding the random effect on alpha (c_i) significantly improves model fit (3-RE vs 2-RE models).

Usage

compare_hurdle_models(model_full, model_reduced)

Arguments

Value

Invisibly returns a list with:

lr_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

P-value from chi-squared distribution

model_comparison

Data frame with model comparison statistics

See Also

fit_demand_hurdle

Examples

data(apt)
fit3 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0", "alpha"))
fit2 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0"))
compare_hurdle_models(fit3, fit2)

Description

Compare multiple demand models using information criteria and likelihood ratio tests (when applicable). Works with all beezdemand model classes.

Usage

compare_models(..., test = c("auto", "lrt", "none"))

Arguments

Details

Models are compared using AIC and BIC. For models from the same statistical backend (e.g., two hurdle models or two NLME models), likelihood ratio tests can be performed if the models are nested.

When comparing models from different backends (e.g., hurdle vs NLME), only information criteria comparisons are possible since the likelihoods are not directly comparable for LRT purposes.

Backend Compatibility

Value

An object of class beezdemand_model_comparison containing:

comparison

Data frame with model fit statistics

test_type

Type of test performed

lrt_results

LRT results if performed (NULL otherwise)

best_model

Index of best model by BIC

notes

Character vector of notes/warnings

nesting_verified

Logical; always FALSE since nesting is not automatically verified. Users must ensure models are properly nested for valid LRT interpretation.

Statistical Notes

The likelihood ratio test (LRT) assumes that:

1. The models are nested (the reduced model is a special case of the full model obtained by constraining parameters).

2. Both models are fit to identical data.

3. Under the null hypothesis, the LR statistic follows a chi-square distribution with degrees of freedom equal to the difference in the number of parameters.

Important caveat for mixed-effects models: When variance components are tested at the boundary (e.g., testing whether a random effect variance is zero), the standard chi-square distribution is not appropriate. The correct null distribution is a mixture of chi-squares (Stram & Lee, 1994). The p-values reported here use the standard chi-square approximation, which is conservative (p-values are too large) for boundary tests.

This function does not automatically verify that models are nested. Users should ensure models are properly nested before interpreting LRT p-values.

References

Stram, D. O., & Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171-1177.

See Also

compare_hurdle_models() for the legacy hurdle-specific comparison

Examples

data(apt)
fit2 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0"))
fit3 <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id",
                          random_effects = c("zeros", "q0", "alpha"))
compare_models(fit2, fit3)

Description

Computes confidence intervals for Q0, alpha, and k parameters from individual demand curve fits. Uses asymptotic normal approximation based on standard errors when available.

Usage

## S3 method for class 'beezdemand_fixed'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

Details

For beezdemand_fixed objects, confidence intervals are computed using the asymptotic normal approximation: estimate +/- z * SE. If standard errors are not available for a parameter, the confidence bounds will be NA.

When the underlying NLS fit objects are available (from detailed = TRUE), this method attempts to use nlstools::confint2() for more accurate profile-based intervals.

Value

A tibble with columns: id, term, estimate, conf.low, conf.high, level.

Examples

fit <- fit_demand_fixed(apt, equation = "hs", k = 2)
confint(fit)
confint(fit, level = 0.90)
confint(fit, parm = "Q0")

Description

Computes confidence intervals for fixed effect parameters from a TMB-based hurdle demand model using the asymptotic normal approximation.

Usage

## S3 method for class 'beezdemand_hurdle'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("internal", "natural"),
  ...
)

Arguments

Details

Confidence intervals are computed using the asymptotic normal approximation based on standard errors from TMB::sdreport(). For parameters estimated on the log scale (Q0, alpha, k), intervals can be back-transformed to the natural scale using report_space = "natural".

The transformation uses:

• For log-scale parameters: exp(estimate +/- z * SE)

Value

A tibble with columns: term, estimate, conf.low, conf.high, level, component, estimate_scale.

Examples

data(apt)
fit <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
confint(fit)

Description

Computes confidence intervals for fixed effect parameters from an NLME-based mixed-effects demand model.

Usage

## S3 method for class 'beezdemand_nlme'
confint(object, parm = NULL, level = 0.95, method = c("wald", "profile"), ...)

Arguments

Details

For Wald intervals, confidence bounds are computed as estimate ± z * SE using standard errors from the model summary.

For profile intervals, nlme::intervals() is called on the underlying nlme model object. This method provides more accurate intervals but can be computationally intensive for complex models.

Value

A tibble with columns: term, estimate, conf.low, conf.high, level, component.

Examples

data(ko)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", equation_form = "zben")
confint(fit)

Description

Confidence Intervals for TMB Model Parameters

Usage

## S3 method for class 'beezdemand_tmb'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("internal", "natural"),
  method = c("wald", "simulate"),
  R = 1000L,
  seed = NULL,
  ...
)

Arguments

Details

method = "simulate" is Monte Carlo simulation from the asymptotic Gaussian sampling approximation to the MLE (neither a Bayesian posterior, nor a data-resampling bootstrap, nor a profile-likelihood interval). Because the sampled distribution is the same Gaussian that Wald assumes, the simulated per-coefficient quantiles converge to the Wald intervals as R -> Inf; the method does not improve on Wald at boundary cases and offers no positivity guarantee on the internal scale (⁠logsigma_*⁠ intervals can be negative). Its value is (a) a diagnostic side-by-side check on the Gaussian approximation, and (b) a shared draw primitive (.tmb_parametric_draws()) for derived-metric confidence intervals.

Value

A tibble with term, estimate, conf.low, conf.high, level.

See Also

confint.beezdemand_nlme(), vcov.beezdemand_tmb().

Examples

data(apt)
fit <- fit_demand_tmb(apt, equation = "exponential", verbose = 0)
confint(fit)
confint(fit, report_space = "natural")
# Diagnostic Monte Carlo intervals (asymptotically Wald-equivalent):
confint(fit, method = "simulate", R = 1000, seed = 42)

Description

Wald confidence intervals via stats::confint.default().

Usage

## S3 method for class 'cp_model_lm'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

Value

A tibble with columns term, estimate, conf.low, conf.high, level, method.

Description

Wraps lme4::confint.merMod(). Defaults to Wald intervals; method = "profile" or method = "boot" use the corresponding lme4 backends and are slower.

Usage

## S3 method for class 'cp_model_lmer'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = c("Wald", "profile", "boot"),
  ...
)

Arguments

Value

A tibble with columns term, estimate, conf.low, conf.high, level, method.

Description

Computes confidence intervals for parameters from a nonlinear cross-price demand model using nlstools::confint2().

Usage

## S3 method for class 'cp_model_nls'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = c("asymptotic", "profile"),
  ...
)

Arguments

Details

This method wraps nlstools::confint2() to provide confidence intervals for the log10-parameterized coefficients (log10_qalone, I, log10_beta).

For back-transformed natural-scale confidence intervals, apply the transformation: 10^conf.low and 10^conf.high for log10-scale parameters.

Value

A tibble with columns: term, estimate, conf.low, conf.high, level, method.

Examples

data(etm)
fit <- fit_cp_nls(etm, equation = "exponentiated")
confint(fit)

Description

A small illustrative dataset of price (x) and consumption (y) target (target), and group (group).

Usage

cp

Format

A data frame with N rows and columns:

id

unique identifier

x

price of the alternative

y

consumption/demand

target

target (e.g., alone, own, alt)

group

e.g., drug or product

Description

Performs pairwise comparisons of intercepts between groups in a cross-price demand model, but only when a significant interaction is present. The emmeans table showing estimated marginal means for intercepts is always returned.

Usage

cp_posthoc_intercepts(object, alpha = 0.05, adjust = "tukey", ...)

Arguments

Value

List containing the emmeans table and optionally pairwise comparisons if interaction is significant

Examples

data(etm)
fit <- fit_cp_linear(etm, type = "mixed", group_effects = TRUE)
cp_posthoc_intercepts(fit)

Description

Performs pairwise comparisons of slopes between groups in a cross-price demand model, but only when a significant interaction is present. The emmeans table showing estimated marginal means for slopes is always returned.

Usage

cp_posthoc_slopes(object, alpha = 0.05, adjust = "tukey", ...)

Arguments

Value

List containing the emmeans table and optionally pairwise comparisons if interaction is significant

Examples

data(etm)
fit <- fit_cp_linear(etm, type = "mixed", group_effects = TRUE)
cp_posthoc_slopes(fit)

Description

A dataset containing ETM data for a small number of participants

Usage

etm

Format

Long-form data.frame with columns: id, x, y, target, group. Participants were asked how many cigarettes, e-cigarettes, combustible, and non-combustible products they would buy at various prices.

Description

A convenience function to extract coefficients from any type of cross-price demand model in a unified format. For mixed effects models, returns a list with different coefficient types.

Usage

extract_coefficients(object, ...)

Arguments

Value

For cp_model_nls and cp_model_lm, returns the model coefficients. For cp_model_lmer, returns a list with fixed, random, and combined coefficients.

Examples

data(etm, package = "beezdemand")
fit <- fit_cp_nls(etm, equation = "exponentiated")
extract_coefficients(fit)

Description

Extra Sum of Squares F-test

Usage

ExtraF(
  dat,
  equation = "hs",
  groups = NULL,
  verbose = FALSE,
  k,
  compare = "alpha",
  idcol = "id",
  xcol = "x",
  ycol = "y",
  groupcol = NULL,
  start_alpha = 0.001
)

Arguments

Details

One alpha better than individual alphas?

Value

List of results and models

Author(s)

Brent Kaplan [email protected]

Examples

## Compare two groups using equation by Koffarnus et al., 2015 and a fixed k of 2

apt$group <- NA
apt[apt$id %in% sample(unique(apt$id), length(unique(apt$id))/2), "group"] <- "a"
apt$group[is.na(apt$group)] <- "b"
ExtraF(apt, "koff", k = 2, groupcol = "group")

Description

Bisection search over n_subjects for the smallest N whose Monte Carlo power estimate from power_demand() reaches target_power. The search accounts for Monte Carlo noise: at each evaluated N, replicates are added in batches (up to n_sim_max) until the Wilson interval for power lies wholly above or below the target; if it still straddles the target at the cap, the decision falls back to the point estimate and the result is flagged uncertain. The selected N and its lower neighbor are then re-evaluated with fresh replicates before minimality is claimed.

The returned n is an estimated minimum under Monte Carlo uncertainty rather than an exact bound. For grant-quality reporting, rerun power_demand() at the returned n with a large n_sim (2000+) and report that estimate with its Monte Carlo confidence interval.

Monotonicity assumption. Bisection presumes that power is non-decreasing in n_subjects. Because every evaluated N is judged from its own independent replicates (and a convergence-conditioned denominator), a Monte Carlo fluctuation can make a lower N read "below" when its true power is above the target, so the search may step past a lower crossing that it never revisits. Evaluated N that contradict the assumption (a lower N reading "above" the selected N) demote the status to "uncertain"; N that were never evaluated cannot be checked. Widen n_sim/n_sim_max when the reported n matters. When the target is already met at n_range[1], that bound is likewise re-evaluated with fresh replicates before "at_lower_bound" is reported; if the second look does not clear the target the bound is treated as below and the bisection proceeds upward.

Usage

find_n_demand(
  target_power = 0.8,
  effect = list(delta_q0 = NULL, delta_alpha = NULL),
  design = list(),
  n_range = c(6, 200),
  n_sim = 200,
  n_sim_max = 4 * n_sim,
  alpha = 0.05,
  df = NULL,
  seed = NULL,
  equation = "simplified",
  random_effects = NULL,
  multi_start = FALSE,
  verbose = TRUE,
  design_type = c("within", "between"),
  ...
)

Arguments

Value

An object of class beezdemand_power_n: a list with

n

Estimated smallest n_subjects reaching target_power; NA when the confirmation pass contradicted the search (status = "unresolved").

target_power

As supplied.

status

"confirmed" (selected N re-confirmed above target and N - 1 below), "uncertain" (a decision relied on a point estimate, N - 1 also cleared the target on reconfirmation, or an evaluated lower N read above the target, so the returned N may not be minimal), "unresolved" (the selected N failed its own reconfirmation; n is NA), or "at_lower_bound" (the target was already met at n_range[1] on two independent looks; smaller N was not explored; widen n_range downward if that matters). These labels describe a heuristic Monte Carlo decision rule (repeated looks at ordinary Wilson intervals across several N) rather than a formal sequential error guarantee.

uncertain

Logical; TRUE when any search decision was made on a point estimate rather than a conclusive Wilson interval, or the status is not "confirmed"/"at_lower_bound".

evaluations

Tibble of every evaluation: n_subjects, n_sim_total, n_used, usable_fraction, power, ci_lower, ci_upper, decision. A warning fires when any evaluation had fewer than 95% usable fits.

alpha, df, effect, design, n_range, n_sim, n_sim_max, seed, settings, call

Echoed inputs and effective settings.

See Also

power_demand() for the Monte Carlo engine.

Other power-analysis: power_demand()

Examples

# Tiny search for demonstration only (use n_sim >= 200 and a wide
# n_range for real planning; see vignette("power-analysis"))
res <- find_n_demand(
  target_power = 0.8,
  effect = list(delta_q0 = log(2.5)),
  n_range = c(4, 8), n_sim = 5, n_sim_max = 10, seed = 1, verbose = FALSE
)
print(res)

Description

Fit a Linear Cross-Price Demand Model

Usage

fit_cp_linear(
  data,
  type = c("fixed", "mixed"),
  x_var = "x",
  y_var = "y",
  id_var = "id",
  group_var = "group",
  target_var = "target",
  filter_target = TRUE,
  target_level = "alt",
  formula = NULL,
  log10x = FALSE,
  group_effects = FALSE,
  random_slope = FALSE,
  return_all = TRUE,
  ...
)

fit_cp_linear.default(
  data,
  formula = NULL,
  log10x = FALSE,
  return_all = FALSE,
  ...
)

fit_cp_linear.mixed(
  data,
  formula = NULL,
  log10x = FALSE,
  return_all = FALSE,
  ...
)

Arguments

Value

Fitted linear model.

Examples

data(etm)
## Fixed-effects linear cross-price model
fit_fixed <- fit_cp_linear(etm, type = "fixed", group_effects = TRUE)
summary(fit_fixed)

## Mixed-effects linear cross-price model
fit_mixed <- fit_cp_linear(etm, type = "mixed", group_effects = TRUE)
summary(fit_mixed)

Description

Fits a cross-price demand curve using log10-parameterization for numerical stability. The optimizer estimates parameters on the log10 scale where applicable, so positivity constraints hold automatically.

Equation forms:

• Exponentiated (default):

  $y = Q_{alone} \cdot 10^{I \cdot \exp(-\beta \cdot x)}$

• Exponential (fits on log10 response scale):

  $\log_{10}(y) = \log_{10}(Q_{alone}) + I \cdot \exp(-\beta \cdot x)$

• Additive (level on $y$):

  $y = Q_{alone} + I \cdot \exp(-\beta \cdot x)$

where $x$ is the alternative product price (or "cross" price) and $y$ is consumption of the target good.

Optimizer parameters (log10 parameterization):

• log10_qalone: $\log_{10}(Q_{alone})$ - baseline consumption when the alternative is effectively absent.

• I: cross-price interaction intensity; sign and magnitude reflect substitution/complementarity. Unconstrained (can be negative for substitutes).

• log10_beta: $\log_{10}(\beta)$ - rate at which cross-price influence decays as $x$ increases.

Natural-scale values are recovered as $Q_{alone} = 10^{log10\_qalone}$ and $\beta = 10^{log10\_beta}$.

The function first attempts a multi-start nonlinear least squares fit (nls.multstart). If that fails (or if explicit start_values are provided), it falls back to minpack.lm::nlsLM. Optionally, it will make a final attempt with nlsr::wrapnlsr. Returns either the fitted model or a structured object with metadata for downstream methods.

Usage

fit_cp_nls(
  data,
  equation = c("exponentiated", "exponential", "additive"),
  x_var = "x",
  y_var = "y",
  start_values = NULL,
  start_vals = lifecycle::deprecated(),
  iter = 100,
  bounds = NULL,
  fallback_to_nlsr = TRUE,
  return_all = TRUE
)

Arguments

Details

Start values. When start_values is missing, the function: (1) estimates a reasonable range for log10_qalone from the observed y, (2) estimates log10_beta from the price range, and (3) launches a multi-start grid in nls.multstart.

Zero handling for exponential equation. Since the exponential equation fits on the $\log_{10}(y)$ scale, observations with $y \le 0$ are automatically removed with a warning. Use the exponentiated or additive forms if you need to retain zero consumption values.

Fitting pipeline (short-circuiting):

1. nls.multstart::nls_multstart() with random starts.

2. If that fails (or if start_values provided): minpack.lm::nlsLM() using start_values (user or internally estimated).

3. If that fails and fallback_to_nlsr = TRUE: nlsr::wrapnlsr().

The returned object has class "cp_model_nls" (when return_all = TRUE) with components: model, method (the algorithm used), equation, start_vals, nlsLM_fit, nlsr_fit, and the data used. This is convenient for custom print/summary/plot methods.

Value

If return_all = TRUE (default): a list of class "cp_model_nls":

• model: the fitted object from the successful backend.

• method: one of "nls_multstart", "nlsLM", or "wrapnlsr".

• equation: the model family used.

• start_vals: named list of starting values (final used; kept for backward compatibility).

• nlsLM_fit, nlsr_fit: fits from later stages (if attempted; NULL when that stage was not reached, or was attempted but did not itself produce a fit).

• data: the 2-column data frame actually fit.

• convergence: list with isConv, finIter, stopCode, stopMessage for the winning backend, read from model$convInfo; populated for nls/nlsLM-class fits, and for nlsr::wrapnlsr() fits too when it returns a plain nls-class object (its usual successful case); isConv = NA only when the winning fit carries no convInfo at all.

If return_all = FALSE: the fitted model object from the successful backend.

Convergence & warnings

• summary() and confint() on the returned object emit a classed warning (beezdemand_cp_nls_nonconverged_warning) when convergence$isConv is explicitly FALSE for the winning fit; unreported convergence (NA, when the winning object carries no convInfo) is silent.

• Poor scaling or extreme y dispersion can make parameters weakly identified.

• For "exponential", the model fits on the $\log_{10}(y)$ scale internally.

See Also

check_unsystematic_cp for pre-fit data screening, validate_cp_data for input validation.

Examples

## --- Example: Real data (E-Cigarettes, id = 1) ---
dat <- structure(list(
  id = c(1, 1, 1, 1, 1, 1),
  x = c(2, 4, 8, 16, 32, 64),
  y = c(3, 5, 5, 16, 17, 13),
  target = c("alt", "alt", "alt", "alt", "alt", "alt"),
  group = c("E-Cigarettes", "E-Cigarettes", "E-Cigarettes",
            "E-Cigarettes", "E-Cigarettes", "E-Cigarettes")
), row.names = c(NA, -6L), class = c("tbl_df", "tbl", "data.frame"))

## Fit the default (exponentiated) cross-price form
fit_ecig <- fit_cp_nls(dat, equation = "exponentiated", return_all = TRUE)
summary(fit_ecig)         # model summary
fit_ecig$method           # backend actually used (e.g., "nls_multstart")
coef(fit_ecig$model)      # parameter estimates: log10_qalone, I, log10_beta

Description

Modern interface for fitting individual demand curves via nonlinear least squares. Returns a structured S3 object with standard methods including summary(), tidy(), and glance().

Usage

fit_demand_fixed(
  data,
  equation = c("hs", "koff", "simplified", "linear", "exponential", "exponentiated"),
  k = 2,
  agg = NULL,
  x_var = "x",
  y_var = "y",
  id_var = "id",
  param_space = c("natural", "log10"),
  by = NULL,
  multistart = TRUE,
  S = NULL,
  ...
)

Arguments

Details

This function is a modern wrapper around the legacy FitCurves() function. It provides the same fitting capabilities but returns a structured S3 object with standardized methods for model interrogation.

Multi-start rescue protocol

fit_demand_fixed() always runs FitCurves()'s existing heuristic start exactly as before (the "production start"). A subject whose production fit is strict-converged (converged_strict: the optimizer's own convergence flag AND finite coefficients/objective AND not sitting on a user-supplied bound) is accepted immediately; no sampled starts are ever run for it, so its row, fitted model, predictions, and data are byte-identical to the multistart = FALSE / S = 1 protocol by construction. Only subjects whose production fit is not strict-converged are re-fit from S - 1 additional starts, sampled log-uniformly in interpretable (Q0, Pmax) coordinates and mapped to each equation's native (Q0, alpha) parameterization via the same closed forms used by beezdemand_calc_pmax_omax(). Among the sampled attempts that themselves strict-converge, the minimum-residual-SS start wins (ties broken by draw order). If none of the sampled starts strict-converge, the original (non-converged) production row is kept. equation = "linear" is a closed-form fit and is never multistarted. FitCurves() itself is unchanged; sampling draws from the ambient RNG stream (call set.seed() before fit_demand_fixed() for reproducibility; the helpers never call set.seed() themselves).

Value

An object of class beezdemand_fixed with components:

results

Data frame of fitted parameters for each subject. Gains n_starts_tried, n_starts_converged, and start_source ("production", "sampled", or "none") from the multi-start protocol; see Details.

fits

List of model fit objects (if detailed = TRUE internally)

predictions

List of prediction data frames

data_used

List of data frames used for each fit

multistart

List describing the multi-start protocol: multistart, S (resolved budget), equation, eligible (whether this equation supports rescue), and summary (per-subject start metadata, or NULL when not applicable)

call

The original function call

equation

The equation form used

k_spec

Description of k specification

agg

Aggregation method used

n_total

Total number of subjects/fits attempted

n_success

Number of successful fits

n_fail

Number of failed fits

See Also

fit_demand_tmb() for TMB mixed-effects models, fit_demand_mixed() for NLME mixed-effects models, fit_demand_hurdle() for hurdle models.

Other demand-fitting: fit_demand_hurdle(), fit_demand_mixed(), fit_demand_tmb()

Examples

data(apt)
fit <- fit_demand_fixed(apt, equation = "hs", k = 2)
print(fit)
summary(fit)
tidy(fit)
glance(fit)

# Grouped analysis: fit separately by gender (subset keeps it fast)
data(apt_full)
ids <- unique(apt_full[c("id", "gender")])
ids <- ids[ids$gender %in% c("Male", "Female"), ]
keep <- unlist(lapply(split(ids$id, ids$gender), head, 40))
dat <- apt_full[apt_full$id %in% keep, ]
fit_g <- fit_demand_fixed(dat, equation = "hs", k = 2, by = "gender")
tidy(fit_g)   # group column prepended
glance(fit_g)  # one row per group

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

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

re_cov_fallback

Logical. TRUE when the estimated random-effects covariance was not positive definite and the subject-level effects were computed from an uncorrelated (diagonal) approximation; the reported correlations do not apply to those outputs. summary() and print() carry a note when it is TRUE.

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

The default part2 = "zhao_exponential" places no $Q_0$ inside the exponent ($\exp(-\alpha_{Zhao} \, P)$), whereas the Hursh & Silberberg form used by fit_demand_fixed(equation = "hs") and fit_demand_tmb(equation = "exponential") uses $\exp(-\alpha_{HS} \, Q_0 \, P)$. The two $\alpha$ values are therefore on different scales and relate as

$\alpha_{HS} = \alpha_{Zhao} / Q_0,$

so a Zhao $\alpha$ is not directly comparable with an HS $\alpha$ unless divided by the subject's (or group's) $Q_0$. Use part2 = "exponential" (alias "hs_stdq0") for an HS-scaled $\alpha$.

Convergence

The 3-random-effect spec (random_effects = c("zeros", "q0", "alpha")) can false-converge on real purchase task data: nlminb() reports a nonzero convergence code (e.g. "false convergence (8)") while the Hessian is not positive definite (fit$hessian_pd == FALSE). This reflects weak identification of the alpha random effect rather than a broken model specification (companion datasets from the same source can converge cleanly with the identical spec). fit$opt$message and fit$opt$convergence retain nlminb's own diagnostic fields alongside fit$converged and fit$hessian_pd; print() and summary() display a prominent warning block (quoting the optimizer message) whenever either signals a problem, so a non-converged or non-PD-Hessian fit is hard to miss even without inspecting these fields directly.

The recommended stability check when the 3RE spec reports non-convergence or a non-PD Hessian is to refit with random_effects = c("zeros", "q0") (dropping the alpha random effect, which converges more readily) and compare the empirical-Bayes subject parameters between the two fits. Broadly similar per-subject estimates support treating the 2RE fit's conclusions as stable across the two specifications. This is a diagnostic comparison rather than a claim that the 2RE spec is generally preferred, and beezdemand does not refit automatically.

See Also

summary.beezdemand_hurdle, predict.beezdemand_hurdle (and its Scoring predictions section: score zero-inclusive consumption with type = "demand", the default since 0.3.0), plot.beezdemand_hurdle, compare_hurdle_models, simulate_hurdle_data

fit_demand_tmb() for continuous mixed-effects models, fit_demand_mixed() for NLME-based fitting, fit_demand_fixed() for individual NLS curves.

Other demand-fitting: fit_demand_fixed(), fit_demand_mixed(), fit_demand_tmb()

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)

Description

Fits a nonlinear mixed-effects model for behavioral economic demand data. The function allows Q0 and alpha parameters to vary by specified factors and supports different demand equation forms.

Usage

fit_demand_mixed(
  data,
  y_var,
  x_var,
  id_var,
  factors = NULL,
  factor_interaction = FALSE,
  equation_form = c("zben", "simplified", "exponentiated"),
  param_space = c("log10", "natural"),
  k = NULL,
  custom_model_formula = NULL,
  fixed_rhs = NULL,
  continuous_covariates = NULL,
  start_value_method = c("heuristic", "pooled_nls"),
  random_effects = Q0 + alpha ~ 1,
  covariance_structure = c("pdDiag", "pdSymm"),
  start_values = NULL,
  collapse_levels = NULL,
  nlme_control = list(msMaxIter = 200, niterEM = 100, maxIter = 200, pnlsTol = 0.001,
    tolerance = 1e-06, apVar = TRUE, minScale = 1e-09, opt = "nlminb", msVerbose = FALSE),
  method = "ML",
  ...
)

Arguments