### Bootstrapping for inference about beta_1 in the toluca example:


toluca.data <- read.table(file = "http://people.stat.sc.edu/Hitchcock/tolucadata.txt", 
header=FALSE, col.names = c('lotsize', 'workhrs'))

# attaching the data frame:

attach(toluca.data)

# fitting the regression model:

toluca.reg <- lm(workhrs ~ lotsize)

# getting the summary regression output:

b1 <- coef(toluca.reg)[2]


# getting a 95% confidence interval for the true slope beta_1 (an indirect way):

alpha <- 0.05
b1 <- summary(toluca.reg)$coef[2,1]
s.b1 <- summary(toluca.reg)$coef[2,2]
error.df <- summary(toluca.reg)$df[2]
lower <- b1 - qt(1-alpha/2, df=error.df)*s.b1
upper <- b1 + qt(1-alpha/2, df=error.df)*s.b1
print(paste(100*(1-alpha), "percent CI for slope:", round(lower,4), round(upper,4)))

# Compare the classical CI above to the bootstrap CI results:

### Fixed X resampling:

B = 1000
b1.star.vec = rep(0,times=B)

for (i in 1:B) {

boot.samp.resids <- resid(toluca.reg)[sample(1:(nrow(toluca.data)), replace=T)]

y.star.vec <- fitted(toluca.reg) + boot.samp.resids

toluca.reg.boot <- lm(y.star.vec ~ lotsize)

b1.star.vec[i] <- coef(toluca.reg.boot)[2]
}

hist(b1.star.vec)

alpha=0.05

# Percentile-method 95% CI for beta_1:

L.p <- quantile(b1.star.vec, alpha/2)
U.p <- quantile(b1.star.vec, 1-alpha/2)
print(paste("Percentile-method 95% CI for beta_1:", round(L.p,4), round(U.p,4)) )

# Reflection-method 95% CI for beta_1:

d1 <- b1 - quantile(b1.star.vec, alpha/2)
d2 <- quantile(b1.star.vec, 1-alpha/2) - b1
L.r <- b1 - d2
U.r <- b1 + d1
print(paste("Reflection-method 95% CI for beta_1:", round(L.r,4), round(U.r,4)) )




#############################################################################

## Bootstrapping for inference about beta_1 in the blood pressure example:

bloodpressure.data <- 
read.table(file = "http://people.stat.sc.edu/Hitchcock/bloodpressuredata.txt", 
header=FALSE, col.names = c('age', 'dbp'))

# attaching the data frame:

attach(bloodpressure.data)

# Regressing the response, dbp, against the predictor, age 

bp.reg <- lm(dbp ~ age) 

# Regressing the absolute residuals against the predictor, age  
 
abs.res <- abs(resid(bp.reg))
bp.reg2 <- lm(abs.res ~ age)

# Defining the weights using the fitted values from this second regression:  

weight.vec <- 1/((fitted(bp.reg2))^2);


# Using the weights option in lm to get the WLS estimates:  
 
bp.wls.reg <- lm(dbp ~ age, weights = weight.vec)

b1 <- coef(bp.wls.reg)[2]

### Random X resampling:

B = 1000
b1.star.vec = rep(0,times=B)

for (i in 1:B) {

boot.samp <- bloodpressure.data[sample(1:(nrow(bloodpressure.data)), replace=T),]

bp.reg.boot <- lm(dbp ~ age, data=boot.samp)

abs.res.boot <- abs(resid(bp.reg.boot))
bp.reg2.boot <- lm(abs.res ~ age, data=boot.samp)
weight.vec.boot <- 1/((fitted(bp.reg2.boot))^2);

bp.wls.reg.boot <- lm(dbp ~ age, weights = weight.vec.boot, data=boot.samp)

b1.star.vec[i] <- coef(bp.wls.reg.boot)[2]
}

hist(b1.star.vec)

alpha=0.05

# Percentile-method 95% CI for beta_1:

L.p <- quantile(b1.star.vec, alpha/2)
U.p <- quantile(b1.star.vec, 1-alpha/2)
print(paste("Percentile-method 95% CI for beta_1:", round(L.p,4), round(U.p,4)) )

# Reflection-method 95% CI for beta_1:

d1 <- b1 - quantile(b1.star.vec, alpha/2)
d2 <- quantile(b1.star.vec, 1-alpha/2) - b1
L.r <- b1 - d2
U.r <- b1 + d1
print(paste("Reflection-method 95% CI for beta_1:", round(L.r,4), round(U.r,4)) )