Version:	0.2.0
Date:	2026-02-27
Title:	Behavioral Economic Easy Demand
Author:	Brent Kaplan [aut, cre, cph], Shawn Gilroy [ctb]
Maintainer:	Brent Kaplan <bkaplan.ku@gmail.com>
Description:	Facilitates many of the analyses performed in studies of behavioral economic demand. The package supports commonly-used options for modeling operant demand including (1) data screening proposed by Stein, Koffarnus, Snider, Quisenberry, & Bickel (2015; <doi:10.1037/pha0000020>), (2) fitting models of demand such as linear (Hursh, Raslear, Bauman, & Black, 1989, <doi:10.1007/978-94-009-2470-3_22>), exponential (Hursh & Silberberg, 2008, <doi:10.1037/0033-295X.115.1.186>) and modified exponential (Koffarnus, Franck, Stein, & Bickel, 2015, <doi:10.1037/pha0000045>), and (3) calculating numerous measures relevant to applied behavioral economists (Intensity, Pmax, Omax). Also supports plotting and comparing data.
Depends:	R (≥ 4.1.0)
Imports:	nlsr, nlstools, nls2, ggplot2, stats, optimx, broom, lme4, emmeans, minpack.lm, nls.multstart, performance, scales, tibble, lifecycle, dplyr, tidyr, nlme, rlang, utils, TMB
Suggests:	broom.mixed, GGally, knitr, tidyverse, rmarkdown, purrr, conflicted, devtools, here, readr, patchwork, testthat (≥ 3.0.0)
LinkingTo:	TMB, RcppEigen
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://brentkaplan.github.io/beezdemand/, https://github.com/brentkaplan/beezdemand
BugReports:	https://github.com/brentkaplan/beezdemand/issues
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.3
VignetteBuilder:	knitr
SystemRequirements:	GNU make
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2026-03-03 01:10:39 UTC; brent
Repository:	CRAN
Date/Publication:	2026-03-03 09:00:19 UTC

beezdemand: Behavioral Economic Easy Demand

Description

Facilitates many of the analyses performed in studies of behavioral economic demand. The package supports commonly-used options for modeling operant demand including (1) data screening proposed by Stein, Koffarnus, Snider, Quisenberry, & Bickel (2015; doi:10.1037/pha0000020), (2) fitting models of demand such as linear (Hursh, Raslear, Bauman, & Black, 1989, doi:10.1007/978-94-009-2470-3_22), exponential (Hursh & Silberberg, 2008, doi:10.1037/0033-295X.115.1.186) and modified exponential (Koffarnus, Franck, Stein, & Bickel, 2015, doi:10.1037/pha0000045), and (3) calculating numerous measures relevant to applied behavioral economists (Intensity, Pmax, Omax). Also supports plotting and comparing data.

Author(s)

Maintainer: Brent Kaplan bkaplan.ku@gmail.com (ORCID) [copyright holder]

Other contributors:

• Shawn Gilroy shawn.gilroy@temple.edu [contributor]

See Also

Useful links:

• https://brentkaplan.github.io/beezdemand/

• https://github.com/brentkaplan/beezdemand

• Report bugs at https://github.com/brentkaplan/beezdemand/issues

Derived Metrics Registry

Description

Defines canonical names for derived metrics and their legacy mappings.

Usage

.beezdemand_derived_registry

Format

An object of class list of length 10.

Equation Registry

Description

Maps equation IDs to their parameter requirements and scales.

Usage

.beezdemand_equation_registry

Format

An object of class list of length 13.

Canonical Parameter Definitions

Description

A list defining all canonical parameters, their constraints, and scale info.

Usage

.beezdemand_param_registry

Format

An object of class list of length 15.

Compute Normalized Alpha (Alpha Star) via the Delta Method

Description

Implements Strategy B normalization of the demand elasticity parameter \alpha so that values are comparable across different k values (Rzeszutek et al., 2025). The formula is \alpha^* = -\alpha / \ln(1 - 1/(k \cdot \ln(b))) where b is the logarithmic base used by the demand equation (10 for HS/Koff, e for hurdle models).

Usage

.calc_alpha_star(params, param_scales, vcov = NULL, base = c("e", "10"))

Arguments

Details

Standard errors are obtained via the delta method when a variance–covariance matrix (or SE vector) is supplied.

Value

A list with elements:

estimate

Numeric scalar; the alpha_star value, or NA.

se

Numeric scalar; delta-method SE, or NA.

note

Character or NULL; diagnostic message if alpha_star could not be computed.

References

Rzeszutek, M. J., Regnier, S. D., Franck, C. T., & Koffarnus, M. N. (2025). Overviewing the exponential model of demand and introducing a simplification that solves issues of span, scale, and zeros. Experimental and Clinical Psychopharmacology.

Calculate Observed Pmax and Omax for a Single Subject

Description

Calculate Observed Pmax and Omax for a Single Subject

Usage

.calc_observed_pmax_omax(price, consumption)

Arguments

Value

List with observed metrics

Check Unit Elasticity Condition

Description

Check Unit Elasticity Condition

Usage

.check_unit_elasticity(elasticity, tol = 0.1)

Arguments

Value

Logical TRUE if |eta + 1| < tol

Calculate Elasticity at a Given Price

Description

Computes eta(p) = d log(Q(p)) / d log(p) via central finite differences

Usage

.elasticity_at_price(demand_fn, price, delta = 1e-05)

Arguments

Value

Numeric elasticity value

Build TMB Data List for Hurdle Model

Description

Internal function to construct the TMB data list for hurdle demand models.

Usage

.hurdle_build_tmb_data(prepared_data, model_name, epsilon)

Arguments

Value

A list suitable for TMB::MakeADFun(data = ...).

Compute Subject-Specific Parameters

Description

