1 Today’s Question

In Lectures 2, 3, and 4, we estimated the adjusted difference in CVD risk between smokers and non-smokers.

Those lectures gave us point estimates:

one best guess from one dataset.

Lecture 5 asks a new question:

How uncertain are those estimates?

We will use the bootstrap to build 95% confidence intervals and do hypothesis testing for:

  1. IPW
  2. Outcome regression
  3. Doubly robust estimation

Our teaching analysis still uses the same variables and does not use survey weights.

Symbol Meaning in this lecture
T Smoking status: 1 = smoker, 0 = non-smoker
Y CVD status: 1 = had CVD, 0 = did not have CVD
X Age, gender, race, education, income-to-poverty ratio, and BMI

2 Why Bootstrap?

A point estimate is like taking one photo.

A confidence interval asks:

If we could repeat the study many times, how much might the estimate move around?

But we usually have only one dataset. The bootstrap uses that one dataset to create many “pretend new datasets.”

3 Bootstrap Sampling

In one bootstrap sample:

  1. We randomly draw the same number of people as the original dataset.
  2. We draw with replacement, so one person can appear more than once.
  3. Some people from the original dataset may not appear in that bootstrap sample.
  4. We refit the models and recalculate the adjusted risk difference.

4 Load and Prepare the Data

data_path <- file.path("..", "Data", "NHANES_data.csv")

if (!file.exists(data_path)) {
  data_path <- file.path("Data", "NHANES_data.csv")
}

nhanes <- read.csv(data_path, stringsAsFactors = FALSE)

baseline_vars <- c(
  "age_yr",
  "gender",
  "race",
  "educ_lvl",
  "inc_to_pov_ratio",
  "bmi"
)

analysis_vars <- c("smoker_indicator", "cvd_indicator", baseline_vars)

analysis_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]

analysis_data$smoker_indicator <- as.integer(analysis_data$smoker_indicator)
analysis_data$cvd_indicator <- as.integer(analysis_data$cvd_indicator)

analysis_data$smoking_group <- ifelse(
  analysis_data$smoker_indicator == 1,
  "Smoker",
  "Non-smoker"
)

analysis_data$smoking_group <- factor(
  analysis_data$smoking_group,
  levels = c("Non-smoker", "Smoker")
)

analysis_data$gender_label <- factor(
  analysis_data$gender,
  levels = c("F", "M"),
  labels = c("Female", "Male")
)

analysis_data$race_label <- factor(
  analysis_data$race,
  levels = c(
    "mex_american",
    "other_hispanic",
    "nh_white",
    "nh_black",
    "nh_asian"
  ),
  labels = c(
    "Mexican American",
    "Other Hispanic",
    "Non-Hispanic White",
    "Non-Hispanic Black",
    "Non-Hispanic Asian"
  )
)

analysis_data$educ_label <- factor(
  analysis_data$educ_lvl,
  levels = c(
    "lt_9th_grade",
    "9_to_12th_grade_no_diploma",
    "hs_grad_or_ged",
    "start_college_to_aa",
    "college_grad_or_above"
  ),
  labels = c(
    "Less than 9th grade",
    "9th-12th, no diploma",
    "High school or GED",
    "Some college or AA",
    "College graduate or above"
  )
)

nrow(analysis_data)
## [1] 6299
sample_summary <- data.frame(
  Quantity = c(
    "Number of people",
    "Number of smokers",
    "Number of non-smokers",
    "Number with CVD"
  ),
  Value = c(
    nrow(analysis_data),
    sum(analysis_data$smoker_indicator == 1),
    sum(analysis_data$smoker_indicator == 0),
    sum(analysis_data$cvd_indicator == 1)
  )
)

knitr::kable(sample_summary)
Quantity Value
Number of people 6299
Number of smokers 2559
Number of non-smokers 3740
Number with CVD 646

5 The Estimators We Will Bootstrap

Each bootstrap sample repeats the same three adjusted analyses.

weighted_mean <- function(outcome, weights) {
  sum(weights * outcome) / sum(weights)
}

