Modelling Within-Subject Variation in psychological trials with Reaction Times

gistfile1.txt

library(easystats)

set.seed(1234)
n <- 100
baseline <- rnorm(n, 5, 1.5)
d <- data.frame(
  RT1 = rnorm(n, baseline, 0.5),
  RT2 = rnorm(n, baseline, 0.5),
  RT3 = rnorm(n, baseline, 0.5),
  RT4 = rnorm(n, baseline, 0.5),
  RT5 = rnorm(n, baseline, 0.5),
  RT6 = rnorm(n, baseline, 0.5),
  RT7 = rnorm(n, baseline, 0.5),
  RT8 = rnorm(n, baseline, 0.5),
  RT9 = rnorm(n, baseline, 0.5),
  RT10 = rnorm(n, baseline, 0.5),
  ID = 1:n,
  time = 1
)

set.seed(1234)
d2 <- data.frame(
  RT1 = rnorm(n, d$RT1, 0.5),
  RT2 = rnorm(n, d$RT2, 0.5),
  RT3 = rnorm(n, d$RT3, 0.5),
  RT4 = rnorm(n, d$RT4, 0.5),
  RT5 = rnorm(n, d$RT5, 0.5),
  RT6 = rnorm(n, d$RT6, 0.5),
  RT7 = rnorm(n, d$RT7, 0.5),
  RT8 = rnorm(n, d$RT8, 0.5),
  RT9 = rnorm(n, d$RT9, 0.5),
  RT10 = rnorm(n, d$RT10, 0.5),
  ID = 1:n,
  time = 2
)

set.seed(1234)
d3 <- data.frame(
  RT1 = rnorm(n, d$RT1, 0.5),
  RT2 = rnorm(n, d$RT2, 0.5),
  RT3 = rnorm(n, d$RT3, 0.5),
  RT4 = rnorm(n, d$RT4, 0.5),
  RT5 = rnorm(n, d$RT5, 0.5),
  RT6 = rnorm(n, d$RT6, 0.5),
  RT7 = rnorm(n, d$RT7, 0.5),
  RT8 = rnorm(n, d$RT8, 0.5),
  RT9 = rnorm(n, d$RT9, 0.5),
  RT10 = rnorm(n, d$RT10, 0.5),
  ID = 1:n,
  time = 3
)

d_final <- rbind(d, d2, d3)

d_final <- reshape_longer(
  d_final,
  select = RT1:RT10,
  names_to = "trial",
  values_to = "ReactionTime"
)

d_median <- rbind(d, d2, d3)
d_median$ReactionTime <- apply(data_select(d_median, select = RT1:RT10), MARGIN = 1, stats::median)

d_final$trialID <- as.numeric(as.factor(d_final$trial))

priors <- data.frame(prior = "gamma(10,1)", class = "ranef")
m1 <- glmmTMB::glmmTMB(
  ReactionTime ~ time + (1 + time | ID),
  data = d_final,
  priors = priors
)
m2 <- glmmTMB::glmmTMB(
  ReactionTime ~ time + (1 + time | ID),
  data = d_median,
  priors = priors
)
m3 <- glmmTMB::glmmTMB(
  ReactionTime ~ time + (1 + time + trialID | ID),
  data = d_final,
  priors = priors
)

model_parameters(m1, ci_random = 0.95)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |        95% CI |     z |      p
#> -----------------------------------------------------------------
#> (Intercept) |        4.76 | 0.15 | [ 4.46, 5.06] | 31.36 | < .001
#> time        |   -6.65e-03 | 0.01 | [-0.03, 0.02] | -0.48 | 0.629 
#> 
#> # Random Effects
#> 
#> Parameter                | Coefficient |       95% CI
#> -----------------------------------------------------
#> SD (Intercept: ID)       |        1.49 | [1.29, 1.72]
#> SD (time: ID)            |        0.04 | [0.02, 0.08]
#> Cor (Intercept~time: ID) |        0.89 | [0.34, 0.96]
#> SD (Residual)            |        0.59 | [0.57, 0.60]
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.

