## 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)) 

plot of chunk unnamed-chunk-1

par(mfrow=c(1,1))

## ss 10.4.2: Estimates of model parameters 
it2.reml <- update(it2.lmer, method="REML") 
## Warning: Argument 'method' is deprecated. Use the REML argument to specify
## ML or REML estimation.
summary(it2.reml) 
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(it) ~ (tint + target + agegp + sex)^2 + (1 | id)
##    Data: tinting
## 
##      AIC      BIC   logLik deviance df.resid 
##     -3.7     50.7     18.9    -37.7      165 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.456 -0.640 -0.107  0.545  2.474 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.1204   0.347   
##  Residual             0.0293   0.171   
## Number of obs: 182, groups:  id, 26
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)             3.61907    0.11987   30.19
## tint.L                  0.16095    0.04265    3.77
## tint.Q                  0.02096    0.04359    0.48
## targethicon            -0.11807    0.04081   -2.89
## agegpOlder              0.47121    0.21449    2.20
## sexm                    0.08213    0.21449    0.38
## tint.L:targethicon     -0.09193    0.04441   -2.07
## tint.Q:targethicon     -0.00722    0.04648   -0.16
## tint.L:agegpOlder       0.13075    0.04742    2.76
## tint.Q:agegpOlder       0.06972    0.05013    1.39
## tint.L:sexm            -0.09794    0.04742   -2.07
## tint.Q:sexm             0.00542    0.05013    0.11
## targethicon:agegpOlder -0.13887    0.05634   -2.46
## targethicon:sexm        0.07785    0.05634    1.38
## agegpOlder:sexm         0.33164    0.29999    1.11
## 
## Correlation of Fixed Effects:
##             (Intr) tint.L tint.Q trgthc aggpOl sexm   tnt.L:t tnt.Q:t
## tint.L       0.000                                                   
## tint.Q       0.000  0.038                                            
## targethicon -0.188  0.066 -0.037                                     
## agegpOlder  -0.551  0.000  0.000  0.062                              
## sexm        -0.551  0.000  0.000  0.062  0.293                       
## tnt.L:trgth  0.000 -0.595  0.000  0.073  0.000  0.000                
## tnt.Q:trgth  0.000  0.000 -0.556 -0.040  0.000  0.000  0.079         
## tnt.L:ggpOl  0.000 -0.342 -0.034 -0.059  0.000  0.000  0.000   0.000 
## tnt.Q:ggpOl  0.000 -0.033 -0.354  0.032  0.000  0.000  0.000   0.000 
## tint.L:sexm  0.000 -0.342 -0.034 -0.059  0.000  0.000  0.000   0.000 
## tint.Q:sexm  0.000 -0.033 -0.354  0.032  0.000  0.000  0.000   0.000 
## trgthcn:ggO  0.080 -0.048  0.027 -0.425 -0.146  0.056  0.000   0.000 
## trgthcn:sxm  0.080 -0.048  0.027 -0.425  0.056 -0.146  0.000   0.000 
## aggpOldr:sx  0.385  0.000  0.000  0.000 -0.699 -0.699  0.000   0.000 
##             tn.L:O tn.Q:O tnt.L:s tnt.Q:s trgt:O trgth:
## tint.L                                                 
## tint.Q                                                 
## targethicon                                            
## agegpOlder                                             
## sexm                                                   
## tnt.L:trgth                                            
## tnt.Q:trgth                                            
## tnt.L:ggpOl                                            
## tnt.Q:ggpOl  0.096                                     
## tint.L:sexm -0.385 -0.037                              
## tint.Q:sexm -0.037 -0.385  0.096                       
## trgthcn:ggO  0.140 -0.076 -0.054   0.029               
## trgthcn:sxm -0.054  0.029  0.140  -0.076  -0.385       
## aggpOldr:sx  0.000  0.000  0.000   0.000   0.000  0.000
 # NB: The final column, giving degrees of freedom, is not in the 
 # summary output for version 0.995-2 of lme4. It is our addition. 

uid <- unique(tinting$id) 
subs <- match(uid, tinting$id) 
with(tinting, table(sex[subs], agegp[subs])) 
##    
##     Younger Older
##   f       9     4
##   m       4     9
## Sec 10.5: A Generalized Linear Mixed Model
## Inclusion of 'Bank' (no A moths) leads to convergence problems
## Treat comparisons with 'Bank' separately (no A moths; 1 transect)
moths40 <- subset(moths, habitat!="Bank")
moths40$transect <- 1:40  # Each row is from a different transect 
moths40$habitat <- relevel(moths40$habitat, ref="Lowerside") 
(A.glmer <-  glmer(A~habitat+log(meters)+(1|transect),  
                   family=poisson, data=moths40))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: A ~ habitat + log(meters) + (1 | transect)
