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)