# Chapter 9

library(MASS, first=T)

attach(wtloss)
# alter margin 4; others are default
par(mar=c(5.1, 4.1, 4.1, 4.1))
plot(Days, Weight, type="p",ylab="Weight (kg)")
Wt.lbs <- pretty(range(Weight*2.205))
axis(side=4, at=Wt.lbs/2.205, lab=Wt.lbs, srt=90)
mtext("Weight (lb)", side=4, line=3)

# 9.1

wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th),
    data = wtloss,
    start = list(b0=90, b1=95, th=120),
    trace = T)
wtloss.fm

# 9.2

parameters(wtloss) <- list(b0=90, b1=95, th=120)
wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th),
    data = wtloss, trace=T)
#param(wtloss, "b0") <- 80

# 9.3

expn <- function(b0, b1, th, x)
{
    temp <- 2^(-x/th)
    model.func <- b0 + b1 * temp
    Z <- cbind(1, temp, (b1 * x * temp * log(2))/th^2)
    dimnames(Z) <- list(NULL, c("b0","b1","th"))
    attr(model.func, "gradient") <- Z
    model.func
}

wtloss.gr <- nls(Weight ~ expn(b0, b1, th, Days),
    data = wtloss, trace = T)
expn1 <- deriv(y ~ b0 + b1 * 2^(-x/th), c("b0", "b1", "th"),
               function(b0, b1, th, x) NULL)

# 9.4

summary(wtloss.gr)
deviance(wtloss.gr)
vcov(wtloss.gr)

# 9.5

parameters(wtloss) <- list(th=120)
wtloss.pl <- nls(Weight ~ cbind(1, 2^(-Days/th)),
             data=wtloss, algorithm="plinear", trace=T)
summary(wtloss.pl)

# 9.6

time.pds <- deriv(~ -th * log((y0 - b0)/b1)/log(2),
    c("b0","b1","th"), function(y0, b0, b1, th) NULL)
est.time <- function(y0, obj) {
    b <- coef(obj)
    tmp <- time.pds(y0, b["b0"], b["b1"], b["th"])
    x0 <- as.vector(tmp)
    lam <- attr(tmp, "gradient")
    v <- (lam %*% vcov(obj) * lam) %*% matrix(1, 3, 1)
    a <- cbind(x0, sqrt(v))
    dimnames(a) <- list(paste(y0, "kg: "), c("x0", "SE"))
    a
}
est.time(c(100, 90, 85), wtloss.gr)

attach(wtloss)
plot(Days, Weight, xlab= "days", ylab ="weight (kg)",
   xlim=c(0,730), ylim=c(70, 200))
xx <- seq(0, 730, 10)
lines(xx, 81.37+ 102.68 * 2^(-xx/141.91))
wtloss.quad <- lm(Weight ~ poly(Days, 2))
lines(xx, predict.gam(wtloss.quad, data.frame(Days=xx)), lty=2)
wtloss.cub <- lm(Weight ~ poly(Days, 3))
lines(xx, predict.gam(wtloss.cub, data.frame(Days=xx)), lty=3)
legend(locator(1), c("exponential", "quadratic", "cubic"),
    lty=1:3)
detach()

fm0 <- lm(Wt*Time ~ Viscosity + Time - 1,  data=stormer)
b0 <- as.list(coef(fm0))
names(b0) <- paste("b",1:2,sep="")
unlist(b0)
parameters(stormer) <- b0
storm.fm <- nls(Time ~ b1*Viscosity/(Wt-b2), data=stormer,
          trace=T)
attach(stormer,1)
z1 <- Wt/Viscosity
z2 <- -1/Viscosity
detach(1,save="stormer")
attach(stormer)
storm.gm <- glm(Time ~ z1 + z2 - 1,
     family=quasi(link=inverse, variance=constant),
     data=stormer, trace=T)
g <- coef(storm.gm)
b <- coef(storm.fm)
b0 <- c(1/g[1], g[2]/g[1])
cbind(b,b0)

J <- matrix(c(-1/g[1]^2, -g[2]/g[1]^2, 0, 1/g[1]), 2, 2)
J %*% vcov(storm.gm) %*% t(J)
vcov(storm.fm)
bc <- coef(storm.fm)
se <- sqrt(diag(vcov(storm.fm)))
dv <- deviance(storm.fm)
par(pty = "s")
b1 <- bc[1] + seq(-3*se[1], 3*se[1], length = 51)
b2 <- bc[2] + seq(-3*se[2], 3*se[2], length = 51)
bv <- expand.grid(b1, b2)
ssq <- function(b)
       sum((Time - b[1] * Viscosity/(Wt-b[2]))^2)
dbeta <- apply(bv, 1, ssq)
cc <- matrix(Time - rep(bv[,1],rep(23,2601)) *
       Viscosity/(Wt - rep(bv[,2], rep(23,2601))), 23)
