I'm trying to fit a GLMM on simulated response time data (continuous, positively skewed) obtained from hypothetical participants (Subject) responding to the same set of hypothetical items (Item), so it's a fully-crossed design. I intend to include several crossed-random effects for Subject and Item, so in lme4 language, it would look like the following:
glmer(y ~ x1*x2*z1 + (1+x1+x2|Subject) + (1|Item), family=Gamma("identity"), data=foo)
However, as I read from Ben Bolker's GLMM FAQ (https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#fn1), the estimation procedure used by glmer (adaptive Gauss-Hermite quadrature) can only handle up to 2-3 random effects. Indeed, running glmer on my simulated data not only results in inevitable non-convergence but also takes such a long time to run.
Someone recommended to me to use MASS::glmmPQL instead because the cases in which penalized quasi-likelihood (PQL) would perform poorly (count/binomial DV, mean DV < 5) doesn't apply to my data (continuous DV, identity link, many items, and many subjects). Moreover, PQL could handle more random effects than GHQ; it could also allow for correlations of random effects to be estimated; and it estimates the model faster than GHQ. (I don't actually know about any of those being accurate characterizations of PQL and GHQ; would be happy to be corrected and pointed to the right direction.)
I initially thought that the response in this thread (Specifying multiple (separate) random effects in lme) gave the answer, but it models the random effect for Item as if it was nested under Subject, but they're crossed (note Item %in% Subject in results of a test model (i.e., doesn't have the random effects I want as in the model I mentioned earlier) below):
> bar <- glmmPQL(y ~ x1*x2*z1, random=list(Subject=~1, Item=~1),
family=Gamma("identity"),data=foo)
> summary(bar)
Linear mixed-effects model fit by maximum likelihood
Data: foo
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 88.78614
Formula: ~1 | Item %in% Subject
(Intercept) Residual
StdDev: 197.4256 0.0005360151
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: y ~ x1 * x2 * z1
Value Std.Error DF t-value p-value
(Intercept) 515.1268 13.382142 3735 38.49360 0.0000
x1 -27.4088 4.163655 3735 -6.58287 0.0000
x2 -18.1321 6.422738 3735 -2.82312 0.0048
z1 -26.6217 11.966715 48 -2.22464 0.0308
x1:x2 2.7502 5.884197 3735 0.46739 0.6402
x1:z1 4.0648 3.657627 3735 1.11132 0.2665
x2:z1 2.2683 5.685076 3735 0.39899 0.6899
x1:x2:z1 1.3081 5.183023 3735 0.25238 0.8008
*snip results*
Number of Observations: 3791
Number of Groups:
Subject Item %in% Subject
50 3791
I was wondering if someone can point me to how I can specify the model such that the random effects of Subject and Item are truly crossed. Thanks!
