plot.loess <- function(x,...){ind <- order(x$x);plot(x$x[ind],x$fitted[ind],...)}
lines.loess <- function(x,...){ind <- order(x$x);lines(x$x[ind],x$fitted[ind],...)}

bmd <- read.delim("../../data/bmd.txt")
m <- bmd[bmd$gender=="male",]
f <- bmd[bmd$gender=="female",]

## Slide 2
x <- c(1:40,60:100)
xx <- seq(1,100,len=101)
y <- x+rnorm(length(x))
plot(x,y,pch=19,col="gray")
fit <- ksmooth(x,y,kernel="normal",bandwidth=30)
lines(fit$x,fit$y,col="red",lwd=3)

## Slide 6
par(mfrow=c(1,2))
plot(x,y,pch=19,col="gray")
fit <- ksmooth(x,y,kernel="normal",bandwidth=30)
lines(fit$x,fit$y,col="red",lwd=2)
abline(v=3,lty=2)
k <- function(x){(abs(x) < 1)*3/4*(1-x^2)}
lines(xx,y[3]+30*k((3-xx)/30),type="l",col="red",lty=2)

B <- cbind(1,x)
W <- diag(k((3-x)/30))
l <- as.numeric(B[3,]%*%solve(crossprod(B,W)%*%B)%*%crossprod(B,W))
L <- NULL
for (i in 1:length(y))
  {
    W <- diag(k((x[i]-x)/30))
    L <- rbind(L,as.numeric(B[i,]%*%solve(crossprod(B,W)%*%B)%*%crossprod(B,W)))
  }
L <- B%*%solve(crossprod(B,W)%*%B)%*%crossprod(B,W)
lines(x,y[3]+350*l,type="l",col="blue",lty=2)
fit <- loess(y~x,span=.3)
lines(fit,col="blue",lwd=2)
legend("top",lty=1,col=c("red","blue"),legend=c("N-W","Loess"))

plot(x,y,pch=19,col="gray")
fit <- ksmooth(x,y,kernel="normal",bandwidth=30)
lines(fit$x,fit$y,col="red",lwd=2)
abline(v=48,lty=2)
k <- function(x){(abs(x) < 1)*3/4*(1-x^2)}
lines(xx,y[40]+30*k((48-xx)/30),type="l",col="red",lty=2)
legend("top",lty=1,col=c("red","blue"),legend=c("N-W","Loess"))

## Slide 8
x <- seq(0,2*pi,len=101)
xx <- seq(min(x),max(x),len=101)
y <- sin(x) + rnorm(length(x),sd=.2)
plot(x,y,pch=19,col="gray")
fit <- loess(y~x,degree=1,span=.5)
lines(xx,predict(fit,newdata=xx),col="blue",lwd=3)
fit <- loess(y~x,degree=2,span=.8)
lines(xx,predict(fit,newdata=xx),col="red",lwd=3)
legend("topright",legend=c("Linear","Quadratic"),lty=1,col=c("blue","red"),lwd=3)

## Slide 10
span <- seq(.1,1,len=31)
degree <- 0:2
n <- nrow(m)
GCV <- p <- matrix(NA,nrow=length(span),ncol=3)
for (i in 1:length(span))
  {
    for (j in 1:length(degree))
      {
        fit <- loess(spnbmd~age,m,span=span[i],degree=degree[j])
        GCV[i,j] <- mean(((fit$y-fit$fitted)/(1-fit$trace.hat/n))^2)
        p[i,j] <- fit$trace.hat
      }
  }
matplot(span,GCV,type="l",lwd=3,lty=1,col=c("black","blue","red"))
legend("topleft",legend=c("degree=0","degree=1","degree=2"),lty=1,col=c("black","blue","red"))

matplot(span,p,type="l",lty=1,xlab="Span",ylab="Effective degrees of freedom",col=c("black","blue","red"),lwd=2)
legend("topright",legend=c("degree=0","degree=1","degree=2"),lty=1,col=c("black","blue","red"))
optimal <- apply(GCV,2,which.min)
points(span[optimal],c(p[optimal[1],1],p[optimal[2],2],p[optimal[3],3]),pch=19,col=c("black","blue","red"))

## Slide 13
fit <- loess(spnbmd~age,m,span=span[optimal[1]],degree=0)
plot(fit,type="l",ylim=range(m$spnbmd),xlab="Age",ylab="Relative change in spinal BMD",main="",lwd=2)
fit <- loess(spnbmd~age,m,span=span[optimal[2]],degree=1)
lines(fit,col="blue",lwd=2)
fit <- loess(spnbmd~age,m,span=span[optimal[3]],degree=2)
lines(fit,col="red",lwd=2)
points(m$age,m$spnbmd,pch=19,cex=0.7,col="gray")
abline(h=0,lty=2)
legend("topright",legend=c("degree=0","degree=1","degree=2"),lty=1,col=c("black","blue","red"))

## Slide 14
span <- seq(.1,1,len=31)
nm <- nrow(m)
nf <- nrow(f)
GCV <- p <- matrix(NA,nrow=length(span),ncol=2)
for (i in 1:length(span))
  {
    fit <- loess(spnbmd~age,m,span=span[i])
    GCV[i,1] <- sum(((fit$y-fit$fitted)/(1-fit$trace.hat/nm))^2)
    p[i,1] <- fit$trace.hat
    fit <- loess(spnbmd~age,f,span=span[i])
    GCV[i,2] <- sum(((fit$y-fit$fitted)/(1-fit$trace.hat/nf))^2)
    p[i,2] <- fit$trace.hat
  }
optimal <- apply(GCV,2,which.min)
p[optimal[1],1]
p[optimal[2],2]
plot(m$age,m$spnbmd,pch=19,col=rgb(.7,.7,1),ylim=range(bmd$spnbmd),xlab="Age",ylab="Relative change in spinal BMD")
points(f$age,f$spnbmd,pch=19,col=rgb(1,.7,.7))
fit <- loess(spnbmd~age,m,span=span[optimal[1]],degree=2)
lines(fit,col="blue",lwd=3)
fit <- loess(spnbmd~age,f,span=span[optimal[2]],degree=2)
lines(fit,col="red",lwd=3)
legend("topright",legend=c("Males","Females"),lty=1,lwd=3,col=c("blue","red"))
abline(h=0,lty=2)