delay_discounting_fit()

Introduction

Delay discounting refers to the decline in value of a reward based on the time to its delivery. The steepness of delay discounting tends to be consistent within individuals, for which this attribute counts as a behavioral trait (Odum, 2011).

Delay discounting data are typically obtained from protocols involving choices between an immediate smaller reward and a delayed larger reward (i.e., intertemporal choice tasks). By dynamically varying the delay or magnitude of some of these alternatives, one can obtain indifference points where the subjective value of both alternatives is assumed to be equivalent.

The function describing the decay of indifference points along the values of the delay for one of the alternatives can be analyzed with many different parametric models. Nonetheless, the most used by far are the exponential and the hyperbolic models. The normalized exponential and hyperbolic models, respectively, are expressed as:

• \(V(d)=exp(-k*d)\)
• \(V(d)=\frac{1}{(1+k*d)}\)

Where \(V(d)\) is the subjective value as a function of delay \(d\), and \(k\) the discount parameter for both models.

However, the data can be similarly described by the different models, which render the problem of model-selection almost impossible to solve only by means of statistical adequacy. This is why some authors have proposed a model-agnostic index for the delay discounting function, the area under the curve (Myerson et al., 2001).

Here we introduce three separate functions for applying both the hyperbolic and exponential fit methods, as well as the AUC method to delay discounting data.

The hyperbolic_fit() and exp_fit() functions take the following parameters:

• value a numeric vector of the subjective values (indifference points).
• delay a numeric vector of the delays used.
• initial_guess a numeric value providing an initial start for \(k\).
• max_iter a positive integer for the maximum number of iterations for the adjusting process (this is optional, the default value is 100000).
• scale_offset A constant to be added if the residuals are close to 0, this is to avoid division by 0, which is know to cause problems of convergence (0 by default).

The trapezoid_auc() function take the following parameters:

• x a numeric vector of delay values.
• y a numeric vector of subjective values (indifference points).

NOTE: Values are not normalized, so if you need normalized results, use unit_normalization()on your data

Example

First let’s load an example data set from an experimental session:

data("DD_data")
norm_sv <- DD_data$norm_sv
delays <- DD_data$Delay

DD_data

##   norm_sv Delay
## 1   0.934     1
## 2   0.746     6
## 3   0.746    12
## 4   0.488    36
## 5   0.684    60
## 6   0.441   120

Now let’s fit lineal, hyperbolic and exponential models for comparison purposes and calculate the Akaike criterion for each one:

# first, fit a linear model
lineal_m <- lm(norm_sv ~ delays)
# hyperbolic model
hyp_m <- hyperbolic_fit(norm_sv, delays, 0.1)
# exponential model
exp_m <- exp_fit(norm_sv, delays, 0.1)
AIC(lineal_m, hyp_m, exp_m)

##          df       AIC
## lineal_m  3 -4.089953
## hyp_m     2 -3.534409
## exp_m     2 -1.427746

Notice we omitted the max_iter value, therefore using the default one.

Now we extract the \(k\) value from each model:

k_hyp <- coef(hyp_m)
k_exp <- coef(exp_m)
k_lin <- coef(lineal_m)

## [1] "K_hyp:  0.0154943251678617"

## [1] "K_exp:  0.00941780689009577"

## [1] "K_lin:  0.796330986258254"    "K_lin:  -0.00314462092574266"

Additionally, we can calculate the AUC using the trapezoid_auc() function like this:

delay_norm <- delays / max(delays) # It is important to normalize the delay values first in order to get a coherent AUC.
AUC_value <- trapezoid_auc(delay_norm, norm_sv)

Finally let’s plot the resulting fitted curves:

Note: to create the data points for the hyperbolic fit curve we used the eq_hyp() function included in the package which simulate hyperbolic delay discounting data given a specific \(k\) value and a delay vector. The function takes two parameters: k and delay. The function is used like this:

y_data <- seq(0, max(delays), len = 200)
x_data <- eq_hyp(k = k_hyp, y_data)

For the exponential fit curve, we just applied the model with the resulting k_exp value to the y_data vector like this:

x_data <- exp(-k_exp * y_data)