Bootstrap vs Standard Error Confidence Intervals

[This article was first published on R on Jason Bryer, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

During a recent class a student asked whether bootstrap confidence intervals were more robust than confidence intervals estimated using the standard error (i.e. \(SE = \frac{s}{\sqrt{n}}\)). In order to answer this question I wrote a function to simulate taking a bunch of random samples from a population, calculate the confidence interval for that sample using the standard error approach (the t distribution is used by default, see the cv parameter. To use the normal distribution, for example, set cv = 1.96.), and then also calculating a confidence interval using the boostrap.

library(dplyr)
library(ggplot2)

#' Simulate random samples to estimate confidence intervals and bootstrap
#' estimates.
#'
#' @param pop a numeric vector representing the population.
#' @param n sample size for each random sample from the population.
#' @param n_samples the number of random samples.
#' @param n_boot number of bootstrap samples to take for each sample.
#' @param seed a seed to use for the random process.
#' @param cv critical value to use for calculating confidence intervals.
#' @return a data.frame with the sample and bootstrap mean and confidence
#'        intervals along with a logical variable indicating whether a Type I
#'        error would have occurred with that sample.
bootstrap_clt_simulation <- function(
		pop,
		n = 30,
		n_samples = 500,
		n_boot = 500,
		cv = abs(qt(0.025, df = n - 1)),
		seed,
		verbose = interactive()
) {
	if(missing(seed)) {
		seed <- sample(100000)
	}
	results <- data.frame(
		seed = 1:n_samples,
		samp_mean = numeric(n_samples),
		samp_se = numeric(n_samples),
		samp_ci_low = numeric(n_samples),
		samp_ci_high = numeric(n_samples),
		samp_type1 = logical(n_samples),
		boot_mean = numeric(n_samples),
		boot_ci_low = numeric(n_samples),
		boot_ci_high = numeric(n_samples),
		boot_type1 = logical(n_samples)
	)
	if(verbose) {
		pb <- txtProgressBar(min = 0, max = n_samples, style = 3)
	}
	for(i in 1:n_samples) {
		if(verbose) {
			setTxtProgressBar(pb, i)
		}
		set.seed(seed + i)
		samp <- sample(pop, size = n)
		boot_samp <- numeric(n_boot)
		for(j in 1:n_boot) {
			boot_samp[j] <- sample(samp, size = length(samp), replace = TRUE) |>
				mean()
		}
		results[i,]$seed <- seed + i
		results[i,]$samp_mean <- mean(samp)
		results[i,]$samp_se <- sd(samp) / sqrt(length(samp))
		results[i,]$samp_ci_low <- mean(samp) - cv * results[i,]$samp_se
		results[i,]$samp_ci_high <- mean(samp) + cv * results[i,]$samp_se
		results[i,]$samp_type1 <- results[i,]$samp_ci_low > mean(pop) |
			mean(pop) > results[i,]$samp_ci_high
		results[i,]$boot_mean <- mean(boot_samp)
		results[i,]$boot_ci_low <- mean(boot_samp) - cv * sd(boot_samp)
		results[i,]$boot_ci_high <- mean(boot_samp) + cv * sd(boot_samp)
		results[i,]$boot_type1 <- results[i,]$boot_ci_low > mean(pop) |
			mean(pop) > results[i,]$boot_ci_high
	}
	if(verbose) {
		close(pb)
	}
	return(results)
}

Uniform distribution for the population

Let’s start with a uniform distribution for our population.

pop_unif <- runif(1e5, 0, 1)
ggplot(data.frame(x = pop_unif), aes(x = x)) + geom_density()

The mean of the population is 0.4999484. We can now simulate samples and their corresponding bootstrap estimates.

results_unif <- bootstrap_clt_simulation(pop = pop_unif, seed = 42, verbose = FALSE)

4% of our samples did not contain the population mean in the confidence interval (i.e. Type I error rate) compared to rmean(results_unif$boot_type1) * 100`% of the bootstrap estimates. The following table compares the Type I errors for each sample compared to the bootstrap estiamted from that sample.

tab <- table(results_unif$samp_type1, results_unif$boot_type1, useNA = 'ifany')
tab
##        
##         FALSE TRUE
##   FALSE   477    3
##   TRUE      0   20

In general committing a type I error is the same regardless of method, though there were 3 instances where the bootstrap would have led to a type I error rate where the standard error approach would not.

The following plots show the relationship between the estimated mean (left) and condifence interval width (right) for each sample and its corresponding bootstrap.

results_unif |>
	ggplot(aes(x = samp_mean, y = boot_mean)) +
	geom_vline(xintercept = mean(pop_unif), color = 'blue') +
	geom_hline(yintercept = mean(pop_unif), color = 'blue') +
	geom_abline() +
	geom_point() +
	ggtitle("Sample mean vs bootstrap mean")

results_unif |>
	dplyr::mutate(samp_ci_width = samp_ci_high - samp_ci_low,
				  boot_ci_width = boot_ci_high - boot_ci_low) |>
	ggplot(aes(x = samp_ci_width, y = boot_ci_width)) +
	geom_abline() +
	geom_point() +
	ggtitle('Sample vs boostrap confidence interval width')

Skewed distribution for the population

We will repeat the same analysis using a positively skewed distribution.

pop_skewed <- rnbinom(1e5, 3, .5)
ggplot(data.frame(x = pop_skewed), aes(x = x)) + geom_density(bw = 0.75)

The mean of the population for this distribution is 2.99792

results_skewed <- bootstrap_clt_simulation(pop = pop_skewed, seed = 42, verbose = FALSE)
mean(results_skewed$samp_type1) # Percent of samples with Type I error
## [1] 0.05
mean(results_skewed$boot_type1) # Percent of bootstrap estimates with Type I error
## [1] 0.052
# CLT vs Bootstrap Type I error rate
table(results_skewed$samp_type1, results_skewed$boot_type1, useNA = 'ifany')
##        
##         FALSE TRUE
##   FALSE   473    2
##   TRUE      1   24
results_skewed |>
	ggplot(aes(x = samp_mean, y = boot_mean)) +
	geom_vline(xintercept = mean(pop_skewed), color = 'blue') +
	geom_hline(yintercept = mean(pop_skewed), color = 'blue') +
	geom_abline() +
	geom_point() +
	ggtitle("Sample mean vs bootstrap mean")

results_skewed |>
	dplyr::mutate(samp_ci_width = samp_ci_high - samp_ci_low,
				  boot_ci_width = boot_ci_high - boot_ci_low) |>
	ggplot(aes(x = samp_ci_width, y = boot_ci_width)) +
	geom_abline() +
	geom_point() +
	ggtitle('Sample vs boostrap confidence interval width')

We can see the results are very similar to that of the uniform distirubtion. Exploring the one case where the bootstrap would have resulted in a Type I error where the standard error approach would not reveals that it is very close with the difference being less than 0.1.

results_differ <- results_skewed |>
	dplyr::filter(!samp_type1 & boot_type1)
results_differ
##   seed samp_mean   samp_se samp_ci_low samp_ci_high samp_type1 boot_mean
## 1  443  3.866667 0.4516466    2.942946     4.790388      FALSE  3.924733
## 2  474  3.933333 0.4816956    2.948155     4.918511      FALSE  3.956800
##   boot_ci_low boot_ci_high boot_type1
## 1    3.044802     4.804665       TRUE
## 2    3.018549     4.895051       TRUE
set.seed(results_differ[1,]$seed)
samp <- sample(pop_skewed, size = 30)
boot_samp <- numeric(500)
for(j in 1:500) {
	boot_samp[j] <- sample(samp, size = length(samp), replace = TRUE) |>
		mean()
}
cv = abs(qt(0.025, df = 30 - 1))
mean(pop_skewed)
## [1] 2.99792
ci <- c(mean(samp) - cv * sd(samp) / sqrt(30), mean(samp) + cv * sd(samp) / sqrt(30))
ci
## [1] 2.942946 4.790388
mean(pop_skewed) < ci[1] | mean(pop_skewed) > ci[2]
## [1] FALSE
ci_boot <- c(mean(boot_samp) - cv * sd(boot_samp), mean(boot_samp) + cv * sd(boot_samp))
ci_boot
## [1] 3.044802 4.804665
mean(pop_skewed) < ci_boot[1] | mean(pop_skewed) > ci_boot[2]
## [1] TRUE

Adding an outlier

Let’s consider a sample that forces the largest value from the population to be in the sample.

set.seed(2112)
samp_outlier <- c(sample(pop_skewed, size = 29), max(pop_skewed))
boot_samp <- numeric(500)
for(j in 1:500) {
	boot_samp[j] <- sample(samp, size = length(samp), replace = TRUE) |>
		mean()
}

ci <- c(mean(samp_outlier) - cv * sd(samp_outlier) / sqrt(30), mean(samp_outlier) + cv * sd(samp_outlier) / sqrt(30))
ci
## [1] 1.647006 4.952994
mean(pop_skewed) < ci[1] | mean(pop_skewed) > ci[2]
## [1] FALSE
ci_boot <- c(mean(boot_samp) - cv * sd(boot_samp), mean(boot_samp) + cv * sd(boot_samp))
ci_boot
## [1] 2.905153 4.781381
mean(pop_skewed) < ci_boot[1] | mean(pop_skewed) > ci_boot[2]
## [1] FALSE

In this example we do see that the presense of the outlier does have a bigger impact on the confidence interval with the bootstrap confidence interval being much smaller.

To leave a comment for the author, please follow the link and comment on their blog: R on Jason Bryer.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)