V The V Lab

Causal Inference in Depth

From Target Trials, Causal Diagrams, and Potential Outcomes to Standardization, Propensity Scores, and Sensitivity Analysis

The V Lab Home Biostatistics Foundations Medical Statistical Tests Linear Regression Logistic Regression Multilevel Modelling Causal Inference Propensity Score Matching Survival Analysis AFT and Cox PH 中文版

AudienceLearners in public health, epidemiology, medicine, and data science

Study timeApproximately 180–240 minutes

PrerequisitesRegression, confidence intervals, and basic R syntax

About the data and findings in this tutorial All individuals, interventions, counterfactual outcomes, and numerical findings were simulated with fixed random seeds solely to demonstrate the methods. In real data, only one potential outcome is observed for each person, and the true causal effect cannot be looked up directly. This tutorial retains the “truth” only so that we can test whether each method recovers a known answer.

How to use this tutorial

This tutorial centers on one question: If everyone in the target population participated in an intensive blood pressure management program, rather than everyone receiving usual care, how much would mean systolic blood pressure differ at 6 months?

The recommended sequence is: “define the question → emulate the target trial → draw the causal diagram → state the identification assumptions → choose an estimation method → examine diagnostics → conduct sensitivity analyses → report transparently.” Code is shown by default; you can run it section by section or collapse it with the page controls.

Learning objectives

After completing this tutorial, you should be able to:

• Distinguish descriptive, predictive, and causal questions, and define a causal question precisely using the target trial framework;
• Define the ATE, ATT, CATE, risk difference, risk ratio, and other causal estimands using potential outcomes;
• Explain consistency, exchangeability, positivity, and no interference;
• Use directed acyclic graphs to identify confounders, mediators, colliders, and variables that should not be adjusted for;
• Understand randomization, intention-to-treat analysis, adherence, and loss to follow-up in randomized trials;
• Estimate causal effects using standardization (g-computation), propensity score weighting, matching, and augmented inverse probability weighting;
• Assess covariate balance, propensity score overlap, extreme weights, model form, and changes in the target population;
• Distinguish total effects, direct effects, effect modification, and exploratory subgroup analyses;
• Describe the core assumptions underlying difference-in-differences, regression discontinuity, and instrumental variable designs;
• Design sensitivity analyses for unmeasured confounding, missingness, selection, and modeling choices;
• Write an auditable causal research report that does not overstate its conclusions.

1 Causal Questions Begin with “What If?”

1.1 Association, prediction, and causation are different tasks

The same data and regression function can serve different goals, but they answer different questions:

Task	Typical question	Primary basis for evaluation
Description	How much does mean blood pressure differ between program participants and nonparticipants?	Accuracy of the sample, measurements, and description
Prediction	Who is likely to still have high blood pressure after 6 months?	Out-of-sample calibration, discrimination, and error
Causation	What would happen to outcomes if the same target population participated in the program rather than receiving usual care?	Design, temporal ordering, identification assumptions, and estimation

An “adjusted regression coefficient” may still be only a conditional association. A causal interpretation is not granted by lm(), glm(), or a statistically significant p-value. It comes from a clearly defined intervention contrast, a credible design, and assumptions sufficient to connect observed data to counterfactual outcomes.

1.2 Start by emulating the ideal target trial

A target trial is not necessarily a trial that will actually be conducted. It is a protocol describing how an ideal randomized trial would answer the question if ethics, time, and resources allowed. An observational study can attempt to emulate this protocol.

Protocol element	Target trial in this tutorial
Eligibility criteria	Members of the target population who begin management at baseline and meet prespecified clinical criteria
Treatment strategies	Enroll immediately in the intensive program or continue usual care
Assignment procedure	Random assignment in the ideal trial; adjustment for recorded factors in the observational cohort
Time zero	The common time at which eligibility is confirmed and the treatment strategy is assigned
Follow-up	From time zero through 6 months
Outcome	Systolic blood pressure at 6 months (mmHg)
Causal contrast	Mean difference if everyone participated in the program versus if everyone received usual care
Analysis principle	First estimate the population-average effect of assignment to a strategy

The target trial makes visible many problems that software can hide. Were periods before and after treatment initiation mixed together? Did people who had to “survive until treatment” receive immortal time? Did the eligibility criteria use future information? Did follow-up begin at the same time in each group?

Time zero must be aligned If eligibility confirmation, treatment assignment, and the start of follow-up occur at different times, selection bias and immortal time bias may arise before any model is fit. Adding covariates usually cannot repair this design misalignment.

1.3 Define treatment versions and comparison strategies

“Receive better care,” “exercise more,” and “control blood pressure” are not sufficiently well-defined interventions. An interpretable causal question must specify, at a minimum, who receives what and when, for how long, which co-interventions are allowed, and what the comparison strategy is.

Different versions of the “program” may involve different follow-up schedules, authority to adjust medications, and adherence support. If these versions affect outcomes differently but are combined into one binary variable, the consistency assumption becomes ambiguous. Whenever possible, write each strategy as a rule that another person could implement or audit.

2 Potential Outcomes and Causal Estimands

2.1 Each person has two potential outcomes

Let $A=1$ denote participation in the program and $A=0$ denote usual care. For individual $i$:

$$Y_i(1)=\text{the individual’s 6-month blood pressure under the program},\qquad Y_i(0)=\text{the individual’s 6-month blood pressure under usual care}.$$

The individual causal effect is $Y_i(1)-Y_i(0)$. In reality, only one of these outcomes can be observed:

$$Y_i=A_iY_i(1)+(1-A_i)Y_i(0).$$

The unobserved potential outcome is not an ordinary missing value; it cannot be recovered by measuring the same person again at the same time. This is the fundamental problem of causal inference. The task of research design and statistical methods is to construct credible counterfactual comparisons at the population level.

2.2 ATE, ATT, and CATE answer different questions

The most common estimand is the average treatment effect (ATE):

$$ATE=E\{Y(1)-Y(0)\}.$$

This tutorial uses the direction “program minus usual care,” so a negative value means that the program lowers blood pressure. Two other common estimands are:

Estimand	Definition	Target population
ATE	$E\{Y(1)-Y(0)\}$	Entire study target population
ATT	$E\{Y(1)-Y(0)\mid A=1\}$	People who actually participated in the program
CATE	$E\{Y(1)-Y(0)\mid X=x\}$	Subgroup with characteristics $X=x$

If the treatment effect varies with baseline blood pressure and program participants more often have high baseline blood pressure, the ATE and ATT may differ. The estimand must be chosen before selecting a weighting or matching method; the target population should not be changed after seeing which result is more statistically significant.

Because the simulation retains both otherwise unobservable potential outcomes, we know the finite-sample truth for this realized simulated cohort: the ATE is -5.06 mmHg, and the ATT among actual program participants is -5.16 mmHg. Subsequent analyses provide the methods only with the ordinarily observable causal_data.

2.3 The effect scale must match the decision

For continuous outcomes, possible scales include mean differences, ratios, and quantile differences. Common scales for binary outcomes include:

$$RD=P\{Y(1)=1\}-P\{Y(0)=1\}, \quad RR=\frac{P\{Y(1)=1\}}{P\{Y(0)=1\}}, \quad OR=\frac{P\{Y(1)=1\}/P\{Y(1)=0\}} {P\{Y(0)=1\}/P\{Y(0)=0\}}.$$

The risk difference directly expresses how many fewer or additional events occur per 100 people, the risk ratio expresses relative change, and the odds ratio generally does not equal the risk ratio. Effect modification also depends on the scale: a variable may interact with treatment on the risk-difference scale but not on the risk-ratio scale.

Write the estimand in one sentence first Recommended template: In [target population], compare [strategy 1] with [strategy 0] starting at [time zero], with respect to the [marginal effect scale] for [outcome] through [time horizon].

3 From Counterfactuals to Observed Data: Identification Assumptions

3.1 Consistency

Consistency requires that a person who actually receives $A=a$ has an observed outcome equal to the corresponding potential outcome: if $A_i=a$, then $Y_i=Y_i(a)$. It also implies that treatment versions are defined with sufficient precision.

Potential threats to consistency include:

• The “program” consists of entirely different services across regions;
• Electronic health records document only a “referral,” with no indication of whether the intervention was actually received;
• Exposure dose, initiation time, or co-interventions are undefined;
• A self-reported exposure is not the same behavior one intends to intervene on.

3.2 Conditional exchangeability: no uncontrolled common causes

Randomized trials use randomization to make $A\perp\{Y(1),Y(0)\}$. An observational study can usually claim only that, conditional on a set of pretreatment common causes $L$,

$$\{Y(1),Y(0)\}\perp A\mid L.$$

This is often called “no unmeasured confounding.” It cannot be established from the data alone; it requires subject-matter knowledge, measurement quality, temporal ordering, and sensitivity analysis. Automatically placing many variables in a model does not guarantee that all important common causes were measured correctly.

3.3 Positivity: every type of person needs comparable strategies

For every $L=l$ that occurs with positive probability in the target population, we require:

$$0<P(A=1\mid L=l)<1.$$

A structural positivity violation occurs when a group of people cannot receive a strategy by design—for example, someone with an absolute contraindication cannot receive a medication. Changing the model cannot create the missing counterfactual information. Practical positivity problems arise when, in a finite sample, some covariate combinations almost always receive only one strategy. They appear as propensity scores near 0 or 1, extreme weights, and unstable estimates.

3.4 No interference and treatment definition

The usual framework also assumes that one person’s outcome is unaffected by other people’s treatment, an assumption known as no interference. Spillover effects are common, however, with vaccines, infectious disease control, peer interventions, and hospital workflows. In those settings, group coverage, network exposure, or cluster-level strategies must be defined; an individual binary treatment no longer describes the question adequately.

3.5 Identification formula

Under consistency, conditional exchangeability, and positivity, the mean potential outcome can be identified by the g-formula:

$$E\{Y(a)\}=E_L[E(Y\mid A=a,L)].$$

The right-hand side contains only observable distributions: compare treatment strategies within each value of $L$, then average over the distribution of $L$ in the target population. Standardization, outcome regression, stratification, and many machine-learning g-computation methods all implement this logic.

4 Causal Diagrams: Deciding What to Adjust For

4.1 Nodes, arrows, and paths in a DAG

A directed acyclic graph (DAG) uses arrows to represent assumed direct causal relationships. A DAG is not discovered automatically from the data; it is an explicit statement of the researcher’s understanding of the data-generating process.

draw_node <- function(x, y, label, fill = "white", width = 0.16) {
  rect(x - width / 2, y - 0.055, x + width / 2, y + 0.055,
       col = fill, border = palette_ci["navy"], lwd = 1.7)
  text(x, y, label, cex = 0.88)
}

draw_arrow <- function(x0, y0, x1, y1, color = palette_ci["gray"]) {
  arrows(x0, y0, x1, y1, length = 0.08, lwd = 1.8, col = color)
}

