Skip to content

plotting/examining residuals of a mixed generalised linear model

3 messages · kath_oreilly, Ben Bolker

Original
kath_oreilly
# 
Dear R users,

I'm hoping that more experienced users will be able to assist me in
examining the model fit of a mixed generalised linear model. The example
using the data 'bacteria' within the MASS package will hopefully illustrate
what I would like to acheive;

library(MASS)
library(nlme)
attach(bacteria) # y being output and the trt - treatment group being an
explanatory variable. There is pseudoreplication as each patient (ID) is
sampled multiple times (week)

bacteria$y<-1*(bacteria$y=="y")  # to make readable in lmer
table(bacteria$y,bacteria$trt)
hs <- groupedData(y~trt|ID,outer=~trt,data=bacteria)  # I don't think this
is really necessary
model <- lmer(y ~ trt + (week|ID),family=binomial,data=hs)
summary(model)

Here I would like to examine the fit of the variable trt by examining the
residuals. In lm (using lme in "the R book, p. 657"), one would be able to
use

plot(model,trt~resid(.))

However it doesn't work.

If some one would explain why, that would be great. I've come across the
package "zelig", which uses simulation to examine model fit. I haven't
gotten my head round this yet, and was hoping some one would advise the
approp path to take.

Many thanks!
Ben Bolker
#
kath_oreilly wrote:
I think you're a little confused. Thanks for the reproducible example,
though.

First of all, you're confusing nlme ("old mixed models", described most
thoroughly
in Pinheiro and Bates 2000, good for correlation and heteroscedasticity
models
but not capable of doing GLMMs) with lme4 ("new mixed models", not yet
thoroughly documented but D. Bates is working on a book ...) Mixing
them is usually a bad idea.

The plotting capabilities aren't as highly developed for lme4 as for
nlme, but you can still get what you want done.

The r-sig-mixed-models list may be better for this kind of question.

Here's my solution to your problem:

library(MASS)
data(bacteria)

with(bacteria,table(y,trt))
library(lme4)
bacteria$pres <- as.numeric(bacteria$y=="y")  ## derive 0/1 variable
model <- lmer(pres ~ trt + (week|ID),family=binomial,data=bacteria)
summary(model)

with(bacteria,boxplot(residuals(model)~trt))

## or boxplot(residuals(model)~bacteria$trt)