[R] Residuals for binomial lmer fits

[Andy Fugard]

Dear all,

I would like to use the residuals in a general linear mixed effect model 
to diagnose model fit.

I know that the resid function has been implemented for linear mixed 
models but not yet for general linear mixed effects.  Is there a way to 
get them out of lmer fit objects?

I tried searching the r-help archive and found nothing.

I tried and failed to replicate what (I guessed would be equivalent to 
what) resid does, starting with one predictor and a random intercept. 
The residuals output by resid had the following properties:

 > summary(resid(lmer1))
      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
-5.92e+01 -1.29e+01 -1.10e+00  2.62e-15  1.31e+01  6.66e+01
 > sd(resid(lmer1))
[1] 19.6

My naive method gave:

 > summary(my.resids)
      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
-1.84e+02 -3.81e+01  3.12e+00 -1.93e-13  4.00e+01  1.71e+02
 > sd(my.resids)
[1] 56.4

Both the mean and the sd are MUCH higher.

All I did was use the fixed effects with random effects flattened to 
predict y's.  Code below.

Best wishes,

Andy


##############################################################

require(lme4)

set.seed(0)

# Data for 100 people, each of whom has 10 observations

people = 100
obs = 10

# Generate person-level variation

sub_noise = data.frame(id = 1:people, rand_int = rnorm(people,0,60))

# And within person noise

resids = rnorm(people*obs, 0, 20)

id = rep(1:people,obs)

thedata = data.frame(id)
thedata = merge(sub_noise, thedata)
thedata$x = rep(1:obs,people)
thedata$y = 23*thedata$x + 20 + thedata$rand_int + resids
thedata$y.flat = 23*thedata$x + 20 + resids

# Fit a random intercept model

lmer1 = lmer(y ~ x + (1|id), data=thedata)
summary(lmer1)

# Get the intercepts

ranef(lmer1)$id[1]

ran_effects = data.frame(rownames(ranef(lmer1)$id), ranef(lmer1)$id[1])
names(ran_effects) = c("id", "b")
ran_effects_data = merge(thedata, ran_effects)

# Calculate the predicted y's using the fixed effects and flattening out
# using the random effects:

predicted.y = fixef(lmer1)[1] + ran_effects_data$x * fixef(lmer1)[2]
                 + ran_effects_data$b

# Now how far off were we?

my.resids = predicted.y - ran_effects_data$y

summary(resid(lmer1))
sd(resid(lmer1))
summary(my.resids)
sd(my.resids)