This is the 3rd and final email of hints on hw 06. 

Note: I would like you to turn in whatever you have for hw 06 by this
Fri, Apr 28.  

For problem #3, which is ex 8.7 of HTF, perhaps a useful thing to do
is to show how to get the cubic spline basis functions and the design
matrix for spline regression out of R.  Here is some sample code...

library(splines) # to get the ns() function for lm()...
                 # ns() can be used to build design matrices
                 # and to evaluate spline functions at new data 
                 # points.

dat <- read.table("bone.data.txt",header=T);

x <- dat$age
y <- dat$spnbmd

mydf <- 20

fit <- lm(y ~ ns(x,df=mydf))

plot(x,y,xlab="Age",ylab="Relative Change in Spinal BMD")

N.new <- 100
new.x <- seq(min(x),max(x),length=N.new)

X <- cbind(1,ns(x,df=mydf)) # the matrix "H" in HTF, sec 8.2-8.3

new.X <- cbind(1,ns(new.x,df=mydf))  # rows are h(x)'s, sec 8.2-8.3

XTXI <- solve(t(X)%*%X)  # not very efficient but it'll do...

betas <- XTXI %*% t(X) %*% y

new.y <- new.X %*% betas # or new.X %*% fit$coef
#
# note: new.y <- predict(fit,data.frame(x=new.x))
# gives a somewhat different set of predictions
# and I'm not sure why...  I trust what I remember
# from linear models though so I will use
#        new.y <- new.X %*% betas
               
res.se <- summary(fit)$sigma 
new.se <- sqrt(diag(new.X %*% XTXI %*% t(new.X))) * res.se

z <- qnorm(.95)  # for 90% pointwise confidence bands...

lines(new.x,new.y)
lines(new.x,new.y + z*new.se,lty=2)
lines(new.x,new.y - z*new.se,lty=2)

# these confidence bands are pretty tight, but I don't think they're
# incorrectly computed---at least they're not any worse than the 
# confidence bands you get by fitting a simple linear regression 
# model lm(y ~ x).  

You can either build your own crossvalidation and boostrap functions
for parts (a) and (c), or use the "crossval" and "bootstrap" functions
in "library(bootstrap)"...  For the Bayesian part, part (b), you might
try a lot of different values of tau; this is analogous to the
"shrinkage" plots we made in the Bayesian part of the course before
spring break.

Hope this helps,

-BJ