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

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. This enables easy model diagnostics and visualization with the tidyverse.

Usage

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

Arguments

Details

For two-part hurdle models:

• .fitted gives predicted demand on the natural consumption scale

• .fitted_prob gives the predicted probability of positive consumption

• .resid is defined only for positive observations as log(y) - .fitted_link

• Observations with zero consumption have .resid = NA since they are explained by Part I (the zero-probability component), not Part II

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 on log scale for positive observations, NA for zeros

.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

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 (handles zeros gracefully)
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

• Uses modern beezdemand styling

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

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

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

• 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

residuals

Summary statistics for residuals

random_effects

Summary of random effects (if applicable)

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

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

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)

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

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

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

This function 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

This function 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

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, ensuring positive constraints are naturally satisfied.

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

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

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

Convergence & warnings

• Check convergence codes and residual diagnostics from the underlying fit.

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

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.

Value

An object of class beezdemand_fixed with components:

results

Data frame of fitted parameters for each subject

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

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

Examples

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

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

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)

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

list(
  Q0 = list(factor_name = list(new_level = c(old_levels), ...)),
  alpha = list(factor_name = list(new_level = c(old_levels), ...))
)

Value

An object of class beezdemand_nlme.

Examples

# Basic mixed-effects demand fit with apt data
# Transform consumption using LL4 for the zben equation
apt_ll4 <- apt |> dplyr::mutate(y_ll4 = ll4(y))

fit <- fit_demand_mixed(
  data = apt_ll4,
  y_var = "y_ll4",
  x_var = "x",
  id_var = "id",
  equation_form = "zben"
)
print(fit)
summary(fit)

Description

Analyzes purchase task data

Usage

FitCurves(
  dat,
  equation,
  k,
  agg = NULL,
  detailed = FALSE,
  xcol = "x",
  ycol = "y",
  idcol = "id",
  groupcol = NULL,
  lobound,
  hibound,
  constrainq0 = NULL,
  startq0 = NULL,
  startalpha = NULL,
  param_space = c("natural", "log10")
)

Arguments

Details

FitCurves() has been superseded by fit_demand_fixed(), which provides a modern S3 interface with standardized methods (summary(), tidy(), glance(), predict()). FitCurves() will continue to work but is no longer recommended for new code. See vignette("migration-guide") for migration instructions.

Value

If detailed == FALSE (default), a dataframe of results. If detailed == TRUE, a 3 element list consisting of (1) dataframe of results, (2) list of model objects, (3) list of individual dataframes used in fitting

Author(s)

Brent Kaplan [email protected] Shawn Gilroy [email protected]

See Also

fit_demand_fixed() for the modern interface

Examples

## Analyze using Hursh & Silberberg, 2008 equation with a k fixed to 2
FitCurves(apt[sample(apt$id, 5), ], "hs", k = 2)

Description

Fits curve to pooled or mean data

Usage

FitMeanCurves(
  dat,
  equation,
  k,
  remq0e = FALSE,
  replfree = NULL,
  rem0 = FALSE,
  nrepl = NULL,
  replnum = NULL,
  plotcurves = FALSE,
  method = NULL,
  indpoints = TRUE,
  vartext = NULL
)

Arguments

Details

FitMeanCurves() has been superseded by fit_demand_fixed() with the agg parameter. FitMeanCurves() will continue to work but is no longer recommended for new code. See vignette("migration-guide") for migration instructions.

Value

Data frame

Author(s)

Brent Kaplan [email protected]

See Also

fit_demand_fixed() for the modern interface with agg parameter

Examples

## Fit aggregated data (mean only) using Hursh & Silberberg, 2008 equation with a k fixed at 2
FitMeanCurves(apt[sample(apt$id, 5), ], "hs", k = 2, method = "Mean")

Description

Modern wrapper for fitting individual demand curves via nonlinear least squares. Returns a structured S3 object with standard methods.

Description

S3 method for fixef for objects of class beezdemand_nlme. Extracts the fixed-effect coefficients from the fitted nlme model.

Usage

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

Arguments

Value

A named numeric vector of fixed-effect coefficients.

See Also

coef.beezdemand_nlme, ranef.beezdemand_nlme

Description

Extract Fixed Effects from Mixed-Effects Cross-Price Model

Usage

## S3 method for class 'cp_model_lmer'
fixef(object, ...)

Arguments

Value

Named vector of fixed effects

Description

