Assessing the fit of GLMM implementation of a Rasch model to binary data using lme4

Question

I'd like to assess the fit of the kinds of models described by de Boeck et al (2011) (http://www.jstatsoft.org/v39/i12). They are GLMM implementations of Rasch family models, e.g.:

all.1.NULL.Rasch <- glmer(score ~ -1 
                          + item
                          + (-1 + X064.part.Part | person), data = IN11c, family=binomial)

As Gelman and Hill (2006) explain, a binned residual plot is a better approach than simply plotting the fitted values against the residuals (also discussed here: How to assess the fit of a binomial GLMM fitted with lme4 (> 1.0)?).

However, my plots (using the arm package) exhibit a striking pattern, which seems to me likely to be due to some aspect of the model which is different from those Gelman and Hill were talking about.

Can anyone explain this, and, preferably also a useful way to assess the fit of the type of model I have?

First, I'm assuming that this is correct:

binnedplot(fitted(all.1.NULL.Rasch),residuals(all.1.NULL.Rasch))

The data actually has 89 variables, so I'm supplying details for those in the quoted model:

> str(IN11c$score)
 Ord.factor w/ 2 levels "0"<"1": 2 2 2 2 2 2 2 2 2 2 ...

> str(IN11c$X064.part.Part)
 Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

> str(IN11c$person)
 Factor w/ 9961 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...

> summary(IN11c$score)
     0      1   NA's 
109222 238188   1225 

> summary(IN11c$X064.part.Part)
     1      2      3      4 
 69727  79688  69727 129493 

> summary(IN11c$person)
      1       2       3       4       5       6       7       8       9      10      11      12      13      14      15 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 16      17      18      19      20      21      22      23      24      25      26      27      28      29      30 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 31      32      33      34      35      36      37      38      39      40      41      42      43      44      45 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 46      47      48      49      50      51      52      53      54      55      56      57      58      59      60 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 61      62      63      64      65      66      67      68      69      70      71      72      73      74      75 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 76      77      78      79      80      81      82      83      84      85      86      87      88      89      90 
 35      35      35      35      35      35      35      35      35      35      35      35      35      35      35 
 91      92      93      94      95      96      97      98      99 (Other) 
 35      35      35      35      35      35      35      35      35  345170

and the output:

> summary(all.1.NULL.Rasch)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)  ['glmerMod']

Family: binomial  ( logit )

Formula: score ~ -1 + item + (-1 + X064.part.Part | person)
Data: IN11c

  AIC       BIC    logLik  deviance  df.resid 
369799.0  370283.1 -184854.5  369709.0    347365 

Scaled residuals: 
Min      1Q  Median      3Q     Max 
-6.4212 -0.6800  0.3506  0.5897  6.6114 

Random effects:
Groups Name            Variance Std.Dev. Corr          
person X064.part.Part1 2.2550   1.5017                 
    X064.part.Part2 0.7472   0.8644   0.71          
    X064.part.Part3 2.4486   1.5648   0.54 0.75     
    X064.part.Part4 0.7124   0.8441   0.63 0.86 0.69
Number of obs: 347410, groups:  person, 9961

Fixed effects:
    Estimate Std. Error z value Pr(>|z|)    
