The binned residual plot in the R arm package is often recommended as a way to check if a logistic model is making any systematic errors. The general idea is that the mean residual for a group of observations with similar fitted values should be close to zero. Or equivalently, for a group of observation with mean fitted value p, the proportion of the group for which the response = 1 should be roughly p.
When calculating fitted values for a GLMM, the lme4 package provides the option of either including or excluding the random effects. For the purpose of a binned residual plot, I would have thought random effects should be included. However, doing so produces a plot that shows a clear sinusoidal pattern, whereas the plot with random effects excluded looks much better. The same pattern is visible for all datasets I've checked; two reproducible examples are included below.
In both of the examples below, a LRT shows that the addition of random effects significantly improves the models relative to equivalent GLMs. Given this is the case, why would including random effects in the fitted value show that the model is making a systematic error?
Examples:
Note that both fitted values and residuals are on the response scale.
Using the verbal aggression data, which is included in the lme4 package
library(lme4)
library(arm)
data(VerbAgg, package = 'lme4')
verb_mod <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 + (1|id) + (1|item), family = binomial, data = VerbAgg)
par(mfcol=c(1, 2))
binnedplot(predict(verb_mod, type="response", re.form=NULL), resid(verb_mod, type="response"), nclass=40, main='With random effects')
binnedplot(predict(verb_mod, type="response", re.form=NA), resid(verb_mod, type="response"), nclass=40, main='Without random effects')
Using the bdf language score data, which is included in the nlme package
library(nlme)
data(bdf, package = "nlme")
bdf <- subset(bdf, select = c(schoolNR, Minority, ses, repeatgr))
bdf$repeatgr[bdf$repeatgr == 2] <- 1
bdf_mod <- glmer(repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR), data = bdf, family = binomial(link = "logit"))
par(mfcol=c(1, 2))
binnedplot(predict(bdf_mod, type="response", re.form=NULL), resid(bdf_mod, type="response"), main='With random effects', nclass=20)
binnedplot(predict(bdf_mod, type="response", re.form=NA), resid(bdf_mod, type="response"), main='Without random effects', nclass=20)
Link to bdf tutorial http://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
