Statistical Inference: confidence intervals
From Description to Inference
library(tidyverse)
set.seed(301)
# Build 100 confidence intervals from 100 different samples
cis <- map_dfr(1:100, function(i) {
s <- rnorm(30, mean = 170, sd = 10)
tt <- t.test(s)
tibble(trial = i, est = mean(s),
lower = tt$conf.int[1], upper = tt$conf.int[2])
}) %>%
mutate(captures = lower <= 170 & upper >= 170)
cat("Intervals containing the true mean:", sum(cis$captures), "out of 100\n")Intervals containing the true mean: 95 out of 100ggplot(cis, aes(trial, est, color = captures)) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
geom_point(size = 1) +
geom_hline(yintercept = 170, color = "black", linewidth = 0.8) +
scale_color_manual(values = c("FALSE" = "#E74C3C", "TRUE" = "#27AE60")) +
labs(title = "100 samples, 100 confidence intervals",
subtitle = "Red intervals missed the true mean of 170",
x = "Sample number", y = "Mean height (cm)", color = "Captured?") +
theme_minimal(base_size = 11)
95% confidence means that if we were to repeat this process many times, about 95% of the resulting confidence intervals would contain the true mean. Here, five of the 100 confidence intervals missed the true mean of 170.
Important
A 95% CI is a range built by a procedure that captures the true parameter 95% of the time in repeated sampling.
Letโs Try It Out
We calculate a 95% confidence interval for the population mean height using a random sample of 30 people. It uses the t-distribution because the population standard deviation is unknown and is estimated using the sample data.
Calculating a confidence interval is a four-step process:
set.seed(301)
# Generate a random sample of 30 heights (in cm) from a normal distribution
heights <- rnorm(30, mean = 170, sd = 10) # one random normal sample
# Compute the sample mean and standard error
# heights refers to the 30 simulated observations
x_bar <- mean(heights) # avg of 30 sample means
s <- sd(heights) # sample standard deviation
n <- length(heights) # Find the sample size
se <- s / sqrt(n) # Find the standard error of the mean
# Note: The standard error estimates how much the
# sample mean would vary across repeated random samples
# of the same size.
# Find the critical value from the t-distribution.
# The sigma value is unknown, so we use the
# t-distribution with n-1 degrees of freedom
# How to find the critical value? For a 95%
# confidence level, we want the value that leaves 2.5% in each tail.
t_crit <- qt(0.975, df = n - 1)
# We calculate the lower and upper bounds of the 95%
# confidence interval
# The interval provides a range of plausible values for the
# population mean, based on the sample data and the assumptions
# of the t-interval method.
ci <- c(x_bar - t_crit * se, x_bar + t_crit * se)
round(c(mean = x_bar, se = se, lower = ci[1], upper = ci[2]), 2) mean se lower upper
169.33 1.78 165.68 172.97 Tip
The standard error (se) estimates how much the sample mean would vary across repeated random samples of the same size.
Some ideas about confidence intervals
Purpose: A confidence interval estimates a plausible range for the unknown population mean height.
A 95% confidence interval means that if we were to repeat this process many times, about 95% of the resulting confidence intervals would contain the true mean. Here, five of the 100 confidence intervals missed the true mean of 170.
Confidence level: A 95% confidence interval is produced by using the 0.975 quantile of the t-distribution.
Important
[1] 165.6761 172.9745
attr(,"conf.level")
[1] 0.95[1] 164.4073 174.2433
attr(,"conf.level")
[1] 0.99[1] 0.3233236 0.5228954
attr(,"conf.level")
[1] 0.95Tip
Wider confidence = wider interval. A 99% CI is more likely to be right and less informative. A measure of certainty of the estimate is traded off against precision of the estimate.
โ Wrong
โThere is a 95% probability that ฮผ is between 167 and 173.โ
The true ฮผ is a fixed number. It is either in your interval or it is not โ no probability involved. ๐ฒ
โ Right
โWe are 95% confidentโฆโ meaning: this method produces intervals that capture ฮผ 95% of the time.
The randomness lives in the interval, not in ฮผ.
The Practical Reading
A CI tells you which values are compatible with your data. If a claimed value falls outside your interval, your data argues against it. ๐
[-0.5, 8.2] points.
What is the best summary?
library(tidyverse)
set.seed(301)
# One trial = sum of five dice; replicate it 10,000 times
totals <- replicate(10000, sum(sample(1:6, 5, replace = TRUE)))
# Simulated probability
p_sim <- mean(totals >= 20)
p_sim
# Normal approximation:
# one die has mean 3.5 and variance 35/12
mu <- 5 * 3.5
sigma <- sqrt(5 * 35 / 12)
# Continuity correction: P(X >= 20) uses 19.5
p_norm <- pnorm(19.5, mu, sigma, lower.tail = FALSE)
c(simulated = p_sim, normal_approx = p_norm)
# The CLT at work: five dice already look bell-shaped
ggplot(tibble(t = totals), aes(t)) +
geom_bar(fill = "#8E44AD") +
theme_minimal()
Note
| Task | Code |
|---|---|
| Random sample | sample(x, n, replace = TRUE) |
| Reproducibility | set.seed(301) |
| Density / mass | dnorm(), dbinom(), dpois() |
| P(X โค q) | pnorm(q, mean, sd) |
| P(X > q) | pnorm(q, ..., lower.tail = FALSE) |
| Quantile | qnorm(p, mean, sd) |
| Random draws | rnorm(n, mean, sd) |
| Repeat a simulation | replicate(n, expr) |
| Standard error | sd(x) / sqrt(length(x)) |
| CI for a mean | t.test(x)$conf.int |
| CI for a proportion | prop.test(successes, n)$conf.int |
| t critical value | qt(0.975, df = n - 1) |