source("http://www.uvm.edu/~rsingle/stat5210/data/scripts-5210.R")
 dat = otherdata("fev2.dat")
 dat$fev = 10*dat$fev
 head(dat)
 set.seed(17)

#Full model
 mod = lm(fev~height+age, data=dat)
 summary(mod)
 vcov(mod)

 L.hat.full = -1*mod$coef[2] + mod$coef[3]
 L.hat.full
 
#Compute a bootstrap estimate of L = -1*Beta_1 + Beta_2 and SE(L)
 
 boot.fn = function(data,index)
 { 
   tmp = lm(fev~height+age, data, subset=index)
   L.hat = -1*tmp$coef[2] + tmp$coef[3] 
   return( unname(L.hat) )
 }
 boot.fn(dat,  1:100) #evaluate L.hat for one sub-sample
 boot.fn(dat,131:215) #evaluate L.hat for another sub-sample

#install.packages(boot)
 library(boot)

#compute L.hat for 1000 bootstrap samples (w/ replacement) of same size as dat 
 temp = boot(dat, boot.fn, 1000)
 temp
 L.hat.full

#--------------------------------
#Attributes of a boot() object

 temp = boot(dat, boot.fn, R=1000)
 print(temp)
 attributes(temp)

 temp$t0                #mean of original data
 hist(temp$t)           #histogram of Bootstrap estimates           
 mean(temp$t)           #Bootstrap estimate (mean of Bootstrap samples)
 mean(temp$t) - temp$t0 #Bias
 text(temp$t0,0,"X")    #place an "X" at the original mean 

 t0   = temp$t0 
 bias = mean(temp$t) - temp$t0
#-------------------------------- 
#Bootstrap Confidence Intervals 
 boot.ci(temp)
 
#CI: Percentile 
 L95 = quantile(sort(temp$t),.025)
 U95 = quantile(sort(temp$t),.975)
 c(L95,U95) 
  
#CI: Normal (Wald)
 c(t0-bias - 1.96*0.1632286, t0-bias + 1.96*0.1632286)

#CI: Basic (Pivotal or Empirical) 
#    Better than Percentile for dists with strange tails 
 c(2*t0 - U95, 2*t0 - L95) 
  
#BCa (Bias-Corrected accelerated)