Conducts pairwise comparisons for Q0 and/or alpha parameters from a beezdemand_nlme model across levels of specified factors. Comparisons are performed on the log10 scale of the parameters. Results include estimates of differences (on log10 scale) and optionally, ratios (on the natural scale by applying 10^difference).

Usage

get_demand_comparisons(
  fit_obj,
  params_to_compare = c("Q0", "alpha"),
  compare_specs = NULL,
  contrast_type = "pairwise",
  contrast_by = NULL,
  adjust = "tukey",
  at = NULL,
  ci_level = 0.95,
  report_ratios = TRUE,
  ...
)

Arguments

Value

A list named by parameter. Each element contains:

S3 class beezdemand_comparison is assigned.

Examples

data(ko, package = "beezdemand")
ko$y_ll4 <- ll4(ko$y, lambda = 4)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
  id_var = "monkey", factors = "dose", equation_form = "zben")
get_demand_comparisons(fit)

Description

Calculates Estimated Marginal Means (EMMs) for Q0 and alpha parameters from a beezdemand_nlme model for all combinations of specified factor levels. Reports parameters on both their estimation scale (log10) and their natural, back-transformed scale. Optionally includes Essential Value (EV).

Usage

get_demand_param_emms(
  fit_obj,
  factors_in_emm = NULL,
  at = NULL,
  ci_level = 0.95,
  include_ev = FALSE,
  ...
)

Arguments

Value

A tibble containing:

Examples

data(ko, package = "beezdemand")
ko$y_ll4 <- ll4(ko$y, lambda = 4)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
  id_var = "monkey", factors = "dose", equation_form = "zben")
get_demand_param_emms(fit)

Description

Computes the trend (slope) of Q0 and/or alpha with respect to one or more continuous covariates using emmeans::emtrends() on a fitted beezdemand_nlme model. Trends are computed on the parameter estimation scale (log10), consistent with how parameters are modeled.

Usage

get_demand_param_trends(
  fit_obj,
  params = c("Q0", "alpha"),
  covariates,
  specs = ~1,
  at = NULL,
  ci_level = 0.95,
  ...
)

Arguments

Value

A tibble combining trends for each requested parameter and covariate, including columns for grouping factors (from specs), parameter, covariate, trend (slope on log10 scale), and its CI (lower.CL, upper.CL).

Examples

data(ko)
ko$dose_num <- as.numeric(as.character(ko$dose))
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", factors = "drug",
                        equation_form = "zben")
trends <- get_demand_param_trends(fit, covariates = "dose_num",
                                  specs = ~ drug)

Description

Calculates summary statistics for consumption data at each price point, including measures of central tendency (mean, median), variability (SD), range (min, max), and data quality (proportion of zeros, missing values).

This is the modern replacement for GetDescriptives(), returning a structured S3 object with dedicated methods for printing, summarizing, and visualizing.

Usage

get_descriptive_summary(data, x_var = "x", y_var = "y", id_var = "id")

Arguments

Details

For each unique price in the dataset, the function calculates:

• Mean - Average consumption across subjects (rounded to 2 decimals)

• Median - Median consumption (rounded to 2 decimals)

• SD - Standard deviation (rounded to 2 decimals)

• PropZeros - Proportion of subjects with zero consumption (0-1)

• NAs - Count of missing values

• Min - Minimum consumption value (rounded to 2 decimals)

• Max - Maximum consumption value (rounded to 2 decimals)

Value

An S3 object of class beezdemand_descriptive containing:

• statistics - Data frame with 8 columns (Price, Mean, Median, SD, PropZeros, NAs, Min, Max) and one row per unique price

• call - The matched call

• data_summary - List with n_subjects, n_prices, and prices vector

See Also

• GetDescriptives() - Legacy function (superseded)

• plot.beezdemand_descriptive() - Visualization method

• summary.beezdemand_descriptive() - Extended summary

Examples

data(apt, package = "beezdemand")

# Calculate descriptive statistics
desc <- get_descriptive_summary(apt)
print(desc)

# View statistics table
desc$statistics

# Create visualization
plot(desc)

# Extended summary with distribution info
summary(desc)

Description

Calculates empirical (model-free) measures of demand from purchase task data. These metrics characterize consumption patterns without fitting a demand curve model.

This is the modern replacement for GetEmpirical(), returning a structured S3 object with dedicated methods for printing, summarizing, and visualizing.

Usage

get_empirical_measures(data, x_var = "x", y_var = "y", id_var = "id")

Arguments

Details

Empirical Measures

