# load weight loss data
data(wtloss,package=MASS)

# estimate theta~
tmp.fun <- function(theta) {
    theta <- theta/100
    ssanova(Weight~exp(-theta*Days),data=wtloss)$score
}
nlm(tmp.fun,1)$estimate
# gcv values of L-spline fit and cubic spline fit
ssanova(Weight~exp(-.004885*Days),data=wtloss)$score
ssanova(Weight~Days,data=wtloss)$score

# calculate and evaluate L-spline fit
wtloss$dd <- exp(-.004885*wtloss$Days)
wtloss.fit1 <- ssanova(Weight~dd,data=wtloss)
tt <- seq(0,250,length=101)
est1 <- predict(wtloss.fit1,data.frame(dd=exp(-.004885*tt)),se=TRUE)
# calculate cubic spline fit
wtloss.fit2 <- ssanova(Weight~Days,data=wtloss)
est2 <- predict(wtloss.fit2,data.frame(Days=tt),se=TRUE)

# draw figure 4.3 left frame
est0 <- 81.374+102.68*2^(-tt/141.91)
plot(wtloss$Days,wtloss$Weight,xlab="days",ylab="weight (kg)")
lines(tt,est0,col=6)
lines(tt,est1$fit,col=4)
lines(tt,est2$fit,col=5)

# draw figure 4.3 right frame
plot(tt,est1$fit-est0,type="l",ylim=c(-1.5,1.5),xlab="days",ylab="")
lines(tt,est1$fit-est0-1.96*est1$se,lty=2)
lines(tt,est1$fit-est0+1.96*est1$se,lty=2)
lines(tt,est2$fit-est0,col=5)
lines(tt,est2$fit-est0-1.96*est2$se,col=5,lty=2)
lines(tt,est2$fit-est0+1.96*est2$se,col=5,lty=2)


# generate simulation data
set.seed(5732)
tt <- ((1:100)-.5)/100
yy <- 5+3*exp(-4*tt)+2*exp(-8*tt)+.5*rnorm(tt)

# estimate theta~
tmp.fun <- function(theta) {
    ssanova(yy~exp(-theta*tt))$score
}
nlm(tmp.fun,4)$estimate
# gcv values of L-spline fit and cubic spline fit
ssanova(yy~exp(-5.2559*tt))$score
ssanova(yy~tt)$score

# calculate and evaluate L-spline and cubic spline fit
ttt <- exp(-5.2559*tt)
fit.L <- ssanova(yy~ttt)
est.L <- predict(fit.L,data.frame(ttt=ttt),se=TRUE)
fit.c <- ssanova(yy~tt)
est.c <- predict(fit.c,data.frame(tt=tt),se=TRUE)

# draw figure 4.4 left frame
plot(tt,yy,xlab="x",ylab="y")
lines(tt,5+3*exp(-4*tt)+2*exp(-8*tt),col=6)
lines(tt,est.L$fit,col=4)
lines(tt,est.L$fit-1.96*est.L$se,col=4,lty=2)
lines(tt,est.L$fit+1.96*est.L$se,col=4,lty=2)
lines(tt,est.c$fit,col=5)
lines(tt,est.c$fit-1.96*est.c$se,col=5,lty=2)
lines(tt,est.c$fit+1.96*est.c$se,col=5,lty=2)

# draw figure 4.4 right frame
est0 <- 5+3*exp(-4*tt)+2*exp(-8*tt)
plot(tt,est.L$fit-est0,type="l",ylim=c(-.75,.75),
xlab="x",ylab="")
lines(c(0,1),c(0,0),col=6)
lines(tt,est.L$fit-est0-1.96*est.L$se,lty=2)
lines(tt,est.L$fit-est0+1.96*est.L$se,lty=2)
lines(tt,est.c$fit-est0,col=5)
lines(tt,est.c$fit-est0-1.96*est.c$se,col=5,lty=2)
lines(tt,est.c$fit-est0+1.96*est.c$se,col=5,lty=2)