[R-sig-ME] residuals from glmm using lmer

[Ben Bolker]

martin highfield <mhighfield106 at ...> writes:

> 
> Dear all
> 
> I am trying to validate a GLMM using plots of residuals versus fitted values
> and other explanatory variables. However, I am unsure what type of residuals
> (or fitted values) are returned when using the residuals function on a lmer
> object. My data and model are similar to the example given below  where I
> have a binary response variable (y) , a continous explanatory (day) ,a two
> level nominal explanatory variable (treat) and a random effect of group ( 5
> groups).
> 
> set.seed(100)
> y <- rbinom(n=100, size=1, prob=0.5)
> day <- sample(1:10, 100, replace=T)
> group <- rep(paste("group", 1:5, sep=""), times=20)
> treat <- sample(c("T1", "T2"), 100, replace=T)
> 
> dataf <- data.frame(y,day,treat, group)
> library(lme4)
> mod.lmer <- lmer(y~day+treat+(1|group), family=binomial, data=dataf)
> summary(mod.lmer)
> 
> I can extract the residuals and fitted values using the residuals and fitted
> extractor functions and plot them.
> 
> res <- residuals(mod.lmer)
> fit <- fitted(mod.lmer)
> 
> plot(fit,res)
> abline(h=0)
> plot(res~dataf$treat)
> plot(res~dataf$day)
> 
> However, what do these residuals represent? Are they simply the observed
> minus the fitted values, Pearsons residuals (y_i - mu_i /sqrt(mu_i)),
> deviance residuals etc? 

  You can convince yourself that they're the Pearson residuals
(which are NOT what you've quoted here -- you've quoted the Pearson
residuals for the Poisson case [and I think you're missing parentheses?])

head(res)
head(y-fit)  ## not the same as res
est.sd <- sqrt(fit*(1-fit))
head((y-fit)/est.sd)   ## the same as res


> Also, I would like to extract both the residuals
> conditional on the random effect and also the residuals based on both the
> fixed and random effects but am unsure how to do this. Any
> advice/clarification would be gratefully appreciated.

  fitted() gives you the unconditional (fixed+random) estimates.
  
  Don't have time to work it out now, but you can construct your
own fixed- and random-effect predictions by constructing the appropriate
model matrices and multiplying by the random effects (ranef())
and fixed effect parameter (fixef()) vectors ...

  By the way, since you didn't add any group-level variation
to your example above, all the random effects are zero ...