In Lecture 5, we used bootstrapping to ask:
How much random sampling uncertainty do we have?
Lecture 6 asks a different question:
What if there is an important confounder we did not measure?
This is called sensitivity analysis.
The goal is not to prove that our answer is perfect. The goal is to ask:
How sensitive is our conclusion to a possible hidden problem?
For this setting, we use the E-value.
Yes, E-value is possible here.
But there is one important detail:
E-values are naturally calculated for a risk ratio, not a risk difference.
In earlier lectures, we focused on the risk difference:
risk among smokers - risk among non-smokers
For the E-value, we also calculate the adjusted risk ratio:
risk among smokers / risk among non-smokers
Then the E-value asks:
How strong would an unmeasured confounder need to be to explain away the adjusted association?
This is exactly the kind of question we want for a sensitivity analysis. We are not pretending that all confounders were measured perfectly. Instead, we ask how powerful a missing confounder would have to be before our adjusted result could disappear.
Imagine our adjusted analysis gives this result:
Adjusted risk ratio = 1.4
E-value is about 2.1
The risk ratio of 1.4 means:
After adjustment, smoker CVD risk is about 1.4 times the non-smoker CVD risk.
The E-value of 2.1 means:
To make this result disappear, a missing confounder would need to have two strong links.
The missing confounder would need to be:
Both links matter. If the hidden factor is strongly related to smoking but barely related to CVD, it cannot explain much. If it strongly raises CVD risk but is equally common among smokers and non-smokers, it also cannot explain much.
So the easiest interpretation is:
Higher E-value = the hidden confounder would need to be stronger.
For our example, E-value = 2.1 is
moderate. The result is not extremely fragile, but it
is not invincible either.
A simple reading guide:
| E-value size | Easy interpretation |
|---|---|
| Close to 1 | The result is not very robust. A small hidden confounder could explain it away. |
| Around 2 | A moderately strong hidden confounder would be needed. |
| Much bigger than 2 | A very strong hidden confounder would be needed. The result is harder to explain away. |
For our smoking and CVD analysis, the E-values are around 2.0 to 2.2. That means:
The adjusted result is moderately robust, but not invincible.
In plain language:
A hidden factor would need to be pretty strongly related to both smoking and CVD to fully erase the adjusted association.
We already adjusted for measured baseline covariates:
But maybe there is another variable we did not measure. Call it
U.
If U affects both smoking and CVD, it could act like a
hidden confounder.
Examples might include long-term diet, family history, stress, physical activity, or other health behaviors not measured well in our teaching dataset.
The E-value is a number on the risk ratio scale. VanderWeele and Ding introduced it as a simple way to summarize sensitivity to possible unmeasured confounding.
For an observed risk ratio above 1:
E-value = RR + sqrt(RR * (RR - 1))
Here, RR means the adjusted risk ratio:
RR = CVD risk if everyone smoked / CVD risk if everyone did not smoke
The E-value tells us the minimum common strength of association an unmeasured confounder would need to have with both:
to explain away the adjusted association, after the measured covariates have already been adjusted for.
The phrase after the measured covariates have already been adjusted for is important. The hidden confounder would need to add information beyond age, gender, race, education, income-to-poverty ratio, and BMI.
The E-value is not a p-value. It is not a probability. It is not a new estimate of CVD risk. It is a robustness benchmark.
data_path <- file.path("..", "Data", "NHANES_data.csv")
if (!file.exists(data_path)) {
data_path <- file.path("Data", "NHANES_data.csv")
}
nhanes <- read.csv(data_path, stringsAsFactors = FALSE)
baseline_vars <- c(
"age_yr",
"gender",
"race",
"educ_lvl",
"inc_to_pov_ratio",
"bmi"
)
analysis_vars <- c("smoker_indicator", "cvd_indicator", baseline_vars)
analysis_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]
analysis_data$smoker_indicator <- as.integer(analysis_data$smoker_indicator)
analysis_data$cvd_indicator <- as.integer(analysis_data$cvd_indicator)
analysis_data$smoking_group <- ifelse(
analysis_data$smoker_indicator == 1,
"Smoker",
"Non-smoker"
)
analysis_data$smoking_group <- factor(
analysis_data$smoking_group,
levels = c("Non-smoker", "Smoker")
)
analysis_data$gender_label <- factor(
analysis_data$gender,
levels = c("F", "M"),
labels = c("Female", "Male")
)
analysis_data$race_label <- factor(
analysis_data$race,
levels = c(
"mex_american",
"other_hispanic",
"nh_white",
"nh_black",
"nh_asian"
),
labels = c(
"Mexican American",
"Other Hispanic",
"Non-Hispanic White",
"Non-Hispanic Black",
"Non-Hispanic Asian"
)
)
analysis_data$educ_label <- factor(
analysis_data$educ_lvl,
levels = c(
"lt_9th_grade",
"9_to_12th_grade_no_diploma",
"hs_grad_or_ged",
"start_college_to_aa",
"college_grad_or_above"
),
labels = c(
"Less than 9th grade",
"9th-12th, no diploma",
"High school or GED",
"Some college or AA",
"College graduate or above"
)
)
nrow(analysis_data)
## [1] 6299
sample_summary <- data.frame(
Quantity = c(
"Number of people",
"Number of smokers",
"Number of non-smokers",
"Number with CVD"
),
Value = c(
nrow(analysis_data),
sum(analysis_data$smoker_indicator == 1),
sum(analysis_data$smoker_indicator == 0),
sum(analysis_data$cvd_indicator == 1)
)
)
knitr::kable(sample_summary)
| Quantity | Value |
|---|---|
| Number of people | 6299 |
| Number of smokers | 2559 |
| Number of non-smokers | 3740 |
| Number with CVD | 646 |
We use the same three methods as before:
weighted_mean <- function(outcome, weights) {
sum(weights * outcome) / sum(weights)
}
fit_sensitivity_methods <- function(input_data) {
treatment_model <- glm(
smoker_indicator ~ age_yr + gender_label + race_label +
educ_label + inc_to_pov_ratio + bmi,
data = input_data,
family = binomial()
)
e_hat <- predict(treatment_model, type = "response")
outcome_model <- glm(
cvd_indicator ~ smoker_indicator + age_yr + gender_label +
race_label + educ_label + inc_to_pov_ratio + bmi,
data = input_data,
family = binomial()
)
data_if_non_smoker <- input_data
data_if_smoker <- input_data
data_if_non_smoker$smoker_indicator <- 0
data_if_smoker$smoker_indicator <- 1
m0_hat <- predict(
outcome_model,
newdata = data_if_non_smoker,
type = "response"
)
m1_hat <- predict(
outcome_model,
newdata = data_if_smoker,
type = "response"
)
smoker_rows <- input_data$smoker_indicator == 1
non_smoker_rows <- input_data$smoker_indicator == 0
ipw_risk_smoker <- weighted_mean(
input_data$cvd_indicator[smoker_rows],
1 / e_hat[smoker_rows]
)
ipw_risk_non_smoker <- weighted_mean(
input_data$cvd_indicator[non_smoker_rows],
1 / (1 - e_hat[non_smoker_rows])
)
or_risk_smoker <- mean(m1_hat)
or_risk_non_smoker <- mean(m0_hat)
T_obs <- input_data$smoker_indicator
Y_obs <- input_data$cvd_indicator
dr_risk_smoker <- mean(m1_hat + T_obs / e_hat * (Y_obs - m1_hat))
dr_risk_non_smoker <- mean(
m0_hat + (1 - T_obs) / (1 - e_hat) * (Y_obs - m0_hat)
)
c(
IPW_risk0 = ipw_risk_non_smoker,
IPW_risk1 = ipw_risk_smoker,
IPW_difference = ipw_risk_smoker - ipw_risk_non_smoker,
IPW_ratio = ipw_risk_smoker / ipw_risk_non_smoker,
Outcome_regression_risk0 = or_risk_non_smoker,
Outcome_regression_risk1 = or_risk_smoker,
Outcome_regression_difference = or_risk_smoker - or_risk_non_smoker,
Outcome_regression_ratio = or_risk_smoker / or_risk_non_smoker,
Doubly_robust_risk0 = dr_risk_non_smoker,
Doubly_robust_risk1 = dr_risk_smoker,
Doubly_robust_difference = dr_risk_smoker - dr_risk_non_smoker,
Doubly_robust_ratio = dr_risk_smoker / dr_risk_non_smoker
)
}
safe_fit_sensitivity_methods <- function(input_data) {
tryCatch(
fit_sensitivity_methods(input_data),
error = function(e) {
rep(NA_real_, 12)
}
)
}
point_results <- fit_sensitivity_methods(analysis_data)
methods <- c("IPW", "Outcome_regression", "Doubly_robust")
method_labels <- c("IPW", "Outcome regression", "Doubly robust")
point_risk0 <- point_results[paste0(methods, "_risk0")]
point_risk1 <- point_results[paste0(methods, "_risk1")]
point_difference <- point_results[paste0(methods, "_difference")]
point_ratio <- point_results[paste0(methods, "_ratio")]
names(point_risk0) <- methods
names(point_risk1) <- methods
names(point_difference) <- methods
names(point_ratio) <- methods
point_table <- data.frame(
Method = method_labels,
Non_smoker_risk_percent = round(100 * point_risk0, 1),
Smoker_risk_percent = round(100 * point_risk1, 1),
Difference_pp = round(100 * point_difference, 2),
Risk_ratio = round(point_ratio, 2)
)
knitr::kable(
point_table,
col.names = c(
"Method",
"Non-smoker risk (%)",
"Smoker risk (%)",
"Difference (pp)",
"Risk ratio"
)
)
| Method | Non-smoker risk (%) | Smoker risk (%) | Difference (pp) | Risk ratio | |
|---|---|---|---|---|---|
| IPW | IPW | 8.7 | 11.7 | 2.99 | 1.34 |
| Outcome_regression | Outcome regression | 8.5 | 12.0 | 3.45 | 1.40 |
| Doubly_robust | Doubly robust | 8.6 | 11.8 | 3.25 | 1.38 |
We still report a 95% confidence interval for the risk ratio, using the same bootstrap idea from Lecture 5.
For the E-value itself, we will keep the lecture simple and use only the point-estimate E-value.
set.seed(20260627)
B <- 300
n <- nrow(analysis_data)
bootstrap_ratios <- matrix(
NA_real_,
nrow = B,
ncol = length(methods)
)
colnames(bootstrap_ratios) <- methods
ratio_names <- paste0(methods, "_ratio")
for (b in seq_len(B)) {
bootstrap_rows <- sample.int(n, size = n, replace = TRUE)
bootstrap_data <- analysis_data[bootstrap_rows, ]
bootstrap_results <- safe_fit_sensitivity_methods(bootstrap_data)
bootstrap_ratios[b, ] <- bootstrap_results[ratio_names]
}
bootstrap_ratios <- bootstrap_ratios[complete.cases(bootstrap_ratios), ]
number_of_successful_bootstraps <- nrow(bootstrap_ratios)
number_of_successful_bootstraps
## [1] 300
ratio_ci <- t(apply(
bootstrap_ratios,
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
colnames(ratio_ci) <- c("Lower", "Upper")
log_ratio_se <- (log(ratio_ci[, "Upper"]) - log(ratio_ci[, "Lower"])) /
(2 * qnorm(0.975))
ratio_z <- log(point_ratio) / log_ratio_se
ratio_p_value <- 2 * pnorm(-abs(ratio_z))
names(ratio_p_value) <- methods
format_p_value <- function(p) {
ifelse(p < 0.001, "<0.001", sprintf("%.3f", p))
}
ratio_table <- data.frame(
Method = method_labels,
Risk_ratio = round(point_ratio, 2),
CI_lower = round(ratio_ci[, "Lower"], 2),
CI_upper = round(ratio_ci[, "Upper"], 2),
P_value = format_p_value(ratio_p_value)
)
knitr::kable(
ratio_table,
col.names = c(
"Method",
"Risk ratio",
"95% CI lower",
"95% CI upper",
"p-value"
)
)
| Method | Risk ratio | 95% CI lower | 95% CI upper | p-value | |
|---|---|---|---|---|---|
| IPW | IPW | 1.34 | 1.17 | 1.58 | <0.001 |
| Outcome_regression | Outcome regression | 1.40 | 1.22 | 1.64 | <0.001 |
| Doubly_robust | Doubly robust | 1.38 | 1.21 | 1.63 | <0.001 |
The p-value tests the null hypothesis:
H0: adjusted risk ratio = 1
Here, 1 means no difference in CVD risk between smokers
and non-smokers after adjustment.
We calculate the p-values on the log risk-ratio scale, using the bootstrap confidence interval to estimate the standard error.
e_value <- function(rr) {
rr_away_from_null <- ifelse(rr < 1, 1 / rr, rr)
ifelse(
rr_away_from_null <= 1,
1,
rr_away_from_null +
sqrt(rr_away_from_null * (rr_away_from_null - 1))
)
}
point_evalues <- e_value(point_ratio)
names(point_evalues) <- methods
evalue_table <- data.frame(
Method = method_labels,
Risk_ratio = round(point_ratio, 2),
RR_95_CI = paste0(
"(",
round(ratio_ci[, "Lower"], 2),
", ",
round(ratio_ci[, "Upper"], 2),
")"
),
P_value = format_p_value(ratio_p_value),
E_value = round(point_evalues, 2)
)
knitr::kable(
evalue_table,
col.names = c(
"Method",
"Risk ratio",
"95% CI",
"p-value",
"E-value"
)
)
| Method | Risk ratio | 95% CI | p-value | E-value | |
|---|---|---|---|---|---|
| IPW | IPW | 1.34 | (1.17, 1.58) | <0.001 | 2.02 |
| Outcome_regression | Outcome regression | 1.40 | (1.22, 1.64) | <0.001 | 2.16 |
| Doubly_robust | Doubly robust | 1.38 | (1.21, 1.63) | <0.001 | 2.10 |
In this lecture, the E-value is calculated from the point estimate of the risk ratio.
The p-value and the E-value answer different questions:
| Number | Question it answers |
|---|---|
| p-value | Is the adjusted association statistically different from no association? |
| E-value | How strong would hidden confounding need to be to explain away the adjusted association? |
In our table, all three p-values are <0.001. This
means the adjusted smoker versus non-smoker risk ratios are
statistically different from 1 in this dataset. However, a small p-value
does not prove causation; that is why we also use the E-value.
The formula is:
E-value = RR + sqrt(RR * (RR - 1))
display_rr <- round(point_ratio, 2)
evalue_calculation_table <- data.frame(
Method = method_labels,
Risk_ratio = display_rr,
Calculation = paste0(
display_rr,
" + sqrt(",
display_rr,
" * (",
display_rr,
" - 1))"
),
E_value = round(point_evalues, 2)
)
knitr::kable(
evalue_calculation_table,
col.names = c("Method", "Risk ratio", "Calculation", "E-value")
)
| Method | Risk ratio | Calculation | E-value | |
|---|---|---|---|---|
| IPW | IPW | 1.34 | 1.34 + sqrt(1.34 * (1.34 - 1)) | 2.02 |
| Outcome_regression | Outcome regression | 1.40 | 1.4 + sqrt(1.4 * (1.4 - 1)) | 2.16 |
| Doubly_robust | Doubly robust | 1.38 | 1.38 + sqrt(1.38 * (1.38 - 1)) | 2.10 |
The table uses the rounded risk ratios for display. R calculates using the exact risk ratios and then rounds the final E-value.
For example, for IPW:
RR = 1.34
E-value = 1.34 + sqrt(1.34 * (1.34 - 1))
= 2.02
So the IPW E-value of 2.02 means:
To fully explain away the IPW point estimate, a hidden confounder would need to be associated with both smoking and CVD by about 2.02-fold each, after adjusting for the measured covariates.
The E-value is a mathematical function of the risk ratio.
When the risk ratio is close to 1, the E-value is close to 1.
As the risk ratio moves farther away from 1, the E-value becomes larger.
This curve shows why a larger risk ratio leads to a larger E-value. Our three methods sit close together because their adjusted risk ratios are close together.
The dotted vertical line at 2 is a visual reference. It helps us see whether the E-value is around 2-fold, below 2-fold, or above 2-fold.
For our classroom dataset, the adjusted risk ratios are around 1.34 to 1.4.
The E-values for the point estimates are around 2.02 to 2.16.
That means an unmeasured confounder would need to be associated with both smoking and CVD by about a 2-fold to 2.2-fold risk ratio, above and beyond the covariates we already adjusted for, to fully explain away the point estimates.
VanderWeele and Ding often write the interpretation as one complete sentence.
The sentence pattern is:
The observed risk ratio of [RR] could be explained away by an unmeasured
confounder that was associated with both the treatment and the outcome by
a risk ratio of [E-value]-fold each, above and beyond the measured covariates.
But if both hidden-confounder links were weaker than this, they could not fully
explain away the point estimate.
Applied to our smoking and CVD example:
IPW. The observed risk ratio of 1.34 could be explained away by an unmeasured confounder that was associated with both smoking and CVD by a risk ratio of 2.02-fold each, above and beyond the measured covariates. If both hidden-confounder links were weaker than 2.02-fold, they could not fully explain away this point estimate.
Outcome regression. The observed risk ratio of 1.40 could be explained away by an unmeasured confounder that was associated with both smoking and CVD by a risk ratio of 2.16-fold each, above and beyond the measured covariates. If both hidden-confounder links were weaker than 2.16-fold, they could not fully explain away this point estimate.
Doubly robust. The observed risk ratio of 1.38 could be explained away by an unmeasured confounder that was associated with both smoking and CVD by a risk ratio of 2.10-fold each, above and beyond the measured covariates. If both hidden-confounder links were weaker than 2.10-fold, they could not fully explain away this point estimate.
This wording is useful because it forces us to say exactly how strong the hidden confounder would need to be.
We compare the E-value with a realistic guess about a possible hidden confounder.
For example, suppose we worry about a missing variable such as physical activity.
So in practice:
If both links could realistically be around 2-fold, then the result may be sensitive to that hidden confounder. If both links are probably much weaker than 2-fold, then the result is more reassuring.
This is why the E-value is useful: it changes a vague worry into a concrete question.
Could a missing variable be this strongly related to both smoking and CVD?
Suppose there is a hidden variable U.
This is closer to Peng Ding’s figure: U sits above
treatment and outcome because it is the hidden variable that may create
extra association between them.
If the E-value is about 2, then a hidden confounder would need to be
related to both smoking and CVD by about 2-fold each, after adjustment
for X, to explain away the point estimate.
This does not mean such a confounder exists. It gives us a way to judge how worried we should be.
One more subtle point from the paper:
The E-value gives the minimum common strength if the two hidden-confounder links are equally strong.
If one link is weaker than the E-value, the other link would need to be even stronger to explain away the result.
E-values are a good first sensitivity analysis because they are simple and visual.
Other useful sensitivity checks include:
For this lecture, the E-value is the best main topic because it directly answers:
How strong would hidden confounding need to be?
The E-value is helpful, but it is not magic.
It does not fix:
The E-value only asks about one type of concern:
possible unmeasured confounding.
A large E-value does not prove causation. A small E-value does not prove there is no effect. The E-value simply tells us how easy or hard it would be for hidden confounding to explain away the adjusted association.
This appendix is optional. It gives the algebra behind the formula.
Suppose the observed adjusted risk ratio is greater than 1:
\[ RR_{\text{obs}} > 1. \]
For us, \(RR_{\text{obs}}\) is the adjusted smoker-versus-non-smoker CVD risk ratio.
Now imagine there is an unmeasured confounder \(U\).
VanderWeele and Ding describe two hidden-confounder strengths:
\[ \begin{aligned} RR_{EU} &= \text{strength of association between exposure and } U, \\ RR_{UD} &= \text{strength of association between } U \text{ and outcome}. \end{aligned} \]
In our lecture, \(RR_{EU}\) means how strongly \(U\) is related to smoking, and \(RR_{UD}\) means how strongly \(U\) is related to CVD. These are considered after the measured covariates \(X\) have already been adjusted for.
The risk-ratio bias factor is:
\[ B = \frac{RR_{EU} \times RR_{UD}} {RR_{EU} + RR_{UD} - 1}. \]
Think of \(B\) as the largest amount by which hidden confounding could multiply the risk ratio. If there is no hidden confounding, \(B = 1\). If hidden confounding is strong, \(B > 1\).
The adjusted risk ratio could be explained away if hidden confounding is strong enough to move the observed risk ratio down to 1:
\[ \frac{RR_{\text{obs}}}{B} \leq 1. \]
This is equivalent to:
\[ B \geq RR_{\text{obs}}. \]
So the question becomes:
How large do \(RR_{EU}\) and \(RR_{UD}\) need to be so that the bias factor is at least as large as the observed risk ratio?
The E-value gives a simple benchmark by setting the two hidden-confounder strengths equal:
\[ RR_{EU} = RR_{UD} = r. \]
Then the bias factor becomes:
\[ B = \frac{r \times r}{r + r - 1} = \frac{r^2}{2r - 1}. \]
To find the smallest equal strength that could explain away the result, set the bias factor equal to the observed risk ratio:
\[ \frac{r^2}{2r - 1} = RR_{\text{obs}}. \]
Now solve for \(r\):
\[ \begin{aligned} r^2 &= RR_{\text{obs}}(2r - 1), \\ r^2 &= 2RR_{\text{obs}}r - RR_{\text{obs}}, \\ r^2 - 2RR_{\text{obs}}r + RR_{\text{obs}} &= 0. \end{aligned} \]
Using the quadratic formula:
\[ r = \frac{2RR_{\text{obs}} \pm \sqrt{(2RR_{\text{obs}})^2 - 4RR_{\text{obs}}}}{2}. \]
Simplify:
\[ \begin{aligned} r &= RR_{\text{obs}} \pm \sqrt{RR_{\text{obs}}^2 - RR_{\text{obs}}} \\ &= RR_{\text{obs}} \pm \sqrt{RR_{\text{obs}}(RR_{\text{obs}} - 1)}. \end{aligned} \]
The smaller solution is below 1, but a risk-ratio strength is defined away from the null and must be at least 1. So we use the larger solution:
\[ \boxed{ \text{E-value} = RR_{\text{obs}} + \sqrt{RR_{\text{obs}}(RR_{\text{obs}} - 1)} } \]
That is the E-value formula.
For example, using the IPW risk ratio:
\[ \begin{aligned} RR_{\text{obs}} &= 1.34, \\ \text{E-value} &= 1.34 + \sqrt{1.34(1.34 - 1)} \\ &= 2.02. \end{aligned} \]
For a protective association where \(RR_{\text{obs}} < 1\), we first invert the risk ratio:
\[ RR_{\text{away from null}} = \frac{1}{RR_{\text{obs}}}. \]
Then we apply the same formula to \(RR_{\text{away from null}}\).
This appendix shows the math behind the three adjusted estimators 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} \]
For IPW and doubly robust estimation, we first estimate the probability of being a smoker given the baseline covariates:
\[ e(X_i) = P(T_i = 1 \mid X_i). \]
In this lecture, we estimate \(e(X_i)\) using logistic regression:
\[ \text{logit}\{e(X_i)\} = \alpha_0 + \alpha^\top X_i. \]
The fitted value is:
\[ \hat e_i = \hat e(X_i). \]
IPW means inverse probability weighting.
The idea is:
Give more weight to people whose observed smoking status was less likely given their covariates.
For smokers, the weight is:
\[ w_i^{(1)} = \frac{T_i}{\hat e_i}. \]
For non-smokers, the weight is:
\[ w_i^{(0)} = \frac{1 - T_i}{1 - \hat e_i}. \]
The normalized IPW estimator of the smoker risk is:
\[ \hat\psi_1^{IPW} = \frac{ \sum_{i=1}^n \dfrac{T_iY_i}{\hat e_i} }{ \sum_{i=1}^n \dfrac{T_i}{\hat e_i} }. \]
The normalized IPW estimator of the non-smoker risk is:
\[ \hat\psi_0^{IPW} = \frac{ \sum_{i=1}^n \dfrac{(1 - T_i)Y_i}{1 - \hat e_i} }{ \sum_{i=1}^n \dfrac{1 - T_i}{1 - \hat e_i} }. \]
Then:
\[ \begin{aligned} \widehat{RD}^{IPW} &= \hat\psi_1^{IPW} - \hat\psi_0^{IPW}, \\ \widehat{RR}^{IPW} &= \frac{\hat\psi_1^{IPW}}{\hat\psi_0^{IPW}}. \end{aligned} \]
Outcome regression models the outcome directly.
Let:
\[ m(t, X_i) = P(Y_i = 1 \mid T_i = t, X_i). \]
In this lecture, we estimate this using logistic regression:
\[ \text{logit}\{m(T_i, X_i)\} = \beta_0 + \beta_1T_i + \beta^\top X_i. \]
After fitting the model, we predict each person’s CVD risk twice:
\[ \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.
The outcome regression estimator is:
\[ \hat\psi_1^{OR} = \frac{1}{n} \sum_{i=1}^n \hat m_1(X_i), \]
and:
\[ \hat\psi_0^{OR} = \frac{1}{n} \sum_{i=1}^n \hat m_0(X_i). \]
Then:
\[ \begin{aligned} \widehat{RD}^{OR} &= \hat\psi_1^{OR} - \hat\psi_0^{OR}, \\ \widehat{RR}^{OR} &= \frac{\hat\psi_1^{OR}}{\hat\psi_0^{OR}}. \end{aligned} \]
Doubly robust estimation combines the two ideas:
For the smoker risk:
\[ \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]. \]
For the non-smoker risk:
\[ \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]. \]
The first part, \(\hat m_t(X_i)\), is the outcome model prediction.
The second part is a correction term:
\[ \frac{\mathbf{1}(T_i = t)}{\hat P(T_i = t \mid X_i)} \{Y_i - \hat m_t(X_i)\}. \]
This correction uses the observed data to adjust the prediction when the model is off.
Then:
\[ \begin{aligned} \widehat{RD}^{DR} &= \hat\psi_1^{DR} - \hat\psi_0^{DR}, \\ \widehat{RR}^{DR} &= \frac{\hat\psi_1^{DR}}{\hat\psi_0^{DR}}. \end{aligned} \]
The doubly robust estimator has a useful property:
It can still be consistent if either the treatment model is correct or the outcome model is correct.
In simple language:
This does not mean it is magic. If both models are badly wrong, the doubly robust estimator can still be biased.
VanderWeele, T. J., & Ding, P. (2017). Sensitivity analysis in observational research: Introducing the E-value. Annals of Internal Medicine, 167(4), 268-274. https://doi.org/10.7326/M16-2607
Ding, P., & VanderWeele, T. J. (2016). Sensitivity analysis without assumptions. Epidemiology, 27(3), 368-377.