Internal function to compute subject-specific demand parameters from fixed effects and random effects.

Usage

.hurdle_compute_subject_pars(
  coefficients,
  random_effects_mat,
  subject_levels,
  n_re,
  part2,
  price,
  subject_id,
  epsilon,
  id_var
)

Arguments

Value

Data frame of subject-specific parameters.

Generate Default Starting Values for Hurdle Model

Description

Internal function to generate sensible default starting values for TMB hurdle model optimization.

Usage

.hurdle_default_starts(consumption, n_re, n_subjects, has_k)

Arguments

Value

A named list of starting values for TMB parameters.

Extract Results from TMB Hurdle Fit

Description

Internal function to extract coefficients, standard errors, and derived quantities from a fitted TMB hurdle model.

Usage

.hurdle_extract_estimates(opt, obj, sdr, param_names, n_subjects, n_re)

Arguments

Value

A list with coefficients, se, variance_components, correlations, and random effects matrix.

Get Parameter Names for Hurdle Model

Description

Internal function to get the names of fixed effect parameters and variance component names based on model configuration.

Usage

.hurdle_get_param_names(n_re, has_k)

Arguments

Value

A list with fixed_names, var_names, and rho_names.

Normalize User-Provided Starting Values

Description

Internal function to handle backwards-compatibility and normalization of user-provided starting values.

Usage

.hurdle_normalize_starts(start_values, n_re, has_k)

Arguments

Value

Normalized starting values list.

Prepare Hurdle Model Data

Description

Internal function to prepare data structures for TMB hurdle model fitting. Converts subject IDs to 0-indexed integers and creates derived variables.

Usage

.hurdle_prepare_data(data, y_var, x_var, id_var)

Arguments

Value

A list containing:

price

Numeric vector of prices

consumption

Numeric vector of consumption values

delta

Integer vector (1 if consumption == 0, else 0)

logQ

Log consumption (0 for zeros)

subject_id

0-indexed subject IDs for C++

subject_levels

Unique subject IDs (original)

n_subjects

Number of unique subjects

Transform Random Effects to Original Scale

Description

Internal function to transform standardized random effects (u) to the original scale using the Cholesky decomposition of the covariance matrix.

Usage

.hurdle_transform_random_effects(u_hat, coefficients, n_re)

Arguments

Value

Matrix of random effects on original scale with named columns.

Analytic Pmax for HS/Exponential Model (Lambert W)

Description

For the exponential model: Q(p) = Q0 * 10^(k * (exp(-alpha * Q0 * p) - 1)) Pmax = -W_0(-1 / (k * ln(10))) / (alpha * Q0)

Usage

.pmax_analytic_hs(alpha_nat, q0_nat, k_nat)

Arguments

Value

List with pmax, method, and notes

Analytic Pmax for Hurdle Model (Natural-parameter exponential form)

Description

For hurdle: Q(p) = Q0 * exp(k * (exp(-alpha * p) - 1)) Note: No Q0 normalization in exponent (unlike HS) Pmax = -W_0(-1/k) / alpha

Usage

.pmax_analytic_hurdle(alpha_nat, k_nat)

Arguments

Value

List with pmax, method, and notes

Analytic Pmax for HS-Standardized Hurdle Part II (Q0 inside exponent)

Description

For hurdle_hs_stdq0: Q(p) = Q0 * exp(k * (exp(-alpha * Q0 * p) - 1)) Pmax = -W_0(-1/k) / (alpha * Q0)

Usage

.pmax_analytic_hurdle_hs_stdq0(alpha_nat, q0_nat, k_nat)

Arguments

Value

List with pmax, method, and notes

Analytic Pmax for Simplified/SND Model

Description

For SND: Q(p) = Q0 * exp(-alpha * Q0 * p) Pmax = 1 / (alpha * Q0)

Usage

.pmax_analytic_snd(alpha_nat, q0_nat)

Arguments

Value

List with pmax, method, and notes

Numerical Pmax via Optimization

Description

Numerical Pmax via Optimization

Usage

.pmax_numerical(expenditure_fn, price_range)

Arguments

Value

List with pmax, method, boundary info, and notes

Validate and Convert Parameters to Natural Scale

Description

Validate and Convert Parameters to Natural Scale

Usage

.standardize_params_to_natural(params, param_scales = NULL)

Arguments

Value

List with natural-scale parameters and conversion notes

Convert Parameter from Specified Scale to Natural Scale

Description

Convert Parameter from Specified Scale to Natural Scale

Usage

.to_natural_scale(value, scale)

Arguments

Value

Numeric value on natural scale

AIC for Hurdle Demand Model

Description

AIC for Hurdle Demand Model

Usage

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

Arguments

Value

A numeric AIC value.

BIC for Hurdle Demand Model

Description

BIC for Hurdle Demand Model

Usage

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

Arguments

Value

A numeric BIC value.

ChangeData

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 bkaplan.ku@gmail.com

Examples

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

Check Column Names

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 bkaplan.ku@gmail.com

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

Systematic Purchase Task Data Checker

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 bkaplan.ku@gmail.com

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)

ExtraF

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 bkaplan.ku@gmail.com

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

FitCurves

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 bkaplan.ku@gmail.com Shawn Gilroy shawn.gilroy@temple.edu

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)

Fit Pooled/Mean Curves

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 bkaplan.ku@gmail.com

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

Get pmax

Description

...

Usage

GetAnalyticPmax(Alpha, K, Q0)

Arguments

Details

...

Value

Numeric

Author(s)

Shawn Gilroy sgilroy1@lsu.edu

Examples

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

Analytic Pmax Fallback

Description

Fallback method for Analytic Pmax

Usage

GetAnalyticPmaxFallback(K_, A_, Q0_)

Arguments

Details

Derivative-based optimization strategy

