## General setup
require(lattice)
require(splines)
require(mgcv)
heart <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/621/data/heart.txt")
bmd <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/621/data/bmd.txt")
m <- subset(bmd, gender=="male")
f <- subset(bmd, gender=="female")

## Slide 5
form <- list(formula(y~x+I(x^2)),formula(y~x+I(x^2)+I(x^3)+I(x^4)),formula(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6)))
for (J in 1:4) {
  f <- function(x){dnorm(x,0.5,0.18)}
  main <- list(expression("Up to "*x^2),expression("Up to "*x^4),expression("Up to "*x^6))
  N <- 50
  n <- 100
  xx <- seq(0,1,len=101)
  par(mfrow=c(1,3))
  for (j in 1:3) {
    plot(xx,f(xx),col="red",lwd=3,type="l",ylim=c(-1,3),xlab="x",ylab="f(x)",main=main[[j]])
    for (i in 1:N) {
      x <- runif(n)
      y <- rnorm(n,f(x))
      fit <- lm(form[[j]])
      lines(xx,predict(fit,data.frame(x=xx)),col=rgb(.7,.7,.7,alpha=.5))
    }
    lines(xx,f(xx),col="red",lwd=3)
  }
}

## Building up to splines
## Slide 8
myPanel <- function(x, y, ...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  c1 <- mean(y[x < knots[1]])
  c2 <- mean(y[x >= knots[1] & x < knots[2]])
  c3 <- mean(y[x >= knots[2]])
  x1 <- seq(min(x),knots[1],len=33)
  x2 <- seq(knots[1],knots[2],len=33)
  x3 <- seq(knots[2],max(x),len=33)
  llines(x1, c1, col="blue", lwd=3)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y),col="black", lty=2)
  llines(x2, c2, col="blue", lwd=3)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y),col="black", lty=2)
  llines(x3, c3, col="blue", lwd=3)
}
xyplot(spnbmd~age|gender, bmd, panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 9
myPanel <- function(x, y, ...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  ind1 <- x < knots[1]
  ind2 <- x >= knots[1] & x < knots[2]
  ind3 <- x >= knots[2]
  x1 <- x[ind1]
  x2 <- x[ind2]
  x3 <- x[ind3]
  fit1 <- lm(y[ind1]~x1)
  fit2 <- lm(y[ind2]~x2)
  fit3 <- lm(y[ind3]~x3)
  x1 <- seq(min(x),knots[1],len=33)
  x2 <- seq(knots[1],knots[2],len=33)
  x3 <- seq(knots[2],max(x),len=33)
  llines(x1,predict(fit1,data.frame(x1=x1)), col="blue", lwd=3)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y), col="black", lty=2)
  llines(x2,predict(fit2,data.frame(x2=x2)),col="blue", lwd=3)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y), col="black", lty=2)
  llines(x3,predict(fit3,data.frame(x3=x3)),col="blue", lwd=3)
}
xyplot(spnbmd~age|gender, bmd, panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 10
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  fit <- lm(y~bs(x,knots=knots,degree=1))
  x0 <- seq(min(x),max(x),len=101)
  llines(x0,predict(fit,data.frame(x=x0)),col="blue", lwd=3)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y), col="black", lty=2)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y), col="black", lty=2)
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 15
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  fit <- lm(y~bs(x,knots=knots,degree=2))
  x0 <- seq(min(x),max(x),len=101)
  llines(x0,predict(fit,data.frame(x=x0)),col="blue", lwd=3)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y),col="black", lty=2)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y),col="black", lty=2)
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 16
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  fit <- lm(y~bs(x,knots=knots,degree=3))
  x0 <- seq(min(x),max(x),len=101)
  llines(x0,predict(fit,data.frame(x=x0)),col="blue", lwd=3)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y),col="black", lty=2)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y),col="black", lty=2)
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 21
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  knots=quantile(x,c(1/3,2/3))
  fit <- lm(y~ns(x,knots=knots))
  x0 <- seq(min(x),max(x),len=101)
  llines(x0,predict(fit,data.frame(x=x0)),col="blue", lwd=3)
  lsegments(x0=min(x),y0=min(y),x1=min(x),y1=max(y),col="black", lty=2)
  lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y),col="black", lty=2)
  lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y),col="black", lty=2)
  lsegments(x0=max(x),y0=min(y),x1=max(x),y1=max(y),col="black", lty=2)
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 22
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)
  fit <- lm(y~ns(x,df=5))
  x0 <- seq(min(x),max(x),len=101)
  llines(x0,predict(fit,data.frame(x=x0)),col="blue", lwd=3)
  knots <- c(min(x),attr(ns(x,df=5),"knots"),max(x))
  for (i in 1:length(knots)) {
    lsegments(x0=knots[i],y0=min(y),x1=knots[i],y1=max(y),col="black", lty=2)
  }
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 23
par(mfrow=c(1,1))
pcol <- c("#FF4E3780", "#008DFF80")
lcol <- c("#FF4E37FF", "#008DFFFF")
with(f, plot(age, spnbmd, col=pcol[1], pch=16, xlab="Age", ylab="Relative change in spinal BMD"))
with(m, points(age, spnbmd, col=pcol[2], pch=16))
fit.f <- lm(spnbmd~ns(age, df=5), f)
fit.m <- lm(spnbmd~ns(age, df=5), m)
age <- seq(min(bmd$age),max(bmd$age),len=101)
Y <- cbind(predict(fit.f, data.frame(age=age)), predict(fit.m, data.frame(age=age)))
matlines(age, Y, col=lcol, lwd=3, lty=1)
legend("topright", legend=c("Female","Male"), col=lcol, lty=1, lwd=3, ncol=2)