draw_path_arrow <- function(x, y, color = palette_ci["gray"]) {
  lines(x[-length(x)], y[-length(y)], lwd = 1.8, col = color)
  arrows(
    x[length(x) - 1], y[length(y) - 1], x[length(x)], y[length(y)],
    length = 0.08, lwd = 1.8, col = color
  )
}

plot.new()
plot.window(xlim = c(0, 1), ylim = c(0, 1))

draw_node(0.13, 0.75, "L\nBaseline common causes", "#E8F2F1", 0.22)
draw_node(0.39, 0.75, "A\nProgram", "#E7F1FA")
draw_node(0.64, 0.75, "M\nEngagement", "#FFF3D8")
draw_node(0.87, 0.75, "Y\n6-month SBP", "#FCE8E4")
draw_arrow(0.24, 0.75, 0.30, 0.75)
draw_arrow(0.47, 0.75, 0.55, 0.75)
draw_arrow(0.72, 0.75, 0.79, 0.75)
draw_path_arrow(c(0.13, 0.13, 0.87, 0.87),
                c(0.81, 0.89, 0.89, 0.81))
draw_path_arrow(c(0.39, 0.39, 0.87, 0.87),
                c(0.69, 0.61, 0.61, 0.69), palette_ci["blue"])

draw_node(0.13, 0.27, "A\nProgram", "#E7F1FA")
draw_node(0.40, 0.27, "S\nPost-treatment visit", "#F5E7F2", 0.22)
draw_node(0.67, 0.27, "U\nUnmeasured health need", "#EEEEEE", 0.20)
draw_node(0.87, 0.27, "Y\n6-month SBP", "#FCE8E4")
draw_arrow(0.21, 0.27, 0.31, 0.27)
draw_arrow(0.59, 0.27, 0.49, 0.27)
draw_arrow(0.75, 0.27, 0.79, 0.27)
text(
  0.50, 0.12,
  "S is the collider A -> S <- U; conditioning on S opens a noncausal path",
  cex = 0.78
)

box()

Simplified causal diagrams of a pretreatment common cause, a mediator, and a collider.

In the first diagram, $L\rightarrow A$ and $L\rightarrow Y$, so $A\leftarrow L\rightarrow Y$ is a backdoor path that must be blocked through design or analysis. $M$ lies on $A\rightarrow M\rightarrow Y$ and represents part of the program’s effect; it generally should not be adjusted for when estimating the total effect.

In the second diagram, the post-treatment visit $S$ is affected by both the program $A$ and unmeasured health need $U$. Restricting the sample to “people with at least one visit” or including $S$ in the model induces a conditional association between $A$ and $U$, thereby opening $A\rightarrow S\leftarrow U\rightarrow Y$.

4.2 The backdoor criterion and minimally sufficient adjustment sets

An adjustment set must block every backdoor path from $A$ to $Y$. When estimating a total effect, a safe introductory rule is to exclude post-treatment mediators and colliders; more general adjustment criteria must be evaluated for each graph. A minimally sufficient set is usually preferable to “every available variable.” A reasonable set for the primary analysis in this tutorial includes the pretreatment variables age, baseline blood pressure, smoking, rural residence, and health literacy.

Variable role	Structure in the graph	Usual treatment when estimating the total effect
Confounder	$A\leftarrow L\rightarrow Y$	Adjust to block the backdoor path
Pure outcome predictor	$A\quad L\rightarrow Y$	May improve precision but is not required to remove confounding
Instrumental variable	$Z\rightarrow A\rightarrow Y$	Usually does not need adjustment in an ordinary outcome regression
Mediator	$A\rightarrow M\rightarrow Y$	Do not adjust when estimating the total effect
Collider	$A\rightarrow S\leftarrow U\rightarrow Y$	Do not condition, stratify, or select the sample on it
Exposure or outcome proxy	Measurement structure depends on the specific process	Cannot be handled based on association alone

“Measured before treatment” is not a sufficient reason A variable occurring before treatment is not automatically a confounder. Common causes, instrumental variables, ancestors of colliders, and variables that affect only precision play different roles. Draw the graph based on time and subject-matter mechanisms before deciding how to use each variable.

4.3 What happens when you adjust for mediators and colliders?

In the simulated data, engagement_score is a post-treatment mediator and clinic_contact is a collider. The following comparison shows the program coefficient after reasonable pretreatment adjustment, after adding the mediator, and after adding the collider. These coefficients are not formal causal direct effects in every setting; the purpose is to show how “adjusting for more variables” can change the question or introduce bias.

model_pre_treatment <- lm(
  six_month_sbp ~ program_num + age_c10 + baseline_sbp_c10 +
    smoking_num + rural_num + health_literacy,
  data = causal_data
)

model_with_mediator <- update(
  model_pre_treatment,
  . ~ . + engagement_score
)

model_with_collider <- update(
  model_pre_treatment,
  . ~ . + clinic_contact_num
)

bad_control_table <- data.frame(
  Model = c(
    "Pretreatment common causes only",
    "Add the post-treatment mediator",
    "Add the post-treatment collider"
  ),
  `Program coefficient` = c(
    coef(model_pre_treatment)["program_num"],
    coef(model_with_mediator)["program_num"],
    coef(model_with_collider)["program_num"]
  ),
  `Question answered` = c(
    "Adjusted total-effect projection; not generally the population ATE",
    "Blocks part of the causal pathway; the question changes",
    "May open a noncausal path"
  ),
  check.names = FALSE
)

knitr::kable(
  bad_control_table,
  digits = 2,
  caption = "Program coefficients under different adjustment strategies"
)

Program coefficients under different adjustment strategies

Model	Program coefficient	Question answered
Pretreatment common causes only	-5.06	Adjusted total-effect projection; not generally the population ATE
Add the post-treatment mediator	-2.65	Blocks part of the causal pathway; the question changes
Add the post-treatment collider	-5.30	May open a noncausal path

The finite-sample true ATE is approximately -5.06 mmHg. Because the pretreatment model omits the known treatment-by-baseline-SBP interaction, its program coefficient is an adjusted model projection rather than the population ATE; it is shown only as a baseline for the bad-control comparison. After engagement is added, the program coefficient primarily retains the portion of the effect that does not operate through this mediator. Adding the visit variable may induce an artificial association between the program and unmeasured health need. Calling either quantity a “direct effect” would additionally require a well-defined intervention on the mediator, further identification assumptions, and appropriate handling of treatment–mediator interaction.

5 Randomized Trials: Design Comes Before Adjustment

5.1 Randomization creates exchangeability

If the program is randomly assigned with a fixed probability in the same target population, pretreatment characteristics will not systematically determine assignment on average. Randomization does not guarantee perfect balance in every finite sample, but it provides the design basis for quantifying random error and interpreting the contrast causally.

Below, the same simulated participants are reassigned at random. To avoid looking at unobservable counterfactuals, the analysis uses only the outcome that would be observed under each participant’s randomized assignment.

set.seed(20260811)
random_program <- rbinom(n, 1, 0.5)
random_outcome <- ifelse(
  random_program == 1,
  causal_truth$y1,
  causal_truth$y0
)

trial_difference <- with(
  data.frame(random_program, random_outcome),
  mean(random_outcome[random_program == 1]) -
    mean(random_outcome[random_program == 0])
)

trial_fit <- lm(random_outcome ~ random_program)
trial_ci <- confint(trial_fit)["random_program", ]

trial_result <- data.frame(
  Measure = c(
    "Randomized-trial mean difference",
    "Lower 95% CI",
    "Upper 95% CI",
    "Finite-sample simulated ATE benchmark"
  ),
  mmHg = c(trial_difference, trial_ci[1], trial_ci[2], true_ate)
)

knitr::kable(
  trial_result,
  digits = 2,
  caption = "Intention-to-treat contrast in one simulated randomized trial"
)

Intention-to-treat contrast in one simulated randomized trial

Measure	mmHg
Randomized-trial mean difference	-5.61
Lower 95% CI	-6.51
Upper 95% CI	-4.71
Finite-sample simulated ATE benchmark	-5.06

An estimate from one randomized trial need not equal the truth exactly; it is affected by random assignment and sampling variation. Adding strong pretreatment outcome predictors can improve precision, but randomization itself is the primary source of exchangeability. The interval shown here is an ordinary linear-model approximation. A formal trial should use design-consistent inference appropriate to its randomization scheme, heteroskedasticity, stratification, or cluster design.

5.2 Intention-to-treat, adherence, and per-protocol effects

An intention-to-treat (ITT) analysis compares participants according to their original randomized assignment, regardless of actual adherence. It preserves randomization and usually answers the question, “What is the effect of implementing the assignment strategy?”

A simple comparison by treatment actually received breaks randomization because adherence may be influenced by health status, preferences, and resources. If the target is a per-protocol effect under full adherence, deviations from the strategy must be defined clearly and appropriate methods must be used—for example, weighting for baseline and time-varying covariates, the g-formula, structural nested models, or, under additional assumptions, randomized assignment as an instrumental variable.

5.3 Loss to follow-up, missingness, and lack of blinding can still cause bias

Randomization does not guarantee that:

• Outcomes are not selectively missing because of loss to follow-up;
• The behavior of participants, clinicians, or outcome assessors is unaffected by knowledge of assignment;
• Outcome definitions and measurements are identical in the two groups;
• The randomized sample represents the target population to which the findings are intended to generalize;
• The treatment version used in the trial is the same as the strategy implemented in practice.

Randomized trials must therefore still report allocation concealment, blinding, adherence, between-group co-interventions, reasons for loss to follow-up, missing outcomes, and the limits of generalizability.

6 Observational Data: Start With the Unadjusted Comparison

6.1 Inspecting the analysis data

stopifnot(
  nrow(causal_data) == nrow(causal_truth),
  all(causal_data$program_num %in% c(0, 1)),
  !anyNA(causal_data),
  all(causal_data$baseline_sbp > 0),
  all(causal_data$six_month_sbp > 0)
)

preview <- head(causal_data[, c(
  "participant_id", "age", "rural", "smoking", "baseline_sbp",
  "health_literacy", "program", "six_month_sbp"
)])

knitr::kable(
  preview,
  col.names = c(
    "Participant", "Age", "Residence", "Smoking", "Baseline SBP",
    "Health literacy", "Management strategy", "6-month SBP"
  ),
  digits = 1,
  caption = "First six rows of the simulated observational cohort"
)

First six rows of the simulated observational cohort

Participant	Age	Residence	Smoking	Baseline SBP	Health literacy	Management strategy	6-month SBP
C0001	36	Rural	Not current	142	-1.9	Usual care	132
C0002	42	Urban	Not current	149	0.1	Usual care	135
C0003	65	Rural	Not current	150	0.4	Usual care	148
C0004	42	Urban	Not current	128	0.1	Usual care	116
C0005	48	Rural	Not current	150	-1.6	Usual care	147
C0006	69	Urban	Not current	137	3.2	Program	116

The unit of analysis is the participant. Pretreatment variables are measured before the program begins, and the outcome is measured at 6 months. The mediator and post-treatment clinic contact remain in the data as cautionary examples, but they are not included in the confounding adjustment set for the primary analysis.

6.2 Why the crude comparison is biased

