277 lines
8.7 KiB
R
277 lines
8.7 KiB
R
#-*- R -*-
|
|
|
|
## Script from Fourth Edition of `Modern Applied Statistics with S'
|
|
|
|
# Chapter 10 Random and Mixed Effects
|
|
|
|
library(MASS)
|
|
library(lattice)
|
|
pdf(file="ch10.pdf", width=8, height=6, pointsize=9)
|
|
options(width=65, digits=5)
|
|
library(nlme)
|
|
|
|
# 10.1 Linear models
|
|
|
|
xyplot(Y ~ EP | No, data = petrol,
|
|
xlab = "ASTM end point (deg. F)",
|
|
ylab = "Yield as a percent of crude",
|
|
panel = function(x, y) {
|
|
panel.grid()
|
|
m <- sort.list(x)
|
|
panel.xyplot(x[m], y[m], type = "b", cex = 0.5)
|
|
})
|
|
|
|
Petrol <- petrol
|
|
names(Petrol)
|
|
Petrol[, 2:5] <- scale(Petrol[, 2:5], scale = FALSE)
|
|
pet1.lm <- lm(Y ~ No/EP - 1, Petrol)
|
|
matrix(round(coef(pet1.lm), 2), 2, 10, byrow = TRUE,
|
|
dimnames = list(c("b0", "b1"), levels(Petrol$No)))
|
|
|
|
pet2.lm <- lm(Y ~ No - 1 + EP, Petrol)
|
|
anova(pet2.lm, pet1.lm)
|
|
|
|
pet3.lm <- lm(Y ~ SG + VP + V10 + EP, Petrol)
|
|
anova(pet3.lm, pet2.lm)
|
|
|
|
pet3.lme <- lme(Y ~ SG + VP + V10 + EP,
|
|
random = ~ 1 | No, data = Petrol)
|
|
summary(pet3.lme)
|
|
|
|
pet3.lme <- update(pet3.lme, method = "ML")
|
|
summary(pet3.lme)
|
|
|
|
anova(pet3.lme, pet3.lm)
|
|
|
|
pet4.lme <- update(pet3.lme, fixed = Y ~ V10 + EP)
|
|
anova(pet4.lme, pet3.lme)
|
|
fixed.effects(pet4.lme)
|
|
coef(pet4.lme)
|
|
|
|
pet5.lme <- update(pet4.lme, random = ~ 1 + EP | No)
|
|
anova(pet4.lme, pet5.lme)
|
|
|
|
nl1 <- nlschools
|
|
attach(nl1)
|
|
classMeans <- tapply(IQ, class, mean)
|
|
nl1$IQave <- classMeans[as.character(class)]
|
|
detach()
|
|
cen <- c("IQ", "IQave", "SES")
|
|
nl1[cen] <- scale(nl1[cen], center = TRUE, scale = FALSE)
|
|
|
|
options(contrasts = c("contr.treatment", "contr.poly"))
|
|
nl.lme <- lme(lang ~ IQ*COMB + IQave + SES,
|
|
random = ~ IQ | class, data = nl1)
|
|
summary(nl.lme)
|
|
|
|
## singular.ok = TRUE is the default in R
|
|
summary(lm(lang ~ IQ*COMB + SES + class, data = nl1, singular.ok = TRUE))
|
|
|
|
nl2 <- cbind(aggregate(nl1[c(1,7)], list(class = nl1$class), mean),
|
|
unique(nl1[c("class", "COMB", "GS")]))
|
|
summary(lm(lang ~ IQave + COMB, data = nl2, weights = GS))
|
|
|
|
sitka.lme <- lme(size ~ treat*ordered(Time),
|
|
random = ~1 | tree, data = Sitka, method = "ML")
|
|
Sitka <- Sitka # make a local copy for S-PLUS
|
|
attach(Sitka)
|
|
Sitka$treatslope <- Time * (treat == "ozone")
|
|
detach()
|
|
sitka.lme2 <- update(sitka.lme,
|
|
fixed = size ~ ordered(Time) + treat + treatslope)
|
|
anova(sitka.lme, sitka.lme2)
|
|
|
|
# fitted curves
|
|
matrix(fitted(sitka.lme2, level = 0)[c(301:305, 1:5)],
|
|
2, 5, byrow = TRUE,
|
|
dimnames = list(c("control", "ozone"), unique(Sitka$Time)))
|
|
|
|
|
|
# 10.2 Classic nested designs
|
|
|
|
if(FALSE) {
|
|
summary(raov(Conc ~ Lab/Bat, data = coop, subset = Spc=="S1"))
|
|
|
|
is.random(coop) <- T
|
|
is.random(coop$Spc) <- F
|
|
is.random(coop)
|
|
|
|
varcomp(Conc ~ Lab/Bat, data = coop, subset = Spc=="S1")
|
|
|
|
varcomp(Conc ~ Lab/Bat, data = coop, subset = Spc=="S1",
|
|
method = c("winsor", "minque0"))
|
|
}
|
|
|
|
#oats <- oats # make a local copy: needed in S-PLUS
|
|
oats$Nf <- ordered(oats$N, levels = sort(levels(oats$N)))
|
|
|
|
oats.aov <- aov(Y ~ Nf*V + Error(B/V), data = oats, qr = TRUE)
|
|
summary(oats.aov)
|
|
summary(oats.aov, split = list(Nf = list(L = 1, Dev = 2:3)))
|
|
|
|
plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]]))
|
|
abline(h = 0, lty = 2)
|
|
oats.pr <- proj(oats.aov)
|
|
qqnorm(oats.pr[[4]][,"Residuals"], ylab = "Stratum 4 residuals")
|
|
qqline(oats.pr[[4]][,"Residuals"])
|
|
|
|
oats.aov <- aov(Y ~ N + V + Error(B/V), data = oats, qr = TRUE)
|
|
model.tables(oats.aov, type = "means", se = TRUE)
|
|
# we can get the unimplemented standard errors from
|
|
se.contrast(oats.aov, list(N == "0.0cwt", N == "0.2cwt"), data=oats)
|
|
se.contrast(oats.aov, list(V == "Golden.rain", V == "Victory"), data=oats)
|
|
|
|
|
|
# is.random(oats$B) <- T
|
|
# varcomp(Y ~ N + V + B/V, data = oats)
|
|
|
|
lme(Conc ~ 1, random = ~1 | Lab/Bat, data = coop,
|
|
subset = Spc=="S1")
|
|
|
|
options(contrasts = c("contr.treatment", "contr.poly"))
|
|
summary(lme(Y ~ N + V, random = ~1 | B/V, data = oats))
|
|
|
|
|
|
# 10.3 Non-linear mixed effects models
|
|
|
|
options(contrasts = c("contr.treatment", "contr.poly"))
|
|
sitka.nlme <- nlme(size ~ A + B * (1 - exp(-(Time-100)/C)),
|
|
fixed = list(A ~ treat, B ~ treat, C ~ 1),
|
|
random = A + B ~ 1 | tree, data = Sitka,
|
|
start = list(fixed = c(2, 0, 4, 0, 100)), verbose = TRUE)
|
|
|
|
summary(sitka.nlme)
|
|
|
|
summary(update(sitka.nlme,
|
|
corr = corCAR1(0.95, ~Time | tree)))
|
|
|
|
Fpl <- deriv(~ A + (B-A)/(1 + exp((log(d) - ld50)/th)),
|
|
c("A","B","ld50","th"), function(d, A, B, ld50, th) {})
|
|
|
|
st <- coef(nls(BPchange ~ Fpl(Dose, A, B, ld50, th),
|
|
start = c(A = 25, B = 0, ld50 = 4, th = 0.25),
|
|
data = Rabbit))
|
|
Rc.nlme <- nlme(BPchange ~ Fpl(Dose, A, B, ld50, th),
|
|
fixed = list(A ~ 1, B ~ 1, ld50 ~ 1, th ~ 1),
|
|
random = A + ld50 ~ 1 | Animal, data = Rabbit,
|
|
subset = Treatment == "Control",
|
|
start = list(fixed = st))
|
|
## The next fails on some R platforms and some versions of nlme
|
|
## Rm.nlme <- update(Rc.nlme, subset = Treatment=="MDL")
|
|
## so update starting values
|
|
st <- coef(nls(BPchange ~ Fpl(Dose, A, B, ld50, th),
|
|
start = c(A = 25, B = 0, ld50 = 4, th = 0.25),
|
|
data = Rabbit, subset = Treatment == "MDL"))
|
|
Rm.nlme <- update(Rc.nlme, subset = Treatment=="MDL",
|
|
start = list(fixed = st))
|
|
|
|
Rc.nlme
|
|
Rm.nlme
|
|
|
|
c1 <- c(28, 1.6, 4.1, 0.27, 0)
|
|
R.nlme1 <- nlme(BPchange ~ Fpl(Dose, A, B, ld50, th),
|
|
fixed = list(A ~ Treatment, B ~ Treatment,
|
|
ld50 ~ Treatment, th ~ Treatment),
|
|
random = A + ld50 ~ 1 | Animal/Run, data = Rabbit,
|
|
start = list(fixed = c1[c(1, 5, 2, 5, 3, 5, 4, 5)]))
|
|
summary(R.nlme1)
|
|
R.nlme2 <- update(R.nlme1,
|
|
fixed = list(A ~ 1, B ~ 1, ld50 ~ Treatment, th ~ 1),
|
|
start = list(fixed = c1[c(1:3, 5, 4)]))
|
|
anova(R.nlme2, R.nlme1)
|
|
summary(R.nlme2)
|
|
|
|
xyplot(BPchange ~ log(Dose) | Animal * Treatment, Rabbit,
|
|
xlab = "log(Dose) of Phenylbiguanide",
|
|
ylab = "Change in blood pressure (mm Hg)",
|
|
subscripts = TRUE, aspect = "xy", panel =
|
|
function(x, y, subscripts) {
|
|
panel.grid()
|
|
panel.xyplot(x, y)
|
|
sp <- spline(x, fitted(R.nlme2)[subscripts])
|
|
panel.xyplot(sp$x, sp$y, type = "l")
|
|
})
|
|
|
|
|
|
# 10.4 Generalized linear mixed models
|
|
|
|
contrasts(bacteria$trt) <- structure(contr.sdif(3),
|
|
dimnames = list(NULL, c("drug", "encourage")))
|
|
summary(glm(y ~ trt * week, binomial, data = bacteria))
|
|
summary(glm(y ~ trt + week, binomial, data = bacteria))
|
|
|
|
summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data = epil))
|
|
|
|
epil2 <- epil[epil$period == 1, ]
|
|
epil2["period"] <- rep(0, 59); epil2["y"] <- epil2["base"]
|
|
epil["time"] <- 1; epil2["time"] <- 4
|
|
epil2 <- rbind(epil, epil2)
|
|
epil2$pred <- unclass(epil2$trt) * (epil2$period > 0)
|
|
epil2$subject <- factor(epil2$subject)
|
|
epil3 <- aggregate(epil2, list(epil2$subject, epil2$period > 0),
|
|
function(x) if(is.numeric(x)) sum(x) else x[1])
|
|
epil3$pred <- factor(epil3$pred, labels = c("base", "placebo", "drug"))
|
|
|
|
contrasts(epil3$pred) <- structure(contr.sdif(3),
|
|
dimnames = list(NULL, c("placebo-base", "drug-placebo")))
|
|
summary(glm(y ~ pred + factor(subject) + offset(log(time)),
|
|
family = poisson, data = epil3))
|
|
|
|
glm(y ~ factor(subject), family = poisson, data = epil)
|
|
|
|
library(survival)
|
|
bacteria$Time <- rep(1, nrow(bacteria))
|
|
coxph(Surv(Time, unclass(y)) ~ week + strata(ID),
|
|
data = bacteria, method = "exact")
|
|
coxph(Surv(Time, unclass(y)) ~ factor(week) + strata(ID),
|
|
data = bacteria, method = "exact")
|
|
coxph(Surv(Time, unclass(y)) ~ I(week > 2) + strata(ID),
|
|
data = bacteria, method = "exact")
|
|
|
|
fit <- glm(y ~ trt + I(week> 2), binomial, data = bacteria)
|
|
summary(fit)
|
|
sum(residuals(fit, type = "pearson")^2)
|
|
|
|
if(FALSE) { # slow
|
|
## package available from http://www.stats.ox.ac.uk/pub/RWin
|
|
## but these examples fail
|
|
library(GLMMGibbs)
|
|
# declare a random intercept for each subject
|
|
epil$subject <- Ra(data = factor(epil$subject))
|
|
glmm(y ~ lbase*trt + lage + V4 + subject, family = poisson,
|
|
data = epil, keep = 100000, thin = 100)
|
|
|
|
epil3$subject <- Ra(data = factor(epil3$subject))
|
|
glmm(y ~ pred + subject, family = poisson,
|
|
data = epil3, keep = 100000, thin = 100)
|
|
}
|
|
|
|
|
|
summary(glmmPQL(y ~ trt + I(week> 2), random = ~ 1 | ID,
|
|
family = binomial, data = bacteria))
|
|
summary(glmmPQL(y ~ lbase*trt + lage + V4,
|
|
random = ~ 1 | subject,
|
|
family = poisson, data = epil))
|
|
summary(glmmPQL(y ~ pred, random = ~1 | subject,
|
|
family = poisson, data = epil3))
|
|
|
|
|
|
# 10.5 GEE models
|
|
|
|
## modified for YAGS >= 3.21-3
|
|
## package available from http://www.stats.ox.ac.uk/pub/RWin
|
|
if(require(yags)) {
|
|
print(yags(y == "y" ~ trt + I(week > 2), family = binomial, alphainit=0,
|
|
id = ID, corstr = "exchangeable", data = bacteria))
|
|
|
|
print(yags(y ~ lbase*trt + lage + V4, family = poisson, alphainit=0,
|
|
id = subject, corstr = "exchangeable", data = epil))
|
|
}
|
|
|
|
options(contrasts = c("contr.sum", "contr.poly"))
|
|
library(gee)
|
|
summary(gee(y ~ pred + factor(subject), family = poisson,
|
|
id = subject, data = epil3, corstr = "exchangeable"))
|
|
|
|
# End of ch10
|