# packrat/lib-R/MASS/scripts/ch10.R In UBC-MDS/Karl:

```#-*- 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
```
UBC-MDS/Karl documentation built on May 22, 2019, 1:53 p.m.