model_parameters(m2, ci_random = 0.95)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |        95% CI |     z |      p
#> -----------------------------------------------------------------
#> (Intercept) |        4.77 | 0.15 | [ 4.48, 5.06] | 31.91 | < .001
#> time        |       -0.01 | 0.01 | [-0.04, 0.01] | -1.20 | 0.230 
#> 
#> # Random Effects
#> 
#> Parameter                | Coefficient |       95% CI
#> -----------------------------------------------------
#> SD (Intercept: ID)       |        1.49 | [1.29, 1.71]
#> SD (time: ID)            |        0.10 | [0.08, 0.12]
#> Cor (Intercept~time: ID) |        0.37 | [0.12, 0.56]
#> SD (Residual)            |        0.10 | [0.08, 0.11]
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.

model_parameters(m3, ci_random = 0.95)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |        95% CI |     z |      p
#> -----------------------------------------------------------------
#> (Intercept) |        4.74 | 0.15 | [ 4.46, 5.03] | 32.65 | < .001
#> time        |   -6.85e-03 | 0.01 | [-0.03, 0.02] | -0.52 | 0.603 
#> 
#> # Random Effects
#> 
#> Parameter                   | Coefficient |        95% CI
#> ---------------------------------------------------------
#> SD (Intercept: ID)          |        1.64 | [ 1.42, 1.89]
#> SD (time: ID)               |        0.04 | [ 0.02, 0.08]
#> SD (trialID: ID)            |        0.06 | [ 0.05, 0.08]
#> Cor (Intercept~time: ID)    |        0.87 | [ 0.19, 0.96]
#> Cor (Intercept~trialID: ID) |       -0.51 | [-0.81, 0.05]
#> Cor (time~trialID: ID)      |       -0.12 | [-0.57, 0.33]
#> SD (Residual)               |        0.56 | [ 0.54, 0.57]
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.


estimate_means(m1, c("time", "ID=c(14,18,79)"))
#> Estimated Marginal Means
#> 
#> time | ID | Mean |   SE |       95% CI |     z
#> ----------------------------------------------
#> 1    | 14 | 4.85 | 0.15 | [4.55, 5.15] | 31.55
#> 2    | 14 | 4.84 | 0.16 | [4.53, 5.15] | 30.88
#> 3    | 14 | 4.84 | 0.16 | [4.52, 5.15] | 30.04
#> 1    | 18 | 3.47 | 0.15 | [3.17, 3.78] | 22.61
#> 2    | 18 | 3.44 | 0.16 | [3.13, 3.74] | 21.92
#> 3    | 18 | 3.40 | 0.16 | [3.08, 3.71] | 21.11
#> 1    | 79 | 5.41 | 0.15 | [5.11, 5.71] | 35.21
#> 2    | 79 | 5.42 | 0.16 | [5.11, 5.73] | 34.57
#> 3    | 79 | 5.43 | 0.16 | [5.11, 5.75] | 33.72
#> 
#> Variable predicted: ReactionTime
#> Predictors modulated: time, ID=c(14,18,79)

estimate_means(m2, c("time", "ID=c(14,18,79)"))
#> Estimated Marginal Means
#> 
#> time | ID | Mean |   SE |       95% CI |     z
#> ----------------------------------------------
#> 1    | 14 | 4.89 | 0.15 | [4.59, 5.19] | 31.97
#> 2    | 14 | 4.79 | 0.16 | [4.49, 5.10] | 30.54
#> 3    | 14 | 4.70 | 0.16 | [4.39, 5.02] | 29.05
#> 1    | 18 | 3.62 | 0.15 | [3.32, 3.92] | 23.71
#> 2    | 18 | 3.60 | 0.16 | [3.29, 3.90] | 22.91
#> 3    | 18 | 3.57 | 0.16 | [3.25, 3.89] | 22.05
#> 1    | 79 | 5.42 | 0.15 | [5.12, 5.72] | 35.46
#> 2    | 79 | 5.49 | 0.16 | [5.18, 5.80] | 34.98
#> 3    | 79 | 5.56 | 0.16 | [5.25, 5.88] | 34.37
#> 
#> Variable predicted: ReactionTime
#> Predictors modulated: time, ID=c(14,18,79)