##    Data: moths40
##      AIC      BIC   logLik deviance df.resid 
##   207.73   222.93   -94.86   189.73       31 
## Random effects:
##  Groups   Name        Std.Dev.
##  transect (Intercept) 0.483   
## Number of obs: 40, groups:  transect, 40
## Fixed Effects:
##      (Intercept)  habitatDisturbed     habitatNEsoak     habitatNWsoak  
##           1.0201           -1.2625           -0.8431            1.5517  
##    habitatSEsoak     habitatSWsoak  habitatUpperside       log(meters)  
##           0.0532            0.2506           -0.1707            0.1544
A.glm <- glm(A~habitat+log(meters), data=moths40, family=quasipoisson)
par(mfrow=c(2,2), pty="s") 
fitglm <- fitted(A.glm) 
fitglmer <- fitted(A.glmer) 
## Plots from quasipoisson analysis 
plot(resid(A.glm) ~ moths40$meters, log="x") 
plot(resid(A.glm) ~ fitglm, log="x") 
## 
## Plots from glmer mixed model analysis 
plot(resid(A.glmer) ~ moths40$meters, log="x") 
plot(resid(A.glmer) ~ fitglmer, log="x") 

plot of chunk unnamed-chunk-1

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

##         Mixed models with a binomial error and logit link 

## Sec 10.6: Repeated Measures in Time 
##             The theory of repeated measures modeling  
##                       *Correlation structure 
##         Different approaches to repeated measures analysis 
## ss 10.6.1: Example -- random variation between profiles 
## Footnote Code
## Plot points and fitted lines (panel A) 
library(lattice) 
xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, 
        panel=function(x,y,subscripts,groups,...){ 
           yhat <- fitted(lm(y ~ groups*x)) 
           panel.superpose(x, y, subscripts, groups, pch=1:5 
) 
           panel.superpose(x, yhat, subscripts, groups, type="l") 
        }, 
        xlab="Watts per kilogram", 
        ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")")) 

plot of chunk unnamed-chunk-1

##               Separate lines for different athletes 
## Calculate intercepts and slopes; plot Slopes vs Intercepts 
## Uses the function lmList() from the lme4 package 
library(lme4) 
hp.lmList <- lmList(o2 ~ wattsPerKg | id, data=humanpower1) 
coefs <- coef(hp.lmList) 
names(coefs) <- c("Intercept", "Slope") 
plot(Slope ~ Intercept, data=coefs) 
abline(lm(Slope~Intercept, data=coefs)) 

plot of chunk unnamed-chunk-1

##                    A random coefficients model 
hp.lmer <- lmer(o2 ~ wattsPerKg + (wattsPerKg | id), 
                data=humanpower1) 
summary(hp.lmer) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: o2 ~ wattsPerKg + (wattsPerKg | id)
##    Data: humanpower1
## 
## REML criterion at convergence: 124.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9947 -0.4878 -0.0835  0.4616  2.1212 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 50.73    7.12          
##           wattsPerKg   7.15    2.67     -1.00
##  Residual              4.13    2.03          
## Number of obs: 28, groups:  id, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)     2.09       3.78    0.55
## wattsPerKg     13.91       1.36   10.23
## 
## Correlation of Fixed Effects:
##            (Intr)
## wattsPerKg -0.992
hat<-fitted(hp.lmer) 
lmhat<-with(humanpower1, fitted(lm(o2 ~ id*wattsPerKg))) 
panelfun <- 
  function(x, y, subscripts, groups, ...){ 
           panel.superpose(x, hat, subscripts, groups, type="l",lty=2) 
           panel.superpose(x, lmhat, subscripts, groups, type="l",lty=1) 
} 
xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, panel=panelfun,  
       xlab="Watts", 
       ylab=expression("Oxygen intake ("*ml.min^{-1}*"."*kg^{-1}*")")) 

plot of chunk unnamed-chunk-1

## Footnote Code
## Plot of residuals 
xyplot(resid(hp.lmer) ~ wattsPerKg, groups=id, type="b", data=humanpower1) 

plot of chunk unnamed-chunk-1

