Additive aging curve

Warning

This post is a draft.

Recently, I helped design an experiment measuring a binary response against a continuous variable. If the user abandoned their cart at time zero, then we delayed for a variable number of minutes before reminding them to finish their purchase. The delay has a non-monotonic relationship to the response: as the delay increases, so does the purchase rate; then the rate peaks; and finally it decreases.

Causally, we may decompose this process into two: as the delay increases, the user 1) becomes more available for and 2) loses interest in purchasing the product. This is a common phenomena in different time-based scenarios. In sports, the “aging curve” refers to how a player’s performance increases with age, then decreases. As the player gets older, they get 1) better at the sport and 2) physically weaker.

Andrew Gelman wrote about this a couple of times in his blog: see his posts from 2018 and 2023, where Gelman suggests modeling these processes using an additive function like:

\[g(x) = g_1(x) + g_2(x),\]

where
\(g_1(x)\) is a monotonically increasing function with a right asymptote; and
\(g_2(x)\) is a monotonically decreasing function with a left asymptote.

In this post, we’ll analyse an experimental dataset by fitting and comparing three different models: a non-parametric bootstrap, a semi-parametric spline and a fully parametric decomposable curve like \(g(x)\).

The Digit Span test

The motivation for Gelman’s 2018 post was a study relating age to peak cognitive functioning (Hartshorne and Germine 2015). According to the study, one of their experiments was a large scale online experimentation platform:

Hartshorne, Joshua K., and Laura T. Germine. 2015. “When Does Cognitive Functioning Peak? The Asynchronous Rise and Fall of Different Cognitive Abilities Across the Life Span.” Psychological Science 26 (4): 433–43. https://doi.org/10.1177/0956797614567339.

Participants in Experiment 2 (N = 10,394; age range = 10–69 years old) […] were visitors to TestMyBrain.org, who took part in experiments in order to contribute to scientific research and in exchange for performance-related feedback. […] We continued data collection for each experiment for approximately 1 year, sufficient to obtain around 10,000 participants, which allowed fine-grained age-of-peak-performance analysis.

The dataset for Experiment 2 is available online (Germine and Hartshorne 2016) and includes results of the Digit Span verbal working memory test, part of the Wechsler Adult Intelligence Scale (WAIS) and Wechsler Memory Scale (WMS) supertests. In the Digit Span test, subjects must repeat lists of digits, either in the same or reversed order.

Germine, Laura, and Joshua K Hartshorne. 2016. “Hartshorne & Germine (2015) When Does Cognitive Functioning Peak?” OSF. osf.io/f2saj.

Let’s plot the relationship between age and Digit Span performance:

Visually, it’s still unclear if this relationship follows an aging curve, but we’ll get back to this matter in the next section.

Bootstrap estimates

In the original paper, the authors describe a bootstrap resampling procedure to estimate the distribution of ages of peak performance:

Estimates and standard errors for age of peak performance were calculated using a bootstrap resampling procedure identical to the one used in Experiment 1 but applied to raw performance data. To dampen noise, we smoothed means for each age using a moving 3-year window prior to identifying age of peak performance in each sample. Other methods of dampening noise provide similar results.

Let’s decompose this method (as I understand it) into steps:

1. With replacement, sample \(n\) observations from the dataset;
2. Calculate the mean performance for each sample and age;
3. Repeat steps 1 and 2 \(m\) times to get multiple samples;
4. Sort each sample by age and smooth age means using a 3-year rolling average;
5. Find the age of peak performance for each sample.

import polars as pl

n = experiment.height
m = 10000
nm = n * m
seed = 37
samples = (
    experiment.sample(nm, with_replacement=True, seed=seed)
    .with_columns(sample=pl.arange(1, nm + 1) % m)
    .group_by("sample", "age")
    .agg(mean=pl.col("y").mean())
    .sort("sample", "age")
    .with_columns(smoothed_mean=pl.col("mean").rolling_mean(3).over("sample"))
)
peak = samples.group_by("sample").agg(
    age=pl.col("age").get(pl.col("smoothed_mean").arg_max())
)

