1 Today’s Question

Lecture 2 used inverse probability weighting:

Model the treatment: who is likely to be a smoker, based on baseline covariates?

Lecture 3 used outcome regression:

Model the outcome: who is likely to have CVD, based on smoking and baseline covariates?

Lecture 4 combines both ideas:

Can we use both the treatment model and the outcome model together?

This is called doubly robust estimation.

Our treatment, outcome, and covariates are the same as before:

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

As before, this teaching analysis does not use survey weights.

2 The Big Picture

Doubly robust estimation uses two helpers.

The name “doubly robust” comes from a powerful idea:

If either the treatment model is good enough or the outcome model is good enough, the estimator can still work well.

But it is not magic. If both models are badly wrong, the estimate can still be wrong.

3 Why Add a Correction?

Outcome regression predicts risk under T = 0 and under T = 1. Doubly robust estimation starts there, then adds a correction using the observed outcomes.

If the outcome model predicts very well, the correction is small. If the outcome model misses something but the treatment model is good, the correction can help repair the estimate.

4 Load and Prepare the Data

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

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

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

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

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

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

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

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

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

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

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

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

nrow(dr_data)
## [1] 6299

We use the same complete-case sample as Lectures 2 and 3.

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

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

5 Step 1: Fit the Treatment Model

The treatment model estimates the propensity score:

Propensity score = P(smoker | age, gender, race, education, income, BMI)

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

dr_data$propensity_score <- predict(treatment_model, type = "response")

summary(dr_data$propensity_score)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05416 0.26826 0.39915 0.40625 0.53934 0.85971

6 Step 2: Fit the Outcome Model

The outcome model estimates CVD risk:

Outcome model = E(CVD | smoking, age, gender, race, education, income, BMI)

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

dr_data$predicted_risk_observed <- predict(outcome_model, type = "response")

data_if_non_smoker <- dr_data
data_if_smoker <- dr_data

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

dr_data$predicted_risk_if_non_smoker <- predict(
  outcome_model,
  newdata = data_if_non_smoker,
  type = "response"
)

dr_data$predicted_risk_if_smoker <- predict(
  outcome_model,
  newdata = data_if_smoker,
  type = "response"
)

outcome_model_summary <- data.frame(
  Quantity = c(
    "Average predicted risk if everyone non-smoker",
    "Average predicted risk if everyone smoker"
  ),
  Percent = round(
    100 * c(
      mean(dr_data$predicted_risk_if_non_smoker),
      mean(dr_data$predicted_risk_if_smoker)
    ),
    1
  )
)

knitr::kable(outcome_model_summary)
Quantity Percent
Average predicted risk if everyone non-smoker 8.5
Average predicted risk if everyone smoker 12.0

7 Step 3: Build the Doubly Robust Estimate

The doubly robust estimator starts with the outcome-regression prediction and adds a correction.

For the smoker world, T = 1:

DR risk under T = 1 = average(predicted risk if T = 1 + smoker correction)

For the non-smoker world, T = 0:

DR risk under T = 0 = average(predicted risk if T = 0 + non-smoker correction)

The correction uses the observed outcome minus the model prediction. That difference is called a residual.

T_obs <- dr_data$smoker_indicator
Y_obs <- dr_data$cvd_indicator
e_hat <- dr_data$propensity_score
m1_hat <- dr_data$predicted_risk_if_smoker
m0_hat <- dr_data$predicted_risk_if_non_smoker

dr_data$correction_if_smoker <- T_obs / e_hat * (Y_obs - m1_hat)
dr_data$correction_if_non_smoker <- (1 - T_obs) / (1 - e_hat) *
  (Y_obs - m0_hat)

dr_risk_smoker <- mean(m1_hat + dr_data$correction_if_smoker)
dr_risk_non_smoker <- mean(m0_hat + dr_data$correction_if_non_smoker)

