---
title: "Using beezdemand"
author: "Brent Kaplan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using beezdemand}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  dpi = 144,
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)
```

# Rationale Behind `beezdemand`

Behavioral economic demand is gaining in popularity. The motivation behind `beezdemand` was to
create an alternative tool to conduct these
analyses. This package is not necessarily meant to be a replacement
for other softwares; rather, it is meant to serve as an
additional tool in the behavioral economist's toolbox. It is meant for
researchers to conduct behavioral economic (be) demand the easy (ez) way.

## Note About Use

`beezdemand` is under active development but aims to be stable for typical applied use.
If you find issues or would like to contribute, please open an issue on my [GitHub page](https://github.com/brentkaplan/beezdemand) or [email me](mailto:bkaplan.ku@gmail.com).
Please check `NEWS.md` for changes in recent versions.

## Installing `beezdemand`

### CRAN Release (recommended method)

The latest stable version of `beezdemand` can be found on [CRAN](https://CRAN.R-project.org/package=beezdemand) and installed using the following command. The first time you install the package, you may be asked to select a CRAN mirror. Simply select the mirror geographically closest to you.

```{r cran-install, eval = FALSE}
install.packages("beezdemand")

library(beezdemand)
```

### GitHub Release

To install a stable release directly from [GitHub](https://github.com/brentkaplan/beezdemand), first install and load the `devtools` package. Then, use `install_github` to install the package and associated vignette. You *don't* need to download anything directly from [GitHub](https://github.com/brentkaplan/beezdemand), as you should use the following instructions:

```{r git-install, eval = FALSE}
install.packages("devtools")

devtools::install_github("brentkaplan/beezdemand", build_vignettes = TRUE)

library(beezdemand)
```

### GitHub Development Version

To install the development version of the package, specify the development branch in `install_github`:

```{r gitdev-install, eval = FALSE}
devtools::install_github("brentkaplan/beezdemand@develop")
```

## Using the Package

```{r packages, include = FALSE, echo = FALSE}
# Package dependencies are specified in DESCRIPTION
# They should already be installed when building vignettes
library(dplyr)
library(tidyr)
library(ggplot2)
library(beezdemand)
```

### Example Dataset

An example dataset of responses on an Alcohol Purchase Task is provided. This object is called `apt` and is located within the `beezdemand` package. These data are a subset from the paper by Kaplan & Reed (2018). Participants (id) reported the number of alcoholic drinks (y) they would be willing to purchase and consume at various prices (x; USD). Note the format of the data, which is called "long format". Long format data are data structured such that repeated observations are stacked in multiple rows, rather than across columns. First, take a look at an extract of the dataset `apt`, where I've subsetted rows 1 through 10 and 17 through 26:

```{r example-data-set, echo=FALSE, results='asis'}
knitr::kable(
    apt[c(1:10, 17:26), ],
    caption = "Example APT (Alcohol Purchase Task) data in long format"
)
```

The first column contains the row number. The second column contains the id number of the series within the dataset. The third column contains the x values (in this specific dataset, price per drink) and the fourth column contains the associated responses (number of alcoholic drinks purchased at each respective price). There are replicates of id because for each series (or participant), several x values were presented.

### Converting from Wide to Long and Vice Versa

Take for example the format of most datasets that would be exported from a data collection software such as Qualtrics or SurveyMonkey or Google Forms:

```{r example-wide}
## the following code takes the apt data, which are in long format, and converts
## to a wide format that might be seen from data collection software
wide <- as.data.frame(tidyr::pivot_wider(apt, names_from = x, values_from = y))
colnames(wide) <- c("id", paste0("price_", seq(1, 16, by = 1)))
knitr::kable(
    wide[1:5, 1:10],
    caption = "Example data in wide format (first 5 participants, first 10 prices)"
)
```

A dataset such as this is referred to as "wide format" because each participant series contains a single row and multiple measurements within the participant are indicated by the columns. This data format is fine for some purposes; however, for `beezdemand`, data are required to be in "long format" (in the same format as the example data described earlier). The `pivot_demand_data()` function makes this conversion easy.

#### Quick conversion with `pivot_demand_data()`

Since our column names (`price_1`, `price_2`, ...) don't encode the actual prices, we supply them via `x_values`:

```{r example-pivot}
long <- pivot_demand_data(
  wide,
  format = "long",
  x_values = c(0, 0.5, 1, 1.50, 2, 2.50, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20)
)
knitr::kable(
    head(long),
    caption = "Wide to long conversion using pivot_demand_data()"
)
```

If your wide data already has numeric column names (e.g., `"0"`, `"0.5"`, `"1"`), `pivot_demand_data()` will auto-detect the prices and no `x_values` argument is needed. See `?pivot_demand_data` for details.

<details>
<summary>Manual approach with tidyr (click to expand)</summary>

For users who want to understand the underlying mechanics, here is the step-by-step approach using `tidyr` directly. First, rename columns to their actual prices:

```{r example-wide-manual}
## make a copy for the manual approach
wide_manual <- wide
newcolnames <- c("id", "0", "0.5", "1", "1.50", "2", "2.50", "3",
                 "4", "5", "6", "7", "8", "9", "10", "15", "20")