dbeta <- matrix(drop(rep(1, 23) %*% cc^2),51)
fstat <- matrix( ((dbeta - dv)/2) / (dv/21), 51, 51)
qf(0.95, 2, 21)
plot(b1, b2, type="n")
lev <- c(1,2,5,7,10,15,20)
contour(b1, b2, fstat, levels=lev, labex=0.75, lty=2, add=T)
contour(b1, b2, fstat, levels=qf(0.95,2,21), add=T, labex=0)
text(31.6,0.3,"95% CR", adj=0, cex=0.75)
points(bc[1], bc[2], pch=3, mkh=0.1)

# 9.7

plist <- enlist(coef(wtloss.gr))
expn2 <- deriv3(~ b0 + b1*2^(-x/th), c("b0","b1","th"),
      function(x, b0, b1, th) NULL)
wtloss.he <- nls(Weight ~ expn2(Days, b0, b1, th),
      start=plist, data=wtloss)
rms.curv(wtloss.he)
stab <- deriv3(~ ((u0*u2-u1^2) +
    (u0-u1)^2 *((u1-u2)/(u0-u1))^((x-40)/80))/(u0-2*u1+u2),
    c("u0","u1","u2"), function(x, u0, u1, u2) NULL)
mu <- predict(lm(Weight~Days+Days^2, data=wtloss),
    newdata=data.frame(Days=c(40,120,200)))
names(mu) <- paste("u", 0:2, sep="")
wtloss.st <- nls(Weight ~ stab(Days, u0, u1, u2),
    start=mu, data=wtloss, trace=T)
rms.curv(wtloss.st)
summary(wtloss.st)$correlation
par(mfrow=c(2,3), pty="s")
plot(profile(wtloss.gr), ask=F)
plot(profile(wtloss.st))

# 9.8

attach(faithful)
hist(waiting, prob=T, xlim=c(35,105), ylim=c(0,0.045),
   style.bar="old")
wait.dns <- density(waiting, 200, width=10)
lines(wait.dns, lty=2)
lmix2 <- deriv3(
     ~ -log(p*dnorm((x-u1)/s1)/s1 + (1-p)*dnorm((x-u2)/s2)/s2),
     c("p", "u1", "s1", "u2", "s2"),
     function(x, p, u1, s1, u2, s2) NULL)
p0 <- list(p=sum(waiting < 70)/length(waiting),
           u1=50, s1=5, u2=80, s2=5)
unlist(p0)
tr.ms <- function(info, theta, grad, scale, flags, fit.pars)
{
    cat(signif(info[3]), ":", signif(theta), "\n")
    invisible(list())
}
wait.mix2 <- ms(~ lmix2(waiting, p, u1, s1, u2, s2),
    start=p0, data = faithful, trace = tr.ms)
dmix2 <- function(x, p, u1, s1, u2, s2)
             p * dnorm(x, u1, s1) + (1-p) * dnorm(x, u2, s2)
cf <- enlist(coef(wait.mix2))
wait.fdns <- list(x = wait.dns$x, y =
    eval(expression(dmix2(wait.dns$x, p, u1, s1, u2, s2)),
    local = cf))
lines(wait.fdns)
par(usr = c(0,1,0,1))
legend(0.1, 0.9, c("Normal mixture", "Nonparametric"),
    lty = c(1,2), bty = "n")
pmix2 <- deriv(~ p*pnorm((x-u1)/s1) + (1-p)*pnorm((x-u2)/s2),
      "x", function(x, p, u1, s1, u2, s2) NULL)
pr0 <- (1:272 - 0.5)/272
x0 <- x1 <- sort(waiting)
for(i in 1:10) {
    pr <- eval(expression(pmix2(x0, p, u1, s1, u2, s2)),
               local = cf)
    del <- (pr - pr0)/attr(pr, "gradient")
    x0 <- x0 - del
    print(del <- max(abs(del)))
    if(del < 0.0005) break
}
par(pty = "s")
plot(x0, x1, xlim = range(x0, x1), ylim = range(x0, x1),
    xlab = "Model quantiles", ylab = "Waiting time")
abline(0,1)
vmat <- wait.mix2$hessian        # lower triangle only
vmat <- solve(vmat + t(vmat) - diag(diag(vmat)))
# S-Plus 3.2 users can use the inverse information
# vmat <- summary(wait.mix2)$Information
cbind(coef(wait.mix2), sqrt(diag(vmat)))
(cf$s1 - cf$s2) /sqrt(c(0,0,1,0,1) %*% vmat %*% c(0,0,1,0,1))
wait.mix2e <- ms(~ -log(p * dnorm(waiting, u1, s1) +
   (1-p) * dnorm(waiting, u2, s1)), start = cf[-5],
   data = faithful)
coef(wait.mix2e)
2*(wait.mix2e$value - wait.mix2$value)
detach()