dr_difference <- dr_risk_smoker - dr_risk_non_smoker
dr_ratio <- dr_risk_smoker / dr_risk_non_smoker

dr_components <- data.frame(
  Quantity = c(
    "Outcome-model risk if everyone non-smoker",
    "Average correction for non-smoker world",
    "Doubly robust risk if everyone non-smoker",
    "Outcome-model risk if everyone smoker",
    "Average correction for smoker world",
    "Doubly robust risk if everyone smoker"
  ),
  Percent = round(
    100 * c(
      mean(m0_hat),
      mean(dr_data$correction_if_non_smoker),
      dr_risk_non_smoker,
      mean(m1_hat),
      mean(dr_data$correction_if_smoker),
      dr_risk_smoker
    ),
    2
  )
)

knitr::kable(dr_components)
Quantity Percent
Outcome-model risk if everyone non-smoker 8.53
Average correction for non-smoker world 0.03
Doubly robust risk if everyone non-smoker 8.56
Outcome-model risk if everyone smoker 11.98
Average correction for smoker world -0.17
Doubly robust risk if everyone smoker 11.82

In this dataset, the correction is small. That means the IPW and outcome-regression pieces are giving fairly similar answers.

8 Picture of the Doubly Robust Calculation

9 Compare All Four Methods

Now we compare:

  1. Simple observed comparison
  2. IPW
  3. Outcome regression
  4. Doubly robust estimation
smoker_rows <- dr_data$smoking_group == "Smoker"
non_smoker_rows <- dr_data$smoking_group == "Non-smoker"

observed_risk_smoker <- mean(dr_data$cvd_indicator[smoker_rows])
observed_risk_non_smoker <- mean(dr_data$cvd_indicator[non_smoker_rows])
observed_difference <- observed_risk_smoker - observed_risk_non_smoker
observed_ratio <- observed_risk_smoker / observed_risk_non_smoker

ipw_weight <- ifelse(
  dr_data$smoker_indicator == 1,
  1 / dr_data$propensity_score,
  1 / (1 - dr_data$propensity_score)
)

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

ipw_risk_smoker <- weighted_risk(
  dr_data$cvd_indicator[smoker_rows],
  ipw_weight[smoker_rows]
)

ipw_risk_non_smoker <- weighted_risk(
  dr_data$cvd_indicator[non_smoker_rows],
  ipw_weight[non_smoker_rows]
)

ipw_difference <- ipw_risk_smoker - ipw_risk_non_smoker
ipw_ratio <- ipw_risk_smoker / ipw_risk_non_smoker

or_risk_smoker <- mean(dr_data$predicted_risk_if_smoker)
or_risk_non_smoker <- mean(dr_data$predicted_risk_if_non_smoker)
or_difference <- or_risk_smoker - or_risk_non_smoker
or_ratio <- or_risk_smoker / or_risk_non_smoker

comparison_table <- data.frame(
  Method = c(
    "Observed comparison",
    "IPW",
    "Outcome regression",
    "Doubly robust"
  ),
  Non_smoker_risk_percent = round(
    100 * c(
      observed_risk_non_smoker,
      ipw_risk_non_smoker,
      or_risk_non_smoker,
      dr_risk_non_smoker
    ),
    1
  ),
  Smoker_risk_percent = round(
    100 * c(
      observed_risk_smoker,
      ipw_risk_smoker,
      or_risk_smoker,
      dr_risk_smoker
    ),
    1
  ),
  Difference_percentage_points = round(
    100 * c(
      observed_difference,
      ipw_difference,
      or_difference,
      dr_difference
    ),
    1
  ),
  Risk_ratio = round(c(observed_ratio, ipw_ratio, or_ratio, dr_ratio), 2)
)

