Semiparametric Efficiency Lower Bound for the Smoking-CVD Risk Difference
Se Yoon Lee
2026-07-06
1 Today’s Question
Lectures 1-6 estimated and interpreted the causal risk difference for cardiovascular disease (CVD) under smoking versus non-smoking.
This note asks the PhD-level version of the same classroom problem:
Among all regular estimators of the identified risk difference, what is the smallest possible large-sample variance?
That smallest possible variance is the semiparametric efficiency lower bound.
The target parameter is still
\[ \psi = E\{Y(1)\} - E\{Y(0)\}, \]
where \(Y(1)\) is the potential CVD outcome under smoking and \(Y(0)\) is the potential CVD outcome under non-smoking.
The bridge from Lecture 4 to this note is the doubly robust estimator. Lecture 4 introduced it as a prediction-plus-correction estimator. Here we show that the same correction is not a trick; it is the empirical version of the efficient influence function.
We will again use the classroom NHANES complete-case dataset from Lectures 1-6. The efficiency bound derived below is the bound for the i.i.d. classroom analysis used throughout these notes.
2 Source Map
The attached material points to one coherent story.
| Source material | Role in this note |
|---|---|
| Kennedy, Semiparametric Theory and Empirical Processes in Causal Inference | Tangent spaces, efficient influence functions, ATE EIF, product-rate remainder |
| Wager, Causal Inference | AIPW intuition, strong double robustness, cross-fitting perspective |
| On the Asymptotic Properties of Debiased Machine Learning Estimators | Neyman orthogonality and cross-fitting details |
| debiased_ml_slides | Pedagogical version of orthogonal scores and sample splitting |
The note does not quote these sources at length. Instead, it adapts their machinery to the same smoking-CVD risk difference used throughout this lecture course.
3 Data Structure and Identification
For each person, let
\[ O = (X,T,Y), \]
where:
| Symbol | Meaning in this course |
|---|---|
| \(T=1\) | smoker |
| \(T=0\) | non-smoker |
| \(Y=1\) | had CVD before interview |
| \(Y=0\) | did not have CVD before interview |
| \(X\) | age, gender, race/ethnicity, education, income-to-poverty ratio, BMI |
The causal estimand is
\[ \psi = \psi_1 - \psi_0, \qquad \psi_t = E\{Y(t)\}. \]
The observed-data identification assumptions are the same ones used earlier:
- Consistency: if \(T=t\), then \(Y=Y(t)\).
- Conditional exchangeability: \(\{Y(1),Y(0)\} \perp T \mid X\).
- Positivity: \(0 < P(T=1 \mid X) < 1\) for covariate values in the target population.
Under these assumptions,
\[ \psi_t = E\{E(Y \mid T=t,X)\}. \]
Define the nuisance functions
\[ e(X) = P(T=1 \mid X), \]
and
\[ Q_t(X) = E(Y \mid T=t,X), \qquad t \in \{0,1\}. \]
Then the identified parameter is
\[ \psi = E\{Q_1(X)-Q_0(X)\}. \]
4 The Semiparametric Model
A parametric model specifies a finite-dimensional family \(p(o;\theta)\). A semiparametric model lets part of the distribution be infinite-dimensional.
For this note, the main model is the nonparametric observed-data model:
\[ P \in \mathcal{M}_{np}, \]
where the distribution of \(X\), the treatment mechanism \(e(X)\), and the outcome regression \(Q_t(X)\) are unrestricted except for positivity and regularity conditions.
This matters because the target parameter \(\psi(P)\) is one-dimensional, but \(P\) contains many nuisance features:
\[ P \quad \Longleftrightarrow \quad \{p_X(x), e(x), Q_0(x), Q_1(x), \text{conditional outcome variation}\}. \]
The semiparametric question is:
If we refuse to fully parametrize these nuisance parts, what is the best possible first-order variance for estimating \(\psi\)?
5 Influence Functions
An estimator \(\hat\psi\) is asymptotically linear if
\[ \hat\psi - \psi = \frac{1}{n}\sum_{i=1}^n \phi(O_i) + o_p(n^{-1/2}), \]
where \(E\{\phi(O)\}=0\) and \(E\{\phi(O)^2\}<\infty\).
Then
\[ \sqrt{n}(\hat\psi-\psi) \rightsquigarrow N\{0, E(\phi^2)\}. \]
The function \(\phi\) is an influence function for the estimator. If \(\phi\) is the efficient influence function, its variance is the smallest possible variance among regular asymptotically linear estimators in the model.
6 Tangent-Space Geometry
Semiparametric efficiency is geometry in
\[ L_0^2(P) = \{f(O):E_P f(O)=0,\ E_P f(O)^2<\infty\}, \]
with inner product
\[ \langle f,h\rangle_P = E_P\{f(O)h(O)\}. \]
The key objects are:
| Object | Notation | Meaning |
|---|---|---|
| Regular submodel | \(\{P_\epsilon:\epsilon\in(-\delta,\delta), P_0=P\}\) | A one-dimensional path through the statistical model \(\mathcal M\). |
| Score direction | \(S(O)=\left.\frac{d}{d\epsilon}\log p_\epsilon(O)\right|_{\epsilon=0}\) | The velocity of the submodel at \(P\), viewed as a vector in \(L_0^2(P)\). |
| Tangent space | \(\mathcal T(P)=\overline{\{S: S \text{ is a regular score at }P\}}\) | All locally possible directions through \(\mathcal M\). |
| Pathwise derivative | \(\dot\psi_P(S)=\left.\frac{d}{d\epsilon}\psi(P_\epsilon)\right|_{\epsilon=0}\) | Directional derivative of the target along score \(S\). |
| Gradient | \(g_\psi\) satisfying \(\dot\psi_P(S)=\langle g_\psi,S\rangle_P\) | A Riesz representer of the derivative. |
| Nuisance tangent space | \(\mathcal T_\eta=\{S\in\mathcal T(P):\dot\psi_P(S)=0\}\) | Directions that change nuisance features but not the target. |
| Canonical gradient / EIF | \(D_P^*\) or \(\phi_{\mathrm{EIF}}\) | The gradient after removing all nuisance variation. |
The projection formula is the central geometric statement:
\[ D_P^* = \Pi_{\mathcal T_\eta^\perp}g_\psi = g_\psi-\Pi_{\mathcal T_\eta}g_\psi, \qquad \langle D_P^*,S_\eta\rangle_P=0 \quad\text{for all }S_\eta\in\mathcal T_\eta. \]
Thus the efficient influence function is not just any influence function. It is the canonical gradient, the part of the target gradient that remains after projecting away directions irrelevant to \(\psi(P)\).
For the ATE/risk-difference problem, the observed data are \(O=(X,T,Y)\) and the likelihood factorizes as
\[ p(o)=p(x)\,p(t\mid x)\,p(y\mid t,x). \]
This induces the orthogonal tangent-space decomposition
\[ \mathcal T(P) = \mathcal T_X \oplus \mathcal T_{T\mid X} \oplus \mathcal T_{Y\mid T,X}. \]
The ATE target is
\[ \psi(P)=E_P\{Q_1(X)-Q_0(X)\}, \qquad Q_t(X)=E_P(Y\mid T=t,X), \qquad e(X)=P(T=1\mid X). \]
The efficient influence function decomposes into interpretable orthogonal pieces:
\[ D_P^*(O) = \underbrace{Q_1(X)-Q_0(X)-\psi}_{\phi_X(O)} + \underbrace{0}_{\phi_{T\mid X}(O)} + \underbrace{H(T,X)\{Y-Q_T(X)\}}_{\phi_{Y\mid T,X}(O)}, \]
where
\[ H(T,X) = \frac{T}{e(X)} - \frac{1-T}{1-e(X)} \]
is the residual-weight function multiplying the outcome residual. This is the same augmentation term that turns outcome standardization into the doubly robust estimator. In estimated form, the DR estimator can be written as the solution of the empirical EIF equation
\[ P_nD_{\hat P,\psi}^*(O)=0, \]
where
\[ D_{\hat P,\psi}^*(O) = \hat Q_1(X)-\hat Q_0(X) + \hat H(T,X)\{Y-\hat Q_T(X)\} - \psi. \]
This geometric view is why the same phrase, orthogonality, appears in semiparametric efficiency, doubly robust estimation, and debiased machine learning.
7 Efficient Influence Function for \(\psi_t\)
For
\[ \psi_t = E\{Q_t(X)\}, \]
the efficient influence function in the nonparametric observed-data model is
\[ D_t^*(O) = \frac{I(T=t)}{P(T=t \mid X)} \{Y-Q_t(X)\} + Q_t(X)-\psi_t. \]
For \(t=1\),
\[ D_1^*(O) = \frac{T}{e(X)}\{Y-Q_1(X)\} + Q_1(X)-\psi_1. \]
For \(t=0\),
\[ D_0^*(O) = \frac{1-T}{1-e(X)}\{Y-Q_0(X)\} + Q_0(X)-\psi_0. \]
Each \(D_t^*\) has two components:
| Component | Formula | Meaning |
|---|---|---|
| Residual correction | \(I(T=t)\{Y-Q_t(X)\}/P(T=t\mid X)\) | uses observed residuals from people actually in treatment arm \(t\) |
| Standardization variation | \(Q_t(X)-\psi_t\) | variation in predicted counterfactual risk across covariates |
8 Efficient Influence Function for the Risk Difference
The risk difference is
\[ \psi = \psi_1-\psi_0. \]
Therefore its efficient influence function is
\[ D^*(O) = D_1^*(O)-D_0^*(O). \]
Equivalently,
\[ D^*(O) = \frac{T}{e(X)} \{Y-Q_1(X)\} - \frac{1-T}{1-e(X)} \{Y-Q_0(X)\} + Q_1(X)-Q_0(X)-\psi. \]
This formula is the mathematical heart of Lecture 4.
The doubly robust estimator is obtained by replacing unknown nuisance functions with estimates and averaging the non-\(-\psi\) part:
\[ \hat\psi_{DR} = \frac{1}{n} \sum_{i=1}^{n} \left[ \hat Q_1(X_i)-\hat Q_0(X_i) + \frac{T_i}{\hat e(X_i)} \{Y_i-\hat Q_1(X_i)\} - \frac{1-T_i}{1-\hat e(X_i)} \{Y_i-\hat Q_0(X_i)\} \right]. \]
The empirical mean of the estimated influence function is exactly zero if we define
\[ \hat D_i = \left[ \hat Q_1(X_i)-\hat Q_0(X_i) + \frac{T_i}{\hat e(X_i)} \{Y_i-\hat Q_1(X_i)\} - \frac{1-T_i}{1-\hat e(X_i)} \{Y_i-\hat Q_0(X_i)\} \right] - \hat\psi_{DR}. \]
9 The Efficiency Lower Bound
The semiparametric efficiency lower bound for estimating \(\psi\) is
\[ \mathcal{I}_{eff}^{-1} = E\{D^*(O)^2\}. \]
Thus any regular estimator \(\hat\psi\) must satisfy, asymptotically,
\[ \mathrm{Var}(\hat\psi) \geq \frac{E\{D^*(O)^2\}}{n} \]
in the nonparametric i.i.d. model.
If an estimator is asymptotically linear with influence function \(D^*\), then it attains the bound:
\[ \sqrt{n}(\hat\psi-\psi) \rightsquigarrow N\left(0, E\{D^*(O)^2\}\right). \]
This is the formal meaning of “efficient” in semiparametric theory.
10 Variance Decomposition
For the risk difference, the efficient variance has a useful decomposition:
\[ E\{D^*(O)^2\} = \mathrm{Var}\{Q_1(X)-Q_0(X)\} + E\left[ \frac{\mathrm{Var}(Y \mid T=1,X)}{e(X)} + \frac{\mathrm{Var}(Y \mid T=0,X)}{1-e(X)} \right]. \]
For binary \(Y\),
\[ \mathrm{Var}(Y \mid T=t,X) = Q_t(X)\{1-Q_t(X)\}. \]
So
\[ E\{D^*(O)^2\} = \mathrm{Var}\{Q_1(X)-Q_0(X)\} + E\left[ \frac{Q_1(X)\{1-Q_1(X)\}}{e(X)} + \frac{Q_0(X)\{1-Q_0(X)\}}{1-e(X)} \right]. \]
This decomposition is extremely informative:
- The first term is heterogeneity in conditional causal risk differences.
- The second term is residual outcome noise under the two treatment worlds.
- Positivity appears directly through \(1/e(X)\) and \(1/\{1-e(X)\}\).
11 Orthogonality and the Product Remainder
The efficient influence function is orthogonal to nuisance directions. In modern DML language, the estimating equation is Neyman orthogonal.
For the DR estimating function, nuisance error enters the second-order remainder as a product. Let \(\hat e\), \(\hat Q_1\), and \(\hat Q_0\) be estimated nuisance functions. A key remainder term is of the form
\[ R_2(\hat\eta,\eta) \approx E\left[ \{\hat e(X)-e(X)\} \left\{ \frac{\hat Q_1(X)-Q_1(X)}{\hat e(X)} + \frac{\hat Q_0(X)-Q_0(X)}{1-\hat e(X)} \right\} \right]. \]
The exact sign convention depends on how the remainder is written, but the important point is the product:
\[ |R_2(\hat\eta,\eta)| \lesssim \|\hat e-e\| \left( \|\hat Q_1-Q_1\|+\|\hat Q_0-Q_0\| \right) \]
under positivity.
This explains the \(n^{-1/4}\) rule:
\[ \|\hat e-e\| = o_p(n^{-1/4}) \quad \text{and} \quad \|\hat Q_t-Q_t\| = o_p(n^{-1/4}) \]
are enough to make the product
\[ o_p(n^{-1/2}). \]
That lets the first-order behavior be governed by the empirical average of \(D^*(O)\), yielding asymptotic normality and efficiency.
12 Double Robustness vs Efficiency
Double robustness and efficiency are related but distinct.
| Property | Meaning |
|---|---|
| Double robustness | consistency can hold if either \(e(X)\) or \(Q_t(X)\) is consistently estimated |
| Local efficiency | if all nuisance functions are consistently estimated well enough, the estimator has influence function \(D^*\) |
| Semiparametric efficiency | among regular estimators in the model, no estimator has smaller asymptotic variance than \(E(D^{*2})/n\) |
So DR is not merely “two chances to be right.” Its deeper role is:
It targets the efficient influence function and therefore can attain the efficiency lower bound.
13 Load and Prepare the NHANES Classroom Data
data_candidates <- c(
file.path("..", "Data", "NHANES_data.csv"),
file.path("Data", "NHANES_data.csv"),
"C:/Users/seyoo/OneDrive/GitHub_repository/Causal_Inference_Survival_Analysis/Causal_Inference_Lecture/Data/NHANES_data.csv"
)
data_path <- data_candidates[file.exists(data_candidates)][1]
if (is.na(data_path)) {
stop("Could not find NHANES_data.csv. Check the data path.")
}
nhanes <- read.csv(data_path, stringsAsFactors = FALSE)
baseline_vars <- c(
"age_yr",
"gender",
"race",
"educ_lvl",
"inc_to_pov_ratio",
"bmi"
)
analysis_vars <- c("smoker_indicator", "cvd_indicator", baseline_vars)
eff_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]
eff_data$smoker_indicator <- as.integer(eff_data$smoker_indicator)
eff_data$cvd_indicator <- as.integer(eff_data$cvd_indicator)
eff_data$smoking_group <- factor(
ifelse(eff_data$smoker_indicator == 1, "Smoker", "Non-smoker"),
levels = c("Non-smoker", "Smoker")
)
eff_data$gender_label <- factor(
eff_data$gender,
levels = c("F", "M"),
labels = c("Female", "Male")
)
eff_data$race_label <- factor(
eff_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"
)
)
eff_data$educ_label <- factor(
eff_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(eff_data)## [1] 6299sample_summary <- data.frame(
Group = c("Non-smoker", "Smoker", "Total"),
N = c(
sum(eff_data$smoker_indicator == 0),
sum(eff_data$smoker_indicator == 1),
nrow(eff_data)
),
CVD_percent = round(100 * c(
mean(eff_data$cvd_indicator[eff_data$smoker_indicator == 0]),
mean(eff_data$cvd_indicator[eff_data$smoker_indicator == 1]),
mean(eff_data$cvd_indicator)
), 2)
)
mkable(sample_summary)| Group | N | CVD_percent |
|---|---|---|
| Non-smoker | 3740 | 6.93 |
| Smoker | 2559 | 15.12 |
| Total | 6299 | 10.26 |
14 Estimate the Nuisance Functions
We use the same classroom covariate set as Lectures 2-6.
treatment_model <- glm(
smoker_indicator ~ age_yr + gender_label + race_label +
educ_label + inc_to_pov_ratio + bmi,
data = eff_data,
family = binomial()
)
outcome_model <- glm(
cvd_indicator ~ smoker_indicator + age_yr + gender_label +
race_label + educ_label + inc_to_pov_ratio + bmi,
data = eff_data,
family = binomial()
)
eff_data$e_hat <- bound_probability(predict(treatment_model, type = "response"))
data_if_smoker <- eff_data
data_if_non_smoker <- eff_data
data_if_smoker$smoker_indicator <- 1
data_if_non_smoker$smoker_indicator <- 0
eff_data$q1_hat <- bound_probability(
predict(outcome_model, newdata = data_if_smoker, type = "response")
)
eff_data$q0_hat <- bound_probability(
predict(outcome_model, newdata = data_if_non_smoker, type = "response")
)
nuisance_summary <- data.frame(
Quantity = c(
"Estimated propensity score \\(\\hat e(X)\\)",
"Predicted CVD risk if smoker \\(\\hat Q_1(X)\\)",
"Predicted CVD risk if non-smoker \\(\\hat Q_0(X)\\)"
),
Mean_SD = c(
mean_sd(eff_data$e_hat),
mean_sd(eff_data$q1_hat),
mean_sd(eff_data$q0_hat)
),
Min = round(c(
min(eff_data$e_hat),
min(eff_data$q1_hat),
min(eff_data$q0_hat)
), 4),
Max = round(c(
max(eff_data$e_hat),
max(eff_data$q1_hat),
max(eff_data$q0_hat)
), 4)
)
names(nuisance_summary) <- c("Quantity", "Mean (SD)", "Min", "Max")
mkable(nuisance_summary)| Quantity | Mean (SD) | Min | Max |
|---|---|---|---|
| Estimated propensity score \(\hat e(X)\) | 0.4063 (0.1769) | 0.0542 | 0.8597 |
| Predicted CVD risk if smoker \(\hat Q_1(X)\) | 0.1198 (0.1205) | 0.0031 | 0.6744 |
| Predicted CVD risk if non-smoker \(\hat Q_0(X)\) | 0.0853 (0.0915) | 0.0020 | 0.5737 |
15 Estimate the DR Effect and Its EIF
T_obs <- eff_data$smoker_indicator
Y_obs <- eff_data$cvd_indicator
e_hat <- eff_data$e_hat
q1_hat <- eff_data$q1_hat
q0_hat <- eff_data$q0_hat
dr_component <- q1_hat - q0_hat +
T_obs / e_hat * (Y_obs - q1_hat) -
(1 - T_obs) / (1 - e_hat) * (Y_obs - q0_hat)
psi_dr <- mean(dr_component)
eif_dr <- dr_component - psi_dr
or_component <- q1_hat - q0_hat
psi_or <- mean(or_component)
ipw_psi1 <- weighted_mean(Y_obs[T_obs == 1], 1 / e_hat[T_obs == 1])
ipw_psi0 <- weighted_mean(Y_obs[T_obs == 0], 1 / (1 - e_hat[T_obs == 0]))
psi_ipw <- ipw_psi1 - ipw_psi0
psi_crude <- mean(Y_obs[T_obs == 1]) - mean(Y_obs[T_obs == 0])
effect_table <- data.frame(
Estimator = c(
"\\(\\hat\\psi_{crude}\\)",
"\\(\\hat\\psi_{IPW}\\)",
"\\(\\hat\\psi_{OR}\\)",
"\\(\\hat\\psi_{DR}\\)"
),
Method = c("Observed comparison", "IPW", "Outcome regression", "Doubly robust"),
Risk_difference_pp = round(100 * c(
psi_crude,
psi_ipw,
psi_or,
psi_dr
), 3)
)
names(effect_table) <- c("Estimator", "Method", "Risk difference (pp)")
mkable(effect_table)| Estimator | Method | Risk difference (pp) |
|---|---|---|
| \(\hat\psi_{crude}\) | Observed comparison | 8.198 |
| \(\hat\psi_{IPW}\) | IPW | 2.987 |
| \(\hat\psi_{OR}\) | Outcome regression | 3.453 |
| \(\hat\psi_{DR}\) | Doubly robust | 3.253 |
The estimated efficient influence function for the DR estimator is centered by construction:
eif_check <- data.frame(
Quantity = c(
"Mean estimated EIF \\(P_n\\hat D\\)",
"SD estimated EIF \\(sd(\\hat D_i)\\)",
"Estimated EIF SE \\(sd(\\hat D_i)/\\sqrt n\\)",
"Estimated EIF SE in percentage points"
),
Value = c(
mean(eif_dr),
sd(eif_dr),
sd(eif_dr) / sqrt(nrow(eff_data)),
100 * sd(eif_dr) / sqrt(nrow(eff_data))
)
)
eif_check$Value <- signif(eif_check$Value, 5)
mkable(eif_check)| Quantity | Value |
|---|---|
| Mean estimated EIF \(P_n\hat D\) | 0.000000 |
| SD estimated EIF \(sd(\hat D_i)\) | 0.604530 |
| Estimated EIF SE \(sd(\hat D_i)/\sqrt n\) | 0.007617 |
| Estimated EIF SE in percentage points | 0.761700 |
The estimated semiparametric efficiency lower bound for \(\mathrm{Var}(\hat\psi)\) is
\[ \frac{\widehat{\mathrm{Var}}(\hat D)}{n}. \]
This is not a theorem that our fitted logistic models are correct. It is a plug-in estimate of the bound based on the nuisance functions used in this classroom analysis.
16 Decompose the Estimated Bound
The empirical analog of the bound decomposition is:
tau_hat_x <- q1_hat - q0_hat
heterogeneity_component <- var(tau_hat_x)
residual_component <- mean(
q1_hat * (1 - q1_hat) / e_hat +
q0_hat * (1 - q0_hat) / (1 - e_hat)
)
total_decomposition <- heterogeneity_component + residual_component
direct_eif_variance <- var(eif_dr)
bound_table <- data.frame(
Component = c(
"Heterogeneity \\(\\widehat{\\mathrm{Var}}\\{\\hat Q_1(X)-\\hat Q_0(X)\\}\\)",
"Residual noise \\(P_n[\\hat Q_1(1-\\hat Q_1)/\\hat e + \\hat Q_0(1-\\hat Q_0)/(1-\\hat e)]\\)",
"Total decomposition \\(\\hat V_{\\mathrm{EIF}}\\)",
"Direct sample variance \\(\\widehat{\\mathrm{Var}}(\\hat D_i)\\)"
),
Variance_scale = round(c(
heterogeneity_component,
residual_component,
total_decomposition,
direct_eif_variance
), 5),
SE_percentage_points = round(100 * sqrt(c(
heterogeneity_component,
residual_component,
total_decomposition,
direct_eif_variance
) / nrow(eff_data)), 3)
)
names(bound_table) <- c("Component", "Variance scale", "SE (pp)")
mkable(bound_table)| Component | Variance scale | SE (pp) |
|---|---|---|
| Heterogeneity \(\widehat{\mathrm{Var}}\{\hat Q_1(X)-\hat Q_0(X)\}\) | 0.00088 | 0.037 |
| Residual noise \(P_n[\hat Q_1(1-\hat Q_1)/\hat e + \hat Q_0(1-\hat Q_0)/(1-\hat e)]\) | 0.37734 | 0.774 |
| Total decomposition \(\hat V_{\mathrm{EIF}}\) | 0.37822 | 0.775 |
| Direct sample variance \(\widehat{\mathrm{Var}}(\hat D_i)\) | 0.36546 | 0.762 |
Interpretation. The estimated efficiency bound is almost entirely driven by the residual-noise component, not by treatment-effect heterogeneity. Numerically, the heterogeneity contribution is about \(0.0009\), while the residual contribution is about \(0.3773\). This means that, after conditioning on the observed covariates, most of the unavoidable uncertainty comes from the binary outcome variation \(Y-Q_T(X)\), amplified by the inverse propensity weights \(1/\hat e(X)\) and \(1/\{1-\hat e(X)\}\). In this example, improving precision is therefore less about explaining heterogeneity in \(\hat Q_1(X)-\hat Q_0(X)\), and more about the hard statistical fact that CVD is a noisy binary outcome and overlap controls how much that noise is inflated.
The total decomposition \(\hat V_{\mathrm{EIF}}=0.3782\) is also close to the direct empirical variance \(\widehat{\mathrm{Var}}(\hat D_i)=0.3655\). This agreement is a useful diagnostic: the closed-form decomposition and the sample variance of the estimated EIF are telling essentially the same first-order uncertainty story. On the standard-error scale, both imply an uncertainty of about \(0.76\)–\(0.78\) percentage points for the DR risk-difference estimate.
17 Bootstrap SE vs EIF SE
Lecture 5 used the bootstrap. Semiparametric theory gives an analytic first-order SE estimate through the EIF:
\[ \widehat{SE}_{EIF} = \sqrt{ \frac{1}{n} \widehat{\mathrm{Var}}(\hat D) }. \]
The bootstrap and EIF SE answer similar first-order questions. The bootstrap reruns the estimator; the EIF SE uses the estimated first-order linear representation.
fit_dr_once <- function(input_data) {
treatment_model_b <- glm(
smoker_indicator ~ age_yr + gender_label + race_label +
educ_label + inc_to_pov_ratio + bmi,
data = input_data,
family = binomial()
)
outcome_model_b <- glm(
cvd_indicator ~ smoker_indicator + age_yr + gender_label +
race_label + educ_label + inc_to_pov_ratio + bmi,
data = input_data,
family = binomial()
)
e_b <- bound_probability(predict(treatment_model_b, type = "response"))
d1 <- input_data
d0 <- input_data
d1$smoker_indicator <- 1
d0$smoker_indicator <- 0
q1_b <- bound_probability(predict(outcome_model_b, newdata = d1, type = "response"))
q0_b <- bound_probability(predict(outcome_model_b, newdata = d0, type = "response"))
t_b <- input_data$smoker_indicator
y_b <- input_data$cvd_indicator
mean(q1_b - q0_b +
t_b / e_b * (y_b - q1_b) -
(1 - t_b) / (1 - e_b) * (y_b - q0_b))
}
set.seed(20260706)
B_boot <- 200
n_eff <- nrow(eff_data)
boot_dr <- numeric(B_boot)
for (b in seq_len(B_boot)) {
rows_b <- sample.int(n_eff, size = n_eff, replace = TRUE)
boot_dr[b] <- tryCatch(
fit_dr_once(eff_data[rows_b, ]),
error = function(e) NA_real_
)
}
boot_dr <- boot_dr[!is.na(boot_dr)]
se_compare <- data.frame(
Method = c("\\(\\widehat{SE}_{EIF}\\)", "\\(\\widehat{SE}_{boot}\\)"),
SE_percentage_points = round(100 * c(
sd(eif_dr) / sqrt(n_eff),
sd(boot_dr)
), 3),
Number_of_bootstrap_samples = c(NA, length(boot_dr))
)
names(se_compare) <- c("SE estimator", "SE (pp)", "Bootstrap samples")
mkable(se_compare)| SE estimator | SE (pp) | Bootstrap samples |
|---|---|---|
| \(\widehat{SE}_{EIF}\) | 0.762 | NA |
| \(\widehat{SE}_{boot}\) | 0.788 | 200 |
Interpretation. The bootstrap distribution is centered near the original DR estimate, shown by the red vertical line. The dashed black lines are the normal 95% interval based on the estimated EIF standard error, not on the bootstrap quantiles. The fact that the bootstrap standard error, \(0.788\) percentage points, is close to the EIF standard error, \(0.762\) percentage points, is reassuring: two different uncertainty calculations are giving nearly the same scale.
This comparison supports the asymptotic linearity approximation
\[ \hat\psi_{DR}-\psi \approx P_n D^*(O). \]
In words, once the nuisance functions are fit, the remaining sampling variability of the DR estimator behaves much like the empirical average of its influence-function contributions. The bootstrap histogram is slightly irregular because it uses only 200 resamples and refits the nuisance models each time, but its width is broadly consistent with the analytic EIF-based interval. That is exactly what we hope to see before trusting Wald-style confidence intervals from semiparametric theory.
18 Simulation 1: Efficiency Bound in a Known DGP
Now we simulate from a known observational data-generating process. This lets us compute the true \(\psi\) and the true efficiency bound.
The data-generating process is:
\[ X_1,X_2 \sim N(0,1), \qquad X_3 \sim \mathrm{Bernoulli}(0.45), \]
\[ e(X) = \mathrm{expit}(-0.25 + 0.70X_1 - 0.55X_2 + 0.35X_3), \]
\[ Q_0(X) = \mathrm{expit}(-2.10 + 0.80X_1 - 0.45X_2 + 0.35X_3), \]
and
\[ Q_1(X) = \mathrm{expit}\{\mathrm{logit}(Q_0(X)) + 0.45 + 0.35X_1\}. \]
Then \(T \sim \mathrm{Bernoulli}\{e(X)\}\), and \(Y \sim \mathrm{Bernoulli}\{Q_T(X)\}\).
generate_dgp <- function(n, overlap_strength = 1) {
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rbinom(n, size = 1, prob = 0.45)
e <- inv_logit(overlap_strength * (-0.25 + 0.70 * x1 - 0.55 * x2 + 0.35 * x3))
q0 <- inv_logit(-2.10 + 0.80 * x1 - 0.45 * x2 + 0.35 * x3)
q1 <- inv_logit(safe_logit(q0) + 0.45 + 0.35 * x1)
t <- rbinom(n, size = 1, prob = e)
y <- rbinom(n, size = 1, prob = ifelse(t == 1, q1, q0))
data.frame(
x1 = x1,
x2 = x2,
x3 = x3,
t = t,
y = y,
e_true = e,
q1_true = q1,
q0_true = q0
)
}
truth_from_large_population <- function(n_pop = 200000, overlap_strength = 1) {
pop <- generate_dgp(n_pop, overlap_strength = overlap_strength)
psi <- mean(pop$q1_true - pop$q0_true)
v_heterogeneity <- var(pop$q1_true - pop$q0_true)
v_residual <- mean(
pop$q1_true * (1 - pop$q1_true) / pop$e_true +
pop$q0_true * (1 - pop$q0_true) / (1 - pop$e_true)
)
c(
psi = psi,
v_eff = v_heterogeneity + v_residual,
heterogeneity = v_heterogeneity,
residual = v_residual
)
}
estimate_sim_methods <- function(dat,
e_spec = c("correct", "wrong"),
q_spec = c("correct", "wrong")) {
e_spec <- match.arg(e_spec)
q_spec <- match.arg(q_spec)
if (e_spec == "correct") {
e_formula <- t ~ x1 + x2 + x3
} else {
e_formula <- t ~ x1
}
if (q_spec == "correct") {
q_formula <- y ~ t + x1 + x2 + x3 + t:x1
} else {
q_formula <- y ~ t + x1
}
e_fit <- glm(e_formula, data = dat, family = binomial())
q_fit <- glm(q_formula, data = dat, family = binomial())
e_hat_s <- bound_probability(predict(e_fit, type = "response"))
d1 <- dat
d0 <- dat
d1$t <- 1
d0$t <- 0
q1_hat_s <- bound_probability(predict(q_fit, newdata = d1, type = "response"))
q0_hat_s <- bound_probability(predict(q_fit, newdata = d0, type = "response"))
crude <- mean(dat$y[dat$t == 1]) - mean(dat$y[dat$t == 0])
ipw <- weighted_mean(dat$y[dat$t == 1], 1 / e_hat_s[dat$t == 1]) -
weighted_mean(dat$y[dat$t == 0], 1 / (1 - e_hat_s[dat$t == 0]))
outcome_regression <- mean(q1_hat_s - q0_hat_s)
dr <- mean(
q1_hat_s - q0_hat_s +
dat$t / e_hat_s * (dat$y - q1_hat_s) -
(1 - dat$t) / (1 - e_hat_s) * (dat$y - q0_hat_s)
)
oracle <- mean(
dat$q1_true - dat$q0_true +
dat$t / dat$e_true * (dat$y - dat$q1_true) -
(1 - dat$t) / (1 - dat$e_true) * (dat$y - dat$q0_true)
)
c(
Crude = crude,
IPW = ipw,
Outcome_regression = outcome_regression,
Doubly_robust = dr,
Oracle_EIF = oracle
)
}
parametric_or_se <- function(q_fit, d1, d0, B = 120) {
beta_hat <- coef(q_fit)
vc <- vcov(q_fit)
if (any(!is.finite(beta_hat)) || any(!is.finite(vc))) {
return(NA_real_)
}
p <- length(beta_hat)
vc <- vc + diag(1e-10, p)
coef_draws <- tryCatch(
{
z <- matrix(rnorm(B * p), nrow = B, ncol = p)
sweep(z %*% chol(vc), 2, beta_hat, "+")
},
error = function(e) NULL
)
if (is.null(coef_draws)) {
return(NA_real_)
}
x1_mat <- model.matrix(delete.response(terms(q_fit)), data = d1)
x0_mat <- model.matrix(delete.response(terms(q_fit)), data = d0)
q1_draws <- inv_logit(x1_mat %*% t(coef_draws))
q0_draws <- inv_logit(x0_mat %*% t(coef_draws))
psi_draws <- colMeans(q1_draws - q0_draws)
sd(psi_draws)
}
estimate_sim_methods_with_ci <- function(dat,
e_spec = c("correct", "wrong"),
q_spec = c("correct", "wrong")) {
e_spec <- match.arg(e_spec)
q_spec <- match.arg(q_spec)
if (e_spec == "correct") {
e_formula <- t ~ x1 + x2 + x3
} else {
e_formula <- t ~ x1
}
if (q_spec == "correct") {
q_formula <- y ~ t + x1 + x2 + x3 + t:x1
} else {
q_formula <- y ~ t + x1
}
e_fit <- glm(e_formula, data = dat, family = binomial())
q_fit <- glm(q_formula, data = dat, family = binomial())
e_hat_s <- bound_probability(predict(e_fit, type = "response"))
d1 <- dat
d0 <- dat
d1$t <- 1
d0$t <- 0
q1_hat_s <- bound_probability(predict(q_fit, newdata = d1, type = "response"))
q0_hat_s <- bound_probability(predict(q_fit, newdata = d0, type = "response"))
n <- nrow(dat)
treated <- dat$t == 1
untreated <- dat$t == 0
crude <- mean(dat$y[treated]) - mean(dat$y[untreated])
crude_se <- sqrt(
var(dat$y[treated]) / sum(treated) +
var(dat$y[untreated]) / sum(untreated)
)
ipw_psi1 <- weighted_mean(dat$y[treated], 1 / e_hat_s[treated])
ipw_psi0 <- weighted_mean(dat$y[untreated], 1 / (1 - e_hat_s[untreated]))
ipw <- ipw_psi1 - ipw_psi0
ipw_den1 <- mean(dat$t / e_hat_s)
ipw_den0 <- mean((1 - dat$t) / (1 - e_hat_s))
ipw_if <- dat$t / e_hat_s * (dat$y - ipw_psi1) / ipw_den1 -
(1 - dat$t) / (1 - e_hat_s) * (dat$y - ipw_psi0) / ipw_den0
ipw_se <- sd(ipw_if) / sqrt(n)
outcome_regression <- mean(q1_hat_s - q0_hat_s)
or_se <- parametric_or_se(q_fit, d1, d0)
if (!is.finite(or_se)) {
or_se <- sd(q1_hat_s - q0_hat_s) / sqrt(n)
}
dr_component <- q1_hat_s - q0_hat_s +
dat$t / e_hat_s * (dat$y - q1_hat_s) -
(1 - dat$t) / (1 - e_hat_s) * (dat$y - q0_hat_s)
dr <- mean(dr_component)
dr_se <- sd(dr_component - dr) / sqrt(n)
oracle_component <- dat$q1_true - dat$q0_true +
dat$t / dat$e_true * (dat$y - dat$q1_true) -
(1 - dat$t) / (1 - dat$e_true) * (dat$y - dat$q0_true)
oracle <- mean(oracle_component)
oracle_se <- sd(oracle_component - oracle) / sqrt(n)
estimates <- c(
Crude = crude,
IPW = ipw,
Outcome_regression = outcome_regression,
Doubly_robust = dr,
Oracle_EIF = oracle
)
ses <- c(
Crude = crude_se,
IPW = ipw_se,
Outcome_regression = or_se,
Doubly_robust = dr_se,
Oracle_EIF = oracle_se
)
data.frame(
Method = names(estimates),
Estimate = unname(estimates),
SE = unname(ses),
Lower = unname(estimates - 1.96 * ses),
Upper = unname(estimates + 1.96 * ses),
row.names = NULL
)
}
estimate_dr_with_eif <- function(dat,
e_spec = c("correct", "wrong"),
q_spec = c("correct", "wrong")) {
e_spec <- match.arg(e_spec)
q_spec <- match.arg(q_spec)
if (e_spec == "correct") {
e_formula <- t ~ x1 + x2 + x3
} else {
e_formula <- t ~ x1
}
if (q_spec == "correct") {
q_formula <- y ~ t + x1 + x2 + x3 + t:x1
} else {
q_formula <- y ~ t + x1
}
e_fit <- glm(e_formula, data = dat, family = binomial())
q_fit <- glm(q_formula, data = dat, family = binomial())
e_hat_s <- bound_probability(predict(e_fit, type = "response"))
d1 <- dat
d0 <- dat
d1$t <- 1
d0$t <- 0
q1_hat_s <- bound_probability(predict(q_fit, newdata = d1, type = "response"))
q0_hat_s <- bound_probability(predict(q_fit, newdata = d0, type = "response"))
dr_component <- q1_hat_s - q0_hat_s +
dat$t / e_hat_s * (dat$y - q1_hat_s) -
(1 - dat$t) / (1 - e_hat_s) * (dat$y - q0_hat_s)
psi_hat <- mean(dr_component)
eif_hat <- dr_component - psi_hat
c(
psi = psi_hat,
se = sd(eif_hat) / sqrt(nrow(dat)),
lower = psi_hat - 1.96 * sd(eif_hat) / sqrt(nrow(dat)),
upper = psi_hat + 1.96 * sd(eif_hat) / sqrt(nrow(dat)),
e_min = min(e_hat_s),
e_max = max(e_hat_s),
max_inverse_weight = max(1 / e_hat_s, 1 / (1 - e_hat_s))
)
}
oracle_eif_with_se <- function(dat) {
oracle_component <- dat$q1_true - dat$q0_true +
dat$t / dat$e_true * (dat$y - dat$q1_true) -
(1 - dat$t) / (1 - dat$e_true) * (dat$y - dat$q0_true)
psi_hat <- mean(oracle_component)
eif_hat <- oracle_component - psi_hat
c(
psi = psi_hat,
se = sd(eif_hat) / sqrt(nrow(dat)),
lower = psi_hat - 1.96 * sd(eif_hat) / sqrt(nrow(dat)),
upper = psi_hat + 1.96 * sd(eif_hat) / sqrt(nrow(dat))
)
}
summarize_ci_performance <- function(results, truth, v_eff = NULL) {
split_results <- split(results, list(results$Method, results$n), drop = TRUE)
out <- lapply(split_results, function(current) {
n_current <- unique(current$n)
method_current <- unique(current$Method)
estimate <- current$Estimate
se <- current$SE
data.frame(
n = n_current,
Method = method_current,
Bias = mean(estimate) - truth,
Monte_Carlo_SD = sd(estimate),
Mean_EIF_SE = mean(se),
Coverage = mean(current$Lower <= truth & current$Upper >= truth),
Efficiency_ratio = if (is.null(v_eff)) {
NA_real_
} else {
var(estimate) / (v_eff / n_current)
},
row.names = NULL
)
})
do.call(rbind, out)
}
evaluate_aipw_with_nuisances <- function(dat, e_hat_s, q1_hat_s, q0_hat_s) {
mean(
q1_hat_s - q0_hat_s +
dat$t / e_hat_s * (dat$y - q1_hat_s) -
(1 - dat$t) / (1 - e_hat_s) * (dat$y - q0_hat_s)
)
}
make_rich_features <- function(dat) {
dat$x1_2 <- dat$x1^2
dat$x2_2 <- dat$x2^2
dat$x1_3 <- dat$x1^3
dat$x2_3 <- dat$x2^3
dat$x1x2 <- dat$x1 * dat$x2
dat
}
rich_e_formula <- t ~ x1 + x2 + x3 + x1_2 + x2_2 + x1_3 + x2_3 + x1x2
rich_q_formula <- y ~ t * (x1 + x2 + x3 + x1_2 + x2_2 + x1_3 + x2_3 + x1x2)
estimate_dr_with_formula <- function(dat, e_formula, q_formula) {
dat_rich <- make_rich_features(dat)
e_fit <- glm(e_formula, data = dat_rich, family = binomial())
q_fit <- glm(q_formula, data = dat_rich, family = binomial())
e_hat_s <- bound_probability(predict(e_fit, type = "response"))
d1 <- dat_rich
d0 <- dat_rich
d1$t <- 1
d0$t <- 0
q1_hat_s <- bound_probability(predict(q_fit, newdata = d1, type = "response"))
q0_hat_s <- bound_probability(predict(q_fit, newdata = d0, type = "response"))
dr_component <- q1_hat_s - q0_hat_s +
dat_rich$t / e_hat_s * (dat_rich$y - q1_hat_s) -
(1 - dat_rich$t) / (1 - e_hat_s) * (dat_rich$y - q0_hat_s)
psi_hat <- mean(dr_component)
eif_hat <- dr_component - psi_hat
c(
psi = psi_hat,
se = sd(eif_hat) / sqrt(nrow(dat_rich)),
lower = psi_hat - 1.96 * sd(eif_hat) / sqrt(nrow(dat_rich)),
upper = psi_hat + 1.96 * sd(eif_hat) / sqrt(nrow(dat_rich))
)
}
estimate_dr_crossfit <- function(dat, e_formula, q_formula, K = 5) {
dat_rich <- make_rich_features(dat)
n <- nrow(dat_rich)
folds <- sample(rep(seq_len(K), length.out = n))
e_hat_s <- q1_hat_s <- q0_hat_s <- rep(NA_real_, n)
for (k in seq_len(K)) {
valid <- which(folds == k)
train <- which(folds != k)
e_fit <- glm(e_formula, data = dat_rich[train, ], family = binomial())
q_fit <- glm(q_formula, data = dat_rich[train, ], family = binomial())
valid_data <- dat_rich[valid, ]
d1 <- valid_data
d0 <- valid_data
d1$t <- 1
d0$t <- 0
e_hat_s[valid] <- bound_probability(
predict(e_fit, newdata = valid_data, type = "response")
)
q1_hat_s[valid] <- bound_probability(
predict(q_fit, newdata = d1, type = "response")
)
q0_hat_s[valid] <- bound_probability(
predict(q_fit, newdata = d0, type = "response")
)
}
dr_component <- q1_hat_s - q0_hat_s +
dat_rich$t / e_hat_s * (dat_rich$y - q1_hat_s) -
(1 - dat_rich$t) / (1 - e_hat_s) * (dat_rich$y - q0_hat_s)
psi_hat <- mean(dr_component)
eif_hat <- dr_component - psi_hat
c(
psi = psi_hat,
se = sd(eif_hat) / sqrt(n),
lower = psi_hat - 1.96 * sd(eif_hat) / sqrt(n),
upper = psi_hat + 1.96 * sd(eif_hat) / sqrt(n)
)
}set.seed(20260706)
sim_truth <- truth_from_large_population()
sim_truth_table <- data.frame(
Quantity = truth_math_label(names(sim_truth)),
Value = round(unname(sim_truth), 5)
)
mkable(sim_truth_table)| Quantity | Value |
|---|---|
| \(\psi\) | 0.08069 |
| \(V_{\mathrm{EIF}} = E\{D^*(O)^2\}\) | 0.65295 |
| \(\mathrm{Var}\{Q_1(X)-Q_0(X)\}\) | 0.00607 |
| \(E[Q_1(1-Q_1)/e + Q_0(1-Q_0)/(1-e)]\) | 0.64688 |
The oracle EIF estimator is not available in practice because it uses the true nuisance functions. It is useful in simulation because its variance should match the efficiency bound.
set.seed(20260706)
R_eff <- 200
n_sim <- 800
sim_method_names <- c(
"Crude",
"IPW",
"Outcome_regression",
"Doubly_robust",
"Oracle_EIF"
)
sim_estimates <- matrix(NA_real_, nrow = R_eff, ncol = length(sim_method_names))
sim_ses <- sim_lower <- sim_upper <- sim_estimates
colnames(sim_estimates) <- colnames(sim_ses) <-
colnames(sim_lower) <- colnames(sim_upper) <- sim_method_names
for (r in seq_len(R_eff)) {
dat_r <- generate_dgp(n_sim)
sim_fit_r <- tryCatch(
estimate_sim_methods_with_ci(dat_r, e_spec = "correct", q_spec = "correct"),
error = function(e) data.frame(
Method = sim_method_names,
Estimate = NA_real_,
SE = NA_real_,
Lower = NA_real_,
Upper = NA_real_
)
)
sim_fit_r <- sim_fit_r[match(sim_method_names, sim_fit_r$Method), ]
sim_estimates[r, ] <- sim_fit_r$Estimate
sim_ses[r, ] <- sim_fit_r$SE
sim_lower[r, ] <- sim_fit_r$Lower
sim_upper[r, ] <- sim_fit_r$Upper
}
complete_sim_rows <- complete.cases(sim_estimates) &
complete.cases(sim_ses) &
complete.cases(sim_lower) &
complete.cases(sim_upper)
sim_estimates <- sim_estimates[complete_sim_rows, , drop = FALSE]
sim_ses <- sim_ses[complete_sim_rows, , drop = FALSE]
sim_lower <- sim_lower[complete_sim_rows, , drop = FALSE]
sim_upper <- sim_upper[complete_sim_rows, , drop = FALSE]
summarize_mc <- function(estimates, ses, lower, upper, truth) {
data.frame(
Method = method_math_label(colnames(estimates)),
Mean = apply(estimates, 2, mean),
Bias = apply(estimates, 2, mean) - truth,
Monte_Carlo_SD = apply(estimates, 2, sd),
Mean_SE = apply(ses, 2, mean),
Coverage_95 = colMeans(lower <= truth & upper >= truth),
RMSE = sqrt(apply((t(t(estimates) - truth))^2, 2, mean)),
row.names = NULL
)
}
mc_eff_table <- summarize_mc(
sim_estimates,
sim_ses,
sim_lower,
sim_upper,
sim_truth["psi"]
)
mc_eff_table[, -1] <- round(mc_eff_table[, -1], 4)
names(mc_eff_table) <- c(
"Estimator",
"\\(E_{MC}(\\hat\\psi)\\)",
"Bias",
"MC SD",
"Mean SE",
"95% CI coverage",
"RMSE"
)
theory_row <- data.frame(
Estimator = "\\(\\sqrt{V_{\\mathrm{EIF}}/n}\\)",
Mean = NA_real_,
Bias = NA_real_,
Monte_Carlo_SD = round(sqrt(sim_truth["v_eff"] / n_sim), 4),
Mean_SE = round(sqrt(sim_truth["v_eff"] / n_sim), 4),
Coverage_95 = NA_real_,
RMSE = NA_real_
)
names(theory_row) <- names(mc_eff_table)
mkable(rbind(mc_eff_table, theory_row))| Estimator | \(E_{MC}(\hat\psi)\) | Bias | MC SD | Mean SE | 95% CI coverage | RMSE |
|---|---|---|---|---|---|---|
| \(\hat\psi_{crude}\) | 0.1913 | 0.1106 | 0.0289 | 0.0283 | 0.025 | 0.1143 |
| \(\hat\psi_{IPW}\) | 0.0810 | 0.0003 | 0.0286 | 0.0317 | 0.965 | 0.0285 |
| \(\hat\psi_{OR}\) | 0.0806 | -0.0001 | 0.0286 | 0.0277 | 0.940 | 0.0285 |
| \(\hat\psi_{DR}\) | 0.0807 | 0.0000 | 0.0288 | 0.0286 | 0.950 | 0.0287 |
| \(\hat\psi_{oracle\ EIF}\) | 0.0809 | 0.0002 | 0.0287 | 0.0285 | 0.945 | 0.0286 |
| \(\sqrt{V_{\mathrm{EIF}}/n}\) | NA | NA | 0.0286 | 0.0286 | NA | NA |
Two comments are important.
First, the Oracle EIF row is the oracle benchmark for the nonparametric observed-data model: it uses the true \(e(X)\), \(Q_1(X)\), and \(Q_0(X)\), and its Monte Carlo SD should be close to \(\sqrt{V_{\mathrm{EIF}}/n}\).
Second, the outcome regression row may have a smaller RMSE in this particular simulation because the simulated outcome model is exactly logistic-linear and the estimator is allowed to exploit that correct parametric restriction. That does not violate the semiparametric lower bound. It means the outcome-regression estimator is being judged in a smaller model than the nonparametric model used to define \(V_{\mathrm{EIF}}\). If that logistic outcome model were wrong, the same estimator could have persistent bias.
When the nuisance models are correctly specified, the DR estimator approaches the oracle EIF benchmark. In finite samples, it may not be exactly equal to the bound because nuisance functions are estimated, models are finite-dimensional, and the Monte Carlo size is limited. The parametric outcome-regression estimator is included as a warning: if the simulation makes its working model exactly true, it can look more efficient than the nonparametric semiparametric benchmark because it is borrowing extra model information.
18.1 Simulation 1A: Sample Size, Coverage, and the Efficiency Ratio
This simulation has one purpose: to check the theorem in the cleanest possible setting.
The theorem does not say that every reasonable estimator must have variance \(V_{\mathrm{EIF}}/n\). It says that the estimator whose first-order expansion is governed by the efficient influence function has asymptotic variance \(V_{\mathrm{EIF}}/n\), and that no regular estimator in the same nonparametric observed-data model can do better.
So here we remove a major source of distraction: we use the true nuisance functions \(e(X)\), \(Q_1(X)\), and \(Q_0(X)\). This gives the oracle one-step estimator
\[ \hat\psi_{\mathrm{EIF}} = P_n\left[ Q_1(X)-Q_0(X) + \frac{T}{e(X)}\{Y-Q_1(X)\} - \frac{1-T}{1-e(X)}\{Y-Q_0(X)\} \right]. \]
If the efficiency-bound calculation is correct, then as \(n\) increases we should see:
- Monte Carlo bias close to 0.
- Monte Carlo standard deviation close to \(\sqrt{V_{\mathrm{EIF}}/n}\).
- Mean EIF-based standard error close to both quantities.
- Wald coverage close to 0.95.
- Efficiency ratio \(\mathrm{Var}_{MC}(\hat\psi_{\mathrm{EIF}})/(V_{\mathrm{EIF}}/n)\) close to 1.
Suppose we repeat the simulation \(R\) times at the same sample size \(n\). Let \(\hat\psi^{(r)}\) denote the estimate from repetition \(r\), and let
\[ \bar{\hat\psi}_{MC} = \frac{1}{R}\sum_{r=1}^R \hat\psi^{(r)}. \]
The Monte Carlo standard deviation is the empirical standard deviation of the repeated estimates:
\[ \widehat{\mathrm{SD}}_{MC}(\hat\psi) = \left[ \frac{1}{R-1} \sum_{r=1}^R \left\{\hat\psi^{(r)}-\bar{\hat\psi}_{MC}\right\}^2 \right]^{1/2}. \]
In this simulation, the theorem predicts
\[ \widehat{\mathrm{SD}}_{MC}(\hat\psi_{\mathrm{EIF}}) \approx \sqrt{\frac{V_{\mathrm{EIF}}}{n}}, \]
up to Monte Carlo error from using a finite number of repetitions.
set.seed(20260706)
sample_size_grid <- c(250, 500, 1000, 2000, 5000)
R_size <- rep(600, length(sample_size_grid))
oracle_theorem_results <- vector("list", length(sample_size_grid))
for (j in seq_along(sample_size_grid)) {
n_j <- sample_size_grid[j]
R_j <- R_size[j]
current_results <- lapply(seq_len(R_j), function(r) {
dat_r <- generate_dgp(n_j)
fit_r <- oracle_eif_with_se(dat_r)
data.frame(
n = n_j,
Method = "Oracle_EIF",
Estimate = fit_r["psi"],
SE = fit_r["se"],
Lower = fit_r["lower"],
Upper = fit_r["upper"],
row.names = NULL
)
})
oracle_theorem_results[[j]] <- do.call(rbind, current_results)
}
oracle_theorem_results <- do.call(rbind, oracle_theorem_results)
oracle_theorem_summary <- summarize_ci_performance(
oracle_theorem_results,
truth = sim_truth["psi"],
v_eff = sim_truth["v_eff"]
)
oracle_theorem_summary <- oracle_theorem_summary[order(oracle_theorem_summary$n), ]
oracle_theorem_summary$MC_Reps <- vapply(
seq_len(nrow(oracle_theorem_summary)),
function(i) {
sum(oracle_theorem_results$n == oracle_theorem_summary$n[i])
},
integer(1)
)
oracle_theorem_summary$Theoretical_SE <- sqrt(sim_truth["v_eff"] / oracle_theorem_summary$n)
oracle_theorem_summary$Ratio_MC_SE <- oracle_theorem_summary$Efficiency_ratio *
sqrt(2 / pmax(oracle_theorem_summary$MC_Reps - 1, 1))
oracle_theorem_display <- oracle_theorem_summary[, c(
"n", "MC_Reps", "Bias", "Monte_Carlo_SD", "Mean_EIF_SE",
"Theoretical_SE", "Coverage", "Efficiency_ratio"
)]
oracle_theorem_display[, c(
"Bias", "Monte_Carlo_SD", "Mean_EIF_SE",
"Theoretical_SE", "Coverage", "Efficiency_ratio"
)] <- round(
oracle_theorem_display[, c(
"Bias", "Monte_Carlo_SD", "Mean_EIF_SE",
"Theoretical_SE", "Coverage", "Efficiency_ratio"
)],
3
)
names(oracle_theorem_display) <- c(
"\\(n\\)",
"MC reps",
"Bias",
"MC SD",
"Mean \\(\\widehat{SE}_{EIF}\\)",
"\\(\\sqrt{V_{\\mathrm{EIF}}/n}\\)",
"95% CI coverage",
"\\(\\mathrm{Var}_{MC}(\\hat\\psi)/(V_{\\mathrm{EIF}}/n)\\)"
)
mkable(oracle_theorem_display)| \(n\) | MC reps | Bias | MC SD | Mean \(\widehat{SE}_{EIF}\) | \(\sqrt{V_{\mathrm{EIF}}/n}\) | 95% CI coverage | \(\mathrm{Var}_{MC}(\hat\psi)/(V_{\mathrm{EIF}}/n)\) |
|---|---|---|---|---|---|---|---|
| 250 | 600 | -0.004 | 0.048 | 0.051 | 0.051 | 0.957 | 0.896 |
| 500 | 600 | 0.001 | 0.037 | 0.036 | 0.036 | 0.945 | 1.021 |
| 1000 | 600 | 0.000 | 0.026 | 0.026 | 0.026 | 0.957 | 1.034 |
| 2000 | 600 | 0.000 | 0.018 | 0.018 | 0.018 | 0.953 | 1.034 |
| 5000 | 600 | 0.000 | 0.012 | 0.011 | 0.011 | 0.952 | 1.015 |
This is the desired pattern. The green Monte Carlo standard deviation, the blue EIF standard error, and the red theoretical curve all move together on the \(n^{-1/2}\) scale. The efficiency ratio stays near 1, up to ordinary Monte Carlo error from estimating a variance with \(R=600\) repeated simulations.
This is why we do the experiment: it connects the abstract lower-bound formula \(V_{\mathrm{EIF}}=E\{D^*(O)^2\}\) to observable repeated-sampling behavior. Once that connection is clear, feasible DR estimation is easier to understand. A feasible DR estimator tries to behave like this oracle EIF estimator after replacing \(e,Q_1,Q_0\) by estimates. If the nuisance estimators are good enough, its first-order behavior is governed by the same EIF; if they are unstable or misspecified, the finite-sample behavior can deviate.
IPW and OR are useful comparisons, but they answer different questions. IPW does not use the outcome regression information in \(Q_t(X)\), so it often pays a variance price. OR can look extremely good when the outcome regression model is exactly correct, but that is because the simulation has given it extra parametric structure. The clean efficiency-bound check is the oracle EIF pattern above.
18.2 Simulation 1B: When the EIF/DR Estimator Beats Single-Robust Estimators
The first simulation intentionally made the basic parametric outcome model correct. That is useful for understanding the efficiency bound, but it can hide the practical advantage of EIF-based estimation. Here we add a more demanding data-generating process with nonlinear confounding and treatment-effect heterogeneity:
\[ \mathrm{logit}\{e(X)\} = -0.20+0.70X_1-0.65X_2+0.35X_3 +0.55X_1X_2-0.30X_1^2+0.25\sin(X_2), \]
and
\[ \mathrm{logit}\{Q_0(X)\} = -2.05+0.70X_1-0.45X_2+0.35X_3 +0.60X_1X_2-0.35X_2^2+0.35\sin(2X_1), \]
with
\[ \mathrm{logit}\{Q_1(X)\} = \mathrm{logit}\{Q_0(X)\} +0.45+0.40X_1-0.25X_2+0.45X_1X_2. \]
We compare:
- Simple IPW using only \(X_1,X_2,X_3\) in the propensity model.
- Simple OR using only \(T,X_1,X_2,X_3\) in the outcome model.
- A richer EIF/DR estimator that includes nonlinear and interaction terms in both nuisance fits, with estimated propensities truncated to \([0.025,0.975]\) to avoid finite-sample inverse-weight explosions.
- The infeasible oracle EIF estimator.
This is a practical stress test, not a theorem that DR always beats OR or IPW. The lesson is sharper: when single-robust working models are under-specified, the EIF/DR construction can be much closer to the oracle benchmark if its nuisance fits capture the important structure.
make_advantage_features <- function(dat) {
dat$x1_2 <- dat$x1^2
dat$x2_2 <- dat$x2^2
dat$x1x2 <- dat$x1 * dat$x2
dat$sin_x2 <- sin(dat$x2)
dat$sin2_x1 <- sin(2 * dat$x1)
dat
}
generate_advantage_dgp <- function(n) {
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rbinom(n, size = 1, prob = 0.45)
x1_2 <- x1^2
x2_2 <- x2^2
x1x2 <- x1 * x2
sin_x2 <- sin(x2)
sin2_x1 <- sin(2 * x1)
e <- bound_probability(
inv_logit(
-0.20 + 0.70 * x1 - 0.65 * x2 + 0.35 * x3 +
0.55 * x1x2 - 0.30 * x1_2 + 0.25 * sin_x2
),
eps = 0.01
)
q0 <- bound_probability(
inv_logit(
-2.05 + 0.70 * x1 - 0.45 * x2 + 0.35 * x3 +
0.60 * x1x2 - 0.35 * x2_2 + 0.35 * sin2_x1
),
eps = 0.01
)
q1 <- bound_probability(
inv_logit(safe_logit(q0) + 0.45 + 0.40 * x1 - 0.25 * x2 + 0.45 * x1x2),
eps = 0.01
)
t <- rbinom(n, size = 1, prob = e)
y <- rbinom(n, size = 1, prob = ifelse(t == 1, q1, q0))
data.frame(
x1 = x1,
x2 = x2,
x3 = x3,
t = t,
y = y,
e_true = e,
q1_true = q1,
q0_true = q0
)
}
truth_advantage_dgp <- function(n_pop = 250000) {
pop <- generate_advantage_dgp(n_pop)
psi <- mean(pop$q1_true - pop$q0_true)
v_eff <- var(pop$q1_true - pop$q0_true) +
mean(
pop$q1_true * (1 - pop$q1_true) / pop$e_true +
pop$q0_true * (1 - pop$q0_true) / (1 - pop$e_true)
)
c(psi = psi, v_eff = v_eff)
}
estimate_advantage_methods_with_ci <- function(dat) {
dat_rich <- make_advantage_features(dat)
n <- nrow(dat_rich)
e_simple <- glm(t ~ x1 + x2 + x3, data = dat_rich, family = binomial())
q_simple <- glm(y ~ t + x1 + x2 + x3, data = dat_rich, family = binomial())
e_rich <- glm(
t ~ x1 + x2 + x3 + x1_2 + x2_2 + x1x2 + sin_x2,
data = dat_rich,
family = binomial()
)
q_rich <- glm(
y ~ t * (x1 + x2 + x3 + x1_2 + x2_2 + x1x2 + sin2_x1),
data = dat_rich,
family = binomial()
)
e_simple_hat <- bound_probability(predict(e_simple, type = "response"), eps = 0.025)
e_rich_hat <- bound_probability(predict(e_rich, type = "response"), eps = 0.025)
d1 <- dat_rich
d0 <- dat_rich
d1$t <- 1
d0$t <- 0
q1_simple <- bound_probability(predict(q_simple, newdata = d1, type = "response"))
q0_simple <- bound_probability(predict(q_simple, newdata = d0, type = "response"))
q1_rich <- bound_probability(predict(q_rich, newdata = d1, type = "response"))
q0_rich <- bound_probability(predict(q_rich, newdata = d0, type = "response"))
treated <- dat_rich$t == 1
untreated <- dat_rich$t == 0
crude <- mean(dat_rich$y[treated]) - mean(dat_rich$y[untreated])
crude_se <- sqrt(
var(dat_rich$y[treated]) / sum(treated) +
var(dat_rich$y[untreated]) / sum(untreated)
)
ipw_psi1 <- weighted_mean(dat_rich$y[treated], 1 / e_simple_hat[treated])
ipw_psi0 <- weighted_mean(dat_rich$y[untreated], 1 / (1 - e_simple_hat[untreated]))
ipw <- ipw_psi1 - ipw_psi0
ipw_den1 <- mean(dat_rich$t / e_simple_hat)
ipw_den0 <- mean((1 - dat_rich$t) / (1 - e_simple_hat))
ipw_if <- dat_rich$t / e_simple_hat * (dat_rich$y - ipw_psi1) / ipw_den1 -
(1 - dat_rich$t) / (1 - e_simple_hat) * (dat_rich$y - ipw_psi0) / ipw_den0
ipw_se <- sd(ipw_if) / sqrt(n)
outcome_regression <- mean(q1_simple - q0_simple)
or_se <- parametric_or_se(q_simple, d1, d0, B = 80)
if (!is.finite(or_se)) {
or_se <- sd(q1_simple - q0_simple) / sqrt(n)
}
dr_component <- q1_rich - q0_rich +
dat_rich$t / e_rich_hat * (dat_rich$y - q1_rich) -
(1 - dat_rich$t) / (1 - e_rich_hat) * (dat_rich$y - q0_rich)
dr <- mean(dr_component)
dr_se <- sd(dr_component - dr) / sqrt(n)
oracle_component <- dat_rich$q1_true - dat_rich$q0_true +
dat_rich$t / dat_rich$e_true * (dat_rich$y - dat_rich$q1_true) -
(1 - dat_rich$t) / (1 - dat_rich$e_true) * (dat_rich$y - dat_rich$q0_true)
oracle <- mean(oracle_component)
oracle_se <- sd(oracle_component - oracle) / sqrt(n)
estimates <- c(
Crude = crude,
IPW = ipw,
Outcome_regression = outcome_regression,
Doubly_robust = dr,
Oracle_EIF = oracle
)
ses <- c(
Crude = crude_se,
IPW = ipw_se,
Outcome_regression = or_se,
Doubly_robust = dr_se,
Oracle_EIF = oracle_se
)
data.frame(
Method = names(estimates),
Estimate = unname(estimates),
SE = unname(ses),
Lower = unname(estimates - 1.96 * ses),
Upper = unname(estimates + 1.96 * ses),
row.names = NULL
)
}set.seed(20260706)
advantage_truth <- truth_advantage_dgp()
R_advantage <- 120
n_advantage <- 1200
advantage_method_names <- c(
"Crude",
"IPW",
"Outcome_regression",
"Doubly_robust",
"Oracle_EIF"
)
advantage_estimates <- matrix(
NA_real_,
nrow = R_advantage,
ncol = length(advantage_method_names),
dimnames = list(NULL, advantage_method_names)
)
advantage_ses <- advantage_lower <- advantage_upper <- advantage_estimates
for (r in seq_len(R_advantage)) {
fit_r <- tryCatch(
estimate_advantage_methods_with_ci(generate_advantage_dgp(n_advantage)),
error = function(e) data.frame(
Method = advantage_method_names,
Estimate = NA_real_,
SE = NA_real_,
Lower = NA_real_,
Upper = NA_real_
)
)
fit_r <- fit_r[match(advantage_method_names, fit_r$Method), ]
advantage_estimates[r, ] <- fit_r$Estimate
advantage_ses[r, ] <- fit_r$SE
advantage_lower[r, ] <- fit_r$Lower
advantage_upper[r, ] <- fit_r$Upper
}
complete_advantage_rows <- complete.cases(advantage_estimates) &
complete.cases(advantage_ses) &
complete.cases(advantage_lower) &
complete.cases(advantage_upper)
advantage_estimates <- advantage_estimates[complete_advantage_rows, , drop = FALSE]
advantage_ses <- advantage_ses[complete_advantage_rows, , drop = FALSE]
advantage_lower <- advantage_lower[complete_advantage_rows, , drop = FALSE]
advantage_upper <- advantage_upper[complete_advantage_rows, , drop = FALSE]
advantage_mc_table <- summarize_mc(
advantage_estimates,
advantage_ses,
advantage_lower,
advantage_upper,
advantage_truth["psi"]
)
advantage_label_map <- c(
Crude = "\\(\\hat\\psi_{crude}\\)",
IPW = "\\(\\hat\\psi_{IPW}^{simple}\\)",
Outcome_regression = "\\(\\hat\\psi_{OR}^{simple}\\)",
Doubly_robust = "\\(\\hat\\psi_{DR/EIF}^{rich}\\)",
Oracle_EIF = "\\(\\hat\\psi_{oracle\\ EIF}\\)"
)
advantage_mc_table$Method <- advantage_label_map[advantage_method_names]
advantage_mc_table[, -1] <- round(advantage_mc_table[, -1], 4)
names(advantage_mc_table) <- c(
"Estimator",
"\\(E_{MC}(\\hat\\psi)\\)",
"Bias",
"MC SD",
"Mean SE",
"95% CI coverage",
"RMSE"
)
advantage_bound_row <- data.frame(
Estimator = "\\(\\sqrt{V_{\\mathrm{EIF}}/n}\\)",
Mean = NA_real_,
Bias = NA_real_,
Monte_Carlo_SD = round(sqrt(advantage_truth["v_eff"] / n_advantage), 4),
Mean_SE = round(sqrt(advantage_truth["v_eff"] / n_advantage), 4),
Coverage_95 = NA_real_,
RMSE = NA_real_
)
names(advantage_bound_row) <- names(advantage_mc_table)
mkable(rbind(advantage_mc_table, advantage_bound_row))| Estimator | \(E_{MC}(\hat\psi)\) | Bias | MC SD | Mean SE | 95% CI coverage | RMSE |
|---|---|---|---|---|---|---|
| \(\hat\psi_{crude}\) | 0.1751 | 0.0804 | 0.0228 | 0.0232 | 0.0500 | 0.0835 |
| \(\hat\psi_{IPW}^{simple}\) | 0.1213 | 0.0266 | 0.0224 | 0.0243 | 0.8250 | 0.0347 |
| \(\hat\psi_{OR}^{simple}\) | 0.1249 | 0.0302 | 0.0214 | 0.0227 | 0.7750 | 0.0370 |
| \(\hat\psi_{DR/EIF}^{rich}\) | 0.0931 | -0.0016 | 0.0205 | 0.0215 | 0.9583 | 0.0205 |
| \(\hat\psi_{oracle\ EIF}\) | 0.0938 | -0.0009 | 0.0216 | 0.0220 | 0.9583 | 0.0215 |
| \(\sqrt{V_{\mathrm{EIF}}/n}\) | NA | NA | 0.0218 | 0.0218 | NA | NA |
In this experiment, the simple IPW and simple OR estimators are biased because the nuisance models omit nonlinear confounding and treatment-effect heterogeneity. The rich EIF/DR estimator is not simply averaging predictions; it uses the efficient influence-function correction
\[ \hat Q_1(X)-\hat Q_0(X) + \frac{T}{\hat e(X)}\{Y-\hat Q_1(X)\} - \frac{1-T}{1-\hat e(X)}\{Y-\hat Q_0(X)\}, \]
so the estimator behaves much closer to the oracle EIF benchmark. This is the simulation where the lecture message should feel practical: the EIF is not just a variance formula; it is also the recipe for a robust estimating equation.
19 Simulation 2: Double Robustness
Now we intentionally misspecify nuisance models.
| Scenario | Propensity model | Outcome model |
|---|---|---|
| both correct | correct | correct |
| only \(e\) correct | correct | wrong |
| only \(Q\) correct | wrong | correct |
| both wrong | wrong | wrong |
set.seed(20260706)
R_dr <- 120
n_dr <- 1000
scenario_table <- data.frame(
Scenario = c(
"Both correct",
"Only e correct",
"Only Q correct",
"Both wrong"
),
e_spec = c("correct", "correct", "wrong", "wrong"),
q_spec = c("correct", "wrong", "correct", "wrong"),
stringsAsFactors = FALSE
)
dr_results <- lapply(seq_len(nrow(scenario_table)), function(s) {
estimates <- matrix(NA_real_, nrow = R_dr, ncol = 5)
colnames(estimates) <- colnames(sim_estimates)
for (r in seq_len(R_dr)) {
dat_r <- generate_dgp(n_dr)
estimates[r, ] <- tryCatch(
estimate_sim_methods(
dat_r,
e_spec = scenario_table$e_spec[s],
q_spec = scenario_table$q_spec[s]
),
error = function(e) rep(NA_real_, 5)
)
}
estimates <- estimates[complete.cases(estimates), , drop = FALSE]
data.frame(
Scenario = scenario_table$Scenario[s],
Method = c("Outcome regression", "IPW", "Doubly robust"),
Bias = c(
mean(estimates[, "Outcome_regression"]) - sim_truth["psi"],
mean(estimates[, "IPW"]) - sim_truth["psi"],
mean(estimates[, "Doubly_robust"]) - sim_truth["psi"]
),
Monte_Carlo_SD = c(
sd(estimates[, "Outcome_regression"]),
sd(estimates[, "IPW"]),
sd(estimates[, "Doubly_robust"])
),
RMSE = c(
sqrt(mean((estimates[, "Outcome_regression"] - sim_truth["psi"])^2)),
sqrt(mean((estimates[, "IPW"] - sim_truth["psi"])^2)),
sqrt(mean((estimates[, "Doubly_robust"] - sim_truth["psi"])^2))
)
)
})
dr_results <- do.call(rbind, dr_results)
dr_display <- dr_results
dr_display[, c("Bias", "Monte_Carlo_SD", "RMSE")] <-
round(dr_display[, c("Bias", "Monte_Carlo_SD", "RMSE")], 4)
dr_display$Method <- method_math_label(gsub(" ", "_", dr_display$Method))
names(dr_display) <- c(
"Scenario",
"Estimator",
"Bias",
"MC SD",
"RMSE"
)
mkable(dr_display)| Scenario | Estimator | Bias | MC SD | RMSE |
|---|---|---|---|---|
| Both correct | \(\hat\psi_{OR}\) | -0.0003 | 0.0269 | 0.0267 |
| Both correct | \(\hat\psi_{IPW}\) | -0.0009 | 0.0289 | 0.0288 |
| Both correct | \(\hat\psi_{DR}\) | -0.0003 | 0.0277 | 0.0276 |
| Only e correct | \(\hat\psi_{OR}\) | 0.0324 | 0.0223 | 0.0393 |
| Only e correct | \(\hat\psi_{IPW}\) | 0.0014 | 0.0254 | 0.0254 |
| Only e correct | \(\hat\psi_{DR}\) | 0.0014 | 0.0247 | 0.0246 |
| Only Q correct | \(\hat\psi_{OR}\) | -0.0020 | 0.0245 | 0.0245 |
| Only Q correct | \(\hat\psi_{IPW}\) | 0.0331 | 0.0236 | 0.0406 |
| Only Q correct | \(\hat\psi_{DR}\) | -0.0017 | 0.0248 | 0.0247 |
| Both wrong | \(\hat\psi_{OR}\) | 0.0308 | 0.0267 | 0.0407 |
| Both wrong | \(\hat\psi_{IPW}\) | 0.0324 | 0.0264 | 0.0417 |
| Both wrong | \(\hat\psi_{DR}\) | 0.0325 | 0.0266 | 0.0419 |
The DR estimator is stable when either nuisance model is correct. When both are wrong, the product remainder no longer vanishes and bias can remain.
19.1 Simulation 2A: A Bias Surface for the Product Remainder
The previous experiment uses four discrete scenarios. The DR remainder is more geometric than that. If \(\hat e-e\) is on one axis and \(\hat Q-Q\) is on the other axis, the bias is approximately flat along the axes and grows when both errors are nonzero.
To isolate that idea, we generate one large dataset from the true DGP and then deliberately perturb the true nuisance functions:
\[ \hat e_\delta(X) = \mathrm{expit}\{\mathrm{logit}(e(X))+\delta_e g_e(X)\}, \]
and
\[ \hat Q_{t,\delta}(X) = \mathrm{expit}\{\mathrm{logit}(Q_t(X))+\delta_Q g_Q(X)\}. \]
If \(\delta_e=0\), the treatment model is correct. If \(\delta_Q=0\), the outcome model is correct. The DR estimating function should be nearly unbiased along either axis and biased away from the axes.
set.seed(20260706)
bias_dat <- generate_dgp(100000)
delta_e_grid <- seq(-1.2, 1.2, by = 0.3)
delta_q_grid <- seq(-1.2, 1.2, by = 0.3)
bias_surface <- expand.grid(
delta_e = delta_e_grid,
delta_q = delta_q_grid
)
g_e <- with(bias_dat, scale(x1 - 0.5 * x2 + 0.25 * x3)[, 1])
g_q <- with(bias_dat, scale(x2 + 0.35 * x1 - 0.40 * x3)[, 1])
bias_surface$Estimate <- NA_real_
for (i in seq_len(nrow(bias_surface))) {
e_tilt <- bound_probability(
inv_logit(safe_logit(bias_dat$e_true) + bias_surface$delta_e[i] * g_e)
)
q1_tilt <- bound_probability(
inv_logit(safe_logit(bias_dat$q1_true) + bias_surface$delta_q[i] * g_q)
)
q0_tilt <- bound_probability(
inv_logit(safe_logit(bias_dat$q0_true) + bias_surface$delta_q[i] * g_q)
)
bias_surface$Estimate[i] <- evaluate_aipw_with_nuisances(
bias_dat,
e_hat_s = e_tilt,
q1_hat_s = q1_tilt,
q0_hat_s = q0_tilt
)
}
bias_surface$Bias <- bias_surface$Estimate - sim_truth["psi"]
axis_slice <- subset(
bias_surface,
abs(delta_e) < 1e-12 | abs(delta_q) < 1e-12 | abs(delta_e - delta_q) < 1e-12
)
axis_slice$Abs_bias <- abs(axis_slice$Bias)
axis_slice_display <- axis_slice[order(axis_slice$delta_e, axis_slice$delta_q), ]
axis_slice_display$Estimate <- round(axis_slice_display$Estimate, 4)
axis_slice_display$Bias <- round(axis_slice_display$Bias, 4)
axis_slice_display$Abs_bias <- round(axis_slice_display$Abs_bias, 4)
names(axis_slice_display) <- c(
"\\(\\delta_e\\)",
"\\(\\delta_Q\\)",
"\\(\\hat\\psi_{DR}(\\delta_e,\\delta_Q)\\)",
"Bias",
"\\(\\lvert\\mathrm{Bias}\\rvert\\)"
)
mkable(head(axis_slice_display, 18))| \(\delta_e\) | \(\delta_Q\) | \(\hat\psi_{DR}(\delta_e,\delta_Q)\) | Bias | \(\lvert\mathrm{Bias}\rvert\) |
|---|---|---|---|---|
| -1.2 | -1.2 | 0.0706 | -0.0101 | 0.0101 |
| -1.2 | 0.0 | 0.0777 | -0.0030 | 0.0030 |
| -0.9 | -0.9 | 0.0711 | -0.0096 | 0.0096 |
| -0.9 | 0.0 | 0.0784 | -0.0023 | 0.0023 |
| -0.6 | -0.6 | 0.0742 | -0.0065 | 0.0065 |
| -0.6 | 0.0 | 0.0787 | -0.0020 | 0.0020 |
| -0.3 | -0.3 | 0.0773 | -0.0034 | 0.0034 |
| -0.3 | 0.0 | 0.0787 | -0.0020 | 0.0020 |
| 0.0 | -1.2 | 0.0787 | -0.0020 | 0.0020 |
| 0.0 | -0.9 | 0.0787 | -0.0020 | 0.0020 |
| 0.0 | -0.6 | 0.0786 | -0.0021 | 0.0021 |
| 0.0 | -0.3 | 0.0785 | -0.0022 | 0.0022 |
| 0.0 | 0.0 | 0.0784 | -0.0023 | 0.0023 |
| 0.0 | 0.3 | 0.0783 | -0.0024 | 0.0024 |
| 0.0 | 0.6 | 0.0784 | -0.0023 | 0.0023 |
| 0.0 | 0.9 | 0.0784 | -0.0022 | 0.0022 |
| 0.0 | 1.2 | 0.0786 | -0.0021 | 0.0021 |
| 0.3 | 0.0 | 0.0771 | -0.0036 | 0.0036 |
This plot is a useful diagnostic picture for Lecture 4. The DR estimator is not magic; it has a second-order remainder. The axes are protected by double robustness. The corners are not.
20 Simulation 3: Positivity and the Bound
This simulation asks a different question from the previous sections. We are no longer asking whether DR is unbiased or whether the oracle EIF reaches its bound. We are asking:
How much statistical information is available when treatment assignment is nearly deterministic for some covariate profiles?
Recall the ATE efficient influence function:
\[ D^*(O) = Q_1(X)-Q_0(X)-\psi + \frac{T}{e(X)}\{Y-Q_1(X)\} - \frac{1-T}{1-e(X)}\{Y-Q_0(X)\}. \]
Its variance is the semiparametric efficiency bound,
\[ V_{\mathrm{EIF}} = E\{D^*(O)^2\}. \]
For binary outcomes, this decomposes into
\[ V_{\mathrm{EIF}} = \underbrace{ \mathrm{Var}\{Q_1(X)-Q_0(X)\} }_{\text{effect heterogeneity}} + \underbrace{ E\left[ \frac{Q_1(X)\{1-Q_1(X)\}}{e(X)} + \frac{Q_0(X)\{1-Q_0(X)\}}{1-e(X)} \right] }_{\text{residual noise amplified by inverse propensity scores}}. \]
The important positivity terms are therefore
\[ \frac{1}{e(X)} \quad \text{and} \quad \frac{1}{1-e(X)}. \]
As the propensity score gets closer to 0 or 1, these inverse factors become large. That makes each observed outcome carry more noise for the counterfactual contrast. This is not a small-sample problem. It is part of the first-order information bound.
In the simulation, the outcome model \(Q_t(X)\) is held fixed. Only the treatment mechanism changes:
\[ e_a(X) = \mathrm{expit} \left[ a\{-0.25+0.70X_1-0.55X_2+0.35X_3\} \right]. \]
The scalar \(a\) is the propensity-score separation strength:
- Small \(a\): \(e_a(X)\) stays closer to 0.5, so treated and untreated subjects overlap well.
- Large \(a\): \(e_a(X)\) moves closer to 0 or 1 for many covariate profiles, so overlap becomes weak.
Because \(Q_1(X)\) and \(Q_0(X)\) do not change with \(a\), the target ATE
\[ \psi = E\{Q_1(X)-Q_0(X)\} \]
should remain approximately constant. What changes is the information bound \(V_{\mathrm{EIF}}\), especially its residual-noise component.
set.seed(20260706)
overlap_grid <- seq(0.25, 2.75, by = 0.25)
positivity_results <- do.call(rbind, lapply(overlap_grid, function(a) {
pop_a <- generate_dgp(n = 100000, overlap_strength = a)
psi_a <- mean(pop_a$q1_true - pop_a$q0_true)
heterogeneity_a <- var(pop_a$q1_true - pop_a$q0_true)
residual_a <- mean(
pop_a$q1_true * (1 - pop_a$q1_true) / pop_a$e_true +
pop_a$q0_true * (1 - pop_a$q0_true) / (1 - pop_a$e_true)
)
data.frame(
overlap_strength = a,
psi = psi_a,
v_eff = heterogeneity_a + residual_a,
heterogeneity = heterogeneity_a,
residual = residual_a,
p05_e = quantile(pop_a$e_true, 0.05),
p95_e = quantile(pop_a$e_true, 0.95),
mean_inverse_factor = mean(1 / pop_a$e_true + 1 / (1 - pop_a$e_true)),
row.names = NULL
)
}))
positivity_table <- data.frame(
Overlap_strength = positivity_results$overlap_strength,
Psi = round(positivity_results$psi, 4),
P05_e = round(positivity_results$p05_e, 3),
P95_e = round(positivity_results$p95_e, 3),
Mean_inverse_factor = round(positivity_results$mean_inverse_factor, 2),
V_eff = round(positivity_results$v_eff, 4),
Heterogeneity = round(positivity_results$heterogeneity, 4),
Residual = round(positivity_results$residual, 4)
)
names(positivity_table) <- c(
"Overlap strength",
"\\(\\psi\\)",
"5th pct. \\(e_a(X)\\)",
"95th pct. \\(e_a(X)\\)",
"\\(E[1/e_a+1/(1-e_a)]\\)",
"\\(V_{\\mathrm{EIF}}\\)",
"\\(\\mathrm{Var}\\{Q_1-Q_0\\}\\)",
"Residual component"
)
mkable(positivity_table)| Overlap strength | \(\psi\) | 5th pct. \(e_a(X)\) | 95th pct. \(e_a(X)\) | \(E[1/e_a+1/(1-e_a)]\) | \(V_{\mathrm{EIF}}\) | \(\mathrm{Var}\{Q_1-Q_0\}\) | Residual component |
|---|---|---|---|---|---|---|---|
| 0.25 | 0.0806 | 0.403 | 0.587 | 4.05 | 0.5408 | 0.0060 | 0.5348 |
| 0.50 | 0.0802 | 0.312 | 0.667 | 4.22 | 0.5582 | 0.0060 | 0.5522 |
| 0.75 | 0.0804 | 0.234 | 0.740 | 4.52 | 0.5943 | 0.0060 | 0.5883 |
| 1.00 | 0.0805 | 0.171 | 0.803 | 5.03 | 0.6522 | 0.0061 | 0.6461 |
| 1.25 | 0.0808 | 0.122 | 0.851 | 5.82 | 0.7383 | 0.0061 | 0.7322 |
| 1.50 | 0.0808 | 0.085 | 0.892 | 7.09 | 0.8785 | 0.0061 | 0.8724 |
| 1.75 | 0.0810 | 0.059 | 0.920 | 9.10 | 1.0890 | 0.0061 | 1.0829 |
| 2.00 | 0.0808 | 0.041 | 0.942 | 12.43 | 1.4122 | 0.0061 | 1.4061 |
| 2.25 | 0.0804 | 0.027 | 0.959 | 18.32 | 1.9706 | 0.0061 | 1.9645 |
| 2.50 | 0.0805 | 0.019 | 0.971 | 28.16 | 2.7967 | 0.0060 | 2.7906 |
| 2.75 | 0.0807 | 0.013 | 0.979 | 46.65 | 4.5981 | 0.0060 | 4.5921 |
The target effect barely changes in this simulation because the outcome model is fixed. The information bound changes sharply because treatment assignment becomes less overlapping.
The table and figure should be read in the following way:
- The 5th percentile of \(e_a(X)\) decreases and the 95th percentile increases as \(a\) grows. This means some covariate profiles are almost always untreated and others are almost always treated.
- The average inverse factor \(E[1/e_a(X)+1/\{1-e_a(X)\}]\) increases. This is a direct numerical summary of positivity stress.
- The heterogeneity component \(\mathrm{Var}\{Q_1(X)-Q_0(X)\}\) stays nearly fixed because \(Q_1,Q_0\), and the distribution of \(X\), are fixed.
- The residual component increases because the same Bernoulli outcome noise is multiplied by \(1/e_a(X)\) or \(1/\{1-e_a(X)\}\).
So the point is not that confounding bias got worse. In this population calculation there is no estimation bias at all. The point is that weak positivity reduces information: even the best regular estimator must have a larger asymptotic variance.
20.1 Simulation 3A: Positivity, Inverse Weights, and Coverage
The previous positivity experiment computes the population bound. Here we ask what the analyst sees in repeated samples.
For each overlap strength \(a\), we generate \(R\) datasets of size \(n=800\). In each dataset we estimate the DR one-step estimator and its EIF standard error,
\[ \widehat{SE}_{EIF} = \left[ \frac{1}{n} P_n\{\hat D_i-\bar{\hat D}\}^2 \right]^{1/2}, \]
where
\[ \hat D_i = \hat Q_1(X_i)-\hat Q_0(X_i) + \frac{T_i}{\hat e(X_i)}\{Y_i-\hat Q_1(X_i)\} - \frac{1-T_i}{1-\hat e(X_i)}\{Y_i-\hat Q_0(X_i)\}. \]
The Wald interval is
\[ \hat\psi_{DR} \pm 1.96\,\widehat{SE}_{EIF}. \]
We also record the largest inverse-probability factor in each sample:
\[ \max_i \left[ \max\left\{ \frac{1}{\hat e(X_i)}, \frac{1}{1-\hat e(X_i)} \right\} \right]. \]
This is a finite-sample warning sign. When this number is large, a small number of observations can dominate the EIF correction term.
As overlap worsens, three things tend to move together:
- The true lower-bound SE \(\sqrt{V_{\mathrm{EIF}}/n}\) increases.
- The empirical EIF SE increases.
- The largest inverse-probability weights become unstable.
set.seed(20260706)
overlap_coverage_grid <- c(0.50, 1.00, 1.75, 2.50)
R_overlap <- 80
n_overlap <- 800
overlap_results <- vector("list", length(overlap_coverage_grid))
for (j in seq_along(overlap_coverage_grid)) {
a_j <- overlap_coverage_grid[j]
truth_j <- truth_from_large_population(n_pop = 100000, overlap_strength = a_j)
current <- vector("list", R_overlap)
for (r in seq_len(R_overlap)) {
dat_r <- generate_dgp(n_overlap, overlap_strength = a_j)
dr_fit <- tryCatch(
estimate_dr_with_eif(dat_r, e_spec = "correct", q_spec = "correct"),
error = function(e) rep(NA_real_, 7)
)
current[[r]] <- data.frame(
Overlap_strength = a_j,
Truth = truth_j["psi"],
Bound_SE = sqrt(truth_j["v_eff"] / n_overlap),
Estimate = dr_fit["psi"],
SE = dr_fit["se"],
Lower = dr_fit["lower"],
Upper = dr_fit["upper"],
Min_e_hat = dr_fit["e_min"],
Max_e_hat = dr_fit["e_max"],
Max_inverse_weight = dr_fit["max_inverse_weight"]
)
}
overlap_results[[j]] <- do.call(rbind, current)
}
overlap_results <- do.call(rbind, overlap_results)
overlap_results <- overlap_results[complete.cases(overlap_results), ]
overlap_summary <- do.call(rbind, lapply(
split(overlap_results, overlap_results$Overlap_strength),
function(current) {
data.frame(
Overlap_strength = unique(current$Overlap_strength),
Truth = unique(current$Truth),
Bound_SE = unique(current$Bound_SE),
Bias = mean(current$Estimate) - unique(current$Truth),
Monte_Carlo_SD = sd(current$Estimate),
Mean_EIF_SE = mean(current$SE),
Coverage = mean(current$Lower <= current$Truth & current$Upper >= current$Truth),
Median_max_inverse_weight = median(current$Max_inverse_weight),
P95_max_inverse_weight = quantile(current$Max_inverse_weight, 0.95),
row.names = NULL
)
}
))
overlap_display <- overlap_summary
overlap_display[, -1] <- round(overlap_display[, -1], 3)
names(overlap_display) <- c(
"Overlap strength",
"\\(\\psi\\)",
"\\(\\sqrt{V_{\\mathrm{EIF}}/n}\\)",
"Bias",
"MC SD",
"Mean \\(\\widehat{SE}_{EIF}\\)",
"95% CI coverage",
"Median max inverse weight",
"95th pct. max inverse weight"
)
mkable(overlap_display)| Overlap strength | \(\psi\) | \(\sqrt{V_{\mathrm{EIF}}/n}\) | Bias | MC SD | Mean \(\widehat{SE}_{EIF}\) | 95% CI coverage | Median max inverse weight | 95th pct. max inverse weight |
|---|---|---|---|---|---|---|---|---|
| 0.50 | 0.081 | 0.026 | 0.005 | 0.026 | 0.026 | 0.975 | 5.983 | 9.077 |
| 1.00 | 0.081 | 0.029 | -0.006 | 0.029 | 0.028 | 0.950 | 23.177 | 54.540 |
| 1.75 | 0.080 | 0.037 | 0.001 | 0.037 | 0.037 | 0.938 | 231.512 | 854.732 |
| 2.50 | 0.081 | 0.061 | 0.000 | 0.053 | 0.045 | 0.975 | 2062.850 | 8177.446 |
The left panel compares three standard-error scales. The red curve is the population lower-bound standard error, \(\sqrt{V_{\mathrm{EIF}}/n}\). The blue curve is the Monte Carlo standard deviation of the repeated DR estimates. The green curve is the average EIF standard error computed by the analyst. Under stable overlap, these should be close. As overlap weakens, all three move upward because the EIF itself becomes more variable.
The right panel explains why. The maximum inverse weight increases as the propensity scores approach 0 or 1. This does not necessarily create bias when the models are correct, but it makes the estimating equation noisier. In a real data analysis, a large maximum inverse weight means the estimate is partly determined by a few observations in covariate regions with little treatment overlap.
This is why positivity is not only an identification assumption. It is also an information condition. Even with correct models and no unmeasured confounding, limited overlap can make the lower bound too large for precise estimation.
21 Simulation 4: The \(n^{-1/4}\) Product-Rate Story
Suppose the treatment model error is approximately \(n^{-\alpha_e}\) and the outcome model error is approximately \(n^{-\alpha_Q}\). The DR second-order remainder is roughly
\[ n^{-(\alpha_e+\alpha_Q)}. \]
For this to be negligible relative to \(n^{-1/2}\), we need
\[ \alpha_e+\alpha_Q > 1/2. \]
n_grid <- round(exp(seq(log(200), log(50000), length.out = 80)))
rate_cases <- data.frame(
Case = c(
"Both n^{-1/4}",
"One n^{-1/2}, one n^{-0.20}",
"Both too slow",
"One n^{-0.35}, one n^{-0.20}"
),
alpha_e = c(0.25, 0.50, 0.15, 0.35),
alpha_q = c(0.25, 0.20, 0.15, 0.20)
)
rate_curves <- lapply(seq_len(nrow(rate_cases)), function(i) {
data.frame(
n = n_grid,
Case = rate_cases$Case[i],
scaled_remainder = n_grid^(0.5 - rate_cases$alpha_e[i] - rate_cases$alpha_q[i])
)
})
rate_curves <- do.call(rbind, rate_curves)
rate_table <- rate_cases
rate_table$alpha_sum <- rate_table$alpha_e + rate_table$alpha_q
rate_table$Root_n_negligible <- ifelse(rate_table$alpha_sum > 0.5, "yes", "borderline/no")
rate_table$Case <- c(
"Both \\(n^{-1/4}\\)",
"One \\(n^{-1/2}\\), one \\(n^{-0.20}\\)",
"Both too slow",
"One \\(n^{-0.35}\\), one \\(n^{-0.20}\\)"
)
names(rate_table) <- c(
"Case",
"\\(\\alpha_e\\)",
"\\(\\alpha_Q\\)",
"\\(\\alpha_e+\\alpha_Q\\)",
"\\(\\sqrt n R_2=o_p(1)\\)?"
)
mkable(rate_table)| Case | \(\alpha_e\) | \(\alpha_Q\) | \(\alpha_e+\alpha_Q\) | \(\sqrt n R_2=o_p(1)\)? |
|---|---|---|---|---|
| Both \(n^{-1/4}\) | 0.25 | 0.25 | 0.50 | borderline/no |
| One \(n^{-1/2}\), one \(n^{-0.20}\) | 0.50 | 0.20 | 0.70 | yes |
| Both too slow | 0.15 | 0.15 | 0.30 | borderline/no |
| One \(n^{-0.35}\), one \(n^{-0.20}\) | 0.35 | 0.20 | 0.55 | yes |
If the curve decreases toward zero, the second-order nuisance remainder is asymptotically negligible on the \(\sqrt{n}\) scale. If it stays flat or increases, first-order inference can fail. The flat green curve is the boundary case \(\alpha_e+\alpha_Q=1/2\). The blue and orange curves both decrease because their exponent sums are larger than \(1/2\), but the blue curve decreases faster because \(0.50+0.20=0.70\), whereas \(0.35+0.20=0.55\).
21.1 Simulation 4A: A Population-Level Product-Rate Perturbation
The previous figure was algebraic. We can make it more concrete by tracking the leading product-remainder magnitude directly.
A population DR expansion has the schematic form
\[ \hat\psi_{DR}-\psi = (P_n-P)D^*(O) + R_2(\hat\eta,\eta) + \text{higher-order terms}, \]
where the second-order drift is approximately a product of nuisance errors:
\[ |R_2(\hat\eta,\eta)| \approx C \|\hat e-e\| \left(\|\hat Q_1-Q_1\|+\|\hat Q_0-Q_0\|\right). \]
If
\[ \|\hat e-e\| \asymp n^{-\alpha_e}, \qquad \|\hat Q_t-Q_t\| \asymp n^{-\alpha_Q}, \]
then
\[ |R_2(\hat\eta,\eta)| \asymp n^{-(\alpha_e+\alpha_Q)} \]
and the root-\(n\) scaled remainder behaves like
\[ \sqrt n\, |R_2(\hat\eta,\eta)| \asymp n^{1/2-(\alpha_e+\alpha_Q)}. \]
This figure plots that leading product magnitude. This is cleaner than plotting the signed bias from arbitrary nuisance perturbation directions, because signed bias can accidentally cross zero and create a misleading U-shape after taking absolute values.
rate_bias_n_grid <- round(exp(seq(log(250), log(30000), length.out = 25)))
rate_bias_cases <- data.frame(
Case = c(
"alpha_e + alpha_Q = 0.70",
"alpha_e + alpha_Q = 0.50",
"alpha_e + alpha_Q = 0.30"
),
alpha_e = c(0.35, 0.25, 0.15),
alpha_q = c(0.35, 0.25, 0.15),
stringsAsFactors = FALSE
)
rate_bias_results <- lapply(seq_len(nrow(rate_bias_cases)), function(i) {
current_case <- rate_bias_cases[i, ]
out <- lapply(rate_bias_n_grid, function(n_nominal) {
nuisance_product <- n_nominal^(
-(current_case$alpha_e + current_case$alpha_q)
)
data.frame(
n = n_nominal,
Case = current_case$Case,
alpha_sum = current_case$alpha_e + current_case$alpha_q,
Product_remainder = nuisance_product,
Root_n_scaled_remainder = sqrt(n_nominal) * nuisance_product
)
})
do.call(rbind, out)
})
rate_bias_results <- do.call(rbind, rate_bias_results)
rate_bias_display <- do.call(rbind, lapply(
split(rate_bias_results, rate_bias_results$Case),
function(current) {
current[c(1, round(nrow(current) / 2), nrow(current)), ]
}
))
row.names(rate_bias_display) <- NULL
rate_bias_display$Case <- paste0(
"\\(\\alpha_e+\\alpha_Q = ",
sprintf("%.2f", rate_bias_display$alpha_sum),
"\\)"
)
rate_bias_display$Product_remainder <- signif(rate_bias_display$Product_remainder, 3)
rate_bias_display$Root_n_scaled_remainder <- signif(rate_bias_display$Root_n_scaled_remainder, 3)
names(rate_bias_display) <- c(
"\\(n\\)",
"Case",
"\\(\\alpha_e+\\alpha_Q\\)",
"\\(n^{-(\\alpha_e+\\alpha_Q)}\\)",
"\\(\\sqrt n\\,n^{-(\\alpha_e+\\alpha_Q)}\\)"
)
mkable(rate_bias_display)| \(n\) | Case | \(\alpha_e+\alpha_Q\) | \(n^{-(\alpha_e+\alpha_Q)}\) | \(\sqrt n\,n^{-(\alpha_e+\alpha_Q)}\) |
|---|---|---|---|---|
| 250 | \(\alpha_e+\alpha_Q = 0.30\) | 0.3 | 0.191000 | 3.020 |
| 2243 | \(\alpha_e+\alpha_Q = 0.30\) | 0.3 | 0.098800 | 4.680 |
| 30000 | \(\alpha_e+\alpha_Q = 0.30\) | 0.3 | 0.045400 | 7.860 |
| 250 | \(\alpha_e+\alpha_Q = 0.50\) | 0.5 | 0.063200 | 1.000 |
| 2243 | \(\alpha_e+\alpha_Q = 0.50\) | 0.5 | 0.021100 | 1.000 |
| 30000 | \(\alpha_e+\alpha_Q = 0.50\) | 0.5 | 0.005770 | 1.000 |
| 250 | \(\alpha_e+\alpha_Q = 0.70\) | 0.7 | 0.021000 | 0.331 |
| 2243 | \(\alpha_e+\alpha_Q = 0.70\) | 0.7 | 0.004510 | 0.214 |
| 30000 | \(\alpha_e+\alpha_Q = 0.70\) | 0.7 | 0.000735 | 0.127 |
The green curve represents a nuisance product rate faster than \(n^{-1/2}\), so the root-\(n\) scaled remainder vanishes. The orange curve is the boundary case \(\alpha_e+\alpha_Q=1/2\), so the root-\(n\) scaled remainder stays constant. The red curve represents nuisance learning that is too slow; the second-order term remains visible, and even grows, after root-\(n\) scaling.
22 Cross-Fitting and DML
When nuisance functions are estimated with flexible machine learning, empirical process conditions such as Donsker restrictions can be hard to verify.
Cross-fitting is the practical workaround:
- Split the sample into \(K\) folds.
- For each fold, train \(\hat e\) and \(\hat Q\) on the other folds.
- Predict nuisance functions in the held-out fold.
- Evaluate the EIF estimating function using held-out predictions.
- Average the fold-specific EIF contributions.
Cross-fitting helps remove “own-observation” bias: the same outcome \(Y_i\) is not used both to train the nuisance function and to evaluate its residual correction for observation \(i\).
In DML language, the efficient influence function is the orthogonal score.
22.1 Simulation 5: Cross-Fitting with Rich Nuisance Models
This simulation has a specific role in the lecture. Up to this point, we have studied the EIF as if the nuisance functions \(e(X)\), \(Q_1(X)\), and \(Q_0(X)\) were either known or estimated by simple parametric models. In modern semiparametric inference, however, nuisance functions are often learned by flexible regression methods. Then the main technical danger is that the same observation can be used twice:
- once to train \(\hat e\) or \(\hat Q\);
- again to evaluate the EIF residual for that same observation.
That reuse can create own-observation feedback. The fitted nuisance function may partially adapt to \(Y_i\) or \(T_i\), making the residual correction look artificially stable.
Cross-fitting changes the construction of the EIF contribution. Let \(I_1,\ldots,I_K\) be a partition of the sample. For observation \(i\in I_k\), estimate the nuisance functions using only the training sample \(I_{-k}\):
\[ \hat e^{(-k)}(X_i), \qquad \hat Q_1^{(-k)}(X_i), \qquad \hat Q_0^{(-k)}(X_i). \]
Then compute the held-out one-step contribution
\[ \hat\phi_i^{CF} = \hat Q_1^{(-k)}(X_i)-\hat Q_0^{(-k)}(X_i) + \frac{T_i}{\hat e^{(-k)}(X_i)} \{Y_i-\hat Q_1^{(-k)}(X_i)\} - \frac{1-T_i}{1-\hat e^{(-k)}(X_i)} \{Y_i-\hat Q_0^{(-k)}(X_i)\}. \]
The cross-fitted DR estimate is
\[ \hat\psi_{CF} = P_n \hat\phi_i^{CF}. \]
The centered estimated EIF contribution used for standard errors is then
\[ \hat D_i^{CF} = \hat\phi_i^{CF}-\hat\psi_{CF}. \]
This is not a different causal estimand. It is a different way to build the same EIF-based estimating equation so that nuisance learning is separated from EIF evaluation. That separation is why DML theory can allow rich nuisance learners under product-rate conditions rather than relying only on classical Donsker empirical-process assumptions.
To see the mechanics, we fit richer nuisance models than the true DGP requires: cubic terms and an \(X_1X_2\) interaction. This is still a GLM, not a full machine-learning learner, but it creates a useful classroom version of the DML workflow.
We compare:
- A rich in-sample DR estimator.
- A rich 5-fold cross-fitted DR estimator.
- The infeasible oracle EIF estimator.
set.seed(20260706)
R_cf <- 60
n_cf <- 600
K_cf <- 5
crossfit_results <- vector("list", 3 * R_cf)
counter <- 1
for (r in seq_len(R_cf)) {
dat_r <- generate_dgp(n_cf)
insample_fit <- tryCatch(
estimate_dr_with_formula(dat_r, rich_e_formula, rich_q_formula),
error = function(e) rep(NA_real_, 4)
)
crossfit_fit <- tryCatch(
estimate_dr_crossfit(dat_r, rich_e_formula, rich_q_formula, K = K_cf),
error = function(e) rep(NA_real_, 4)
)
oracle_fit <- oracle_eif_with_se(dat_r)
fit_matrix <- rbind(
`Rich in-sample DR` = insample_fit,
`Rich cross-fitted DR` = crossfit_fit,
`Oracle EIF` = oracle_fit
)
for (m in rownames(fit_matrix)) {
crossfit_results[[counter]] <- data.frame(
n = n_cf,
Method = m,
Estimate = fit_matrix[m, "psi"],
SE = fit_matrix[m, "se"],
Lower = fit_matrix[m, "lower"],
Upper = fit_matrix[m, "upper"]
)
counter <- counter + 1
}
}
crossfit_results <- do.call(rbind, crossfit_results)
crossfit_results <- crossfit_results[complete.cases(crossfit_results), ]
crossfit_summary <- summarize_ci_performance(
crossfit_results,
truth = sim_truth["psi"],
v_eff = sim_truth["v_eff"]
)
crossfit_display <- crossfit_summary
crossfit_display[, c("Bias", "Monte_Carlo_SD", "Mean_EIF_SE", "Coverage", "Efficiency_ratio")] <-
round(crossfit_display[, c("Bias", "Monte_Carlo_SD", "Mean_EIF_SE", "Coverage", "Efficiency_ratio")], 3)
crossfit_display$Method <- method_math_label(crossfit_display$Method)
names(crossfit_display) <- c(
"\\(n\\)",
"Estimator",
"Bias",
"MC SD",
"Mean \\(\\widehat{SE}_{EIF}\\)",
"95% CI coverage",
"\\(\\mathrm{Var}_{MC}(\\hat\\psi)/(V_{\\mathrm{EIF}}/n)\\)"
)
mkable(crossfit_display)| \(n\) | Estimator | Bias | MC SD | Mean \(\widehat{SE}_{EIF}\) | 95% CI coverage | \(\mathrm{Var}_{MC}(\hat\psi)/(V_{\mathrm{EIF}}/n)\) |
|---|---|---|---|---|---|---|
| 600 | \(\hat\psi_{oracle\ EIF}\) | -0.001 | 0.035 | 0.033 | 0.917 | 1.131 |
| 600 | \(\hat\psi_{DR,cf}\) | 0.051 | 0.544 | 0.136 | 0.950 | 272.333 |
| 600 | \(\hat\psi_{DR,rich}\) | -0.002 | 0.035 | 0.032 | 0.917 | 1.149 |
The table compares the repeated-sampling behavior of three estimators:
- \(\hat\psi_{\mathrm{oracle}}\): uses the true nuisance functions. This is the unattainable benchmark.
- \(\hat\psi_{\mathrm{CF}}\): estimates nuisances on training folds and evaluates EIF contributions on held-out folds.
- \(\hat\psi_{\mathrm{in}}\): estimates nuisances and evaluates the EIF on the same observations.
The comparison is not meant to prove that cross-fitting always dominates in every small simulation. With simple GLMs, the difference can be modest because GLMs are relatively stable. The point is structural: cross-fitting creates the empirical version of the theoretical condition that the nuisance error is separated from the score evaluation.
The DML paper often visualizes this idea with histograms: a conventional or non-orthogonal construction can be shifted or poorly matched to its Gaussian approximation, while an orthogonal score is much closer to the reference normal distribution. The exact partially linear model in that paper is different from our ATE setting, so the analogous diagnostic here is:
- same-sample EIF evaluation: fit nuisance functions and evaluate the score on the same observations;
- cross-fitted orthogonal EIF evaluation: fit nuisance functions away from observation \(i\), then evaluate the EIF contribution for observation \(i\).
The cleanest scale for this comparison is the studentized repeated-sampling error
\[ Z_r = \frac{\hat\psi^{(r)}-\psi}{\widehat{SE}^{(r)}}. \]
If the first-order EIF approximation is working well, the histogram of \(Z_r\) should be close to the red \(N(0,1)\) curve.
The point of this figure is not that every small Monte Carlo run must look perfectly normal. Here \(R=60\), so the histograms are necessarily rough. The main diagnostic is whether the studentized errors are centered near zero and have roughly unit spread. Cross-fitting is designed to make this approximation more credible by removing the direct path from \(Y_i\) into both nuisance training and EIF evaluation for the same observation.
The boxplot displays estimation error in percentage points over repeated simulations. The oracle EIF distribution is the benchmark target. The cross-fitted estimator should be read as the feasible estimator that most closely follows the DML construction. The in-sample estimator can perform well with stable parametric regressions, but it is not the construction that protects us from adaptive overfitting.
The role of this simulation is therefore conceptual and technical:
- It shows how the EIF is evaluated after nuisance learning.
- It explains why DML uses out-of-fold nuisance predictions.
- It connects the product-rate condition to an implementable algorithm.
- It clarifies that cross-fitting is not a new estimand; it is a device for making the same EIF-based estimator behave like its asymptotic expansion.
23 Practical Checklist
When using an EIF-based lower bound for the smoking-CVD risk difference, check:
- Estimand: Is the target \(\psi=E\{Y(1)\}-E\{Y(0)\}\), a risk ratio, or another scale?
- Observed-data model: Are we using the same i.i.d. complete-case model assumed by this lecture note?
- Identification: Are consistency, exchangeability, and positivity plausible?
- Nuisance functions: Are \(e(X)\), \(Q_1(X)\), and \(Q_0(X)\) estimated well enough?
- Overlap: Are estimated propensities bounded away from 0 and 1?
- Estimator: Does the estimator solve or approximate the EIF equation?
- Inference: Are SEs based on the EIF, bootstrap, or both?
- Sensitivity: What does Lecture 6 say about unmeasured confounding?
24 Key Takeaways
- The DR estimator from Lecture 4 is the sample analog of the efficient influence function for the risk difference.
- The semiparametric efficiency lower bound is \(E\{D^*(O)^2\}/n\).
- The bound decomposes into treatment-effect heterogeneity plus residual outcome noise amplified by inverse propensity scores.
- Weak positivity inflates the bound even if there is no bias.
- DR estimators can be consistent if either the treatment model or outcome model is correct, but they are efficient only when the full first-order expansion is controlled.
- The \(n^{-1/4}\) product-rate condition explains why flexible nuisance estimation can still yield root-\(n\) inference.
- The simulation experiments show the same theory from several angles: efficiency ratios, Wald coverage, bias surfaces, positivity stress tests, and cross-fitted nuisance learning.
- Cross-fitting separates nuisance learning from EIF evaluation and makes the DR estimating equation usable with richer learners.
25 Appendix 1: Pathwise Derivative Sketch
Let \(P_\epsilon\) be a regular parametric submodel through \(P\) with score
\[ S(O) = \left. \frac{\partial}{\partial \epsilon} \log p_\epsilon(O) \right|_{\epsilon=0}. \]
A parameter \(\psi(P)\) is pathwise differentiable if
\[ \left. \frac{\partial}{\partial \epsilon} \psi(P_\epsilon) \right|_{\epsilon=0} = E\{D^*(O)S(O)\} \]
for every regular submodel.
For
\[ \psi_t(P) = E\{Q_t(X)\}, \]
the candidate gradient is
\[ D_t^*(O) = \frac{I(T=t)}{P(T=t \mid X)} \{Y-Q_t(X)\} + Q_t(X)-\psi_t. \]
One can verify:
- \(E\{D_t^*(O)\}=0\).
- The weighted residual term captures perturbations of the conditional outcome law \(Y \mid T,X\).
- The \(Q_t(X)-\psi_t\) term captures perturbations of the covariate distribution \(P_X\).
- The treatment mechanism perturbation contributes zero because the residual term has conditional mean zero given \(T,X\).
Thus \(D_t^*\) represents the pathwise derivative and is the canonical gradient in the nonparametric model.
26 Appendix 2: Mean-Zero Verification
For \(t=1\):
\[ E\left[ \frac{T}{e(X)}\{Y-Q_1(X)\} \mid X \right] = \frac{E(T\mid X)}{e(X)} E\{Y-Q_1(X)\mid T=1,X\} = 0. \]
Also,
\[ E\{Q_1(X)-\psi_1\}=0. \]
Therefore \(E\{D_1^*(O)\}=0\). The proof for \(D_0^*\) is analogous.
27 Appendix 3: Variance Decomposition
Write
\[ D^*(O) = R_1(O)-R_0(O)+\tau(X)-\psi, \]
where
\[ R_1(O)=\frac{T}{e(X)}\{Y-Q_1(X)\}, \]
\[ R_0(O)=\frac{1-T}{1-e(X)}\{Y-Q_0(X)\}, \]
and
\[ \tau(X)=Q_1(X)-Q_0(X). \]
Conditional on \(X\), \(R_1\) and \(R_0\) have mean zero. Also \(R_1R_0=0\) because \(T\) cannot equal both 1 and 0. Therefore cross-terms vanish, and
\[ E(D^{*2}) = E\{R_1^2\} + E\{R_0^2\} + \mathrm{Var}\{\tau(X)\}. \]
Now
\[ E(R_1^2 \mid X) = \frac{1}{e(X)^2} E[T\{Y-Q_1(X)\}^2 \mid X] = \frac{\mathrm{Var}(Y \mid T=1,X)}{e(X)}. \]
Similarly,
\[ E(R_0^2 \mid X) = \frac{\mathrm{Var}(Y \mid T=0,X)}{1-e(X)}. \]
This gives the decomposition used above.
28 Appendix 4: DR Remainder Sketch
Let
\[ m(O;\eta) = Q_1(X)-Q_0(X) + \frac{T}{e(X)}\{Y-Q_1(X)\} - \frac{1-T}{1-e(X)}\{Y-Q_0(X)\}. \]
The DR estimator is \(P_n m(O;\hat\eta)\). Add and subtract terms:
\[ \hat\psi_{DR}-\psi = (P_n-P)m(O;\eta) + (P_n-P)\{m(O;\hat\eta)-m(O;\eta)\} + P\{m(O;\hat\eta)-m(O;\eta)\}. \]
The first term gives the efficient influence function. The second term is an empirical process term controlled by Donsker conditions or cross-fitting plus regularity. The third term is the drift or second-order remainder.
Under the true distribution,
\[ P\{m(O;\hat\eta)-m(O;\eta)\} \]
is a product of treatment-model and outcome-model errors. This is the analytic source of double robustness and the \(n^{-1/4}\) product-rate rule.
29 Appendix 5: Relation to Lecture 5 Bootstrap
If
\[ \hat\psi-\psi = P_nD^*(O)+o_p(n^{-1/2}), \]
then the first-order variance is
\[ \frac{E(D^{*2})}{n}. \]
The nonparametric bootstrap approximates the sampling distribution by resampling observations and rerunning the estimator. When the estimator is regular and the bootstrap is valid, the bootstrap SE estimates the same first-order target.
For highly adaptive, non-smooth, or design-based estimators, bootstrap validity may require additional care. This is why influence-function-based inference and bootstrap inference should be understood as related but not interchangeable.
30 References
Bang, H., & Robins, J. M. (2005). Doubly robust estimation in missing data and causal inference models. Biometrics, 61(4), 962-973.
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), C1-C68.
Hahn, J. (1998). On the role of the propensity score in efficient semiparametric estimation of average treatment effects. Econometrica, 66(2), 315-331.
Hernan, M. A., & Robins, J. M. (2020). Causal Inference: What If. Chapman & Hall/CRC.
Hirano, K., Imbens, G. W., & Ridder, G. (2003). Efficient estimation of average treatment effects using the estimated propensity score. Econometrica, 71(4), 1161-1189.
Kennedy, E. H. (2016). Semiparametric theory and empirical processes in causal inference. arXiv:1510.04740.
Newey, W. K. (1994). The asymptotic variance of semiparametric estimators. Econometrica, 62(6), 1349-1382.
Newey, W. K., & Robins, J. R. (2018). Cross-fitting and fast remainder rates for semiparametric estimation. arXiv:1801.09138.
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.
Tsiatis, A. A. (2006). Semiparametric Theory and Missing Data. Springer.
van der Laan, M. J., & Robins, J. M. (2003). Unified Methods for Censored Longitudinal Data and Causality. Springer.
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.
Velez, A. On the Asymptotic Properties of Debiased Machine Learning Estimators.
Wager, S. Causal Inference lecture notes.