crude_summary <- aggregate(
  cbind(baseline_sbp, six_month_sbp) ~ program,
  data = causal_data,
  FUN = mean
)

crude_difference <- with(
  causal_data,
  mean(six_month_sbp[program_num == 1]) -
    mean(six_month_sbp[program_num == 0])
)

knitr::kable(
  crude_summary,
  col.names = c("Management strategy", "Mean baseline SBP", "Mean 6-month SBP"),
  digits = 1,
  caption = "Crude description by observed program participation"
)

Crude description by observed program participation

Management strategy	Mean baseline SBP	Mean 6-month SBP
Usual care	145	136
Program	150	136

The unadjusted mean outcome difference is -0.65 mmHg, whereas the true ATE is -5.06 mmHg. Participants already had higher baseline blood pressure before entering the program, so a higher outcome or a smaller apparent decrease in the program group does not by itself show that the program was ineffective. This is a classic example of confounding by indication: people with greater need for intervention are more likely to receive it.

6.3 Using standardized differences to describe baseline imbalance

The standardized mean difference (SMD) does not increase mechanically with sample size, making it useful for describing baseline differences between groups:

$$SMD=\frac{\bar L_1-\bar L_0}{\sqrt{(s_1^2+s_0^2)/2}}.$$

It is not a test for confounding, and no threshold is definitive. An absolute SMD below 0.10 is often used as a rough diagnostic reference, but distributional shape, extreme values, important interactions, and nonlinear terms should also be examined.

balance_variables <- list(
  "Age" = causal_data$age,
  "Baseline SBP" = causal_data$baseline_sbp,
  "Current smoking" = causal_data$smoking_num,
  "Rural residence" = causal_data$rural_num,
  "Health literacy" = causal_data$health_literacy
)

unadjusted_smd <- vapply(
  balance_variables,
  standardized_difference,
  numeric(1),
  z = causal_data$program_num
)

balance_unadjusted <- data.frame(
  Covariate = names(unadjusted_smd),
  SMD = unname(unadjusted_smd),
  Absolute_SMD = abs(unname(unadjusted_smd)),
  check.names = FALSE
)

knitr::kable(
  balance_unadjusted,
  col.names = c("Covariate", "SMD", "Absolute SMD"),
  digits = 2,
  caption = "Unadjusted balance in pretreatment covariates"
)

Unadjusted balance in pretreatment covariates

Covariate	SMD	Absolute SMD
Age	0.46	0.46
Baseline SBP	0.64	0.64
Current smoking	0.32	0.32
Rural residence	0.00	0.00
Health literacy	0.08	0.08

Do not use p-values for baseline variables to select adjustment covariates. In large samples, trivial differences can be statistically significant; in small samples, important imbalances may not be. What matters more is each variable’s role in the causal structure and the practical magnitude of the imbalance.

7 Standardization and the Parametric g-Formula

7.1 Constructing two counterfactual worlds from an outcome model

Standardization first estimates $E(Y\mid A,L)$, then places every member of the target population under both $A=1$ and $A=0$, and finally averages the individual predictions. It produces a marginal effect in the target population, rather than a regression coefficient for a single “average” person.

This example allows the program effect to vary linearly with baseline blood pressure because both the data-generating mechanism and the scientific question support this heterogeneity.

outcome_model <- lm(
  six_month_sbp ~
    program_num * baseline_sbp_c10 +
    age_c10 + smoking_num + rural_num + health_literacy,
  data = causal_data
)

data_program <- transform(causal_data, program_num = 1)
data_usual <- transform(causal_data, program_num = 0)

predicted_y1 <- predict(outcome_model, newdata = data_program)
predicted_y0 <- predict(outcome_model, newdata = data_usual)

gcomp_means <- c(
  program = mean(predicted_y1),
  usual_care = mean(predicted_y0)
)
gcomp_ate <- unname(gcomp_means["program"] - gcomp_means["usual_care"])

gcomp_table <- data.frame(
  Quantity = c(
    "Everyone participates in the program",
    "Everyone receives usual care",
    "ATE: program minus usual care"
  ),
  Estimate = unname(c(gcomp_means, gcomp_ate)),
  Unit = "mmHg",
  row.names = NULL
)

knitr::kable(
  gcomp_table,
  col.names = c("Quantity", "Estimate", "Unit"),
  digits = 2,
  caption = "Standardized counterfactual population means and their ATE contrast"
)

Standardized counterfactual population means and their ATE contrast

Quantity	Estimate	Unit
Everyone participates in the program	133.27	mmHg
Everyone receives usual care	138.32	mmHg
ATE: program minus usual care	-5.05	mmHg

The standardized ATE is -5.05 mmHg, close to the finite-sample simulated truth of -5.06 mmHg. Each computational step is worth checking:

1. The model uses only the pretreatment adjustment set and prespecified effect modification.
2. data_program and data_usual retain the same participants and their values of $L$; only the treatment strategy changes.
3. Each participant’s outcome is predicted under both strategies.
4. Predictions are averaged over the target-population distribution and then contrasted.

7.2 Why the treatment regression coefficient may not equal the ATE

The model contains the program_num * baseline_sbp_c10 interaction. The program_num coefficient therefore represents the conditional contrast at a baseline SBP of 145 mmHg, whereas the ATE averages predicted contrasts across the target population’s distribution of baseline blood pressure.

Even without an interaction, conditional ratios from nonlinear models such as logistic or Cox regression generally do not equal marginal risks or population-average causal effects. A conditional OR, HR, or coefficient at a particular reference value should therefore not automatically be called an ATE.

conditional_program_coefficient <- coef(outcome_model)["program_num"]

coefficient_comparison <- data.frame(
  Quantity = c("Program main-effect coefficient", "Standardized ATE"),
  Estimate = unname(c(conditional_program_coefficient, gcomp_ate)),
  Interpretation = c(
    "Conditional effect at baseline SBP = 145 mmHg",
    "Marginal average effect in the current target population"
  )
)

knitr::kable(
  coefficient_comparison,
  digits = 2,
  caption = "Conditional regression coefficient and marginal standardized effect"
)

Conditional regression coefficient and marginal standardized effect

Quantity	Estimate	Interpretation
Program main-effect coefficient	-5.03	Conditional effect at baseline SBP = 145 mmHg
Standardized ATE	-5.05	Marginal average effect in the current target population

7.3 Using the nonparametric bootstrap to quantify sampling uncertainty

The standard error should reflect both outcome-model fitting and standardization. A nonparametric bootstrap provides an intuitive implementation by resampling participants, refitting the model, and repeating the standardization in every replicate.

estimate_gcomp <- function(data, index) {
  d <- data[index, , drop = FALSE]
  fit <- lm(
    six_month_sbp ~
      program_num * baseline_sbp_c10 +
      age_c10 + smoking_num + rural_num + health_literacy,
    data = d
  )
  d1 <- transform(d, program_num = 1)
  d0 <- transform(d, program_num = 0)
  individual_contrast <-
    predict(fit, newdata = d1) - predict(fit, newdata = d0)
  high_baseline <- d$baseline_sbp >= 150
  c(
    ATE = mean(individual_contrast),
    low_baseline = mean(individual_contrast[!high_baseline]),
    high_baseline = mean(individual_contrast[high_baseline])
  )
}

set.seed(20260812)
n_boot <- 250
gcomp_bootstrap <- replicate(
  n_boot,
  estimate_gcomp(causal_data, sample.int(nrow(causal_data), replace = TRUE))
)
gcomp_ci <- unname(
  quantile(gcomp_bootstrap["ATE", ], c(0.025, 0.975))
)

data.frame(
  Method = "Standardization",
  ATE = gcomp_ate,
  CI_lower = gcomp_ci[1],
  CI_upper = gcomp_ci[2],
  Bootstrap_replicates = n_boot,
  check.names = FALSE
) |>
  knitr::kable(
    col.names = c(
      "Method", "ATE", "Lower 95% CI", "Upper 95% CI",
      "Bootstrap replicates"
    ),
    digits = 2,
    caption = "Percentile bootstrap interval for the standardized ATE"
  )

Percentile bootstrap interval for the standardized ATE

Method	ATE	Lower 95% CI	Upper 95% CI	Bootstrap replicates
Standardization	-5.05	-5.63	-4.43	250

To keep the tutorial practical to render, this example uses only 250 bootstrap replicates. A formal analysis would generally use at least 1,000 replicates and assess Monte Carlo error. The interval adopts a superpopulation interpretation under independent and identically distributed sampling of participants. It quantifies statistical uncertainty from the current sampling and modeling process, but does not include identification uncertainty arising from unmeasured confounding, an incorrect DAG, treatment misclassification, selection bias, or an ambiguously defined intervention.

7.4 What diagnostics does the outcome model need?

Standardization depends on correctly estimating the conditional mean over the relevant data range. Check whether:

• continuous covariates have reasonable functional forms or require splines or other nonlinear representations;
• scientifically important treatment-by-covariate interactions have been omitted;
• predictions require substantial extrapolation beyond observed treatment-covariate combinations;
• residuals, unusual observations, clustering, or measurement error affect the model;
• the outcome type and link function align with the target effect scale; and
• complex models use cross-validation or sample splitting to limit overfitting.

Good outcome-model fit does not establish exchangeability. The model can address only confounders that were measured and represented appropriately.

8 Propensity Scores and Inverse Probability Weighting

8.1 Propensity scores describe the treatment assignment mechanism

The propensity score is defined as:

$$e(L)=P(A=1\mid L).$$

It is the probability of treatment given pretreatment covariates, not an outcome risk or an individual treatment effect. Its principal uses are to construct a pseudo-population, matched sample, or strata with comparable covariate distributions.

propensity_model <- glm(
  program_num ~
    age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
    health_literacy,
  family = binomial(),
  data = causal_data
)

propensity_hat <- clamp_probability(
  predict(propensity_model, type = "response"),
  epsilon = 1e-6
)

propensity_summary <- aggregate(
  propensity_hat,
  by = list(Strategy = causal_data$program),
  FUN = function(x) c(
    min = min(x), q25 = quantile(x, 0.25),
    median = median(x), q75 = quantile(x, 0.75), max = max(x)
  )
)

propensity_quantiles <- propensity_summary[[2]]
propensity_summary <- data.frame(
  Strategy = propensity_summary$Strategy,
  propensity_quantiles,
  row.names = NULL,
  check.names = FALSE
)

knitr::kable(
  propensity_summary,
  digits = 2,
  caption = "Distribution of estimated propensity scores in the two groups"
)

Distribution of estimated propensity scores in the two groups

Strategy	min	q25.25%	median	q75.75%	max
Usual care	0.05	0.24	0.34	0.46	0.88
Program	0.09	0.35	0.47	0.60	0.92

The goal of a propensity-score model is not to maximize classification accuracy or the AUC. A model that almost perfectly separates the treatment groups may instead signal poor overlap. Variable selection should follow the causal structure; after fitting, the focus should be on overlap and weighted balance.

8.2 Examine overlap before estimating the effect

