STA 506 Homework 7: Bootstrap Confidence Intervals

1 Introduction

For this assignment, we will be building different types of bootstrapped confidence intervals to estimate the mean of a lognormal distribution. Specifically, we will use four different types: normal, percentile, bias-corrected and accelerated (BCa), and studentized (asymptotic) at confidence level of 95%. For all of them, we will first build them manually and then build them using the boot() and boot.ci() functions from the boot package.

To compute these confidence intervals, we will need data. The data we will be using consists of soil lead (Pb) concentrations (mg/kg) from 55 urban garden sites. In order to estimate the mean of a lognormal, we will assume this data follows a lognormal distribution.

lead <- c(0.85, 1.23, 0.92, 3.45, 2.11, 1.56, 4.89, 2.34, 1.78, 6.72, 0.95, 1.34, 8.91, 2.67, 1.89, 5.43, 1.12, 3.78, 2.45, 7.65, 1.05, 1.45, 12.34, 2.89, 2.01, 4.56, 1.23, 4.32, 2.67, 9.87, 0.99, 1.56, 15.23, 3.12, 2.34, 3.89, 1.34, 5.67, 2.89, 11.45, 1.12, 1.67, 18.90, 3.45, 2.56, 3.45, 1.45, 6.78, 3.12, 14.56, 1.23, 1.78, 22.34, 3.78, 2.78)

2 Part A

For Part A, we will be constructing a normal (or standard) bootstrapped confidence interval. For this interval, we will assume that our bootstrap sampling distribution will follow a normal distribution.

2.1 Manual Construction of the Confidence Interval

First, we will manually build the confidence interval.

## Manually Building a 95% Normal Bootstrap Confidence Interval

set.seed(1) # controls randomness of sampling

# Manual Bootstrap
B <- 5000
n <- length(lead)
boot_means <- numeric(B) # Empty Vector to Store Bootstrap Sample Means

# Loop to Create Samples and Calculate Sample Means
for(b in 1:B) {
  boot_sample <- sample(lead, n, replace=TRUE)  # Generating bootstrap sample
  boot_means[b] <- mean(boot_sample) # Calculate Bootstrap Sample Means
}


mean.samp <- mean(lead) # mean from the original sample
boot.se <- sd(boot_means) # standard error of bootstrap 
z <- qnorm(0.975) # critical value for confidence interval
moe = z*boot.se # margin of error for confidence interval

# The Confidence Interval 
ci_normal_manual <- c(LCL = mean.samp - moe, UCL = mean.samp + moe)
ci_normal_manual

     LCL      UCL 
3.173717 5.621919 

Our resuling confidence interval is (3.1737, 5.6219).

2.2 Constructing the Confidence Interval with the boot Package

Next, we will use the boot() and boot.ci() functions from the boot package to build the confidence interval.

## Calculating Bootstrap Confidence Interval Via Boot Package

set.seed(1) # still controls randomness of sampling

# creates statistic function 
mean_stat <- function(data, indices) {
  sample_data<- data[indices]
  return(mean(sample_data))
}

# Boot function generating bootstrap sampling distribution of the mean
boot.means <- boot(lead, mean_stat, R = 5000) # generating 5000 samples
# 95% Confidence Interval for Bootstrap Sample Mean
ci_boot_normal1 <- boot.ci(boot.means, type = "norm")
ci_boot_normal1

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.means, type = "norm")

Intervals : 
Level      Normal        
95%   ( 3.165,  5.655 )  
Calculations and Intervals on Original Scale

2.3 Histogram of the Sampling Distribution

And finally, we will plot the actual sampling distribution itself.

# plot histogram of sampling distribution of means
hist(boot.means$t, prob = TRUE, col = "lightblue", main = "Bootstrap Sampling Distribution of Sample Means", xlab = "Sample Mean") 

# Add Gaussian Density Curve of Sampling Distribution
lines(density(boot.means$t, kernel = "gaussian"), col = "red", lwd = 2)

There does appear to be some positive skew. However the skew looks a bit exaggerated due to the graph not being centered. Other than that, nothing looks out of the ordinary. The density curve definitely takes shape of a bell curve thanks to the 5000 sample means we have.

2.4 TO RECAP

Our manually constructed normal bootstrap confidence interval is (3.1737, 5.6219), and the interval constructed via the boot package is (3.165, 5.655). I expected the intervals to differ slightly because the sample mean differs from the bootstrap mean.

3 Part B

For this part, we will construct a percentile bootstrap confidence interval. These are the simplest confidence intervals, as they are built by directly taking the corresponding quantiles of the empirical distribution of bootstrap sampling distribution. This method works well when the bootstrap distribution is approximately symmetric and the estimator has a low amount of bias.

3.1 Manually Building the Confidence Interval

First, We will manually build this confidence interval. We will use the same sampling distribution from earlier. Since the sampling distribution looks roughly symmetric and we are building a 95% confidence interval, we will calculate the 2.5% and 97.5% percentiles and use those values as the lower and upper limits of this confidence interval.

boot.mean.values <- boot.means$t
alpha <- 0.05 # 1 - alpha = level of confidence

# The COnfidence Interval
ci_percentile_manual <- quantile(boot.mean.values, probs=c(alpha/2, 1-alpha/2))
ci_percentile_manual

    2.5%    97.5% 
3.245577 5.729459 

Our resulting confidence interval is (3.2456, 5.7295).

3.2 Constructing the Interval via boot Package

Then, we will construct our confidence interval via the functions from the boot package.

## Generate 95% Bootstrap Percentile Confidence Interval of Sample Mean
set.seed(1) # still controls randomness of sampling

