19.5 Hierarchical sampling and variance components analysis

Hierarchical data are often encountered in observational studies where information is collected at a range of

different spatial scales. 

The principal aim is to discover the scale at which most of the variation is generated.

This information would then allow a closer focus on mechanisms operating at this scale in subsequent more

detailed studies

The following study involves a test with a mean score of 100 administered to children

in four British towns. 

Each town was divided into districts by postcodes, and 

six districts were selected at random. 

Within districts, 10 streets were selected at random, and 

within streets, four households were

Naturally, different households had different numbers of children (childless households

were excluded from the study) and there was no control over sex ratio of children within household

library(lme4)
data <- read.table("c:\\temp\\childfull.txt",header=T)
attach(data)

head(data)

d <- town:district
s <- town:district:factor(street)

h <- town:district:factor(street):house

model <- lmer(response~gender+(1|town)+(1|d)+(1|s)+(1|h))

summary(model)

results <- read.table("c:\\temp\\fertilizer.txt",header=T)

attach(results)

we have a single fixed effect (a two-level categorical variable, with fertilizer added

or not) and

 six replicate plants in each treatment, with each plant measured on five occasions (after 2, 4, 6,

8 or 10 weeks of growth). 

The response variable is root length, measured non-destructively through a glass

panel, which is opened to the light only when the root length measurements are being taken. The fixed-effect
formula looks like this:

names(results)

library(nlme)

library(lattice)

results <- groupedData(root~week|plant,outer = ~ fertilizer,results)

plot(results)

plot(results,outer=T)

model <- lme(root~fertilizer,random=~week|plant)

summary(model)

model2 <- aov(root~fertilizer,subset=(week==10))

summary(model2)

summary.lm(model2)

19.7 Time series analysis in mixed-effects models

model <- lme(follicles~sin(2*pi*Time)+cos(2*pi*Time),
data=Ovary,random=~ 1| Mare)

summary(model)

plot(ACF(model),alpha=0.05)

model2 <- update(model,correlation=corARMA(q=2))

anova(model,model2)

model3 <- update(model2,correlation=corAR1())

anova(model2,model3)

plot(model3,resid(.,type="p")~fitted(.)|Mare)

qqnorm(model3,~resid(.)|Mare)