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.
R
is an open-source statistical programming language. It is powerful and allows for nearly endless
customizability.
This package is a work in progress. I welcome suggestions, feedback, questions, and comments regarding what other researchers might be interested in seeing. If you encounter bugs, errors, or other discrepancies please either open an issue on the package's GitHub page or contact me and I will do my best to fix the problem.
Right now the package can be obtained from my GitHub page.
There are plans to make it available
properly on CRAN. In any case, to install the package first install
Hadley Wickham's devtools
package:
Then simply install the package using the following command:
```devtools::install_github("brentkaplan/beezdemand", build_vignettes
= TRUE)```
By indicating `build_vignettes = TRUE`, the installation will compile this vignette.
## Using the Package
### Example Dataset
I include an example dataset to demonstrate how data should be
entered and how to use the functions. This example dataset consists
of participant responses on an alcohol purchase task. Participants (id)
reported the number of alcoholic drinks (y) they would be willing to
purchase and consume at various prices (x; USD). Note the long format:
| | id| x| y|
|:--|--:|---:|--:|
|1 | 19| 0.0| 10|
|2 | 19| 0.5| 10|
|3 | 19| 1.0| 10|
|4 | 19| 1.5| 8|
|5 | 19| 2.0| 8|
|6 | 19| 2.5| 8|
|7 | 19| 3.0| 7|
|8 | 19| 4.0| 7|
|17 | 30| 0.0| 3|
|18 | 30| 0.5| 3|
|19 | 30| 1.0| 3|
|20 | 30| 1.5| 3|
|21 | 30| 2.0| 2|
|22 | 30| 2.5| 2|
|23 | 30| 3.0| 2|
|24 | 30| 4.0| 2|
### Converting from Wide to Long and vice versa
Some datasets read into `R` will be in a "wide" format, where column names indicate dataset identifiers:
| x| 19| 30| 38| 60| 68| 106| 113| 142| 156| 188|
|----:|--:|--:|--:|--:|--:|---:|---:|---:|---:|---:|
| 0.0| 10| 3| 4| 10| 10| 5| 6| 8| 7| 5|
| 0.5| 10| 3| 4| 10| 10| 5| 6| 8| 7| 5|
| 1.0| 10| 3| 4| 8| 9| 5| 6| 8| 7| 5|
| 1.5| 8| 3| 4| 8| 9| 5| 6| 6| 7| 5|
| 2.0| 8| 2| 4| 6| 8| 4| 5| 6| 6| 4|
| 2.5| 8| 2| 4| 6| 8| 4| 5| 5| 6| 4|
| 3.0| 7| 2| 4| 5| 7| 4| 5| 5| 5| 4|
| 4.0| 7| 2| 3| 5| 6| 3| 5| 4| 5| 3|
| 5.0| 7| 2| 3| 4| 5| 3| 5| 3| 4| 3|
| 6.0| 6| 2| 3| 4| 5| 2| 5| 3| 3| 2|
| 7.0| 6| 2| 3| 3| 5| 2| 4| 3| 3| 2|
| 8.0| 5| 2| 2| 3| 4| 0| 4| 3| 2| 1|
| 9.0| 5| 1| 2| 2| 4| 0| 4| 3| 2| 1|
| 10.0| 4| 1| 2| 2| 3| 0| 4| 3| 2| 1|
| 15.0| 3| 1| 0| 0| 0| 0| 3| 3| 1| 0|
| 20.0| 2| 1| 0| 0| 0| 0| 2| 3| 0| 0|
The functions in this package primarily deal with data that is in long format, for example the provided dataset `apt` described initially. In order to convert from wide to long formats and vice versa, I recommend using the following commands.
__Wide to Long__
```r
long <- tidyr::gather(wide, id, y, -x)
Long to Wide
wide <- tidyr::spread(long, id, y)
Descriptive values of responses at each price. Includes mean, standard
deviation, proportion of zeros, numer of NAs, and minimum and maximum values. If bwplot = TRUE
, a box-and-whisker plot is also created and saved. Notice the red crosses indicate the mean. User may additionally specify the directory that the plot should save into, the type of file (either "png"
or "pdf"
), and the filename. Defaults are shown here:
GetDescriptives(apt, bwplot = TRUE, outdir = "../plots/", device = "png", filename = "bwplot")
Price | Mean | Median | SD | PropZeros | NAs | Min | Max |
---|---|---|---|---|---|---|---|
0 | 6.8 | 6.5 | 2.62 | 0.0 | 0 | 3 | 10 |
0.5 | 6.8 | 6.5 | 2.62 | 0.0 | 0 | 3 | 10 |
1 | 6.5 | 6.5 | 2.27 | 0.0 | 0 | 3 | 10 |
1.5 | 6.1 | 6.0 | 1.91 | 0.0 | 0 | 3 | 9 |
2 | 5.3 | 5.5 | 1.89 | 0.0 | 0 | 2 | 8 |
2.5 | 5.2 | 5.0 | 1.87 | 0.0 | 0 | 2 | 8 |
3 | 4.8 | 5.0 | 1.48 | 0.0 | 0 | 2 | 7 |
4 | 4.3 | 4.5 | 1.57 | 0.0 | 0 | 2 | 7 |
5 | 3.9 | 3.5 | 1.45 | 0.0 | 0 | 2 | 7 |
6 | 3.5 | 3.0 | 1.43 | 0.0 | 0 | 2 | 6 |
7 | 3.3 | 3.0 | 1.34 | 0.0 | 0 | 2 | 6 |
8 | 2.6 | 2.5 | 1.51 | 0.1 | 0 | 0 | 5 |
9 | 2.4 | 2.0 | 1.58 | 0.1 | 0 | 0 | 5 |
10 | 2.2 | 2.0 | 1.32 | 0.1 | 0 | 0 | 4 |
15 | 1.1 | 0.5 | 1.37 | 0.5 | 0 | 0 | 3 |
20 | 0.8 | 0.0 | 1.14 | 0.6 | 0 | 0 | 3 |
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"
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
ChangeData(apt, nrepl = 1, replnum = 0.01, rem0 = FALSE, remq0e = FALSE, replfree = NULL)
Examine consistency of demand data using Stein et al.'s (2015) alogrithm for identifying unsystematic responses. Default values shown, but they can be customized.
CheckUnsystematic(apt, deltaq = 0.025, bounce = 0.1, reversals = 0, ncons0 = 2)
id | TotalPass | DeltaQ | DeltaQPass | Bounce | BouncePass | Reversals | ReversalsPass | NumPosValues |
---|---|---|---|---|---|---|---|---|
19 | 3 | 0.2112 | Pass | 0 | Pass | 0 | Pass | 16 |
30 | 3 | 0.1437 | Pass | 0 | Pass | 0 | Pass | 16 |
38 | 3 | 0.7885 | Pass | 0 | Pass | 0 | Pass | 14 |
60 | 3 | 0.9089 | Pass | 0 | Pass | 0 | Pass | 14 |
68 | 3 | 0.9089 | Pass | 0 | Pass | 0 | Pass | 14 |
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, and exponentiated model. Soon, I will be including the nonlinear mixed effects model, mixed effects versions of the exponential and exponentiated model, and others.
Empirical measures can be obtained separately on their own:
GetEmpirical(apt)
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
19 | 10 | NA | 20 | 45 | 15 |
30 | 3 | NA | 20 | 20 | 20 |
38 | 4 | 15 | 10 | 21 | 7 |
60 | 10 | 15 | 10 | 24 | 8 |
68 | 10 | 15 | 10 | 36 | 9 |
FitCurves()
has several important arguments that can be passed:
equation
can accept hs
or koff
, two of the contemporary equations proposed by Hursh & Silberberg (2008) and Koffarnus et al. (2015), respectively.
k
by default will be calculated based on the maximum and minimum y values of the entire sample and adding .5. Adding this amount was originally proposed by Steven R. Hursh in an early iteration of a Microsoft Excel spreadsheet used to calculate demand metrics. This adjustment was adopted for two reasons. First, when fitting \(Q_0\) as a derived parameter, the value may exceed the empirically observed intensity value. Thus, a k value calculated based only on the observed range of data may underestimate the full fitted range of the curve. Second, we have found that values of \(\alpha\) (as well as values that rely on \(\alpha\), i.e. approximate \(P_{max}\)) display greater discrepancies when smaller values of k are used compared to larger values of k. Other options include "ind"
, which will calculate k based on individual basis, "fit"
, which will fit k as a free parameter on an individual basis, "share"
, which will fit k as a single shared parameter across all data sets (while fitting individual \(Q_0\) and \(\alpha\)).
agg = NULL
is the default. When agg = "Mean"
, data are fit averaged data disregarding any error. When agg = "Pooled"
, all data are used and clustering within individual is ignored.
detailed = FALSE
is the default. This will output a single dataframe of results, as shown below. When detailed = TRUE
, the output is a 3 element list that includes (1) dataframe of results, (2) list of nonlinear regression model objects, (3) list of dataframes containing predicted x and y values (to be used in subsequent plotting), and (4) list of individual dataframes used in fitting.
lobound
and hibound
can accept named vectors that will be used as lower and upper bounds, respectively during fitting. If k = "fit"
, then it should look as follows: lobound = c("q0" = 0, "k" = 0, "alpha" = 0)
and hibound = c("q0" = 25, "k" = 10, "alpha" = 1)
. If k
is not being fit as a parameter, then only "q0"
and "alpha"
should be used in bounding.
Note: Fitting with an equation 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:
Warning message:
Zeros found in data not compatible with equation! Dropping zeros!
The simplest use of FitCurves()
is shown here, only needing to specify dat
and equation
. All other arguments shown are set to their default values.
FitCurves(dat = apt, equation = "hs", k, agg = NULL, detailed = FALSE, xcol = "x", ycol = "y", idcol = "id", groupcol = NULL, lobound, hibound)
Note that this ouput returns a message (No k value specified. Defaulting to empirical mean range +.5
) and the aforementioned warning (Warning message: Zeros found in data not compatible with equation! Dropping zeros!
). With detailed = FALSE
, the only output is the dataframe of results (broken up to show the different types of results). This example fits the exponential equation proposed by Hursh & Silberberg (2008):
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
19 | 10 | NA | 20 | 45 | 15 |
30 | 3 | NA | 20 | 20 | 20 |
38 | 4 | NA | 10 | 21 | 7 |
60 | 10 | NA | 10 | 24 | 8 |
68 | 10 | NA | 10 | 36 | 9 |
Equation | Q0d | K | Alpha | R2 |
---|---|---|---|---|
hs | 10.475734 | 1.031479 | 0.0046571 | 0.9660008 |
hs | 2.932407 | 1.031479 | 0.0134557 | 0.7922379 |
hs | 4.523155 | 1.031479 | 0.0087935 | 0.8662632 |
hs | 10.492134 | 1.031479 | 0.0102231 | 0.9664814 |
hs | 10.651760 | 1.031479 | 0.0061262 | 0.9699408 |
Q0se | Alphase | N | AbsSS | SdRes | Q0Low | Q0High | AlphaLow | AlphaHigh |
---|---|---|---|---|---|---|---|---|
0.4159581 | 0.0002358 | 16 | 0.0193354 | 0.0371632 | 9.583592 | 11.367875 | 0.0041515 | 0.0051628 |
0.2506946 | 0.0017321 | 16 | 0.0978350 | 0.0835955 | 2.394720 | 3.470093 | 0.0097408 | 0.0171706 |
0.2357693 | 0.0008878 | 14 | 0.0259083 | 0.0464653 | 4.009458 | 5.036852 | 0.0068592 | 0.0107277 |
0.6219725 | 0.0005118 | 14 | 0.0236652 | 0.0444083 | 9.136972 | 11.847296 | 0.0091080 | 0.0113382 |
0.3841063 | 0.0002713 | 14 | 0.0109439 | 0.0301992 | 9.814865 | 11.488656 | 0.0055350 | 0.0067173 |
EV | Omaxd | Pmaxd | Notes |
---|---|---|---|
2.0496979 | 45.49394 | 14.393110 | converged |
0.7094189 | 15.74586 | 17.796221 | converged |
1.0855466 | 24.09418 | 17.654534 | converged |
0.9337418 | 20.72481 | 6.546546 | converged |
1.5581899 | 34.58471 | 10.760891 | converged |
Here, the simplest form is shown specifying another equation, "koff"
. This fits the modified exponential equation proposed by Koffarnus et al. (2015):
FitCurves(apt, "koff")
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
19 | 10 | NA | 20 | 45 | 15 |
30 | 3 | NA | 20 | 20 | 20 |
38 | 4 | 15 | 10 | 21 | 7 |
60 | 10 | 15 | 10 | 24 | 8 |
68 | 10 | 15 | 10 | 36 | 9 |
Equation | Q0d | K | Alpha | R2 |
---|---|---|---|---|
koff | 10.131767 | 1.429419 | 0.0029319 | 0.9668576 |
koff | 2.989613 | 1.429419 | 0.0093716 | 0.8136932 |
koff | 4.607551 | 1.429419 | 0.0070562 | 0.8403625 |
koff | 10.371088 | 1.429419 | 0.0068127 | 0.9659117 |
koff | 10.703627 | 1.429419 | 0.0044361 | 0.9444897 |
Q0se | Alphase | N | AbsSS | SdRes | Q0Low | Q0High | AlphaLow | AlphaHigh |
---|---|---|---|---|---|---|---|---|
0.2438729 | 0.0001663 | 16 | 2.908243 | 0.4557758 | 9.608712 | 10.654822 | 0.0025752 | 0.0032886 |
0.1721284 | 0.0013100 | 16 | 1.490454 | 0.3262837 | 2.620434 | 3.358792 | 0.0065620 | 0.0121812 |
0.3078231 | 0.0010631 | 16 | 4.429941 | 0.5625161 | 3.947336 | 5.267766 | 0.0047761 | 0.0093362 |
0.4069382 | 0.0004577 | 16 | 5.010982 | 0.5982703 | 9.498292 | 11.243884 | 0.0058310 | 0.0077945 |
0.4677467 | 0.0003736 | 16 | 8.350830 | 0.7723263 | 9.700410 | 11.706844 | 0.0036349 | 0.0052373 |
EV | Omaxd | Pmaxd | Notes |
---|---|---|---|
1.9957818 | 46.56622 | 15.140905 | converged |
0.6243741 | 14.56810 | 16.052915 | converged |
0.8292622 | 19.34861 | 13.833934 | converged |
0.8588915 | 20.03993 | 6.365580 | converged |
1.3190322 | 30.77608 | 9.472147 | converged |
FitCurves(apt, "hs", agg = "Mean")
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
mean | 6.8 | NA | 20 | 23.1 | 7 |
Equation | Q0d | K | Alpha | R2 |
---|---|---|---|---|
hs | 7.637437 | 1.429419 | 0.0066817 | 0.9807508 |
Q0se | Alphase | N | AbsSS | SdRes | Q0Low | Q0High | AlphaLow | AlphaHigh |
---|---|---|---|---|---|---|---|---|
0.3258955 | 0.0002218 | 16 | 0.02187 | 0.039524 | 6.93846 | 8.336413 | 0.0062059 | 0.0071574 |
EV | Omaxd | Pmaxd | Notes |
---|---|---|---|
0.8757419 | 20.43309 | 8.813583 | converged |
FitCurves(apt, "hs", agg = "Pooled")
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
pooled | 6.8 | NA | 20 | 23.1 | 7 |
Equation | Q0d | K | Alpha | R2 |
---|---|---|---|---|
hs | 6.592488 | 1.031479 | 0.0085032 | 0.460412 |
Q0se | Alphase | N | AbsSS | SdRes | Q0Low | Q0High | AlphaLow | AlphaHigh |
---|---|---|---|---|---|---|---|---|
0.4260508 | 0.0007125 | 146 | 4.677846 | 0.1802361 | 5.750367 | 7.434609 | 0.0070949 | 0.0099115 |
EV | Omaxd | Pmaxd | Notes |
---|---|---|---|
1.122607 | 24.91675 | 12.52644 | converged |
As mentioned earlier, in the function FitCurves
, 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.
FitCurves(apt, "hs", k = "share")
id | Intensity | BP0 | BP1 | Omaxe | Pmaxe |
---|---|---|---|---|---|
19 | 10 | NA | 20 | 45 | 15 |
30 | 3 | NA | 20 | 20 | 20 |
38 | 4 | NA | 10 | 21 | 7 |
60 | 10 | NA | 10 | 24 | 8 |
68 | 10 | NA | 10 | 36 | 9 |
Equation | Q0d | K | Alpha | R2 |
---|---|---|---|---|
hs | 10.014577 | 3.318335 | 0.0011616 | 0.9820968 |
hs | 2.766313 | 3.318335 | 0.0033331 | 0.7641766 |
hs | 4.485810 | 3.318335 | 0.0024580 | 0.8803145 |
hs | 9.721378 | 3.318335 | 0.0024219 | 0.9705985 |
hs | 10.293139 | 3.318335 | 0.0015879 | 0.9722310 |
Q0se | Alphase | N | AbsSS | SdRes | Q0Low | Q0High | AlphaLow | AlphaHigh |
---|---|---|---|---|---|---|---|---|
0.2429150 | 0.0000308 | 16 | 0.0101816 | 0.0269677 | 9.493576 | 10.535578 | 0.0010955 | 0.0012277 |
0.2192797 | 0.0003739 | 16 | 0.1110490 | 0.0890622 | 2.296005 | 3.236621 | 0.0025312 | 0.0041350 |
0.2074990 | 0.0001963 | 14 | 0.0231862 | 0.0439566 | 4.033709 | 4.937912 | 0.0020302 | 0.0028858 |
0.4371060 | 0.0000778 | 14 | 0.0207584 | 0.0415916 | 8.769006 | 10.673751 | 0.0022523 | 0.0025914 |
0.3179670 | 0.0000523 | 14 | 0.0101100 | 0.0290259 | 9.600348 | 10.985930 | 0.0014740 | 0.0017018 |
EV | Omaxd | Pmaxd | Notes |
---|---|---|---|
1.4241851 | 44.55169 | 13.160535 | converged |
0.4963278 | 15.52624 | 16.603781 | converged |
0.6730390 | 21.05416 | 13.884782 | converged |
0.6830742 | 21.36808 | 6.502491 | converged |
1.0418467 | 32.59129 | 9.366896 | converged |
When one has multiple groups, it may be beneficial to compare whether separate curves are preferred over a single curve. This is accomplished by the Extra Sum-of-Squares F-test. This function (using the argument compare
) will determine whether a single \(\alpha\) or a single \(Q_0\) is better than multiple $\alpha$s or $Q0$s. A single curve will be fit, the residual deviations calculated and those residuals are compared to residuals obtained from multiple curves. A resulting _F statistic will be reporting along with a p value.
Example forthcoming.
ExtraF(apt, "hs")
Plots can be created using the PlotCurves
function. This function takes the output from FitCurves
when the argument from FitCurves
, detailed = TRUE
. The default will be to save figures into a plots folder created one directory above the current working directory. Figures can be saved as either PNG or PDF. If the argument ask = TRUE
, then plots will be shown interactively and not saved (ask = FALSE
is the default). Graphs can automatically be created at both an aggregate and individual level.
To learn more about a function and what arguments it takes, type “?” in front of the function name.
?CheckUnsystematic
CheckUnsystematic package:beezdemand R Documentation
Systematic Purchase Task Data Checker
Description:
Applies Stein, Koffarnus, Snider, Quisenberry, & Bickels (2015)
criteria for identification of nonsystematic purchase task data.
Usage:
CheckUnsystematic(dat, deltaq = 0.025, bounce = 0.1, reversals = 0,
ncons0 = 2)
Arguments:
dat: Dataframe in long form. Colums are id, x, y.
deltaq: Numeric vector of length equal to one. The criterion by which
the relative change in quantity purchased will be compared.
Relative changes in quantity purchased below this criterion
will be flagged. Default value is 0.025.
bounce: Numeric vector of length equal to one. The criterion by which
the number of price-to-price increases in consumption that
exceed 25% of initial consumption at the lowest price,
expressed relative to the total number of price increments,
will be compared. The relative number of price-to-price
increases above this criterion will be flagged. Default value
is 0.10.
reversals: Numeric vector of length equal to one. The criterion by
which the number of reversals from number of consecutive (see
ncons0) 0s will be compared. Number of reversals above this
criterion will be flagged. Default value is 0.
ncons0: Numer of consecutive 0s prior to a positive value is used to
flag for a reversal. Value can be either 1 (relatively more
conservative) or 2 (default; as recommended by Stein et al.,
(2015).
Details:
This function applies the 3 criteria proposed by Stein et al.,
(2015) for identification of nonsystematic purchase task data. The
three criteria include trend (deltaq), bounce, and reversals from
0. Also reports number of positive consumption values.
Value:
Dataframe
Author(s):
Brent Kaplan <bkaplan.ku@gmail.com>
Examples:
## Using all default values
CheckUnsystematic(apt, deltaq = 0.025, bounce = 0.10, reversals = 0, ncons0 = 2)
## Specifying just 1 zero to flag as reversal
CheckUnsystematic(apt, deltaq = 0.025, bounce = 0.10, reversals = 0, ncons0 = 1)
Derek D. Reed, Applied Behavioral Economics Laboratory www.behavioraleconlab.com
Paul E. Johnson, Center for Research Methods and Data Analysis, University of Kansas www.crmda.ku.edu
Michael Amlung, Cognitive Neuroscience of Addictions Laboratory www.cnalab.weebly.com
Peter G. Roma, Institutes for Behavior Resources, Inc. www.ibrinc.org
Steven R. Hursh, Institutes for Behavior Resources, Inc. www.ibrinc.org
Shawn P. Gilroy, GitHub
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.
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.
Hursh, S. R., & Silberberg, A. (2008). Economic demand and essential value. Psychological Review, 115 (1), 186-198.
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.
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