colnames(wide_manual) <- newcolnames
```

Then pivot to long format, rename columns, and coerce to numeric:

```{r example-w2l-manual}
long_manual <- tidyr::pivot_longer(wide_manual, -id,
                                   names_to = "price", values_to = "consumption")
long_manual <- dplyr::arrange(long_manual, id)
colnames(long_manual) <- c("id", "x", "y")
long_manual$x <- as.numeric(long_manual$x)
long_manual$y <- as.numeric(long_manual$y)
```

</details>

The dataset is now "tidy" because: (1) each variable forms a column, (2) each observation forms a row, and (3) each type of observational unit forms a table (in this case, our observational unit is the Alcohol Purchase Task data). To learn more about the benefits of tidy data, readers are encouraged to consult Hadley Wikham's essay on [Tidy Data](https://vita.had.co.nz/papers/tidy-data.html).

### Obtain Descriptive Data

Descriptive statistics at each price point can be obtained using
`get_descriptive_summary()`, which returns an S3 object with `print()`,
`summary()`, and `plot()` methods. The legacy `GetDescriptives()` is also
available for backward compatibility.

```{r descriptive}
desc <- get_descriptive_summary(apt)
knitr::kable(
    desc$statistics,
    caption = "Descriptive statistics by price point",
    digits = 2
)
```

The `plot()` method creates a box-and-whisker plot with mean values shown
as red crosses:

```{r descriptive-plot, fig.width=7, fig.height=5}
plot(desc)
```

### Change Data

There are certain instances in which data are to be modified before fitting, for example when using an equation that logarithmically transforms y values. The following function can help with modifying data:

* `nrepl` indicates number of replacement 0 values, either as an integer or `"all"`. If this value is an integer, `n`, then the first `n` 0s will be replaced

* `replnum` indicates the number that should replace 0 values

* `rem0` removes all zeros

* `remq0e` removes y value where x (or price) equals 0

* `replfree` replaces where x (or price) equals 0 with a specified number

```{r change-data, eval = FALSE}
ChangeData(apt, nrepl = 1, replnum = 0.01, rem0 = FALSE, remq0e = FALSE, replfree = NULL)
```

### Identify Unsystematic Responses

Use `check_systematic_demand()` to examine the consistency of purchase task data using Stein et al.'s (2015)
algorithm for identifying unsystematic responses. Default values are shown, but they can be customized.

```{r unsystematic, eval=FALSE}
check_systematic_demand(
  data = apt,
  trend_threshold = 0.025,
  bounce_threshold = 0.1,
  max_reversals = 0,
  consecutive_zeros = 2
)
```

```{r unsystematic-output, echo=FALSE, results='asis'}
unsys <- check_systematic_demand(
  data = apt,
  trend_threshold = 0.025,
  bounce_threshold = 0.1,
  max_reversals = 0,
  consecutive_zeros = 2
)
knitr::kable(
    head(unsys$results, 5),
    caption = "Systematicity check results (first 5 participants)"
)
```

### Analyze Demand Data

Results of the analysis return both empirical and derived measures for
use in additional analyses and model specification. Equations include
the linear model, exponential model, exponentiated model, and
simplified exponential model (Rzeszutek et al., 2025). `beezdemand`
also supports mixed-effects and hurdle demand models (see the dedicated
vignettes for those workflows).

#### Obtaining Empirical Measures

Empirical measures can be obtained separately on their own. The modern
`get_empirical_measures()` returns a tibble with consistent column naming;
the legacy `GetEmpirical()` is also available:

```{r empirical, eval=FALSE}
get_empirical_measures(apt)
```

```{r empirical-output, echo=FALSE, results='asis'}
knitr::kable(
    head(get_empirical_measures(apt)$measures, 5),
    caption = "Empirical demand measures (first 5 participants)",
    digits = 3
)
```

#### Obtaining Derived Measures

Starting with `beezdemand` version `0.2.0`, the recommended interface for fitting
demand curves is `fit_demand_fixed()`. It returns a structured S3 object with
consistent methods like `summary()`, `tidy()`, `glance()`, `predict()`,
`confint()`, `augment()`, and `plot()`.

Key arguments:

* `equation` can be `"linear"`, `"hs"`, `"koff"`, or `"simplified"` (Hursh & Silberberg, 2008; Koffarnus et al., 2015; Rzeszutek et al., 2025).

* `k` can be a fixed numeric value (e.g., `2`) or one of the helper modes:
  `"ind"`, `"fit"`, or `"share"`.

* `agg = NULL` fits per-subject curves. Use `agg = "Mean"` or `agg = "Pooled"`
  for group-level curves.

* `param_space` controls whether optimization happens on the natural scale or
  log10 scale (see `?fit_demand_fixed`).

Note: Fitting with an equation (e.g., `"linear"`, `"hs"`) that doesn't work happily with zero consumption values results in the following. One, a message will appear saying that zeros are incompatible with the equation. Two, because zeros are removed prior to finding empirical (i.e., observed) measures, resulting BP0 values will be all NAs (reflective of the data transformations). The warning message will look as follows:

```{r zero-warning, eval=FALSE}
Warning message:
Zeros found in data not compatible with equation! Dropping zeros!
```

The simplest use of `fit_demand_fixed()` is shown here. This example fits the
exponential equation proposed by Hursh & Silberberg (2008):

```{r hs, eval=FALSE}
fit_demand_fixed(data = apt, equation = "hs", k = 2)
```

```{r hs-setup, include=FALSE}
fit_hs <- fit_demand_fixed(apt, equation = "hs", k = 2)
hs_tidy <- tidy(fit_hs)
hs_glance <- glance(fit_hs)
hs_confint <- confint(fit_hs)
hs_aug <- augment(fit_hs)
hs_diag <- check_demand_model(fit_hs)
```

```{r hs-output, echo=FALSE, results='asis'}
knitr::kable(head(hs_tidy, 10), caption = "Parameter estimates (`tidy()`, first 10 rows)")
knitr::kable(hs_glance, caption = "Model summary (`glance()`)")
knitr::kable(head(hs_confint, 10), caption = "Confidence intervals (`confint()`, first 10 rows)")
knitr::kable(head(hs_aug, 10), caption = "Fitted values and residuals (`augment()`, first 10 rows)")
```

```{r hs-diagnostics}
hs_diag
```

```{r hs-plot, fig.width=7, fig.height=4}
plot(fit_hs)
```

```{r hs-residuals, fig.width=7, fig.height=4}
plot_residuals(fit_hs)$fitted
```

#### Normalized Alpha ($\alpha^*$)

When `k` varies across participants or studies, raw $\alpha$ values are not
directly comparable. The `alpha_star` column in `tidy()` output provides a
normalized version (Strategy B; Rzeszutek et al., 2025) that adjusts for the
scaling constant:

$$\alpha^* = \frac{-\alpha}{\ln\!\left(1 - \frac{1}{k \cdot \ln(b)}\right)}$$

where $b$ is the logarithmic base (10 for HS/Koff equations). Standard errors
are computed via the delta method. `alpha_star` requires $k \cdot \ln(b) > 1$;
otherwise `NA` is returned.

```{r alpha-star, echo=TRUE}
## alpha_star is included in tidy() output for HS and Koff equations
hs_tidy[hs_tidy$term == "alpha_star", c("id", "term", "estimate", "std.error")]
```

Here is the same idea specifying the `"koff"` equation (Koffarnus et al., 2015):

```{r koff, eval=FALSE}
fit_demand_fixed(data = apt, equation = "koff", k = 2)
```

The `"simplified"` equation (Rzeszutek et al., 2025), also known as the Simplified
Exponential with Normalized Decay (SND), handles zeros natively without requiring
data transformation and does not require a `k` parameter:

```{r simplified, eval=FALSE}
fit_demand_fixed(data = apt, equation = "simplified")
```

For a more detailed treatment of the simplified equation, see `vignette("fixed-demand")`.

By specifying `agg = "Mean"`, y values at each x value are aggregated and a single curve is fit to the data (disregarding error around each averaged point):

```{r agg-mean, eval = FALSE}
fit_demand_fixed(data = apt, equation = "hs", k = 2, agg = "Mean")
```

By specifying `agg = "Pooled"`, y values at each x value are aggregated and a single curve is fit to the data and error around each averaged point (but disregarding within-subject clustering):

```{r agg-pooled, eval = FALSE}
fit_demand_fixed(data = apt, equation = "hs", k = 2, agg = "Pooled")
```

### Share k Globally; Fit Other Parameters Locally

As mentioned earlier, in `fit_demand_fixed()`, when `k = "share"` this parameter will be a shared parameter across all datasets (globally) while estimating $Q_0$ and $\alpha$ locally. While this works, it may take some time with larger sample sizes.

```{r share, eval=FALSE}
fit_demand_fixed(data = apt, equation = "hs", k = "share")
```

```{r share-setup, include=FALSE}
fit_share <- fit_demand_fixed(apt, equation = "hs", k = "share")
share_tidy <- tidy(fit_share)
share_glance <- glance(fit_share)
```

```{r share-output, echo=FALSE, results='asis'}
knitr::kable(head(share_tidy, 10), caption = "Parameter estimates (`tidy()`, first 10 rows)")
knitr::kable(share_glance, caption = "Model summary (`glance()`)")
```

### Learn More About Functions

To learn more about a function and what arguments it takes, type "?" in front of the function name.

```{r learn, eval=FALSE}
## Modern interface (recommended)
?check_systematic_demand

