R - Confidence Intervals | Application Architect

Key Insights

• Confidence intervals quantify estimation uncertainty by providing a range of plausible values for population parameters, with the 95% level being standard practice in most fields • R provides multiple approaches for calculating confidence intervals: analytical methods using t.test() and prop.test(), bootstrap resampling via boot package, and Bayesian credible intervals through rstanarm or brms • The choice between methods depends on sample size, distribution assumptions, and parameter type—bootstrap methods offer robust alternatives when normality assumptions fail

Understanding Confidence Intervals

A confidence interval provides a range of values that likely contains an unknown population parameter. A 95% confidence interval means that if you repeated your sampling process infinitely, 95% of the constructed intervals would contain the true parameter value.

The general form for a confidence interval is:

estimate ± (critical value × standard error)

For a mean with known population standard deviation, this becomes:

# Manual calculation for population mean
sample_mean <- 25.3
pop_sd <- 5.2
n <- 50
z_critical <- qnorm(0.975)  # 1.96 for 95% CI

se <- pop_sd / sqrt(n)
margin_error <- z_critical * se
ci_lower <- sample_mean - margin_error
ci_upper <- sample_mean + margin_error

cat(sprintf("95%% CI: [%.2f, %.2f]\n", ci_lower, ci_upper))

Confidence Intervals for Means

When the population standard deviation is unknown (the typical case), use the t-distribution instead of the normal distribution. R’s t.test() function handles this automatically:

# Generate sample data
set.seed(123)
data <- rnorm(30, mean = 100, sd = 15)

# Calculate CI for mean
result <- t.test(data)
print(result)

# Extract just the confidence interval
ci <- result$conf.int
cat(sprintf("95%% CI: [%.2f, %.2f]\n", ci[1], ci[2]))

For different confidence levels, adjust the conf.level parameter:

# 90% confidence interval
t.test(data, conf.level = 0.90)$conf.int

# 99% confidence interval
t.test(data, conf.level = 0.99)$conf.int

When comparing two groups, t.test() provides the CI for the difference in means:

# Two independent samples
group1 <- rnorm(40, mean = 75, sd = 10)
group2 <- rnorm(35, mean = 80, sd = 12)

# Two-sample t-test
comparison <- t.test(group1, group2)
cat(sprintf("Difference CI: [%.2f, %.2f]\n", 
            comparison$conf.int[1], 
            comparison$conf.int[2]))

# Paired samples
before <- c(120, 135, 128, 142, 130)
after <- c(115, 128, 122, 138, 125)

paired_result <- t.test(before, after, paired = TRUE)
print(paired_result$conf.int)

Confidence Intervals for Proportions

For binary outcomes, use prop.test() to calculate confidence intervals for proportions:

# Single proportion
successes <- 45
total <- 100

prop_result <- prop.test(successes, total)
cat(sprintf("Proportion: %.3f\n", prop_result$estimate))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", 
            prop_result$conf.int[1], 
            prop_result$conf.int[2]))

The default method uses the Wilson score interval, which performs better than the Wald interval for small samples:

# Comparing two proportions
successes <- c(30, 45)
totals <- c(100, 120)

two_prop <- prop.test(successes, totals)
print(two_prop$conf.int)  # CI for difference in proportions

For exact binomial confidence intervals with very small samples:

# Exact binomial test
binom.test(12, 20, conf.level = 0.95)$conf.int

Bootstrap Confidence Intervals

Bootstrap resampling provides confidence intervals without distributional assumptions. This approach works well for medians, trimmed means, or any custom statistic:

library(boot)

# Define statistic function
median_func <- function(data, indices) {
  return(median(data[indices]))
}

# Generate data
set.seed(456)
sample_data <- rexp(50, rate = 0.1)

# Bootstrap
boot_results <- boot(sample_data, median_func, R = 10000)

# Calculate different types of bootstrap CIs
boot.ci(boot_results, type = c("norm", "basic", "perc", "bci"))

