## Chapter 10: Multi-level Models, and Repeated Measures 

## Sec 10.1: A One-Way Random Effects Model 
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

## ss 10.1.1: Analysis with {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.2: A More Formal Approach 
##       Relations between variance components and mean squares 
##               Interpretation of variance components 
##                      Intra-class correlation 
## ss 10.1.3: Analysis using {lmer()} 
library(lme4) 
## Loading required package: Matrix
## Loading required package: Rcpp
ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b) 

## Output is from version 0.999375-28 of lmer 
## Note that there is no degrees of freedom information. 
ant111b.lmer 
## Linear mixed model fit by REML ['lmerMod']
## Formula: harvwt ~ 1 + (1 | site)
##    Data: ant111b
## REML criterion at convergence: 94.42
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 1.54    
##  Residual             0.76    
## Number of obs: 32, groups:  site, 8
## Fixed Effects:
## (Intercept)  
##        4.29
##              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 
CI95 <- confint(profile(ant111b.lmer), level=0.95) 
rbind("sigmaL^2"=CI95[1,]^2, "sigma^2"=exp(CI95[2,])^2, 
      "(Intercept)"=CI95[3,]) 
##              2.5 % 97.5 %
## sigmaL^2    0.7965  6.936
## sigma^2     3.2337  7.981
## (Intercept) 3.1277  5.456
##         Handling more than two levels of random variation 
## house.lmer <- lmer(price ~ 1 + (1|city) + (1|city:suburb)) 


## Sec 10.2: 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") 
attach(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

detach(classmeans) 

## ss 10.2.1: Alternative models 
science.lmer <- lmer(like ~ sex + PrivPub + (1 | school) + 
                     (1 | school:class), data = science,  
                     na.action=na.exclude) 
print(summary(science.lmer), show.resid=FALSE) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: like ~ sex + PrivPub + (1 | school) + (1 | school:class)
##    Data: science
## 
## REML criterion at convergence: 5546
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  school:class (Intercept) 0.321    0.566   
##  school       (Intercept) 0.000    0.000   
##  Residual                 3.052    1.747   
## Number of obs: 1383, groups:  school:class, 66; school, 41
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)     4.7218     0.1624   29.07
## sexm            0.1823     0.0982    1.86
## PrivPubpublic   0.4117     0.1857    2.22
## 
## Correlation of Fixed Effects:
##             (Intr) sexm  
## sexm        -0.309       
## PrivPubpblc -0.795  0.012
## Footnote Code
## The variances are included in the output from VarCorr() 
VarCorr(science.lmer)  # Displayed output differs slightly 
##  Groups       Name        Std.Dev.
##  school:class (Intercept) 0.566   
##  school       (Intercept) 0.000   
##  Residual                 1.747
 # The between students (Residual) component of variance is 
 # attr(VarCorr(science.lmer),"sc")^2 

science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class), 
                      data = science, na.action=na.exclude) 
print(summary(science1.lmer), show.resid=FALSE) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: like ~ sex + PrivPub + (1 | school:class)
##    Data: science
## 
## REML criterion at convergence: 5546
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  school:class (Intercept) 0.321    0.566   
##  Residual                 3.052    1.747   
## Number of obs: 1383, groups:  school:class, 66
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)     4.7218     0.1624   29.07
## sexm            0.1823     0.0982    1.86
## PrivPubpublic   0.4117     0.1857    2.22
## 
## Correlation of Fixed Effects:
##             (Intr) sexm  
## sexm        -0.309       
## PrivPubpblc -0.795  0.012
# mcmcsamp() did noi work reliably, & is no longer available
# set.seed(41) 
# science1.mcmc <- mcmcsamp(science1.lmer, n=1000) 
# samps <- VarCorr(science1.mcmc, type="varcov") 
# colnames(samps) <- c("sigma_Class^2", "sigma^2") 
# signif(HPDinterval(samps, prob=0.95), 3) 

## Use profile likelihood instead
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) 
par(mfrow=c(2,2), pty="s") 
## 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) 
## Within plot 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) 
## Normal probability plot of site effects     
qqnorm(ranf, ylab="Ordered site effects", main="") 
## Normal probability plot of residuals 
qqnorm(res, ylab="Ordered w/i class residuals", main="") 

plot of chunk unnamed-chunk-1

