# [R] Residuals for binomial lmer fits

Andy Fugard a.fugard at ed.ac.uk
Mon Oct 8 13:19:28 CEST 2007

```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))
 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)
 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

ran_effects = data.frame(rownames(ranef(lmer1)\$id), ranef(lmer1)\$id)
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) + ran_effects_data\$x * fixef(lmer1)
+ 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)

```