Quantitative Models of Operant Demand

Brent A. Kaplan

Innovation Department, Advocates for Human Potential, Inc.

codedbx.com

Author Note

Brent A. Kaplan http://orcid.org/0000-0002-3758-6776

The code to recreate this chapter is located at github.com/brentkaplan/quant-demand. Please visit the GitHub repository or https://codedbx.com for the most current version of this chapter. The author has no conflicts of interest to declare.

Correspondence concerning this article should be addressed to Brent A. Kaplan, Innovation Department, Advocates for Human Potential, Inc., 490-B Boston Post Road, Sudbury, MA, Email: [email protected]; [email protected]

Abstract

This chapter provides a guide to quantitative methods for behavioral economic demand data. I begin with a primer on regression approaches, providing the reader with intuition behind modeling. I then offer seven practical steps for conducting demand analyses using R Statistical Software, from initial data preparation through model fitting and interpretation. I demonstrate these concepts using published data from an Alcohol Purchase Task, illustrating various modeling approaches, including fit-to-group, two-stage, and mixed-effects methods. To maximize accessibility for implementing these techniques, the chapter is written in Quarto, an open-source technical publishing framework. I provide detailed code examples and visualizations while highlighting helpful tools explicitly created to aid demand curve analyses. I discuss some challenges one might encounter when conducting these analyses and provide recommendations and useful resources for further learning. This chapter serves as a practical tutorial and a theoretical overview, making it valuable for researchers at various skill levels.

Keywords: behavioral economics, demand, nonlinear regression, R Statistical Software, mixed effects models

Quantitative Models of Operant Demand

The field of behavioral economic demand focuses on quantifying the value of a reinforcer or outcome. As we will see, methods for quantifying value are vast and there is no single “correct” approach. This chapter is written assuming you, the reader, have some cursory knowledge of behavioral economic demand methods. If you have picked up this chapter without knowing much about demand, I encourage you to peruse Chapter 2 Introduction to Operant Demand and either Chapter 3 Nonhuman Research Methods/Procedures for Studying Operant Demand or Chapter 4 Human Research Methods for Studying Operant Demand. These chapters will provide the background knowledge necessary to understand how the quantitative approaches fit within the broader backdrop and framework of demand methodology.

The chapter is organized as follows. I will first provide a brief primer on regression methods, including linear and nonlinear regression, and the various approaches to analyzing demand data. This primer should give you some intuitive ideas about what happens when we try to fit a model to the data and a better understanding of the differences between regression approaches. I will then highlight a few of the more popular and frequently used demand models (don’t worry, I will sprinkle in citations to direct further reading). Then, similar to Chapter 10 Quantitative Models of Discounting, I will outline and walk through the steps of analyzing demand data using the R Statistical Software and demonstrate several approaches to modeling demand data using freely available data. Then, I will finish the chapter by discussing some issues and considerations when modeling demand data and providing suggestions for resources.

A quick note about this chapter before moving on: the chapter is written using Quarto (Allaire et al., 2024), an “open-source scientific and technical publishing system built on Pandoc” and the source code and supporting files are available in my repository: https://github.com/brentkaplan/quant-demand. I highly encourage you to look at the source file, especially when navigating through some of the R code (while I expose some of the R code used in this document, there is plenty more that you can not see). This chapter would be much more tedious to write if it were not for the developers of the apaquarto extension (Schneider, n.d.), which allows for a seamless process of creating documents in accordance with American Psychological Association’s 7th Edition style requirements.

Important Terms

Before diving into the numerous quantitative models proposed to describe demand data, there are a handful of terms with which familiarity will be helpful while navigating the chapter. Adapted in part from Koffarnus et al. (2022), Table 1 lists some of these terms and their definitions. These terms will be used throughout the chapter when referencing different aspects of demand and its analysis, and they are helpful generally when working with behavioral economic demand data and analyses.

Table 1 ABOUT HERE

A Note on Regression

Regression is the statistical process of investigating the relationship between two or more variables. Typically, regression aims to understand the effect of some independent variable (in other words, some variable that the experimenter manipulates or some variable that may causally affect the outcome) on a dependent variable (i.e., the outcome variable). Regression can take many forms, but for this chapter, we will distinguish between linear and nonlinear regression and spend most of the focus on the latter.

Linear regression versus nonlinear regression

If you have learned about or conducted any regression modeling, you have probably heard of or performed linear regression. General linear regression models the relationship between the independent and dependent variables and, importantly, assumes this relationship is linear and the residuals (i.e., the difference between what the model predicts and the observed data) follow a conditionally normal distribution. This means that no matter what, one unit increase in the value of the independent variable (e.g., age) will result in the same amount of change in the dependent variable (e.g., salary). The general linear model is simply a specific case of the class of models known as generalized linear models, where the assumptions about the residual structure are not so strict. Nonlinear regression models, on the other hand, can hold the assumption of conditionally normally distributed residuals, but is not necessarily inherent and, potentially more importantly, the distinction between the aforementioned linear regression and nonlinear regression lies in the relation of the predictor variables - they are expressed as a nonlinear combination. As such, nonlinear models can be highly flexible in their ability to bend and curve. Because behavioral economic demand methods attempt to model some consummatory behavior (or the self-report of such behavior) by capturing both the inelastic and elastic relationship with price, nonlinear models are especially useful for describing this pattern. Historically, numerous models have been proposed to describe how consumption decreases as a function of increasing costs, starting as far back as 1988 (Hursh et al., 1988).

Regression Approaches for Demand

I want to highlight three regression approaches commonly used to analyze demand data. I will briefly describe each approach but the reader is encouraged to read Kaplan et al. (2021) for a more in-depth discussion of each.

Fit-to-Group Approach

This first approach can be conducted in two ways. The first way is fitting a model to pre-aggregated averaged data. In this approach, data at each price is first averaged across individuals (or organisms) and then a line is fit to these averaged data. Because all the variability in the data is compressed into a singular point for each price, this approach ignores the intra-price variability and is typically not appropriate for statistical inference. However, this approach may be useful for descriptive, graphical, or theoretical equation testing. The second approach is fitting a curve to all data without aggregation, also known as a pooled approach. Under typical circumstances, the resulting point estimates from the regression will be the same regardless of either of these fit-to-group approaches. However, in the case of this pooled approach, the standard errors are likely to be incorrect because all data points are treated as independent observations when there is dependence within an individual or condition.

Two-Stage Approach

This two-stage approach is quite common for analyzing demand data. In this approach, a demand model is fit to each individual’s data separately. This approach results in parameter estimates for each individual and condition. We have called this the two-stage approach because the modeling serves as the first stage, and the resulting demand metrics are used in the second stage of statistical analysis to make some sort of inference or comparison (e.g., a t-test, correlation). One downside to this approach is that the model effectively overfits to the individual’s data and ignores any information about how other individuals in the sample respond to the same experimental conditions.

Mixed-Effects Modeling

Both of the approaches mentioned above are considered “fixed-effects” only approaches. Fixed effects parameter estimates that are constant or unchanging across different observations (i.e., they are “fixed” across individuals or conditions). In the approaches above, these fixed effects are Q₀ and \alpha; for example, in the fit-to-means approach only one sample level Q₀ and \alpha is estimated so each individual inherits those values. Random effects, on the other hand, account for variability that differs across individuals or experimental sessions. They allow the model to recognize that there may be differences between subjects or conditions that are not captured by the fixed effects. You can think of mixed-effects models broadly as the ability to simultaneously estimate “population-level” effects (akin to the fit-to-group approaches) and “subject-level” effects (akin to the two-stage approach), with the model taking into account all the data provided. As such, mixed-effects models combine both fixed and random effects to capture both the within-subject variability and between-subject variability. This approach allows the researcher to properly account for the nested and crossed nature that is typical of demand data. For example, prices are often nested within an individual and the same individuals typically participate in multiple conditions, which is a classic example of crossed random effects. For a much more thorough and detailed discussion of these regression approaches including the numerous benefits of mixed-effects models, see Kaplan et al. (2021).

Brief Overview of Models of Demand

This section will provide a brief overview of some of the more frequently used demand models, and a new model that appears promising in simplifying demand analyses. There are certainly many other models not discussed here and while they could be afforded an entire chapter, that focus would detract from that of the current chapter on implenting the quantitative techniques of fitting these models. As such, a comprehensive review of all models is beyond the scope of the current chapter. Instead, I direct you toward Koffarnus et al. (2022) for a thorough discussion on the various models, including the advantages and disadvantages of each. For a quick reference to the seminal papers and model formulations, see Table 2.

