7  Confidence Intervals

---

The sampling distribution describes how a statistic varies across samples; the confidence interval turns knowledge about that sampling distribution into a statement about the unknown parameter.

This chapter asks one main question: Which population values are reasonable, and which are unreasonable?

• Bootstrap percentile interval: this takes the middle quantiles of the bootstrap distribution of a statistic as the reasonable values.
• Normal approximation: this gives the interval \(\hat{m} \pm E\), where the margin of error \(E\) comes from the estimated standard error and a mathematical formula for the Normal distribution.

ImportantKey Definition

A confidence interval is a range of values, computed from the sample, intended to contain the unknown population parameter.

The coverage level \(1-\alpha\) is the long-run proportion of intervals (constructed from repeated samples by the same procedure) that contain the true parameter.

Along with the estimate \(\hat{m}\), we report a range and a coverage level as a “\(Z\%\) CI”.

Confidence intervals are also useful for reporting uncertainty alongside an estimate. The range shows how much sampling noise could plausibly move the estimate.

7.1 Resampling Intervals

If we have a sampling distribution for a statistic, the simplest way to read off a confidence interval is to take its middle quantiles.

The critical values of a sampling distribution are the lower and upper quantiles that bracket the middle \(Z\%\) of the distribution. Critical values are useful because they translate any coverage level into a pair of numbers that can be looked up or simulated, without committing to a parametric model for the sampling distribution. Often, we have only one sample, in which case we can construct a bootstrap or jackknife sampling distribution (as in Sampling & Resampling) and compute the critical values directly from it:

• The \(5^{\text{th}}\) and \(95^{\text{th}}\) percentiles are the critical values for the \(90\%\) confidence interval.
• The \(2.5^{\text{th}}\) and \(97.5^{\text{th}}\) percentiles are the critical values for the \(95\%\) confidence interval.
• The \(0.5^{\text{th}}\) and \(99.5^{\text{th}}\) percentiles are the critical values for the \(99\%\) confidence interval.

The histogram below shows the bootstrap distribution of the mean with the \(95\%\) interval marked.

## MAIN EXAMPLE
# Pretend these 3294 hourly wages are the whole population
# Source: US workers, with hourly wages in 1980 dollars
data(Wages1, package='Ecdat')
population_wages <- Wages1[, 'wage']

# Lab data: one random sample of 50 workers, as if it were all we observed
set.seed(1) # to be replicable
x_hat <- sample(population_wages, 50, replace=FALSE)
m_hat <- mean(x_hat)
m_hat
## [1] 5.708014
# The sample mean is about 5.71 dollars an hour.

# Bootstrap Distribution
boot_means <- rep(NA, 9999)
for(b in seq_along(boot_means)){
    x_boot <- sample(x_hat, replace=TRUE) # c.f. jackknife
    m_b_hat <- mean(x_boot)
    boot_means[b] <-m_b_hat
}

# 95% CI
boot_ci <- quantile(boot_means, probs=c(.025, .975), na.rm=TRUE)
boot_ci
##    2.5%   97.5% 
## 4.96511 6.46076
# The interval runs from about 4.97 to 6.46 dollars.
# We are 95% confident it contains the population mean wage.
# The pretend population mean is about 5.76, so this interval does contain it.

# Visualize
hist(boot_means, freq=FALSE,
    main=NA,
    border=NA, breaks=50,
    xlab='Bootstrap means of hourly wages (1980 dollars)')
title('Bootstrap CI', font.main=1)
abline(v=boot_ci, col=rgb(0, 0, 0, .8), lwd=2)

NoteMust Know

Using the bootstrap distribution above, here are the \(90\%\) and \(99\%\) confidence intervals.

# 90% CI: 5th and 95th percentiles
quantile(boot_means, probs=c(.05, .95), na.rm=TRUE)
##       5%      95% 
## 5.086016 6.335946

# 99% CI: 0.5th and 99.5th percentiles
quantile(boot_means, probs=c(.005, .995), na.rm=TRUE)
##     0.5%    99.5% 
## 4.758652 6.695011

The \(99\%\) interval is wider than the \(90\%\) interval: demanding higher coverage forces a longer interval.

NoteMust Know

Here is the same bootstrap \(95\%\) confidence interval for a second sample, x2_hat.

# Hourly wages of the 1569 women in Wages1 (1980 dollars)
# Wages1 is really a sample of US workers, so these women are a sample too.
x2_hat <- Wages1[Wages1[, 'sex']=='female', 'wage']
m2_hat <- mean(x2_hat)
m2_hat
## [1] 5.146924

# Bootstrap Distribution
boot_means2 <- rep(NA, 9999)
for(b in seq_along(boot_means2)){
    x_boot <- sample(x2_hat, replace=TRUE)
    m_b_hat <- mean(x_boot)
    boot_means2[b] <- m_b_hat
}

# 95% CI
boot_ci2 <- quantile(boot_means2, probs=c(.025, .975))
boot_ci2
##     2.5%    97.5% 
## 5.007022 5.289265
# The interval runs from about 5.01 to 5.29 dollars.
# We are 95% confident it contains the population mean wage for women.

# Visualize
hist(boot_means2, freq=FALSE,
    main=NA,
    border=NA, breaks=50,
    xlab='Bootstrap means of hourly wages (1980 dollars)')
title('Bootstrap CI', font.main=1)
abline(v=boot_ci2, col=rgb(0, 0, 0, .8), lwd=2)

TipTest Yourself

Repeat the bootstrap \(95\%\) confidence interval with the murder arrest rates in USArrests.

# These are one sample with no known population.
# So there is no true mean to compare with.
x3_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
m3_hat <- mean(x3_hat)
m3_hat
## [1] 7.788

