This file contains final MLMs for reaction times on the flanker task for the PD Inhibition DDM project.

Note: I have adopted the common convention of only analyzing RTs for correct responses. This assumes that whatever process gives rise to the RT for an incorrect response is not commensurable. The alternative would be to allow for correct to enter into the model as a predict of the current trial’s RT. This could support moderation analyses – for example, whether error RTs are faster in the mostly congruent condition for incongruent trials. The downside is all of the ensuing multi-way interactions, many of which are not of central interest.

N.B. For the AIC comparisons to be valid, we need to use method=‘ML’ throughout. Once we have a winning model, or one where the parameters matter, re-estimate with REML. update(mm, method='REML')

To reduce the number of models, I don’t consider block diagonal matrices for L2 random effects. This is achieved with pdDiag() and basically estimates a unique variance component for each level of a factor. But the fits of these models are extremely close to the more conventional block within id id/block parameterization for random effects. Given that the latter is much more conventional, let’s stick with that when we want to check for heterogeneity in L2 by condition.

LME approach

These models permit heterogeneous L1 residual variability

m1: design effects, random intercept of subject

minfo[["m1"]] <- c(fixed="Congruency, Block", l1="Homogeneous", l2="Subject: Intercept")
m1 <- lme(rt_inv ~ cond*block, 
          random = ~ 1 | id,
          na.action = na.exclude, 
          data = flanker, method='ML')

m1_reml <- update(m1, method='REML')
flextable(car::Anova(m1, type=3)) %>% autofit()

Chisq

Df

Pr(>Chisq)

8196.32

1

0.00e+00

162.18

1

3.78e-37

7.38

1

6.58e-03

14.82

1

1.18e-04

#Get model-predicted means in original RT scale. Force back onto original RT using custom function
rg <- update(ref_grid(m1_reml), tran=list(linkinv=function(x) { -1000/x }, mu.eta=function(x) { -1000/x }), predict.type="response")
res_trans <- emmeans(rg, ~cond|block)

flextable(as.data.frame(res_trans)) %>% colformat_num(j=c("response", "SE", "lower.CL", "upper.CL"), digits=2) %>% autofit()

cond

block

response

SE

df

lower.CL

upper.CL

congruent

most_incon

370.13

11.09

106

362.17

378.46

incongruent

most_incon

386.56

11.38

106

378.02

395.48

congruent

most_con

373.50

11.00

106

365.54

381.83

incongruent

most_con

397.91

11.95

106

388.70

407.56

m2: design effects, allow for L2 random intercept to vary by block

Use nested RE specification.

minfo[["m2"]] <- c(fixed="Congruency, Block", l1="Homogeneous", l2="Block Within Subject: Intercept")
m2 <- lme(rt_inv ~ cond*block,
          random=~1|id/block,
          na.action = na.exclude,
          data = flanker, method='ML')
flextb_vif(m2)

Predictor

vif

cond

2.08

block

1.39

cond:block

2.09

flextable(car::Anova(m2, type=3)) %>% autofit()

Chisq

Df

Pr(>Chisq)

7997.73

1

0.00e+00

166.21

1

4.98e-38

3.62

1

5.71e-02

14.51

1

1.39e-04

m3: design effects, allow for L1 residual variance to vary by block

minfo[["m3"]] <- c(fixed="Congruency, Block", l1="Per Block", l2="Block within subject: Intercept")
m3 <- lme(rt_inv ~ cond*block,
          random=~1|id/block,
          weights=varIdent(form=~1|block),
          na.action = na.exclude,
          data = flanker, method='ML')
