## Chapter 10: Multi-level Models, and Repeated Measures
##                  Corn yield measurements example
library(lattice); library(DAAG)
Site <- with(ant111b, reorder(site, harvwt, FUN=mean))
stripplot(Site ~ harvwt, data=ant111b, scales=list(tck=0.5),
          xlab="Harvest weight of corn")

plot of chunk unnamed-chunk-1

## Sec 10.1: Corn Yield Data --- Analysis Using {aov()}
library(DAAG)
ant111b.aov <- aov(harvwt ~ 1 + Error(site), data=ant111b)

summary(ant111b.aov)
## 
## Error: site
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  7   70.3    10.1               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24   13.9   0.578
##                   Interpreting the mean squares
##                    Details of the calculations
##         Practical use of the analysis of variance results
##                  Random effects vs. fixed effects
##            Nested factors -- a variety of applications
## ss 10.1.1: A More Formal Approach
##       Relations between variance components and mean squares
##               Interpretation of variance components
##                      Intra-class correlation

## Sec 10.2: Analysis using {lmer()}, from the { lme4} package
library(lme4)
## Loading required package: Matrix
## Loading required package: Rcpp
ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b)

## Note that there is no degrees of freedom information.
print(ant111b.lmer, ranef.comp="Variance", digits=3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: harvwt ~ 1 + (1 | site)
##    Data: ant111b
## REML criterion at convergence: 94.42
## Random effects:
##  Groups   Name        Variance
##  site     (Intercept) 2.368   
##  Residual             0.578   
## Number of obs: 32, groups:  site, 8
## Fixed Effects:
## (Intercept)  
##        4.29
##              The processing of  output from {lmer()}
print(coef(summary(ant111b.lmer)), digits=3)
##             Estimate Std. Error t value
## (Intercept)     4.29       0.56    7.66
##              Fitted values and residuals in {lmer()}
s2W <- 0.578; s2L <- 2.37; n <- 4
sitemeans <- with(ant111b, sapply(split(harvwt, site), mean))
grandmean <- mean(sitemeans)
shrinkage <- (n*s2L)/(n*s2L+s2W)
grandmean + shrinkage*(sitemeans - grandmean)
##  DBAN  LFAN  NSAN  ORAN  OVAN  TEAN  WEAN  WLAN 
## 4.851 4.212 2.217 6.764 4.801 3.108 5.455 2.925
##
## More directly, use fitted() with the lmer object
unique(fitted(ant111b.lmer))
## [1] 4.851 4.212 2.217 6.764 4.801 3.108 5.455 2.925
##
## Compare with site means
sitemeans
##  DBAN  LFAN  NSAN  ORAN  OVAN  TEAN  WEAN  WLAN 
## 4.885 4.207 2.090 6.915 4.833 3.036 5.526 2.841
##  *Uncertainty in the parameter estimates --- profile  likelihood and alternatives
prof.lmer <- profile(ant111b.lmer)
CI95 <- confint(prof.lmer, level=0.95)
rbind("sigmaL^2"=CI95[1,]^2, "sigma^2"=CI95[2,]^2)
##           2.5 % 97.5 %
## sigmaL^2 0.7965  6.936
## sigma^2  0.3444  1.079
CI95[3,]
##  2.5 % 97.5 % 
##  3.128  5.456
library(lattice)
print(xyplot(prof.lmer, conf=c(50, 80, 95, 99)/100, 
             aspect=0.8, between=list(x=0.35)))

plot of chunk unnamed-chunk-1

##         Handling more than two levels of random variation

## Sec 10.3: Survey Data, with Clustering
## Footnote Code
## Means of like (data frame science: DAAG), by class
classmeans <- with(science,
                   aggregate(like, by=list(PrivPub, Class), mean) )
  # NB: Class identifies classes independently of schools
  #     class identifies classes within schools
names(classmeans) <- c("PrivPub", "Class", "avlike")
with(classmeans, {
     ## Boxplots: class means by Private or Public school
     boxplot(split(avlike, PrivPub), horizontal=TRUE, las=2,
                   xlab = "Class average of score", boxwex = 0.4)
     rug(avlike[PrivPub == "private"], side = 1)
     rug(avlike[PrivPub == "public"], side = 3)
})

plot of chunk unnamed-chunk-1

## ss 10.3.1: Alternative models
science.lmer <- lmer(like ~ sex + PrivPub + (1 | school) +
                     (1 | school:class), data = science,
                     na.action=na.exclude)

print(VarCorr(science.lmer), comp="Variance", digits=3)
##  Groups       Name        Variance
##  school:class (Intercept) 0.321   
##  school       (Intercept) 0.000   
##  Residual                 3.052
print(coef(summary(science.lmer)), digits=2)
##               Estimate Std. Error t value
## (Intercept)       4.72      0.162    29.1
## sexm              0.18      0.098     1.9
## PrivPubpublic     0.41      0.186     2.2
summary(science.lmer)$ngrps
## school:class       school 
##           66           41
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
                      data = science, na.action=na.exclude)
