library(bootstrap) x = diabetes$age y = diabetes$logCpeptide n = length(y) plot (x,y) sm.0 = smooth.spline(x,y,all.k=T,spar=0) lines(predict(sm.0),col="green") sm.05 = smooth.spline(x,y,all.k=T,spar=0.5) lines(predict(sm.05),col="red") sm.1 = smooth.spline(x,y,all.k=T,spar=1) lines(predict(sm.1),col="blue") sm.2 = smooth.spline(x,y,all.k=T,spar=2) lines(predict(sm.2),col="black") #plug-in estimate of lambda spars = seq(0,2,by=0.01) lambs = loss= NULL for (s in spars){ sm.now = smooth.spline(x,y,all.kn=T,spar=s) lambs = c(lambs, sm.now$lambda) pred.now = predict(sm.now,x)$y loss = c(loss, mean((y-pred.now)^2)) } boot.l = boot.o = boot.b632 = numeric (length(spars)) boot.n632=NULL B=20 for (b in 1:B){ ist = sample (n,rep=T) xst = x[ist] yst = y[ist] not.ist = (1:n)[!((1:n)%in%ist)] loss.now = opt.now = b632.now= NULL for (spar in spars){ sm.now = smooth.spline(xst,yst,all.kn=T,spar=spar) pred.now = predict(sm.now,x)$y loss.now = c(loss.now, mean((y-pred.now)^2)) pred.st.now = predict(sm.now,xst)$y opt.now = c(opt.now, mean((y-pred.now)^2)-mean((yst-pred.st.now)^2)) pred.notst.now = predict(sm.now,x[not.ist])$y b632.now = c(b632.now,sum((y[not.ist]-pred.notst.now)^2)) } boot.l = boot.l + loss.now boot.o = boot.o + opt.now boot.b632 = boot.b632 + b632.now boot.n632 = c(boot.n632, length(not.ist)) } ll = lambs<0.5 par (mfrow=c(2,2)) plot (lambs[ll], boot.l[ll]/B,ylim=c(0.33,0.45),xlim=c(0,0.03),type="l") plot (lambs[ll], loss[ll]+boot.o[ll]/B,ylim=c(0.33,0.45),xlim=c(0,0.03),type="l") plot (lambs[ll], 0.368*loss[ll]+0.632*boot.b632[ll]/sum(boot.n632),xlim=c(0,0.03),ylim=c(0.33,0.45),type="l")