hist(
  propensity_hat[causal_data$program_num == 0],
  breaks = seq(0, 1, by = 0.04), probability = TRUE,
  col = grDevices::adjustcolor(palette_ci["orange"], alpha.f = 0.50),
  border = "white", xlim = c(0, 1),
  xlab = "Estimated propensity score", ylab = "Density",
  main = "Propensity score overlap"
)
hist(
  propensity_hat[causal_data$program_num == 1],
  breaks = seq(0, 1, by = 0.04), probability = TRUE,
  col = grDevices::adjustcolor(palette_ci["teal"], alpha.f = 0.50),
  border = "white", add = TRUE
)
legend(
  "topright",
  legend = c("Usual care", "Program"),
  fill = grDevices::adjustcolor(
    c(palette_ci["orange"], palette_ci["teal"]), 0.50
  ),
  border = NA, bty = "n"
)

Distribution of estimated propensity scores by observed management strategy.

Within the region of common support, both groups contain participants with similar values of $L$. If a region contains only one group, estimating the ATE requires model extrapolation. Appropriate responses include redefining the target population, restricting the analysis to common support, changing the estimand, or acknowledging that the data cannot answer the original question. Any restriction must be accompanied by a description of the new target population.

8.3 Constructing stabilized inverse probability weights

The stabilized inverse probability of treatment weight for the ATE is:

$$SW_i=A_i\frac{P(A=1)}{e(L_i)}+ (1-A_i)\frac{P(A=0)}{1-e(L_i)}.$$

The denominator reweights treatment within covariate patterns, while the numerator stabilizes the scale of the weights. In the ideal weighted pseudo-population, treatment is approximately independent of the measured variables in $L$.

treatment_prevalence <- mean(causal_data$program_num)

stabilized_weight <- ifelse(
  causal_data$program_num == 1,
  treatment_prevalence / propensity_hat,
  (1 - treatment_prevalence) / (1 - propensity_hat)
)

effective_sample_size <- function(w) {
  sum(w)^2 / sum(w^2)
}

weight_summary <- data.frame(
  Metric = c(
    "Minimum", "1st percentile", "Median", "99th percentile", "Maximum",
    "Mean weight", "Overall effective sample size",
    "Program-group effective sample size",
    "Usual-care-group effective sample size"
  ),
  Value = unname(c(
    min(stabilized_weight),
    quantile(stabilized_weight, 0.01),
    median(stabilized_weight),
    quantile(stabilized_weight, 0.99),
    max(stabilized_weight),
    mean(stabilized_weight),
    effective_sample_size(stabilized_weight),
    effective_sample_size(stabilized_weight[causal_data$program_num == 1]),
    effective_sample_size(stabilized_weight[causal_data$program_num == 0])
  ))
)

knitr::kable(
  weight_summary,
  digits = 2,
  caption = "Diagnostics for stabilized inverse probability weights"
)

Diagnostics for stabilized inverse probability weights

Metric	Value
Minimum	0.44
1st percentile	0.50
Median	0.89
99th percentile	2.78
Maximum	4.84
Mean weight	1.00
Overall effective sample size	1848.07
Program-group effective sample size	696.86
Usual-care-group effective sample size	1157.13

The mean of stabilized weights should generally be close to 1. The effective sample size (ESS) is:

$$ESS=\frac{(\sum_i w_i)^2}{\sum_i w_i^2}.$$

An ESS substantially below the original sample size indicates that a small number of highly weighted observations dominate the estimate. ESS is a diagnostic summary, not the actual amount of independent information, and it does not replace a variance method appropriate for weighted estimation.

8.4 Reassess balance after weighting

weighted_smd <- vapply(
  balance_variables,
  standardized_difference,
  numeric(1),
  z = causal_data$program_num,
  w = stabilized_weight
)

balance_table <- data.frame(
  Covariate = names(unadjusted_smd),
  Unadjusted_SMD = unname(unadjusted_smd),
  Weighted_SMD = unname(weighted_smd),
  check.names = FALSE
)

knitr::kable(
  balance_table,
  col.names = c("Covariate", "Unadjusted SMD", "Weighted SMD"),
  digits = 2,
  caption = "Pretreatment covariate balance before and after inverse probability weighting"
)

Pretreatment covariate balance before and after inverse probability weighting

Covariate	Unadjusted SMD	Weighted SMD
Age	0.46	-0.02
Baseline SBP	0.64	-0.02
Current smoking	0.32	0.00
Rural residence	0.00	-0.01
Health literacy	0.08	0.00

old_margin <- par("mar")
par(mar = c(5.1, 8.2, 4.1, 2.1))
plot(
  abs(unadjusted_smd), seq_along(unadjusted_smd),
  pch = 16, col = palette_ci["orange"],
  xlim = c(0, max(0.35, abs(unadjusted_smd), abs(weighted_smd))),
  ylim = c(0.5, length(unadjusted_smd) + 0.5),
  yaxt = "n", xlab = "Absolute SMD", ylab = "",
  main = "Covariate balance"
)
balance_plot_labels <- c(
  "Age", "Baseline SBP", "Current smoking", "Rural residence",
  "Health literacy"
)
axis(2, at = seq_along(unadjusted_smd), labels = balance_plot_labels, las = 1)
points(abs(weighted_smd), seq_along(weighted_smd),
       pch = 17, col = palette_ci["teal"])
abline(v = 0.10, lty = 2, col = palette_ci["gray"])
legend(
  "topright", legend = c("Unadjusted", "Weighted", "0.10 reference"),
  pch = c(16, 17, NA), lty = c(NA, NA, 2),
  col = c(palette_ci["orange"], palette_ci["teal"], palette_ci["gray"]),
  bty = "n"
)

Absolute standardized mean differences before and after inverse probability weighting.

par(mar = old_margin)

The balance plot evaluates only variables that were measured and displayed. Excellent measured balance cannot rule out unmeasured confounding or repair an incorrect time zero, a problematic selection mechanism, or measurement error.

Diagnose propensity-score models by balance If important covariates remain imbalanced after weighting, revisit their coding, nonlinear terms, interactions, overlap, and the treatment assignment mechanism rather than first choosing the model that produces a more favorable outcome effect.

8.5 Estimating the ATE with Hájek-type weighted means

Normalizing the weighted means within treatment groups prevents finite-sample departures of the weight totals from their target sizes from carrying directly into the estimated means:

ipw_y1 <- weighted_mean(
  causal_data$six_month_sbp[causal_data$program_num == 1],
  stabilized_weight[causal_data$program_num == 1]
)
ipw_y0 <- weighted_mean(
  causal_data$six_month_sbp[causal_data$program_num == 0],
  stabilized_weight[causal_data$program_num == 0]
)
ipw_ate <- ipw_y1 - ipw_y0

ipw_table <- data.frame(
  Quantity = c("Program-strategy weighted mean", "Usual-care-strategy weighted mean", "IPW ATE"),
  Estimate = c(ipw_y1, ipw_y0, ipw_ate),
  Unit = "mmHg"
)

knitr::kable(
  ipw_table,
  digits = 2,
  caption = "Marginal mean outcomes estimated by inverse probability weighting"
)

Marginal mean outcomes estimated by inverse probability weighting

Quantity	Estimate	Unit
Program-strategy weighted mean	133.06	mmHg
Usual-care-strategy weighted mean	138.33	mmHg
IPW ATE	-5.27	mmHg

estimate_ipw <- function(data, index) {
  d <- data[index, , drop = FALSE]
  ps_fit <- glm(
    program_num ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    family = binomial(),
    data = d
  )
  ps <- clamp_probability(predict(ps_fit, type = "response"), 1e-6)
  prevalence <- mean(d$program_num)
  w <- ifelse(
    d$program_num == 1,
    prevalence / ps,
    (1 - prevalence) / (1 - ps)
  )
  weighted_mean(d$six_month_sbp[d$program_num == 1], w[d$program_num == 1]) -
    weighted_mean(d$six_month_sbp[d$program_num == 0], w[d$program_num == 0])
}

set.seed(20260813)
ipw_bootstrap <- replicate(
  n_boot,
  estimate_ipw(causal_data, sample.int(nrow(causal_data), replace = TRUE))
)
ipw_ci <- unname(quantile(ipw_bootstrap, c(0.025, 0.975)))

data.frame(
  Method = "Stabilized IPW",
  ATE = ipw_ate,
  CI_lower = ipw_ci[1],
  CI_upper = ipw_ci[2],
  Bootstrap_replicates = n_boot,
  check.names = FALSE
) |>
  knitr::kable(
    col.names = c(
      "Method", "ATE", "Lower 95% CI", "Upper 95% CI",
      "Bootstrap replicates"
    ),
    digits = 2,
    caption = "Percentile bootstrap interval for the IPW ATE"
  )

Percentile bootstrap interval for the IPW ATE

Method	ATE	Lower 95% CI	Upper 95% CI	Bootstrap replicates
Stabilized IPW	-5.27	-6.01	-4.54	250

Each bootstrap replicate re-estimates both the propensity scores and the weights. Treating weights from a single fitted model as fixed and then using the default standard error from an ordinary weighted lm() generally fails to capture the uncertainty introduced by estimating those weights.

9 Propensity Score Matching

9.1 Matching changes the comparison and may also change the target population

Matching finds controls whose pretreatment characteristics are similar to those of treated participants. Common 1:1 nearest-neighbor matching more naturally estimates the ATT among treated participants who were successfully matched, rather than the ATE in the original population. Once unmatched participants are excluded, the new analytic population must be described.

The following base R example demonstrates greedy matching without replacement on the logit of the propensity score. Applied analyses should use validated software so that processing order, ties, replacement, calipers, weights, and variance estimation are handled explicitly.

Continue with the complete PSM workflow This section connects causal identification to the idea of matching. For a complete MatchIt case study covering sample flow, overlap, SMDs, Love plots, matching weights, ATT inference, and design sensitivity, continue to Propensity Score Matching in Depth.

greedy_match <- function(z, logit_ps, caliper, seed = 1) {
  treated <- which(z == 1)
  available_controls <- which(z == 0)
  set.seed(seed)
  treated <- sample(treated)

  matched_treated <- integer(0)
  matched_control <- integer(0)

  for (treated_id in treated) {
    if (!length(available_controls)) break
    distance <- abs(logit_ps[available_controls] - logit_ps[treated_id])
    nearest_position <- which.min(distance)
    if (distance[nearest_position] <= caliper) {
      matched_treated <- c(matched_treated, treated_id)
      matched_control <- c(
        matched_control, available_controls[nearest_position]
      )
      available_controls <- available_controls[-nearest_position]
    }
  }

  data.frame(treated = matched_treated, control = matched_control)
}

logit_propensity <- qlogis(propensity_hat)
match_caliper <- 0.20 * sd(logit_propensity)
matched_pairs <- greedy_match(
  causal_data$program_num,
  logit_propensity,
  caliper = match_caliper,
  seed = 20260814
)

pair_difference <-
  causal_data$six_month_sbp[matched_pairs$treated] -
  causal_data$six_month_sbp[matched_pairs$control]