This yields the following bootstrap distribution of ages of peak performance:

This distribution suggests two important things:

1. The most probable age of peak performance is 33;
2. There is a non-negligible probability that the age of peak performance happens in the early 20s, but a negligible probability that it happens in the late 20s.

Thing 2 certainly deserves attention. This is possibly caused by a confound variable or some measuring error, but I won’t investigate this any further. Instead, let’s get back to estimating curves. We will use the samples from step 4 to summarize the distribution of mean performances. For each age, we calculate the mean and 90% interquantile range, yielding a nonparametric curve:

This figure is analogue to figure 3a in the paper. This curve is entirely empirical, so there isn’t much to interpret here. It does indicate, however, the rising and falling behavior of an aging curve, with a steep increase from age 10 to 20, followed by a plateau between 20 and 30, and a slow decline beginning at 40.

Penalized splines

Splines are wiggly curves…

\[ \begin{align} g(x) &= \alpha + Z \bf{b} \\ y &\sim \mathrm{Normal}(g(x), \sigma) \\ \alpha &\sim \mathrm{Student}(3, 0, 0.1) \\ \sigma &\sim \mathrm{HalfCauchy}(1) \end{align} \]

Polynomials… (Kolassa 2017)
Adaptive smoothing…

Kolassa, Stephan. 2017. “B-Splines VS High Order Polynomials in Regression.” Cross Validated. https://stats.stackexchange.com/q/314878.

\[ \begin{align} b &= \tau \bf{z} \\ \tau &\sim \mathrm{HalfCauchy}(1) \\ \bf{z} &\sim \mathrm{Normal}(0, 1) \end{align} \]

https://www.pymc.io/projects/examples/en/latest/howto/spline.html
https://www.tjmahr.com/random-effects-penalized-splines-same-thing/

import pymc as pm

with pm.Model() as spline:
    Z = pm.ConstantData("Z", Z)
    α = pm.StudentT("α", 3, 0, sigma=0.1)
    τ = pm.HalfCauchy("τ", 1)
    z = pm.Normal("z", 0, 1, size=Z.shape[1])
    b = pm.Deterministic("b", τ * z)
    μ = pm.Deterministic("μ", α + pm.math.dot(Z, b.T))
    σ = pm.HalfCauchy("σ", 1)
    pm.Normal("y", μ, σ, observed=y)

When, where and how things change… https://www.youtube.com/watch?v=Zxokd_Eqrcg&t=506s

Twin peaks!

Additive function

\[ \begin{align} g_1(x) &= \alpha + \beta_1 \exp(-\lambda_1 x) \\ g_2(x) &= \beta_2 \exp(\lambda_2 x) \\ \\ g(x) &= g_1(x) + g_2(x) \\ &= \alpha + \beta_1 \exp(-\lambda_1 x) + \beta_2 \exp(\lambda_2 x) \\ \\ y &\sim \mathrm{Normal}(g(x), \sigma) \\ \alpha &\sim \mathrm{Normal}(0, 2) \\ \lambda &\sim \mathrm{Exponential}(0.01) \\ \sigma &\sim \mathrm{Exponential}(1) \\ \end{align} \]

import pymc as pm


def g_1(x):
    return α + β[0] * pm.math.exp(-λ[0] * x)


def g_2(x):
    return β[1] * pm.math.exp(λ[1] * x)


def g(x):
    return g_1(x) + g_2(x)


with pm.Model() as additive:
    x = pm.ConstantData("x", x)
    α = pm.Normal("α", 0, 1)
    β = pm.Normal("β", 0, 1, size=2)
    λ = pm.HalfNormal("λ", 0.004, size=2)
    μ = pm.Deterministic("μ", g(x))
    σ = pm.HalfNormal("σ", 1)
    pm.Normal("y", mu=μ, sigma=σ, observed=y)