Nothing
#-*- R -*-
# initialization
library(nlme)
library(lattice)
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 = "ch08.pdf")
# Chapter 8 Fitting Nonlinear Mixed-Effects Models
# 8.1 Fitting Nonlinear Models in S with nls and nlsList
## outer = ~1 is used to display all five curves in one panel
plot(Orange, outer = ~1)
logist <-
function(x, Asym, xmid, scal) Asym/(1 + exp(-(x - xmid)/scal))
logist <- deriv(~Asym/(1+exp(-(x-xmid)/scal)),
c("Asym", "xmid", "scal"), function(x, Asym, xmid, scal){})
Asym <- 180; xmid <- 700; scal <- 300
logist(Orange$age[1:7], Asym, xmid, scal)
fm1Oran.nls <- nls(circumference ~ logist(age, Asym, xmid, scal),
data = Orange, start = c(Asym = 170, xmid = 700, scal = 500))
summary(fm1Oran.nls)
plot(fm1Oran.nls)
plot(fm1Oran.nls, Tree ~ resid(.), abline = 0)
Orange.sortAvg <- sortedXyData("age", "circumference", Orange)
Orange.sortAvg
NLSstClosestX(Orange.sortAvg, 130)
logistInit <- function(mCall, LHS, data) {
xy <- sortedXyData(mCall[["x"]], LHS, data)
if(nrow(xy) < 3) {
stop("Too few distinct input values to fit a logistic")
}
Asym <- max(abs(xy[,"y"]))
if (Asym != max(xy[,"y"])) Asym <- -Asym # negative asymptote
xmid <- NLSstClosestX(xy, 0.5 * Asym)
scal <- NLSstClosestX(xy, 0.75 * Asym) - xmid
value <- c(Asym, xmid, scal)
names(value) <- mCall[c("Asym", "xmid", "scal")]
value
}
logist <- selfStart(logist, initial = logistInit)
class(logist)
logist <- selfStart(~ Asym/(1 + exp(-(x - xmid)/scal)),
initial = logistInit, parameters = c("Asym", "xmid", "scal"))
getInitial(circumference ~ logist(age, Asym, xmid, scal), Orange)
nls(circumference ~ logist(age, Asym, xmid, scal), Orange)
getInitial(circumference ~ SSlogis(age,Asym,xmid,scal), Orange)
nls(circumference ~ SSlogis(age, Asym, xmid, scal), Orange)
fm1Oran.lis <-
nlsList(circumference ~ SSlogis(age, Asym, xmid, scal) | Tree,
data = Orange)
fm1Oran.lis <- nlsList(SSlogis, Orange)
fm1Oran.lis.noSS <-
nlsList(circumference ~ Asym/(1+exp(-(age-xmid)/scal)),
data = Orange,
start = c(Asym = 170, xmid = 700, scal = 500))
fm1Oran.lis
summary(fm1Oran.lis)
plot(intervals(fm1Oran.lis), layout = c(3,1))
plot(fm1Oran.lis, Tree ~ resid(.), abline = 0)
Theoph[1:4,]
fm1Theo.lis <- nlsList(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
data = Theoph)
fm1Theo.lis
plot(intervals(fm1Theo.lis), layout = c(3,1))
pairs(fm1Theo.lis, id = 0.1)
# 8.2 Fitting Nonlinear Mixed-Effects Models with nlme
## no need to specify groups, as Orange is a groupedData object
## random is omitted - by default it is equal to fixed
(fm1Oran.nlme <-
nlme(circumference ~ SSlogis(age, Asym, xmid, scal),
data = Orange,
fixed = Asym + xmid + scal ~ 1,
start = fixef(fm1Oran.lis)))
summary(fm1Oran.nlme)
summary(fm1Oran.nls)
pairs(fm1Oran.nlme)
fm2Oran.nlme <- update(fm1Oran.nlme, random = Asym ~ 1)
anova(fm1Oran.nlme, fm2Oran.nlme)
plot(fm1Oran.nlme)
## level = 0:1 requests fixed (0) and within-group (1) predictions
plot(augPred(fm2Oran.nlme, level = 0:1),
layout = c(5,1))
qqnorm(fm2Oran.nlme, abline = c(0,1))
(fm1Theo.nlme <- nlme(fm1Theo.lis))
## IGNORE_RDIFF_BEGIN
try( intervals(fm1Theo.nlme, which="var-cov") ) ## could fail: Non-positive definite...
## IGNORE_RDIFF_END
(fm2Theo.nlme <- update(fm1Theo.nlme,
random = pdDiag(lKe + lKa + lCl ~ 1)))
fm3Theo.nlme <-
update(fm2Theo.nlme, random = pdDiag(lKa + lCl ~ 1))
anova(fm1Theo.nlme, fm3Theo.nlme, fm2Theo.nlme)
plot(fm3Theo.nlme)
qqnorm(fm3Theo.nlme, ~ ranef(.))
CO2
plot(CO2, outer = ~Treatment*Type, layout = c(4,1))
(fm1CO2.lis <- nlsList(SSasympOff, CO2))
## IGNORE_RDIFF_BEGIN
(fm1CO2.nlme <- nlme(fm1CO2.lis))
## IGNORE_RDIFF_END
(fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1))
anova(fm1CO2.nlme, fm2CO2.nlme)
plot(fm2CO2.nlme,id = 0.05,cex = 0.8,adj = -0.5)
fm2CO2.nlmeRE <- ranef(fm2CO2.nlme, augFrame = TRUE)
fm2CO2.nlmeRE
class(fm2CO2.nlmeRE)
plot(fm2CO2.nlmeRE, form = ~ Type * Treatment)
contrasts(CO2$Type)
contrasts(CO2$Treatment)
fm3CO2.nlme <- update(fm2CO2.nlme,
fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),
start = c(32.412, 0, 0, 0, -4.5603, 49.344))
summary(fm3CO2.nlme)
anova(fm3CO2.nlme, Terms = 2:4)
fm3CO2.nlmeRE <- ranef(fm3CO2.nlme, aug = TRUE)
plot(fm3CO2.nlmeRE, form = ~ Type * Treatment)
fm3CO2.fix <- fixef(fm3CO2.nlme)
fm4CO2.nlme <- update(fm3CO2.nlme,
fixed = list(Asym + lrc ~ Type * Treatment, c0 ~ 1),
start = c(fm3CO2.fix[1:5], 0, 0, 0, fm3CO2.fix[6]))
## IGNORE_RDIFF_BEGIN
summary(fm4CO2.nlme)
## IGNORE_RDIFF_END
fm5CO2.nlme <- update(fm4CO2.nlme, random = Asym ~ 1)
anova(fm4CO2.nlme, fm5CO2.nlme)
CO2$type <- 2 * (as.integer(CO2$Type) - 1.5)
CO2$treatment <- 2 * (as.integer(CO2$Treatment) - 1.5)
fm1CO2.nls <- nls(uptake ~ SSasympOff(conc, Asym.Intercept +
Asym.Type * type + Asym.Treatment * treatment +
Asym.TypeTreatment * type * treatment, lrc.Intercept +
lrc.Type * type + lrc.Treatment * treatment +
lrc.TypeTreatment * type * treatment, c0), data = CO2,
start = c(Asym.Intercept = 32.371, Asym.Type = -8.0086,
Asym.Treatment = -4.2001, Asym.TypeTreatment = -2.7253,
lrc.Intercept = -4.5267, lrc.Type = 0.13112,
lrc.Treatment = 0.093928, lrc.TypeTreatment = 0.17941,
c0 = 50.126))
anova(fm5CO2.nlme, fm1CO2.nls)
# plot(augPred(fm5CO2.nlme, level = 0:1), ## FIXME: problem with levels
# layout = c(6,2)) ## Actually a problem with contrasts.
## This fit just ping-pongs.
#fm1Quin.nlme <-
# nlme(conc ~ quinModel(Subject, time, conc, dose, interval,
# lV, lKa, lCl),
# data = Quinidine, fixed = lV + lKa + lCl ~ 1,
# random = pdDiag(lV + lCl ~ 1), groups = ~ Subject,
# start = list(fixed = c(5, -0.3, 2)),
# na.action = NULL, naPattern = ~ !is.na(conc), verbose = TRUE)
#fm1Quin.nlme
#fm1Quin.nlmeRE <- ranef(fm1Quin.nlme, aug = TRUE)
#fm1Quin.nlmeRE[1:3,]
# plot(fm1Quin.nlmeRE, form = lCl ~ Age + Smoke + Ethanol + ## FIXME: problem in max
# Weight + Race + Height + glyco + Creatinine + Heart,
# control = list(cex.axis = 0.7))
#fm1Quin.fix <- fixef(fm1Quin.nlme)
#fm2Quin.nlme <- update(fm1Quin.nlme,
# fixed = list(lCl ~ glyco, lKa + lV ~ 1),
# start = c(fm1Quin.fix[3], 0, fm1Quin.fix[2:1]))
fm2Quin.nlme <-
nlme(conc ~ quinModel(Subject, time, conc, dose, interval,
lV, lKa, lCl),
data = Quinidine, fixed = list(lCl ~ glyco, lV + lKa ~ 1),
random = pdDiag(diag(c(0.3,0.3)), form = lV + lCl ~ 1),
groups = ~ Subject,
start = list(fixed = c(2.5, 0, 5.4, -0.2)),
na.action = NULL, naPattern = ~ !is.na(conc))
summary(fm2Quin.nlme) # wrong values
options(contrasts = c("contr.treatment", "contr.poly"))
fm2Quin.fix <- fixef(fm2Quin.nlme)
## subsequent fits don't work
#fm3Quin.nlme <- update(fm2Quin.nlme,
# fixed = list(lCl ~ glyco + Creatinine, lKa + lV ~ 1),
# start = c(fm2Quin.fix[1:2], 0.2, fm2Quin.fix[3:4]))
#summary(fm3Quin.nlme)
#fm3Quin.fix <- fixef(fm3Quin.nlme)
#fm4Quin.nlme <- update(fm3Quin.nlme,
# fixed = list(lCl ~ glyco + Creatinine + Weight, lKa + lV ~ 1),
# start = c(fm3Quin.fix[1:3], 0, fm3Quin.fix[4:5]))
#summary(fm4Quin.nlme)
## This fit just ping-pongs
##fm1Wafer.nlmeR <-
## nlme(current ~ A + B * cos(4.5679 * voltage) +
## C * sin(4.5679 * voltage), data = Wafer,
## fixed = list(A ~ voltage + I(voltage^2), B + C ~ 1),
## random = list(Wafer = A ~ voltage + I(voltage^2),
## Site = pdBlocked(list(A~1, A~voltage+I(voltage^2)-1))),
### start = fixef(fm4Wafer), method = "REML", control = list(tolerance=1e-2))
## start = c(-4.255, 5.622, 1.258, -0.09555, 0.10434),
## method = "REML", control = list(tolerance = 1e-2))
##fm1Wafer.nlmeR
##fm1Wafer.nlme <- update(fm1Wafer.nlmeR, method = "ML")
(fm2Wafer.nlme <-
nlme(current ~ A + B * cos(w * voltage + pi/4),
data = Wafer,
fixed = list(A ~ voltage + I(voltage^2), B + w ~ 1),
random = list(Wafer = pdDiag(list(A ~ voltage + I(voltage^2), B + w ~ 1)),
Site = pdDiag(list(A ~ voltage+I(voltage^2), B ~ 1))),
start = c(-4.255, 5.622, 1.258, -0.09555, 4.5679)))
plot(fm2Wafer.nlme, resid(.) ~ voltage | Wafer,
panel = function(x, y, ...) {
panel.grid()
panel.xyplot(x, y)
panel.loess(x, y, lty = 2)
panel.abline(0, 0)
})
## anova(fm1Wafer.nlme, fm2Wafer.nlme, test = FALSE)
# intervals(fm2Wafer.nlme)
# 8.3 Extending the Basic nlme Model
#fm4Theo.nlme <- update(fm3Theo.nlme,
# weights = varConstPower(power = 0.1))
# this fit is way off
#fm4Theo.nlme
#anova(fm3Theo.nlme, fm4Theo.nlme)
#plot(fm4Theo.nlme)
## xlim used to hide an unusually high fitted value and enhance
## visualization of the heteroscedastic pattern
# plot(fm4Quin.nlme, xlim = c(0, 6.2))
#fm5Quin.nlme <- update(fm4Quin.nlme, weights = varPower())
#summary(fm5Quin.nlme)
#anova(fm4Quin.nlme, fm5Quin.nlme)
#plot(fm5Quin.nlme, xlim = c(0, 6.2))
var.nlme <- nlme(follicles ~ A + B * sin(2 * pi * w * Time) +
C * cos(2 * pi * w *Time), data = Ovary,
fixed = A + B + C + w ~ 1, random = pdDiag(A + B + w ~ 1),
# start = c(fixef(fm5Ovar.lme), 1))
start = c(12.18, -3.298, -0.862, 1))
##fm1Ovar.nlme
##ACF(fm1Ovar.nlme)
##plot(ACF(fm1Ovar.nlme, maxLag = 10), alpha = 0.05)
##fm2Ovar.nlme <- update(fm1Ovar.nlme, correlation = corAR1(0.311))
##fm3Ovar.nlme <- update(fm1Ovar.nlme, correlation = corARMA(p=0, q=2))
##anova(fm2Ovar.nlme, fm3Ovar.nlme, test = FALSE)
##intervals(fm2Ovar.nlme)
##fm4Ovar.nlme <- update(fm2Ovar.nlme, random = A ~ 1)
##anova(fm2Ovar.nlme, fm4Ovar.nlme)
##if (interactive()) fm5Ovar.nlme <- update(fm4Ovar.nlme, correlation = corARMA(p=1, q=1))
# anova(fm4Ovar.nlme, fm5Ovar.nlme)
# plot(ACF(fm5Ovar.nlme, maxLag = 10, resType = "n"),
# alpha = 0.05)
# fm5Ovar.lmeML <- update(fm5Ovar.lme, method = "ML")
# intervals(fm5Ovar.lmeML)
# fm6Ovar.lmeML <- update(fm5Ovar.lmeML, random = ~1)
# anova(fm5Ovar.lmeML, fm6Ovar.lmeML)
# anova(fm6Ovar.lmeML, fm5Ovar.nlme)
# intervals(fm5Ovar.nlme, which = "fixed")
fm1Dial.lis <-
nlsList(rate ~ SSasympOff(pressure, Asym, lrc, c0) | QB,
data = Dialyzer)
fm1Dial.lis
plot(intervals(fm1Dial.lis))
fm1Dial.gnls <- gnls(rate ~ SSasympOff(pressure, Asym, lrc, c0),
data = Dialyzer, params = list(Asym + lrc ~ QB, c0 ~ 1),
start = c(53.6, 8.6, 0.51, -0.26, 0.225))
fm1Dial.gnls
Dialyzer$QBcontr <- 2 * (Dialyzer$QB == 300) - 1
fm1Dial.nls <-
nls(rate ~ SSasympOff(pressure, Asym.Int + Asym.QB * QBcontr,
lrc.Int + lrc.QB * QBcontr, c0), data = Dialyzer,
start = c(Asym.Int = 53.6, Asym.QB = 8.6, lrc.Int = 0.51,
lrc.QB = -0.26, c0 = 0.225))
## IGNORE_RDIFF_BEGIN
summary(fm1Dial.nls)
## IGNORE_RDIFF_END
logLik(fm1Dial.nls)
plot(fm1Dial.gnls, resid(.) ~ pressure, abline = 0)
fm2Dial.gnls <- update(fm1Dial.gnls,
weights = varPower(form = ~ pressure))
anova(fm1Dial.gnls, fm2Dial.gnls)
ACF(fm2Dial.gnls, form = ~ 1 | Subject)
plot(ACF(fm2Dial.gnls, form = ~ 1 | Subject), alpha = 0.05)
fm3Dial.gnls <-
update(fm2Dial.gnls, corr = corAR1(0.716, form = ~ 1 | Subject))
fm3Dial.gnls
intervals(fm3Dial.gnls)
anova(fm2Dial.gnls, fm3Dial.gnls)
# restore two fitted models
fm2Dial.lme <-
lme(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,
Dialyzer, ~ pressure + I(pressure^2),
weights = varPower(form = ~ pressure))
fm2Dial.lmeML <- update(fm2Dial.lme, method = "ML")
fm3Dial.gls <-
gls(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,
Dialyzer, weights = varPower(form = ~ pressure),
corr = corAR1(0.771, form = ~ 1 | Subject))
fm3Dial.glsML <- update(fm3Dial.gls, method = "ML")
anova( fm2Dial.lmeML, fm3Dial.glsML, fm3Dial.gnls, test = FALSE)
# cleanup
summary(warnings())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.