# pdf("F:\\sta4210\\CH11TA01.pdf")

bpdata <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH11TA01.txt",
           header=F, col.names=c("age","dbp"))

attach(bpdata)
n <- length(dbp)

#### Fit Ordinary least squares, Regression of |e| vs age

bp.mod1 <- lm(dbp ~ age)
summary(bp.mod1)
e <- residuals(bp.mod1)
abs.e <- abs(e)
b1w <- coef(bp.mod1)[2]


abs_e.mod1 <- lm(abs.e ~ age)     # Regression of |e| on age
summary(abs_e.mod1)


s.hat <- predict(abs_e.mod1)    # Save the fitted SD's for W
w.hat <- 1/s.hat^2              # Obtain the estimated weights in W


### Fit weighted least squares using lm function and w.hat as weight

bp.modwls <- lm(dbp ~ age, weight=w.hat)
summary(bp.modwls)

###################### BEGIN BOOTSTRAP  (random pairs of (X,Y))

num.boot <- 1000
set.seed(34567)

b1w.boot <- rep(0,num.boot)

for (i in 1:num.boot)  {

boot.sample <- sample(1:n,size=n,replace=TRUE)

dbp.b <- dbp[boot.sample]
age.b <- age[boot.sample]

bp.mod1b <- lm(dbp.b ~ age.b)
abs.e <- abs(residuals(bp.mod1b))



abs_e.mod1b <- lm(abs.e ~ age.b)     # Regression of |e| on age


s.hat <- predict(abs_e.mod1b)    # Save the fitted SD's for W
w.hat <- 1/s.hat^2              # Obtain the estimated weights in W


### Fit weighted least squares using lm function and w.hat as weight

bp.modwlsb <- lm(dbp.b ~ age.b, weight=w.hat)
b1w.boot[i] <- coef(bp.modwlsb)[2]
}

hist(b1w.boot,breaks=seq(0.235,1.015,0.0325))

(b1w.boot_025 <- quantile(b1w.boot,.025))
(b1w.boot_975 <- quantile(b1w.boot,.975))

(b1.boot.sd <- sd(b1w.boot))

(d1 <- b1w-b1w.boot_025)
(d2 <- b1w.boot_975-b1w)

(beta1w.95CI <- c(b1w-d2,b1w+d1))



# dev.off()