Evaluation of different double-exponential models

# load the necessary libraries
library(tidyverse)
library(minpack.lm)

Models

The following model appears in Van den Bos & McClure (2013):

\[ V = \omega\delta_1^D + (1-\omega)\delta_2^D. \tag{1} \label{eq1} \]

Another model referenced as the double-exponential model proposed by Van den Bos & McClure appears in Green et al. (2013) and Ballard et al. (2023):

\[ V = A[\omega e^{(-\delta_1D)} + (1-\omega)e^{(-\delta_2D)}]. \tag{2} \label{eq2} \]

I noticed that the model I fit to the AO data was similar to Equation \(\ref{eq1}\), but I accidentally included A where it is not actually included by Van den Bos & McClure. So the model I fit in the lab meeting was:

\[ V = A[\omega\delta_1^D + (1-\omega)\delta_2^D]. \tag{3} \label{eq3} \]

I think the exclusion of A in Equation \(\ref{eq1}\) by Van den Bos & McClure is assuming the data is transformed such that A is 1 (This seems to be supported by how bad the fit is).

Data

I used the average AO data to do some more testing of whether these models are equivalent. The data and the analysis is shown below.

# get the data
analysis_df <- data.frame(
  Delay = c(10, 20, 30, 40, 50, 60, 70, 80),
  RawIndiffPt = c(
    65.03900390039004,
    99.46794679467946,
    117.46524652465247,
    134.9199919991999,
    139.46194619461946,
    146.65466546654665,
    150.31303130313032,
    154.2814281428143
  )
)

# create a new column with the transformed values
analysis_df$TransformedIndiffPt <- min(analysis_df$RawIndiffPt) + max(analysis_df$RawIndiffPt) - analysis_df$RawIndiffPt

# show the data
knitr::kable(analysis_df)

Delay	RawIndiffPt	TransformedIndiffPt
10	65.03900	154.28143
20	99.46795	119.85249
30	117.46525	101.85519
40	134.91999	84.40044
50	139.46195	79.85849
60	146.65467	72.66577
70	150.31303	69.00740
80	154.28143	65.03900

Fitting the models to the average AO data

# Set A by subtracting 15 (FDF mean of LL) from the min + max
A <- min(analysis_df$TransformedIndiffPt) + max(analysis_df$TransformedIndiffPt) - 15

# double-exponential formulas (used for fitting)
equation_1_formula <- TransformedIndiffPt ~ w * d1^Delay + (1 - w) * d2^Delay

equation_2_formula <- TransformedIndiffPt ~ A * (w * exp(-d1 * Delay) + (1 - w) * exp(-d2 * Delay))

equation_3_formula <- TransformedIndiffPt ~ A * (w * d1^Delay + (1 - w) * d2^Delay)

# fit the models
equation_1_fit <- nlsLM(
  formula = equation_1_formula, # forumula used for fitting
  data = analysis_df, # data to fit (column names must match the formula)
  start = list(w = 0.257, d1 = 0.395, d2 = 0.965) # initial values for the parameters
)

equation_2_fit <- nlsLM(
  formula = equation_2_formula,
  data = analysis_df,
  start = list(w = 0.257, d1 = 0.395, d2 = 0.965)
)

equation_3_fit <- nlsLM(
  formula = equation_3_formula,
  data = analysis_df,
  start = list(w = 0.257, d1 = 0.395, d2 = 0.965)
)

# calculate R-squared
sst <- sum((analysis_df$TransformedIndiffPt - mean(analysis_df$TransformedIndiffPt))^2)
equation_1_ssr <- sum(residuals(equation_1_fit)^2)
equation_2_ssr <- sum(residuals(equation_2_fit)^2)
equation_3_ssr <- sum(residuals(equation_3_fit)^2)

equation_1_r2 <- 1 - (equation_1_ssr / sst)
equation_2_r2 <- 1 - (equation_2_ssr / sst)
equation_3_r2 <- 1 - (equation_3_ssr / sst)

# make a data frame with the results
results_df <- data.frame(
  Model = c("Equation 1", "Equation 2", "Equation 3"),
  R2 = c(equation_1_r2, equation_2_r2, equation_3_r2),
  w = c(coef(equation_1_fit)["w"], coef(equation_2_fit)["w"], coef(equation_3_fit)["w"]),
  d1 = c(coef(equation_1_fit)["d1"], coef(equation_2_fit)["d1"], coef(equation_3_fit)["d1"]),
  d2 = c(coef(equation_1_fit)["d2"], coef(equation_2_fit)["d2"], coef(equation_3_fit)["d2"])
)

# show the results in a table
knitr::kable(results_df)