summary(m3)
## Linear mixed-effects model fit by maximum likelihood
##  Data: flanker 
##     AIC   BIC logLik
##   14536 14597  -7260
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)
## StdDev:       0.295
## 
##  Formula: ~1 | block %in% id
##         (Intercept) Residual
## StdDev:      0.0666    0.358
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | block 
##  Parameter estimates:
## most_incon   most_con 
##       1.00       1.07 
## Fixed effects: rt_inv ~ cond * block 
##                                Value Std.Error    DF t-value p-value
## (Intercept)                   -2.702   0.03015 16039   -89.6  0.0000
## condincongruent                0.115   0.00864 16039    13.4  0.0000
## blockmost_con                  0.024   0.01265   106     1.9  0.0581
## condincongruent:blockmost_con  0.049   0.01275 16039     3.8  0.0001
##  Correlation: 
##                               (Intr) cndncn blckm_
## condincongruent               -0.198              
## blockmost_con                 -0.244  0.471       
## condincongruent:blockmost_con  0.134 -0.677 -0.478
## 
## Standardized Within-Group Residuals:
##     Min      Q1     Med      Q3     Max 
## -4.4977 -0.5985  0.0157  0.5960  4.1934 
## 
## Number of Observations: 16255
## Number of Groups: 
##            id block %in% id 
##           107           214
#stargazer(m3, type="html")
flextable(car::Anova(m3, type=3)) %>% autofit()

Chisq

Df

Pr(>Chisq)

8031.10

1

0.00e+00

178.53

1

1.01e-40

3.67

1

5.54e-02

14.49

1

1.41e-04

flextable(render_aictab(minfo=minfo, m1, m2, m3)) %>% autofit()

model

fixed

l1

l2

logLik

AIC

dAIC

df

weight

m1

Congruency, Block

Homogeneous

Subject: Intercept

-7343

14698

162.6

6

4.98e-36

m2

Congruency, Block

Homogeneous

Block Within Subject: Intercept

-7279

14571

35.9

7

1.59e-08

m3

Congruency, Block

Per Block

Block within subject: Intercept

-7260

14536

0.0

8

1.00e+00

m4: Allow L1 residual variability to depend on both condition and block

minfo[["m4"]] <- c(fixed="Congruency, Block", l1="Per Block and Condition", l2="Block within subject: Intercept")
m4 <- lme(rt_inv ~ cond*block,
          random=~1|id/block,
          weights=varIdent(form=~1|block*cond),
          na.action = na.exclude,
          data = flanker, method='ML')
summary(m4)
## Linear mixed-effects model fit by maximum likelihood
##  Data: flanker 
##     AIC   BIC logLik
##   14530 14607  -7255
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)
## StdDev:       0.296
## 
##  Formula: ~1 | block %in% id
##         (Intercept) Residual
## StdDev:      0.0667    0.361
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | block * cond 
##  Parameter estimates:
## most_incon*incongruent   most_incon*congruent     most_con*congruent   most_con*incongruent 
##                  1.000                  0.982                  1.049                  1.104 
## Fixed effects: rt_inv ~ cond * block 
##                                Value Std.Error    DF t-value p-value
## (Intercept)                   -2.702   0.03014 16039   -89.6  0.0000
## condincongruent                0.115   0.00858 16039    13.5  0.0000
## blockmost_con                  0.024   0.01258   106     1.9  0.0569
## condincongruent:blockmost_con  0.048   0.01286 16039     3.8  0.0002
##  Correlation: 
##                               (Intr) cndncn blckm_
## condincongruent               -0.194              
## blockmost_con                 -0.242  0.465       
## condincongruent:blockmost_con  0.130 -0.667 -0.464
## 
## Standardized Within-Group Residuals:
##     Min      Q1     Med      Q3     Max 
## -4.4718 -0.6017  0.0136  0.5967  4.1736 
## 
## Number of Observations: 16255
## Number of Groups: 
##            id block %in% id 
##           107           214
flextable(render_aictab(minfo=minfo, m1, m2, m3, m4)) %>% autofit()

model

fixed

l1

l2

logLik

AIC

dAIC

df

weight

m1

Congruency, Block

Homogeneous

Subject: Intercept

-7343

14698

168.26