print(VarCorr(science1.lmer), comp="Variance", digits=3)
##  Groups       Name        Variance
##  school:class (Intercept) 0.321   
##  Residual                 3.052
print(coef(summary(science1.lmer)), digits=2)
##               Estimate Std. Error t value
## (Intercept)       4.72      0.162    29.1
## sexm              0.18      0.098     1.9
## PrivPubpublic     0.41      0.186     2.2
library(afex)
## Loading required package: car
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:DAAG':
## 
##     vif
## 
## Loading required package: pbkrtest
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:DAAG':
## 
##     hills
## 
## Loading required package: parallel
## Loading required package: reshape2
## ************
## Welcome to afex. Important notes:
## 
## Due to popular demand, afex doesn't change the contrasts globally anymore.
## To set contrasts globally to contr.sum run set_sum_contrasts().
## To set contrasts globally to the default (treatment) contrasts run set_default_contrasts().
## 
## All afex functions are unaffected by global contrasts and use contr.sum as long as check.contr = TRUE (which is the default).
## ************
mixed(like ~ sex + PrivPub + (1 | school:class),
      data = na.omit(science), method="KR")
## Contrasts set to contr.sum for the following variables: sex, PrivPub, school, class
## Fitting 3 (g)lmer() models:
## [...]
## Obtaining 2 p-values:
## [..]
##    Effect    F ndf     ddf F.scaling p.value
## 1     sex 3.44   1 1379.49      1.00     .06
## 2 PrivPub 4.91   1   60.44      1.00     .03
##              More detailed examination of the output
## Use profile likelihood
pp <- profile(science1.lmer, which="theta_")
  # which="theta_": all random parameters
  # which="beta_": fixed effect parameters
var95 <- confint(pp, level=0.95)^2
  # Square to get variances in place of SDs
rownames(var95) <- c("sigma_Class^2", "sigma^2")
signif(var95, 3)
##               2.5 % 97.5 %
## sigma_Class^2 0.178  0.511
## sigma^2       2.830  3.300
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
                      data = science, na.action=na.omit)
ranf <- ranef(obj = science1.lmer, drop=TRUE)[["school:class"]]
flist <- science1.lmer@flist[["school:class"]]
privpub <- science[match(names(ranf), flist), "PrivPub"]
num <- unclass(table(flist)); numlabs <- pretty(num)
opar <- par(mfrow=c(2,2), pty="s", mgp=c(2.25,0.5,0), mar=c(3.6,3.6,2.1, 0.6))
## Panel A: Plot effect estimates vs numbers
plot(sqrt(num), ranf, xaxt="n", pch=c(1,3)[as.numeric(privpub)],
     xlab="# in class (square root scale)",
     ylab="Estimate of class effect")
lines(lowess(sqrt(num[privpub=="private"]),
             ranf[privpub=="private"], f=1.1), lty=2)
lines(lowess(sqrt(num[privpub=="public"]),
             ranf[privpub=="public"], f=1.1), lty=3)
axis(1, at=sqrt(numlabs), labels=paste(numlabs))
res <- residuals(science1.lmer)
vars <- tapply(res, INDEX=list(flist), FUN=var)*(num-1)/(num-2)
## Panel B: Within class variance estimates vs numbers
plot(sqrt(num), vars, pch=c(1,3)[unclass(privpub)])
lines(lowess(sqrt(num[privpub=="private"]),
             as.vector(vars)[privpub=="private"], f=1.1), lty=2)
lines(lowess(sqrt(num[privpub=="public"]),
             as.vector(vars)[privpub=="public"], f=1.1), lty=3)
## Panel C: Normal probability plot of site effects
qqnorm(ranf, ylab="Ordered site effects", main="")
## Panel D: Normal probability plot of residuals
qqnorm(res, ylab="Ordered w/i class residuals", main="")

plot of chunk unnamed-chunk-1

par(opar)

## ss 10.3.2: Instructive, though faulty, analyses
##                Ignoring class as the random effect
science2.lmer <- lmer(like ~ sex + PrivPub + (1 | school),
                      data = science, na.action=na.exclude)