## Footnote Code
## Derive the sd from the data frame coefs that was calculated above 
sd(coefs$Slope) 
## [1] 3.278
## ss 10.6.2: Orthodontic measurements on chlldren 
##                    Preliminary data exploration 
## Footnote Code
## Plot showing pattern of change for each of the 25 individuals 
library(MEMSS, quietly=TRUE) 
## 
## Attaching package: 'MEMSS'
## 
## The following objects are masked from 'package:datasets':
## 
##     CO2, Orange, Theoph
xyplot(distance ~ age | Subject, groups=Sex, data=Orthodont, 
       scale=list(y=list(log=2)), type=c("p","r"), layout=c(11,3)) 

plot of chunk unnamed-chunk-1

## Use lmList() to find the slopes 
ab <- coef(lmList(distance ~ age | Subject, Orthodont)) 
names(ab) <- c("a", "b") 
## Obtain the intercept at x=mean(x)  
## (For each subject, this is independent of the slope) 
ab$ybar <- ab$a + ab$b*11  # mean age is 11, for each subject. 
sex <- substring(rownames(ab), 1 ,1) 
plot(ab[, 3], ab[, 2], col=c(F="gray40", M="black")[sex],              
     pch=c(F=1, M=3)[sex], xlab="Intercept", ylab="Slope") 
extremes <- ab$ybar %in% range(ab$ybar) |  
            ab$b %in% range(ab$b[sex=="M"]) | 
            ab$b %in% range(ab$b[sex=="F"]) 
text(ab[extremes, 3], ab[extremes, 2], rownames(ab)[extremes], pos=4, xpd=TRUE) 

plot of chunk unnamed-chunk-1

## The following makes clear M13's difference from other points 
qqnorm(ab$b) 

plot of chunk unnamed-chunk-1

Orthodont$logdist <- log(Orthodont$distance) 
## Now repeat, with logdist replacing distance 

## Footnote Code
## Compare males slopes with female slopes 
Orthodont$logdist <- log(Orthodont$distance) 
ablog <- coef(lmList(logdist ~ age | Subject, Orthodont)) 
names(ablog) <- c("a", "b") 
## Obtain the intercept at mean age (= 11), for each subject 
## (For each subject, this is independent of the slope) 
ablog$ybar <- with(ablog, a + b*11) 
extreme.males <- rownames(ablog) %in% c("M04","M13") 
sex <- substring(rownames(ab), 1, 1) 
with(ablog,  
  t.test(b[sex=="F"], b[sex=="M" & !extreme.males], var.equal=TRUE)) 
## 
##  Two Sample t-test
## 
## data:  b[sex == "F"] and b[sex == "M" & !extreme.males]
## t = -2.32, df = 23, p-value = 0.02957
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0160529 -0.0009191
## sample estimates:
## mean of x mean of y 
##   0.02115   0.02963
  # Specify var.equal=TRUE, to allow comparison with anova output 

##                    A random coefficients model 
keep <- !(Orthodont$Subject%in%c("M04","M13")) 
orthdiff.lmer <- lmer(logdist ~ Sex * I(age-11) + (I(age-11) | Subject),  
                      data=Orthodont, subset=keep, method="ML")  
## Warning: Argument 'method' is deprecated. Use the REML argument to specify
## ML or REML estimation.
orthdiff.lmer 
## Linear mixed model fit by REML ['lmerMod']
## Formula: logdist ~ Sex * I(age - 11) + (I(age - 11) | Subject)
##    Data: Orthodont
##  Subset: keep
## REML criterion at convergence: -232.2
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 0.079581     
##           I(age - 11) 0.000918 1.00
##  Residual             0.048764     
## Number of obs: 100, groups:  Subject, 25
## Fixed Effects:
##         (Intercept)              SexMale          I(age - 11)  
##             3.11451              0.09443              0.02115  
## SexMale:I(age - 11)  
##             0.00849
orthsame.lmer <-  
  lmer(logdist ~ Sex + I(age - 11) + (I(age - 11) | Subject), 
       data=Orthodont, method="ML", subset=keep) 
