#-*- 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())