knitr::kable(
  comparison_table,
  col.names = c(
    "Method",
    "Non-smoker risk (%)",
    "Smoker risk (%)",
    "Difference (percentage points)",
    "Risk ratio"
  )
)
Method Non-smoker risk (%) Smoker risk (%) Difference (percentage points) Risk ratio
Observed comparison 6.9 15.1 8.2 2.18
IPW 8.7 11.7 3.0 1.34
Outcome regression 8.5 12.0 3.5 1.40
Doubly robust 8.6 11.8 3.3 1.38

The doubly robust estimate is close to the IPW and outcome-regression estimates in this example. That is reassuring because two different adjustment strategies tell a similar story.

10 Interpretation

The simple observed comparison estimates a CVD risk difference of 8.2 percentage points.

After adjustment:

Method Adjusted risk difference
IPW 3 percentage points
Outcome regression 3.5 percentage points
Doubly robust 3.3 percentage points

The doubly robust estimate says:

In this complete-case classroom sample, the adjusted predicted CVD risk is 8.6% if everyone were set to non-smoker and 11.8% if everyone were set to smoker.

The adjusted difference is 3.3 percentage points.

This is smaller than the raw observed difference. That suggests that some of the raw smoker versus non-smoker gap was related to baseline covariate differences between the groups.

11 What Can Go Wrong?

Doubly robust does not mean “always correct.”

We still need important causal assumptions:

  1. No important unmeasured confounders.
  2. Good overlap: people with similar X values could plausibly be smokers or non-smokers.
  3. At least one of the two models is reasonable.
  4. The data variables are measured well enough for the teaching question.

12 Key Takeaways

  1. IPW uses a treatment model.
  2. Outcome regression uses an outcome model.
  3. Doubly robust estimation uses both models.
  4. It starts with model-predicted risks and adds a correction from observed outcomes.
  5. It can work well if at least one of the two models is reasonable, but it is not magic.

13 Appendix: Math of the Doubly Robust Estimator

This appendix gives the mathematical notation for the doubly robust estimator used in this lecture.

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

\[ \begin{aligned} T_i &= \text{smoking indicator}, \\ Y_i &= \text{CVD indicator}, \\ X_i &= \text{baseline covariates}. \end{aligned} \]

Here \(T_i = 1\) means smoker, \(T_i = 0\) means non-smoker, and \(Y_i = 1\) means the person has CVD.

Our goal is to estimate two adjusted risks:

\[ \begin{aligned} \psi_1 &= \text{CVD risk if everyone smoked}, \\ \psi_0 &= \text{CVD risk if everyone did not smoke}. \end{aligned} \]

Then the adjusted risk difference and adjusted risk ratio are:

\[ \begin{aligned} RD &= \psi_1 - \psi_0, \\ RR &= \frac{\psi_1}{\psi_0}. \end{aligned} \]

13.1 Appendix 1: Treatment Model

The treatment model estimates the propensity score:

\[ e(X_i) = P(T_i = 1 \mid X_i). \]

In this lecture, we fit a logistic regression:

\[ \text{logit}\{e(X_i)\} = \alpha_0 + \alpha^\top X_i. \]

The fitted propensity score is:

\[ \hat e_i = \hat e(X_i). \]

13.2 Appendix 2: Outcome Model

The outcome model estimates CVD risk given smoking and covariates:

\[ m(t, X_i) = P(Y_i = 1 \mid T_i = t, X_i). \]

In this lecture, we fit another logistic regression:

\[ \text{logit}\{m(T_i, X_i)\} = \beta_0 + \beta_1 T_i + \beta^\top X_i. \]

Then we predict each person’s CVD risk under both treatment choices:

\[ \begin{aligned} \hat m_1(X_i) &= \hat m(1, X_i), \\ \hat m_0(X_i) &= \hat m(0, X_i). \end{aligned} \]

Here \(\hat m_1(X_i)\) means predicted CVD risk if person \(i\) were set to smoker, and \(\hat m_0(X_i)\) means predicted CVD risk if person \(i\) were set to non-smoker.

13.3 Appendix 3: Doubly Robust Risk Under Smoking