Table 2 ABOUT HERE

Linear Elasticity

The linear elasticity model (Equation 1), proposed in 1988 by Hursh et al. (1988), jointly describes changes in elasticity as a function of a (change in slope with increases in price) and b (the initial slope of the curve). L is analogous to the intercept, or estimated consumption when price is infinitesimally small. In this model, changes in elasticity are assumed to be linear. Later equations improved upon the primary limitations of this model, where the idea of elasticity is jointly described by two parameters (instead of one) and where the model sometimes predicts increases in consumption at very small prices.

log(Q_{j}) = log(L) + b \cdot (log(P)) - a \cdot P_{j} +\varepsilon_{j},j=1,...,k \tag{1}

Exponential

The exponential model of demand (Equation 2), introduced in 2008 by Hursh and Silberberg (2008), utilizes an exponential decay function to describe the relationship between price and consumption. Contrasted with earlier models, this model describes the rate at which elasticity changes across the demand curve in the free parameter \alpha and the “intercept” or consumption when there is no price is Q₀. In this model, the parameter k governs how high and low the curve is allowed to go; generally, smaller values of k will restrict the range of consumption, while larger values of k will allow for wider ranges of consumption. Due to the logarithmic transformations, especially for the dependent variable (y values), this model cannot incorporate zeros as data unless those zero values are manually transformed into a small positive value. Changing these values is not recommended as the decision of what “small” value to use may quite substantially affect the resulting parameter estimates (Koffarnus et al., 2015; Yu et al., 2014).

log(Q_{j}) = log(Q_{0}) + k(e^{-\alpha \cdot Q_{0} \cdot C_{j}}-1)+\varepsilon_{j},j=1,...,k \tag{2}

Exponentiated

The exponentiated model (Equation 3), introduced in 2015 by Koffarnus et al. (2015), serves as a simple mathematical transformation of Equation 2 (both sides of the model are exponentiated), where all the parameters from Equation 2 have the same interpretations. This model has been used frequently in the literature and is often preferred to the original formulation because no logarithmic transformations of y-values are required - the data can be fit directly.

Q_{j}=Q_{0}\cdot10^{k(e^{-\alpha \cdot Q_{0} \cdot C_{j}}-1)}+\varepsilon_{j},j=1,...,k \tag{3}

Simplified Exponential with Normalized Decay

The simplified exponential with normalized decay (Equation 4) is a modification of Equation 3 that functionally removes the need to include the span parameter k. One benefit of this model is that some of the derived metrics (e.g., essential value) have a closed-form solution compared to other models that require approximations. This is a newly proposed model, and studies may benefit from demonstrating its utility compared to other models (e.g., Equation 2, Equation 3).

Q_{j}=Q_{0}\cdot e^{(-\alpha \cdot Q_{0} \cdot C_{j})}+\varepsilon_{j},j=1,...,k \tag{4}

Tools for Demand Analysis

Several tools and software exist for modeling demand. GraphPad Prism (GraphPad Software, n.d.) is one popular program to conduct demand analysis. While Prism is also a popular software for statistical analysis and data visualization (especially in the fied of biostatistics), it is proprietary and costs money to license and use. Other free software languages you can use to run demand analyses (essentially any program that can perform nonlinear regression) include, but are not limited to, software such as R and Python. This chapter will primarily focus on using the R Statistical Programming Language (R Core Team, 2024) to analyze demand data. R is a free, open-source programming language and environment for statistical computing and graphics and is a popular language for statistical analysis and data visualization. R has “packages,” which are collections of functions that can be used to perform specific tasks. For example, the nlmrt and nlme packages (Nash, 2016; Pinheiro et al., 2023) contain functions that are used to run nonlinear fixed effects and mixed effects models, respectively. For this chapter, we will use collections of packages such as tidyverse (Wickham et al., 2019) and easystats (Lüdecke et al., 2022), as well as individual packages such as nlmrt and nlme. I have developed a specific R package called beezdemand [Behavioral Economic Easy Demand; Kaplan et al. (2019); Kaplan (2023)] that contains many helper functions that can be used to conduct various aspects related to demand analyses. Although I will not rely on the beezdemand package for everything in this chapter, I will use some functions and highlight others that may make analyses easier.

One other tool that I have made recently available is shinybeez (Kaplan & Reed, n.d.), which is an R Shiny application that can be used to conduct behavioral economic demand and discounting analyses. A Shiny application is essentially a web application that can be run locally or on the web and allows the user to conduct analyses without having to write any code. The shinybeez application is (currently) available at https://brentkaplan.shinyapps.io/shinybeez/, and the source code is available at https://github.com/brentkaplan/shinybeez (this repository will always contain up-to-date web URLs). In addition to being completely free, shinybeez relies on beezdemand to perform the core functionality of the demand analyses. This allows transparency and replicability in conducting demand analyses regardless of whether you are using the beezdemand package as part of a broader R workflow or using the shinybeez application. Please consult Kaplan and Reed (n.d.) for more detailed information about shinybeez. Finally, you will find many of the papers I have published and referenced in this chapter, as well as tools that I have created at my website: https://codedbx.com. Visit that website to always find the most up-to-date versions of the tools and code that I create.

General Steps for Demand Analysis

Here, I will outline the general steps that can be used to conduct demand analysis. These steps are not necessarily the only way to conduct demand analysis, but they are a good starting point for those new to these analyses, and I will generally follow these steps in this chapter. An important note for these steps is that I highly recommend using them to outline the decision rules used to conduct the analyses a priori. This is important because it will help you to be more transparent and replicable in your analyses and, importantly, decrease the likelihood of getting lost in the garden of forking paths (Gelman & Loken, 2013).

• Step 0: Analysis Development. This step involves developing a plan for the analysis, including the specific demand model and regression approach to be used, the criteria for model selection, and the criteria for model validation. This step should include laying out the decision rules you will use for each of the following seven steps, including contingency plans for if something does not go as planned (e.g., if a model you choose does not converge, what changes will you make?).

• Step 1: Data Preparation and Initial Examination. This step involves preparing the data for analysis, including checking for missing values, outliers, and data transformations. Practical approaches include creating a table of descriptive statistics and potentially creating visualizations such as scatterplots or histograms.

• Step 2: Use Criteria to Examine Systematicity. This step involves using systematic criteria to examine how data conform to our preexisting expectations of demand functions.

• Step 3: Model Fitting. Select the demand model and use nonlinear regression to plot the demand curve. This step involves analyzing the data according to the chosen regression method specified in Step 0 (i.e., fit-to-group, two-stage, mixed effects models).

• Step 4: Visual Analysis and Interpretation. Plot the demand curve using semi- or double-logarithmic scales to create demand (and/or work) functions visually depicting the model. Overlay model fit lines and data points to visually examine how closely the model accurately represents the data points.

• Step 5: Model Comparison and Refinement. If more than one demand model is applicable, compare their fit indices (e.g., R2 values) to select the most accurate data representations. Adjust parameters (e.g., k) to improve the fit and interpretability of the demand curve. Repeat steps 3 through 5 as needed or based on the a priori decision rules.

• Step 6: Calculate Important Metrics. Calculate metrics from the demand curve, including Q₀ (model-based intensity), elasticity, and model-based O_max and P_max. Other metrics that may not be reliant on fitting the demand curve include observed intensity, observed O_max and P_max, and breakpoint.

• Step 7: Interpret and Report Findings. Explain what key measures mean in the analysis context and broader research question. For example, whether specific reinforcers maintain behavior longer or are more resistant to price increases. Provide a detailed summary of the model parameters, including the goodness-of-fit, and discuss any observed deviations or patterns in the demand curve. Report all decisions made during this seven-step process and, ideally, provide the code and data to facilitate replication.

Step 1. Data Preparation and Initial Examination of Example Demand Data

In this chapter, we will work with a dataset I have previously made publicly available in the GitHub repository: https://github.com/brentkaplan/mixed-effects-demand. The data we will use was used previously by Kaplan et al. (2021), where we introduced concepts about mixed-effects models for behavioral economic demand; however, here, I will provide a more step-by-step approach for conducting demand analyses in general according to the steps outlined before.