Value

numeric

Author(s)

Shawn Gilroy sgilroy1@lsu.edu

Examples

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

Get Purchase Task Descriptive Summary

Description

Calculates descriptive statistics from purchase task data.

Usage

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

Arguments

Details

GetDescriptives() has been superseded by get_descriptive_summary(), which provides a modern S3 interface with standardized methods (print(), summary(), plot()). GetDescriptives() will continue to work but is no longer recommended for new code.

Provides the following descriptive statistics from purchase task data at each price: mean consumption, median consumption, standard deviation of consumption, proportion of 0 values, number of NAs, minimum consumption, and maximum consumption.

Value

Dataframe with descriptive statistics

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

See Also

get_descriptive_summary() for the modern interface

Examples

GetDescriptives(apt)

GetEmpirical

Description

Calculates empirical measures for purchase task data

Usage

GetEmpirical(dat, xcol = "x", ycol = "y", idcol = "id")

Arguments

Details

GetEmpirical() has been superseded by get_empirical_measures(), which provides a modern S3 interface with standardized methods (print(), summary(), plot()). GetEmpirical() will continue to work but is no longer recommended for new code.

Will calculate and return the following empirical measures: Intensity, BP0, BP1, Omax, and Pmax

Value

Data frame of empirical measures

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

See Also

get_empirical_measures() for the modern interface

Examples

## Obtain empirical measures
GetEmpirical(apt)

Get K

Description

Calculates a k value by looking for the max/min consumption across entire dataset and adds .5 to that range

Usage

GetK(dat, mnrange = TRUE)

Arguments

Details

GetK() has been superseded by get_k(), which provides explicit parameters for the adjustment value and optional verbose output for better transparency. GetK() will continue to work but is no longer recommended for new code.

Will look for maximum/minimum greater zero

Value

Numeric

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

See Also

get_k() for the modern interface

Examples

GetK(apt)

Get Shared K

Description

Finds shared k among selected datasets using global regression

Usage

GetSharedK(dat, equation, sharecol = "group")

Arguments

Details

Uses global regression to fit a shared k among datasets. Assumes the dataset is in its final form. Used within FitCurves

Value

Numeric value of shared k

Author(s)

Brent Kaplan bkaplan.ku@gmail.com Shawn P Gilroy shawn.gilroy@temple.edu

Examples

## Find a shared k value across datasets indicated by id

GetSharedK(apt, "hs", sharecol = "id")

Get Values for SimulateDemand

Description

Gets values used in SimulateDemand

Usage

GetValsForSim(dat)

Arguments

Details

Gets values used in SimulateDemand

Value

List of 3: setaparams, sdindex, x

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

Examples

GetValsForSim(apt)

Plot Curve

Description

Creates a single plot object

Usage

PlotCurve(adf, dfrow, newdats, yscale = "log", style = c("modern", "apa"))

Arguments

Details

Creates individual demand curves

Value

ggplot2 graphical object

Author(s)

Shawn Gilroy shawn.gilroy@temple.edu

Examples

## Creates a single plot from elements of an object created by FitCurves
if (interactive()) {
  fc <- FitCurves(apt, "hs", k = 2, detailed = TRUE)
  PlotCurve(fc$adfs[[1]], fc$dfres[1, ], fc$newdats[[1]])
}

Plot Curves

Description

Creates plots

Usage

PlotCurves(dat, outdir = NULL, device = "png", ending = NULL, ask = TRUE, ...)

Arguments

Details

Creates and saves plots of individual demand curves

Value

Nothing

Author(s)

Brent Kaplan bkaplan.ku@gmail.com, Shawn Gilroy shawn.gilroy@temple.edu

Examples

## Interactively view plots from output from FitCurves
if (interactive()) {
  fc <- FitCurves(apt, "hs", k = 2, detailed = TRUE)
  PlotCurves(fc, ask = TRUE)
}

Recode Outliers

Description

Recodes outliers

Usage

RecodeOutliers(df, outval = 3.29, unitshigher = 0)

Arguments

Details

Recodes outliers using Tabachnick and Fidell's (2013) criteria. A variable is standardized and values that are greater/less than or equal to a specified outlier value (specified in standard deviations; default 3.29SD) are recoded to a certain number of units (default 0) higher/lower than the greatest nonoutlier value. Disregards 'NA' values.

Value

Invisibly, a dataframe with original and recoded (if any) values

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

Examples

## If any outliers are detected, they would be coded as 1 unit higher

emp <- GetEmpirical(apt)
RecodeOutliers(emp[, c(2:6)], unitshigher = 1)

Replace Zeros

Description

Replaces 0 values

Usage

ReplaceZeros(dat, nrepl = 1, replnum = 0.01)

Arguments

Details

Replaces specified number of 0s with replacement value.

Value

Dataframe (long form)

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

Examples

## Replace all zeros with .01
ReplaceZeros(apt, nrepl = "all", replnum = .01)

Simulate Demand Data

Description

Simulate demand data

Usage

SimulateDemand(nruns = 10, setparams, sdindex, x, outdir = NULL, fn = NULL)

Arguments

Details

Generates and saves simulated datasets in the manner specified in Koffarnus, Franck, Stein, & Bickel (2015).

Value

Invisibly a list consisting of: rounded consumption values, unrounded consumption values, simulation parameters, and inState and outState of seeds.

Author(s)

Brent Kaplan bkaplan.ku@gmail.com

Examples

## set values
setparams <- vector(length = 4)
setparams <- c(-2.5547, .702521, 1.239893, .320221, 3.096, 1.438231)
names(setparams) <- c("alphalm", "alphalsd", "q0lm", "q0lsd", "k", "yvalssd")
sdindex <- c(2.1978, 1.9243, 1.5804, 1.2465, 0.8104, 0.1751, 0.0380, 0.0270)
x <- c(.1, 1, 3, 10, 30, 100, 300, 1000)
set.seed(1234)
sim <- SimulateDemand(nruns = 1, setparams = setparams, sdindex = sdindex, x = x)
sim