## Warning: Argument 'method' is deprecated. Use the REML argument to specify
## ML or REML estimation.
anova(orthsame.lmer, orthdiff.lmer) 
## refitting model(s) with ML (instead of REML)
## Data: Orthodont
## Subset: keep
## Models:
## orthsame.lmer: logdist ~ Sex + I(age - 11) + (I(age - 11) | Subject)
## orthdiff.lmer: logdist ~ Sex * I(age - 11) + (I(age - 11) | Subject)
##               Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## orthsame.lmer  7 -245 -227    130     -259                        
## orthdiff.lmer  8 -247 -226    132     -263  3.71      1      0.054
orthdiffr.lmer <- update(orthdiff.lmer, method="REML") 
## Warning: Argument 'method' is deprecated. Use the REML argument to specify
## ML or REML estimation.
summary(orthdiffr.lmer) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: logdist ~ Sex * I(age - 11) + (I(age - 11) | Subject)
##    Data: Orthodont
##  Subset: keep
## 
## REML criterion at convergence: -232.2
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.325 -0.477  0.001  0.522  3.939 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 6.33e-03 0.079581     
##           I(age - 11) 8.42e-07 0.000918 1.00
##  Residual             2.38e-03 0.048764     
## Number of obs: 100, groups:  Subject, 25
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          3.11451    0.02510   124.1
## SexMale              0.09443    0.03354     2.8
## I(age - 11)          0.02115    0.00330     6.4
## SexMale:I(age - 11)  0.00849    0.00441     1.9
## 
## Correlation of Fixed Effects:
##             (Intr) SexMal I(-11)
## SexMale     -0.748              
## I(age - 11)  0.080 -0.060       
## SxMl:I(-11) -0.060  0.080 -0.748
## At the time of writing, this is not possible 
## orth.mcmc <- mcmcsamp(orthdiffr.lmer, n=1000) 
## HPDinterval(orth.mcmc) 

## Use confint() instead
confint(orthdiffr.lmer, method="boot", quiet=TRUE) 
## Warning: Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
##                                          2.5 %  97.5 %
## sd_(Intercept)|Subject               5.573e-02 0.10489
## cor_I(age - 11).(Intercept)|Subject -1.000e+00 1.00000
## sd_I(age - 11)|Subject               9.968e-05 0.00939
## sigma                                3.982e-02 0.05629
## (Intercept)                          3.066e+00 3.16544
## SexMale                              2.334e-02 0.16773
## I(age - 11)                          1.470e-02 0.02802
## SexMale:I(age - 11)                 -8.840e-04 0.01737
  # A few of the bootstrap fits may fail. Providing there
  # are not too many failures, it seems safe to ignore these
  # Run the above with 'quiet=FALSE' a few times to check.

##                Correlation between successive times 
res <- resid(orthdiffr.lmer) 
Subject <- factor(Orthodont$Subject[keep]) 
orth.acf <- t(sapply(split(res, Subject),  
                     function(x)acf(x, lag=4, plot=FALSE)$acf)) 
## Calculate respective proportions of Subjects for which  
## autocorrelations at lags 1, 2 and 3 are greater than zero. 
apply(orth.acf[,-1], 2, function(x)sum(x>0)/length(x)) 
## [1] 0.20 0.32 0.36
##             *The variance for the difference in slopes 

## Sec 10.7: Further Notes on Multi-level and Other Models with CorrelatedErrors 
## ss 10.7.1: Different sources of variance -- complication or focus of interest? 
## ss 10.7.2: Predictions from models with a complex error structure 
##  Consequences from assuming an overly simplistic error structure 
## ss 10.7.3: An historical perspective on multi-level models 
## ss 10.7.4: Meta-analysis 
## ss 10.7.5: Functional data analysis 
## ss 10.7.6:  Error structure in explanatory variables 

## Sec 10.8: Recap 

## Sec 10.9: Further Reading 
##                             References 
##                   References for further reading 
##           Analysis of variance with multiple error terms 
##              Multi-level models and repeated measures 
##                           Meta-analysis 
##                             References 
##                   References for further reading 
n.omit <- 2 
take <- rep(TRUE, 48) 
take[sample(1:48,2)] <- FALSE 
kiwishade.lmer <- lmer(yield ~ shade + (1|block) + (1|block:plot), 
                      data = kiwishade,subset=take) 
vcov <- VarCorr(kiwishade.lmer) 
print(vcov) 
##  Groups     Name        Std.Dev.
##  block:plot (Intercept) 1.40    
##  block      (Intercept) 1.93    
##  Residual               3.51
## Or, print(vcov^2)

cult.lmer <- lmer(ct ~ Cultivar + Dose + factor(year) +  
                  (-1 + Dose | gp), data = sorption,  
                  REML=TRUE) 
cultdose.lmer <- lmer(ct ~ Cultivar/Dose + factor(year) + 
                      (-1 + Dose | gp), data = sorption,  
                      REML=TRUE)