The dataset is human self-report response data on the Alcohol Purchase Task (APT) from Kaplan and Reed (2018). This study used “happy hour” alcohol pricing to explore how different promotional price frames (e.g., buy-one-get-one-free versus half-price deals) influence alcohol demand. We found that these frames significantly affected consumption patterns, with the buy-one-get-one-free format generating a higher demand. The data we will use here is a subset of the data from this study from one experimental condition.

Alcohol Purchase Task.

The Alcohol Purchase Task consisted of 17 prices: $0, $0.25, $0.5, $1, $1.5, $2, $2.5, $3, $4, $5, $6, $7, $8, $9, $10, $15, $20 and had a standard vignette:

In the questionnaire that follows we would like you to pretend to purchase and consume alcohol during a 5-hour period. Imagine that you and your friends are at a bar on a weekend night from 9:00 p.m. until 2:00 a.m. to see a band. Imagine that you do not have any obligations the next day (i.e., no work or classes). The following questions ask how many drinks you would purchase at various prices. The available drinks are standard size domestic beers (12 oz.), wine (5 oz.), shots of hard liquor (1.5 oz.), or mixed drinks containing one shot of liquor. Assume that you did not drink alcohol or use drugs before you went to the bar, and that you will not drink or use drugs after leaving the bar. You cannot bring your own alcohol or drugs to the bar. Also, assume that the alcohol you are about to purchase is for your consumption only. In other words, you cannot sell the drinks or give them to anyone else. You also cannot bring the drinks home and you have no other alcohol at home. Everything you buy is, therefore, for your own personal use within the 5 hour period that you are at the bar. Please respond to these questions honestly, as if you were actually in this situation. To verify you understand the pretend scenario, you must correctly answer the next three questions before moving on in the questionnaire.

Then, participants were asked to correctly answer three multiple-choice questions before going forward: (a) “In this pretend scenario, how many hours do you have to consume the drinks?”, (b) “In this pretend scenario, how much did you drink before the bar?”, and (c) “In this pretend scenario, what is the drink special?”

Importance of exploring data before fitting

apt |>
  # select only x and y columns
  select(x, y) |>
  # for each unique x value
  group_by(x) |>
  # calculate the following statistics
  summarise(
    Mean = mean(y),
    `Std. Dev` = sd(y),
    Min = min(y),
    `25%` = quantile(y, .25),
    Median = median(y),
    `75%` = quantile(y, .75),
    Max = max(y)
  ) |>
  ungroup() |>
  mutate(x = paste0("$", x)) |>
  rename(
    "Price" = "x"
  )

When collecting demand data, I recommend visualizing or plotting consumption values as a function of the prices. This is done to examine the range of the data and determine whether specific values are too extreme. Table 3 shows descriptive statistics such as mean, standard deviation, and quantiles at each price. As Table 3 shows, the maximum number of drinks a participant reported was 50, which corresponds to 10 alcholic drinks per hour for five hours. These numbers, and numbers that you see, may seem unrealistic and should be investigated. At minimum, the responses by these participant(s) should be examined and at most removed if there is a justifiable and compelling reason for doing so. These decisions should be set forth in Step 0 when creating the analysis plan. Table A1 shows the descriptive table ouput from the GetDescriptives() function in the beezdemand package¹. This table provides similar information to Table 3 that I manually created.

Table 3 ABOUT HERE

Sometimes, it is helpful to examine the data by plotting the relative distribution of responses for each price. Figure 1 shows such data for the APT. In this figure, the gray boxplots represent the 25th to 75th percentile, with the horizontal lines representing the median values, the red line representing the mean, and the green and red dots representing the minimum and maximum values, respectively. You can see that there is quite a bit of variability in the first couple of prices. Overall, the average and median number of self-reported drinks purchased decreases as the price per drink increases, typical of demand data. Plotting data in this way can help you see big-picture trends in the data to identify extreme values or patterns counter to expectations. However, it is also important to attempt to quantify the extent to which demand data are considered “systematic”.

Figure 1 ABOUT HERE

Step 2: Use Criteria to Examine Systematicity

Stein et al. (2015) proposed a set of three criteria for examining and identifying the extent to which demand data (e.g., responses) are deemed “systematic.” These three criteria include trend, bounce, and reversals from zero. Trend assumes that consumption from the first price to the last price should go down (with a suggested default value of at least 0.025 log-unit reduction per log-unit range in price). Bounce is a measure of price-to-price increases/decreases in consumption, as the general expectation is that consumption should decrease overall and not result in frequent price increases. Stein and colleagues suggest a detection limit of 0.10, meaning the criteria is flagged when increases occur for more than 10% of the price increments. Finally, reversals from zero occur when a non-zero response follows two or more consecutive zero responses. Readers are encouraged to consult Stein et al. (2015) for a more detailed discussion of these criteria. Important to note is that these criteria are for identification and are not to be used automatically or thoughtlessly for exclusion.

Although I recommend researchers screen for systematicity and report these numbers, the decision to retain these response sets in subsequent analyses is in the researcher’s hands. If the researcher decides to remove or exclude some response sets based on these criteria, I recommend analyzing both the full dataset and the restricted dataset and reporting the extent to which excluding those response sets affected the results and interpretations.

# pass the dataframe `apt` to the `CheckSystematic` function from the
# beezdemand package. default values are used for the three criteria.
apt_unsys <- beezdemand::CheckUnsystematic(apt)

Here, I use the checkUnsystematic function from the beezdemand package. This function outputs a data frame with the total number of criteria passed, each criteria’s value and whether that criterion was met, and the number of positive values. The number of positive values is helpful because, in some regression approaches (e.g., two-stage), these datasets cannot be fit. Table 4 shows the output from this function². You can see that the criteria values all differ and that only among this small subset, id = 5 passes only two criteria. Numbers and percentages of the number of response sets that pass the three crtieria, as well as how these respondents are treated, should be reported. For example, Table A2 shows the number and percentage of response sets associated with the number of criteria passed.

Table 4 ABOUT HERE

Step 3: Model Fitting

Approaches for Handling k Values

The k value in Equation 2 and Equation 3 generally reflects a the “range” of the data and functionally constrains the upper and lower bounds of the demand model, such that a higher k value will result in a relatively higher and lower asymptote and a lower k value will result in a relatively lower upper bound and a relatively higher lower bound. There are several different ways of calculating k values. One approach does not rely on solving it as a free parameter and is instead calculated based on the observed range of the data. Then, this value is inserted as a constant into the model during the fitting process. Sometimes, a value of 0.5 is added to the calculated observed range because of the potential that an estimated Q₀ value might exceed the range of the observed data. Such changes should be reported when describing the method of calculating k. Another approach is to treat k as a free parameter in the model process. In the two-stage approach, this can be done by solving for k differently for each participant or by sharing a singular (i.e., global) k value across all participants. The former approach typically results in estimates highly specific to each participant, whereas the latter approach considers the entire range of the data.

# manual k calculation
apt_k <- log10(max(apt$y[apt$y > 0]) / min(apt$y[apt$y > 0])) + 0.5

# beezdemand k calculation specifying using the full range of data
apt_k_beez <- beezdemand::GetK(apt, mnrange = FALSE)

# identical results
apt_k

[1] 2.19897

apt_k_beez

[1] 2.19897

In any case, treating k as a free parameter can cause issues with model convergence (i.e., finding an optimal solution based on a set of reasonable parameters). In my experience, when this occurs, the model estimates a very high k value that is many magnitudes higher than the actual range of the data in logarithmic units. For this chapter, I will use the first approach by calculating the observed range in logarithmic units and then adding 0.5 to that range. This approach is equivalent to the beezdemand function GetK with mnrange = FALSE³. The resulting k value is 2.199. See the code chunk.

A Note on Model Specification

Entering the parameters to be solved into the model can be done in one of two ways. In the first, the parameters are entered directly into the model without any “transformation” such that Q₀ and \alpha are optimized in their natural units (e.g., Q₀ in the units of consumption and \alpha in the inverse units of price). Essentially, these units are on different scales. Small changes in Q₀ and \alpha can lead to disproportionately large changes in the predicted demand curve, especially at extreme values (e.g., very high prices where consumption drops sharply). The other method - reexpressing Q₀ and \alpha as logarithms of Q₀ and \alpha during the fitting process - can help make parameter estimation more stable and less susceptible to outliers or noise in the data. For this chapter, I will express the models and the Q₀ and \alpha parameters in the model as logarithms.

Fit to Group

Two “fit-to-group” approaches will be mentioned here briefly. The first is what we have termed the “fit to means,” and the second is what we have termed “pooling data” (or “fit-to-pooled”).