annotation_logticks2

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 shawn.gilroy@temple.edu

ANOVA Method for Hurdle Demand Models

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)

ANOVA Method for NLME Demand Models

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)

Example alcohol purchase task data

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.

Full alcohol purchase task dataset

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)

Augment a beezdemand_fixed Model with Fitted Values and Residuals

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

Augment a beezdemand_hurdle Model with Fitted Values and Residuals

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

Augment a beezdemand_nlme Model with Fitted Values and Residuals

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

Calculate Pmax and Omax with Method Reporting and Parameter-Space Safety

Description

Unified internal engine for pmax/omax computation. Supports analytic solutions (Lambert W for HS/hurdle, closed-form for SND), numerical fallback, and observed (row-wise) metrics. Handles parameter-space conversions transparently.

Usage

beezdemand_calc_pmax_omax(
  model_type = NULL,
  params = NULL,
  param_scales = NULL,
  expenditure_fn = NULL,
  demand_fn = NULL,
  price_obs = NULL,
  consumption_obs = NULL,
  tol = 0.1,
  compute_observed = NULL
)

Arguments

Value

A list with snake_case fields:

pmax_model

Model-based pmax

omax_model

Model-based omax

q_at_pmax_model

Quantity at pmax

method_model

Method used: "analytic_lambert_w", "analytic_snd", "numerical_optimize_observed_domain"

domain_model

Price domain used for computation

is_boundary_model

Logical; is pmax at domain boundary?

elasticity_at_pmax_model

Elasticity evaluated at pmax

unit_elasticity_pass_model

Logical; is elasticity near -1?

note_model

Any notes about model computation

pmax_obs

Observed pmax

omax_obs

Observed omax

method_obs

Method for observed: "row_wise_max"

tie_break_obs

Tie-break rule: "min_price"

n_obs_rows

Number of observation rows

n_unique_prices

Number of unique prices

has_duplicate_prices

Logical; duplicate prices detected?

n_max_ties

Number of rows achieving omax

note_obs

Notes about observed computation

alpha_scale_in

Input scale for alpha

q0_scale_in

Input scale for Q0

k_scale_in

Input scale for k (if applicable)

note_param_space

Notes about parameter conversions

Examples

result <- beezdemand_calc_pmax_omax(
  model_type = "hs",
  params = list(alpha = 0.001, q0 = 10, k = 3),
  price_obs = c(0, 0.5, 1, 2, 4, 8, 16),
  consumption_obs = c(10, 9, 8, 6, 3, 1, 0)
)
result$pmax_model
result$omax_model

Calculate Pmax/Omax for Multiple Subjects

Description

Calculate Pmax/Omax for Multiple Subjects

Usage

beezdemand_calc_pmax_omax_vec(
  params_df,
  model_type,
  param_scales = NULL,
  price_list = NULL,
  consumption_list = NULL,
  ...
)

Arguments

Value

Data frame with pmax/omax results for each subject

Examples

params_df <- data.frame(
  alpha = c(0.001, 0.002),
  q0 = c(10, 15),
  k = c(3, 3)
)
beezdemand_calc_pmax_omax_vec(params_df, model_type = "hs")

S3 Methods for beezdemand_descriptive Objects

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

S3 Methods for beezdemand_empirical Objects

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

Build Fixed-Effects RHS Formula String

Description

Internal helper to construct the right-hand side of a fixed-effects formula from factors, interaction flag, and continuous covariates.

Usage

build_fixed_rhs(
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  data = NULL
)

Arguments

Value

A character string representing the RHS (e.g., "~ 1 + dose + drug").

Calculate Group-Level Demand Metrics

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)

Calculate Observed Pmax/Omax Grouped by ID

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

Calculate Omax and Pmax for Demand Curves

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)

Calculate Omax and Pmax for Multiple Subjects

Description

Vectorized calculation of Omax and Pmax for multiple subjects with individual-specific parameters.

Usage

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

Arguments

Value

A data frame with columns Pmax, Omax, Qmax.

Calculate Amplitude and Persistence

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)

Cannabis/cigarette cross-price responses

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)

Check Demand Model Diagnostics

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)

Check Cross-Price Data for Unsystematic Responding

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)

Check Demand Data for Unsystematic Responding

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)

Check for Unsystematic Patterns in Cross-Price Data

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)

Extract Coefficients from Cross-Price Demand Models

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

Extract Coefficients from Fixed-Effect Demand Fit

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.

Extract Coefficients from Hurdle Demand Model

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.

Extract Coefficients from a beezdemand_nlme Model

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

Collapse Factor Levels for a Specific Parameter

Description

Internal helper to apply level collapsing for a single parameter (Q0 or alpha). Creates new columns with suffix to avoid modifying original factor columns.

Usage

collapse_factor_levels(data, collapse_spec, factors, suffix)

Arguments

Value

A list with:

• data: Modified data frame with new collapsed factor columns

• new_factor_names: Character vector of new factor column names to use

• info: List with original and new levels for each collapsed factor

Compare Nested Hurdle Demand Models

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)

Compare Demand Models

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)

Confidence Intervals for Fixed-Effect Demand Model Parameters

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

Confidence Intervals for Hurdle Demand Model Parameters

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)

Confidence Intervals for Mixed-Effects Demand Model Parameters

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)

Confidence Intervals for Cross-Price NLS Model Parameters

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)

Example cross‐price dataset

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

Run pairwise intercept comparisons for cross-price demand model

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)

Run pairwise slope comparisons for cross-price demand model

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)

Example Experimental Tobacco Marketplace data

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.

