Poor Coverage for Bootstrap Confidence Interval on Mean

Question

Based on a highly skewed population, I was hoping to obtain good coverage with bootstrap (vs standard) confidence intervals. I’m not having much luck.

I specify my population in R as a Gamma distribution with shape = 0.3 and scale = 0.5. I think it’s easiest to see the skew of this population with a graph of its CDF:

plot(1:1000/1000, pgamma(1:1000/1000, shape = 0.3, scale = 0.5), 
     type = "l", yaxs = "i", xaxs = "i", xlim = c(0, 1), ylim = c(0, 1), col = "red", 
     xlab = "Value", ylab = "CDF", main = "Population CDF")

Note that the true mean for this gamma distribution is a function of its parameters: $0.3\times0.5 = 0.15$.

In each of 10,000 simulations, I sample 24 values from this distribution. In each simulation, I calculate the sample mean. I also calculate a standard confidence interval (based on assumed normality) and three bootstrap confidence intervals (using the boot package in R). For each simulation, I determine whether the true mean is within the confidence interval for each method. Note that these simulations took my computer about 25 minutes to run.

library(boot)
set.seed(1776)
(trueMean <- 0.3*0.5)

mean.boot <- function(x, i) mean(x[i]) #To be used with boot package for bootstrapping 

sampSize <- 24
numSims <- 10000
bootstrapIterations <- 10000
resultsDF <- data.frame(sampMean = numeric(numSims), withinStandard = logical(numSims),
                        withinNorm = logical(numSims), withinPerc = logical(numSims), 
                        withinBca = logical(numSims))

for (i in 1:numSims) {
  temp1 <- rgamma(sampSize, shape = 0.3, scale = 0.5)
  sampMean <- mean(temp1)
  sampVariance <- var(temp1)/sampSize
  standard_low <- sampMean - qt(0.95, sampSize - 1)*sqrt(sampVariance)
  standard_high <- sampMean + qt(0.95, sampSize - 1)*sqrt(sampVariance)
  
  boot1 <- boot(temp1, mean.boot, bootstrapIterations)
  bootCI1 <- boot.ci(boot1, conf = 0.9, type = c("norm", "perc", "bca"))
  normal_low <- bootCI1[[4]][2]
  normal_high <- bootCI1[[4]][3]
  perc_low <- bootCI1[[5]][4]
  perc_high <- bootCI1[[5]][5]
  bca_low <- bootCI1[[6]][4]
  bca_high <- bootCI1[[6]][5]
  
  resultsDF$sampMean[i] <- sampMean
  resultsDF$withinStandard[i] <- standard_low <= trueMean & standard_high >= trueMean
  resultsDF$withinNorm[i] <- normal_low <= trueMean & normal_high >= trueMean
  resultsDF$withinPerc[i] <- perc_low <= trueMean & perc_high >= trueMean
  resultsDF$withinBca[i] <- bca_low <= trueMean & bca_high >= trueMean
  if (i %% 100 == 0) print(paste0("Simulation: ", i))
}

Based on the simulations, the following histogram approximates the sampling distribution of the mean.

hist(resultsDF$sampMean, breaks = 39, xlim = c(0, 0.4), main = "Approximation of Sampling Distribution",
     xlab = "Sample Mean")

The histogram’s shape is far from normal, so it’s not surprising that the coverage for the standard confidence interval is poor. Unfortunately, the bootstrap intervals also have poor coverage.

mean(resultsDF$withinStandard)
mean(resultsDF$withinNorm)
mean(resultsDF$withinPerc)
mean(resultsDF$withinBca)

Specifically, I’m estimating the following coverages for my 90% confidence intervals:

Standard: 0.8314, Bootstrap Normal: 0.8128, Bootstrap Percentile: 0.8219, Bootstrap BCA: 0.8452.

I’m wondering if I’m doing something wrong, or if my data are just too skewed for a sample size of 24? Was there something I should have known about bootstrapping that would have made it obvious that my results would be poor? I figured that my results, even if imperfect, would be notably better than the standard results, which actually require a normal sampling distribution. Also, any guidance on what I can do to improve my results would be greatly appreciated.