• Intensity - The consumption value at the lowest price point. Reflects unrestricted demand or preferred level of consumption.

• BP0 (Breakpoint 0) - The first price at which consumption drops to zero. If consumption never reaches zero, BP0 = NA. Indicates the price threshold where the commodity becomes unaffordable or undesirable.

• BP1 (Breakpoint 1) - The highest price at which consumption is still non-zero. If all consumption is zero, BP1 = NA. Represents the upper limit of the commodity's value to the consumer.

• Omaxe (Empirical Omax) - The maximum observed expenditure across all prices (max of price × consumption). Represents peak spending on the commodity.

• Pmaxe (Empirical Pmax) - The price at which maximum expenditure occurs. If maximum expenditure is zero, Pmaxe = 0. If multiple prices have the same maximum expenditure, the highest price is returned.

Value

An S3 object of class beezdemand_empirical containing:

• measures - Data frame with one row per subject and columns:

  • id - Subject identifier

  • Intensity - Consumption at lowest price (demand intensity)

  • BP0 - Breakpoint 0: first price where consumption = 0

  • BP1 - Breakpoint 1: last price with non-zero consumption

  • Omaxe - Empirical Omax: maximum total expenditure (price × consumption)

  • Pmaxe - Empirical Pmax: price at which maximum expenditure occurs

• call - The matched call

• data_summary - List with n_subjects, has_zeros, and complete_cases

Note

• Data must not contain duplicate prices within a subject (will error)

• Missing values (NA) in consumption are preserved in calculations where applicable

• Breakpoints require at least one zero consumption value to be meaningful

See Also

• GetEmpirical() - Legacy function (superseded)

• plot.beezdemand_empirical() - Visualization method

• summary.beezdemand_empirical() - Extended summary

Examples

data(apt, package = "beezdemand")

# Calculate empirical measures
emp <- get_empirical_measures(apt)
print(emp)

# View measures table
emp$measures

# Extended summary with distribution info
summary(emp)

# Visualize distribution of measures
plot(emp)  # histogram by default
plot(emp, type = "matrix")  # scatterplot matrix

Description

Provides summary statistics for subject-level demand parameters from a hurdle demand model. This is analogous to EMMs but based on empirical Bayes estimates of subject-specific parameters.

Usage

get_hurdle_param_summary(fit_obj, ci_level = 0.95)

Arguments

Value

A data frame with summary statistics for each parameter:

parameter

Parameter name

mean

Mean across subjects

sd

Standard deviation across subjects

median

Median across subjects

lcl

Lower confidence limit (based on percentiles)

ucl

Upper confidence limit (based on percentiles)

min

Minimum value

max

Maximum value

See Also

fit_demand_hurdle, get_subject_pars

Examples

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

Description

This function extracts and combines fixed and random effects to calculate individual-level predicted coefficients for all parameter-factor combinations from a beezdemand_nlme model object. It automatically detects the factor structure and calculates coefficients for each individual and factor level.

Usage

get_individual_coefficients(
  fit_obj,
  params = c("Q0", "alpha"),
  format = c("wide", "long")
)

Arguments

Details

Individual-level coefficients represent the predicted parameter values for each subject in the study. For models with factors, these coefficients combine:

1. The baseline intercept effect (fixed + random)

2. The factor-specific effect (fixed + random) for each factor level

This is equivalent to manually calculating: coefficient = intercept_fixed + intercept_random + factor_fixed + factor_random

The function automatically handles:

• Models with or without factors

• Any number of factor levels

• Missing random effects (defaults to 0)

• Complex factor structures with multiple factors

For models without factors, only intercept coefficients are calculated. For models with factors, both intercept and factor-level coefficients are provided.

Value

A data frame with individual-level predicted coefficients.

• In "wide" format: rows are individuals, columns are parameter-factor combinations

• In "long" format: columns are id, parameter, condition, coefficient_value

Column naming convention for wide format:

• ⁠estimated_\{param\}_intercept⁠: Baseline/reference level coefficient

• ⁠estimated_\{param\}_\{factor\}\{level\}⁠: Factor level-specific coefficient

All coefficients are on the log10 scale (same as model estimation scale). To convert to natural scale, use 10^coefficient.

See Also

fit_demand_mixed for fitting the original model coef.beezdemand_nlme for extracting model coefficients get_demand_param_emms for estimated marginal means

Examples

data(ko)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", factors = "drug",
                        equation_form = "zben")
