> n <- length(dbp)
> 
> #### Fit Ordinary least squares, Regression of |e| vs age
> 
> bp.mod1 <- lm(dbp ~ age)
> summary(bp.mod1)

Call:
lm(formula = dbp ~ age)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.4786  -5.7877  -0.0784   5.6117  19.7813 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 56.15693    3.99367  14.061  < 2e-16 ***
age          0.58003    0.09695   5.983 2.05e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 8.146 on 52 degrees of freedom
Multiple R-squared: 0.4077,     Adjusted R-squared: 0.3963 
F-statistic: 35.79 on 1 and 52 DF,  p-value: 2.05e-07 

> 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)

Call:
lm(formula = abs.e ~ age)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.7639 -2.7882 -0.1587  3.0757 10.0350 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.54948    2.18692  -0.709  0.48179    
age          0.19817    0.05309   3.733  0.00047 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.461 on 52 degrees of freedom
Multiple R-squared: 0.2113,     Adjusted R-squared: 0.1962 
F-statistic: 13.93 on 1 and 52 DF,  p-value: 0.0004705 

> 
> 
> 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)

Call:
lm(formula = dbp ~ age, weights = w.hat)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.0230 -0.9939 -0.0327  0.9250  2.2008 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 55.56577    2.52092  22.042  < 2e-16 ***
age          0.59634    0.07924   7.526 7.19e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.213 on 52 degrees of freedom
Multiple R-squared: 0.5214,     Adjusted R-squared: 0.5122 
F-statistic: 56.64 on 1 and 52 DF,  p-value: 7.187e-10 

> 
> ###################### 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))
Error in hist.default(b1w.boot, breaks = seq(0.235, 1.015, 0.0325)) : 
  some 'x' not counted; maybe 'breaks' do not span range of 'x'
> 
> (b1w.boot_025 <- quantile(b1w.boot,.025))
     2.5% 
0.4277667 
> (b1w.boot_975 <- quantile(b1w.boot,.975))
    97.5% 
0.7537555 
> 
> (b1.boot.sd <- sd(b1w.boot))
[1] 0.08637445
> 
> (d1 <- b1w-b1w.boot_025)
      age 
0.1522641 
> (d2 <- b1w.boot_975-b1w)
    97.5% 
0.1737247 
> 
> (beta1w.95CI <- c(b1w-d2,b1w+d1))
      age       age 
0.4063061 0.7322949 
>