# Bootstrap Distribution
boot_means3 <- rep(NA, 9999)
for(b in seq_along(boot_means3)){
    x_boot <- sample(x3_hat, replace=TRUE) # c.f. jackknife
    m_b_hat <- mean(x_boot)
    boot_means3[b] <- m_b_hat
}

# 95% CI
boot_ci3 <- quantile(boot_means3, probs=c(.025, .975), na.rm=TRUE)
boot_ci3
##   2.5%  97.5% 
## 6.6099 8.9720
# The interval runs from about 6.6 to 9.0 arrests per 100,000.

# Visualize
hist(boot_means3, freq=FALSE,
    main=NA,
    border=NA, breaks=50)
title('Bootstrap CI', font.main=1)
abline(v=boot_ci3, col=rgb(0, 0, 0, .8), lwd=2)

This approach can be easily extended to other statistics. For example, a \(Z\%\) confidence interval for the median is an interval where \(Z\%\) of those generated will contain the population median.

TipTest Yourself

The following code creates a bootstrap CI for the median. Edit it to create a jackknife CI instead. Then edit it to create a bootstrap CI for the upper quartile.

med_hat <- quantile(x_hat, probs=0.5) #median
med_hat
##      50% 
## 5.274198

# Bootstrap Distribution
boot_medians <- rep(NA, 9999)
for(b in seq_along(boot_medians)){
    x_boot <- sample(x_hat, replace=TRUE) # c.f. jackknife
    med_b_hat <- quantile(x_boot, probs=0.5)
    boot_medians[b] <-med_b_hat
}

# 95% CI
boot_ci <- quantile(boot_medians, probs=c(.025, .975), na.rm=TRUE)
boot_ci
##     2.5%    97.5% 
## 4.312633 6.736894

# Visualize
hist(boot_medians, freq=FALSE,
    main=NA,
    border=NA, breaks=50)
title('Bootstrap CI', font.main=1)
abline(v=boot_ci, col=rgb(0, 0, 0, .8), lwd=2)

7.2 Normal Approximation

When the sampling distribution is well approximated by a Normal curve, the confidence interval has a simple closed form built from the standard error. The formula lets us see that the CI width depends on

• the sample size
• the spread of the data
• the coverage level \(1-\alpha\)

ImportantKey Definition

The margin of error \(E\) is the half-width of a Normal-approximation confidence interval, so the interval is \(\hat{m} \pm E\). It equals the critical value times the standard error, \(E = q(\alpha/2) \cdot \text{SE}(M) \approx q(\alpha/2) \hat{s}/\sqrt{n}\).

The Normal-approximation interval is useful because it requires only an estimate \(\hat{m}\), a standard error \(\text{SE}(M)\), and a single Normal quantile (no resampling loop, no full sampling distribution). I.e., the classical formula uses the adjusted sample standard deviation \(\hat{s}_{\text{adj}}\); compute \(\text{SE}(M) \approx \hat{s}_{\text{adj}}/\sqrt{n}\). While \(\text{SE}(M)\) typically comes from the classical formula, it can also come from the jackknife (or even the bootstrap).

The \(95\%\) interval is the most common and uses critical value \(q=1.96\) (the \(97.5\)th percentile of the standard Normal); the \(90\%\) interval uses \(q=1.645\). In any case, a coverage level of \(1-\alpha\) means that if the same sampling procedure were repeated \(100\) times from the same population, approximately \(100(1-\alpha)\) percent of the intervals are expected to contain the true population mean. ¹

NoteMust Know

For example, suppose we sample \(n=36\) monthly electricity bills (in dollars) and find \(\hat{m}=142\) and \(\hat{s}_{\text{adj}}=30\). The classical standard error is \(\text{SE}(M) \approx \hat{s}_{\text{adj}}/\sqrt{n} = 30/\sqrt{36} = 5\). From theory, we also know that \(95\%\) confidence interval has upper and lower critical value of \[\begin{aligned} F^{-1}_{\text{Normal}}\left(0.975, \text{SE}(M) \right) &= +1.96 \cdot \text{SE}(M) \\ F^{-1}_{\text{Normal}}\left(0.025, \text{SE}(M) \right) &= -1.96 \cdot \text{SE}(M) \end{aligned}\] Together, we can calculate \[\begin{aligned} E &= q(\alpha/2) \cdot \text{SE}(M) = 1.96 \cdot 5 = 9.8 \\ [\hat{m} - E, \hat{m} + E] &= [142 - 9.8, 142 + 9.8] = [132.2,~ 151.8]. \end{aligned}\] We are \(95\%\) confident that the population mean electricity bill is between \(\$132.20\) and \(\$151.80\).

m_bill_hat <- 142
s_bill_adj_hat <- 30
n <- 36
SE <- s_bill_adj_hat / sqrt(n)
SE
## [1] 5

# 95% CI Manual Calculation
m_bill_hat + c(-1.96, 1.96) * SE
## [1] 132.2 151.8

## 95% CI Intuition
qnorm( c(0.025, 0.975), m_bill_hat, SE)
## [1] 132.2002 151.7998

The main advantages of the Normal approximation are that it

• can be computed formulaically as above
• works well for estimating extreme probabilities, where resampling methods tend to be worse.

The main disadvantage is that

• the sampling distribution might be far from normal.

The Normal approximation is least reliable when the sample is small or the data are strongly skewed: the sampling distribution is then asymmetric, and a symmetric \(\hat{m}\pm E\) interval cannot capture that. The percentile bootstrap interval reads quantiles directly off the resampling distribution, so it can be lopsided when the data call for it. The figure below compares the percentile and Normal intervals, which are not always similar.

