1
$\begingroup$

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

CDF of Population

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

Sampling Distribution

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.

$\endgroup$
2
  • 2
    $\begingroup$ One key thing to know about bootstrapping is that its theoretical justification is purely asymptotic: that is, in the limit of large sample sizes. Samples of a few dozen (or sometimes, even tens of thousands) drawn from highly skewed parent distributions usually are not well characterized by asymptotic results. Intuitively, when the inclusion of a single observation in a bootstrap sample can radically change its characteristics, bootstrapping has to be unreliable -- and this condition characterizes extreme skewness. $\endgroup$ – whuber Aug 19 '20 at 19:45
  • 1
    $\begingroup$ Thanks for the reply and link. From a practical point of view, we just get one sample. Do you know of a rule of thumb for comparing the skewness of the sample with the sample size to determine if a bootstrap interval is likely to be "acceptable" (i.e, low probability of having "inadequate" coverage)? $\endgroup$ – Kendal Aug 19 '20 at 20:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.