Fit to Means.

# nlmrt::wrapnls version
apt_averaged <- apt |>
  # for each x value (i.e., price)
  group_by(x) |>
  # calculate the mean of y
  summarise(mn = mean(y)) |>
  # fit the nonlinear model
  nlmrt::wrapnls(
    # specify the model by Koffarnus et al.
    mn ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    # find starting values that appear feasible based on the data
    start = list(q0 = 1, alpha = -1),
    data = _,
    control = list(maxiter = 1000)
  )

# minpack.lm::nlsLM version
apt_averaged_nlsLM <- apt |>
  # for each x value (i.e., price)
  group_by(x) |>
  # calculate the mean of y
  summarise(mn = mean(y)) |>
  # fit the nonlinear model
  minpack.lm::nlsLM(
    mn ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    start = list(q0 = 1, alpha = -1),
    data = _,
    control = list(maxiter = 1000)
  )

# beezdemand::FitCurves version
apt_averaged_beez <- beezdemand::FitCurves(
  dat = apt,
  # specify the equation
  equation = "koff",
  # specify the k value
  k = apt_k,
  # specify how the data should be aggregated
  agg = "Mean"
)

Here, I have conducted the fit-to-means approach in three different ways. The first two use nonlinear regression packages nlmrt and minpack.lm. The functions from these two packages, compared to nls from base R, tend to be more robust in convergence and identify sensible initial values to start the optimization. The fit-to-means approach typically results in the easiest convergence and optimization (i.e., this approach tends to work well and does not frequently throw errors). The third approach I use is the FitCurves function from the beezdemand package. Instead of having to specify an equation, this function allows you to specify the equation to fit the data (currently supports Equation 2 and Equation 3) and the aggregation method (NULL for two-stage, “Mean” for fit-to-means, and “Pooled” for the fit-to-pooled approach).

We can see the \hat{Q_{0}} and \hat{\alpha} (the hats indicated predicted values) estimates are expressed in the logarithmic units. \hat{Q_{0}} = 0.777 and \hat{\alpha} = -2.207. When transformed back onto the parameter’s original scale, these values are 5.977 and 0.006. Notice how these are very similar to the results of the beezdemand package (identical to three decimal places): 5.977 and 0.006.

Figure 2 graphically displays the fit-to-means approach. As the figure shows, the line tracks the data points reasonably well. This is a typical, normal-looking demand curve when fit to averaged data. This approach is a simple way to examine how closely the predicted line tracks the data. The beezdemand package has functions to plot the data after fitting the data using the FitCurves function⁴.

Figure 2 ABOUT HERE

Pooling Data.

# nlmrt::wrapnls version
apt_pooled <- apt |>
  nlmrt::wrapnls(
    y ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    start = list(q0 = 1, alpha = -1),
    data = _,
    control = list(maxiter = 1000)
  )

# minpack.lm::nlsLM version
apt_pooled_nlsLM <- apt |>
  # fit the nonlinear model
  minpack.lm::nlsLM(
    y ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    start = list(q0 = 1, alpha = -1),
    data = _,
    control = list(maxiter = 1000)
  )

# beezdemand::FitCurves version
apt_pooled_beez <- beezdemand::FitCurves(
  dat = apt,
  equation = "koff",
  k = apt_k,
  agg = "Pooled"
)

Like the fit-to-means approach, I fit the data using the fit-to-pooled approach via three different methods. As mentioned earlier, both fit-to-group approaches result in highly similar parameter estimates. The standard errors of the estimates and goodness-of-fit statistics differ between the two approaches. Table 5 shows the differences in the model estimates.

Table 5 ABOUT HERE

Figure 3 graphically displays the fit-to-pooled approach. Although all data points were fit to this model, I only show the data points between 0 and 30 on the y-axis. Furthermore, given the transparency of the points, more points in the x-y coordinate reflect a darker mass of points. As is often the case, the curve is nearly indistinguishable from the fit-to-mean approach.

Figure 3 ABOUT HERE

Two-Stage Approach to Modeling

# fit a single wrapnls function for use with apt data
nls_fit <- function(df, kval) {
  fit <- try(nlmrt::wrapnls(
    y ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    start = list(q0 = 1, alpha = -1),
    data = df,
    control = list(maxiter = 1000)
  ), silent = T)
  return(fit)
}

# fit a single nlsLM function for use with apt data
nlsLM_fit <- function(df, kval) {
  fit <- try(minpack.lm::nlsLM(
    y ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
    start = list(q0 = 1, alpha = -1),
    data = df,
    control = list(maxiter = 1000)
  ), silent = T)
  return(fit)
}

# needs workaround because any datasets that do not converge results
# in a logical NA instead of a tibble
# create two replacement tibbles for any NAs; tidied and glanced
# and use these for the `otherwise` argument in the next code block
tidied_na <- tibble(
  `term` = c("q0", "alpha"),
  `estimate` = c(NA, NA),
  `std.error` = c(NA, NA),
  `statistic` = c(NA, NA),
  `p.value` = c(NA, NA)
)

glanced_na <- tibble(
  `sigma` = NA,
  `isConv` = NA,
  `finTol` = NA,
  `logLik` = NA,
  `AIC` = NA,
  `BIC` = NA,
  `deviance` = NA,
  `df.residual` = NA,
  `nobs` = NA
)

# nlmrt::wrapnls version
apt_twostage <- apt |>
  nest(data = c(x, y)) |>
  group_by(id) |>
  mutate(
    # we fit the model here using the function we defined above
    fit = map(data, nls_fit, kval),
    # broom-specific functions to extract the model estimates and
    tidied = map(fit, possibly(broom::tidy, otherwise = tidied_na)),
    # goodness-of-fit statistics
    glanced = map(fit, possibly(broom::glance, otherwise = glanced_na))
  ) |>
  unnest(cols = c("data")) |>
  filter(x %in% 0) |>
  unnest(cols = "tidied")

# minpack.lm::nlsLM version
apt_twostage_nlsLM <- apt |>
  nest(data = c(x, y)) |>
  group_by(id) |>
  mutate(
    fit = map(data, nlsLM_fit, kval),
    tidied = map(fit, possibly(broom::tidy, otherwise = tidied_na)),
    glanced = map(fit, possibly(broom::glance, otherwise = glanced_na))
  ) |>
  unnest(cols = c("data")) |>
  filter(x %in% 0) |>
  unnest(cols = "tidied")

# beezdemand version
apt_twostage_beez <- beezdemand::FitCurves(
    dat = apt,
    equation = "koff",
    k = apt_k,
  )

Table 6 shows the estimates from the two-stage approach for the first 10 participants in the sample (ids 1-10). You can see that each participant has specific \hat{Q_{0}} and \hat{\alpha} estimates, and for participant id = 5, the estimates result in NA values because this participant has no positive consumption values so the model cannot be fit.

Table 6 ABOUT HERE

Figure 4 shows the two-stage approach to modeling the APT data for the first ten response sets in the sample (ids 1-10). The figure shows how well each curve tracks each response set’s data. The curve closely tracks the data because each model is optimized to that response set. Response set five (id = 5) has no line because there are no positive consumption values.

Figure 4 ABOUT HERE

Mixed-Effects Modeling

apt_nlme <- nlme::nlme(
  # specify the equation
  y ~ 10^(q0) * 10^(apt_k * (exp(-10^(alpha) * 10^(q0) * x) - 1)),
  # add column for the k value constant being used
  data = apt |>
    mutate(apt_k = apt_k),
  # specify fixed effects
  fixed = list(
    q0 ~ 1,
    alpha ~ 1
  ),
  # specify random effects
  random = list(nlme::pdDiag(q0 + alpha ~ 1)),
  # starting values from the fixed effects only model
  start = list(fixed = c(
    coef(apt_averaged)[1],
    coef(apt_averaged)[2]
  )),
  groups = ~id,
  method = "ML",
  verbose = 2,
  # specify control parameters for fitting
  control = list(
    msMaxIter = 5000,
    niterEM = 5000,
    maxIter = 5000,
    pnlsTol = .001,
    tolerance = .01,
    apVar = T,
    minScale = .0000001,
    opt = "optim"
  )
)

Figure 5 shows the results of the mixed-effects modeling approach. This figure shows that the mixed-effect model can predict curves for each response set (the semi-transparent black lines; analogous to the two-stage approach) while also predicting the population fixed effects (the solid red line; analogous to the fit-to-group approach).

