O'Reilly logo

The R Book by Michael J. Crawley

Stay ahead with the world's most comprehensive technology and business learning platform.

With Safari, you learn the way you learn best. Get unlimited access to videos, live online training, learning paths, books, tutorials, and more.

Start Free Trial

No credit card required

Error Plots from a Hierarchical Analysis

You will find the syntax of model checking in lme both complicated and difficult to remember.

library(nlme)
library(lattice)
trellis.par.set(col.whitebg())

If you use the standard plot(model) with lme you get a single panel showing the standardized residuals as a function of the fitted values. For more comprehensive model checking, it is useful to make the dataframe into a groupedData object, then refit the model. Here we investigate the REML model with gender as a fixed effect:

hs<-groupedData(subject~gender|town/district/street/family/gender/replicate,
                                       outer=~gender,data=hierarchy)
model<-
lme(subject~gender,random=~1|town/district/street/family/gender,data=hs)
plot(model,gender~resid(.))

images

To inspect the constancy of variance across towns and check for heteroscedasticity:

plot(model,resid(.,type="p")~fitted(.)|town)

It should be clear that this kind of plot only makes sense for those variables with informative factor levels such as gender and town; it would make no sense to group together the streets labelled s1 or s3 or the families labelled f1, f2 or f3.

Tests for normality use the familiar QQ plots, but applied to panels:

qqnorm(model,~resid(.)|gender)

images

The residuals are normally distributed for both genders, and within each town: ...

With Safari, you learn the way you learn best. Get unlimited access to videos, live online training, learning paths, books, interactive tutorials, and more.

Start Free Trial

No credit card required