For the smoker world, \(T = 1\), the doubly robust estimator is:

\[ \hat\psi_1^{DR} = \frac{1}{n} \sum_{i=1}^{n} \left[ \hat m_1(X_i) + \frac{T_i}{\hat e_i} \{Y_i - \hat m_1(X_i)\} \right]. \]

This has two pieces:

\[ \text{prediction} = \hat m_1(X_i), \]

and:

\[ \text{correction} = \frac{T_i}{\hat e_i} \{Y_i - \hat m_1(X_i)\}. \]

The correction compares the observed outcome \(Y_i\) with the predicted outcome \(\hat m_1(X_i)\), but it only contributes for people who were actually smokers.

13.4 Appendix 4: Doubly Robust Risk Under Non-Smoking

For the non-smoker world, \(T = 0\), the doubly robust estimator is:

\[ \hat\psi_0^{DR} = \frac{1}{n} \sum_{i=1}^{n} \left[ \hat m_0(X_i) + \frac{1 - T_i}{1 - \hat e_i} \{Y_i - \hat m_0(X_i)\} \right]. \]

Again, this has two pieces:

\[ \text{prediction} = \hat m_0(X_i), \]

and:

\[ \text{correction} = \frac{1 - T_i}{1 - \hat e_i} \{Y_i - \hat m_0(X_i)\}. \]

The correction only contributes for people who were actually non-smokers.

13.5 Appendix 5: Risk Difference and Risk Ratio

Once we estimate the two adjusted risks, the doubly robust risk difference is:

\[ \widehat{RD}^{DR} = \hat\psi_1^{DR} - \hat\psi_0^{DR}. \]

The doubly robust risk ratio is:

\[ \widehat{RR}^{DR} = \frac{\hat\psi_1^{DR}} {\hat\psi_0^{DR}}. \]

In the main lecture, the risk difference is reported in percentage points:

\[ 100 \times \widehat{RD}^{DR}. \]

13.6 Appendix 6: Why the Estimator Is Called Doubly Robust

The doubly robust estimator combines two models:

\[ \hat e(X_i) \quad \text{and} \quad \hat m(t, X_i). \]

It is called doubly robust because it can still estimate the target risk correctly if either one of these models is correct.

For this statement, we still need the usual causal assumptions: consistency, no important unmeasured confounding after adjustment for \(X_i\), and positivity. Let \(Y_i^1\) be the potential CVD outcome if person \(i\) smoked, and let \(Y_i^0\) be the potential CVD outcome if person \(i\) did not smoke.

If the outcome model is correct, then:

\[ E\{Y_i - m(1, X_i) \mid T_i = 1, X_i\} = 0. \]

So the average correction term for the smoker world is centered around zero:

\[ E\left[ \frac{T_i}{e(X_i)} \{Y_i - m(1, X_i)\} \right] = 0. \]

In that case, the outcome-model prediction is already correct enough.

If the treatment model is correct, then the inverse-probability correction reweights the observed smokers so they represent the full covariate distribution:

\[ E\left[ \frac{T_i}{e(X_i)}Y_i \right] = E(Y_i^{1}). \]

Similarly:

\[ E\left[ \frac{1 - T_i}{1 - e(X_i)}Y_i \right] = E(Y_i^{0}). \]

So if the treatment model is correct, the correction can protect us even when the outcome model is imperfect.

The important classroom message is:

Doubly robust means either the treatment model or the outcome model can be correct enough. It does not mean both models can be wrong.

14 References

Bang, H., & Robins, J. M. (2005). Doubly robust estimation in missing data and causal inference models. Biometrics, 61(4), 962-973.

Hernan, M. A., & Robins, J. M. (2020). Causal Inference: What If. Chapman & Hall/CRC.

Robins, J. M., Rotnitzky, A., & Zhao, L. P. (1994). Estimation of regression coefficients when some regressors are not always observed. Journal of the American Statistical Association, 89(427), 846-866.