116 lines
4.1 KiB
R
116 lines
4.1 KiB
R
|
#-*- R -*-
|
||
|
|
|
||
|
|
# initialization
|
||
|
|
|
||
|
|
library(nlme)
|
||
|
|
options(width = 65,
|
||
|
|
## reduce platform dependence in printed output when testing
|
||
|
|
digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)
|
||
|
|
options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))
|
||
|
|
pdf(file = "ch06.pdf")
|
||
|
|
|
||
|
|
# Chapter 6 Nonlinear Mixed-Effects Models:
|
||
|
|
# Basic Concepts and Motivating Examples
|
||
|
|
|
||
|
|
# 6.2 Indomethicin Kinetics
|
||
|
|
|
||
|
|
plot(Indometh)
|
||
|
|
fm1Indom.nls <- nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2),
|
||
|
|
data = Indometh)
|
||
|
|
summary(fm1Indom.nls)
|
||
|
|
plot(fm1Indom.nls, Subject ~ resid(.), abline = 0)
|
||
|
|
(fm1Indom.lis <- nlsList(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2),
|
||
|
|
data = Indometh))
|
||
|
|
plot(intervals(fm1Indom.lis))
|
||
|
|
## IGNORE_RDIFF_BEGIN
|
||
|
|
(fm1Indom.nlme <- nlme(fm1Indom.lis,
|
||
|
|
random = pdDiag(A1 + lrc1 + A2 + lrc2 ~ 1),
|
||
|
|
control = list(tolerance = 0.0001)))
|
||
|
|
## IGNORE_RDIFF_END
|
||
|
|
fm2Indom.nlme <- update(fm1Indom.nlme,
|
||
|
|
random = pdDiag(A1 + lrc1 + A2 ~ 1))
|
||
|
|
anova(fm1Indom.nlme, fm2Indom.nlme)
|
||
|
|
(fm3Indom.nlme <- update(fm2Indom.nlme, random = A1+lrc1+A2 ~ 1))
|
||
|
|
fm4Indom.nlme <-
|
||
|
|
update(fm3Indom.nlme,
|
||
|
|
random = pdBlocked(list(A1 + lrc1 ~ 1, A2 ~ 1)))
|
||
|
|
## IGNORE_RDIFF_BEGIN
|
||
|
|
anova(fm3Indom.nlme, fm4Indom.nlme)
|
||
|
|
## IGNORE_RDIFF_END
|
||
|
|
anova(fm2Indom.nlme, fm4Indom.nlme)
|
||
|
|
plot(fm4Indom.nlme, id = 0.05, adj = -1)
|
||
|
|
qqnorm(fm4Indom.nlme)
|
||
|
|
plot(augPred(fm4Indom.nlme, level = 0:1))
|
||
|
|
summary(fm4Indom.nlme)
|
||
|
|
|
||
|
|
# 6.3 Growth of Soybean Plants
|
||
|
|
|
||
|
|
head(Soybean)
|
||
|
|
plot(Soybean, outer = ~ Year * Variety)
|
||
|
|
(fm1Soy.lis <- nlsList(weight ~ SSlogis(Time, Asym, xmid, scal),
|
||
|
|
data = Soybean,
|
||
|
|
## in R >= 3.4.3, more iterations are needed for "1989P5"
|
||
|
|
## due to a change of initial values in SSlogis();
|
||
|
|
## control is passed to getInitial() only since R 4.1.0
|
||
|
|
control = list(maxiter = 60)))
|
||
|
|
## IGNORE_RDIFF_BEGIN
|
||
|
|
(fm1Soy.nlme <- nlme(fm1Soy.lis))
|
||
|
|
## IGNORE_RDIFF_END
|
||
|
|
fm2Soy.nlme <- update(fm1Soy.nlme, weights = varPower())
|
||
|
|
anova(fm1Soy.nlme, fm2Soy.nlme)
|
||
|
|
plot(ranef(fm2Soy.nlme, augFrame = TRUE),
|
||
|
|
form = ~ Year * Variety, layout = c(3,1))
|
||
|
|
soyFix <- fixef(fm2Soy.nlme)
|
||
|
|
options(contrasts = c("contr.treatment", "contr.poly"))
|
||
|
|
## IGNORE_RDIFF_BEGIN
|
||
|
|
(fm3Soy.nlme <-
|
||
|
|
update(fm2Soy.nlme,
|
||
|
|
fixed = Asym + xmid + scal ~ Year,
|
||
|
|
start = c(soyFix[1], 0, 0, soyFix[2], 0, 0, soyFix[3], 0, 0)))
|
||
|
|
## IGNORE_RDIFF_END
|
||
|
|
anova(fm3Soy.nlme)
|
||
|
|
# The following line is not in the book but needed to fit the model
|
||
|
|
fm4Soy.nlme <-
|
||
|
|
nlme(weight ~ SSlogis(Time, Asym, xmid, scal),
|
||
|
|
data = Soybean,
|
||
|
|
fixed = list(Asym ~ Year*Variety, xmid ~ Year + Variety, scal ~ Year),
|
||
|
|
random = Asym ~ 1,
|
||
|
|
start = c(17, 0, 0, 0, 0, 0, 52, 0, 0, 0, 7.5, 0, 0),
|
||
|
|
weights = varPower(0.95), control = list(verbose = TRUE))
|
||
|
|
# FIXME: An update doesn't work for the fixed argument when fixed is a list
|
||
|
|
## p. 293-4 :
|
||
|
|
summary(fm4Soy.nlme)
|
||
|
|
plot(augPred(fm4Soy.nlme))# Fig 6.14, p. 295
|
||
|
|
|
||
|
|
# 6.4 Clinical Study of Phenobarbital Kinetics
|
||
|
|
|
||
|
|
(fm1Pheno.nlme <-
|
||
|
|
nlme(conc ~ phenoModel(Subject, time, dose, lCl, lV),
|
||
|
|
data = Phenobarb, fixed = lCl + lV ~ 1,
|
||
|
|
random = pdDiag(lCl + lV ~ 1), start = c(-5, 0),
|
||
|
|
na.action = NULL, naPattern = ~ !is.na(conc)))
|
||
|
|
fm1Pheno.ranef <- ranef(fm1Pheno.nlme, augFrame = TRUE)
|
||
|
|
# (These plots used to encounter difficulties, now fine):
|
||
|
|
plot(fm1Pheno.ranef, form = lCl ~ Wt + ApgarInd)
|
||
|
|
plot(fm1Pheno.ranef, form = lV ~ Wt + ApgarInd)
|
||
|
|
|
||
|
|
options(contrasts = c("contr.treatment", "contr.poly"))
|
||
|
|
if(FALSE)## This fit just "ping-pongs" until max.iterations error
|
||
|
|
fm2Pheno.nlme <-
|
||
|
|
update(fm1Pheno.nlme,
|
||
|
|
fixed = list(lCl ~ Wt, lV ~ Wt + ApgarInd),
|
||
|
|
start = c(-5.0935, 0, 0.34259, 0, 0),
|
||
|
|
control = list(pnlsTol = 1e-4, maxIter = 500,
|
||
|
|
msVerbose = TRUE, opt = "nlm"))
|
||
|
|
##summary(fm2Pheno.nlme)
|
||
|
|
##fm3Pheno.nlme <-
|
||
|
|
## update(fm2Pheno.nlme,
|
||
|
|
## fixed = lCl + lV ~ Wt,
|
||
|
|
## start = fixef(fm2Pheno.nlme)[-5])
|
||
|
|
##plot(fm3Pheno.nlme, conc ~ fitted(.), abline = c(0,1))
|
||
|
|
|
||
|
|
# cleanup
|
||
|
|
|
||
|
|
summary(warnings())
|
||
|
|
|