Behavior Research Methods, Instruments, & Computers

2004, 36 (1), 41-51

Delay discountingdescribes an individuals preference
for a smaller immediate reward over a delayed but larger
reward because the later reward is perceived as less valu-
able, due to the passage of time (Bickel, Madden, &
Petry, 1998). Delay-discountingmeasures have been used
as indices of impulsivity, on the basis of the premise that
impulsiveindividualsprefer more immediate gratification
(Crean, de Wit, & Richards, 2000). A higher level of dis-
counting has been shown to be associated with a greater
liability to substance use and abuse (Bickel, Odum, &
Madden, 1999; Madden, Petry, Badger, & Bickel, 1997;
Petry, 2001a, 2001b).
Traditional assessment of discounting involves fitting
summary models to data concerning indifference with
respect to the choice between immediate and delayed re-
wards (Crean et al., 2000). Indifference experiments
offer choices of immediate rewards (value V
I
) versus de-
layed rewards (value V
D
) and ascertain the indifference
point at which V
I
is perceived to be equivalent to V
D
(in-
difference is denoted V
I
*, and the reward ratio R 5V
I
/V
D
is denoted R* at the indifference point). These preference
questions are of the form, Would you prefer to receive
$10 today or $25 in 7 days? (Here, V
D
5 25, and V
I
5
10.) A series of these indifference points is evaluated
across a range of delay times, and nonlinear regression is
used to fit a function describing the decline in perceived
value as a function of time.
The summary model used most commonly is a hyper-
bolic decline in value with increased delay time (Johnson
& Bickel, 2002; Kirby, 1997; Kirby & Marakovi c, 1996;
Kirby, Petry, & Bickel, 1999; Madden, Bickel, & Jacobs,
1999), which has the form
(1)
This allows us to summarize a subjects tendency to dis-
count as a single value, k, the scale parameter for the
hyperbolic function of delay time. The possible range of
k is 0 to `, but the practical range falls between 0 and 0.5.
Small values of k (10
24
, 10
23
, . . .) indicate a lack of dis-
counting and a general preference for delayed rewards,
whereas larger values (0.11) indicate heavy discounting
and a preference for immediate rewards. The response
does seem to vary with the magnitude of V
D
(Green, My-
erson, & McFadden, 1997; Johnson & Bickel, 2002;
Kirby, 1997; Kirby & Marakovi c, 1995), so separate k val-
ues are often assessed for small, medium, and large V
D
.
R
V
V kt
*
*
. = =
+
I
D
1
1
41 Copyright 2004 Psychonomic Society, Inc.
This work was supportedby Transdisciplinary Tobacco Use Research
Center Grant P50 CA/DA84718to C.L. The authors extend their thanks
to K. P. Burnham and J. L. Norris for discussions on assessing confi-
dence intervals and to K. Kirby, W. Bilker, and three anonymous re-
viewers for reviews of an earlier version of this article. Correspondence
concerning this article should be sent to E. P. Wileyto, Tobacco Use Re-
search Center, University of Pennsylvania, 3535 Market Street, Suite
4100, Philadelphia, PA 19104 (e-mail: [email protected]).
Using logistic regression to estimate
delay-discounting functions
E. PAUL WILEYTO and JANET AUDRAIN-MCGOVERN
University of Pennsylvania, Philadelphia, Pennsylvania
LEONARD H. EPSTEIN
State University of New York, Buffalo, New York
and
CARYN LERMAN
University of Pennsylvania, Philadelphia, Pennsylvania
The monetary choice questionnaire (MCQ) and similar computer tasks ask preference questions in
order to ascertainindifference, the perceivedequivalenceof immediate versus larger delayedrewards.
Indifferencedata are thenfittedwith a hyperbolic function, summarizingthe decline in perceivedvalue
with delay time. We present a fitting method that estimates the hyperbolic parameter k directly from
surveyresponses. Binary preferencesare modeled as a function of time (X
2
) and a transformedreward
ratio(X
1
), yielding logistic regressioncoefficientsb
2
and b
1
. The hyperbolic parameter emerges as

k 5
b
2
/b
1
, where the logistic predicted p 5 .5 (the definition of indifference). The MCQ was administered
to 1,073 adolescents and was scored using both standard and logistic methods. The means for ln(

k)
were similar (standard, 24.53; logistic, 24.51), and the results were highly correlated(r 5.973). Sim-
ulated MCQ data showed that

k was unbiased, except where b
1
$ 21, indicating a vague survey re-
sponse. Jackknife standard errors provided excellent coverage.
42 WILEYTO, AUDRAIN-MCGOVERN, EPSTEIN, AND LERMAN
Although fitting a model to indifference points is
straightforward, determining the indifference points is
cumbersome. To keep the number of questions down to
a reasonable level, computer task assessments given in
the laboratory adjust V
I
depending on prior answers
(Crean et al., 2000; Richards, Mitchell, de Wit, & Sei-
den, 1997), yet even this quasi-optimal sampling still
leaves us with an indeterminate number of questions,
usually well over 100 (Lerman, unpublished data;
1
Richards et al., 1997).
Moving from the laboratory venue to more routine as-
sessment requires a shorter uniform series of questions.
Kirby et al. (1999) introduced the 27-item monetary
choice questionnaire (MCQ) as such an alternative to the
computer task (requiring only nine questions devoted to
each of small, medium, and large delayed rewards). The
survey is scored by inspection, assigning a k value for
each block of nine questions by looking for a single in-
difference point at which the subjects behavior switched
from preferring the delayed reward to preferring the im-
mediate reward. Each block is assigned the k value most
consistent with the array of realized choices and the geo-
metric mean of possible values, if more than one k value
turned out to be most consistent with the data block. Kir-
bys (2000) method reduces the scoring of the MCQ to
entering survey data into a spreadsheet and matching be-
havioral breakpoints in a discrete lookup table.
All of the approaches described thus far depend on
treating indifference points as the random variables of
interest, an approach that is ideal for determining what
type of function to use in summarizing discounting, in
addition to estimating the parameters. Although Kirbys
(2000) method is crude, it does highlight the fact that it
takes fewer data to estimate the parameters when it has
been decided ahead of time that the function is hyper-
bolic than it does to determine what type of function fits
best. This focus on assessing indifference also does not
allowus to simulate survey data; thus, the approach can-
not provide the usual measures of bias that accompany
standard estimation procedures. In addition, the standard
MCQ scoring procedure does not yield any measures of
uncertainty.
The purpose of this article is threefold. First, we pre-
sent a new approach to assessing discounting and esti-
mating parameters with logistic regression, to provide a
direct estimate of the hyperbolic parameter k and obtain
confidence intervals, using the individual choices made
on each preference question as the random variables of
interest. Second, we compare the new approach with
standard scoring on MCQ data, to determine how well
the new measure corresponds to the standard. Finally, we
use the logistic model to simulate MCQ survey data and
examine the performance of the estimation procedure
against known values of k. The procedure we develop
here has the potential to streamline analysis of MCQ
data while giving a better account for people who do not
have a perfectly defined indifference point, as well as to
reduce the amount of data needed for the computer task
assessment of discounting, all using standard statistical
software.
A Logistic Probability Model of Discounting
We examined only a simple model that was indepen-
dent of V
D
. Logistic regression (Hosmer & Lemeshow,
1989) uses a transformation of a linear score variable
(SbX) to predict a probability of a binary event:
(2)
Using logistic regression (with no intercept), we fit the
choice responses (0 for immediate, 1 for delayed), with
a logistic model using time (t) and a transformed ratio
variable (1 2 1/R) as predictors, to give us
(3)
We were interested in the value of the reward ratio
(R 5 R*) that sets the probability of taking the delayed
reward to .5. From Equation 2, p 5 .5 where SbX 5 0,
making R*
(4)
and our estimate of k is
The approach is similar philosophicallyto probit analysis,
except that we are using a logistic model of an event prob-
ability, rather than a cumulative normal distribution. Stan-
dard errors are readily estimated for

k and ln(

k), using the

delta method, after obtaining the covariance matrix(V)
for b (Hogg & Craig, 1995). We focused on ln(

k) to pre-
vent the lower confidence limit around

k from dropping
below 0. The variance for our estimate ln(

k) is then
Var[ln(

k)] 5 gVg
T
, where g 5 ln(

k)/b 5 [21/b
1
,
1/b
2
], and the standard error is SE
D
5 Var[ln(

k)]. Stan-
dard errors for ln(

k) may also be obtained using the jack-

knife procedure (SE
JK
), by systematic deletion of one
survey question (Efron & Tibshirani, 1993). (See the ap-
pendices for the code for estimating k, SE, and SE
JK
.)
METHOD
Our purpose was to examine the performance of a new
estimation procedure for rapidly fitting the hyperbolic
delay-discounting function to smaller quantities of uni-
form survey data. The data came from two sources. First,
we utilized data that were collected as part of a longitu-
dinal study of health-relatedbehavior among high school
students in the Washington, D.C. metropolitan area. The
MCQ was administered in the spring of the sophomore
year, and 1,073 students completed the survey. Estimates
of k were generated using both the standard procedure

. k =
b
b
2
1
R
t
* , =
+
1
1
2
1
b
b
SbX = -

+ b b
1 2
1
1
R
t.
p
e
e
=
+
S
S
b
b
X
X
1
.
ESTIMATING DELAY-DISCOUNTING FUNCTIONS 43
and logistic regression. Finally, the estimates from the
two methods were compared using correlation methods.
Although real data allow us to compare the output of
the new technique with a standard approach, real data
cannot tell us how well the procedure works in an ab-
solute way. To address that question, we simulated ran-
dom MCQ data, using Monte Carlo methods based on
Equations 2 and 3, testing the procedures ability to iden-
tify the original parameters used in generating the simu-
lated data. We generated 1,200 surveys for each combi-
nation of parameters, where k 5 [0.002 0.01 0.02 0.1
0.25] and b
1
5 [20.5 21 23 25 28], for a total of
30,000 simulated surveys. We assumed that k was uni-
form across all V
D
levels and estimated only a single
value. We applied the logistic regression procedure to
each survey and calculated an upper and a lower confi-
dence limit for ln(

k), using SE
D
and SE
JK
. We also
recorded all related events, such as success or failure to
converge on an estimate and inclusion of the true value
of ln(k) inside each of the confidence intervals (CIs). Fi-
nally, we calculated summary statistics on our output
data, includingthe bias for

k and ln(

k), the average length

of the CI, the probability of successful convergence, and
the coverage probability for each of the CIs.
Software
We used two software packages. First, we wrote a rou-
tine for standard scoring of the MCQ survey, using the
Matlab programming language (Mathworks, Natick,
MA). The routine duplicated the logic in a spreadsheet
that had been written for scoring the survey (Kirby,
2000), but we included loop statements for f iling
through a large data set and recording output values.
We used Stata (Stata Corporation, College Station,
TX) for estimation of the k parameter, includingthe logit
and the generalized linear models (GLM) procedures
(see Appendix A). We implementedthe delta method for
estimating CIs around k, utilizing the covariance matrix
available from the system variables. We also used the
jackknife procedure to obtain an alternate SE and CI.
Stata also provides a flexible facility for maximum-
likelihood(ML) estimation, allowingus to tailor likelihood
functions that fit our needs. We used the Stata ML to place
boundaries on parameter values, using a penalty function
that came into play any time parameters stepped beyond
the boundaries during optimization (see Appendix B).
RESULTS AND DISCUSSION
Whole Population Analysis
The standard scoring procedure was applied to the stu-
dent survey data, yielding a full 1,073 k estimates for
each of the three V
D
levels. We averaged the ln(

k) values
across the V
D
levels, to yield a mean ln(

k) of 24.53
(SE 5 0.046; CI, 24.623 to 24.443). This is approxi-
mately halfway across the range of the standard assess-
ment (28.740336 to 21.39938). An analysis of the
whole population with the Stata logit procedure yielded
b
1
51.67 and b
2
520.0184, for a population-averaged

k of 0.011, or ln(

k) of 24.51. Using the delta method, we

obtained the variance of ln(

k), VAR
ln(

k)
5 0.00173, for
an SE of 0.042, yielding a 95% CI of 24.593 to 24.427.
Overall, the whole population measures of ln(

k) appear
very similar.
Comparing Logistic With Standard Scoring
We used the Stata GLM procedure to score the student
MCQ survey data, this time analyzing by subject. (Im-
plementation of logistic regression is different for logit
and GLM, and GLM works better on the small samples
of the MCQ.) The GLM procedure determined estimates
for 986 subjects, with a mean ln(k) of 24.51 and a range
of 210.58 to 20.95. (The range for standard scoring is
artificially restricted to 28.74 to 21.40. For extreme
questionnaires, almost all one choice, the procedure can
project an estimate outside of this range.) A number of k
estimates (87) were out of bounds (negative) or failed to
converge. For

k to be positive, both b
1
and b
2
must be
negative.
Logistic-regression based estimates are plottedagainst
their corresponding standard estimates in Figure 1A.
Most of the points lie along the 45 line. Some small ex-
ceptions occur at the extremes, where our point estimate
is free to reach more extreme values than are allowed
under Kirbys (2000) scoring, giving the relationship a
slightly sigmoid appearance. The overall correlation
between the measures is .944. However, if we restrict the
analysis to those with a b
1
, 21, the correlation be-
comes .972. As we shall discuss later, b
1
,21 indicates
a better defined indifference response, and as is indi-
cated by the hollowpoints in Figure 1A, most of the out-
liers were b
1
$ 21.
It is useful nowto discuss the (optimization) pathology
that prevented 87 surveys (8% of the total) from being
scored by Stata GLM. To examine whether the subjects
were responding sensibly, we ordered the questions by
ratio and recorded a regression slope (choice: delayed 5
1 vs. order). We expected a negative trend across the
order. All of the 986 surveys fitted by Stata showed a
negative trend. Of the 87 not fitted, 13 showed a slight
positive trend. This implies that the subjects choices
were changing in the wrong direction, so that they
were choosing the immediate reward more often when it
was small and rejecting it more often when it was closer
in value to the delayed reward. (This is like choosing the
delayed reward when asked, Would you like $100 right
now or in 30 days?) We cannot think of any explanation
for this, other than the fact that answering questions at
random, in the sense of a coin toss, will occasionally lead
to a positive trend. Another 48 of the surveys that failed
to converge showed a zero trend, having been answered
uniformly all the way across. Kirbys (2000) own method
of scoring deals with those 48 uniform surveys in a sen-
sible way, assigning them the high or the low extreme
values from his available range. We cannot obtain a point
estimate for these subjects without placing bounds on
the parameter values, yielding an at least or at most
k value. The remaining 26 showed a slight or vague neg-
44 WILEYTO, AUDRAIN-MCGOVERN, EPSTEIN, AND LERMAN
ative trend. When the actual patterns displayed by those
last 26 surveys were examined, the subjects had made ei-
ther random (coin toss) choices all the way across the
order or uniform choices, with one or two exceptions.
Kirbys (2000) own method determines a consistency
score for each block of nine questions, and he rejects ob-
servations with consistency scores less than 0.77. We ob-
served that over 2% of blocks have low consistency, and
Figure 1. (A) Results from the Stata generalized linear models (GLM)
procedure, comparing logistic and standard scoring. The x-axis repre-
sents the average of the three ln(k) values obtained through standard
scoring, and the y-axis represents the result of the logistic procedure.
Open circles indicate an estimate of b
1
. 21, indicating a vague or less
consistent survey response on the part of the subject. Results lie gener-
ally along the 45 line (gray), indicating a one-to-one correspondence.
The correlation is .947 if the open circles are included, and .972 if the
open circles are discarded. The GLM procedure scored 986 of 1,073 sur-
veys. (B) Similar results using Stata maximum likelihood (ML), with
bounds placed on the parameters (b
1
, b
2
,

k), using a penalty function.
Constraints allow all 1,073 surveys to be scored, and uniform surveys
(all 0 or all 1) are correctly placed at the extremes. The correlation is
.944 if the open circles are included and .972 if the open circles are dis-
carded. Standard scoring for the survey produces results for log(k) be-
tween 28.74 and 21.40, and this restricted range is indicated by gray
dashed lines in both parts of the figure. Our logistic scoring falls outside
of those lines when a set of extreme survey results causes the model to
project p 5 .5 past the standard range. This happened more often for the
constrained procedure, since the constraints made it possible to fit all of
the surveys, giving (B) a sigmoid appearance. Our constraints were
wider than standard scoring, because our purpose in constraining the
estimate was to keep k in the legal bounds for the hyperbolic function,
rather than to restrict it to range of the monetary choice questionnaire
(MCQ).
ESTIMATING DELAY-DISCOUNTING FUNCTIONS 45
3.7% (39) are inconsistent in at least one block. Exclud-
ing the uniform responses, our nonidentifiable rate is
similar to his, and the very surveys that are difficult to
evaluate using the standard method create optimization
problems for the logistic method.
We repeated the logistic regression procedure with
Stata ML, applying penalty function bounds to the pa-
rameters, in order to determine whether a constrained
optimization routine would reduce our failure rate.
Bounds were applied so that 210 , b
1
, 20.2, b
2
,
20.00001, and 0.00002 ,

k ,0.5. The ML logistic pro-

cedure generated a

k for all surveys (Figure 1B). Most of
the distribution looks the same. The exceptions are the
observations with b
1
estimated to be . 21 and the ad-
dition of those at the extreme ends of the distribution.
The bounded procedure correctly assigned the 48 uni-
form responses the most extreme values it could, subject
to the constraints, and matched the behavior of the pro-
cedure outlined by Kirby (2000).
What b
1
and b
2
Mean in the Context of the MCQ
The MCQ was designed perfectly to estimate k with a
small sample of questions. The layout of data points in
the MCQ is shown in Figure 2A. The scatter points at the
lower left correspond to an immediate reward that is very
small versus a comparatively large reward after only a
small delay. At the upper right extent of the scatter, the
reward values are almost equal, but one must wait a very
long time for the slightly larger one. Somewhere across
the range of scatter points, we would expect to see a shift
in behavior between choosing the delayed reward and
choosing the immediate reward. Kirby et al. (1999) set
up the array of questions so that the placement of the best
curve would be apparent. This arrangement works for
the logistic method, because the no-intercept (a forced 0
intercept) logistic model forces the hyperbolic y-intercept
through 1. Like Kirby, we identify the hyperbolic func-
tion from two points: the y-intercept and the point at
which it crosses the picket line.
On reviewing a large number of MCQ results, we
found that some individuals showed very clean shifts in
preference from delayed to immediate rewards, so that
when the ratio was smaller than the indifference point,
the delayed reward was always chosen, and when the
ratio was higher, the immediate reward was always cho-
sen. Others appeared to be answering at random or,
more properly, tossing a hypothetical coin for each sur-
vey item in the mid-range, sometimes showing a clean
response at the extremes, sometimes not. There were
many surveys in between these extremes.
The parameter b
1
can be thought of as a vagueness pa-
rameter, related to the clarity of the indifference response.
It represents the slope of the logistic probabilityfor choos-
ing delayed as the reward ratio increases as the probabil-
ity crosses .5 (at R*). Figure 2B shows a sample indiffer-
ence model, indicated by the bold line with y-intercept of
1. Above and below the bold line we see contours for p 5
.75 and p 5.25 that arise from two different values of b
1
.
The closer those contours are to the bold line, the steeper
the probabilitygradient across the indifference level of .5.
(The profile view is in Figure 2C.) The possible range of
b
1
is (2` , b
1
, 0), although a more realistic range is
(220 , b
1
, 20.1). Large negative values of b
1
(b
1
,
25) are associated with clear, steep behavior changes in
crossing the indifference curve (Figure 2C), and individu-
als who show clean breaks in behavior (the kind that one
expects to see with standard scoring) are represented by
those values. Small negative values are associated with
vague, flat changes (Figures 2B and 2C), and individuals
who appear to be answering at random are represented
by those small negative values. The b
1
parameter allows
us to simulate individual variation in survey behavior and
to examine the impact of a vague or clean response on our
ability to identify the model from the data.
The parameter b
2
is a time coefficient, or a scale pa-
rameter. It can be thought of as scaling the k value to
match the units on the x-axis in Figure 2A or 2B, which
are days but could as easily be hours, weeks, or months.
Performance of the Procedure
Thus far, we have demonstrated that our procedure
based on logistic regression can match the estimates
made with Kirbys own method on his discounting sur-
vey. Here, we present an examinationof the performance
of the procedure, using simulated survey data generated
from known parameter values. As was expected, not all
simulated surveys could be scored. Figure 3A reveals the
probability of obtaining an estimate based on the logis-
tic method. The method failed frequently where the true
value of b
1
was 20.5. Our earlier discussion on the
meaning of b
1
makes sense of this, because this is where
the location of the p 5 .5 contour would be the most
vague. Failures also occurred often where the true b
1
was
21, and the true k was 0.002 (small). The gradient is still
vague where b
1
is 21, and the indifference curve is very
close to the upper edge of the range of the data shown in
Figure 2A. A smaller number of failures occurred at the
other end of the k range (true k 5 0.25).
The meaningful values of