6

2.74e-37

m2

Congruency, Block

Homogeneous

Block Within Subject: Intercept

-7279

14571

41.59

7

8.78e-10

m3

Congruency, Block

Per Block

Block within subject: Intercept

-7260

14536

5.68

8

5.51e-02

m4

Congruency, Block

Per Block and Condition

Block within subject: Intercept

-7255

14530

0.00

10

9.45e-01

Interim conclusions about random effects structure:

  1. Allow for per-block random intercept (nested random effect)
  2. Allow for L1 residual variance to vary by block and trial condition

But!! Obtaining p-values and denominator degrees of freedom is a somewhat tortured process because neither Satterthwaite nor Kenward-Roger corrections for degrees of freedom are provided. Moreover, we don’t have any conceptual interest in L1 heterogeneity. Likewise, when we switch to glmer for accuracy, we no longer have the option for heterogeneous L1 variances since nlme doesn’t provide a binary logistic model as far as I know.

Thus, switch to lmer() since the best-fitting model (m30) depends primarily on nested L2 random effects. This makes it easier to switch to Mplus, too.

LMER-based modeling

m1: design effects, random intercept of subject

minfo <- list() #re-initialize now that we're eliminating l1 heterogeneity
minfo[["m1"]] <- c(fixed="Congruency, Block", l2="Subject: Intercept")
m1 <- lmer(rt_inv ~ cond*block + (1|id),
          na.action = na.exclude, 
          data = flanker, REML=FALSE)
m1_reml <- update(m1, REML=TRUE)
flextable(anova(m1) %>% rownames_to_column("Effect")) %>% autofit() #lmerTest reports Satterthwaite degrees of freedom

Effect

Sum Sq

Mean Sq

NumDF

DenDF

F value

Pr(>F)

cond

66.10

66.10

1

16149

471.4

4.70e-103

block

8.21

8.21

1

16148

58.5

2.12e-14

cond:block

2.08

2.08

1

16148

14.8

1.19e-04

#flextable(car::Anova(m1, type=3)) %>% autofit()

#Get model-predicted means in original RT scale. Force back onto original RT using custom function
rg <- update(ref_grid(m1_reml), tran=list(linkinv=function(x) { -1000/x }, mu.eta=function(x) { -1000/x }), predict.type="response")
res_trans <- emmeans(rg, ~cond|block, lmer.df = "satterthwaite")

# emmeans(m1_reml, ~cond|block, lmer.df = "satterthwaite") #in inverse units

flextable(as.data.frame(res_trans)) %>% 
  colformat_num(j=c("response", "SE", "lower.CL", "upper.CL"), digits=2) %>% 
  autofit()

cond

block

response

SE

df

lower.CL

upper.CL

congruent

most_incon

370.13

11.09

118

362.18

378.45

incongruent

most_incon

386.56

11.38

110

378.03

395.48

congruent

most_con

373.50

11.00

110

365.54

381.82

incongruent

most_con

397.91

11.95

119

388.71

407.55

m2: design effects, allow for L2 random intercept to vary by block

Use nested RE specification.

minfo[["m2"]] <- c(fixed="Congruency, Block", l2="Block within subject: Intercept")
m2 <- lmer(rt_inv ~ cond*block + (1|id/block),
          na.action = na.exclude,
          data = flanker, REML=FALSE)
m2_reml <- update(m2, REML=TRUE)
flextable(anova(m2) %>% rownames_to_column("Effect")) %>% autofit()

Effect

Sum Sq

Mean Sq

NumDF

DenDF

F value

Pr(>F)

cond

66.21

66.21

1

16052

479.9

7.75e-105

block

2.63

2.63

1

120

19.0

2.73e-05

cond:block

2.00

2.00

1

16073

14.5

1.40e-04

#car::Anova(m2, type=3)
flextable(render_aictab(minfo=minfo, m1, m2)) %>% autofit()

model

fixed

l2