Figure 5 ABOUT HERE

Figure 6 shows the results of the mixed-effects modeling approach except for each participant. The black line, again, shows the predicted lines generated from the random effects and the red line, which is identical across plots, shows the group-level effects. Notice how the black lines deviate higher and lower compared to the red line. The mixed-effects model takes into account all the data simultaneously to estimate both individual and group-level effects.

Figure 6 ABOUT HERE

Step 4: Visual Analysis and Interpretation

We have analyzed the models visually and examined how well they describe the data using the three regression approaches. Overall, we have seen that Equation 3 describes the data well in all approaches. Provided our mixed effects regression approach and current model adequately fit the data, we can use the results of this model to conduct any model comparisons (i.e., fit a different model using the mixed effects regression approach) going forward and according to whatever predefined analysis steps are outlined from Step 0.

Step 5: Model Comparison and Refinement

This step involves comparing the model(s), refining the model(s) (e.g., choosing a different k value, using more informative initial values, choosing which response sets to include), and ultimately choosing the model that will be used to calculate key metrics and interpret the results. Therefore, Steps 3 through 5 may be repeated as necessary, and remember, all these decisions should be considered and outlined in Step 0. If decisions are made at this point sans any analysis plan, then decisions could be made based on what you think might be most likely to produce the results you want, not necessarily based on what is most appropriate and justifiable based on the existing literature and research question. For example, I have fit and plotted both Equation 3 and Equation 4 using the different approaches discussed in this chapter in Figure A1, Figure A2, and Figure A3.

Step 6: Calculate Important Metrics

Step 6 involves calculating critical metrics related to the data and model chosen. These metrics typically fall into two buckets: observed metrics (those metrics calculated directly from the data) and derived metrics (those metrics calculated based on model fitting). I will explain the metrics in each section and provide helpful functions from the beezdemand package to help calculate them so they do not have to be calculated manually. For ease, I will demonstrate obtaining metrics based on the two-stage approach (using Equation 3) because often times we will want to use the metrics in some second stage analysis (e.g., t-tests, correlations). Both the beezdemand R package and shinybeez will support automatic calculation of these metrics for mixed-effects models in the future.

Observed Metrics

apt_beez_obs <- beezdemand::GetEmpirical(dat = apt)

Depending on the analysis plan, this step involves generating the important metrics of interest and reporting them in a digestible format, such as a table or a figure. Table 7 shows the observed metrics calculated for the first ten ids using the beezdemand package.

The GetEmpirical() function from the beezdemand package calculates the following observed metrics directly from the data⁵:

• Intensity: Intensity is the consumption (i.e., y value) associated with the lowest price (i.e., x value) in the data. Depending on the task, this could be when price is free or when price is some small value.

• Breakpoint 0 (BP~0): The first and lowest price at which consumption is zero.

• Breakpoint 1 (BP~1): The last and highest price at which consumption is nonzero.

• Omax_e: Represents the maximum expenditure (e.g., effort or cost) an organism is willing to exert to obtain the reinforcer, indicating the peak of “spending” for that commodity. This is calculated by the peak of the expenditure curve by multiplying each x value by its corresponding y value

• Pmax_e: The price (i.e., x value) associated with where Omax_e occurs.

Table 7 ABOUT HERE

Derived Metrics

The FitCurves() function from the beezdemand package calculates the following derived metrics based on the model chosen⁶. Table 8 shows the derived metrics calculated for the first ten ids using the beezdemand package.

• Q0_d: The estimated consumption (i.e., y value) when price is zero. This is analagous to the observed metric of Intensity, but has an associated standard error provided in the output.

• \alpha: The estimated rate of change in elasticity across the entire demand curve.

• Essential Value (EV): Essential value quantifies the relative reinforcing efficacy of a commodity, representing how strongly an organism maintains consumption of a commodity even as price increases or access is restricted, and is inversely related to elasticity. The goal of EV is to be comparable across all commodities (and models).

• Omax_d/a: Similar to Observed Omax, except this is (d)erived/(a)nalytic because it is computed after calculating (d)erived/(a)nalytic Pmax.

• Pmax_d/a: (D)erived/(A)nalytic Pmax is the price at which the slope of the tangent line of the demand curve is equal to -1 (also known as unit elasticity).

Table 8 ABOUT HERE

Step 7: Interpret and Report Findings

The final step is to interpret and report the findings of the demand curve analysis. This step involves taking the observed and/or derived metrics, reporting them either in text or a table, and interpreting them in the context of the research question. This step is crucial for communicating the demand curve analysis findings to others and ensuring that the findings are accurate, reliable, and understandable. This step should also involve reporting the goodness-of-fit metrics used to assess the adequacy of the model reported (e.g., R2, Root mean square error, AIC, BIC) depending on what was specified in Step 0.

Issues and Considerations for Demand Curve Analyses

k

In my experience, one of the most common (and frustrating) issues in conducting demand curve analyses is dealing with the k parameter. Conceptually, we may want always to try and model k as a free parameter to let the model determine the best “range” of the data. However, in practice, this can lead to issues, with the most common problem being that the model attempts to estimate an astronomically high k value. While the model may technically converge, estimates of \alpha tend to be quite different (because in Equation 2 and Equation 3, k is intricately linked with \alpha) than if using a more “reasonable” k value that more closely reflects the range of the data. As I have done here, one straightforward way of dealing with this is to determine a reasonable k value from the observed data and enter that value as a constant into the model. This approach tends to work quite well, and I believe it can be justifiable during the analysis. Another consideration is to test out Equation 4 recently proposed by Rzeszutek et al. (In Press), where k is removed from the model entirely. I encourage the reader to consult Rzeszutek et al. (In Press) for a technical discussion of this model if you want to use it. Both beezdemand and shinybeez will have the functionality to fit this model along with Equation 2 and Equation 3.

Modeling Probabilities

Another methodological approach to assessing “demand” is using probability-based purchase tasks Roma et al. (2016). In these purchase tasks, instead of the respondent responding with a quantity of the outcome being purchased, the respondent responds with a probability of purchasing the outcome (e.g., 0% to 100% likelihood of buying the outcome or engaging in the behavior). On the face of it, these tasks appear to resemble typical quantity-based purchase tasks. However, there are several aspects to consider, not only about the task itself but also about the analysis methods and the interpretations of the model results.

The first consideration is whether probability-based purchase tasks are measuring demand in the conventional sense of what quantity-based purchase tasks attempt to measure. While self-report quantity-based purchase tasks are in and of themselves a slight departure from traditional experiential paradigms (e.g., self-administration), literature has shown adequate correspondence between these types of tasks Strickland et al. (2020). Probability-based purchase tasks, while also being self-report, appear relatively newer, and evidence for adequate correspondence between these types of tasks and conventional quantity-based purchase tasks has not yet been firmly established. Taking a step back, the general concept of demand quantifies the extent to which an organism defends its baseline consumption in the face of increasing costs. Translating this to probability-based purchase tasks, should the interpretation be the extent to which an organism (i.e., human) defends its baseline likelihood of engaging in the behavior? Also, how sensitive are respondents to differences in probabilities? For example, is there a meaningful difference between someone responding with a likelihood of 62% versus 69% or 18% versus 23%? This will depend somewhat on the task, but may inject an artificial degree of accuracy into the data.

The second consideration is when modeling these types of data, what should the interpretation of Q₀ be if it is estimated to exceed 1 (or 100%)? Is it possible for someone to be more than 100% likely to engage in a behavior? This parameter could be constrained at a group or individual level basis leaving only one (i.e., \alpha) or two (i.e., \alpha and k) parameters to be estimated, but this assumption may be too restrictive depending on the research question.

Two considerations for addressing these issues are to reconceptualize the purchase task as a dichotomous choice task (i.e., “yes” or “no” at each price point) and consider using a logistic regression model. When conceptualizing a probability-based purchase task, it may be helpful to consider that, in reality, the respondent can only make one choice: they either engage in the behavior or they do not. Therefore, a dichotomous choice task may more accurately reflect the individual’s actual behavior in the real world. Data from such a task lends itself well to logistic regression models where the outcome is a binary variable. The logistic regression model, therefore, is used to estimate the probability of engaging in the behavior at each price point rather than relying on the respondent to report these probabilities with an unknown degree of accuracy.

Recommended Resources and Readings

