library(lme4) library(RePsychLing) library(knitr) opts_chunk$set(comment=NA) options(width=92,show.signif.stars = FALSE)
Some of the data from @barrseyfedd2010 are available as the data frame bs10
in the RePsychLing
package.
str(bs10)
As with other data frames in this package, the subject and item factors are called subj
and item
. The response being modelled, dif
, is the difference in two response times.
The two experimental factors S
and F
, both at two levels, are represented in the -1/+1 encoding, as is their interaction, SF
. The S
factor is the speaker condition with levels -1 for the same speaker in both trials and +1 for different speakers. The F
factor is the filler condition with levels -1 for NS
and +1 for FP
.
The maximal model has a full factorial design 1+S+F+SF
for the fixed-effects and for potentially correlated vector-valued random effects for subj
and for item
. We use the parameter estimates from an lmm
fit using the MixedModels package package for Julia, which can be much faster than fitting this model with lmer
. (In addition to being faster, the fit from lmm
produced a significantly lower deviance.)
m0 <- lmer(dif ~ 1+S+F+SF + (1+S+F+SF|subj) + (1+S+F+SF|item), bs10, REML=FALSE, start=thcvg$bs10$m0, control=lmerControl(optimizer="Nelder_Mead",optCtrl=list(maxfun=1L), check.conv.grad="ignore",check.conv.hess="ignore")) summary(m0)
The model converges with warnings. Six of 12 correlation parameters are estimated at the +/-1 boundary.
Notice that there are only 12 items. Expecting to estimate 4 variances and 6 covariances from 12 items is optimistic.
A summary of a principal components analysis (PCA) of the random-effects variance-covariance matrices shows
summary(rePCA(m0))
For both the by-subject and the by-item random effects the estimated variance-covariance matrices are singular. There are at least 2 dimensions with no variation in the subject-related random-effects and 3 directions with no variation for the item-related random-effects.
Clearly the model is over-specified.
A zero-correlation-parameter model fits independent random effects for the intercept, the experimental factors and their interaction for each of the subj
and item
grouping factors.
It can be conveniently specified using the ||
operator in the random-effects terms.
print(summary(m1 <- lmer(dif ~ 1+S+F+SF + (1+S+F+SF||subj) + (1+S+F+SF||item), bs10, REML=FALSE)), corr=FALSE) anova(m1, m0)
The zcpLMM fits significantly worse than the maxLMM. However, the results strongly suggest that quite a few of the variance components are not supported by the data.
summary(rePCA(m1))
The random effect for filler, F
, by subject has essentially zero variance and the random effect for the interaction, SF
, accounts for less than 15% of the total variation in the random effects. There is no evidence for item-related random effects in m1
.
Remove all variance components estimated with a value of zero.
print(summary(m2 <- lmer(dif ~ 1+S+F+SF + (1+S+SF||subj), bs10, REML=FALSE)), corr=FALSE)
Naturally, the fit for this model is equivalent to that for m1
because it is only the variance components with zero estimates that are eliminated.
anova(m2,m1)
Next we check whether the variance component for the interaction, SF
, could reasonably be zero.
print(summary(m3 <- lmer(dif ~ 1+S+F+SF + (1+S||subj), bs10, REML=FALSE)), corr=FALSE) anova(m3, m2)
Not quite significant, but could be considered. The fit is still worse than for the maxLMM m0
. We now reintroduce a correlation parameters in the vector-valued random effects for subj
.
print(summary(m4 <- lmer(dif ~ 1+S+F+SF + (1+S|subj), bs10, REML=FALSE)), corr=FALSE) anova(m3, m4, m0)
Looks like we have evidence for a significant correlation parameter. Moreover, LMM m4
fits as well as maxLMM.
LMM m4
is the optimal model. It might be worth while to check the theoretical contribution of the correlation parameter.
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.