individual_coefs <- get_individual_coefficients(fit)
head(individual_coefs)

Description

Calculates the k scaling parameter used in demand curve equations to normalize consumption across different units or ranges. The k value is derived from the logarithmic range of consumption values.

This is the modern replacement for GetK(), with explicit parameters for the adjustment value and optional verbose output.

Usage

get_k(
  data,
  use_means = TRUE,
  adjustment = 0.5,
  x_var = "x",
  y_var = "y",
  verbose = FALSE
)

Arguments

Details

The k parameter is calculated as:

$k = \log_{10}(\text{max}) - \log_{10}(\text{min}) + \text{adjustment}$

where max and min are the maximum and minimum non-zero consumption values.

Use in Demand Equations

The k parameter appears in several demand curve equations:

• Hursh & Silberberg (2008): Scales the exponential term

• Koffarnus et al. (2015): Normalizes the exponentiated model

• Ensures numerical stability during model fitting

Calculation Modes

• use_means = TRUE (default): Calculates k from mean consumption at each price point. Recommended when data has multiple subjects, as it reduces influence of individual outliers.

• use_means = FALSE: Calculates k from the full range of individual consumption values. May be preferable for single-subject data or when individual variability is theoretically important.

Value

A single numeric value representing the k scaling parameter

Note

• Only non-zero consumption values are used (zero values are excluded)

• Missing values (NA) are automatically removed via na.rm = TRUE

• The default adjustment of 0.5 is conventional but can be modified

See Also

• GetK() - Legacy function (superseded)

• FitCurves() - Uses k parameter in demand curve fitting

Examples

data(apt, package = "beezdemand")

# Calculate k using default settings (mean range + 0.5)
k_val <- get_k(apt)

# Calculate k from individual values
k_ind <- get_k(apt, use_means = FALSE)

# Calculate with custom adjustment
k_custom <- get_k(apt, adjustment = 1.0)

# Show calculation details
k_verbose <- get_k(apt, verbose = TRUE)

Description

This function is a wrapper around get_demand_param_emms. It first calls get_demand_param_emms to calculate Estimated Marginal Means (EMMs) for Q0 and alpha parameters over all combinations of the specified factor levels. It then filters these results to return EMMs only for the combinations of factor levels that were actually present in the original dataset used to fit the beezdemand_nlme model.

Usage

get_observed_demand_param_emms(
  fit_obj,
  factors_in_emm = NULL,
  at = NULL,
  ci_level = 0.95,
  include_ev = FALSE,
  ...
)

Arguments

Value

A tibble similar to the output of get_demand_param_emms, but filtered to include only rows corresponding to factor level combinations that were observed in the original fit_obj$data. Contains:

See Also

get_demand_param_emms

Examples

data(ko, package = "beezdemand")
ko$y_ll4 <- ll4(ko$y, lambda = 4)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
  id_var = "monkey", factors = "dose", equation_form = "zben")
get_observed_demand_param_emms(fit)

Description

Convenience function to extract subject-specific demand parameters from a fitted hurdle demand model. Equivalent to accessing object$subject_pars.

Usage

get_subject_pars(object)

Arguments

Value

Data frame with subject-specific parameters including:

id

Subject identifier

a_i

Random effect for Part I (zeros)

b_i

Random effect for Part II (Q0)

c_i

Random effect for alpha (3-RE model only)

Q0

Subject-specific intensity (consumption at price 0)

alpha

Subject-specific elasticity

breakpoint

Price where P(quit) = 0.5

Pmax

Price at maximum expenditure

Omax

Maximum expenditure

See Also

fit_demand_hurdle

Examples

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

Description

...

Usage

GetAnalyticPmax(Alpha, K, Q0)

Arguments

Details

...

Value

Numeric

Author(s)

Shawn Gilroy [email protected]

Examples

GetAnalyticPmax(Alpha = 0.001, K = 3, Q0 = 10)

Description

Fallback method for Analytic Pmax

Usage

GetAnalyticPmaxFallback(K_, A_, Q0_)

Arguments

Details

Derivative-based optimization strategy

Value

numeric

Author(s)

Shawn Gilroy [email protected]

Examples

GetAnalyticPmaxFallback(K_ = 1, A_ = 0.001, Q0_ = 10)

Description

Calculates descriptive statistics from purchase task data.

Usage

GetDescriptives(
  dat,
  bwplot = FALSE,
  outdir = "../plots/",
  device = "png",
  filename = "bwplot"
)

Arguments