matched_att_ci <- mean_ci(pair_difference)
true_matched_att <- mean(
  causal_truth$individual_effect[matched_pairs$treated]
)

matched_result <- data.frame(
  Metric = c(
    "Number in the original program group",
    "Number of successfully matched pairs",
    "Number of unmatched program participants",
    "Matched ATT (demonstration)",
    "Lower bound of the illustrative 95% CI",
    "Upper bound of the illustrative 95% CI",
    "Finite-sample simulated ATT truth among successfully matched program participants",
    "Finite-sample simulated ATT truth among all program participants"
  ),
  Value = c(
    as.character(sum(causal_data$program_num == 1)),
    as.character(nrow(matched_pairs)),
    as.character(sum(causal_data$program_num == 1) - nrow(matched_pairs)),
    formatC(
      unname(c(
        matched_att_ci["estimate"], matched_att_ci["lower"],
        matched_att_ci["upper"], true_matched_att, true_att
      )),
      digits = 2, format = "f"
    )
  ),
  Unit = c("participants", "pairs", "participants", rep("mmHg", 5))
)

knitr::kable(
  matched_result,
  align = c("l", "r", "l"),
  caption = "Results of 1:1 nearest-neighbor propensity score matching"
)

Results of 1:1 nearest-neighbor propensity score matching

Metric	Value	Unit
Number in the original program group	888	participants
Number of successfully matched pairs	743	pairs
Number of unmatched program participants	145	participants
Matched ATT (demonstration)	-5.08	mmHg
Lower bound of the illustrative 95% CI	-5.92	mmHg
Upper bound of the illustrative 95% CI	-4.24	mmHg
Finite-sample simulated ATT truth among successfully matched program participants	-5.11	mmHg
Finite-sample simulated ATT truth among all program participants	-5.16	mmHg

matched_indices <- c(matched_pairs$treated, matched_pairs$control)
matched_z <- c(
  rep(1, nrow(matched_pairs)),
  rep(0, nrow(matched_pairs))
)

matched_smd <- vapply(
  balance_variables,
  function(x) standardized_difference(x[matched_indices], matched_z),
  numeric(1)
)

data.frame(
  Covariate = names(matched_smd),
  `SMD before matching` = unname(unadjusted_smd),
  `SMD after matching` = unname(matched_smd),
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "Balance of pretreatment covariates before and after matching"
  )

Balance of pretreatment covariates before and after matching

Covariate	SMD before matching	SMD after matching
Age	0.46	0.03
Baseline SBP	0.64	0.00
Current smoking	0.32	-0.03
Rural residence	0.00	0.03
Health literacy	0.08	0.02

The postmatching paired mean difference applies only to participants who were successfully matched. The paired t interval in the table treats the matched pairs as fixed and does not fully account for propensity score estimation or the discontinuity of the matching algorithm. It is therefore a teaching demonstration, not a general variance-estimation template for matched analyses. Different matching algorithms can produce different samples; investigators should prespecify the distance measure, caliper, use of replacement, matching ratio, processing order, number excluded, balance assessment, and a variance estimator appropriate for the design.

10 Doubly Robust Estimation: AIPW

10.1 Using both treatment and outcome models

Augmented inverse probability weighting (AIPW) combines standardization with IPW. Define $m_a(L)=E(Y\mid A=a,L)$. The ATE estimator is:

$$\widehat\psi_{AIPW}=\frac{1}{n}\sum_{i=1}^n\left[ \hat m_1(L_i)-\hat m_0(L_i) +\frac{A_i\{Y_i-\hat m_1(L_i)\}}{\hat e(L_i)} -\frac{(1-A_i)\{Y_i-\hat m_0(L_i)\}}{1-\hat e(L_i)} \right].$$

The first component is the contrast in outcome-model predictions; the final two terms use weighted residuals to correct prediction errors. Under standard regularity conditions, the estimator is consistent if either the propensity score model or the outcome model is correctly specified, which gives the method its “doubly robust” designation.

10.2 Using cross-fitting to reduce overfitting bias

Cross-fitting trains nuisance models in one subset of the data and generates propensity score and outcome predictions in a fold that was not used for fitting. This sample separation is especially important with complex machine-learning models; the example uses simple regressions to demonstrate the computational structure.

set.seed(20260815)
k_folds <- 5
fold_id <- integer(nrow(causal_data))
for (treatment_level in 0:1) {
  level_index <- which(causal_data$program_num == treatment_level)
  fold_id[level_index] <- sample(
    rep(seq_len(k_folds), length.out = length(level_index))
  )
}

crossfit_ps <- crossfit_m1 <- crossfit_m0 <- rep(NA_real_, nrow(causal_data))

for (fold in seq_len(k_folds)) {
  test_index <- which(fold_id == fold)
  train_data <- causal_data[fold_id != fold, , drop = FALSE]
  test_data <- causal_data[test_index, , drop = FALSE]

  ps_fit <- glm(
    program_num ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    family = binomial(),
    data = train_data
  )

  outcome_fit_1 <- lm(
    six_month_sbp ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    data = train_data[train_data$program_num == 1, , drop = FALSE]
  )
  outcome_fit_0 <- lm(
    six_month_sbp ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    data = train_data[train_data$program_num == 0, , drop = FALSE]
  )

  crossfit_ps[test_index] <- clamp_probability(
    predict(ps_fit, newdata = test_data, type = "response"),
    1e-6
  )
  crossfit_m1[test_index] <- predict(outcome_fit_1, newdata = test_data)
  crossfit_m0[test_index] <- predict(outcome_fit_0, newdata = test_data)
}

stopifnot(
  !anyNA(crossfit_ps),
  !anyNA(crossfit_m1),
  !anyNA(crossfit_m0)
)

aipw_score <-
  crossfit_m1 - crossfit_m0 +
  causal_data$program_num *
    (causal_data$six_month_sbp - crossfit_m1) / crossfit_ps -
  (1 - causal_data$program_num) *
    (causal_data$six_month_sbp - crossfit_m0) / (1 - crossfit_ps)

aipw_ate <- mean(aipw_score)
aipw_se <- sd(aipw_score) / sqrt(length(aipw_score))
aipw_ci <- aipw_ate + qnorm(c(0.025, 0.975)) * aipw_se

data.frame(
  Method = "Treatment-stratified 5-fold cross-fitted AIPW",
  ATE = aipw_ate,
  SE = aipw_se,
  `Lower 95% CI` = aipw_ci[1],
  `Upper 95% CI` = aipw_ci[2],
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "AIPW estimate of the population average treatment effect"
  )

AIPW estimate of the population average treatment effect

Method	ATE	SE	Lower 95% CI	Upper 95% CI
Treatment-stratified 5-fold cross-fitted AIPW	-5.04	0.35	-5.74	-4.35

This interval is a normal-approximation interval based on the empirical influence function for the cross-fitted AIPW estimator, under independent and identically distributed sampling of participants and standard regularity conditions. Clustering, repeated measurements, or complex sampling requires a corresponding variance method.

Doubly robust does not mean “double insurance” AIPW cannot repair unmeasured confounding, incorrect temporal ordering, a failure of consistency, a serious positivity violation, or simultaneous misspecification of both models. Extreme propensity scores can also amplify the residual terms, and cross-fitting cannot create treatment contrasts that do not exist in the data.

11 Bringing the Main ATE Estimates Together

11.1 Crude and adjusted estimates

crude_fit <- lm(six_month_sbp ~ program_num, data = causal_data)
crude_ci <- confint(crude_fit)["program_num", ]

ate_comparison <- data.frame(
  Method = c(
    "Crude mean difference", "Standardization", "Stabilized IPW",
    "Cross-fitted AIPW", "Finite-sample simulated ATE benchmark"
  ),
  Estimate = c(crude_difference, gcomp_ate, ipw_ate, aipw_ate, true_ate),
  `Lower 95% CI` = c(
    crude_ci[1], gcomp_ci[1], ipw_ci[1], aipw_ci[1], NA
  ),
  `Upper 95% CI` = c(
    crude_ci[2], gcomp_ci[2], ipw_ci[2], aipw_ci[2], NA
  ),
  Estimand = c(
    "Crude observed-group contrast (not a causal estimand)", rep("ATE", 4)
  ),
  check.names = FALSE
)

ate_comparison_display <- ate_comparison
for (column in c("Estimate", "Lower 95% CI", "Upper 95% CI")) {
  values <- ate_comparison_display[[column]]
  ate_comparison_display[[column]] <- ifelse(
    is.na(values), "—", formatC(values, digits = 2, format = "f")
  )
}

knitr::kable(
  ate_comparison_display,
  caption = "Crude observed association and three confounding-adjusted ATE estimates"
)

Crude observed association and three confounding-adjusted ATE estimates

Method	Estimate	Lower 95% CI	Upper 95% CI	Estimand
Crude mean difference	-0.65	-1.55	0.25	Crude observed-group contrast (not a causal estimand)
Standardization	-5.05	-5.63	-4.43	ATE
Stabilized IPW	-5.27	-6.01	-4.54	ATE
Cross-fitted AIPW	-5.04	-5.74	-4.35	ATE
Finite-sample simulated ATE benchmark	-5.06	—	—	ATE

plot_rows <- 4:1
comparison_plot_labels <- c(
  "Crude difference", "Standardization", "Stabilized IPW", "Cross-fitted AIPW"
)
old_margin <- par("mar")
par(mar = c(5.1, 8.2, 4.1, 2.1))
plot(
  ate_comparison$Estimate[1:4], plot_rows,
  pch = 16,
  col = c(palette_ci["orange"], rep(palette_ci["teal"], 3)),
  xlim = range(
    ate_comparison[["Lower 95% CI"]][1:4],
    ate_comparison[["Upper 95% CI"]][1:4],
    true_ate
  ),
  ylim = c(0.5, 4.5), yaxt = "n",
  xlab = "Program minus usual care: mean SBP difference (mmHg)", ylab = "",
  main = "Crude and adjusted estimates"
)
segments(
  ate_comparison[["Lower 95% CI"]][1:4], plot_rows,
  ate_comparison[["Upper 95% CI"]][1:4], plot_rows,
  col = c(palette_ci["orange"], rep(palette_ci["teal"], 3)),
  lwd = 2
)
axis(2, at = plot_rows, labels = comparison_plot_labels, las = 1)
abline(v = true_ate, lty = 2, lwd = 2, col = palette_ci["navy"])
abline(v = 0, lty = 3, col = palette_ci["gray"])
legend(
  "bottomright", legend = "Finite-sample ATE benchmark",
  lty = 2, lwd = 2, col = palette_ci["navy"], bty = "n"
)

Crude and adjusted estimates with the finite-sample simulated ATE benchmark.

par(mar = old_margin)

The three adjusted methods use different modeling strategies but rely on the same core identification assumptions. Similar results increase confidence in the implementations, but they do not prove that unmeasured confounding is absent. The crude estimate is biased toward zero because people with higher baseline blood pressure are more likely to participate in the program. The dashed benchmark is the finite-sample ATE in this realized simulated cohort, whereas the displayed intervals reflect the model-based, bootstrap, and influence-function sampling procedures described above; it is a teaching reference rather than the exact parameter those intervals are designed to cover.