## Legacy interface (still available)
?CheckUnsystematic
```

## Next Steps

Now that you are familiar with the basics, explore the other vignettes for more
advanced workflows:

- `vignette("model-selection")` -- Choosing the right model class for your data
- `vignette("fixed-demand")` -- In-depth fixed-effect demand modeling
- `vignette("group-comparisons")` -- Extra sum-of-squares F-test for group comparisons
- `vignette("mixed-demand")` -- Mixed-effects nonlinear demand models (NLME)
- `vignette("mixed-demand-advanced")` -- Advanced mixed-effects topics (factors, EMMs, covariates)
- `vignette("hurdle-demand-models")` -- Two-part hurdle models via TMB
- `vignette("cross-price-models")` -- Cross-price demand analysis
- `vignette("migration-guide")` -- Migrating from `FitCurves()` to `fit_demand_fixed()`

## Free `beezdemand` Alternative with Graphical User Interface

If you are interested in using open source software for analyzing demand curve data but don't know `R`, please check out [shinybeez](https://github.com/brentkaplan/shinybeez), a free open source Shiny app for demand curve analysis (see the companion article [here](https://pubmed.ncbi.nlm.nih.gov/40103003/)).

## Acknowledgments

- Shawn P. Gilroy, Contributor [GitHub](https://github.com/miyamot0)

- Derek D. Reed, Applied Behavioral Economics Laboratory

- Mikhail N. Koffarnus, Addiction Recovery Research Center

- Steven R. Hursh, Institutes for Behavior Resources, Inc.

- Paul E. Johnson, Center for Research Methods and Data Analysis, University of Kansas

- Peter G. Roma, Institutes for Behavior Resources, Inc.

- W. Brady DeHart, Addiction Recovery Research Center

- Michael Amlung, Cognitive Neuroscience of Addictions Laboratory


## Recommended Readings

- 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., Kaplan, B. A., & Becirevic, A. (2015). Basic research on the behavioral economics of reinforcer value. In *Autism Service Delivery* (pp. 279-306). Springer New York. https://doi.org/10.1007/978-1-4939-2656-5_10

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

- 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://dx.doi.org/10.1037/pha0000045

- 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://dx.doi.org/10.1037/pha0000020

- 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

- 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

- 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

- Reed, D. D., Kaplan, B. A., & Gilroy, S. P. (2025). *Handbook of Operant Behavioral Economics: Demand, Discounting, Methods, and Applications* (1st ed.). Academic Press. https://shop.elsevier.com/books/handbook-of-operant-behavioral-economics/reed/978-0-323-95745-8

- Kaplan, B. A. (2025). Quantitative models of operant demand. In D. D. Reed, B. A. Kaplan, & S. P. Gilroy (Eds.), *Handbook of Operant Behavioral Economics: Demand, Discounting, Methods, and Applications* (1st ed.). Academic Press. https://shop.elsevier.com/books/handbook-of-operant-behavioral-economics/reed/978-0-323-95745-8

- Kaplan, B. A., & Reed, D. D. (2025). shinybeez: A Shiny app for behavioral economic easy demand and discounting. *Journal of the Experimental Analysis of Behavior*. https://doi.org/10.1002/jeab.70000

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

- Rzeszutek, M. J., Regnier, S. D., Kaplan, B. A., Traxler, H. K., Stein, J. S., Tomlinson, D., & Koffarnus, M. N. (2025). Identification and management of nonsystematic cross-commodity data: Toward best practice. *Experimental and Clinical Psychopharmacology*. In press.