# creates statistic function 
mean_stat <- function(data, indices) {
  sample_data<- data[indices]
  return(mean(sample_data))
}

# Boot function generating bootstrap sampling distribution of the mean
boot.means <- boot(lead, mean_stat, R = 5000) # generating 5000 samples
# 95% Confidence Interval for Bootstrap Sample Mean
# Percentile CI
ci_boot_percent <- boot.ci(boot.means, type="perc")
ci_boot_percent

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.means, type = "perc")

Intervals : 
Level     Percentile     
95%   ( 3.243,  5.730 )  
Calculations and Intervals on Original Scale

The Resulting Confidence Interval is (3.2430, 5.7300).

3.3 TO RECAP

The manually constructed percentile confidence interval is (3.2456, 5.7295), and the percentile confidence interval constructed via the boot package is (3.2530, 5.7300).

4 Part C

For this part, we will build the Bias-Corrected and Accelerated (BCa) Bootstrap Confidence Interval. This method is essentially an enhanced version of the percentile method: it provides adjusted quantiles by adjusting for bias and skewness in the sampling distribution. This occurs with two adjustment parameters: the bias correction factor, \(z_0\), which measures median bias by comparing the bootstrap samples to the original sample, and the acceleration factor, \(a\), which measures how the standard error changes.

To construct this interval, we need to rely on a specific type of resampling method called a jackknife. This is a resampling method in which you take n subsamples of your original sample, leaving one observation out of each of them and then calculating any statistic(s) of interest. Then, we get \(z_0\) by getting the proportion of observations in our bootstrap sample that is smaller than our statistic, and then we get \(a\) by seeing how far the confidence interval moves based on the bias adjustment.

## Manually Building A 95% BCa Confidence Interval

# Creates function
ci_bca <- function(data, B, alpha, original_estimate){ 
  n = length(data)
  boot_means <- numeric(B)
  
# Bootstrapping and Finding the Likelihood Estimators for Mean and Variance
for (i in 1:B) {
     boot_sample = sample(data, size = n, replace= TRUE)
     mean.lognorm[i] = sum(log(boot_sample))/n
     var.lognorm[i] = sum((log(boot_sample) - mean.lognorm[i])^2)/n
     bootstrap_means[i] = exp(mean.lognorm[i] +
                                      (var.lognorm[i]/2))
     
     # bias correction factor below
     ## can use normal if statement syntax, just incorporate the [i] to indicate that the statement is evaluating element by element
         if (bootstrap_ex_estimates[i] < original_estimate ){

         }
     } # END BOOTSTRAP LOOP
     
      #### START CORRECTION CALCULATION  
        prop_less = prop_less/b
        z0 = qnorm(prop_less)
        # now have bias correction
      #### END BIAS CORRECTION CALCULATION
  
        
      #### START ACCELERATION CALCULATION
          # good practice to use different "counting variable than was earlier used in function, I'll use j instead of i
     jack_vals = numeric(n)
     
     for (j in 1:n){ # START JACK LOOP
       jack_sample = data[-j]
       
       mu_jack = sum(log(jack_sample)) / (n-1)
        # the jackknife estimate is computed by leaving OUT the ith observation
       theta_squared_jack = sum((log(jack_sample)-mu_jack)^2) / (n-1)
       jack_vals[j] = exp(mu_jack + (theta_squared_jack/2))
     } # END JACK LOOP
     
       theta_jack_mean = mean(jack_vals)
             # we want the difference between the total average and each jackknife estimate that left one out
       numerator = sum((theta_jack_mean - jack_vals)^3)
       denominator = sum((theta_jack_mean - jack_vals)^2)
       
       a_hat = numerator/ (6 * denominator^(3/2))
     #### END ACCELERATION CALCULATION
       
       
    ### FINAL OUTPUT/INTERVAL
       z_score1 = qnorm(alpha/2)
       z_score2 = qnorm(1-alpha/2)
       alpha1 = pnorm(z0 + (z0 + z_score1)/(1-a_hat*(z0 + z_score1)))
       alpha2 = pnorm(z0 + (z0 + z_score2)/(1-a_hat*(z0 + z_score2)))
       
       lower = quantile(bootstrap_ex_estimates, alpha1)
       upper = quantile(bootstrap_ex_estimates, alpha2)
       results = rbind(lower,upper)
       rownames(results) = c("Lower Estimate", "Higher Estimate")
       colnames(results) = "BCA"
       
       return(results)
            
}

4.1 Using the boot Package

set.seed(1)

mean_stat <- function(data, indices) {
  sample_data <- data[indices]
  return(mean(data[indices]))
}

boot.means<- boot(lead, mean_stat, R=5000)

ci_bc_acc <- boot.ci(boot.means, type="bca")
ci_bc_acc 

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.means, type = "bca")

Intervals : 
Level       BCa          
95%   ( 3.391,  5.991 )  
Calculations and Intervals on Original Scale

The resulting confidence interval is (3.391, 5.991).

4.2 TO RECAP

The manually constructed Bias-Corrected and Accelerated interval is (,) and the confidence interval constructed via the boot package is (3.391, 5.991)

5 Part D

For Part D, we will be constructing a studentized (asymptotic) bootstrap confidence interval. This method accounts for variability across bootstrap samples.

## Using the boot Package

set.seed(1)

mean_stat <- function(data, indices) {
  sample_data <- data[indices]
  return(mean(data[indices]))

}
  
#boot_result <- boot(lead, mean_stat, R=5000)
#boot.ci(boot_result, type="stud")