The matching result is not included in this ATE plot because 1:1 matching targets successfully matched treated participants and is therefore closer to an ATT. Placing estimates for different target populations or effect scales in the same column would create a misleading appearance of disagreement among methods.

12 Effect Modification and Heterogeneity

12.1 “For whom does it work better?” is different from “Who is more likely to receive it?”

A confounder affects both treatment and outcome; an effect modifier makes the magnitude of the treatment effect differ across populations. A variable can play both roles or only one of them.

In the simulation mechanism, the program’s effect on blood pressure through engagement varies continuously with baseline SBP. Below, a prespecified 150 mmHg cutoff is used solely to demonstrate subgroup standardization. In an applied study, retaining the continuous information and presenting heterogeneity with curves and intervals would generally be preferable.

high_baseline <- causal_data$baseline_sbp >= 150

estimate_subgroup_effect <- function(index) {
  mean(predicted_y1[index] - predicted_y0[index])
}

subgroup_ci <- t(apply(
  gcomp_bootstrap[c("low_baseline", "high_baseline"), , drop = FALSE],
  1,
  quantile,
  probs = c(0.025, 0.975)
))

subgroup_effects <- data.frame(
  `Baseline group` = c("Baseline SBP < 150", "Baseline SBP ≥ 150"),
  N = c(sum(!high_baseline), sum(high_baseline)),
  `Standardized CATE` = c(
    estimate_subgroup_effect(!high_baseline),
    estimate_subgroup_effect(high_baseline)
  ),
  `Lower 95% CI` = subgroup_ci[, 1],
  `Upper 95% CI` = subgroup_ci[, 2],
  `Finite-sample simulated truth` = c(
    mean(causal_truth$individual_effect[!high_baseline]),
    mean(causal_truth$individual_effect[high_baseline])
  ),
  row.names = NULL,
  check.names = FALSE
)

knitr::kable(
  subgroup_effects,
  digits = 2,
  caption = "Conditional average treatment effects by baseline blood pressure group"
)

Conditional average treatment effects by baseline blood pressure group

Baseline group	N	Standardized CATE	Lower 95% CI	Upper 95% CI	Finite-sample simulated truth
Baseline SBP < 150	1375	-5.00	-5.73	-4.30	-4.89
Baseline SBP ≥ 150	825	-5.12	-5.88	-4.31	-5.34

Subgroup comparisons should report the number of participants, event counts or outcome distributions, effects, and uncertainty in each group. An interaction should not be claimed merely because one subgroup result is statistically significant and another is not. A formal test should directly evaluate the difference in effects, with the scale, functional form, and subgroups specified in advance.

12.2 Effect modification depends on the scale and target population

A treatment may produce a larger absolute benefit for high-risk people on the risk-difference scale while remaining approximately constant on the risk-ratio scale. The population ATE also changes with the proportions of subgroups in the target population. When transporting results, ask:

• Does the new population include covariate ranges that were not represented in the original study?
• Is the distribution of effect modifiers different?
• Are the treatment version, co-interventions, and outcome measurement the same?
• Are the mechanisms governing study participation or treatment access related to the target population?

A subgroup estimate is not automatically an individualized treatment rule. Individual treatment effects are generally unobservable, and using the same data to discover and report “best responders” is highly prone to overfitting. External validation and an evaluation of decision consequences are required.

13 Sensitivity Analyses for Positivity, Models, and Weights

13.1 Truncation must not be done quietly

Extreme weights increase variance and amplify measurement error in a small number of observations. Weight truncation can reduce variance but may reintroduce confounding bias; restricting the analysis to common support explicitly changes the target population. Both choices should be prespecified or transparently reported as sensitivity analyses.

truncation_limits <- quantile(stabilized_weight, c(0.01, 0.99))
truncated_weight <- pmin(
  pmax(stabilized_weight, truncation_limits[1]),
  truncation_limits[2]
)

truncated_ipw <-
  weighted_mean(
    causal_data$six_month_sbp[causal_data$program_num == 1],
    truncated_weight[causal_data$program_num == 1]
  ) -
  weighted_mean(
    causal_data$six_month_sbp[causal_data$program_num == 0],
    truncated_weight[causal_data$program_num == 0]
  )

support_index <- propensity_hat >= 0.10 & propensity_hat <= 0.90
support_data <- causal_data[support_index, , drop = FALSE]
support_ps_model <- glm(
  program_num ~
    age_c10 + baseline_sbp_c10 + smoking_num + rural_num + health_literacy,
  family = binomial(), data = support_data
)
support_ps <- clamp_probability(
  predict(support_ps_model, type = "response"), 1e-6
)
support_prevalence <- mean(support_data$program_num)
support_weight <- ifelse(
  support_data$program_num == 1,
  support_prevalence / support_ps,
  (1 - support_prevalence) / (1 - support_ps)
)
support_ipw <-
  weighted_mean(
    support_data$six_month_sbp[support_data$program_num == 1],
    support_weight[support_data$program_num == 1]
  ) -
  weighted_mean(
    support_data$six_month_sbp[support_data$program_num == 0],
    support_weight[support_data$program_num == 0]
  )

outcome_model_no_interaction <- lm(
  six_month_sbp ~
    program_num + age_c10 + baseline_sbp_c10 + smoking_num +
    rural_num + health_literacy,
  data = causal_data
)
no_interaction_ate <- coef(outcome_model_no_interaction)["program_num"]

sensitivity_table <- data.frame(
  Analysis = c(
    "Primary analysis: stabilized IPW",
    "Weight truncation at the 1st and 99th percentiles",
    "Restricted to 0.10 ≤ estimated PS ≤ 0.90",
    "Standardization with an outcome model without interactions"
  ),
  `Effect estimate` = c(
    ipw_ate, truncated_ipw, support_ipw, no_interaction_ate
  ),
  `Analysis N` = c(
    nrow(causal_data), nrow(causal_data), nrow(support_data),
    nrow(causal_data)
  ),
  `Target description` = c(
    "ATE in the original target population",
    "Same nominal ATE; modified weighting estimator and bias-variance tradeoff",
    "ATE in a restricted empirical population defined by the current data",
    "Misspecified-model sensitivity estimate for the original ATE"
  ),
  check.names = FALSE
)

knitr::kable(
  sensitivity_table,
  digits = 2,
  caption = "Sensitivity analyses for weights, propensity-score restriction, and outcome-model specification"
)

Sensitivity analyses for weights, propensity-score restriction, and outcome-model specification

Analysis	Effect estimate	Analysis N	Target description
Primary analysis: stabilized IPW	-5.27	2200	ATE in the original target population
Weight truncation at the 1st and 99th percentiles	-4.99	2200	Same nominal ATE; modified weighting estimator and bias-variance tradeoff
Restricted to 0.10 ≤ estimated PS ≤ 0.90	-5.21	2176	ATE in a restricted empirical population defined by the current data
Standardization with an outcome model without interactions	-5.06	2200	Misspecified-model sensitivity estimate for the original ATE

Similar numerical results do not establish that every assumption is correct, but they do indicate whether an arbitrary analytic choice dominates the conclusion. Propensity score thresholds define a data-dependent, restricted empirical population, not a common-support population in a strict theoretical sense. If many people are excluded, describe them separately and do not extrapolate the result to them. When the true effect is heterogeneous, the no-interaction OLS coefficient is generally an overlap-related model projection rather than the ATE in the original population; the final row intentionally demonstrates sensitivity to model misspecification.

13.2 A simplified bias grid for unmeasured confounding

Suppose there is an unmeasured binary factor $U$. After adequate adjustment for the measured variables, let $\Delta_U$ be the conditional difference in the prevalence of $U$ between the treatment and control groups, and let $\gamma$ be the additive mean difference in the outcome associated with $U$. In a simplified linear setting without complex interactions, the bias from omitting $U$ is approximately $\gamma\Delta_U$, giving the following bias-corrected estimate:

$$\widehat{ATE}_{corrected}\approx \widehat{ATE}_{observed}-\gamma\Delta_U.$$

delta_values <- c(-0.30, -0.15, 0, 0.15, 0.30)
gamma_values <- c(5, 10, 15)

bias_grid <- expand.grid(
  prevalence_difference = delta_values,
  outcome_mean_difference_for_U = gamma_values,
  KEEP.OUT.ATTRS = FALSE
)
bias_grid$bias_corrected_ATE <-
  aipw_ate -
  bias_grid$prevalence_difference *
  bias_grid$outcome_mean_difference_for_U

bias_display <- reshape(
  bias_grid,
  idvar = "prevalence_difference",
  timevar = "outcome_mean_difference_for_U",
  direction = "wide"
)
names(bias_display) <- c(
  "Treatment-control difference in U prevalence",
  "γ = 5", "γ = 10", "γ = 15"
)

knitr::kable(
  bias_display,
  digits = 2,
  caption = "Bias-corrected ATE under simplified unmeasured-confounding parameters (mmHg)"
)

Bias-corrected ATE under simplified unmeasured-confounding parameters (mmHg)

Treatment-control difference in U prevalence	γ = 5	γ = 10	γ = 15
-0.30	-3.54	-2.04	-0.54
-0.15	-4.29	-3.54	-2.79
0.00	-5.04	-5.04	-5.04
0.15	-5.79	-6.54	-7.29
0.30	-6.54	-8.04	-9.54

This grid is neither a test for unmeasured confounding nor a general correction formula. It forces investigators to state how imbalanced an omitted factor would need to be, how strongly and in what direction it would need to relate to the outcome, and whether those values would materially change the conclusion. A more formal quantitative bias analysis should be tailored to the outcome type, effect scale, interactions, measurement error, and available external information.

13.3 Negative controls and multiple lines of evidence

A negative-control exposure should not, in theory, affect the target outcome; a negative-control outcome should not, in theory, be affected by the target treatment. An observed association may indicate shared confounding, selection, or measurement problems. Negative controls themselves, however, depend on an explicit assumption that the relevant effect should be absent and therefore cannot automatically quantify or remove bias.

A persuasive causal argument usually comes from multiple lines of evidence: different designs, biases expected to operate in different directions, prespecified sensitivity analyses, negative controls, natural experiments, mechanistic evidence, and external replication—not from the p value of a single model.

14 Other Important Causal Designs

14.1 Different designs solve different problems

When conditional exchangeability given measured covariates is not credible, researchers may exploit policies, thresholds, changes over time, or external encouragement to form quasi-experimental comparisons. Each design replaces the assumptions of conventional covariate adjustment with a different set of strong assumptions.

