## 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")
## 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)
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="")
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))
## 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)
## 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)
## 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))
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")
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}*")"))
## 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))
## 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}*")"))
## Footnote Code
## Plot of residuals
xyplot(resid(hp.lmer) ~ wattsPerKg, groups=id, type="b", data=humanpower1)
## 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))
## 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)
## The following makes clear M13's difference from other points
qqnorm(ab$b)
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)