This chapter is not exhaustive of the demand curve literature (including demand curve analysis), regression modeling, nor working in R. There are numerous other resources that serve as valuable resources for these topics. With respect to the demand curve literature, I recommend the reader look at papers by some of the following authors (in no specific order): Reed et al. (2013), Strickland et al. (2020), Roma et al. (2016), Reed et al. (2022), Aston and Meshesha (2020), Reed et al. (2020), Zvorsky et al. (2019), Weinsztok et al. (2023), Kiselica et al. (2016), and Kaplan et al. (2018).

With respect to statistics and the intersection of statistics and behavioral economics, I recommend the reader look at resources by some of the following authors: Young (2018), Yu et al. (2014), Kaplan et al. (2021), Ho et al. (2018), McElreath (2018), and James et al. (2023).

With respect to learning more about R, a number of fantastic resources exist online such as Wickham and Grolemund (2016) and Wickham (2019), https://swirlstats.com, and the R journal. Another good resource for new users of R who are interested in conducting demand curve analyses is the paper by Kaplan et al. (2019) and the associated document “Introduction to R and beezdemand” available at: https://github.com/brentkaplan/beezdemand/tree/master/pobs along with Gilroy and Kaplan (2019) (which introduces and explains GitHub). The former document provides beginner-friendly steps for using R and recommends resources for learning its fundamental functionalities.

Conclusion

Quantitative analysis of behavioral economic demand data continues to evolve with new models, methods, and tools becoming available to researchers. In this chapter, I provided a general approach to conducting demand analyses through seven key steps, from initial planning through final interpretation. I demonstrated various regression approaches including fit-to-group, two-stage, and mixed-effects modeling, and their implementation using freely available tools like R and the beezdemand package. By following these systematic approaches and leveraging modern analytical tools, researchers can produce more robust and replicable demand analyses. The field continues to advance with new methodological developments, offering exciting opportunities for future research. As these methods become more accessible through open-source tools, the quality and sophistication of behavioral economic demand analyses will continue to improve. Understanding how these models work is just a first step in advancing the field.

References

Allaire, J. J., Teague, C., Scheidegger, C., Xie, Y., & Dervieux, C. (2024). Quarto. https://doi.org/10.5281/zenodo.5960048

Aston, E. R., & Meshesha, L. Z. (2020). Assessing Cannabis Demand: A Comprehensive Review of the Marijuana Purchase Task. Neurotherapeutics, 17(1), 87–99. https://doi.org/10.1007/s13311-019-00819-z

Gelino, B. W., Kaplan, B. A., & Reed, D. D. (2023). A Behavioral Economic Analysis of Carbon-Neutral Home Energy Enrollment. Behavior and Social Issues, 32(2), 517–533. https://doi.org/10.1007/s42822-023-00143-4

Gelman, A., & Loken, E. (2013). The garden of forking paths: Why multiple comparisons can be a problem, even when there is no “fishing expedition” or “p-hacking” and the research hypothesis was posited ahead of time.

Gilroy, S. P., & Kaplan, B. A. (2019). Furthering Open Science in Behavior Analysis: An Introduction and Tutorial for Using GitHub in Research. Perspectives on Behavior Science, 42(3), 565–581. https://doi.org/10.1007/s40614-019-00202-5

Gilroy, S. P., Kaplan, B. A., Schwartz, L. P., Reed, D. D., & Hursh, S. R. (2021). A zero-bounded model of operant demand. Journal of the Experimental Analysis of Behavior, 115(3), 729–746. https://doi.org/10.1002/jeab.679

GraphPad Software. (n.d.). GraphPad software.

Harsin, J. D., Gelino, B. W., Strickland, J. C., Johnson, M. W., Berry, M. S., & Reed, D. D. (2021). Behavioral economics and safe sex: Examining condom use decisions from a reinforcer pathology framework. Journal of the Experimental Analysis of Behavior, 116(2), 149–165. https://doi.org/10.1002/jeab.706

Ho, Y.-Y., Nhu Vo, T., Chu, H., Luo, X., & Le, C. T. (2018). A Bayesian hierarchical model for demand curve analysis. Statistical Methods in Medical Research, 27(7), 2038–2049. https://doi.org/10.1177/0962280216673675

Hursh, S. R., Raslear, T. G., Shurtleff, D., Bauman, R., & Simmons, L. (1988). A COST-BENEFIT ANALYSIS OF DEMAND FOR FOOD. Journal of the Experimental Analysis of Behavior, 50(3), 419–440. https://doi.org/10.1901/jeab.1988.50-419

Hursh, S. R., & Silberberg, A. (2008). Economic demand and essential value. Psychological Review, 115(1), 186–198. https://doi.org/10.1037/0033-295X.115.1.186

James, G., Witten, D., Hastie, T., Tibshirani, R., & Taylor, J. (2023). An Introduction to Statistical Learning: With Applications in Python. Springer International Publishing. https://doi.org/10.1007/978-3-031-38747-0

Kaplan, B. A. (2023). Beezdemand: Behavioral economic easy demand.

Kaplan, B. A., Foster, R. N. S., Reed, D. D., Amlung, M., Murphy, J. G., & MacKillop, J. (2018). Understanding alcohol motivation using the alcohol purchase task: A methodological systematic review. Drug and Alcohol Dependence, 191, 117–140. https://doi.org/10.1016/j.drugalcdep.2018.06.029

Kaplan, B. A., Franck, C. T., McKee, K., Gilroy, S. P., & Koffarnus, M. N. (2021). Applying mixed-effects modeling to behavioral economic demand: An introduction. Perspectives on Behavior Science, 44(2), 333–358. https://doi.org/10.1007/s40614-021-00299-7

Kaplan, B. A., Gilroy, S. P., Reed, D. D., Koffarnus, M. N., & Hursh, S. R. (2019). The R package beezdemand: Behavioral Economic Easy Demand. Perspectives on Behavior Science, 42(1), 163–180. https://doi.org/10.1007/s40614-018-00187-7

Kaplan, B. A., & Reed, D. D. (n.d.). Shinybeez: A shiny app for behavioral economic easy demand and discounting.

Kaplan, B. A., & Reed, D. D. (2018). Happy hour drink specials in the Alcohol Purchase Task. Experimental and Clinical Psychopharmacology, 26(2), 156–167. https://doi.org/10.1037/pha0000174

Kiselica, A. M., Webber, T. A., & Bornovalova, M. A. (2016). Validity of the alcohol purchase task: A meta-analysis. Addiction (Abingdon, England), 111(5), 806–816. https://doi.org/10.1111/add.13254

Koffarnus, M. N., Franck, C. T., Stein, J. S., & Bickel, W. K. (2015). A modified exponential behavioral economic demand model to better describe consumption data. Experimental and Clinical Psychopharmacology, 23(6), 504–512. https://doi.org/10.1037/pha0000045

Koffarnus, M. N., Kaplan, B. A., Franck, C. T., Rzeszutek, M. J., & Traxler, H. K. (2022). Behavioral economic demand modeling chronology, complexities, and considerations: Much ado about zeros. Behavioural Processes, 199, 104646. https://doi.org/10.1016/j.beproc.2022.104646

Lüdecke, D., Ben-Shachar, M. S., Patil, I., Wiernik, B. M., Bacher, E., Thériault, R., & Makowski, D. (2022). Easystats: Framework for easy statistical modeling, visualization, and reporting. CRAN. https://doi.org/10.32614/CRAN.package.easystats

Martínez-Loredo, V., González-Roz, A., Secades-Villa, R., Fernández-Hermida, J. R., & MacKillop, J. (2021). Concurrent validity of the Alcohol Purchase Task for measuring the reinforcing efficacy of alcohol: An updated systematic review and meta-analysis. Addiction (Abingdon, England), 116(10), 2635–2650. https://doi.org/10.1111/add.15379