Method	Typical estimand	Core basis for identification	Key diagnostic or threat
Regression/standardization	ATE or CATE in the target population	No unmeasured confounding given $L$	Functional form, overlap, adjustment set
IPW	Marginal effect in the target population	Correct treatment model plus the same identification assumptions	Balance, extreme weights, model error
AIPW	Marginal effect in the target population	Either the treatment or outcome model is correct, plus the same identification assumptions	Both models, extreme scores, cross-fitting
Matching	Often the ATT among matchable treated participants	Exchangeability on measured variables	Match quality, discarded observations, variance
Difference-in-differences	ATT for the treated group	Parallel trends in the absence of treatment	Pretrends, concurrent shocks, anticipation
Regression discontinuity	Local effect near a cutoff	Continuity of potential outcomes at the cutoff	Score manipulation, bandwidth, functional form
Instrumental variables	Often a local average treatment effect among compliers	Relevance, independence, exclusion restriction, and related assumptions	Weak instrument, direct pathways, monotonicity
Interrupted time series	Level or trend change at a policy date	No concurrent shock that also affects the outcome	Seasonality, historical events, autocorrelation

14.2 Difference-in-differences

Difference-in-differences (DiD) compares the pre-to-post change in the treated group with the corresponding change in a comparison group:

$$DiD=(\bar Y_{treated,post}-\bar Y_{treated,pre})- (\bar Y_{control,post}-\bar Y_{control,pre}).$$

The key parallel-trends assumption states that, without the intervention, the two groups would have followed the same average outcome trend. A prepolicy trend plot can reveal some evidence against this assumption, but it cannot prove the unobserved postpolicy counterfactual trend. Changes in group composition, differential concurrent shocks, anticipation, spillovers, and changes in outcome coding also require attention.

With staggered adoption and treatment effects that vary across groups or time, a traditional two-way fixed-effects regression can combine inappropriate comparisons. Use modern group-time estimators designed for staggered treatment, and state whether not-yet-treated or never-treated groups form the comparison.

14.3 Regression discontinuity

When a rule assigns treatment according to whether a continuous score $R$ crosses a cutoff $c$, units immediately on either side of the cutoff may be compared:

$$\tau_{RD}=\lim_{r\downarrow c}E(Y\mid R=r)- \lim_{r\uparrow c}E(Y\mid R=r).$$

This effect usually applies only near the cutoff. Identification requires potential outcomes and important baseline characteristics to be continuous at the cutoff, apart from the jump in treatment probability, and requires that individuals cannot precisely manipulate their score. Show raw observations or binned means, the score density, and covariate continuity. Report local-linear estimates across defensible bandwidths and polynomial orders; avoid high-order global polynomials.

14.4 Instrumental variables

An instrument $Z$ affects treatment $A$ but influences the outcome only through treatment. With a binary instrument and binary treatment, the Wald ratio is:

$$\widehat{LATE}= \frac{E(Y\mid Z=1)-E(Y\mid Z=0)} {E(A\mid Z=1)-E(A\mid Z=0)}.$$

Under instrument relevance, independence, the exclusion restriction, monotonicity, and related conditions, this ratio identifies the local average treatment effect among compliers. It is a complier mean difference for a continuous outcome and a complier risk difference for a binary outcome; no separate “linear risk-difference model” is required. A weak first stage produces instability and severe finite-sample bias, while strong relevance does not establish the exclusion restriction.

Clinician preference, distance, and policy eligibility are sometimes proposed as instruments, but each may affect outcomes directly through care quality, transportation, regional resources, or other services. An instrument is credible because of institutional and mechanistic knowledge, not because a variable has been entered into a two-stage regression.

14.5 Time-varying treatment and confounding

During long-term treatment, prior treatment can affect subsequent health status, which then affects the next treatment decision. If a time-varying covariate is both a consequence of earlier treatment and a confounder of later treatment and outcome, ordinary regression adjustment can block part of the treatment pathway and introduce bias.

Marginal structural models with inverse-probability treatment and censoring weights, the longitudinal g-formula, g-estimation, and related g-methods address this structure. They require treatment strategies, covariate histories, time zero, follow-up, censoring, and positivity to be defined at each time point. Treating multiple longitudinal rows as independent cross-sectional observations does not solve the problem.

14.6 Mediation analysis poses a new intervention question

A total effect can be decomposed into pathways through and outside a mediator, but a treatment coefficient “adjusted for the mediator” is not automatically a causal direct effect. Mediation analysis must additionally address treatment-mediator interaction, unmeasured mediator-outcome confounding, mediator-outcome confounding affected by treatment, and whether an intervention on the mediator is well defined.

Rather than asking vaguely what percentage is mediated, ask a more explicit question: what would happen to the outcome if the mediator distribution were shifted by an implementable intervention while treatment strategy remained fixed? Natural, controlled, and interventional direct and indirect effects have different counterfactual definitions and assumptions.

15 Missingness, Selection, Measurement, and Generalizability

15.1 Complete-case analysis changes who is analyzed

If outcome missingness is affected by both treatment and prognosis, analyzing only participants with observed outcomes conditions on a potential collider. When pretreatment confounders are missing, a complete-case analysis can also change the target population and its covariate distribution.

Options to consider include:

• improve follow-up and data linkage so that missingness is reduced by design;
• report the amount, reasons, and between-group pattern of missingness at every stage;
• use multiple imputation under an explicit missing-data mechanism;
• use inverse-probability-of-censoring weights for loss to follow-up;
• assess sensitivity to MNAR assumptions with extreme scenarios or pattern-mixture models;
• compare the populations and estimands represented by complete-case, imputed, and weighted analyses.

Multiple imputation does not automatically remove unmeasured confounding, and it cannot reliably recreate a key construct that was never collected.

15.2 Measurement error can leave residual confounding

Collapsing smoking to current versus not current, using a single measurement to represent long-term blood pressure, or using billing codes as a proxy for disease severity can all leave adjustment incomplete. Error in a confounder, treatment misclassification, and outcome error have different consequences; differential error is particularly difficult to predict.

Report each variable’s source, timing, repeated measurements, validation evidence, and threshold rules. When external information on sensitivity, specificity, or reliability is available, use probabilistic bias analysis rather than mentioning measurement error only in the limitations.

15.3 Generalizing from the study sample to a target population

Internal validity asks whether the causal contrast within the study is credible; external validity asks whether it applies to another population. If study participation $S$ is jointly affected by effect modifiers and outcome risk, the study-sample ATE may not represent the target-population ATE.

Generalizability or transportability methods can restandardize or reweight to the target population’s covariate distribution, but they require measurement of the important common causes of selection and outcome as well as effect modifiers, with support in the target population. When target-population data are unavailable, state the scope of applicability clearly.

Average effectiveness does not imply equitable access An average causal effect may hide group differences in benefit, harm, access, and measurement quality. Rurality, income, ethnicity, and disability often represent structures and resources rather than fixed biological attributes. When reporting heterogeneity, state the mechanistic assumptions, data support, and risk of stigmatization, and include “who can access the intervention?” in the decision.

16 An Auditable Causal-Analysis Workflow

16.1 Step 1: Specify the target trial before looking at outcomes

Define the target population, eligibility criteria, strategies, assignment, time zero, follow-up, outcome, causal contrast, effect scale, and analysis principle. Record when each variable is measured so that future information cannot enter baseline definitions.

16.2 Step 2: Draw the graph and register the assumptions

Invite subject-matter experts, data stewards, and representatives of the study population to review the DAG. List the minimally sufficient adjustment set, post-treatment variables that should not be adjusted for, potentially unmeasured common causes, selection mechanisms, and measurement proxies. Retain analysis plans for alternative credible DAGs.

16.3 Step 3: Assess feasibility before analyzing the outcome

Inspect treatment-group sizes, covariate distributions, missingness, propensity-score overlap, extreme combinations, follow-up, and outcome counts. If positivity clearly fails, redefine the question rather than hiding extrapolation with a more complex algorithm.

16.4 Step 4: Estimate and diagnose

Choose a method aligned with the estimand. For standardization, diagnose the outcome model; for IPW, diagnose the treatment model, weights, and balance; for matching, inspect the matched sample; and for AIPW, inspect both nuisance models. Use variance methods appropriate for clustering, repeated measures, and estimated weights.

16.5 Step 5: Probe analysis choices and identification fragility

At minimum, consider defensible DAGs, nonlinearity and interactions, overlap restrictions, weight truncation, missing-data handling, treatment and outcome error, unmeasured confounding, negative controls, and alternative effect scales. Sensitivity parameters should have scientifically defensible ranges rather than being searched until a preferred answer appears.

16.6 Step 6: Report absolute results, boundaries, and decision relevance

Report the standardized outcome under each strategy, the effect contrast, uncertainty interval, target population, diagnostics, and limitations. Separate identification assumptions from statistical-model assumptions, and do not imply that a confidence interval contains every source of uncertainty.

16.6.1 Reporting template

We emulated a comparison of [strategy 1] versus [strategy 0], initiated at [a common time zero], among [the target population]. The outcome was [a precisely defined outcome over a stated period], and the primary estimand was [ATE/ATT/CATE on a stated scale]. Based on [the DAG, subject-matter knowledge, and measurement timing], we adjusted for [variables]; identification relied on [consistency, exchangeability, positivity, no interference, and relevant selection/missingness conditions]. We estimated the effect with [standardization/IPW/AIPW/design-based method]; diagnostics showed [balance, overlap, weights, model behavior, and sample support]. Standardized outcomes under strategies 1 and 0 were [values], for an effect of [estimate, 95% CI]. Results were [robust/sensitive] to [sensitivity analyses]; [specific limitations involving unmeasured confounding, measurement, or transportability] may still affect interpretation.

16.7 Four-sentence result for this tutorial

Among 2200 simulated participants, we compared the 6-month mean systolic blood pressure if everyone entered the intensive blood pressure management program with the mean if everyone received usual care. Based on the pretreatment causal structure, we adjusted for age, baseline blood pressure, smoking, rural residence, and health literacy; cross-fitted AIPW estimated an ATE of -5.04 mmHg (95% CI -5.74 to -4.35). The standardization estimate was -5.05 mmHg, stabilized IPW gave -5.27 mmHg, and the largest absolute SMD among measured covariates after weighting was 0.02. This interval does not include identification uncertainty from unmeasured confounding or related biases; all results come from a simulation with a known mechanism and are not evidence about a real program.

17 Common Errors at a Glance

Common statement or practice	Why it is a problem	Better practice
“It is significant after adjustment, so it is causal.”	Significance does not validate identification assumptions	Define the estimand, design, DAG, and identification conditions first
Select confounders by univariable p-values	A p-value does not determine a causal role	Use temporal order, mechanisms, and a DAG
Adjust for every baseline and post-treatment variable	Mediators and colliders can change the question or induce bias	State each variable’s role and adjust for an appropriate set
A high propensity-score AUC proves a good model	Strong separation may indicate poor overlap	Inspect common support, weights, and weighted balance
Skip covariate checks after matching	An algorithm does not guarantee balance in the realized sample	Report distributions and SMDs before and after matching
Drop unmatched participants but still call the result a population ATE	The analysis population has changed	Describe exclusions and name the ATT or overlap-population effect correctly
Truncate weights without reporting the threshold	Obscures how the estimator and its bias-variance tradeoff changed	Report thresholds, counts, balance, and sensitivity results
“Doubly robust” means either model may be arbitrary	Bias remains when both nuisance models are wrong	Diagnose both models and the shared identification assumptions
One subgroup is significant and another is not	This does not establish a difference between subgroup effects	Estimate the interaction or effect difference and its interval directly
Report only a relative effect	Omits baseline risk and absolute decision relevance	Also report strategy-specific outcomes and an absolute effect
A confidence interval covers every uncertainty	It usually reflects sampling under the specified procedure	Analyze identification and measurement uncertainty separately
Similar results across methods prove no confounding	Methods may share the same incorrect assumptions	Seek evidence with different bias structures and external replication