par(mfrow=c(1,1), pty="m") 

## ss 10.2.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
## The numerical values can be extracted from 
VarCorr(science2.lmer)  # The within schools (Residual) component of variance 
##  Groups   Name        Std.Dev.
##  school   (Intercept) 0.407   
##  Residual             1.794
                        # is the square of the scale parameter 
 

##             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.2.3: Predictive accuracy 

## Sec 10.3: A Multi-level Experimental Design 
## ss 10.3.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.3.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.3.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.3.4: The variance components 
## ss 10.3.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 

kiwishade.lmer 
## 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        Std.Dev.
##  block:plot (Intercept) 1.48    
##  block      (Intercept) 2.02    
##  Residual               3.49    
## 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",  
       scales=list(tck=0.5), aspect=1) 

plot of chunk unnamed-chunk-1

## Footnote Code
## Overlaid normal probability plots of 2 sets of simulated effects 
## To do more simulations, change nsim as required, and re-execute 
simvals <- simulate(kiwishade.lmer, nsim=2) 
simeff <- apply(simvals, 2, function(y) scale(ranef(refit(kiwishade.lmer, y), 
                                                drop=TRUE)[[1]])) 
simeff <- data.frame(v1=simeff[,1], v2=simeff[,2]) 
qqmath(~ v1+v2, data=simeff, xlab="Normal quantiles", 
       ylab="Simulated plot effects\n(2 sets, standardized)", 
              scales=list(tck=0.5), aspect=1) 

plot of chunk unnamed-chunk-1

## ss 10.3.6: Predictive accuracy 

## Sec 10.4: Within and Between Subject Effects 
##                       Model fitting criteria 
## ss 10.4.1: Model selection 
## Change initial letters of levels of tinting$agegp to upper case 
library(R.utils)
## Loading required package: R.oo
## Loading required package: R.methodsS3
## R.methodsS3 v1.6.1 (2014-01-04) successfully loaded. See ?R.methodsS3 for help.
## R.oo v1.18.0 (2014-02-22) successfully loaded. See ?R.oo for help.
## 
## Attaching package: 'R.oo'
## 
## The following objects are masked from 'package:methods':
## 
##     getClasses, getMethods
## 
## The following objects are masked from 'package:base':
## 
##     attach, detach, gc, load, save
## 
## R.utils v1.32.4 (2014-05-14) successfully loaded. See ?R.utils for help.
## 
## Attaching package: 'R.utils'
## 
## The following object is masked from 'package:utils':
## 
##     timestamp
## 
## The following objects are masked from 'package:base':
## 
##     cat, commandArgs, getOption, inherits, isOpen, parse, warnings
levels(tinting$agegp) <- capitalize(levels(tinting$agegp)) 
## Fit all interactions: data frame tinting (DAAG) 
it3.lmer <- lmer(log(it) ~ tint*target*agegp*sex + (1 | id), 
                 data=tinting, REML=FALSE) 

it2.lmer <- lmer(log(it) ~ (tint+target+agegp+sex)^2 + (1 | id),  
                 data=tinting, REML=FALSE) 

it1.lmer <- lmer(log(it)~(tint+target+agegp+sex) + (1 | id),  
                 data=tinting, REML=FALSE) 

anova(it1.lmer, it2.lmer, it3.lmer) 
## Data: tinting
## Models:
## it1.lmer: log(it) ~ (tint + target + agegp + sex) + (1 | id)
## it2.lmer: log(it) ~ (tint + target + agegp + sex)^2 + (1 | id)
## it3.lmer: log(it) ~ tint * target * agegp * sex + (1 | id)
##          Df   AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## it1.lmer  8  1.14 26.8   7.43    -14.9                        
## it2.lmer 17 -3.74 50.7  18.87    -37.7 22.88      9     0.0065
## it3.lmer 26  8.15 91.5  21.93    -43.9  6.11      9     0.7288
## Footnote Code
## Code that gives the first four such plots, for the observed data 
par(mfrow=c(2,2))
interaction.plot(tinting$tint, tinting$agegp, log(tinting$it)) 
interaction.plot(tinting$target, tinting$sex, log(tinting$it)) 
interaction.plot(tinting$tint, tinting$target, log(tinting$it)) 
interaction.plot(tinting$tint, tinting$sex, log(tinting$it))