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("http://web.as.uky.edu/statistics/users/pbreheny/621/data/bmd.txt")
m <- subset(bmd, gender=="male")
f <- subset(bmd, gender=="female")

## Slide 9
bw <- seq(1,6,len=21)
R <- matrix(NA,nrow=length(bw),ncol=nrow(m))
for (i in 1:length(bw)) {
  for (j in 1:nrow(m)) {
    yhat <- ksmooth(m$age[-j], m$spnbmd[-j], kernel="normal", bandwidth=bw[i], x.points=m$age[j])$y
    R[i,j] <- m$spnbmd[j]-yhat
  }
}
CV <- apply(R, 1, crossprod)/nrow(m)
plot(bw, 1000*CV, type="l", xlab="Bandwidth", ylab="CV", lwd=3, las=1)
bw[which.min(CV)]

## Nadaraya-Watson smooths
fit <- with(m, ksmooth(age,spnbmd,kernel="normal",bandwidth=1))
plot(m$age, m$spnbmd, pch=19, cex=0.7, col="gray", ylim=range(m$spnbmd), xlab="Age",ylab="Relative change in spinal BMD", main="h = 1")
lines(fit$x, fit$y, type="l", lwd=3, col="red")
abline(h=0,lty=2)
fit <- with(m, ksmooth(age,spnbmd,kernel="normal",bandwidth=2.75))
plot(m$age, m$spnbmd, pch=19, cex=0.7, col="gray", ylim=range(m$spnbmd), xlab="Age",ylab="Relative change in spinal BMD", main="h = 2.75")
lines(fit$x, fit$y, type="l", lwd=3, col="red")
abline(h=0,lty=2)
fit <- with(m, ksmooth(age,spnbmd,kernel="normal",bandwidth=10))
plot(m$age, m$spnbmd, pch=19, cex=0.7, col="gray", ylim=range(m$spnbmd), xlab="Age",ylab="Relative change in spinal BMD", main="h = 10")
lines(fit$x, fit$y, type="l", lwd=3, col="red")
abline(h=0,lty=2)

## Bias illustration
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)

## Kernel carpentry
col <- c("#FF4E37FF","#008DFFFF")
plot(x, y, pch=19, col="gray", cex=0.7)
fit <- ksmooth(x,y,kernel="normal",bandwidth=30)
lines(fit$x, fit$y, col=col[1], lwd=3)
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=col[1],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))
lines(x,y[3]+350*l,type="l",col=col[2],lty=2)
fit <- loess(y~x,span=.3)
lines(fit, col=col[2], lwd=3)

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

B <- cbind(1,x)
W <- diag(k((42-x)/30))
D <- cbind(0,42-x)
l <- as.numeric(B[42,]%*%solve(crossprod(B,W)%*%B)%*%crossprod(B,W))
ll <- predict(splines:::interpSpline(x,y[40]+350*l))
lines(ll$x, ll$y, col=col[2], lty=2)
fit <- loess(y~x,span=.3)
lines(fit,col=col[2], lwd=3)

## Slide 25
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=col[1], lwd=3)
fit <- loess(y~x,degree=2,span=.8)
lines(xx,predict(fit,newdata=xx),col=col[2], lwd=3)