Extract All Coefficient Types from Cross-Price Demand Models

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)

Fit a Linear Cross-Price Demand Model

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)

Fit cross-price demand with NLS (+ robust fallbacks)

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

Fit Fixed-Effect Demand Curves

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)

Fit Two-Part Mixed Effects Hurdle Demand Model

Description

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

Usage

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

Arguments

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)

Fit Nonlinear Mixed-Effects Demand Model

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)

Fixed-Effect Demand Curve Fitting

Description

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

Extract Fixed Effects from a beezdemand_nlme Model

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

Extract Fixed Effects from Mixed-Effects Cross-Price Model

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

Format Parameter Name with Scale Prefix

Description

Formats a canonical parameter name with the appropriate scale prefix.

Usage

format_param_name(param, scale = c("natural", "log", "log10"))

Arguments

Value

Character. Formatted parameter name (e.g., "log10_alpha").

Get Canonical Derived Metric Name

Description

Maps legacy derived metric names to canonical names.

Usage

get_canonical_metric(name)

Arguments

Value

Character. Canonical metric name or original if no mapping found.

Get Canonical Parameter Name

Description

Maps legacy or variant parameter names to canonical names.

Usage

get_canonical_param(name)

Arguments

Value

Character. Canonical parameter name or original if no mapping found.

Get Pairwise Comparisons for Demand Parameters

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)

Get Estimated Marginal Means for Demand Parameters

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)

Get Trends (Slopes) of Demand Parameters with respect to Continuous Covariates

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)

Calculate Descriptive Statistics by Price

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)

Calculate Empirical Demand Measures

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

Get Equation Specification

Description

Returns the full specification for an equation ID.

Usage

get_equation_spec(equation_id)

Arguments

Value

List. Equation specification or NULL if not found.

Get Hurdle Model Parameter Summary

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)

Calculate Individual-Level Predicted Coefficients from beezdemand_nlme Model

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)

Calculate K Scaling Parameter for Demand Curve Fitting

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)

Legacy to Canonical Mapping Table

Description

Returns a data frame mapping legacy column names to canonical names.

Usage

get_legacy_mapping()

Value

Data frame with columns: legacy, canonical, type.

Get Estimated Marginal Means for Observed Factor Combinations

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)

Get Starting Values from a Pooled NLS Model (Internal Helper)

Description

Fits a simpler, pooled NLS model (ignoring random effects and fixed effect factors) to derive initial estimates for global Q0 and alpha parameters. These are then used as starting values for the main NLME model intercepts.

Usage

get_pooled_nls_starts(data, y_var, x_var, equation_form)

Arguments

Value

A named list with Q0 and alpha starting values if successful, else NULL.

Get Subject-Specific Parameters

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)

Glance Method for beezdemand_fixed

Description

Glance Method for beezdemand_fixed

Usage

## S3 method for class 'beezdemand_fixed'
glance(x, ...)

Arguments

Value

A one-row tibble of model statistics

Glance at a beezdemand_hurdle Model

Description

Returns a one-row tibble of model-level statistics.

Usage

## S3 method for class 'beezdemand_hurdle'
glance(x, ...)

Arguments

Value

A one-row tibble with columns:

model_class

"beezdemand_hurdle"

backend

"TMB"

nobs

Number of observations

n_subjects

Number of subjects

n_random_effects

Number of random effects

converged

Convergence status

logLik

Log-likelihood

AIC

Akaike Information Criterion

BIC

Bayesian Information Criterion

Glance method for beezdemand_nlme

Description

Glance method for beezdemand_nlme

Usage

## S3 method for class 'beezdemand_nlme'
glance(x, ...)

Arguments

Value

A one-row tibble of model statistics with columns:

• model_class: "beezdemand_nlme"

• backend: "nlme"

• equation_form: The equation form used

• nobs: Number of observations

• n_subjects: Number of subjects

• converged: Convergence status

• logLik, AIC, BIC: Model fit statistics

• sigma: Residual standard error

Glance Method for beezdemand_systematicity

Description

Glance Method for beezdemand_systematicity

Usage

## S3 method for class 'beezdemand_systematicity'
glance(x, ...)

Arguments

Value

A one-row tibble of overall statistics

Get model summaries from a linear cross-price model

Description

Get model summaries from a linear cross-price model

Usage

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

Arguments

Value

A tibble with model summary statistics.

Get model summaries from a mixed-effects cross-price model

Description

Get model summaries from a mixed-effects cross-price model

Usage

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

Arguments

Value

A tibble with model summary statistics.

Get model summaries from a cross-price model

Description

This function extracts model summary statistics from a cross-price demand model into a single-row data frame, following the conventions of the broom package. It returns goodness-of-fit measures and other model information.

Usage

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

Arguments

Value

A one-row data frame with model summary statistics:

Examples

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

Test for significant interaction in a cross-price demand model

Description

Test for significant interaction in a cross-price demand model

Usage

has_significant_interaction(object, alpha = 0.05)

Arguments

Value

Logical indicating whether interaction is significant

Example nonhuman demand data with drug and dose

Description

Data from Ko et al. (2002)

Usage

ko

Format

A tibble with 135 rows and 6 columns:

monkey

unique identifier

x

fixed-ratio requirement

y

consumption

y_ll4

consumption transformed via ll4

drug

drug

dose

dose

Lambert W

Description

Ben Bolker's port of Lambert W from GNU Scientific Library (GPLV3)

Usage

lambertW(z, b = 0, maxiter = 10, eps = .Machine$double.eps, min.imag = 1e-09)

Arguments

Details

Ben Bolker's port of Lambert W from GNU Scientific Library

Value

numeric

Author(s)

Benjamin Bolker (port)

Examples

## Principal branch: W(1) ~ 0.5671
lambertW(1)