estimate_means(m3, c("time", "ID=c(14,18,79)"))
#> Estimated Marginal Means
#> 
#> time | ID | Mean |   SE |       95% CI |     z
#> ----------------------------------------------
#> 1    | 14 | 4.85 | 0.15 | [4.56, 5.14] | 32.85
#> 2    | 14 | 4.84 | 0.15 | [4.55, 5.14] | 32.07
#> 3    | 14 | 4.84 | 0.16 | [4.53, 5.14] | 31.11
#> 1    | 18 | 3.47 | 0.15 | [3.18, 3.76] | 23.54
#> 2    | 18 | 3.44 | 0.15 | [3.14, 3.73] | 22.75
#> 3    | 18 | 3.40 | 0.16 | [3.09, 3.70] | 21.85
#> 1    | 79 | 5.41 | 0.15 | [5.12, 5.70] | 36.66
#> 2    | 79 | 5.42 | 0.15 | [5.12, 5.72] | 35.90
#> 3    | 79 | 5.43 | 0.16 | [5.13, 5.74] | 34.94
#> 
#> Variable predicted: ReactionTime
#> Predictors modulated: time, ID=c(14,18,79)
#> Predictors averaged: trialID (5.5)

R Code Explanation: Aggregation vs. Trial-Level Data in GLMMs

Overview

This script performs a simulation study to compare different approaches to analyzing repeated-measures reaction time (RT) data. Specifically, it contrasts analyzing raw, trial-level data (every single reaction time) against analyzing aggregated data (medians per participant per time point).

1. Data Generation (Simulation)

The code creates synthetic data for 100 participants (n <- 100) across 3 distinct time points.

Time Point 1 (d)

• Baseline: A latent "true" speed is generated for each participant (rnorm(n, 5, 1.5)).
• Trials: 10 separate trials (RT1...RT10) are generated for each person. They are centered around that person's baseline with noise (SD = 0.5).

Time Point 2 (d2) & Time Point 3 (d3)

• Dependency: Interestingly, the code generates Time 2 and Time 3 based on the observed values of Time 1 (mean = d$RT1, etc.), not the latent baseline.
• Structure: This creates 10 trials for Time 2 and 10 trials for Time 3.

2. Data Formatting

Two different datasets are prepared for analysis:

1. d_final (Long/Raw Format):
   • The data is reshaped so that every single trial is a row.
   • Dimensions: 100 participants ![][image1] 3 time points ![][image1] 10 trials = 3,000 observations.
   • Includes a trialID variable (1-10) to track which trial number it was within the block.
2. d_median (Aggregated Format):
   • The code calculates the median of the 10 trials for each person at each time point.
   • Dimensions: 100 participants ![][image1] 3 time points = 300 observations.
   • This is a common strategy in cognitive psychology to reduce noise.

3. Modeling (glmmTMB)

The script fits three Linear Mixed Models (LMMs).

Model 1 (m1): Trial-Level Analysis

ReactionTime ~ time + (1 + time | ID)
data = d_final

• Level: Uses all 3,000 observations.
• Structure: Calculates the effect of time on RT, allowing every participant (ID) to have their own starting point (intercept) and their own rate of change over time (slope).
• Result:
  • Residual SD: 0.59. (Reflects the noise between individual trials).
  • SD (time random slope): 0.04.

Model 2 (m2): Aggregated Analysis

ReactionTime ~ time + (1 + time | ID)
data = d_median

• Level: Uses the 300 aggregated medians.
• Result:
  • Residual SD: 0.10. (Much lower, because the trial-to-trial noise has been averaged out).
  • SD (time random slope): 0.10. (This estimate is higher than m1).
  • Implication: Aggregating data makes the model "cleaner" (lower residual error), but can affect the estimation of random effects variance.

Model 3 (m3): Trial-Level with Trial Effects

ReactionTime ~ time + (1 + time + trialID | ID)
data = d_final

• Purpose: Adds trialID as a random slope.
• Meaning: It acknowledges that reaction times might change within a session (e.g., fatigue or practice effects from trial 1 to 10), and that this within-session change varies by participant.
• Result:
  • The SD (trialID: ID) is 0.06, suggesting there is some variation in how people change from trial 1 to 10.

4. Marginal Means (estimate_means)

Finally, the script predicts the specific Reaction Time values for three specific participants (IDs 14, 18, and 79) at times 1, 2, and 3.

• Comparison: The predictions for m1 (raw) and m2 (median) are very similar (e.g., ID 14 at Time 1 is ~4.85 in m1 vs ~4.89 in m2).
• Takeaway: While aggregation changes the variance estimates (SDs), it often yields similar fixed effect estimates (means) if the data is balanced.