Specification of crossed random effects model in R

Question

I have an experiment with a design in which subjects answer four items that are of four different types based on two factors (lets call the factors letter: "a" X "b" and big: "A" X "a", resulting in four types of questions A, a, B, b). The order of items (called here 1-4) is held constant and each subject answers one item of each type. The types are randomized. A subject can for example get question-type combinations: 1-a, 2-B, 3-b, 4-A; or 1-B, 2-b, 3-a, 4-A; etc.

I am interested in effects of question types, but expect that the random effects may play a role as well. I tried to use the following model:

glmer(answer ~ (1|subject) + (big*letter|item) + big*letter, data = data, family = binomial(link = "logit"))  

When I compare this model with one without random slopes:

glmer(answer ~ (1|subject) + (1|item) + big*letter, data = data, family = binomial(link = "logit"))

... the first model is not better in any way than the second:

   Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
m2  6 1242.1 1272.1 -615.04   1230.1                        
m1 15 1261.2 1336.2 -615.60   1231.2     0      9          1

So, my first question is whether the model is specified correctly given the design I have. The second question would be, why is it that including random slopes does not improve the model, even though it is possible to see from the data, that the effect of question type obviously differs between the items.

Edit: Summary table for m1:

Generalized linear mixed model fit by maximum likelihood ['glmerMod']
 Family: binomial ( logit )
Formula: answer ~ (1 | subject) + (big * letter | item) + big * letter 
   Data: data 

      AIC       BIC    logLik  deviance 
1261.2010 1336.2061 -615.6005 1231.2010 

Random effects:
 Groups  Name               Variance Std.Dev. Corr             
 subject (Intercept)        0.71862  0.8477                    
 item    (Intercept)        0.00000  0.0000                    
         bigTRUE            0.04241  0.2059     NaN            
         letterTRUE         0.10219  0.3197     NaN  1.00      
         bigTRUE:letterTRUE 0.05749  0.2398     NaN -1.00 -1.00
Number of obs: 1097, groups: subject, 275; item, 4

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          1.8297     0.1798  10.176  < 2e-16 ***
bigTRUE             -0.9339     0.2413  -3.870 0.000109 ***
letterTRUE          -0.7073     0.2734  -2.587 0.009679 ** 
bigTRUE:letterTRUE   0.7458     0.3159   2.361 0.018212 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) bgTRUE ltTRUE
bigTRUE     -0.683              
letterTRUE  -0.602  0.698       
bgTRUE:TRUE  0.521 -0.786 -0.792

Answer 1

You have a within-item design with multiple observations at each level of each factor for each item, and so it is indeed possible in principle to estimate by-item random effects parameters. However, it is important to consider that the estimation algorithm is trying to estimate nine random effects parameters (three slopes and six correlations) based only only four items. It is unlikely that N=4 gives enough information to get stable estimates of these parameters, even having numerous observations for each item in each cell of your design.

If you want to investigate differences among your items then I would suggest including them in your model as fixed effects rather than as random effects.

A few additional points:

The package lme4 has been undergoing rapid development, and I would recommend upgrading to the most recent version on CRAN, which is 1.1-7. There is a helpful blog post by Klinton Bicknell about the various versions of lme4, and 1.1-7 seems to be the best one currently available (especially for fitting logit models).

Also, given that you have a factorial design, I would also recommend centering your predictor variables, e.g.:

data$big.c <- data$big-mean(data$big)

data$letter.c <- data$letter-mean(data$letter)

This will not only give you the typical ANOVA-style interpretation for your main effects, but can sometimes help with convergence.

Finally, the blog post I mentioned suggests using control=glmerControl(optimizer='bobyqa') as an argument to lmer for best performance with categorical data.