Lectures 1-6 built one causal analysis in stages.
| Lecture | Main idea | How this note uses it |
|---|---|---|
| Lecture 1 | Simple observed risk difference | Start with the crude association |
| Lecture 2 | Inverse probability weighting | Test an IPW-adjusted risk difference |
| Lecture 3 | Outcome regression standardization | Test a standardized risk difference |
| Lecture 4 | Doubly robust estimation | Test a doubly robust risk difference |
| Lecture 5 | Bootstrap confidence intervals | Use the bootstrap to estimate uncertainty |
| Lecture 6 | E-values and sensitivity analysis | Keep statistical evidence separate from causal-validity evidence |
The new question is:
Once we have a risk-difference estimate, how much statistical evidence is there against no difference?
The causal question stays the same:
\[ \text{Is CVD risk different under smoking versus non-smoking?} \]
The target causal risk difference is
\[ RD = E\{Y(1)\} - E\{Y(0)\}. \]
Here \(Y(1)\) is the potential CVD outcome under smoking and \(Y(0)\) is the potential CVD outcome under non-smoking.
We will again use the classroom NHANES complete-case dataset from Lectures 1-6 and we will again not use survey weights. The goal is to understand causal estimands, uncertainty, and hypothesis testing in the same teaching example, not to produce a full NHANES survey analysis.
Hypothesis testing does not create a causal estimand and it does not make confounding disappear.
It answers a narrower statistical question:
If the true target contrast were equal to a chosen null value, would our estimate look unusually far away from that value?
In this course, there are three layers:
| Layer | Question | Main lectures |
|---|---|---|
| Causal question | What would happen under smoking versus non-smoking? | Lecture 1 |
| Identification and estimation | How do we estimate \(E\{Y(1)\}\) and \(E\{Y(0)\}\)? | Lectures 2-4 |
| Statistical uncertainty | How variable is the estimate across samples? | Lecture 5 and this note |
| Causal-validity uncertainty | How sensitive is the answer to unmeasured confounding? | Lecture 6 |
So a p-value belongs to the statistical uncertainty layer. A p-value can be small even when the causal assumptions are wrong. A p-value can also be large even when an effect exists but the data are noisy.
For the risk difference, the most common point null hypothesis is
\[ H_0: RD = 0. \]
The two-sided alternative is
\[ H_A: RD \neq 0. \]
In words, the null says that the marginal CVD risk under smoking equals the marginal CVD risk under non-smoking:
\[ E\{Y(1)\} = E\{Y(0)\}. \]
The same idea can be written on other effect scales. For a risk ratio, the no-effect null is
\[ H_0: \frac{E\{Y(1)\}}{E\{Y(0)\}} = 1. \]
For a log risk ratio, the same null is
\[ H_0: \log\left( \frac{E\{Y(1)\}}{E\{Y(0)\}} \right) = 0. \]
This note focuses on the risk difference because Lectures 1-5 used risk differences as the main classroom contrast.
Suppose an estimator satisfies the large-sample approximation
\[ \hat\theta \approx N(\theta, SE^2). \]
Here \(\theta\) is a generic estimand. In our main analysis, \(\theta\) will be the risk difference \(RD\). To test \(H_0: \theta = \theta_0\), the Wald statistic is
\[ Z = \frac{\hat\theta - \theta_0}{\widehat{SE}(\hat\theta)}. \]
Under the null hypothesis and regularity conditions,
\[ Z \approx N(0,1). \]
The two-sided p-value is
\[ p = 2\{1 - \Phi(|Z_{\mathrm{obs}}|)\}. \]
The p-value is not \(P(H_0 \mid \text{data})\). It is the probability, computed under the null reference distribution, of obtaining a test statistic at least as extreme as the one observed.
For a two-sided 5 percent Wald test, the null value \(\theta_0\) is rejected when
\[ \left| \frac{\hat\theta - \theta_0} {\widehat{SE}(\hat\theta)} \right| > 1.96. \]
Equivalently, \(\theta_0\) is rejected if it lies outside the 95 percent Wald confidence interval
\[ \hat\theta \pm 1.96 \widehat{SE}(\hat\theta). \]
Bootstrap percentile intervals use a different construction, but the teaching interpretation is similar: if the 95 percent interval for the risk difference excludes zero, the data are not very compatible with \(RD=0\) under the assumptions behind the analysis.
Confidence intervals are often more informative than p-values because they show the direction and magnitude of the estimated effect.
We use the same classroom NHANES complete-case analytic sample as the previous lectures.
data_candidates <- c(
file.path("..", "Data", "NHANES_data.csv"),
file.path("Data", "NHANES_data.csv"),
"C:/Users/seyoo/OneDrive/GitHub_repository/Causal_Inference_Survival_Analysis/Causal_Inference_Lecture/Data/NHANES_data.csv"
)
data_path <- data_candidates[file.exists(data_candidates)][1]
if (is.na(data_path)) {
stop("Could not find NHANES_data.csv. Check the data path.")
}
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)
ht_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]
ht_data$smoker_indicator <- as.integer(ht_data$smoker_indicator)
ht_data$cvd_indicator <- as.integer(ht_data$cvd_indicator)
ht_data$smoking_group <- factor(
ifelse(ht_data$smoker_indicator == 1, "Smoker", "Non-smoker"),
levels = c("Non-smoker", "Smoker")
)
ht_data$gender_label <- factor(
ht_data$gender,
levels = c("F", "M"),
labels = c("Female", "Male")
)
ht_data$race_label <- factor(
ht_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"
)
)
ht_data$educ_label <- factor(
ht_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(ht_data)
## [1] 6299
sample_summary <- data.frame(
Group = c("Non-smoker", "Smoker", "Total"),
N = c(
sum(ht_data$smoker_indicator == 0),
sum(ht_data$smoker_indicator == 1),
nrow(ht_data)
),
CVD_percent = round(100 * c(
mean(ht_data$cvd_indicator[ht_data$smoker_indicator == 0]),
mean(ht_data$cvd_indicator[ht_data$smoker_indicator == 1]),
mean(ht_data$cvd_indicator)
), 2)
)
knitr::kable(sample_summary)
| Group | N | CVD_percent |
|---|---|---|
| Non-smoker | 3740 | 6.93 |
| Smoker | 2559 | 15.12 |
| Total | 6299 | 10.26 |
First we rebuild the point estimates from the earlier causal lectures.
t1 <- ht_data$smoker_indicator == 1
t0 <- ht_data$smoker_indicator == 0
n1 <- sum(t1)
n0 <- sum(t0)
p1 <- mean(ht_data$cvd_indicator[t1])
p0 <- mean(ht_data$cvd_indicator[t0])
crude_estimates <- c(
psi1 = p1,
psi0 = p0,
rd = p1 - p0,
rr = p1 / p0
)
adjusted_estimates <- fit_course_estimators(ht_data)
all_estimates <- rbind(
Observed_comparison = crude_estimates,
adjusted_estimates
)
estimate_table <- data.frame(
Lecture_connection = c(
"Lecture 1",
"Lecture 2",
"Lecture 3",
"Lecture 4"
),
Method = c(
"Observed comparison",
"IPW",
"Outcome regression",
"Doubly robust"
),
Non_smoker_risk_percent = round(100 * all_estimates[, "psi0"], 2),
Smoker_risk_percent = round(100 * all_estimates[, "psi1"], 2),
Risk_difference_pp = round(100 * all_estimates[, "rd"], 2),
Risk_ratio = round(all_estimates[, "rr"], 3),
row.names = NULL
)
knitr::kable(
estimate_table,
col.names = c(
"Connection",
"Method",
"Non-smoker risk (%)",
"Smoker risk (%)",
"Risk difference (pp)",
"Risk ratio"
)
)
| Connection | Method | Non-smoker risk (%) | Smoker risk (%) | Risk difference (pp) | Risk ratio |
|---|---|---|---|---|---|
| Lecture 1 | Observed comparison | 6.93 | 15.12 | 8.20 | 2.184 |
| Lecture 2 | IPW | 8.71 | 11.70 | 2.99 | 1.343 |
| Lecture 3 | Outcome regression | 8.53 | 11.98 | 3.45 | 1.405 |
| Lecture 4 | Doubly robust | 8.56 | 11.82 | 3.25 | 1.380 |
Here pp means percentage points.
These four rows answer related but not identical questions:
For the Lecture 1 observed comparison, the usual two-sample standard error for a difference in proportions is
\[ \widehat{SE}_{\mathrm{crude}} = \sqrt{ \frac{\hat p_1(1-\hat p_1)}{n_1} + \frac{\hat p_0(1-\hat p_0)}{n_0} }. \]
This is an association test, not a fully adjusted causal test.
rd_crude <- crude_estimates["rd"]
rr_crude <- crude_estimates["rr"]
se_crude <- sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)
z_crude <- rd_crude / se_crude
p_crude <- 2 * pnorm(-abs(z_crude))
ci_crude <- rd_crude + c(-1, 1) * qnorm(0.975) * se_crude
crude_table <- data.frame(
Quantity = c(
"Observed non-smoker CVD risk",
"Observed smoker CVD risk",
"Observed risk difference",
"Observed risk ratio",
"Standard error for risk difference",
"Wald Z statistic for H0: RD = 0",
"Two-sided p-value",
"95% CI lower",
"95% CI upper"
),
Value = c(
p0,
p1,
rd_crude,
rr_crude,
se_crude,
z_crude,
p_crude,
ci_crude[1],
ci_crude[2]
)
)
crude_table$Value <- round(crude_table$Value, 4)
knitr::kable(crude_table)
| Quantity | Value |
|---|---|
| Observed non-smoker CVD risk | 0.0693 |
| Observed smoker CVD risk | 0.1512 |
| Observed risk difference | 0.0820 |
| Observed risk ratio | 2.1838 |
| Standard error for risk difference | 0.0082 |
| Wald Z statistic for H0: RD = 0 | 9.9861 |
| Two-sided p-value | 0.0000 |
| 95% CI lower | 0.0659 |
| 95% CI upper | 0.0981 |
This test is useful because it is transparent. But it tests the raw smoker-versus-non-smoker comparison. It has a causal interpretation only under a very strong assumption: the smoker and non-smoker groups would have been comparable even without adjustment.
For IPW, outcome regression, and the doubly robust estimator, the standard error should reflect the full procedure:
That is why Lecture 5 used 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(20260706)
B <- 300
n <- nrow(ht_data)
bootstrap_rd <- matrix(
NA_real_,
nrow = B,
ncol = nrow(adjusted_estimates)
)
colnames(bootstrap_rd) <- rownames(adjusted_estimates)
for (b in seq_len(B)) {
bootstrap_rows <- sample.int(n, size = n, replace = TRUE)
bootstrap_data <- ht_data[bootstrap_rows, ]
bootstrap_rd[b, ] <- safe_fit_adjusted_rd(bootstrap_data)
}
bootstrap_rd <- bootstrap_rd[complete.cases(bootstrap_rd), , drop = FALSE]
number_of_successful_bootstraps <- nrow(bootstrap_rd)
number_of_successful_bootstraps
## [1] 300
bootstrap_se <- apply(bootstrap_rd, 2, sd)
bootstrap_ci <- t(apply(
bootstrap_rd,
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
colnames(bootstrap_ci) <- c("Lower", "Upper")
adjusted_rd <- adjusted_estimates[, "rd"]
z_adjusted <- adjusted_rd / bootstrap_se[rownames(adjusted_estimates)]
p_adjusted <- 2 * pnorm(-abs(z_adjusted))
adjusted_inference_table <- data.frame(
Method = c("IPW", "Outcome regression", "Doubly robust"),
Estimate_pp = round(100 * adjusted_rd, 2),
Bootstrap_SE_pp = round(
100 * bootstrap_se[rownames(adjusted_estimates)],
2
),
Z_statistic = round(z_adjusted, 2),
P_value = signif(p_adjusted, 3),
CI_lower_pp = round(100 * bootstrap_ci[, "Lower"], 2),
CI_upper_pp = round(100 * bootstrap_ci[, "Upper"], 2),
Conclusion_at_0_05 = ifelse(
p_adjusted < 0.05,
"Reject H0",
"Do not reject H0"
),
row.names = NULL
)
knitr::kable(
adjusted_inference_table,
col.names = c(
"Method",
"Estimate (pp)",
"Bootstrap SE (pp)",
"z",
"p-value",
"95% CI lower (pp)",
"95% CI upper (pp)",
"Conclusion at 0.05"
)
)
| Method | Estimate (pp) | Bootstrap SE (pp) | z | p-value | 95% CI lower (pp) | 95% CI upper (pp) | Conclusion at 0.05 |
|---|---|---|---|---|---|---|---|
| IPW | 2.99 | 0.79 | 3.78 | 1.60e-04 | 1.36 | 4.55 | Reject H0 |
| Outcome regression | 3.45 | 0.78 | 4.43 | 9.50e-06 | 1.84 | 4.98 | Reject H0 |
| Doubly robust | 3.25 | 0.78 | 4.17 | 3.02e-05 | 1.69 | 4.78 | Reject H0 |
The bootstrap standard error measures how much the adjusted estimate moves when we repeatedly rebuild the full estimator in resampled datasets.
Now we place the crude analytic test and the three adjusted bootstrap tests in one table.
combined_inference <- rbind(
data.frame(
Method = "Observed comparison",
Estimate_pp = 100 * rd_crude,
SE_pp = 100 * se_crude,
Z_statistic = z_crude,
P_value = p_crude,
CI_lower_pp = 100 * ci_crude[1],
CI_upper_pp = 100 * ci_crude[2],
SE_source = "Two-sample proportion formula"
),
data.frame(
Method = c("IPW", "Outcome regression", "Doubly robust"),
Estimate_pp = 100 * adjusted_rd,
SE_pp = 100 * bootstrap_se[rownames(adjusted_estimates)],
Z_statistic = z_adjusted,
P_value = p_adjusted,
CI_lower_pp = 100 * bootstrap_ci[, "Lower"],
CI_upper_pp = 100 * bootstrap_ci[, "Upper"],
SE_source = "Nonparametric bootstrap"
)
)
combined_display <- combined_inference
combined_display$Estimate_pp <- round(combined_display$Estimate_pp, 2)
combined_display$SE_pp <- round(combined_display$SE_pp, 2)
combined_display$Z_statistic <- round(combined_display$Z_statistic, 2)
combined_display$P_value <- signif(combined_display$P_value, 3)
combined_display$CI_lower_pp <- round(combined_display$CI_lower_pp, 2)
combined_display$CI_upper_pp <- round(combined_display$CI_upper_pp, 2)
knitr::kable(
combined_display,
col.names = c(
"Method",
"Estimate (pp)",
"SE (pp)",
"z",
"p-value",
"95% CI lower (pp)",
"95% CI upper (pp)",
"SE source"
)
)
| Method | Estimate (pp) | SE (pp) | z | p-value | 95% CI lower (pp) | 95% CI upper (pp) | SE source | |
|---|---|---|---|---|---|---|---|---|
| rd | Observed comparison | 8.20 | 0.82 | 9.99 | 0.00e+00 | 6.59 | 9.81 | Two-sample proportion formula |
| IPW | IPW | 2.99 | 0.79 | 3.78 | 1.60e-04 | 1.36 | 4.55 | Nonparametric bootstrap |
| Outcome_regression | Outcome regression | 3.45 | 0.78 | 4.43 | 9.50e-06 | 1.84 | 4.98 | Nonparametric bootstrap |
| Doubly_robust | Doubly robust | 3.25 | 0.78 | 4.17 | 3.02e-05 | 1.69 | 4.78 | Nonparametric bootstrap |
The vertical dashed line is \(RD=0\), the null value. Estimates farther from zero relative to their standard errors produce larger absolute \(z\)-statistics and smaller p-values.
The crude test and the adjusted tests can lead to different numbers because they estimate different contrasts.
The crude comparison tests:
\[ P(Y=1 \mid T=1) - P(Y=1 \mid T=0) = 0. \]
The adjusted analyses test an identified causal contrast under the assumptions from Lectures 2-4:
\[ E\{Y(1)\} - E\{Y(0)\} = 0. \]
Those two are the same only if the crude comparison is already unconfounded. In observational smoking data, that is usually not a comfortable assumption.
A small adjusted p-value means:
Under the fitted adjustment strategy and its assumptions, the estimated adjusted risk difference is far from zero relative to its estimated sampling uncertainty.
It does not mean:
Conversely, a large p-value does not prove no effect. It may reflect a small effect, a noisy estimator, limited sample size, poor overlap, or an estimand whose standard error is large.
Lecture 6 asked a different question:
How strong would unmeasured confounding have to be to explain away the result?
That is not the same as a p-value.
| Tool | Main question | What it does not solve |
|---|---|---|
| p-value | Is the estimate far from the null relative to sampling uncertainty? | Unmeasured confounding |
| Confidence interval | What range of effect sizes is compatible with sampling uncertainty? | Bias from wrong assumptions |
| E-value | How strong would unmeasured confounding need to be on the risk-ratio scale? | Random sampling variability by itself |
So the course message is:
Report statistical uncertainty and causal-validity concerns side by side.
The best causal report does not stop at “p < 0.05.” It states the estimand, method, confidence interval, p-value, and key assumptions.
If the scientific question is directional, such as whether smoking increases CVD risk, the alternative can be written
\[ H_A: RD > 0. \]
The one-sided p-value is
\[ p_{\mathrm{one-sided}} = P\{N(0,1) \geq Z_{\mathrm{obs}}\} = 1 - \Phi(Z_{\mathrm{obs}}) \]
when large positive values are evidence against the null.
A one-sided test should be chosen before seeing the data. Switching from two-sided to one-sided after observing the estimate is not a valid way to reduce a p-value.
A concise classroom report should include:
For example, a doubly robust result can be reported as follows.
dr_report <- combined_inference[combined_inference$Method == "Doubly robust", ]
cat(
"Using the doubly robust estimator, the adjusted risk difference was ",
round(dr_report$Estimate_pp, 2),
" percentage points, with a 95% interval from ",
round(dr_report$CI_lower_pp, 2),
" to ",
round(dr_report$CI_upper_pp, 2),
" percentage points and a two-sided p-value of ",
signif(dr_report$P_value, 3),
".",
sep = ""
)
## Using the doubly robust estimator, the adjusted risk difference was 3.25 percentage points, with a 95% interval from 1.69 to 4.78 percentage points and a two-sided p-value of 3.02e-05.
This sentence is only complete if it is followed by the causal assumptions and the sensitivity-analysis caveats.
Before interpreting a hypothesis test in this course, check the following:
This appendix collects the formulas behind the calculations in the main text.
The observed risks are
\[ \hat p_1 = \frac{\sum_{i=1}^{n} T_iY_i} {\sum_{i=1}^{n} T_i} \]
and
\[ \hat p_0 = \frac{\sum_{i=1}^{n} (1-T_i)Y_i} {\sum_{i=1}^{n} (1-T_i)}. \]
The observed risk difference is
\[ \widehat{RD}_{obs} = \hat p_1 - \hat p_0. \]
The standard error used in the Lecture 1-style test is
\[ \widehat{SE}(\widehat{RD}_{obs}) = \sqrt{ \frac{\hat p_1(1-\hat p_1)}{n_1} + \frac{\hat p_0(1-\hat p_0)}{n_0} }. \]
Let
\[ e(X_i) = P(T_i = 1 \mid X_i) \]
be the propensity score.
The IPW risk under smoking is
\[ \hat\psi_1^{IPW} = \frac{ \sum_{i=1}^{n} \dfrac{T_iY_i}{\hat e(X_i)} }{ \sum_{i=1}^{n} \dfrac{T_i}{\hat e(X_i)} }. \]
The IPW risk under non-smoking is
\[ \hat\psi_0^{IPW} = \frac{ \sum_{i=1}^{n} \dfrac{(1-T_i)Y_i}{1-\hat e(X_i)} }{ \sum_{i=1}^{n} \dfrac{1-T_i}{1-\hat e(X_i)} }. \]
Then
\[ \widehat{RD}^{IPW} = \hat\psi_1^{IPW} - \hat\psi_0^{IPW}. \]
Let
\[ m(t, X_i) = E(Y_i \mid T_i=t, X_i). \]
After fitting the outcome model, standardization estimates
\[ \hat\psi_t^{OR} = \frac{1}{n} \sum_{i=1}^{n} \hat m(t, X_i), \qquad t \in \{0,1\}. \]
The outcome-regression risk difference is
\[ \widehat{RD}^{OR} = \hat\psi_1^{OR} - \hat\psi_0^{OR}. \]
The doubly robust risk under smoking is
\[ \hat\psi_1^{DR} = \frac{1}{n} \sum_{i=1}^{n} \left[ \hat m(1,X_i) + \frac{T_i}{\hat e(X_i)} \{Y_i - \hat m(1,X_i)\} \right]. \]
The doubly robust risk under non-smoking is
\[ \hat\psi_0^{DR} = \frac{1}{n} \sum_{i=1}^{n} \left[ \hat m(0,X_i) + \frac{1-T_i}{1-\hat e(X_i)} \{Y_i - \hat m(0,X_i)\} \right]. \]
The doubly robust risk difference is
\[ \widehat{RD}^{DR} = \hat\psi_1^{DR} - \hat\psi_0^{DR}. \]
Let \(\widehat{RD}^{*(1)}, \ldots, \widehat{RD}^{*(B)}\) be the bootstrap estimates.
The bootstrap standard error is
\[ \widehat{SE}_{boot} = \sqrt{ \frac{1}{B-1} \sum_{b=1}^{B} \left( \widehat{RD}^{*(b)} - \overline{RD}^{*} \right)^2 }, \]
where
\[ \overline{RD}^{*} = \frac{1}{B} \sum_{b=1}^{B} \widehat{RD}^{*(b)}. \]
The Wald statistic using the bootstrap standard error is
\[ Z_{boot} = \frac{\widehat{RD}}{\widehat{SE}_{boot}} \]
for testing \(H_0: RD=0\).
The percentile bootstrap 95 percent confidence interval is
\[ \left[ Q_{0.025}\{\widehat{RD}^{*}\}, Q_{0.975}\{\widehat{RD}^{*}\} \right], \]
where \(Q_a\{\widehat{RD}^{*}\}\) is the \(a\)-quantile of the bootstrap estimates.
Austin, P. C., & Stuart, E. A. (2015). Moving towards best practice when using inverse probability of treatment weighting using the propensity score to estimate causal treatment effects in observational studies. Statistics in Medicine, 34(28), 3661-3679.
Cole, S. R., & Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology, 168(6), 656-664.
Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall/CRC.
Hernan, M. A., & Robins, J. M. (2020). Causal Inference: What If. Chapman & Hall/CRC.
Imbens, G. W., & Rubin, D. B. (2015). Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction. Cambridge University Press.
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.