## Verify: W(z) * exp(W(z)) == z
w <- lambertW(2)
w * exp(w)

Log-Logistic Transformation (LL4-like)

Description

Applies a log-logistic like transformation, specifically log_base(x^lambda + 1) / lambda. This transformation is useful for compressing data that spans several orders of magnitude while handling zero values gracefully (as x=0 yields 0). It's a variation related to the Box-Cox transformation or a generalized logarithm.

Usage

ll4(x, lambda = 4, base = 10)

Arguments

Value

A numeric vector or scalar of the transformed values. Returns NaN for x < 0 if lambda results in non-real numbers (e.g., even root of negative). However, the intended domain is x >= 0.

Examples

ll4(0)
ll4(1)
ll4(10)
ll4(100)
ll4(c(0, 1, 10, 100, 1000))

# Using a different lambda or base
ll4(10, lambda = 2)
ll4(10, base = exp(1)) # Natural log base

Inverse Log-Logistic Transformation (Inverse LL4-like)

Description

Applies the inverse of the ll4 transformation. Given y = ll4(x), this function calculates x = (base^(y * lambda) - 1)^(1/lambda).

Usage

ll4_inv(y, lambda = 4, base = 10)

Arguments

Value

A numeric vector or scalar of the original, untransformed values. May return NaN if (base^(y * lambda) - 1) is negative and 1/lambda implies an even root (e.g., if lambda is 2 or 4).

Examples

original_values <- c(0, 1, 10, 100, 1000)
transformed_values <- ll4(original_values)
back_transformed_values <- ll4_inv(transformed_values)
print(data.frame(original_values, transformed_values, back_transformed_values))
all.equal(original_values, back_transformed_values) # Should be TRUE or very close

# Example with negative y (log-transformed value)
# If y_ll4 = -0.5 (meaning original value was between 0 and 1 for log10)
ll4_inv(-0.5, lambda = 4, base = 10) # (10^(-0.5*4) - 1)^(1/4) = (0.01 - 1)^(1/4) -> NaN
# The ll4_inv function as provided will return NaN here.
# A more robust version for demand might floor at 0 if NaN occurs.

Extract Log-Likelihood from Hurdle Demand Model

Description

Extracts the log-likelihood from a fitted hurdle demand model. Useful for likelihood ratio tests comparing nested models.

Usage

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

Arguments

Value

An object of class logLik with the log-likelihood value and attributes for degrees of freedom and number of observations.

Examples

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

Low-nicotine cigarette purchase task

Description

Long-form purchase task data across nicotine conditions.

Usage

lowNicClean

Format

A data frame with 5 columns: id, condition, x, y, commodity.

Examples

# data(lowNicClean)

Create a beezdemand_systematicity Object

Description

Create a beezdemand_systematicity Object

Usage

new_beezdemand_systematicity(results, type, call)

Arguments

Experimental Tobacco Marketplace (ETM) data

Description

ETM data across price points with product quantities.

Usage

ongoingETM

Format

A data frame with 6 columns: id, x, AdjCig, FixCig, ECig, flavor.

Examples

# data(ongoingETM)

beezdemand Color Palette

Description

beezdemand Color Palette

Usage

palette_beezdemand()

Value

Named vector of brand colors.

Examples

palette_beezdemand()

Parameter Naming Registry for beezdemand

Description

This file defines canonical parameter naming conventions across all model families. It serves as the single source of truth for parameter names, scales, and mappings.

Naming Convention

Parameters follow snake_case with scale prefixes:

• log10_<param>: parameter in log10 space (e.g., log10_alpha)

• log_<param>: parameter in natural log space (e.g., log_alpha)

• natural_<param>: parameter in natural/linear space (e.g., natural_alpha)

• <param>: canonical name without prefix (context determines scale)

Equation Registry

Equations are grouped by model family:

• NLS (FitCurves): hs, koff, linear. The k parameter defaults to fixed but supports multiple modes: "fit", "ind", "share", "range".

• NLME (fit_demand_mixed): zben, nlme_simplified, nlme_exponentiated. Only nlme_exponentiated uses k (always fixed).

• Hurdle Part II (fit_demand_hurdle): hurdle_zhao, hurdle_hs, hurdle_snd. Parameters are in log space. hurdle_snd has no k parameter. Each variant has 2RE and 3RE sub-models with different variance/correlation parameters.

• Cross-price (fit_cp_nls): cp_exponential, cp_exponentiated, cp_additive. All use log10_qalone, I, log10_beta.

Derived Metrics Naming

• pmax_model: Pmax derived from fitted model parameters

• pmax_obs: Pmax from observed/empirical data (max expenditure price)

• omax_model: Omax derived from fitted model parameters

• omax_obs: Omax from observed/empirical data

• alpha_star: Normalized alpha (Strategy B; Rzeszutek et al., 2025). Makes alpha comparable across different k values by adjusting for the scaling constant. Computed as \alpha^* = -\alpha / \ln(1 - 1/(k \cdot \ln(b))) where b is the logarithmic base (10 for HS/Koff, e for hurdle models). Returned by FitCurves, tidy() on fixed fits, and tidy() on hurdle fits. Standard error (alpha_star_se) is obtained via the delta method. Requires k \cdot \ln(b) > 1; otherwise NA is returned.

Legacy Column Mappings

Legacy FitCurves output uses different naming. The mappings are:

• Pmaxd: derived Pmax from model (approximate formula)

• Pmaxa: analytic Pmax from Lambert W

• Pmaxe: empirical/observed Pmax from data

• Omaxd: derived Omax from model

• Omaxa: analytic Omax from Lambert W

• Omaxe: empirical/observed Omax from data

• Q0d: fitted Q0 (intensity)

• Alpha: fitted alpha (elasticity parameter)

