## BMD
bmd <- read.delim("../../data/bmd.txt")

## Slide 12
ll <- seq(0.2,1.5,by=.01)
criteria <- array(NA,dim=c(length(ll),2,2),dimnames=list(ll,c("male","female"),c("CV","GCV")))
deg <- numeric(length(ll))
for (j in 1:2)
  {
    Data <- bmd[bmd$gender==dimnames(criteria)[[2]][j],]
    x <- Data$age
    y <- Data$spnbmd
    for (i in 1:length(ll))
      {
        fit <- smooth.spline(x,y,spar=ll[i])
        l <- fit$lev[match(x,fit$x)]
        yhat <- predict(fit,x)$y
        criteria[i,j,1] <- mean(((y-yhat)/(1-l))^2)
        criteria[i,j,2] <- mean(((y-yhat)/(1-mean(fit$lev)))^2)
        deg[i] <- sum(fit$lev)
      }
  }
dimnames(criteria)[[1]] <- deg
df <- array2df(criteria,vars=c("ll","gender","method","criterion"))
df$ll <- factor2num(df$ll)
trellis.par.set(superpose.line=list(col=c("black","blue")))
print(xyplot(criterion~ll|gender,df,group=method,type="l",auto.key=list(points=FALSE,lines=TRUE,columns=2),xlab="Degrees of freedom",ylab="",scale=list(x=list(log=2)),xscale.components=xscale.components.log))

## Slide 13
ll <- c(0.45,NA,1.15)
par(mfrow=c(2,3),mar=c(2,2,3,0))
main <- rbind(c("","Males",""),c("","Females",""))
for (j in 1:2)
  {
    Data <- bmd[bmd$gender==dimnames(criteria)[[2]][j],]
    x <- Data$age
    y <- Data$spnbmd
    xx <- seq(min(x),max(x),length=101)
    for (i in 1:3)
      {
        if (i==2) fit <- smooth.spline(x,y)
        else fit <- smooth.spline(x,y,spar=ll[i])
        yhat <- predict(fit,xx)$y
        plot(x,y,col="gray85",xlab="",ylab="",main=main[j,i],pch=19)
        lines(xx,yhat,lwd=2,col="red")
      }
  }

## Slide 17
require(mgcv)
myPanel <- function(...)
  {
    panel.xyplot(...,alpha=.3)
    panel.smooth(...)
  }
panel.smooth <- function(x,y,...)
  {
    fit <- gam(y~s(x))
    xx <- seq(min(x),max(x),len=101)
    yy <- predict(fit,newdata=data.frame(x=xx),se.fit=T)
    lpolygon(c(xx,rev(xx)),c(yy$fit-1.96*yy$se.fit,rev(yy$fit+1.96*yy$se.fit)),col=rgb(.6,.6,.6,alpha=.4),border=F,...)
    llines(xx,yy$fit,col="black",lwd=2)
  }
trellis.par.set(plot.symbol=list(pch=19))
xyplot(spnbmd~age|gender,bmd,panel=myPanel)