The percentile method often performs best:

# Extract percentile CI
percentile_ci <- boot.ci(boot_results, type = "perc")$percent[4:5]
cat(sprintf("Bootstrap 95%% CI: [%.2f, %.2f]\n", 
            percentile_ci[1], 
            percentile_ci[2]))

Bootstrap for more complex statistics like correlation coefficients:

# Correlation coefficient CI
correlation_func <- function(data, indices) {
  d <- data[indices, ]
  return(cor(d$x, d$y))
}

# Create bivariate data
df <- data.frame(
  x = rnorm(60),
  y = rnorm(60)
)
df$y <- df$y + 0.5 * df$x  # Add correlation

boot_cor <- boot(df, correlation_func, R = 5000)
boot.ci(boot_cor, type = "perc")

Confidence Intervals for Regression

Linear models provide confidence intervals for coefficients and predictions:

# Linear regression
set.seed(789)
x <- 1:50
y <- 3 + 2*x + rnorm(50, sd = 5)
model <- lm(y ~ x)

# Coefficient CIs
confint(model)

# Prediction intervals vs confidence intervals
new_data <- data.frame(x = c(25, 30, 35))

# Confidence interval for mean response
predict(model, new_data, interval = "confidence", level = 0.95)

# Prediction interval for individual observations
predict(model, new_data, interval = "prediction", level = 0.95)

The prediction interval is wider because it accounts for both estimation uncertainty and individual variation:

# Visualize both intervals
plot_data <- data.frame(x = seq(min(x), max(x), length.out = 100))
conf_int <- predict(model, plot_data, interval = "confidence")
pred_int <- predict(model, plot_data, interval = "prediction")

plot(x, y, pch = 16, col = "gray50")
abline(model, col = "blue", lwd = 2)
lines(plot_data$x, conf_int[, "lwr"], col = "red", lty = 2)
lines(plot_data$x, conf_int[, "upr"], col = "red", lty = 2)
lines(plot_data$x, pred_int[, "lwr"], col = "green", lty = 3)
lines(plot_data$x, pred_int[, "upr"], col = "green", lty = 3)
legend("topleft", 
       legend = c("Regression", "Confidence", "Prediction"),
       col = c("blue", "red", "green"),
       lty = c(1, 2, 3))

Non-Parametric Confidence Intervals

When data violates normality assumptions, use non-parametric methods:

# Wilcoxon signed-rank test for median
set.seed(101)
skewed_data <- rexp(40, rate = 0.5)

wilcox_result <- wilcox.test(skewed_data, conf.int = TRUE)
cat(sprintf("Median CI: [%.2f, %.2f]\n", 
            wilcox_result$conf.int[1], 
            wilcox_result$conf.int[2]))

For custom percentiles, bootstrap remains the most flexible approach:

# 90th percentile CI
percentile_90 <- function(data, indices) {
  return(quantile(data[indices], 0.90))
}

boot_p90 <- boot(skewed_data, percentile_90, R = 10000)
boot.ci(boot_p90, type = "perc")

Practical Considerations

Sample size directly affects interval width. Demonstrate this relationship:

# CI width vs sample size
sample_sizes <- c(10, 30, 50, 100, 200, 500)
ci_widths <- sapply(sample_sizes, function(n) {
  sample <- rnorm(n, mean = 50, sd = 10)
  ci <- t.test(sample)$conf.int
  return(ci[2] - ci[1])
})

plot(sample_sizes, ci_widths, type = "b", 
     xlab = "Sample Size", ylab = "CI Width",
     main = "Confidence Interval Width vs Sample Size")

Always verify assumptions before choosing a method. For means, check normality with small samples:

# Check normality
shapiro.test(data)

# Q-Q plot
qqnorm(data)
qqline(data, col = "red")

When reporting confidence intervals, include the confidence level, sample size, and method used. Confidence intervals provide more information than p-values alone and support better decision-making in applied settings.