• K: k parameter (range in log units)

Reshape Demand Data Between Wide and Long Formats

Description

Converts demand data between wide format (one row per subject, prices as columns) and long format (one row per observation with id, x, y columns). This is a convenience wrapper around tidyr::pivot_longer() and tidyr::pivot_wider() tailored for behavioral economic purchase task data.

Usage

pivot_demand_data(
  data,
  format = c("long", "wide"),
  id_var = "id",
  x_var = "x",
  y_var = "y",
  price_cols = NULL,
  x_values = NULL,
  drop_na = TRUE
)

Arguments

Details

Wide to Long (format = "long")

The function determines which columns are price columns by:

1. If price_cols is provided, those columns are used directly.

2. If price_cols is NULL, column names are tested with as.numeric(). Columns whose names successfully parse as numbers (e.g., "0", "0.5", "10") are treated as price columns. Remaining non-id_var columns are preserved as extra identifiers.

3. If no column names parse as numeric, all non-id_var columns become price columns and x_values must be supplied.

Actual numeric price values come from x_values if supplied, or from parsing column names. If column names cannot be parsed and x_values is not supplied, an error is raised.

Long to Wide (format = "wide")

Pivots long data so that each unique value in x_var becomes a column, with values from y_var. All columns except x_var and y_var are used as identifiers.

Value

A tibble. When format = "long": columns id, any extra identifier columns, x (numeric price), and y (numeric consumption). When format = "wide": one row per subject with prices as column names.

Examples

# --- Wide to long ---
# Columns named as prices (auto-parsed)
wide_num <- data.frame(
  id = 1:3,
  "0" = c(10, 8, 12), "0.5" = c(9, 7, 11), "1" = c(8, 6, 9),
  check.names = FALSE
)
pivot_demand_data(wide_num, format = "long")

# Columns with non-numeric names require x_values
wide_named <- data.frame(id = 1:2, price_1 = c(10, 8), price_2 = c(5, 4))
pivot_demand_data(wide_named, format = "long",
                  x_values = c(0, 0.5))

# --- Long to wide ---
data(apt, package = "beezdemand")
wide_apt <- pivot_demand_data(apt, format = "wide")
head(wide_apt)

beezdemand Plot Theme and Color Palette

Description

Central plotting style API for beezdemand, with a modern default based on the codedbx "Refined Contemporary" brand and an APA alternative.

Plot Method for beezdemand_fixed

Description

Plot Method for beezdemand_fixed

Usage

## S3 method for class 'beezdemand_fixed'
plot(
  x,
  type = c("demand", "population", "individual", "both"),
  ids = NULL,
  style = c("modern", "apa"),
  show_observed = TRUE,
  show_pred = NULL,
  x_trans = c("log10", "log", "linear", "pseudo_log"),
  y_trans = NULL,
  free_trans = 0.01,
  x_limits = NULL,
  y_limits = NULL,
  n_points = 200,
  x_lab = NULL,
  y_lab = NULL,
  xlab = NULL,
  ylab = NULL,
  facet = NULL,
  observed_point_alpha = 0.5,
  observed_point_size = 1.8,
  pop_line_alpha = 0.9,
  pop_line_size = 1,
  ind_line_alpha = 0.35,
  ind_line_size = 0.7,
  subtitle = NULL,
  ...
)

Arguments

Value

A ggplot2 object.

Plot Demand Curves from Hurdle Demand Model

Description

Creates visualizations of fitted demand curves from a hurdle demand model.

Usage

## S3 method for class 'beezdemand_hurdle'
plot(
  x,
  type = c("demand", "population", "probability", "parameters", "individual", "both"),
  ids = NULL,
  subjects = NULL,
  parameters = c("Q0", "alpha", "breakpoint", "Pmax", "Omax"),
  prices = NULL,
  show_population = TRUE,
  show_pred = NULL,
  show_observed = TRUE,
  x_trans = c("log10", "log", "linear", "pseudo_log"),
  y_trans = NULL,
  free_trans = 0.01,
  facet = NULL,
  x_limits = NULL,
  y_limits = NULL,
  x_lab = NULL,
  y_lab = NULL,
  xlab = NULL,
  ylab = NULL,
  style = c("modern", "apa"),
  observed_point_alpha = 0.5,
  observed_point_size = 1.8,
  pop_line_alpha = 0.9,
  pop_line_size = 1,
  ind_line_alpha = 0.35,
  ind_line_size = 0.7,
  ...
)

Arguments

Value

A ggplot2 object.

Examples

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

# Plot mean demand curve
plot(fit)

# Plot parameter distributions
plot(fit, type = "parameters")

Plot Method for beezdemand_nlme Objects

Description

Creates a ggplot2 visualization of a fitted beezdemand_nlme model, showing observed data points and/or model prediction lines.

Usage

## S3 method for class 'beezdemand_nlme'
plot(
  x,
  type = c("demand", "population", "individual", "both"),
  ids = NULL,
  show_observed = TRUE,
  observed_point_alpha = 0.6,
  show_pred = "population",
  n_points = 200,
  inv_fun = identity,
  facet = NULL,
  at = NULL,
  color_by = NULL,
  linetype_by = NULL,
  shape_by = NULL,
  x_trans = c("log10", "log", "linear", "pseudo_log"),
  y_trans = NULL,
  free_trans = 0.01,
  x_limits = NULL,
  y_limits = NULL,
  style = c("modern", "apa"),
  title = NULL,
  subtitle = NULL,
  x_lab = NULL,
  y_lab = NULL,
  xlab = NULL,
  ylab = NULL,
  observed_point_size = 2,
  pop_line_size = 1,
  ind_line_size = 0.6,
  pop_line_alpha = 0.9,
  ind_line_alpha = 0.3,
  ...
)

Arguments

Value

A ggplot2 object.

