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)