## Slide 6
f <- function(x){dnorm(x,0.5,0.18)}
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)))
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)
  }

## Slide 9
bmd <- read.delim("../../data/bmd.txt")
trellis.par.set(plot.symbol=list(pch=19))
xyplot(spnbmd~age|gender,bmd)

## Remaining slides require these lines
require(splines)
myPanel <- function(...)
  {
    panel.smooth(...)
    panel.xyplot(...,alpha=.5)
  }

## Slide 11
panel.smooth <- function(x,y,...)
  {
    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="black")
    lsegments(x0=knots[1],y0=min(y),x1=knots[1],y1=max(y),col="black",lty=2)
    llines(x2,c2,col="black")
    lsegments(x0=knots[2],y0=min(y),x1=knots[2],y1=max(y),col="black",lty=2)
    llines(x3,c3,col="black")
  }
trellis.par.set(plot.symbol=list(pch=19),plot.line=list(lwd=2))
xyplot(spnbmd~age|gender,bmd,panel=myPanel)

## Slide 12
panel.smooth <- function(x,y,...)
  {
    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="black")
    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="black")
    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="black")
  }
xyplot(spnbmd~age|gender,bmd,panel=myPanel)

## Slide 13
panel.smooth <- function(x,y,...)
  {
    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="black")
    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)

## Slide 18
panel.smooth <- function(x,y,...)
  {
    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="black")
    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)
  }
trellis.par.set(plot.symbol=list(pch=19),plot.line=list(lwd=2))
xyplot(spnbmd~age|gender,bmd,panel=myPanel)

## Slide 19
panel.smooth <- function(x,y,...)
  {
    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="black")
    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)

## Slide 23
panel.smooth <- function(x,y,...)
  {
    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="black")
    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)

## Slide 24
panel.smooth <- function(x,y,...)
  {
    fit <- lm(y~ns(x,df=5))
    fitp <- lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5))
    x0 <- seq(min(x),max(x),len=101)
    llines(x0,predict(fit,data.frame(x=x0)),col="black",lwd=3)
    llines(x0,predict(fitp,data.frame(x=x0)),col="red",lwd=3)
  }
xyplot(spnbmd~age|gender,bmd,panel=myPanel)

## Slide 25
plot(bmd$age[bmd$gender=="female"],bmd$spnbmd[bmd$gender=="female"],col=rgb(1,0,0,alpha=.5),pch=20,xlab="Age",ylab="Relative change in spinal BMD")
points(bmd$age[bmd$gender=="male"],bmd$spnbmd[bmd$gender=="male"],col=rgb(0,0,1,alpha=.5),pch=20)
age.f <- bmd$age[bmd$gender=="female"]
age.m <- bmd$age[bmd$gender=="male"]
fit.f <- lm(spnbmd[bmd$gender=="female"]~ns(age.f,df=5),bmd)
fit.m <- lm(spnbmd[bmd$gender=="male"]~ns(age.m,df=5),bmd)
age <- seq(min(bmd$age),max(bmd$age),len=101)
lines(age,predict(fit.f,data.frame(age.f=age)),col="red",lwd=2)
lines(age,predict(fit.m,data.frame(age.m=age)),col="blue",lwd=2)
legend("topright",legend=c("female","male"),col=c("red","blue"),lty=1)

## Slide 30
panel.smooth <- 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="gray",border=F,...)
    llines(x0,y0$fit,col="black")
  }
xyplot(spnbmd~age|gender,bmd,panel=myPanel)

## Slide 31
plot(bmd$age[bmd$gender=="female"],bmd$spnbmd[bmd$gender=="female"],col=rgb(1,0,0,alpha=.5),pch=20,xlab="Age",ylab="Relative change in spinal BMD")
points(bmd$age[bmd$gender=="male"],bmd$spnbmd[bmd$gender=="male"],col=rgb(0,0,1,alpha=.5),pch=20)
age.f <- bmd$age[bmd$gender=="female"]
age.m <- bmd$age[bmd$gender=="male"]
fit.f <- lm(spnbmd[bmd$gender=="female"]~ns(age.f,df=5),bmd)
fit.m <- lm(spnbmd[bmd$gender=="male"]~ns(age.m,df=5),bmd)
age <- seq(min(bmd$age),max(bmd$age),len=101)
y.f <- predict(fit.f,newdata=data.frame(age.f=age),se.fit=TRUE)
y.m <- predict(fit.m,newdata=data.frame(age.m=age),se.fit=TRUE)
polygon(c(age,rev(age)),c(y.f$fit-2*y.f$se.fit,rev(y.f$fit+2*y.f$se.fit)),col=rgb(1,0,0,alpha=.5),border=F)
polygon(c(age,rev(age)),c(y.m$fit-2*y.m$se.fit,rev(y.m$fit+2*y.m$se.fit)),col=rgb(0,0,1,alpha=.5),border=F)
lines(age,predict(fit.f,data.frame(age.f=age)),col="red",lwd=3)
lines(age,predict(fit.m,data.frame(age.m=age)),col="blue",lwd=3)
legend("topright",legend=c("female","male"),col=c("red","blue"),lty=1,lwd=3)