Model	R2	w	d1	d2
Equation 1	-7.9136890	0.5563225	1.0587799	1.0587708
Equation 2	0.9978066	0.6549981	0.0466837	0.0015134
Equation 3	0.9978066	0.3450024	0.9984878	0.9543892

Graphs

# double-exponential functions (used for plotting)
equation_1 <- function(x, w, d1, d2) {
  w * d1**x + (1 - w) * d2**x
}

equation_2 <- function(x, A, w, d1, d2) {
  A * (w * exp(-d1 * x) + (1 - w) * exp(-d2 * x))
}

equation_3 <- function(x, A, w, d1, d2) {
  A * (w * d1**x + (1 - w) * d2**x)
}

# create a data frame for graphing
graphing_df <- data.frame(
  Delay = seq(10, 80, 0.1),
  Equation1 = equation_1(seq(10, 80, 0.1), coef(equation_1_fit)["w"], coef(equation_1_fit)["d1"], coef(equation_1_fit)["d2"]),
  Equation2 = equation_2(seq(10, 80, 0.1), A, coef(equation_2_fit)["w"], coef(equation_2_fit)["d1"], coef(equation_2_fit)["d2"]),
  Equation3 = equation_3(seq(10, 80, 0.1), A, coef(equation_3_fit)["w"], coef(equation_3_fit)["d1"], coef(equation_3_fit)["d2"])
)

Equation 1

ggplot(analysis_df, aes(x = Delay, y = TransformedIndiffPt)) +
  geom_point(shape = 1, size = 3) +
  geom_line(data = graphing_df, aes(x = Delay, y = Equation1)) +
  labs(x = "Delay", y = "Value") +
  ylim(0, 160) +
  theme_classic() +
  theme(
    text = element_text(size = 20),
    axis.ticks.length = unit(0.25, "cm")
  )

Equation 2

ggplot(analysis_df, aes(x = Delay, y = TransformedIndiffPt)) +
  geom_point(shape = 1, size = 3) +
  geom_line(data = graphing_df, aes(x = Delay, y = Equation2)) +
  labs(x = "Delay", y = "Value") +
  ylim(0, 160) +
  theme_classic() +
  theme(
    text = element_text(size = 20),
    axis.ticks.length = unit(0.25, "cm")
  )

Equation 3

ggplot(analysis_df, aes(x = Delay, y = TransformedIndiffPt)) +
  geom_point(shape = 1, size = 3) +
  geom_line(data = graphing_df, aes(x = Delay, y = Equation3)) +
  labs(x = "Delay", y = "Value") +
  ylim(0, 160) +
  theme_classic() +
  theme(
    text = element_text(size = 20),
    axis.ticks.length = unit(0.25, "cm")
  )

Thoughts

It looks like Equation \(\ref{eq2}\) will fit the AO data assuming that I choose good initial parameter guesses. These guesses are different than what I used for the fit of Equation \(\ref{eq3}\) that I showed in the lab meeting. However, I think Equaiton \(\ref{eq2}\) will be best for the residual analysis since that appears with A in the literature. Also, since the \(R^2\) values are the same for Equations \(\ref{eq2}\) and \(\ref{eq3}\), I’m thinking that Equation \(\ref{eq2}\) might be a transformation of Equation \(\ref{eq1}\) that includes A. I’ll have to look into this more.

Another interesting thing I noticed about these fits is that Van den Bos & McClure (2013) stipulate that (1) both \(\delta_1\) and \(\delta_2\) should be between 0 and 1, and (2) \(\delta_1\) should be less than \(\delta_2\). The second stipulation isn’t true for the fits of Equation \(\ref{eq2}\). I’ll have to look into whether this holds for Equation \(\ref{eq2}\) or if it is just a stipulation for Equation \(\ref{eq1}\). But, if it does hold then even if it produces random residuals, I don’t think it would support the double exponential model as anything more than an empirical description.

References

Ballard, T., Luckman, A., & Konstantinidis, E. (2023). A systematic investigation into the reliability of inter-temporal choice model parameters. Psychonomic Bulletin & Review, 30(4), 1294-1322. https://doi.org/10.3758/s13423-022-02241-7

Green, L., Myerson, J., Oliveira, L., & Chang, S. E. (2013). Delay discounting of monetary rewards over a wide range of amounts. Journal of the experimental analysis of behavior, 100(3), 269-281. https://doi.org/10.1002/jeab.45

Van den Bos, W., & McClure, S. M. (2013). Towards a general model of temporal discounting. Journal of the experimental analysis of behavior, 99(1), 58-73. https://doi.org/10.1002/jeab.6