fit_three_methods <- function(input_data) {
  treatment_model <- glm(
    smoker_indicator ~ age_yr + gender_label + race_label +
      educ_label + inc_to_pov_ratio + bmi,
    data = input_data,
    family = binomial()
  )

  e_hat <- predict(treatment_model, type = "response")

  outcome_model <- glm(
    cvd_indicator ~ smoker_indicator + age_yr + gender_label +
      race_label + educ_label + inc_to_pov_ratio + bmi,
    data = input_data,
    family = binomial()
  )

  data_if_non_smoker <- input_data
  data_if_smoker <- input_data

  data_if_non_smoker$smoker_indicator <- 0
  data_if_smoker$smoker_indicator <- 1

  m0_hat <- predict(
    outcome_model,
    newdata = data_if_non_smoker,
    type = "response"
  )

  m1_hat <- predict(
    outcome_model,
    newdata = data_if_smoker,
    type = "response"
  )

  smoker_rows <- input_data$smoker_indicator == 1
  non_smoker_rows <- input_data$smoker_indicator == 0

  ipw_risk_smoker <- weighted_mean(
    input_data$cvd_indicator[smoker_rows],
    1 / e_hat[smoker_rows]
  )

  ipw_risk_non_smoker <- weighted_mean(
    input_data$cvd_indicator[non_smoker_rows],
    1 / (1 - e_hat[non_smoker_rows])
  )

  ipw_difference <- ipw_risk_smoker - ipw_risk_non_smoker

  or_risk_smoker <- mean(m1_hat)
  or_risk_non_smoker <- mean(m0_hat)
  or_difference <- or_risk_smoker - or_risk_non_smoker

  T_obs <- input_data$smoker_indicator
  Y_obs <- input_data$cvd_indicator

  dr_risk_smoker <- mean(m1_hat + T_obs / e_hat * (Y_obs - m1_hat))
  dr_risk_non_smoker <- mean(
    m0_hat + (1 - T_obs) / (1 - e_hat) * (Y_obs - m0_hat)
  )

  dr_difference <- dr_risk_smoker - dr_risk_non_smoker

  c(
    IPW = ipw_difference,
    Outcome_regression = or_difference,
    Doubly_robust = dr_difference
  )
}

safe_fit_three_methods <- function(input_data) {
  tryCatch(
    fit_three_methods(input_data),
    error = function(e) {
      c(IPW = NA_real_,
        Outcome_regression = NA_real_,
        Doubly_robust = NA_real_)
    }
  )
}

6 Point Estimates From the Original Dataset

original_estimates <- fit_three_methods(analysis_data)

original_table <- data.frame(
  Method = c("IPW", "Outcome regression", "Doubly robust"),
  Risk_difference_percentage_points = round(100 * original_estimates, 2)
)

knitr::kable(
  original_table,
  col.names = c("Method", "Risk difference: smoker minus non-smoker")
)
Method Risk difference: smoker minus non-smoker
IPW IPW 2.99
Outcome_regression Outcome regression 3.45
Doubly_robust Doubly robust 3.25

These are the point estimates. Now we will use bootstrapping to estimate their uncertainty.

7 Run the Bootstrap

For this lecture note, we use B = 300 bootstrap samples so the file can knit quickly. In a research project, it is common to use more, such as 1000 or 2000.

set.seed(20260627)

B <- 300
n <- nrow(analysis_data)

bootstrap_estimates <- matrix(
  NA_real_,
  nrow = B,
  ncol = length(original_estimates)
)

colnames(bootstrap_estimates) <- names(original_estimates)

for (b in seq_len(B)) {
  bootstrap_rows <- sample.int(n, size = n, replace = TRUE)
  bootstrap_data <- analysis_data[bootstrap_rows, ]
  bootstrap_estimates[b, ] <- safe_fit_three_methods(bootstrap_data)
}

bootstrap_estimates <- bootstrap_estimates[complete.cases(bootstrap_estimates), ]

number_of_successful_bootstraps <- nrow(bootstrap_estimates)
number_of_successful_bootstraps
## [1] 300

8 Build 95% Confidence Intervals

We use a simple percentile bootstrap interval.

For each method:

  1. Sort the bootstrap estimates.
  2. Take the 2.5th percentile.
  3. Take the 97.5th percentile.
  4. The middle 95% is the confidence interval.
bootstrap_se <- apply(bootstrap_estimates, 2, sd)

bootstrap_ci <- t(apply(
  bootstrap_estimates,
  2,
  quantile,
  probs = c(0.025, 0.975),
  na.rm = TRUE
))

colnames(bootstrap_ci) <- c("Lower", "Upper")

inference_table <- data.frame(
  Method = c("IPW", "Outcome regression", "Doubly robust"),
  Estimate_pp = round(100 * original_estimates, 2),
  Bootstrap_SE_pp = round(100 * bootstrap_se, 2),
  CI_lower_pp = round(100 * bootstrap_ci[, "Lower"], 2),
  CI_upper_pp = round(100 * bootstrap_ci[, "Upper"], 2)
)

knitr::kable(
  inference_table,
  col.names = c(
    "Method",
    "Estimate (pp)",
    "Bootstrap SE (pp)",
    "95% CI lower (pp)",
    "95% CI upper (pp)"
  )
)
Method Estimate (pp) Bootstrap SE (pp) 95% CI lower (pp) 95% CI upper (pp)
IPW IPW 2.99 0.79 1.62 4.56
Outcome_regression Outcome regression 3.45 0.79 2.11 5.02
Doubly_robust Doubly robust 3.25 0.80 1.94 4.92