Plot Method for Linear Cross-Price Demand Models

Description

Creates a ggplot2 visualization of a fitted linear cross-price demand model (of class cp_model_lm). The plot overlays a prediction line over the observed data points. Axis transformations (e.g., "log10") can be applied. If the model includes group effects, separate lines will be drawn for each group.

Usage

## S3 method for class 'cp_model_lm'
plot(
  x,
  data = NULL,
  inv_fun = identity,
  n_points = 100,
  title = NULL,
  xlab = "Price",
  ylab = "Consumption",
  x_trans = "identity",
  y_trans = "identity",
  point_size = 3,
  ...
)

Arguments

Value

A ggplot2 object displaying the fitted model predictions and observed data.

Plot Method for Mixed-Effects Cross-Price Demand Models

Description

Creates a ggplot2 visualization of a fitted mixed-effects cross-price demand model (of class cp_model_lmer). This function allows you to plot:

Usage

## S3 method for class 'cp_model_lmer'
plot(
  x,
  data = NULL,
  inv_fun = identity,
  n_points = 100,
  title = NULL,
  xlab = "Price",
  ylab = "Consumption",
  x_trans = "identity",
  y_trans = "identity",
  point_size = 3,
  pred_type = c("fixed", "random", "all"),
  ...
)

Arguments

Details

"fixed"

Only the population-level (fixed-effects) prediction.

"random"

Only the subject-specific predictions.

"all"

Both: the fixed-effects and the subject-specific predictions.

If the model includes group effects, separate lines will be drawn for each group.

Value

A ggplot2 object displaying the observed data points along with the prediction curves.

Plot a Cross-Price Demand Model (Nonlinear)

Description

Plot a Cross-Price Demand Model (Nonlinear)

Usage

## S3 method for class 'cp_model_nls'
plot(
  x,
  data = NULL,
  inv_fun = identity,
  n_points = 100,
  title = NULL,
  xlab = "Price",
  ylab = "Consumption",
  x_trans = "identity",
  y_trans = "identity",
  point_size = 3,
  inverse_fun = deprecated(),
  ...
)

Arguments

Value

A ggplot2 object.

Plot Random Effects Q-Q

Description

Creates Q-Q plots for random effects to assess normality assumptions.

Usage

plot_qq(object, which = NULL, ...)

## S3 method for class 'beezdemand_hurdle'
plot_qq(object, which = NULL, ...)

## S3 method for class 'beezdemand_nlme'
plot_qq(object, which = NULL, ...)

Arguments

Value

A ggplot2 object.

Examples

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

Plot Residual Diagnostics

Description

Creates diagnostic plots for model residuals including residuals vs fitted, scale-location, and histogram of residuals.

Usage

plot_residuals(object, type = c("all", "fitted", "histogram", "qq"), ...)

Arguments

Value

A ggplot2 object or list of ggplot2 objects.

Examples

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

Plot Demand Curve for a Single Subject

Description

Creates a demand curve plot for a single subject with optional observed data and population reference curve.

Usage

plot_subject(
  object,
  subject_id,
  prices = NULL,
  show_data = TRUE,
  show_population = TRUE,
  style = c("modern", "apa")
)

Arguments

Value

A ggplot2 object.

Examples

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

Pmax/Omax Engine

Description

Internal engine for calculating pmax, omax, and observed metrics with parameter-space safety and method reporting.

Predict Method for beezdemand_fixed

Description

Predict Method for beezdemand_fixed

Usage

## S3 method for class 'beezdemand_fixed'
predict(
  object,
  newdata = NULL,
  type = c("response", "link"),
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A tibble containing the original newdata columns, plus .fitted and, when requested, .se.fit and .lower/.upper. If newdata does not include an id column, predictions are returned for all subjects (cross product of newdata × subjects) unless k is subject-specific (k = "ind").

Predict Method for Hurdle Demand Models

Description

Returns predictions from a fitted hurdle demand model.

Usage

## S3 method for class 'beezdemand_hurdle'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "parameters", "probability", "demand"),
  prices = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

For type = "parameters", a tibble of subject-level parameters. Otherwise, a tibble containing the newdata columns plus .fitted and helper columns predicted_log_consumption, predicted_consumption, prob_zero, and expected_consumption. When requested, .se.fit and .lower/.upper are included.

Examples

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

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

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

Predict Method for beezdemand_nlme Objects

Description

Generates point predictions from a fitted beezdemand_nlme model. Predictions can be made at the population level (fixed effects only) or group/subject level (fixed + random effects). The output scale depends on the equation_form used during model fitting and whether inv_fun is applied.

Usage

## S3 method for class 'beezdemand_nlme'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "population", "individual"),
  level = 0,
  inv_fun = identity,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  interval_level = 0.95,
  ...
)

Arguments

Value

A tibble containing the original newdata columns plus .fitted. When requested, .se.fit and .lower/.upper are included (currently NA).

See Also

predict.nlme

Examples

data(ko)
fit <- fit_demand_mixed(ko, y_var = "y_ll4", x_var = "x",
                        id_var = "monkey", equation_form = "zben")
# Population-level predictions
preds <- predict(fit, level = 0)

# Subject-level predictions
preds_subj <- predict(fit, level = 1)

Predict method for cp_model_lm objects.

Description

Predict method for cp_model_lm objects.

Usage

## S3 method for class 'cp_model_lm'
predict(object, newdata = NULL, ...)

Arguments

Value

Data frame with predictions.

Predict from a Mixed-Effects Cross-Price Demand Model

Description

Generates predictions from a mixed-effects cross-price demand model (of class cp_model_lmer). The function supports two modes:

Usage

## S3 method for class 'cp_model_lmer'
predict(object, newdata = NULL, pred_type = c("fixed", "random"), ...)

Arguments