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.
These models permit heterogeneous L1 residual variability
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 |
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 |
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 |
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:
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.
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 |
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 |