science2.lmer
## Linear mixed model fit by REML ['lmerMod']
## Formula: like ~ sex + PrivPub + (1 | school)
##    Data: science
## REML criterion at convergence: 5584
## Random effects:
##  Groups   Name        Std.Dev.
##  school   (Intercept) 0.407   
##  Residual             1.794   
## Number of obs: 1383, groups:  school, 41
## Fixed Effects:
##   (Intercept)           sexm  PrivPubpublic  
##         4.738          0.197          0.417
## Footnote Code
##             Ignoring the random structure in the data
## Faulty analysis, using lm
science.lm <- lm(like ~ sex + PrivPub, data=science)
summary(science.lm)$coef
##               Estimate Std. Error t value   Pr(>|t|)
## (Intercept)     4.7402    0.09955  47.616 1.545e-293
## sexm            0.1509    0.09860   1.531  1.261e-01
## PrivPubpublic   0.3951    0.10511   3.759  1.779e-04
## ss 10.3.3: Predictive accuracy

## Sec 10.4: A Multi-level Experimental Design
## ss 10.4.1: The anova table
## Analysis of variance: data frame kiwishade (DAAG)
kiwishade.aov <- aov(yield ~ shade + Error(block/shade),
                     data=kiwishade)

summary(kiwishade.aov)
## 
## Error: block
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  2    172    86.2               
## 
## Error: block:shade
##           Df Sum Sq Mean Sq F value Pr(>F)
## shade      3   1395     465    22.2 0.0012
## Residuals  6    126      21               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 36    439    12.2
## ss 10.4.2: Expected values of mean squares
model.tables(kiwishade.aov, type="means")
## Tables of means
## Grand mean
##       
## 96.53 
## 
##  shade 
## shade
##    none Aug2Dec Dec2Feb Feb2May 
##  100.20  103.23   89.92   92.77
## Footnote Code
## Calculate treatment means
with(kiwishade, sapply(split(yield, shade), mean))
##    none Aug2Dec Dec2Feb Feb2May 
##  100.20  103.23   89.92   92.77
## ss 10.4.3: * The analysis of variance sums of squares  breakdown
## Footnote Code
## For each plot, calculate mean, and differences from the mean
vine <- paste("vine", rep(1:4, 12), sep="")
vine1rows <- seq(from=1, to=45, by=4)
kiwivines <- unstack(kiwishade, yield ~ vine)
kiwimeans <- apply(kiwivines, 1, mean)
kiwivines <- cbind(kiwishade[vine1rows,  c("block","shade")],
                   Mean=kiwimeans, kiwivines-kiwimeans)
kiwivines <- with(kiwivines, kiwivines[order(block, shade), ])
mean(kiwimeans)      # the grand mean
## [1] 96.53
## ss 10.4.4: The variance components
## ss 10.4.5: The mixed model analysis
kiwishade.lmer <- lmer(yield ~ shade + (1|block) + (1|block:plot),
                         data=kiwishade)
  # block:shade is an alternative to block:plot

print(kiwishade.lmer, ranef.comp="Variance", digits=3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: yield ~ shade + (1 | block) + (1 | block:plot)
##    Data: kiwishade
## REML criterion at convergence: 252
## Random effects:
##  Groups     Name        Variance
##  block:plot (Intercept)  2.19   
##  block      (Intercept)  4.08   
##  Residual               12.18   
## Number of obs: 48, groups:  block:plot, 12; block, 3
## Fixed Effects:
##  (Intercept)  shadeAug2Dec  shadeDec2Feb  shadeFeb2May  
##       100.20          3.03        -10.28         -7.43
##                  Residuals and estimated effects
## Footnote Code
## Simplified version of plot
xyplot(residuals(kiwishade.lmer) ~ fitted(kiwishade.lmer)|block, data=kiwishade,
                groups=shade, layout=c(3,1), par.strip.text=list(cex=1.0),
                xlab="Fitted values (Treatment + block + plot effects)",
                ylab="Residuals", pch=1:4, grid=TRUE, aspect=1,
                scales=list(x=list(alternating=FALSE), tck=0.5),
                key=list(space="top", points=list(pch=1:4),
                         text=list(labels=levels(kiwishade$shade)),columns=4))

plot of chunk unnamed-chunk-1

## Footnote Code
## Simplified version of graph that shows the plot effects
ploteff <- ranef(kiwishade.lmer, drop=TRUE)[[1]]
qqmath(ploteff, xlab="Normal quantiles", ylab="Plot effect estimates",
       aspect=1, scales=list(tck=0.5))