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)