k span four orders of magni-
tude. Thus, we calculated bias [E(

u) 2 q] for log10(

k),
using successful estimates, so that bias would be pre-
sented in orders of magnitude (Figure 3B).

k appears to
have been relatively unbiased, except for a positive bias
where k was small (true k 5 0.002) and b
1
was a small
negative number (20.5, 21). The bias in those cases
may have been a direct result of the values omitted due
to estimation failures; the positive bias reflects the sur-
veys that are, by chance, more like the middle of the k
range and easier to estimate. Bias was reduced (Fig-
ure 3B) using the constrained routine in Stata ML. The
most likely explanationis that the estimates from the dif-
ficult data sets were included, rather than excluded, and
were important in centering the distribution.
The delta method 95% CI did not perform well over
much of the parameter range tested (Figure 4A). The prob-
ability of including the true value of k (coverage probabil-
46 WILEYTO, AUDRAIN-MCGOVERN, EPSTEIN, AND LERMAN
Figure 2. (A) How standard scoring of the monetary choice questionnaire (MCQ)
identifies k. Points indicate the reward ratio values in the survey versus delay times,
and they correspond to hyperbolic functions that have specific values of k. The esti-
mate of k is identified by looking for a shift in behavior. At the lower left, the imme-
diate reward is worth very little, as compared with the delayed reward, and subjects
choose the delayed reward almost exclusively, whereas at the upper right, the imme-
diate reward is almost as valuable as the delayed reward, and they choose the im-
mediate reward almost exclusively. The shift in behavior indicated here corresponds
to a k of approximately 0.01. Standard scoring has a more difficult time with re-
sponses that do not break cleanly, which are quite likely because the choices are given
out of order and respondents shoot from the hip. (B) Our approachdoes not require
a clean break, because we model the probability of choosing the delayed reward and
identify where the probability is .5. The contours reflect a gradual change in the
probability that the participant will choose delayed rewards, and the slope of this
gradual shift in behavior is reflected by parameter b
1
. If the change in behavior is
sharper, b
1
is a large negative number. If the change is more gradual, b
1
is a smaller
negative number. The bold line represents a probability contour of p 5 .5 for choos-
ing the delayed reward. Additional lines show the .75 and .25 probability contours
under two different values of b
1
. (C) Here, we display the same model of a gradual
shift in behavior in profile, indicating elevations along the dashed line in contour
view. It shows gradients crossing the p 5 .5 contour for various values of b
1
. Large
negative values may not be separated easily with small samples. Statas generalized
linear models procedure seems to work better than the logit procedure for these ex-
treme values.
ESTIMATING DELAY-DISCOUNTING FUNCTIONS 47
ity) was closest to .95 for b
1
520.5 and, generally, .4 .5
for lower values of b
1
. This may be related to the small
sample obtained by a single survey and the excessively
small standard errors generated by the delta method under
some parameter combinations. However, the jackknife SE
(SE
JK
), used in conjunctionwith the GLM procedure, pro-
duced CIs with consistently good coverage (Figure 4B),
generally falling in the 89% 97% range. We also applied
the jackknife procedure to the boundedMLprocedure. The
initial assessment of coverage included values as low as
70%, higher than the worst given by the delta method but
not close to the coverage provided when the jackknife pro-
cedure was applied to Stata GLM. This was remedied by
removing from considerationCIs calculatedfor point esti-
mates made with constraints active, which left the remain-
ing coverage close to 95% across the board (Figure 4C).
SE
JK
can be exceedingly large, and the distribution of
SE
JK
can be highly skewed. For that reason, Figure 4D
displays the median length of the 95% CI resulting from
SE
JK
, revealing two trends. First, the median size of the CI
is about the same magnitude as true k, and correspond-
ingly, the CI generally increases as true k increases. The
Figure 3. Simulation results. The x-axis shows the true value of ln(k), with separate
lines showing the response for different values of b
1
. (A) Proportion of simulated sur-
veys that could be scored using Stata generalized linear models. Two things con-
tributed to failures: k
true
values at either extreme were near the edge of the surveys
ability to detect them, and small negative numbers for b
1
(solid line) produced a vague
response. (B) Bias for a parameter u is [E(

u) 2 u]. A bias of 0 is perfect. Because the
range of k can be several orders of magnitude, we calculated bias as {E[log
10
(

k)] 2
log
10
(k)}. It was negligible except where b
1
was a small negative number and the prob-
ability of deleting one of the more extreme observations was high. (C) Bias for the
boundedestimate (from Stata maximumlikelihood) was reduced because extreme ob-
servations were included in the average.
48 WILEYTO, AUDRAIN-MCGOVERN, EPSTEIN, AND LERMAN
second trend is that high values of b
1
tend to have larger
CIs than do low values of b
1
. The high values of b
1
are
more vague responses to the series of survey questions
(Figure 2A). One exception, where the CI is unusually
large, occurs where b
1
is 20.5 and k is 0.002. This com-
bination implies a vague response on the part of the sub-
ject, and a true k value near the edge of the surveys abil-
ity to detect.
Figure 4. Performance of the confidence interval. The x-axis shows the true
value of ln(k), with separate lines showing the response for different values of
b
1
. (A) Coverage of the confidence interval derived via the delta method falls
far short of the purported 95%. The worst performance comes from what
should be the best data (clean behavioral shifts and k
true
in the mid-range).
(B) Coverage of the confidence interval derived via the jackknife remains close
to 95% across all parameter values tested. We recommend the jackknife pro-
cedure. (C) Coverage of the confidence interval derived via the jackknife when
bounds are applied to the parameters remains close to 95% across all param-
eter values tested, as long as standard errors are discarded when the constraint
is active. Coverage reached as low as 70% when observations with active con-
straints were included. We recommend ignoring the standard error estimate
when any constraint is active for the point estimate. (D) The median width of
jackknife confidence intervals for the constrained procedure (corresponding to
that in panel C). Intervals are widest for large values of k
true
, and for a small
negative b
1
.
ESTIMATING DELAY-DISCOUNTING FUNCTIONS 49
Limitations of Our Analysis
Althoughour methodof analysis can be appliedto com-
puter task data, as well as to any fixed format survey de-
signed to assess discounting, we restricted our analysis to
data that would be found in the MCQ (Kirby et al., 1999),
becauseit is incurrent use and amenableto simulation. Our
restriction to MCQ data was also appropriate because this
survey is growing in popularity; it is the only standard sur-
vey for assessing delay discounting and is being placed in
a large number of current protocols because it is compact
(Audrain-McGovern et al., 2004; K. N. Kirby, personal
communication; N. M. Petry, personal communication;
M. Seligman, personal communication). What we report
is probably the most extreme test of the fitting technique,
and better results would be expected with larger formats.
For simplicity in this initial analysis, we treated dis-
countingas if it were uniform, regardless of the magnitude
of the delayed reward. Nevertheless, we know that k is
smaller if the delayed reward is more valuable (Green
et al., 1997; Johnson & Bickel, 2002; Kirby, 1997; Kirby
& Marakovi c, 1995) and that future efforts will need to
account for this. Kirbys (2000) own method accounts for
this by determining an estimate for k over each of three
blocks of nine questions. Our failure rate when attempt-
ing to estimate k for each block of questions with the
simpler unconstrained routines was too high to support
that approach, and for that reason, we focused on aver-
ages. However, the constrained routine was reliable in
providing an estimate by block, with coefficients of de-
termination for regressions of our results on those from
Kirbys (2000) own method at .88, .89, and .90 for small,
medium, and large reward blocks, respectively. With an
expanded format for this survey, we expect that this ap-
proach will enable us to fit more complex models that in-
corporate V
D
as a covariate and to model the effect of the
magnitudeof the delayed reward. (The multiplicativefor-
mulation described by Ho, Mobini, Chiang, Bradshaw, &
Szabadi, 1999, appears to work with the same approach.
Another suggestion would be to replace the hyperbolic
parameter k with a linear model of k: k 5 k
0
1bV
D
.)
Finally, we recommend using the constrained proce-
dure when possible and the jackknife procedure for ob-
taining CIs. The interval estimates are based on SEs that
are inaccurate if the constraints are active. The user
should examine estimates of b
1
, b
2
, and k, determine
whether any of them are at the constrained values, and
ignore SEs and CIs calculated for those observations.
Further Research
The technique described here will probably see its
greatest use in scoring larger fixed discounting surveys
and reducing the size of computer task assessments of
discounting. Many of those larger surveys and tasks were
constructed to ascertain indifference, to be used in de-
termining what function best fits the data, as well as in
estimating the parameters. Most authors now assume
that the best function to fit is hyperbolic, so gathering
data sets large enough to determine the best-fitting func-
tion is not necessary.
We have conductedsome preliminary trials (by simula-
tion) of slightly expanded formats for the survey, includ-
ing simple range extensions of the current layout (includ-
ing smaller ratios at the left end and longer times at the
right end of the scheme in Figure 2C and different lay-
outs). We have seen some initial improvements in the abil-
ity to identify the model at the edges of the usual k range,
but a full analysis is outside the scope of this article.
REFERENCES
Audrain-McGovern, J., Rodriguez, D., Tercyak, K. P., Epstein,
L. H., Goldman, P., & Wileyto, E. P. (2004). Applyinga behavioral
economic framework to understanding adolescent smoking. Psy-
chology of Addictive Behaviors, 18, 64-73.
Bickel, W. K., Madden, G. J., & Petry, N. M. (1998). The price of
change: The behavioral economics of drug dependence. Behavior
Therapy, 29, 545-565.
Bickel, W. K., Odum, A. L., & Madden, G. J. (1999). Impulsivity and
cigarette smoking: Delay discountingin current, never, and ex-smokers.
Psychopharmacology, 146, 447-454.
Crean, J. P., de Wit, H., & Richards, J. B. (2000). Reward discount-
ing as a measure of impulsive behavior in a psychiatric outpatient
population. Experimental & Clinical Psychopharmacology, 8, 155-
162.
Efron, B., & Tibshirani, R. J. (1993). An introductionto the bootstrap.
New York: Chapman & Hall.
Green, L., Myerson, J., & McFadden, E. (1997). Rate of temporal
discounting decreases with amount of reward. Memory & Cognition,
25, 715-723.
Ho, M. Y., Mobini, S., Chiang, T. J., Bradshaw, C. M., & Szabadi, E.
(1999). Theory and method in the quantitative analysis of impulsive
choice behaviour: Implications for psychopharmacology. Psycho-
pharmacology, 146, 362-372.
Hogg, R. V., & Craig, A. T. (1995). Introduction to mathematical sta-
tistics. Englewood Cliffs, NJ: Prentice-Hall.
Hosmer, D. W., & Lemeshow, S. (1989). Applied logistic regression
(1st ed.). New York: Wiley.
Johnson, M. W., & Bickel, W. K. (2002). Within-subject comparison
of real and hypothetical money rewards in delay discounting. Jour-
nal of the Experimental Analysis of Behavior, 77, 129-146.
Kirby, K. N. (1997). Bidding on the future: Evidence against normative
discounting of delayed rewards. Journal of Experimental Psychol-
ogy: General, 126, 54-70.
Kirby, K. N. (2000). Instructions for inferring discount rates from
choices between immediate and delayed rewards. Unpublished man-
uscript.
Kirby, K. N., & Marakovic, N. N. (1995). Modelingmyopic decisions:
Evidence for hyperbolic delay-discounting within subjects and
amounts. Organizational Behavior & Human Decision Processes,
64, 22-30.
Kirby, K. N., & Marakovic, N. N. (1996). Delay-discounting proba-
bilistic rewards: Rates decrease as amounts increase. Psychonomic
Bulletin & Review, 3, 100-104.
Kirby, K. N., Petry, N. M., & Bickel, W. K. (1999). Heroin addicts
have higher discount rates for delayed rewards than nondrug-using
controls. Journal of Experimental Psychology: General, 128, 78-87.
Madden, G. J., Bickel, W. K., & Jacobs, E. A. (1999). Discounting of
delayed rewards in opioid-dependent outpatients: Exponential or
hyperbolic discounting functions? Experimental & Clinical Psycho-
pharmacology, 7, 284-293.
Madden, G. J., Petry, N. M., Badger, G. J., & Bickel, W. K. (1997).
Impulsive and self-control choices in opioid-dependent patients and
non-drug-using control participants: Drug and monetary rewards.
Experimental & Clinical Psychopharmacology, 5, 256-262.
50 WILEYTO, AUDRAIN-MCGOVERN, EPSTEIN, AND LERMAN
Petry, N. M. (2001a). Delay discounting of money and alcohol in ac-
tively using alcoholics, currently abstinent alcoholics, and controls.
Psychopharmacology, 154, 243-250.
Petry, N. M. (2001b). Pathological gamblers, with and without sub-
stance use disorders, discount delayed rewards at high rates. Journal
of Abnormal Psychology, 110, 482-487.
Richards, J. B., Mitchell, S. H., de Wit, H., & Seiden, L. S. (1997).
Determination of discount functions in rats with an adjusting-amount
procedure. Journal of the Experimental Analysis of Behavior, 67,
353-366.
NOTE
1. Our own use of this software has generally required approximately
120 questions to identify all indifference points. The actual number is
indeterminate and may be different for each subject.
APPENDIX A
Stata Code for Estimating k Using Logit and GLM
Data must be in long format, with each question item on a separate line, 27 lines per person. We use inte-
ger variables idn for subject, choice for outcome (0,1), and (continuous) variables ratio as reward ratio
V
I
/V
D
, time as delay time, and ip1 as the constructedvariable (1 21/Ratio). The best routine in Stata for
our purposes is glm, which fails to converge less often than logit. The entry for estimating the model across
the whole sample is:
glm choice ip1 time, nocons f(binomial) l(logit)
To collect individual estimates frommany subjects, each identified by idn, we use the STATSBY com-
mand (/* and */ continue a command across lines):
statsby glm choice ip1 time, nocons f(binomial) l(logit) /*
*/ _b lkhat 5 (log(_b[time]/ _b[ip1])), by(idn) clear
So far, our code generates only estimates of k for each subject. In order to use the jackknife procedure to
estimate standard errors, more coding is required. First, install the levels command (from the ssc Web
archive) if you have not already done so.
ssc install levels
We then GENERATE empty variables to receive our results, create a local list of the idn values, use a loop
(FOREACH) to apply the JACKKNIFE procedure to each subject, saving results for each successful run,
and COLLAPSE to one line per subject.
gen b1 5 .
gen b2 5 .
gen lkhat 5 .
gen lkerr 5.
levels idn, local(idnlevs)
foreach i of local idnlevs {
capture jknife glm choice ip1 time, nocons f(binomial) l(logit) /*
*/ _b lkjk 5(log(_b[time]/_b[ip1])) if idn 55 `i,e
scalar rc3 5_rc
replace b1 5 r(stat1) if idn 55 `i & rc3 55 0
replace b2 5 r(stat2) if idn 55 `i & rc3 55 0
replace lkhat 5r(stat3) if idn 5 5`i & rc3 5 50
replace lkerr 5 r(se3) if idn 5 5`i & rc3 550
}
collapse b1 b2 lk*, by(idn)
The confidence interval is then the estimate lkhat, 61.96 lkerr.
ESTIMATING DELAY-DISCOUNTING FUNCTIONS 51
APPENDIX B
Routines for Logistic Regression Using Stata ML, Constraining Parameter Values
First specify the likelihood function, including a penalty for violating bounds (copy exactly).
program define hyplogmlec
version 7
args todo b lnf
tempvar theta1 tmp1 myb2 d1 d2 d3 d4 d5
capture matrix `myb2 5`b
mleval `theta1 5 `b,eq(1)
quietly gen double `tmp1 5exp(`theta1)/(1 1 (exp(`theta1)))
mlsum `lnf 5$ML_y1*ln(`tmp1) 1(1-$ML_y1)*ln(1-`tmp1)
note: Calculation of penalties.
scalar `d1 5 ((`myb2[1,2]/`myb2[1,1]) /*
*/ , 0.0001)*((`myb2[1,2]/`myb2[1,1])- 0.0001)^2
scalar `d2 5 ((`myb2[1,2]/`myb2[1,1]) /*
*/ . 1)*((`myb2[1,2]/`myb2[1,1])-1)^2
scalar `d3 5 ((`myb2[1,1]) , -10)*((`myb2[1,1] ) 110)^2
scalar `d4 5 ((`myb2[1,1]) . -.2)*((`myb2[1,1]) 1 .2)^2
scalar `d5 5 ((`myb2[1,2]) . -.001)*((`myb2[1,2]) 1 .001)^2
scalar `lnf 5`lnf -(`d2 1 `d4 1`d1 1 `d3 1`d5)*5000
if `lnf 55 . { exit}
end
Again, we record the estimates, and the jackknifestandard error of log(

k), using a routine that calls the like-

lihood function. The ID variable idn may represent subject or subject by delayed reward size block (small,
medium, large). The constrained procedure also provides estimates for the smaller nine-item sample. Com-
plete routines for STATA may be found in the Statistical Software Components Archive (maintainedby Boston
College): http://ideas.repec.org/s/boc/bocode.html.
gen b1 5 .
gen b2 5 .
gen lkhat 5 .
gen lkerr 5.
levels idn, local(idnlevs)
foreach i of local idnlevs {
capture jknife /*
*/ ml model d0 hyplogmlec (choicevar 5ip1 time, nocons),max /*
*/ _b bias 5(log(_b[time]/_b[ip1])) if idn 55 `i,e
scalar rc3 5_rc
replace b1 5 r(stat1) if idn 55 `i & rc3 55 0
replace b2 5 r(stat2) if idn 55 `i & rc3 55 0
replace lkhat 5r(stat3) if idn 5 5`i & rc3 5 50
replace lkerr 5 r(se3) if idn 5 5`i & rc3 550
}
collapse b1 b2 lk*, by(idn)
(Manuscript received October 1, 2002;
revision accepted for publicationSeptember 9, 2003.)