Here, pp means percentage points.

For example, a risk difference of 3.0 percentage points means:

CVD risk is estimated to be 3.0 percentage points higher under smoking than under non-smoking.

9 Picture of the Bootstrap Distributions

Each histogram shows the bootstrap estimates from many pretend datasets.

If the bootstrap estimates are tightly packed, the method has a smaller standard error. If they are more spread out, the estimate is more uncertain.

10 Picture of the 95% Confidence Intervals

The dashed vertical line means zero difference.

If a confidence interval does not cross zero, that suggests the adjusted difference is unlikely to be exactly zero.

11 Hypothesis Testing

Now we formally ask:

Is there evidence that the adjusted CVD risk is different for smokers and non-smokers?

For each method, our hypotheses are:

H0: adjusted risk difference = 0

H1: adjusted risk difference is not 0

The null hypothesis H0 means:

After adjustment, smoking and non-smoking have the same CVD risk.

The alternative hypothesis H1 means:

After adjustment, smoking and non-smoking have different CVD risk.

12 Bootstrap Standard Error Test

We estimate the standard error using the bootstrap:

z_statistic <- original_estimates / bootstrap_se
p_value <- 2 * (1 - pnorm(abs(z_statistic)))

hypothesis_table <- data.frame(
  Method = c("IPW", "Outcome regression", "Doubly robust"),
  Estimate_pp = round(100 * original_estimates, 2),
  Bootstrap_SE_pp = round(100 * bootstrap_se, 2),
  Z_statistic = round(z_statistic, 2),
  P_value = signif(p_value, 3),
  Conclusion_at_0_05 = ifelse(
    p_value < 0.05,
    "Reject H0",
    "Do not reject H0"
  )
)

knitr::kable(
  hypothesis_table,
  col.names = c(
    "Method",
    "Estimate (pp)",
    "Bootstrap SE (pp)",
    "z",
    "p-value",
    "Conclusion at 0.05"
  )
)
Method Estimate (pp) Bootstrap SE (pp) z p-value Conclusion at 0.05
IPW IPW 2.99 0.79 3.76 1.71e-04 Reject H0
Outcome_regression Outcome regression 3.45 0.79 4.37 1.22e-05 Reject H0
Doubly_robust Doubly robust 3.25 0.80 4.08 4.57e-05 Reject H0

The p-value is a way to measure how surprising our estimate would be if the true adjusted risk difference were really zero.

Small p-values give evidence against H0.

13 Hypothesis Test Picture

14 Interpretation

For this complete-case classroom sample:

interpretation_table <- data.frame(
  Method = c("IPW", "Outcome regression", "Doubly robust"),
  Estimate = paste0(round(100 * original_estimates, 1), " pp"),
  CI_95 = paste0(
    "(",
    round(100 * bootstrap_ci[, "Lower"], 1),
    ", ",
    round(100 * bootstrap_ci[, "Upper"], 1),
    ") pp"
  ),
  P_value = signif(p_value, 3)
)

knitr::kable(
  interpretation_table,
  col.names = c("Method", "Estimate", "95% CI", "p-value")
)
Method Estimate 95% CI p-value
IPW IPW 3 pp (1.6, 4.6) pp 1.71e-04
Outcome_regression Outcome regression 3.5 pp (2.1, 5) pp 1.22e-05
Doubly_robust Doubly robust 3.3 pp (1.9, 4.9) pp 4.57e-05

The IPW, outcome-regression, and doubly robust estimates all compare the predicted CVD risk if everyone were a smoker with the predicted CVD risk if everyone were a non-smoker, after adjusting for measured baseline covariates.

If a method’s 95% confidence interval is entirely above zero, the bootstrap analysis suggests evidence that the adjusted smoker risk is higher than the adjusted non-smoker risk.

But this still depends on causal assumptions:

  1. We adjusted for the important confounders.
  2. Smokers and non-smokers have enough overlap in covariates.
  3. The models used by each method are reasonable enough.
  4. The data were measured well enough for this teaching question.

15 What the Bootstrap Does Not Fix

Bootstrap helps with random sampling uncertainty.

It does not automatically fix:

  1. Unmeasured confounding
  2. Bad measurement
  3. Poor model choices
  4. Lack of overlap
  5. Survey-design issues that we intentionally skipped for this high-school-level lesson

So the bootstrap tells us:

How much our estimate moves around because the sample could have been different.

It does not prove:

Smoking caused the exact difference we estimated.