18 Exercises and Answers

18.1 Exercise 1: Write the question as an estimand

What is missing from “Does the health program work?” Rewrite it as an estimable causal question.

Show answer

At minimum, the target population, well-defined program version, comparison strategy, common time zero, follow-up duration, outcome, effect scale, and target-population contrast are missing. Example: “Among adults newly diagnosed with hypertension in 2026 who meet the clinical eligibility criteria, what is the population mean difference in 6-month systolic blood pressure if all are offered the specified 6-month intensive program at diagnosis rather than usual care?”

18.2 Exercise 2: Identify variable roles

Baseline disease severity affects treatment receipt and the outcome; post-treatment adherence is affected by treatment and affects the outcome; follow-up visits are affected by both treatment and symptom worsening. What role does each variable play, and how should each be handled when estimating the total effect?

Show answer

Baseline severity is a confounder and should be controlled by design or analysis. Post-treatment adherence is a mediator and is generally not adjusted for when estimating the total effect. Follow-up visits are a collider of treatment and symptom worsening; restricting or conditioning on visits may open a noncausal pathway.

18.3 Exercise 3: Positivity or sample size?

A clinical rule prohibits a medication during pregnancy, but the study seeks a population ATE that includes pregnant people. Can collecting more data under the same rule solve the problem?

Show answer

No. This is a structural positivity violation: the subgroup will never have observations under treatment. Redefine the target population or treatment-strategy question, use a different design, or acknowledge that these data cannot identify that population ATE. A larger sample or more complex extrapolation does not create the missing counterfactual support.

18.4 Exercise 4: What does standardization do?

Why does standardization create predictions under both $A=1$ and $A=0$ for every participant rather than setting all continuous covariates to their sample means?

Show answer

The ATE is an average counterfactual contrast over the target-population distribution. Preserving each participant’s covariates, changing only treatment, and then averaging correctly integrates nonlinearity and effect modification. A prediction for the “average covariate profile” may describe no actual person and generally is not the population average.

18.5 Exercise 5: Interpret weight diagnostics

An IPW analysis has a maximum weight of 85, an ESS that falls from 3,000 to 420, and an absolute weighted SMD of 0.24 for an important covariate. Is the effect ready to report?

Show answer

Not without further work. These findings suggest poor overlap, model misspecification, or data-quality problems, and balance has not been achieved. Check nonlinearities and interactions in the treatment model, coding, common support, extreme covariate combinations, and the target population. If needed, redefine the estimand and report transparent truncation or restriction sensitivity analyses.

18.6 Exercise 6: ATE versus ATT

One-to-one matching without replacement excludes 35% of treated participants and 70% of controls. Can the matched mean difference still be called the original population ATE?

Show answer

Usually not. It more closely targets the ATT among successfully matched treated participants, depending on the matching weights and algorithm. Describe the exclusions, matched population, balance, and estimand. If many treated participants are excluded, the result may not even represent the ATT among all treated participants.

18.7 Exercise 7: Limits of double robustness

Both the propensity-score and outcome models in an AIPW analysis fit well. May an unmeasured disease-severity common cause now be ignored?

Show answer

No. Double robustness concerns correct specification of at least one nuisance model. If both models omit an unmeasured confounder, the core exchangeability condition still fails. Better measurement, a different design, negative controls, or quantitative bias analysis is needed.

18.8 Exercise 8: Subgroup “significance”

The effect has p = 0.03 in an urban subgroup and p = 0.20 in a rural subgroup. Can we conclude that only urban participants benefit?

Show answer

No. Two separate significance decisions do not test the difference between effects. Estimate the urban-rural effect difference or treatment interaction and its interval on a prespecified scale, while examining sample size, overlap, measurement, and multiplicity in both groups.

18.9 Exercise 9: Parallel trends in difference-in-differences

The prepolicy trends appear parallel. Does this prove that the postpolicy counterfactual trends would also have been parallel?

Show answer

No. Prepolicy data can reveal evidence against parallel trends and inform plausibility, but the untreated postpolicy counterfactual remains unobserved. Also justify the absence of differential concurrent shocks, anticipation, changing group composition, and spillovers.

18.10 Exercise 10: Interpret a confidence interval

The AIPW 95% CI excludes zero. Does this mean unmeasured confounding, selection bias, and misclassification cannot reverse the conclusion?

Show answer

No. A conventional interval mainly quantifies sampling uncertainty given the data, models, and identification conditions. Systematic bias can shift the entire interval. Report statistical uncertainty separately from sensitivity to identification, measurement, selection, and transportability.

19 Quick Reference

19.1 Core formulas

Concept	Formula	Interpretation
Individual effect	$Y_i(1)-Y_i(0)$	Counterfactual difference between two strategies for one unit
ATE	$E\{Y(1)-Y(0)\}$	Average effect in the target population
ATT	$E\{Y(1)-Y(0)\mid A=1\}$	Average effect among treated participants
Consistency	$A=a\Rightarrow Y=Y(a)$	The observed strategy maps to a well-defined potential outcome
Exchangeability	$\{Y(1),Y(0)\}\perp A\mid L$	No uncontrolled common causes conditional on $L$
Positivity	$0<P(A=a\mid L)<1$	Every target covariate pattern supports both strategies
g-formula	$E[Y(a)]=E_L[E(Y\mid A=a,L)]$	Standardize conditional outcomes to the target population
Propensity score	$e(L)=P(A=1\mid L)$	Probability of treatment given pretreatment covariates
Stabilized weight	$P(A=A_i)/P(A=A_i\mid L_i)$	Create a pseudo-population balanced on measured variables
ESS	$(\sum w)^2/\sum w^2$	Summary of weight concentration

19.2 Common base R patterns

Goal	Code pattern
Crude mean difference	with(d, mean(y[a == 1]) - mean(y[a == 0]))
Outcome model	lm(y ~ a * x1 + x2, data = d)
Two counterfactual worlds	transform(d, a = 1); transform(d, a = 0)
Standardized ATE	mean(predict(fit, d1) - predict(fit, d0))
Propensity score	glm(a ~ x1 + x2, family = binomial(), data = d)
ATE weights	ifelse(a == 1, mean(a)/ps, (1-mean(a))/(1-ps))
Strategy-specific weighted mean	mu1 <- sum(w[a == 1] * y[a == 1]) / sum(w[a == 1]); calculate mu0 analogously for a == 0
IPW ATE	mu1 - mu0
Weight ESS	sum(w)^2 / sum(w^2)
Nonparametric bootstrap	Resample with sample.int(nrow(d), replace = TRUE) and refit
Linear-model coefficient interval	confint(fit)

19.3 Glossary

Term	Plain-language meaning
Counterfactual/potential outcome	The outcome that would occur for the same unit under another strategy
Estimand	A quantity jointly defined by the target population, strategies, outcome, time, and effect scale
Confounding	Treatment and outcome share causes, so the observed groups are not directly comparable
Backdoor path	A noncausal path between treatment and outcome that begins with an arrow into treatment
Mediator	A variable on a causal pathway from treatment to outcome
Collider	A variable receiving two arrows; conditioning on it can open a path
Standardization	Average conditional counterfactual predictions over the target covariate distribution
Propensity score	Probability of treatment given pretreatment covariates
Positivity	Every type of target individual has some chance of receiving each strategy
Overlap	Both groups contain participants with similar pretreatment characteristics
Double robustness	The estimator remains consistent when either the treatment or outcome nuisance model is correct, under the shared identification assumptions
Cross-fitting	Train nuisance models and generate predictions in separate data folds to reduce overfitting bias
Local effect	An effect applying only near a threshold or to a subgroup such as compliers

19.4 Pre-analysis checklist

• Are the target population, eligibility criteria, and treatment strategies implementable and explicit?
• Are eligibility assessment, treatment assignment, and follow-up initiation aligned at the same time zero?
• Are the outcome, follow-up window, competing events, and censoring prespecified?
• Is the main estimand an ATE, ATT, CATE, or local effect, and on what scale?
• Is the measurement time clear for every covariate, mediator, selection indicator, and outcome?
• Is the DAG supported by subject-matter knowledge, with plausible alternatives considered?
• Does the adjustment set block all known backdoor paths without including inappropriate post-treatment variables?
• Are key confounders measured with sufficient quality and completeness?
• Is there structural or practical nonpositivity in the sample?
• Were the analysis plan, subgroups, and sensitivity ranges set before inspecting results?

19.5 Post-analysis checklist

• Does every reported estimate correspond to the stated target population and estimand?
• Are the marginal outcomes under both strategies reported along with their contrast?
• Were functional form, interactions, residuals, and extrapolation assessed for standardization?
• Were full distributions and balance of important covariates checked after weighting or matching?
• Were propensity-score overlap, weight quantiles, the maximum weight, and ESS reported?
• How did exclusions, failed matches, or support restrictions change who was analyzed, and how did truncation change each participant’s contribution?
• Does the variance method account for estimated weights, matching, clustering, and repeated measures?
• Were missingness, loss to follow-up, misclassification, and selection mechanisms addressed directly?
• Do sensitivity parameters for unmeasured confounding have a scientific basis and explicit direction?
• Are sampling uncertainty and identification uncertainty distinguished?
• Is the average result kept within populations supported by the data rather than overextended to individual decisions?
• Are code, seeds, data processing, and the software environment sufficient for reproduction?

Final knowledge check

Before publishing a causal conclusion, answer each question in one sentence:

1. What exactly would it mean for every target individual to receive strategy 1?
2. What is the comparison strategy, and where is time zero?
3. Does the effect target an ATE, ATT, CATE, or another population?
4. Which pathways create confounding, and why does the adjustment set block them?
5. Which key variables should not be adjusted for, and why?
6. Which identification assumptions cannot be tested directly with these data?
7. In which parts of the population do the data lack treatment overlap?
8. Does the main estimate depend on a few extreme observations or model extrapolation?
9. How strong and in what direction would an unmeasured bias need to be to alter the decision?
10. To whom can the result be generalized, and to whom can it not?

If any answer is only “the software handled it” or “because p < 0.05,” the analysis is not finished.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60
##  [5] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31
##  [9] lifecycle_1.0.5 cli_3.6.6       sass_0.4.10     jquerylib_0.1.4
## [13] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [17] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

Teaching use: This tutorial demonstrates the design, identification, estimation, diagnostics, and reporting of causal questions. All data and ground-truth effects are simulated and are not evidence about an intensive blood pressure management program in any real population. A real study should involve subject-matter experts, data auditing, a prespecified analysis plan, validated software, and ethical review.