McElreath, R. (2018). Statistical Rethinking: A Bayesian Course with Examples in R and Stan (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315372495

Nash, J. C. (2016). Nlmrt: Functions for nonlinear least squares solutions [Manual].

Newman, M., & Ferrario, C. R. (2020). An improved demand curve for analysis of food or drug consumption in behavioral experiments. Psychopharmacology, 237(4), 943–955. https://doi.org/10.1007/s00213-020-05491-2

Pinheiro, J., Bates, D., & R Core Team. (2023). Nlme: Linear and nonlinear mixed effects models [Manual].

R Core Team. (2024). R: A language and environment for statistical computing [Manual]. R Foundation for Statistical Computing.

Reed, D. D., Kaplan, B. A., Becirevic, A., Roma, P. G., & Hursh, S. R. (2016). Toward quantifying the abuse liability of ultraviolet tanning: A behavioral economic approach to tanning addiction. Journal of the Experimental Analysis of Behavior, 106(1), 93–106. https://doi.org/10.1002/jeab.216

Reed, D. D., Naudé, G. P., Salzer, A. R., Peper, M., Monroe-Gulick, A. L., Gelino, B. W., Harsin, J. D., Foster, R. N. S., Nighbor, T. D., Kaplan, B. A., Koffarnus, M. N., & Higgins, S. T. (2020). Behavioral economic measurement of cigarette demand: A descriptive review of published approaches to the cigarette purchase task. Experimental and Clinical Psychopharmacology, 28(6), 688–705. https://doi.org/10.1037/pha0000347

Reed, D. D., Niileksela, C. R., & Kaplan, B. A. (2013). Behavioral economics: A tutorial for behavior analysts in practice. Behavior Analysis in Practice, 6(1), 34–54. https://doi.org/10.1007/BF03391790

Reed, D. D., Strickland, J. C., Gelino, B. W., Hursh, S. R., Jarmolowicz, D. P., Kaplan, B. A., & Amlung, M. (2022). Applied behavioral economics and public health policies: Historical precedence and translational promise. Behavioural Processes, 198, 104640. https://doi.org/10.1016/j.beproc.2022.104640

Roma, P. G., Hursh, S. R., & Hudja, S. (2016). Hypothetical purchase task questionnaires for behavioral economic assessments of value and motivation. Managerial and Decision Economics, 37(4-5), 306–323. https://doi.org/10.1002/mde.2718

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

Schneider, W. J. (n.d.). Apaquarto.

Stein, J. S., Koffarnus, M. N., Snider, S. E., Quisenberry, A. J., & Bickel, W. K. (2015). Identification and management of nonsystematic purchase task data: Toward best practice. Experimental and Clinical Psychopharmacology, 23(5), 377–386. https://doi.org/10.1037/pha0000020

Strickland, J. C., Campbell, E. M., Lile, J. A., & Stoops, W. W. (2020). Utilizing the commodity purchase task to evaluate behavioral economic demand for illicit substances: A review and meta-analysis. Addiction, 115(3), 393–406. https://doi.org/10.1111/add.14792

Weinsztok, S. C., Reed, D. D., & Amlung, M. (2023). Substance-related cross-commodity purchase tasks: A systematic review. Psychology of Addictive Behaviors, 37(1), 72–86. https://doi.org/10.1037/adb0000851

Wickham, H. (2019). Advanced R (Second edition). CRC Press/Taylor & Francis Group.

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686

Wickham, H., & Grolemund, G. (2016). R for data science: Import, tidy, transform, visualize, and model data (1st ed.). O’Reilly Media.

Wilson, A. G., Franck, C. T., Koffarnus, M. N., & Bickel, W. K. (2015). Behavioral economics of cigarette purchase tasks: Within-subject comparison of real, potentially real, and hypothetical cigarettes. Nicotine & Tobacco Research, 18(5), 524–530. https://doi.org/10.1093/ntr/ntv154

Winger, G., Galuska, C. M., Hursh, S. R., & Woods, J. H. (2006). Relative Reinforcing Effects of Cocaine, Remifentanil, and Their Combination in Rhesus Monkeys. Journal of Pharmacology and Experimental Therapeutics, 318(1), 223–229. https://doi.org/10.1124/jpet.105.100461

Young, M. E. (2018). A place for statistics in behavior analysis. Behavior Analysis: Research and Practice, 18(2), 193–202. https://doi.org/10.1037/bar0000099

Yu, J., Liu, L., Collins, R. L., Vincent, P. C., & Epstein, L. H. (2014). Analytical Problems and Suggestions in the Analysis of Behavioral Economic Demand Curves. Multivariate Behavioral Research, 49(2), 178–192. https://doi.org/10.1080/00273171.2013.862491

Zvorsky, I., Nighbor, T. D., Kurti, A. N., DeSarno, M., Naudé, G., Reed, D. D., & Higgins, S. T. (2019). Sensitivity of hypothetical purchase task indices when studying substance use: A systematic literature review. Preventive Medicine, 128, 105789. https://doi.org/10.1016/j.ypmed.2019.105789

Table 1

Important behavioral economic demand terms

Term	Similar Terms	Definition
Cost	Price, Unit Price, Effort	The response effort or monetary expenditure associated with obtaining one unit of a commodity. It represents the price, effort, or number of responses required, or a combination of these factors.
Consumption		The amount of a commodity earned or obtained and then consumed at a specific cost.
Purchasing		The quantity of a commodity obtained or earned, typically equal to consumption but particularly relevant in hypothetical scenarios or tasks where not all purchased items are consumed.
Demand Intensity	Q₀	The level of consumption associated with minimal or zero cost. This may be measured directly (such as free consumption) or derived from demand models.
Demand Elasticity		Reflects the reduction in consumption when cost increases by a single unit. Elasticity often changes with cost, following a curvilinear demand curve. It is generally derived from model parameters, making interpretation dependent on the specific model used and often requiring multiple parameters for precision.
Change in Elasticity	α	As elasticity is cost-dependent and demand is curvilinear, some models measure the rate of elasticity change as cost rises. Like elasticity, this parameter is derived from model-based parameters, which influences interpretation and makes it contextually dependent.
Omax	Max Expenditure	The peak amount spent or maximum effort exerted for any given price. This value can be computed by identifying Pmax and then multiplying it by the consumption level at that price point.
Pmax		The cost associated with the Omax value. Calculated by determining where the demand curve has a slope of -1 in log-log space, irrespective of the specific demand model used.
Breakpoint		The cost threshold that brings consumption down to zero. It can be derived from raw data or model-based equations, though models with exponential decay functions may not reach zero, requiring specialized estimation methods.

Table 2

Quantitative Models of Demand

Citation	Model Formulation
Hursh et al. (1988)	\log(Q) = \log(L) + b(\log(P)) - aP
Winger et al. (2006)	\log(Q) = \log(L)e^{-aP}
Hursh and Silberberg (2008)	\log(Q) = \log(Q_{0}) + k(e^{-\alpha Q_{0}C} - 1)
Yu et al. (2014)	Q = L \cdot P^{b} \cdot e^{-aP}
Koffarnus et al. (2015)	Q = Q_{0} \cdot 10^{k(e^{-\alpha Q_{0}C} - 1)}
Newman and Ferrario (2020)	Q = Q_{0} \left[ 1 + \left( \frac{P}{P_{0}} \right)^b \right]^{-a/b}
Gilroy et al. (2021)	\text{IHS}(Q) = \text{IHS}(Q_{0}) + \text{IHS}(Q_{0}) \cdot e^{-\alpha Q_{0}P} - 1 where \text{IHS}(Q_{0}) = \log_{10} \left(0.5 \cdot Q_{0} + \sqrt{0.25 \cdot Q_{0}^2 + 1} \right)
Rzeszutek et al. (In Press)	Q_{j}=Q_{0}\cdot e^{(-\alpha \cdot Q_{0} \cdot C_{j})}

Note. As noted in previous works, such as Kaplan et al. (2021), many original model formulations did not provide explicit error terms. Error terms describe deviations from a regression line, so they have not been included here.

Table 3

Descriptive statistics for the Alcohol Purchase Task data.

Price	Mean	Std. Dev	Min	25%	Median	75%	Max
$0	5.86	4.71	0	3	5	8	50
$0.25	5.43	4.31	0	3	5	8	40
$0.5	5.26	4.20	0	2	5	8	30
$1	5.05	3.94	0	2	5	7	30
$1.5	4.72	3.71	0	2	4	7	30
$2	4.40	3.42	0	2	4	6	30
$2.5	4.08	3.27	0	2	4	6	30
$3	3.69	3.02	0	1	3	5	30
$4	3.16	2.81	0	1	3	5	30
$5	2.66	2.41	0	1	2	4	20
$6	2.17	2.26	0	0	2	3	25
$7	1.80	2.35	0	0	1	3	44
$8	1.43	1.78	0	0	1	2	19
$9	1.12	1.74	0	0	0	2	30
$10	0.92	1.70	0	0	0	1	30
$15	0.49	1.19	0	0	0	1	25
$20	0.35	0.90	0	0	0	0	16

Table 4

First ten rows from the checkUnsystematic function from the beezdemand package.

id	TotalPass	DeltaQ	DeltaQPass	Bounce	BouncePass	Reversals	ReversalsPass	NumPosValues
1	3	0.30	Pass	0.00	Pass	0	Pass	17
2	3	0.82	Pass	0.00	Pass	0	Pass	12
3	3	0.91	Pass	0.00	Pass	0	Pass	16
4	3	0.84	Pass	0.00	Pass	0	Pass	16
5	2	0.00	Fail	0.00	Pass	0	Pass	0
6	3	0.88	Pass	0.00	Pass	0	Pass	14
7	3	0.82	Pass	0.00	Pass	0	Pass	15
8	3	0.84	Pass	0.00	Pass	0	Pass	15
9	3	0.70	Pass	0.06	Pass	0	Pass	5
10	3	0.88	Pass	0.00	Pass	0	Pass	13

Table 5

Parameter estimates from the fit-to-group approaches for the Alcohol Purchase Task data.

Parameter	Estimate	Std. error	t value	p value	Approach
q0	0.7765	0.0083	93.33	0.0000	Fit-to-Mean
alpha	-2.2069	0.0160	-137.74	0.0000	Fit-to-Mean
q0	0.7765	0.0036	214.14	0.0000	Fit-to-Pooled
alpha	-2.2069	0.0070	-316.03	0.0000	Fit-to-Pooled

Note. Estimates and standard errors are expressed in log(10) units.

Table 6

Two-stage approach to modeling the Alcohol Purchase Task data showing the first 5 participants’ estimates.

id	Parameter	Estimate	Std. error	t value	p value
1	q0	1.0515	0.0173	60.72	0.0000
1	alpha	-2.7218	0.0401	-67.80	0.0000
2	q0	0.7779	0.0367	21.18	0.0000
2	alpha	-2.1152	0.0676	-31.28	0.0000
3	q0	1.0638	0.0253	42.09	0.0000
3	alpha	-2.5226	0.0495	-50.97	0.0000
4	q0	0.8307	0.0208	39.98	0.0000
4	alpha	-2.2015	0.0388	-56.71	0.0000
5	q0	NA	NA	NA	NA
5	alpha	NA	NA	NA	NA

Table 7

Observed metrics for the Alcohol Purchase Task data generated from the beezdemand package. Table shows the first ten response sets.

id	Intensity	BP0	BP1	Omaxe	Pmaxe
1	10	NA	20.00	42.00	7.00
2	5	8.00	7.00	12.50	2.50
3	10	20.00	15.00	30.00	6.00
4	6	20.00	15.00	15.00	15.00
5	0	NA	NA	0.00	0.00
6	8	10.00	9.00	24.00	4.00
7	5	15.00	10.00	12.00	6.00
8	6	15.00	10.00	16.00	8.00
9	2	2.00	1.50	3.00	1.50
10	8	9.00	8.00	20.00	4.00

Table 8

Derived metrics for the Alcohol Purchase Task data generated from the beezdemand package. Table shows the first ten response sets and a subset of the derived metrics reported.

id	Equation	Q0d	K	Alpha	R2	Q0se	Alphase	N	AbsSS	SdRes	EV	Omaxd	Pmaxd	Omaxa	Pmaxa	Notes
1	koff	11.26	2.20	0.00190	0.92	0.45	0.00018	17	13.00	0.93	1.62	42.99	11.95	42.99	11.93	converged
2	koff	6.00	2.20	0.00767	0.86	0.51	0.00119	17	11.56	0.88	0.40	10.64	5.55	10.64	5.54	converged
3	koff	11.58	2.20	0.00300	0.90	0.67	0.00034	17	23.50	1.25	1.02	27.18	7.34	27.18	7.33	converged
4	koff	6.77	2.20	0.00629	0.94	0.32	0.00056	17	4.91	0.57	0.49	12.98	6.00	12.98	5.98	converged
5	koff	-0.00	2.20	0.02000	NA	NA	NA	17	0.00	0.00	0.15	4.08	-89,554,659,598,886.59	4.08	-89,383,368,470,915.56	wrapnls failed to converge, reverted to nlxb
6	koff	9.32	2.20	0.00379	0.89	0.59	0.00047	17	17.90	1.09	0.81	21.53	7.23	21.53	7.22	converged
7	koff	5.53	2.20	0.00834	0.95	0.25	0.00070	17	2.87	0.44	0.37	9.78	5.54	9.78	5.53	converged
8	koff	6.91	2.20	0.00512	0.92	0.36	0.00052	17	6.53	0.66	0.60	15.92	7.21	15.92	7.20	converged
9	koff	3.60	2.20	0.04074	0.74	0.58	0.01139	17	8.47	0.75	0.08	2.00	1.74	2.00	1.74	converged
10	koff	9.52	2.20	0.00462	0.87	0.73	0.00066	17	24.59	1.28	0.66	17.65	5.80	17.65	5.79	converged

Figure 1

Descriptive statistics for the Alcohol Purchase Task data. Boxplots represent the 25% to 75% interquartile range, horizontal lines within the boxplots represent the median, high and low points represent the maximum and minimum values, and the red line represents the mean.

Figure 2

Fit-to-means approach for the Alcohol Purchase Task data.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure 3

Fit-to-pooled approach for the Alcohol Purchase Task data. All data points within 0 - 30 on the y-axis are shown, even though all data are used to fit the model. Data points are transparent so the darker the mass of points, the more data points are located at that price - consumption combination.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure 4

Two-stage approach to modeling the Alcohol Purchase Task data. Figure shows the first ten response sets. The fifth response set (id = 5) does not have a curve because the model cannot be fit to all zeros.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure 5

Predicted values from the mixed-effects model for the Alcohol Purchase Task data. Semi-transparent black lines indicate predicted values (random effects) for each participant while the solid red line indicates the predictions (fixed effects) for the entire sample.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure 6

Mixed-effects modeling approach showing the individual predicted lines from the random effects in black and the predicted line from the population fixed effects in red. Note that the fifth response set (id = 5) now has an individual predicted line because the model leverages information from the entire sample to assign a “best guess” for that response set.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Appendix

Table A1

Descriptive statistics for the Alcohol Purchase Task data generated from the beezdemand package.

Price	Mean	Median	SD	PropZeros	NAs	Min	Max
0	5.86	5	4.71	0.08	0	0	50
0.25	5.43	5	4.31	0.13	0	0	40
0.5	5.26	5	4.20	0.15	0	0	30
1	5.05	5	3.94	0.14	0	0	30
1.5	4.72	4	3.71	0.17	0	0	30
2	4.40	4	3.42	0.17	0	0	30
2.5	4.08	4	3.27	0.19	0	0	30
3	3.69	3	3.02	0.19	0	0	30
4	3.16	3	2.81	0.22	0	0	30
5	2.66	2	2.41	0.25	0	0	20
6	2.17	2	2.26	0.31	0	0	25
7	1.80	1	2.35	0.37	0	0	44
8	1.43	1	1.78	0.43	0	0	19
9	1.12	0	1.74	0.50	0	0	30
10	0.92	0	1.70	0.56	0	0	30
15	0.49	0	1.19	0.71	0	0	25
20	0.35	0	0.90	0.76	0	0	16

Table A2

Breakdown of number of response sets with the number of criteria passed.

TotalPass	n	percentage
0	36	3.27%
1	66	6%
2	52	4.73%
3	946	86%

Figure A1

Model comparison using the fit-to-means approach.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure A2

Model comparison using the fit-to-pooled approach for the Alcohol Purchase Task data. All data points within 0 - 30 on the y-axis are shown, even though all data are used to fit the model. Data points are transparent so the darker the mass of points, the more data points are located at that price - consumption combination.

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Figure A3

Model comparison using the mixed effects model approach

Note. Zero is undefined on the log scale, so this plot uses a pseudo-log scale to visualize zero.

Footnotes

1. shinybeez automatically creates a descriptive table for you when you load your data.

2. shinybeez automatically creates a table of unsystematic criteria similar to Table 4. Criteria values can be specified in the application.

3. shinybeez and beezdemand allow several different options for calculating k values, including treating k as a free parameter.

4. shinybeez automatically creates a plot of the data depending on the regression method used (e.g., fit-to-group or two-stage).

5. shinybeez automatically creates a table of observed metrics for you when you load your data.

6. shinybeez returns a table of derived metrics for you automatically when you run the demand analysis.