16 Key Takeaways

  1. A point estimate gives one best guess.
  2. A bootstrap confidence interval shows uncertainty around that guess.
  3. Bootstrapping repeats the whole analysis on many resampled datasets.
  4. We used bootstrapping for IPW, outcome regression, and doubly robust estimation.
  5. Hypothesis testing asks whether zero difference is plausible.
  6. Statistical evidence is not the same as causal proof; causal assumptions still matter.

17 Appendix: Bootstrap Math Used in This Lecture

This appendix gives the mathematical notation for the bootstrap procedure used in this lecture.

For each person \(i = 1, \ldots, n\), let:

\[ O_i = (Y_i, T_i, X_i) \]

where \(Y_i\) is CVD status, \(T_i\) is smoking status, and \(X_i\) is the vector of measured baseline covariates.

For each method \(m\), where:

\[ m \in \{\text{IPW}, \text{OR}, \text{DR}\}, \]

let the estimated adjusted risk difference from the original dataset be:

\[ \hat\Delta_m. \]

In words, \(\hat\Delta_m\) estimates:

\[ \Delta_m = \text{CVD risk if everyone smoked} - \text{CVD risk if everyone did not smoke}. \]

In this lecture, we report this risk difference in percentage points:

\[ 100 \times \hat\Delta_m. \]

17.1 Bootstrap Resampling

The empirical distribution of the observed data is:

\[ \hat F_n = \frac{1}{n} \sum_{i=1}^{n} \delta_{O_i}, \]

where \(\delta_{O_i}\) means a point mass at person \(i\)’s observed data.

A bootstrap sample is created by drawing:

\[ O_1^{*(b)}, O_2^{*(b)}, \ldots, O_n^{*(b)} \overset{iid}{\sim} \hat F_n, \]

for bootstrap repetition \(b = 1, \ldots, B\).

In this lecture:

\[ B = 300. \]

Because we sample from \(\hat F_n\), a person can appear more than once in a bootstrap sample, and another person may not appear at all.

17.2 Bootstrap Estimates

For each bootstrap sample \(b\), we repeat the whole analysis and recompute the risk difference:

\[ \hat\Delta_m^{*(b)} = \text{risk difference estimate from bootstrap sample } b. \]

So for each method \(m\), the bootstrap gives many estimates:

\[ \hat\Delta_m^{*(1)}, \hat\Delta_m^{*(2)}, \ldots, \hat\Delta_m^{*(B)}. \]

These values form the bootstrap distribution.

17.3 Bootstrap Standard Error

The bootstrap mean is:

\[ \bar{\Delta}_m^* = \frac{1}{B} \sum_{b=1}^{B} \hat\Delta_m^{*(b)}. \]

The bootstrap standard error is:

\[ \widehat{SE}_{boot}(\hat\Delta_m) = \sqrt{ \frac{1}{B - 1} \sum_{b=1}^{B} \left( \hat\Delta_m^{*(b)} - \bar{\Delta}_m^* \right)^2 }. \]

This number measures how much the estimate moves across the bootstrap samples.

17.4 Percentile Bootstrap Confidence Interval

The lecture uses a percentile bootstrap confidence interval.

Let:

\[ q_{0.025,m}^* = \text{2.5th percentile of } \left\{ \hat\Delta_m^{*(1)}, \ldots, \hat\Delta_m^{*(B)} \right\}, \]

and:

\[ q_{0.975,m}^* = \text{97.5th percentile of } \left\{ \hat\Delta_m^{*(1)}, \ldots, \hat\Delta_m^{*(B)} \right\}. \]

Then the 95% bootstrap confidence interval is:

\[ \left[ q_{0.025,m}^*, q_{0.975,m}^* \right]. \]

In percentage points, this becomes:

\[ \left[ 100q_{0.025,m}^*, 100q_{0.975,m}^* \right]. \]

17.5 Bootstrap-Based Hypothesis Test

For hypothesis testing, the lecture tests:

\[ H_0: \Delta_m = 0 \]

against:

\[ H_1: \Delta_m \neq 0. \]

The test statistic is:

\[ z_m = \frac{\hat\Delta_m} {\widehat{SE}_{boot}(\hat\Delta_m)}. \]

The two-sided p-value is:

\[ p_m = 2 \left\{ 1 - \Phi\left(|z_m|\right) \right\}, \]

where \(\Phi(\cdot)\) is the cumulative distribution function of the standard normal distribution.

If \(p_m < 0.05\), we reject \(H_0\) at the 5% significance level.

In plain language:

If zero is far away from the estimate compared with the bootstrap standard error, we get stronger evidence against no adjusted risk difference.

18 References

Efron, B. (1979). Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1), 1-26.

Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall/CRC.