## Slide 24
myPanel <- function(x,y,...) {
  panel.xyplot(x, y, ..., alpha=.5, pch=16, col="gray", cex=0.8)  
  fit.spline <- lm(y~ns(x,df=5))
  fit.loess <- loess(y~x,span=.61,degree=2)
  x0 <- seq(min(x), max(x), len=101)
  llines(x0, predict(fit.loess, data.frame(x=x0)), col=lcol[1], lwd=3)
  llines(x0, predict(fit.spline, data.frame(x=x0)), col=lcol[2], lwd=3)
}
trellis.par.set(superpose.line=list(lwd=3, col=lcol))
xyplot(spnbmd~age|gender, bmd, panel=myPanel, auto.key=list(columns=2,points=FALSE,lines=TRUE,text=c("loess","spline")), xlab="Age", ylab="Relative change in spinal BMD")

## Slide 27
myPanel <- function(x,y,...) {
  fit <- lm(y~ns(x,df=5))
  x0 <- seq(min(x),max(x),len=101)
  y0 <- predict(fit,newdata=data.frame(x=x0),se.fit=T)
  lpolygon(c(x0,rev(x0)), c(y0$fit-2*y0$se.fit,rev(y0$fit+2*y0$se.fit)), col="gray75", border=F,...)
  lpoints(x, y, col="gray30", cex=0.5, pch=16)
  llines(x0, y0$fit, col="blue", lwd=3)
}
xyplot(spnbmd~age|gender,bmd,panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")

## Slide 28
par(mar=c(5,5,1,1),mfrow=c(1,2))
fit <- glm(chd~ns(sbp,df=3), data=heart, family="binomial")
visreg(fit, partial=FALSE, xlab="Systolic blood pressure", ylab=expression(hat(f)(x)))
visreg(fit, partial=FALSE, scale="response", xlab="Systolic blood pressure", ylab=expression(hat(pi)(x)))

## CV, GCV for splines
Data <- list(f,m)
ll <- seq(0.2,1.5,by=.01)
criteria <- array(NA,dim=c(length(ll),2,2),dimnames=list(ll,c("female","male"),c("CV","GCV")))
deg <- numeric(length(ll))
for (j in 1:2) {
  x <- Data[[j]]$age
  y <- Data[[j]]$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)
col <- c("#FF4E37FF", "#008DFFFF")
trellis.par.set(superpose.line=list(col=col, lwd=3))
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))
matrix(deg[apply(criteria,2:3,which.min)], ncol=2, dimnames=dimnames(criteria)[2:3])

## Slide 40
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")
  }
}

## Confidence bands
myPanel <- function(x,y,...) {
  fit <- gam(y~s(x))
  x0 <- seq(min(x),max(x),len=101)
  y0 <- predict(fit,newdata=data.frame(x=x0),se.fit=T)
  lpolygon(c(x0,rev(x0)), c(y0$fit-2*y0$se.fit,rev(y0$fit+2*y0$se.fit)), col="gray75", border=F,...)
  lpoints(x, y, col="gray30", cex=0.5, pch=16)
  llines(x0, y0$fit, col="blue", lwd=3)
}
xyplot(spnbmd~age|gender, bmd, panel=myPanel, xlab="Age", ylab="Relative change in spinal BMD")