## MAIN EXAMPLE
# Bootstrap Distribution with Percentile CI
hist(boot_means, breaks=25,
    main=NA, border=NA,
    freq=FALSE, ylim=c(0, 1.2),
    xlab=expression(hat(m)^group('(', b, ')')))
title('Percentile vs Normal 95% CIs', font.main=1)
boot_ci_percentile <- quantile(boot_means, probs=c(.025, .975))
abline(v=boot_ci_percentile, col=rgb(0, 0, 0, .8), lwd=2)

# Normal Approximation with IID Theory SEs
classic_se <- sd(x_hat)/sqrt(length(x_hat))
x <- seq(3, 9, by=0.01)
fx <- dnorm(x, m_hat, classic_se)
lines(x, fx, col=rgb(0, 0, 1, .8), lty=2)
ci_normal <- qnorm(c(.025, .975), m_hat, classic_se)
#ci_normal <- m_hat+c(-1.96, +1.96)*classic_se
abline(v=ci_normal, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Percentile', 'Normal'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

# The percentile and Normal intervals are quite similar for this sample.

NoteMust Know

Here is the same comparison for the second sample, x2_hat.

# x2_hat (hourly wages of the 1569 women), m2_hat, and boot_means2 defined above
c(mean=m2_hat, median=median(x2_hat))
##     mean   median 
## 5.146924 4.693262
# Hourly wages are strongly right-skewed: the median is well below the mean.

# Bootstrap Distribution with Percentile CI
hist(boot_means2, breaks=25,
    main=NA, border=NA,
    freq=FALSE, ylim=c(0, 5.5),
    xlab=expression(hat(m)^group('(', b, ')')))
title('Percentile vs Normal 95% CIs', font.main=1)
boot_ci2_percentile <- quantile(boot_means2, probs=c(.025, .975))
abline(v=boot_ci2_percentile, col=rgb(0, 0, 0, .8), lwd=2)

# Normal Approximation with IID Theory SEs
classic_se2 <- sd(x2_hat)/sqrt(length(x2_hat))
x <- seq(4.5, 5.8, by=0.005)
fx2 <- dnorm(x, m2_hat, classic_se2)
lines(x, fx2, col=rgb(0, 0, 1, .8), lty=2)
ci2_normal <- qnorm(c(.025, .975), m2_hat, classic_se2)
abline(v=ci2_normal, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Percentile', 'Normal'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

# Compare the two intervals
rbind(percentile=boot_ci2_percentile, normal=ci2_normal)
##                2.5%    97.5%
## percentile 5.007022 5.289265
## normal     5.004605 5.289242
# With 1569 women, the bootstrap distribution of the mean is close to Normal.
# So the two intervals nearly coincide, even though wages are skewed.

The Normal approximation needs the sampling distribution of the mean to be close to Normal, not the data themselves. Skewed data are mainly a problem when the sample is also small.

TipTest Yourself

Repeat the comparison with the murder arrest rates in x3_hat.

# x3_hat, m3_hat, and boot_means3 defined above
# Bootstrap Distribution with Percentile CI
hist(boot_means3, breaks=25,
    main=NA, border=NA,
    freq=FALSE, ylim=c(0, 0.7),
    xlab=expression(hat(m)^group('(', b, ')')))
title('Percentile vs Normal 95% CIs', font.main=1)
boot_ci3_percentile <- quantile(boot_means3, probs=c(.025, .975))
abline(v=boot_ci3_percentile, col=rgb(0, 0, 0, .8), lwd=2)

# Normal Approximation with IID Theory SEs
classic_se3 <- sd(x3_hat)/sqrt(length(x3_hat))
x <- seq(5, 10, by=0.01)
fx3 <- dnorm(x, m3_hat, classic_se3)
lines(x, fx3, col=rgb(0, 0, 1, .8), lty=2)
ci3_normal <- qnorm(c(.025, .975), m3_hat, classic_se3)
abline(v=ci3_normal, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Percentile', 'Normal'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

# The percentile and Normal intervals are quite similar for these data.

TipTest Yourself

Note that you can also plug in Bootstrap or Jackknife SEs to the Normal Approximation.

## Bootstrap Distribution
hist(boot_means, breaks=25,
    main=NA, border=NA,
    freq=FALSE, ylim=c(0, 1.2),
    xlab=expression(hat(m)^group('(', b, ')')))
abline(v=boot_ci_percentile, col=rgb(0, 0, 0, .8), lwd=2)
title('Percentile vs Normal 95% CIs (with Jackknife SE)', font.main=1)

# Jackknife SE Estimate
n <- length(x_hat)
jack_means <- rep(NA, n)
for(i in seq_along(jack_means)){
    x_jack <- x_hat[-i]
    m_i_hat <- mean(x_jack)
    jack_means[i] <- m_i_hat
}
se_jack <- sd(jack_means)*sqrt(n)

# Normal Approx with Jack SE
x <- seq(3, 9, by=0.01)
fx <- dnorm(x, m_hat, se_jack)
lines(x, fx, col=rgb(0, 0, 1, .8), lty=2)
ci_normal_alt <- qnorm(c(.025, .975), m_hat, se_jack)
abline(v=ci_normal_alt, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Percentile', 'Normal'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

Update the example to use se_boot <- sd(boot_means) in place of se_jack

7.3 Theoretical Properties

A \(Z\%\) confidence interval for the mean is an interval where \(Z\%\) of the intervals generated on different samples will contain the population mean \(\mathbb{E}[X_{i}]\). As in Sampling, pretend the full dataset that the lab sample came from is the whole population, so we know the population mean. We draw \(100\) samples of \(n=30\) observations, compute a \(90\%\) Normal confidence interval for each, and check which intervals contain the population mean.

## MAIN EXAMPLE
# population_wages, the pretend population, defined above
population_wage_mean <- mean(population_wages)

n <- 30
n_samples <- 100
sample_CIs <- matrix(nrow=n_samples, ncol=2)
for(i in seq(1, n_samples)){
    # Random sample of workers
    x_sim <- sample(population_wages, n, replace=FALSE)
    # Confidence Interval for each sample
    m_sim_hat <- mean(x_sim)
    SE <- sd(x_sim)/sqrt(n)
    CI <- qnorm( c(0.05, 0.95), m_sim_hat, SE)
    sample_CIs[i, ] <- CI
}

# Explicit calculation: whether the population mean is in each CI
covered <- population_wage_mean >= sample_CIs[, 1] & population_wage_mean <= sample_CIs[, 2]
sum(covered) # number of intervals that contain the population mean
## [1] 89
# The count is typically a bit below 90.
# Wages are right-skewed, so with n=30 the Normal approximation is not exact.
# The long-run coverage of this procedure is about 88%, not 90%.

# Visualize confidence intervals
plot.new()
plot.window(xlim = range(sample_CIs), ylim = c(0, n_samples))
for(i in seq(1, n_samples) ) {
  col_i <- if (covered[i]) rgb(0, 0, 0, 0.3) else rgb(1, 0, 0, 0.5)
  segments(
    sample_CIs[i, 1], i,
    sample_CIs[i, 2], i,
    col = col_i, lwd = 2)
}
abline(v = population_wage_mean, col = rgb(0, 0, 1, .8), lty = 2)
axis(1)
title('90% Coverage (Red = Missed)', font.main=1)

When the data are skewed and the sample is small, the Normal approximation is not exact, so the long-run coverage can fall below the chosen level.

NoteMust Know

Here is the same coverage check for uniform random data, where the true mean is \(0.5\).

n <- 30
n_samples <- 100
sample_CIs <- matrix(nrow=n_samples, ncol=2)
for(i in seq(1, n_samples)){
    # Random Sample
    x_sim <- runif(n)
    # Confidence Interval for each sample
    m_sim_hat <- mean(x_sim)
    SE <- sd(x_sim)/sqrt(n)
    CI <- qnorm( c(0.05, 0.95), m_sim_hat, SE)
    sample_CIs[i, ] <- CI
}


# Explicit calculation: whether the true mean is in each CI
true_mean <- 0.5 # theoretical result for uniform samples
covered <- true_mean >= sample_CIs[, 1] & true_mean <= sample_CIs[, 2]

# Visualize confidence intervals
plot.new()
plot.window(xlim = range(sample_CIs), ylim = c(0, n_samples))
for(i in seq(1, n_samples) ) {
  col_i <- if (covered[i]) rgb(0, 0, 0, 0.3) else rgb(1, 0, 0, 0.5)
  segments(
    sample_CIs[i, 1], i,
    sample_CIs[i, 2], i,
    col = col_i, lwd = 2)
}
abline(v = true_mean, col = rgb(0, 0, 1, .8), lty = 2)
axis(1)
title('90% Coverage (Red = Missed)', font.main=1)

We typically have only one sample and compute a single confidence interval with the data we have. The interval you computed either contains the true mean or it does not. ² So “\(Z\%\) confident” translates to “\(Z\%\) of intervals constructed on different samples will contain \(\mathbb{E}[X_{i}]\)”.

NoteMust Know

A confidence interval contains many values, and every one of them is a hypothesis the data cannot rule out. Suppose a \(95\%\) confidence interval for a population mean comes out as \([6.6, 9.0]\). The data are consistent with \(\mathbb{E}[X_{i}]=7\), and also with \(\mathbb{E}[X_{i}]=8\), and also with \(\mathbb{E}[X_{i}]=8.9\): none of these competing values can be ruled out at the \(5\%\) level. You would, however, rule out \(\mathbb{E}[X_{i}]=5\) and \(\mathbb{E}[X_{i}]=10\), since they fall outside the interval. “Failing to reject” a value is therefore a weak statement: it says only that the value is one of many still compatible with the data, not that it is correct.

Sampling Distribution Quantiles

To better understand the link between our theoretical interpretation and our calculations, we will compute the quantiles of the sampling distribution. E.g., We construct a \(90\%\) confidence interval of the mean by taking the \(5^{\text{th}}\) and \(95^{\text{th}}\) percentiles of the sampling distribution of the mean. We draw \(300\) samples of \(n=50\) observations from the population, and then construct a confidence interval using quantiles.

# Create 300 samples, each with 50 workers
wage_samples <- matrix(nrow=300, ncol=50)
for(i in seq(1, nrow(wage_samples))){
    wage_samples[i, ] <- sample(population_wages, 50, replace=FALSE)
}
wage_means <- apply(wage_samples, 1, mean) # mean for each sample (row)

# Middle 90%
mq <- quantile(wage_means, probs=c(.05, .95))
paste0('we are 90% confident that the mean is between ', 
    round(mq[1], 2), ' and ', round(mq[2], 2) )
## [1] "we are 90% confident that the mean is between 5.07 and 6.5"

hist(wage_means,
    breaks=seq(3, 9, by=.1),
    border=NA, freq=FALSE,
    main=NA)
title('90% Confidence Interval for the Mean', font.main=1)
abline(v=mq, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=population_wage_mean, col=rgb(0, 0, 1, .8), lty=2)

NoteMust Know

Here is the same calculation for \(300\) uniform random samples of size \(n=1000\).

# Create 300 samples, each with 1000 random uniform variables
X_samples <- matrix(nrow=300, ncol=1000)
for(i in seq(1, nrow(X_samples))){
    X_samples[i, ] <- runif(1000)
}
sample_means <- apply(X_samples, 1, mean) # mean for each sample (row)

# Middle 90%
mq <- quantile(sample_means, probs=c(.05, .95))
paste0('we are 90% confident that the mean is between ',
    round(mq[1], 2), ' and ', round(mq[2], 2) )
## [1] "we are 90% confident that the mean is between 0.49 and 0.51"
# The interval runs from about 0.486 to 0.515, a width of about 0.03.

hist(sample_means,
    breaks=seq(.4, .6, by=.001),
    border=NA, freq=FALSE,
    main=NA)
title('90% Confidence Interval for the Mean', font.main=1)
abline(v=mq, col=rgb(0, 0, 0, .8), lwd=2)

TipTest Yourself

For another example, consider the median. We now repeat the above process to estimate the median for each sample, instead of the mean.

## Sample Quantiles (medians)
sample_quants <- apply(wage_samples, 1, quantile, probs=0.5) #quantile for each sample (row)

# Middle 90% of estimates
mq <- quantile(sample_quants, probs=c(.05, .95))
paste0('we are 90% confident that the median is between ', 
    round(mq[1], 2), ' and ', round(mq[2], 2) )
## [1] "we are 90% confident that the median is between 4.56 and 6.11"

hist(sample_quants,
    breaks=seq(3, 9, by=.1),
    border=NA, freq=FALSE,
    main=NA)
title('90% Confidence Interval for the Median', font.main=1)
abline(v=mq, col=rgb(0, 0, 0, .8), lwd=2)

Interval Size

Confidence intervals shrink with more data, as averaging washes out random fluctuations.

## MAIN EXAMPLE
# Create 300 samples of workers, each of size n
par(mfrow=c(1, 3))
for(n in c(25, 100, 250)){
    # Draw 300 samples
    wage_samples <- matrix(nrow=300, ncol=n)
    for(i in seq(1, nrow(wage_samples))){
        wage_samples[i, ] <- sample(population_wages, n, replace=FALSE)
    }

    # Compute means for each sample (row)
    wage_means <- apply(wage_samples, 1, mean)
    # 90% Confidence Interval
    mq <- quantile(wage_means, probs=c(.05, .95))
    hist(wage_means,
        breaks=seq(3, 9, by=.05),
        border=NA, freq=FALSE,
        main=NA)
    title(paste0('n=', n), font.main=1)
    abline(v=mq, col=rgb(0, 0, 0, .8), lwd=2)
}

NoteMust Know

Here is the same comparison for uniform random samples.

# Create 300 samples, each of size n
par(mfrow=c(1, 3))
for(n in c(25, 100, 250)){
    # Simulate 300 samples
    X_samples <- matrix(nrow=300, ncol=n)
    for(i in seq(1, nrow(X_samples))){
        X_samples[i, ] <- runif(n)
    }

    # Compute means for each sample (row)
    sample_means <- apply(X_samples, 1, mean)
    # 90% Confidence Interval
    mq <- quantile(sample_means, probs=c(.05, .95))
    hist(sample_means,
        breaks=seq(.1, .9, by=.005),
        border=NA, freq=FALSE,
        main=NA)
    title(paste0('n=', n), font.main=1)
    abline(v=mq, col=rgb(0, 0, 0, .8), lwd=2)
}

Confidence intervals also shrink with lower probabilities of containing the theoretical value. For a fixed sample size \(n\), there is a trade-off between precision: the width of a confidence interval, and accuracy: the probability that a confidence interval contains the theoretical value.

TipTest Yourself

Here is an example showing what a Normal CI estimates. Theoretically, \(\pm 1 \cdot \text{SE}(M)\) has 2/3 coverage. Change to \(\pm 2 \cdot \text{SE}(M)\) to see Precision-Accuracy tradeoff.

# Normal CI for each sample
# Uses the last wage_samples above (n=250)
xq <- apply(wage_samples, 1, function(x_row_hat){ #theoretical se's
    mean(x_row_hat) + c(-1, 1)*sd(x_row_hat)/sqrt(length(x_row_hat))
})
# First 3 interval estimates
xq[, c(1, 2, 3)]
##          [,1]     [,2]     [,3]
## [1,] 5.386284 5.803504 5.227066
## [2,] 5.762652 6.249468 5.671207

# Explicit calculation
# Logical vector: whether the population mean is in each CI
covered <- population_wage_mean >= xq[1, ] & population_wage_mean <= xq[2, ]
# Empirical coverage rate
coverage_rate <- mean(covered)
cat(sprintf('Estimated coverage probability: %.2f%%\n', 100 * coverage_rate))
## Estimated coverage probability: 72.00%

7.4 Invert a CI

One way to test a hypothesis is to look at whether a confidence interval contains the hypothesized value: if the value lies outside the interval, sampling variability alone cannot easily explain the gap. This gives a decision rule:

• reject the hypothesis if the hypothesized value falls outside the interval
• fail to reject the hypothesis if it falls inside

A \(95\%\) confidence interval leaves \(5\%\) of the sampling distribution outside it, split between the two tails. So a test at level \(\alpha\) uses an interval with coverage \(1-\alpha\).

For example, suppose you hypothesize the mean equals \(m_{0}\). You then construct a bootstrap distribution with \(95\%\) confidence interval, and find your hypothesized value falls outside of the confidence interval. Then, after accounting for sampling variability (which you estimate), it still seems extremely unlikely that the theoretical mean actually equals \(m_{0}\), so you reject that hypothesis. (If the theoretical value landed in the interval, you would “fail to reject” the theoretical mean equals \(m_{0}\).)

## MAIN EXAMPLE
# x_hat (the sample of 50 workers) and boot_means defined above
ci_95 <- quantile(boot_means, probs=c(.025, .975))
ci_95
##    2.5%   97.5% 
## 4.96511 6.46076

hist(boot_means, breaks=25,
    border=NA,
    freq=FALSE,
    main=NA,
    xlab='Bootstrap means of hourly wages (1980 dollars)')
abline(v=ci_95, col=rgb(0, 0, 0, .8), lwd=2)
# Hypothesized mean: 6.5
abline(v=6.5, col=rgb(0, 0, 1, .8), lty=2)

# The 95% interval runs from about 4.97 to 6.46 dollars an hour.
# The hypothesized 6.5 falls just above it, so we reject this hypothesis.
# The population mean is about 5.76, so this hypothesis is in fact false.

TipTest Yourself

Repeat this test with the murder arrest rates in x3_hat, for the hypothesized mean of \(9\) arrests per \(100{,}000\).

# x3_hat and boot_means3 defined above
ci3_95 <- quantile(boot_means3, probs=c(.025, .975))
ci3_95
##   2.5%  97.5% 
## 6.6099 8.9720

hist(boot_means3, breaks=25,
    border=NA,
    freq=FALSE,
    main=NA,
    xlab='Bootstrap Samples')
abline(v=ci3_95, col=rgb(0, 0, 0, .8), lwd=2)
# Hypothesized mean: 9
abline(v=9, col=rgb(0, 0, 1, .8), lty=2)

# The hypothesized 9 falls just above the 95% interval, so we reject it.

TipTest Yourself

The above procedure also generalizes to many other statistics. Perhaps the most informative are those for spread or shape. E.g., you can conduct hypothesis tests for sd and IQR, or skew and kurtosis.

# Bootstrap Distribution for SD
s_adj_hat <- sd(x_hat)
boot_sds <- rep(NA, 999)
for(b in seq_along(boot_sds)){
    x_boot <- sample(x_hat, replace=TRUE)
    s_adj_b_hat <- sd(x_boot)
    boot_sds[b] <- s_adj_b_hat
}

# Test for SD Differences (Invert CI)
sd_0 <- 2.2 # hypothesized SD
hist(boot_sds, freq=FALSE, breaks=50,
    border=NA, xlab='Bootstrap',
    main=NA)
title('Standard Deviations (Invert CI)', font.main=1)
sd_ci <- quantile(boot_sds, probs=c(0.025, .975) )
abline(v=sd_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=sd_0, col=rgb(0, 0, 1, .8), lty=2)

# The hypothesized 2.2 falls just inside the 95% interval, so we fail to reject it.

To better your understanding, try redoing the above for any function (such as IQR(x_boot)/median(x_boot))

NoteMust Know

Suppose a student scores \(77\%\) on a \(44\)-question multiple-choice exam but believes they are really an \(85\%\) student who had a bad day. The professor is skeptical: “where’s your evidence?” The student bootstraps a confidence interval for their “true” proportion correct, resampling from the \(44\) right/wrong answers on that one exam.

This creates a conundrum. If the confidence interval comes out narrow, it may exclude \(85\%\), and the \(85\%\) hypothesis is rejected. But an interval wide enough to contain \(85\%\) is also wide enough to contain \(70\%\), so the same evidence fails to rule out a worse student who simply got lucky. A single exam cannot easily distinguish the two stories.

We compute the student’s interval and mark both \(70\%\) and \(85\%\) on it.

# One exam: 34 correct out of 44 questions
n_questions <- 44
n_correct <- 34
outcomes <- c(rep(1, n_correct), rep(0, n_questions - n_correct)) # 1=correct, 0=incorrect

# Bootstrap Distribution for the proportion correct
boot_props <- rep(NA, 999)
for(b in seq_along(boot_props)){
    x_boot <- sample(outcomes, replace=TRUE)
    boot_props[b] <- mean(x_boot)
}

hist(boot_props, breaks=25,
    border=NA, freq=FALSE,
    main=NA, xlab='Bootstrap Proportion Correct')
ci_95 <- quantile(boot_props, probs=c(.025, .975))
abline(v=ci_95, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=c(.70, .85), col=rgb(0, 0, 1, .8), lty=2)

# The 95% interval runs from about 0.64 to 0.89.
# It comfortably contains both 0.70 and 0.85.
# So the student can reject neither story.

The deeper issue is what the bootstrap resample actually represents. Resampling only reshuffles the same \(44\) questions, so it approximates how the score would vary if the student retook that exact exam. It says nothing about how the student would do on a different exam with different questions. The student needs that sampling distribution to make their case. Estimating that sampling distribution properly needs new, independent samples: more exams, not more resamples of one exam, and that is a cost the student is unlikely to accept for the rest of the semester.

7.5 Misc. Topics

Prediction Intervals

Sometimes we want to know what range a single new observation will likely fall in.

ImportantKey Definition

A prediction interval is a range that aims to contain a future single observation, not a statistic. We compute it from the lower and upper quantiles of the data themselves, so its width reflects the spread of individual observations rather than the spread of a sampling distribution.

Prediction intervals are useful when a decision depends on individual outcomes. For example, a firm may want to forecast whether a single warranty claim will fall in a tolerable range. Note that \(Z\%\) confidence intervals do not generally cover \(Z\%\) of the data. Compare the \(90\%\) bootstrap confidence interval for the mean with the middle \(90\%\) of the observations themselves.

## MAIN EXAMPLE
# 90% CI for the mean, from the bootstrap distribution above
boot_ci_90 <- quantile(boot_means, probs=c(.05, .95))
boot_ci_90
##       5%      95% 
## 5.086016 6.335946
# Share of sampled workers with a wage inside the CI for the mean
mean(x_hat >= boot_ci_90[1] & x_hat <= boot_ci_90[2])
## [1] 0.14
# Only about 14% of these workers have a wage in the CI, roughly [5.09, 6.34].

# Middle 90% of values: prediction interval
xq <- quantile(x_hat, probs=c(.05, .95))
paste0('we are 90% confident that a new observation will be between ',
    round(xq[1], 2), ' and ', round(xq[2], 2) )
## [1] "we are 90% confident that a new observation will be between 1.79 and 10.75"
# The prediction interval, roughly [1.79, 10.75], is much wider than that CI.

hist(x_hat,
    breaks=seq(0, 16, by=1), border=NA,
    freq=FALSE, main=NA,
    xlab='Hourly wage (1980 dollars)')
rug(x_hat, col=grey(0, .5))
title('Prediction Interval vs CI for the Mean', font.main=1)
abline(v=xq, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=boot_ci_90, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Prediction interval', 'CI for the mean'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

The prediction interval is much wider because it reflects the spread of individual observations rather than the spread of the sample mean.

NoteMust Know

Here is the same comparison for the second sample, x2_hat.

# x2_hat (hourly wages of the 1569 women), m2_hat, and boot_means2 defined above
# 90% CI for the mean
boot_ci2_90 <- quantile(boot_means2, probs=c(.05, .95))
boot_ci2_90
##       5%      95% 
## 5.027807 5.263843
# Share of women with a wage inside the CI for the mean
mean(x2_hat >= boot_ci2_90[1] & x2_hat <= boot_ci2_90[2])
## [1] 0.03696622
# Only about 4% of women earn a wage inside the CI, roughly [5.03, 5.26].

# Middle 90% of values: prediction interval
xq2 <- quantile(x2_hat, probs=c(.05, .95))
xq2
##        5%       95% 
##  1.508883 10.105341
# The prediction interval is roughly [1.51, 10.11].

# Distance from the mean to each end of the prediction interval
c(below=m2_hat - xq2[[1]], above=xq2[[2]] - m2_hat)
##    below    above 
## 3.638041 4.958417
# It reaches about 3.64 dollars below the mean but 4.96 above it.
# It is not centered on the mean because wages are right-skewed.

hist(x2_hat,
    breaks=seq(0, 40, by=.5), border=NA,
    freq=FALSE, main=NA,
    xlab='Hourly wage (1980 dollars)')
title('Prediction Interval vs CI for the Mean', font.main=1)
abline(v=xq2, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=boot_ci2_90, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Prediction interval', 'CI for the mean'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

When the data are skewed, a prediction interval is also not centered on the mean.

TipTest Yourself

Repeat the comparison with the murder arrest rates in x3_hat.

# x3_hat and boot_means3 defined above
# 90% CI for the mean, from the bootstrap distribution above
boot_ci3_90 <- quantile(boot_means3, probs=c(.05, .95))
boot_ci3_90
##     5%    95% 
## 6.8098 8.7760
# Share of states with a murder rate inside the CI for the mean
mean(x3_hat >= boot_ci3_90[1] & x3_hat <= boot_ci3_90[2])
## [1] 0.12
# Only about 12% of states have a murder rate inside the CI, roughly [6.8, 8.8].

# Middle 90% of values: prediction interval
xq3 <- quantile(x3_hat, probs=c(.05, .95))
paste0('we are 90% confident that a new observation will be between ',
    round(xq3[1], 2), ' and ', round(xq3[2], 2) )
## [1] "we are 90% confident that a new observation will be between 2.14 and 15.4"
# The prediction interval, roughly [2.1, 15.4], is much wider than that CI.
# It reflects the spread of individual states, not of the sample mean.

hist(x3_hat,
    breaks=seq(0, 20, by=1), border=NA,
    freq=FALSE, main=NA,
    xlab='Murder arrests per 100,000')
rug(x3_hat, col=grey(0, .5))
title('Prediction Interval vs CI for the Mean', font.main=1)
abline(v=xq3, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=boot_ci3_90, col=rgb(0, 0, 0, .8), lwd=2, lty=2)
legend('topright', legend=c('Prediction interval', 'CI for the mean'),
    lty=c(1, 2), lwd=2, col=rgb(0, 0, 0, .8), bty='n')

NoteMust Know

For uniform random data, the \(90\%\) confidence interval for the mean from the quantile example above theoretically covers far less than \(0.9\) of the data.

x_future_hat <- runif(1000)
# Middle 90% of values
xq0 <- quantile(x_future_hat, probs=c(.05, .95))
paste0('we are 90% confident that a future data point will be between ',
    round(xq0[1], 2), ' and ', round(xq0[2], 2) )
## [1] "we are 90% confident that a future data point will be between 0.06 and 0.95"
# By contrast, the 90% CI for the mean above has a width of about 0.03.
# So it theoretically covers only about 0.03 of the data.

hist(x_future_hat,
    breaks=seq(0, 1, by=.05), border=NA,
    freq=FALSE, main=NA)
rug(x_future_hat, col=grey(0, .5))
title('Prediction Interval', font.main=1)
abline(v=xq0, col=rgb(0, 0, 0, .8), lwd=2)

One-Sided Intervals

Above, our confidence intervals were two-sided: they contained the middle \(Z\%\) of the sampling distribution. We can also construct one-sided intervals that extend to infinity in one direction.

A one-sided interval is shifted to one side, containing one tail rather than the middle. For example, an upper-bounded interval uses \((-\infty, q_{0.95}]\), where \(q_{0.95}\) is the \(95^{\text{th}}\) percentile of the bootstrap distribution: we are \(95\%\) confident the true value is at most \(q_{0.95}\). A lower-bounded interval uses \([q_{0.05}, \infty)\), where \(q_{0.05}\) is the \(5^{\text{th}}\) percentile of the bootstrap distribution: we are \(95\%\) confident the true value is at least \(q_{0.05}\).

TipTest Yourself

For the electricity bill example above with \(\hat{m}=142\) and \(\text{SE}=5\), a one-sided \(95\%\) upper-bounded confidence interval is \[(-\infty,~ \hat{m} + 1.645 \times \text{SE}] = (-\infty,~ 142 + 1.645 \times 5] = (-\infty,~ 150.2].\] We are \(95\%\) confident the population mean bill is at most \(\$150.20\). Note that \(1.645\) is the \(95^{\text{th}}\) percentile of the standard normal (compared to \(1.96\) for the two-sided case).

m_bill_hat <- 142
SE <- 5
# One-sided 95% upper bound
m_bill_hat + qnorm(0.95) * SE
## [1] 150.2243

TipTest Yourself

Here is a bootstrap right-tail example

# Bootstrap Distribution
boot_means <- rep(NA, 9999)
for(b in seq_along(boot_means)){
    x_boot <- sample(x_hat, replace=TRUE) # c.f. jackknife
    m_b_hat <- mean(x_boot)
    boot_means[b] <-m_b_hat
}

# One-sided 95% lower bound
# We are 95% confident the mean is at least q_05
hist(boot_means, border=NA, breaks=50,
    freq=FALSE, main=NA, xlab='Bootstrap')
ci_95 <- quantile(boot_means, probs=c(0.05, 1) )
abline(v=ci_95, col=rgb(0, 0, 0, .8), lwd=2)

7.6 Exercises

1. Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.

2. A 95% confidence interval does not mean there is a 95% probability the true parameter lies inside that particular interval. Explain what “95% confident” actually means in terms of repeated sampling.

3. Suppose a sample of \(n = 49\) light bulbs has a mean lifetime of \(\hat{m} = 1200\) hours with an adjusted sample standard deviation of \(\hat{s}_{\text{adj}} = 140\) hours. Compute the 90% and 99% Normal-approximation confidence intervals for the population mean lifetime. Which interval is wider, and why?

4. Using the USArrests dataset in R, construct a bootstrap 95% confidence interval for the median of the UrbanPop variable. Use at least 9999 bootstrap resamples and plot the bootstrap distribution with the interval boundaries.

Further Reading

• https://learningstatisticswithr.com/book/estimation.html#interval-estimates – Confidence intervals explained with intuition and R code.
• https://www.khanacademy.org/math/statistics-probability/confidence-intervals-one-sample – Video walkthroughs of confidence interval construction and interpretation.
• https://www.econometrics-with-r.org/3.1-hypothesis-tests-and-confidence-intervals-in-the-simple-linear-regression-model.html – Confidence intervals applied to regression, bridging descriptive and inferential statistics.

Recall

We built confidence intervals in two ways, with the bootstrap percentile method and the Normal approximation. Being “\(Z\%\) confident” means that \(Z\%\) of intervals built on different samples will contain the true value. For a fixed sample size, there is a trade-off between precision and accuracy. Inverting a CI tests a hypothesized value: reject it if it falls outside the interval, and fail to reject it otherwise. A prediction interval aims to contain a single new observation, so it uses quantiles of the data themselves and is much wider than a confidence interval for the mean. In the electricity-bill example, \(\hat{m}=142\), \(\hat{s}_{\text{adj}}=30\), and \(n=36\). The classical standard error is \(\text{SE}(M) = 30/\sqrt{36} = 5\). The \(95\%\) margin of error is \(E = 1.96 \times 5 = 9.8\), giving a confidence interval of \([\$132.20, \$151.80]\).

---

1. Generally, a coverage level of \(1-\alpha\) means \(\text{Prob}( M - E < \mathbb{E}[X_{i}] < M + E)=1-\alpha\). Notice that \(\text{Prob}( M - E < \mathbb{E}[X_{i}] < M + E) = \text{Prob}( - E < \mathbb{E}[X_{i}] - M < + E) = \text{Prob}( \mathbb{E}[X_{i}] + E > M > \mathbb{E}[X_{i}] - E)\). So if the interval \([\mathbb{E}[X_{i}] - 10, \mathbb{E}[X_{i}] + 10]\) contains \(95\%\) of all \(M\), then the interval \([M-10, M+10]\) will also contain \(\mathbb{E}[X_{i}]\) in \(95\%\) of the samples. This is because whenever \(M\) is within \(10\) of \(\mathbb{E}[X_{i}]\), the value \(\mathbb{E}[X_{i}]\) is also within \(10\) of \(M\). But for any particular sample, the interval \([\hat{m}-10, \hat{m}+10]\) either does or does not contain \(\mathbb{E}[X_{i}]\). Similarly, if you compute \(\hat{m}=9\) for your particular sample, a coverage level of \(1-\alpha=95\%\) does not mean \(\text{Prob}(9 - E < \mathbb{E}[X_{i}] < 9 + E)=95\%\).

2. As such, a \(Z\%\) confidence interval does not imply a \(Z\%\) probability that the true parameter lies within a particular calculated interval. In practice, people often interpret confidence intervals informally as “showing the uncertainty around our estimate”: wider intervals correspond to higher sampling variability and less precise information about \(\mathbb{E}[X_{i}]\).