itemX1   3.01316    0.04289   70.25  < 2e-16 ***  
itemX2  -0.60039    0.02900  -20.71  < 2e-16 ***  
itemX3   2.06994    0.03466   59.72  < 2e-16 ***  
itemX4   0.05922    0.02830    2.09  0.03635 *  
itemX5   2.69565    0.03966   67.96  < 2e-16 ***  
itemX6   1.67934    0.03239   51.85  < 2e-16 ***  
itemX7   2.27617    0.03612   63.02  < 2e-16 ***  
itemX8   1.61074    0.02833   56.86  < 2e-16 ***  
itemX9   0.60380    0.02391   25.26  < 2e-16 ***  
itemX10  0.46635    0.02362   19.74  < 2e-16 ***  
itemX11  1.60557    0.02830   56.73  < 2e-16 ***  
itemX12  0.63788    0.02398   26.60  < 2e-16 ***  
itemX13  0.87095    0.02463   35.36  < 2e-16 ***  
itemX14  1.20185    0.02602   46.18  < 2e-16 ***  
itemX15  0.92694    0.02484   37.32  < 2e-16 ***  
itemX16  0.30951    0.02872   10.78  < 2e-16 ***  
itemX17 -0.08539    0.02869   -2.98  0.00292 **  
itemX18  0.13537    0.02866    4.72 2.32e-06 ***  
itemX19  1.34097    0.03091   43.39  < 2e-16 ***  
itemX20  1.03142    0.03000   34.38  < 2e-16 ***  
itemX21  0.91493    0.02967   30.84  < 2e-16 ***  
itemX22 -0.14124    0.02874   -4.91 8.94e-07 ***  
itemX23 -0.36678    0.02346  -15.63  < 2e-16 ***  
itemX24  1.89063    0.03027   62.45  < 2e-16 ***  
itemX25  2.38406    0.03510   67.93  < 2e-16 ***  
itemX26  1.06347    0.02527   42.08  < 2e-16 ***  
itemX27  0.59852    0.02378   25.17  < 2e-16 ***  
itemX28  0.38406    0.02341   16.40  < 2e-16 ***  
itemX29  1.43106    0.02709   52.83  < 2e-16 ***  
itemX30  1.22191    0.02600   46.99  < 2e-16 ***  
itemX31  0.97085    0.02489   39.00  < 2e-16 ***  
itemX32  1.16746    0.02568   45.46  < 2e-16 ***  
itemX33 -0.06049    0.02322   -2.60  0.00919 **  
itemX34  0.44368    0.02351   18.87  < 2e-16 ***  
itemX35  1.65890    0.02853   58.15  < 2e-16 ***  

Answer 1

OK, the only thing I can think of (in the absence of any continuous predictors) is link misspecification. Tried this using the logit-power link code from the linked question, but I can't actually get it to look like your data (it's much noisier).

Logit-power (uses (fragile) internal C calls for speed; could use plogis(), qlogis() for readability and stability instead)

logitpower <- function(lambda) {
    L <- list(linkfun=function(mu)
              .Call(stats:::C_logit_link,mu^(1/lambda),PACKAGE="stats"),
              linkinv=function(eta)
              .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")^lambda,
              mu.eta=function(eta) {
                  mu <-  .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")
                  mu.eta <-  .Call(stats:::C_logit_mu_eta,eta,PACKAGE="stats")
                  lambda*mu^(lambda-1)*mu.eta
              },
              valideta = function(eta) TRUE ,
              name=paste0("logit-power(",lambda,")"))
    class(L) <- "link-glm"
    L
}

Comparing the shape of logit-power links, tweaked to have approximately the same x-intercept:

par(las=1,bty="l")
curve(logitpower(1)$linkinv(x),from=-4,to=4)
curve(logitpower(0.5)$linkinv(x-1),add=TRUE,col=2)
curve(logitpower(2)$linkinv(x+0.8),add=TRUE,col=4)
legend("bottomright",bty="n",lty=1,col=c(1,2,4),
       paste0("lambda=",c(1,0.5,2)))
abline(v=0,lty=2)

This suggests that the pattern observed above (underpredict for values < 0.5, overpredict for values > 0.5) could come from a logit-power link with $\lambda>1$

However, when we try simulating and fitting with regular logit and logit-power links, it doesn't reproduce your pattern (and doesn't seem to make much difference):

ff <- binomial(link=logitpower(4))
dd <- expand.grid(f1=factor(1:10),f2=factor(1:20),rep=1:20)
library("lme4")
betavec <- seq(-2,2,length=10)
dd$y <- simulate(~-1+f1+(1|f2),family=ff,
         newdata=dd,
         newparams=list(theta=1,beta=betavec),
         seed=101)[[1]]
fit1 <- glmer(y~-1+f1+(1|f2),family=ff,
              data=dd)
fit2 <- update(fit1,family=binomial)

Binned residual plots:

library("arm")
par(mfrow=c(1,2))
binnedplot(predict(fit1,type="response"),resid(fit1),nclass=100)
binnedplot(predict(fit2,type="response"),resid(fit2),nclass=100)