1 Why Nuisance Misspecification Matters

1.1 Today’s Question

Many causal estimators depend on functions that are not themselves the final scientific target. For example, an estimator of a causal risk difference may require:

  • a model for treatment assignment; and
  • a model for the conditional outcome.

These supporting functions are called nuisance functions. The word nuisance means that they are not the target parameter; it does not mean that they are unimportant.

This lecture asks:

How accurately must nuisance functions be learned, and why can a doubly robust estimator tolerate some nuisance error?

Learning Goals

By the end of the lecture, you should be able to distinguish estimation error from permanent model misspecification, read parametric and nonparametric convergence rates, derive the product-of-errors remainder for a doubly robust estimator, and explain why \(n^{-1/4}\) is a critical balanced benchmark rather than a universal optimal nuisance rate.

1.2 Target Parameters and Nuisance Functions

For person \(i=1,\ldots,n\), let the observed record be

\[ O_i=(Y_i,T_i,X_i)\stackrel{\mathrm{iid}}{\sim}P, \]

where \(X_i\) is a vector of baseline covariates, \(T_i\in\{0,1\}\) is treatment or exposure, and \(Y_i\) is an outcome. For the smoking-CVD example, \(T_i=1\) denotes smoking and \(Y_i=1\) denotes observed CVD. A generic observation is written \(O=(Y,T,X)\).

Two central nuisance functions are

\[ e(x) = P(T=1\mid X=x) \]

and

\[ m(t,x) = E(Y\mid T=t,X=x), \qquad t\in\{0,1\}. \]

For binary \(Y\), this conditional mean is also \(P(Y=1\mid T=t,X=x)\). We use the course shorthand \(m_t(x)=m(t,x)\). Let \(Y_i(t)\) denote person \(i\)’s potential outcome under treatment level \(t\). Under consistency, conditional exchangeability, and positivity, the treatment-specific risks satisfy

\[ \psi_t=E\{Y(t)\}=E\{m_t(X)\},\qquad t\in\{0,1\}, \]

and the causal risk difference is

\[ RD=\psi_1-\psi_0=E\{Y(1)-Y(0)\}=E\{m_1(X)-m_0(X)\}. \]

The nuisance vector is

\[ \eta=(e,m_0,m_1). \]

Its fitted counterpart is \(\hat\eta=(\hat e,\hat m_0,\hat m_1)\).

Nuisance estimation cannot repair a failure of causal identification. Even perfectly learned \(e\) and \(m_t\) do not eliminate bias from unmeasured confounding, ill-defined interventions, interference, selection bias, or structural positivity failure.

2 What Does Misspecification Mean?

2.1 Truth, Working Class, Pseudo-Truth, and Estimate

Let \(f_0\in L_2(P_X)\) denote a generic true conditional function, such as \(e\) or \(m_t\), where \(P_X\) is the marginal distribution of \(X\) under \(P\). For a measurable function \(g\), define

\[ \|g\|_{2,P_X} = \left\{ \int g(x)^2\,dP_X(x) \right\}^{1/2}. \]

Suppose estimation is restricted to a nonempty working class \(\mathcal F\subset L_2(P_X)\). For the squared-\(L_2(P_X)\) geometry used in this section, define the population best approximation by

\[ f^\dagger\in\arg\min_{g\in\mathcal F}\|g-f_0\|_{2,P_X}^2. \]

The function \(f^\dagger\) is the pseudo-truth: the best member of the chosen class under the stated criterion, whether or not \(\mathcal F\) contains \(f_0\). A nonempty closed convex class has a unique \(L_2(P_X)\) projection; otherwise we explicitly assume that a measurable minimizer has been selected. The fitted rule \(\hat f\) is estimated from a finite training sample.

For a different learning loss, one may instead define \(f_\ell^\dagger\in\arg\min_{g\in\mathcal F}E_P[\ell\{O,g\}]\). That loss-based pseudo-truth agrees with the projection above only when the loss and target measure induce the same \(L_2(P_X)\) risk. The triangle inequality below remains valid for any selected \(f^\dagger\in L_2(P_X)\), whereas the later Pythagorean identity requires the stated least-squares projection structure.

The triangle inequality gives

\[ \|\hat f-f_0\|_{2,P_X} \le \underbrace{\|\hat f-f^\dagger\|_{2,P_X}}_{ \text{estimation and optimization error}} + \underbrace{\|f^\dagger-f_0\|_{2,P_X}}_{ \text{approximation or misspecification error}}. \]

Why there are no squares here. This is the triangle inequality for the \(L_2(P_X)\) norm, so the unsquared display is the generally valid statement. Squaring it gives

\[ \|\hat f-f_0\|_{2,P_X}^2 \le \left\{\|\hat f-f^\dagger\|_{2,P_X} +\|f^\dagger-f_0\|_{2,P_X}\right\}^2, \]

which contains a cross term and is not generally equal to the sum of the two squared errors. In the special case where \(f^\dagger\) is the least-squares \(L_2(P_X)\) projection onto a linear subspace and both \(f^\dagger\) and \(\hat f\) lie in that subspace,

\[ \langle \hat f-f^\dagger,\,f^\dagger-f_0\rangle_{P_X}=0, \]

so the Pythagorean identity applies:

\[ \|\hat f-f_0\|_{2,P_X}^2 = \|\hat f-f^\dagger\|_{2,P_X}^2 + \|f^\dagger-f_0\|_{2,P_X}^2. \]

How to read the geometry. This is a two-dimensional schematic cross-section of the Hilbert space \(L_2(P_X)\); the shaded plane represents the closed linear working subspace \(\mathcal F\). Its distances and right angle are induced by the \(L_2(P_X)\) inner product, not by pointwise function differences. The pseudo-truth \(f^\dagger\) is the \(L_2(P_X)\) projection of \(f_0\) onto \(\mathcal F\), and \(\hat f\in\mathcal F\). Therefore the red approximation component is orthogonal to every direction in the plane, including the blue estimation component \(\hat f-f^\dagger\). The purple diagonal is the total error \(\hat f-f_0\). This projection geometry permits the squared-error Pythagorean identity; without these conditions, only the triangle inequality is guaranteed.

The two errors answer different questions:

  • Estimation error: How far is the fitted rule from the best rule available in the chosen class?
  • Approximation error: How far is the best rule in that class from the true function?

With a fixed misspecified class, the first term can shrink rapidly while the second remains permanently positive.

2.2 Pictorial Illustration

Consider

\[ f_0(x)=1+1.25x+1.5x^2, \qquad X\sim\operatorname{Uniform}(-1,1). \]

If the working class contains only linear functions \(a+bx\), let \(q(x)=(1,x)^\top\). The least-squares projection coefficient solves the population normal equations,

\[ \begin{aligned} \theta^\dagger &=\left[E\{q(X)q(X)^\top\}\right]^{-1}E\{q(X)f_0(X)\} \\ &= \begin{pmatrix}1&0\\0&1/3\end{pmatrix}^{-1} \begin{pmatrix}3/2\\5/12\end{pmatrix} = \begin{pmatrix}3/2\\5/4\end{pmatrix}. \end{aligned} \]

Thus, the \(L_2(P_X)\) projection is

\[ f^\dagger(x)=1.5+1.25x. \]

The pointwise approximation difference is

\[ f^\dagger(x)-f_0(x)=\frac12-\frac32x^2. \]

For \(X\sim\operatorname{Uniform}(-1,1)\),

\[ E(X^{2k})=\frac12\int_{-1}^{1}x^{2k}\,dx=\frac{1}{2k+1}, \qquad E(X^2)=\frac13, \qquad E(X^4)=\frac15. \]

Therefore, the squared approximation error is

\[ \begin{aligned} \Delta_{\mathrm{app}}^2 &:=\|f^\dagger-f_0\|_{2,P_X}^2 \\ &=E\left(\frac12-\frac32X^2\right)^2 \\ &=\frac14-\frac32E(X^2)+\frac94E(X^4) \\ &=\frac14-\frac32\left(\frac13\right) +\frac94\left(\frac15\right) \\ &=\frac15. \end{aligned} \]

Taking the square root gives the permanent approximation floor

\[ \Delta_{\mathrm{app}} =\|f^\dagger-f_0\|_{2,P_X} =\frac{1}{\sqrt5} \approx0.447. \]

Single-sample interpretation. The blue fitted line can move closer to the red pseudo-truth as the sample grows. It does not move toward the black truth because the fixed linear class cannot represent curvature. Fast convergence to \(f^\dagger\) is not consistency for \(f_0\). The two vertical brackets are pointwise gaps at \(x=0.72\); the simulation below instead summarizes integrated \(L_2(P_X)\) function errors.

2.2.1 Monte Carlo Confirmation

For this least-squares projection, \(f_0-f^\dagger\) is orthogonal to every linear function. Therefore, every fitted linear rule satisfies the exact Pythagorean identity

\[ \|\hat f-f_0\|_{2,P_X}^2 =\|\hat f-f^\dagger\|_{2,P_X}^2 +\underbrace{\|f^\dagger-f_0\|_{2,P_X}^2}_{\Delta_{\mathrm{app}}^2=1/5}. \]

Consequently, if \(\|\hat f-f^\dagger\|_{2,P_X}\to0\), then \(\|\hat f-f_0\|_{2,P_X}\to\Delta_{\mathrm{app}}=1/\sqrt5\). The error relative to the truth cannot fall below this value asymptotically, which is why it is called the approximation floor.

We now generate \(X_i\sim\operatorname{Uniform}(-1,1)\) and \(Y_i=f_0(X_i)+\varepsilon_i\), where \(\varepsilon_i\sim N(0,0.7^2)\), fit the same misspecified linear regression in each sample, and repeat the experiment 120 times at every sample size.

For this data-generating process, random-design least-squares theory gives

\[ nE\|\hat f-f^\dagger\|_{2,P_X}^2\longrightarrow\frac{523}{350}, \qquad \left\{E\|\hat f-f^\dagger\|_{2,P_X}^2\right\}^{1/2} \sim \sqrt{\frac{523}{350}}\,n^{-1/2}\approx1.222\,n^{-1/2}. \]

Let \(B=120\) be the number of Monte Carlo repetitions and write the fitted rule in repetition \(b\) as

\[ \hat f_n^{(b)}(x)=\hat a_n^{(b)}+\hat b_n^{(b)}x, \qquad b=1,\ldots,B. \]

The table uses the following notation:

\[ \overline{\hat a}_n=\frac1B\sum_{b=1}^B\hat a_n^{(b)}, \qquad \overline{\hat b}_n=\frac1B\sum_{b=1}^B\hat b_n^{(b)}, \]

For either reference function \(g\in\{f^\dagger,f_0\}\), define the root mean integrated squared error

\[ \operatorname{RMISE}_n(g) =\left\{\frac1B\sum_{b=1}^B \|\hat f_n^{(b)}-g\|_{2,P_X}^2\right\}^{1/2}. \]

simulate_pseudo_truth_errors <- function(n) {
  x <- runif(n, -1, 1)
  y <- 1 + 1.25 * x + 1.5 * x^2 + rnorm(n, sd = 0.70)
  fit <- lm(y ~ x)
  fitted_coefficients <- unname(coef(fit))

  if (length(fitted_coefficients) != 2L ||
      any(!is.finite(fitted_coefficients))) {
    stop("The misspecified linear fit produced invalid coefficients.")
  }

  estimation_error_squared <-
    (fitted_coefficients[1] - 1.5)^2 +
      (fitted_coefficients[2] - 1.25)^2 / 3
  total_error_squared <- estimation_error_squared + 1 / 5

  c(
    `Estimated intercept` = fitted_coefficients[1],
    `Estimated slope` = fitted_coefficients[2],
    `Estimation error squared` = estimation_error_squared,
    `Total error squared` = total_error_squared
  )
}

set.seed(20260719)
pseudo_truth_sample_sizes <- c(50L, 100L, 250L, 500L, 1000L, 2500L, 5000L)
pseudo_truth_results <- lapply(
  pseudo_truth_sample_sizes,
  function(n) {
    replicate(rate_simulation_reps, simulate_pseudo_truth_errors(n))
  }
)

pseudo_truth_summary <- do.call(
  rbind,
  Map(
    function(n, errors) {
      estimation_error <- sqrt(errors["Estimation error squared", ])
      total_error <- sqrt(errors["Total error squared", ])
      rm_estimation_error <- sqrt(
        mean(errors["Estimation error squared", ])
      )
      rm_total_error <- sqrt(mean(errors["Total error squared", ]))

      data.frame(
        `Sample size` = n,
        `Mean intercept` = mean(errors["Estimated intercept", ]),
        `Mean slope` = mean(errors["Estimated slope", ]),
        `RM error to pseudo-truth` = rm_estimation_error,
        `Scaled RM error` = sqrt(n) * rm_estimation_error,
        `RM error to truth` = rm_total_error,
        `Approximation floor` = 1 / sqrt(5),
        Estimation_q10 = unname(quantile(estimation_error, 0.10)),
        Estimation_q90 = unname(quantile(estimation_error, 0.90)),
        Total_q10 = unname(quantile(total_error, 0.10)),
        Total_q90 = unname(quantile(total_error, 0.90)),
        check.names = FALSE
      )
    },
    pseudo_truth_sample_sizes,
    pseudo_truth_results
  )
)

if (any(!is.finite(as.matrix(pseudo_truth_summary)))) {
  stop("The pseudo-truth simulation produced a non-finite result.")
}

pseudo_truth_display <- pseudo_truth_summary[
  ,
  c(
    "Sample size",
    "Mean intercept",
    "Mean slope",
    "RM error to pseudo-truth",
    "Scaled RM error",
    "RM error to truth",
    "Approximation floor"
  )
]
names(pseudo_truth_display) <- c(
  "$n$",
  "$\\overline{\\hat a}_n$",
  "$\\overline{\\hat b}_n$",
  "$\\operatorname{RMISE}_n(f^\\dagger)$",
  "$\\sqrt n\\,\\operatorname{RMISE}_n(f^\\dagger)$",
  "$\\operatorname{RMISE}_n(f_0)$",
  "$\\Delta_{\\mathrm{app}}$"
)
for (column in names(pseudo_truth_display)[-1]) {
  pseudo_truth_display[[column]] <- sprintf(
    "%.3f",
    pseudo_truth_display[[column]]
  )
}

knitr::kable(
  pseudo_truth_display,
  align = rep("r", ncol(pseudo_truth_display)),
  escape = FALSE
)
\(n\) \(\overline{\hat a}_n\) \(\overline{\hat b}_n\) \(\operatorname{RMISE}_n(f^\dagger)\) \(\sqrt n\,\operatorname{RMISE}_n(f^\dagger)\) \(\operatorname{RMISE}_n(f_0)\) \(\Delta_{\mathrm{app}}\)
50 1.492 1.248 0.166 1.174 0.477 0.447
100 1.501 1.260 0.134 1.342 0.467 0.447
250 1.488 1.249 0.076 1.209 0.454 0.447
500 1.500 1.248 0.053 1.180 0.450 0.447
1000 1.498 1.250 0.039 1.219 0.449 0.447
2500 1.500 1.251 0.024 1.225 0.448 0.447
5000 1.500 1.253 0.016 1.111 0.447 0.447

Simulation interpretation. The root-mean integrated function error relative to \(f^\dagger\) decreases from 0.166 at \(n=50\) to 0.016 at \(n=5000\). Meanwhile, the root-mean integrated function error relative to \(f_0\) approaches 0.447, close to the permanent approximation floor \(1/\sqrt 5\approx0.447\). The approximately stable \(\sqrt n\)-scaled column is empirical evidence of root-\(n\) convergence to \(f^\dagger\), not convergence to \(f_0\). The shaded regions show the 10th–90th percentiles of the repetition-specific integrated function errors. These are integrated function errors; they are not outcome-prediction RMSEs, which would also contain the irreducible noise variance \(0.7^2\).

2.3 Misspecification Is Not the Same as Other Problems

Problem Meaning Does more data automatically solve it?
Sampling error The fitted rule varies across samples Usually decreases
Optimization error The fitting algorithm has not found its intended solution Only with better optimization
Overfitting The fitted rule adapts too closely to training noise Requires regularization or honest evaluation
Fixed-class misspecification The truth is outside the working class No
Positivity failure Some treatment levels are unavailable at relevant \(X\) No
Unmeasured confounding Important common causes are absent from \(X\) No

3 Convergence Rates of Nuisance Estimators

3.1 Rate Notation

We write

\[ \|\hat f-f_0\|_{2,P_X}=O_p(r_n) \]

when the stochastic error is of order \(r_n\). If

\[ r_n=n^{-\alpha}, \]

then larger \(\alpha\) means faster convergence. Throughout this lecture, a rate such as \(n^{-\alpha}\) refers to the root \(L_2\) error. Its square has rate \(n^{-2\alpha}\).

Always identify what quantity has the stated rate. The \(L_2\) norm, its square, parameter error, excess risk, and pointwise error generally have different rates.

3.2 Correctly Specified Parametric Models

For a fixed-dimensional, correctly specified regular parametric model \(f_\theta\),

\[ \sqrt n(\hat\theta-\theta_0)=O_p(1). \]

If the map \(\theta\mapsto f_\theta\) is sufficiently smooth, then

\[ \|f_{\hat\theta}-f_{\theta_0}\|_{2,P_X} = O_p(n^{-1/2}). \]

This conclusion requires correct specification and fixed effective dimension.

If a fixed parametric class is misspecified, an estimator may still satisfy

\[ \|\hat f-f^\dagger\|_{2,P_X}=O_p(n^{-1/2}), \]

while

\[ \|\hat f-f_0\|_{2,P_X} \xrightarrow{p} \|f^\dagger-f_0\|_{2,P_X}>0. \]

3.3 Nonparametric Rates and the Curse of Dimensionality

Suppose \(f_0\) has smoothness \(s>0\) over \(d\) continuous predictors. For a kernel or local-polynomial estimator of adequate order, under regular-design and moment conditions with \(h\to0\) and \(nh^d\to\infty\), a generic order bound has the form

\[ \|\hat f_h-f_0\|_{2,P_X} = \underbrace{O(h^s)}_{\text{smoothing bias}} + \underbrace{O_p\!\left\{(nh^d)^{-1/2}\right\}}_{\text{sampling variability}}. \]

Choose \(h\) so that neither component dominates. Equating their orders gives

\[ \begin{aligned} h^s &\asymp (nh^d)^{-1/2} \\ \Longleftrightarrow\quad h^{2s} &\asymp (nh^d)^{-1} \\ \Longleftrightarrow\quad nh^{2s+d} &\asymp 1 \\ \Longleftrightarrow\quad h^{2s+d} &\asymp n^{-1} \\ \Longleftrightarrow\quad h &\asymp n^{-1/(2s+d)}. \end{aligned} \]

At this bandwidth, both components have the same order,

\[ \underbrace{h^s}_{\text{bias}} \asymp \underbrace{(nh^d)^{-1/2}}_{\text{variability}} \asymp n^{-s/(2s+d)}. \]

Therefore,

\[ \|\hat f_h-f_0\|_{2,P_X} = O_p\!\left\{n^{-s/(2s+d)}\right\}. \]

Squaring the root-\(L_2\) rate gives \(n^{-2s/(2s+d)}\). Equivalently, under regular-design and moment conditions, the minimax squared \(L_2\) risk over a Hölder-type smoothness class behaves as

\[ \inf_{\tilde f_n} \sup_{f_0\in\mathcal H^s(L)} E\|\tilde f_n-f_0\|_{2,P_X}^2 \asymp n^{-2s/(2s+d)}. \]

This rate is a benchmark for a specified smoothness class, not a universal rate for random forests, boosting, neural networks, or every function called “nonparametric.”

How to read the plot. The vertical coordinate is the squared-\(L_2\) rate benchmark \(n^{-2s/(2s+d)}\) itself, evaluated at \(n=1,000\); the unknown multiplicative constants are omitted. Larger values indicate a worse squared-error benchmark. For fixed \(s\), increasing \(d\) raises the curve because \(2s/(2s+d)\) becomes smaller, so convergence is slower. The red \(n^{-1/2}\) line is the squared-error equivalent of the later root-\(L_2\) \(n^{-1/4}\) boundary. A curve crosses it at \(d=2s\); \(d<2s\) gives a faster rate, whereas equality is only the boundary and does not by itself establish the required little-\(o\) condition.

The complementary figure below fixes \(d=4\) and instead varies the sample size \(n\).

How to read the sample-size plot. Here the horizontal axis is \(n\), the dimension is fixed at \(d=4\), and lower curves indicate faster convergence. On the log–log scale, the slopes are \(-1/3\), \(-1/2\), and \(-2/3\) for \(s=1,2,4\), respectively. Thus, the \(s=1\) benchmark is slower than the squared-risk boundary \(n^{-1/2}\), the \(s=2\) benchmark lies exactly on it, and the \(s=4\) benchmark is faster. The plotted rates again omit unknown multiplicative constants.

3.4 Flexible Learning Does Not Eliminate Assumptions

A growing sieve or adaptive learner can reduce approximation error as \(n\) grows. Its actual rate still depends on structure such as:

  • smoothness or bounded variation;
  • sparsity;
  • low-dimensional or additive structure;
  • compositional structure;
  • regularization and tuning; and
  • the distribution of \(X\).

Calling a method “machine learning” does not prove consistency, a particular convergence rate, or causal validity.

No learner is uniformly best. Performance depends on how its structure and tuning match the data-generating process. A defensible library retains a simple parametric baseline alongside prespecified flexible learners; any stacking weights, hyperparameters, scaling, and feature selection must be learned inside the training folds. Repeating cross-fitting can reveal split instability, but it does not remove systematic bias. These practical lessons are emphasized by Ahrens et al. (2024).

4 Simulation: Parametric, Nonparametric, and Misspecified Rates

4.1 Simulation Design

We now estimate the same one-dimensional conditional mean

\[ f_0(x)=1+1.25x+1.5x^2 \]

from data generated by

\[ X_i\sim\operatorname{Uniform}(-1,1), \qquad Y_i=f_0(X_i)+\varepsilon_i, \qquad \varepsilon_i\sim N(0,0.7^2). \]

We compare three learners:

  1. Correct quadratic regression: a fixed, correctly specified parametric model;
  2. Histogram regression: a piecewise-constant nonparametric learner with \(K_n\asymp n^{1/3}\) bins; and
  3. Misspecified linear regression: a fixed model that omits the quadratic term.

The histogram is deliberately a first-order learner: its root-\(L_2\) bias is \(O(K_n^{-1})\) and its variability is \(O_p\{\sqrt{K_n/n}\}\). Balancing these terms gives \(K_n\asymp n^{1/3}\) and error \(O_p(n^{-1/3})\). Although the truth is twice differentiable, a piecewise-constant histogram does not exploit all of that smoothness; a suitable local-linear method could attain the earlier \(n^{-2/5}\) benchmark for \(s=2\) and \(d=1\).

For each fitted learner, error is evaluated against the true conditional function on a fixed dense midpoint grid that approximates integration under \(X\sim\operatorname{Uniform}(-1,1)\):

\[ \widehat{\operatorname{Err}} = \left[ \frac{1}{M} \sum_{j=1}^M \{\hat f(X_j^{\mathrm{eval}})-f_0(X_j^{\mathrm{eval}})\}^2 \right]^{1/2}. \]

This is function-estimation error. It does not include the irreducible noise \(\varepsilon\). The midpoint quadrature grid is fixed independently of every training sample.

truth_regression <- function(x) {
  1 + 1.25 * x + 1.50 * x^2
}

histogram_predict <- function(x, y, new_x, bins) {
  breaks <- seq(-1, 1, length.out = bins + 1L)
  training_bin <- findInterval(x, breaks, all.inside = TRUE)
  evaluation_bin <- findInterval(new_x, breaks, all.inside = TRUE)

  bin_mean <- tapply(y, training_bin, mean)
  prediction <- as.numeric(bin_mean[as.character(evaluation_bin)])
  prediction[!is.finite(prediction)] <- mean(y)
  prediction
}

simulate_rate_once <- function(n, evaluation_x) {
  x <- runif(n, -1, 1)
  y <- truth_regression(x) + rnorm(n, sd = 0.70)
  truth <- truth_regression(evaluation_x)

  correct_fit <- lm(y ~ x + I(x^2))
  misspecified_fit <- lm(y ~ x)

  correct_prediction <- predict(
    correct_fit,
    newdata = data.frame(x = evaluation_x)
  )
  misspecified_prediction <- predict(
    misspecified_fit,
    newdata = data.frame(x = evaluation_x)
  )
  histogram_prediction <- histogram_predict(
    x = x,
    y = y,
    new_x = evaluation_x,
    bins = max(2L, as.integer(round(n^(1 / 3))))
  )

  c(
    `Correct quadratic` = sqrt(
      mean((correct_prediction - truth)^2)
    ),
    `Histogram regression` = sqrt(
      mean((histogram_prediction - truth)^2)
    ),
    `Misspecified linear` = sqrt(
      mean((misspecified_prediction - truth)^2)
    )
  )
}

set.seed(20260719)
rate_sample_sizes <- c(200L, 400L, 800L, 1600L, 3200L, 6400L)
rate_evaluation_size <- 801L
rate_evaluation_x <- -1 +
  2 * (seq_len(rate_evaluation_size) - 0.5) / rate_evaluation_size
rate_method_names <- c(
  "Correct quadratic",
  "Histogram regression",
  "Misspecified linear"
)

rate_errors <- array(
  NA_real_,
  dim = c(
    length(rate_sample_sizes),
    rate_simulation_reps,
    length(rate_method_names)
  ),
  dimnames = list(
    n = rate_sample_sizes,
    replicate = seq_len(rate_simulation_reps),
    method = rate_method_names
  )
)

for (sample_index in seq_along(rate_sample_sizes)) {
  for (replicate_index in seq_len(rate_simulation_reps)) {
    rate_errors[sample_index, replicate_index, ] <-
      simulate_rate_once(
        n = rate_sample_sizes[sample_index],
        evaluation_x = rate_evaluation_x
      )
  }
}

if (any(!is.finite(rate_errors)) || any(rate_errors <= 0)) {
  stop("The rate simulation produced an invalid function error.")
}

rate_mean <- apply(rate_errors, c(1, 3), mean)
rate_lower <- apply(
  rate_errors,
  c(1, 3),
  quantile,
  probs = 0.10,
  names = FALSE
)
rate_upper <- apply(
  rate_errors,
  c(1, 3),
  quantile,
  probs = 0.90,
  names = FALSE
)

rate_slope <- vapply(
  rate_method_names,
  function(method) {
    unname(
      coef(
        lm(log(rate_mean[, method]) ~ log(rate_sample_sizes))
      )[2]
    )
  },
  numeric(1)
)

rate_slope_table <- data.frame(
  Method = rate_method_names,
  `Expected asymptotic slope` = c(-0.50, -1 / 3, 0),
  `Fitted simulation slope` = unname(round(rate_slope, 3)),
  check.names = FALSE,
  row.names = NULL
)

knitr::kable(
  rate_slope_table,
  align = c("l", "r", "r"),
  row.names = FALSE
)
Method Expected asymptotic slope Fitted simulation slope
Correct quadratic -0.5000000 -0.499
Histogram regression -0.3333333 -0.340
Misspecified linear 0.0000000 -0.004

The teaching render uses 120 independent repetitions at each sample size. The slopes are descriptive finite-sample summaries; they do not prove a minimax theorem.

4.2 One Realized Training Sample

The quadratic model is smooth and accurate because its working class contains the truth. The histogram is rougher but becomes more detailed as \(n\) increases. The linear fit is stable and simple, but its class cannot represent the curved truth.

4.3 Empirical Log-Log Convergence Plot

Interpretation. In this run, the fitted slopes are approximately -0.50 for the correct parametric learner, -0.34 for the nonparametric histogram, and -0.00 for the misspecified linear learner. The red curve approaches the theoretical misspecification plateau \(1/\sqrt{5}\approx 0.447\). The shaded regions show the pointwise 10th–90th percentiles of the repeated-sample function error at each displayed sample size.

A learner can have a very small training error or a fast parameter-estimation rate while remaining inconsistent for the true nuisance function. The relevant error for the later doubly robust argument is distance from the true nuisance function, not distance from the pseudo-truth.

5 From Nuisance Error to Double Robustness

5.1 Pathwise Differentiability, the EIF, and the Remainder

Let \(\mathcal M\) denote the observed-data model and, for each \(P\in\mathcal M\), define the distribution-indexed nuisance functions

\[ e_P(x)=P(T=1\mid X=x), \qquad m_{t,P}(x)=E_P(Y\mid T=t,X=x),\quad t\in\{0,1\}, \]

where the conditional probability and expectation are taken under \(P\).

After identification, write the risk difference as the observed-data functional

\[ \boxed{\Psi:\mathcal M\to\mathbb R,\qquad \Psi(P)=E_P\{m_{1,P}(X)-m_{0,P}(X)\}.} \]

Under positivity, \(m_{t,P}\) is defined on the relevant \(P_X\) support and \(\Psi(P)\) is a version-independent standardized observed-data contrast. Consistency and conditional exchangeability then give it the causal interpretation

\[ RD=E\{Y(1)-Y(0)\}=\Psi(P). \]

If consistency or exchangeability fails while positivity holds, \(\Psi(P)\) remains an associational standardized contrast but need not equal the causal risk difference. If positivity fails, the integral generally depends on how \(m_{t,P}\) is extrapolated into unsupported covariate regions and is not nonparametrically identified without an additional convention or model.

Let \(\{P_\varepsilon:\varepsilon\in(-\delta,\delta)\}\) be a regular one-dimensional submodel through \(P=P_0\), with score \(S(O)=\left.\partial_\varepsilon\log p_\varepsilon(O)\right|_{\varepsilon=0}\) satisfying \(E_P\{S(O)\}=0\) and \(E_P\{S(O)^2\}<\infty\). Pathwise differentiability means that there is a mean-zero, square-integrable function \(D_{\Psi}^*(\,\cdot\,;P)\) such that, for every such submodel,

\[ \left.\frac{d}{d\varepsilon}\Psi(P_\varepsilon)\right|_{\varepsilon=0} =E_P\{D_{\Psi}^*(O;P)S(O)\}. \]

For the risk-difference functional, suppressing the \(P\) subscripts on \(e_P\) and \(m_{t,P}\) on the right-hand side,

\[ D_{\Psi}^*(O;P)=m_1(X)-m_0(X)+\frac{T}{e(X)}\{Y-m_1(X)\}-\frac{1-T}{1-e(X)}\{Y-m_0(X)\}-\Psi(P). \]

In the nonparametric observed-data model, whenever the displayed gradient belongs to \(L_2(P)\), it is the efficient influence function (EIF), and \(E_P\{D_{\Psi}^*(O;P)^2\}\) is the semiparametric efficiency bound. Strong overlap together with \(E_P(Y^2)<\infty\) is a transparent sufficient condition; weak positivity alone is not. For any \(Q\) for which the terms are well-defined, define the von Mises residual algebraically by

\[ \Psi(Q)-\Psi(P)=\int D_{\Psi}^*(o;P)\,d(Q-P)(o)+R_2(Q,P), \]

Along a regular submodel \(Q=P_\varepsilon\), pathwise differentiability implies \(R_2(P_\varepsilon,P)=o(\varepsilon)\). Uniform control over more general nearby \(Q\) requires additional smoothness and topology conditions. The exact product drift derived below is generated by the risk-difference DR score; not every semiparametric target has the same product form.

An efficient influence function is the canonical gradient that represents a parameter’s pathwise derivative; it is not itself a derivative operator. It becomes the first-order influence function of a particular estimator only after that estimator’s empirical-process and remainder terms are controlled. Merely computing a numerical Gateaux derivative or writing down an orthogonal score does not by itself establish consistency or efficiency. Here pathwise differentiability with a finite efficiency bound requires \(D_{\Psi}^*(O;P)\in L_2(P)\); strong overlap and finite outcome second moments are sufficient. Nonsmooth or infinite-dimensional targets may require different theory. Kennedy’s review (Chapter 10, pp. 213-230) develops these distinctions; Mai (2019, pp. 59-60) gives concrete warnings about numerical influence-function construction without the required justification.

5.2 The DR (AIPW) Score

The course notation for the true nuisance functions is \(\eta=(e,m_0,m_1)\). To distinguish truth from an arbitrary working choice, write

\[ \bar\eta=(\bar e,\bar m_0,\bar m_1). \]

The augmented inverse-probability-weighted (AIPW) form of Lecture 4’s doubly robust (DR) score for the risk difference is

\[ \phi_{\Psi}^{DR}(O;\bar\eta)=\bar m_1(X)-\bar m_0(X)+\frac{T}{\bar e(X)}\{Y-\bar m_1(X)\}-\frac{1-T}{1-\bar e(X)}\{Y-\bar m_0(X)\}. \]

This \(\phi_{\Psi}^{DR}\) is an uncentered estimating score, not the EIF itself. At the true nuisance functions,

\[ E_P\{\phi_{\Psi}^{DR}(O;\eta)\}=\Psi(P)=RD, \qquad D_{\Psi}^*(O;P)=\phi_{\Psi}^{DR}(O;\eta)-\Psi(P). \]

To state Neyman orthogonality precisely, define the population moment \(M(\theta,\bar\eta)=E_P\{\phi_{\Psi}^{DR}(O;\bar\eta)-\theta\}\). At \(\theta=\Psi(P)\), for every regular nuisance path \(\eta_r\) through the truth \(\eta\),

\[ \left.\frac{d}{dr}M\{\Psi(P),\eta_r\}\right|_{r=0}=0. \]

Thus, the population moment has no first-order drift along nuisance directions at the truth. Orthogonality does not make an inconsistent nuisance learner consistent; it makes the leading local drift second order.

A DR estimator averages this score after replacing the nuisance functions by fitted values. With cross-fitted nuisances, observation \(i\) is evaluated using functions trained without observation \(i\). Let \(k(i)\) denote the fold containing observation \(i\). We use \(\widehat\Psi_{\mathrm{CF}}^{DR}\) for the statistical-functional notation and retain \(\widehat{RD}_{\mathrm{CF}}^{DR}\) as the equivalent course label:

\[ \widehat\Psi_{\mathrm{CF}}^{DR} \equiv \widehat{RD}_{\mathrm{CF}}^{DR} =\frac1n\sum_{i=1}^n\phi_{\Psi}^{DR}\{O_i;\hat\eta^{(-k(i))}\}. \]

5.3 Exact Population Drift

The central identity is exact, not merely a heuristic. Consider candidate functions \((\bar e,\bar m_0,\bar m_1)\) such that \(0<\bar e(X)<1\) \(P_X\)-almost surely and the displayed score is integrable. Condition on \(X\). Since

\[ E(TY\mid X)=e(X)m_1(X) \]

and

\[ E\{(1-T)Y\mid X\}=\{1-e(X)\}m_0(X), \]

we obtain

\[ E\{\phi_{\Psi}^{DR}(O;\bar\eta)\mid X\}=\bar m_1(X)-\bar m_0(X)+\frac{e(X)}{\bar e(X)}\{m_1(X)-\bar m_1(X)\}-\frac{1-e(X)}{1-\bar e(X)}\{m_0(X)-\bar m_0(X)\}. \]

Subtracting \(m_1(X)-m_0(X)\) and simplifying gives

\[ \boxed{\mathcal B_P(\bar\eta)\equiv E_P\{\phi_{\Psi}^{DR}(O;\bar\eta)\}-\Psi(P)=E_P\!\left[(\bar e-e)\left\{\frac{\bar m_1-m_1}{\bar e}+\frac{\bar m_0-m_0}{1-\bar e}\right\}\right].} \]

The quantity \(\mathcal B_P(\bar\eta)\) is the population drift of the nuisance-indexed score. Its product form makes it a second-order nuisance remainder in the estimator expansion. The distinct notation is intentional: \(R_2(Q,P)\) above is indexed by two full distributions, whereas \(\mathcal B_P(\bar\eta)\) is indexed by candidate nuisance functions while the evaluation law \(P\) is fixed.

5.4 Model Double Robustness

The exact drift is zero if either:

  • \(\bar e=e\), even when the outcome regressions are misspecified; or
  • \(\bar m_1=m_1\) and \(\bar m_0=m_0\), even when the propensity score is misspecified.

For a risk difference, saying “the outcome model is correct” means that both treatment-specific outcome regressions are correct.

Zero population drift is the algebraic core of model double robustness. To turn it into estimator consistency, also require fitted overlap and integrability, a law of large numbers for the evaluated score, and convergence of the fitted nuisance functions to stable limits. Under those regularity conditions:

Propensity nuisance Outcome nuisances Outcome regression IPW DR (AIPW)
Correct Correct Consistent Consistent Consistent
Correct Misspecified Generally biased Consistent Consistent
Misspecified Correct Consistent Generally biased Consistent
Misspecified Misspecified Generally biased Generally biased Generally biased in general

This is model double robustness for consistency. It does not say that an estimator is always unbiased in a finite sample, that it is automatically efficient, or that a standard influence-function confidence interval remains valid under every one-correct-model configuration.

5.5 Positivity and the Product Remainder

Three conditions should not be conflated:

  • Causal positivity: \(0<e(X)<1\) almost surely is part of identification.
  • Strong overlap: \(\delta\le e(X)\le1-\delta\) almost surely is a transparent sufficient condition for a square-integrable EIF when \(Y\) is binary.
  • Fitted overlap: \(\epsilon\le\hat e^{(-k)}(X)\le1-\epsilon\) with probability approaching one, uniformly over folds, controls the fitted denominators. In practice this is an assumption about the fitted rule, not proof of population overlap.

More formally, define the fitted-overlap event

\[ \mathcal A_n = \left\{ \epsilon\le \hat e^{(-k)}(x)\le1-\epsilon \text{ for every }k\text{ and }P_X\text{-almost every }x \right\}, \qquad P(\mathcal A_n)\longrightarrow1. \]

Let \(I_k\) be test fold \(k\), let \(\mathcal D_{-k}\) contain the observations used to train its nuisance functions, and let \(E_P\) integrate over a fresh \(O\sim P\) while holding \(\mathcal D_{-k}\) fixed. Define the fold-specific drift and fold-uniform errors by

\[ \mathcal B_k =E_P\!\left[\phi_{\Psi}^{DR}\{O;\hat\eta^{(-k)}\}-\phi_{\Psi}^{DR}(O;\eta)\mid\mathcal D_{-k}\right], \]

\[ r_{e,n}=\max_{1\le k\le K}\|\hat e^{(-k)}-e\|_{2,P_X},\qquad r_{m,n}=\max_{1\le k\le K}\sum_{t=0}^1\|\hat m_t^{(-k)}-m_t\|_{2,P_X}. \]

By the exact drift identity and Cauchy–Schwarz, on \(\mathcal A_n\) the fitted-overlap condition gives the one-line bound

\[ \boxed{\max_{1\le k\le K}|\mathcal B_k|\le\epsilon^{-1}r_{e,n}r_{m,n}.} \]

All nuisance functions in the compact drift formulas are evaluated at \(X\).

5.6 Exact Cross-Fitted Expansion

Define the centered foldwise empirical-process term and the fold-averaged population drift by

\[ \mathcal E_n=\frac1n\sum_{k=1}^K\sum_{i\in I_k}\!\left[\phi_{\Psi}^{DR}\{O_i;\hat\eta^{(-k)}\}-\phi_{\Psi}^{DR}(O_i;\eta)-E_P\!\left\{\phi_{\Psi}^{DR}(O;\hat\eta^{(-k)})-\phi_{\Psi}^{DR}(O;\eta)\mid\mathcal D_{-k}\right\}\right], \]

\[ \bar{\mathcal B}_n=\sum_{k=1}^K\frac{|I_k|}{n}\mathcal B_k. \]

The estimator therefore has the exact decomposition

\[ \widehat\Psi_{\mathrm{CF}}^{DR}-\Psi(P)=\frac1n\sum_{i=1}^nD_{\Psi}^*(O_i;P)+\mathcal E_n+\bar{\mathcal B}_n. \]

The three terms are, respectively, the first-order EIF average, the foldwise empirical-process error, and the second-order nuisance drift.

Assume the observations are iid, fold assignment is independent of the data, \(K\) is fixed and finite, and \(\min_k|I_k|/n\) is bounded away from zero. Every data-dependent operation—including preprocessing, feature selection, tuning, and stopping—must use only \(\mathcal D_{-k}\). If

\[ \max_{1\le k\le K}\|\phi_{\Psi}^{DR}(\,\cdot\,;\hat\eta^{(-k)})-\phi_{\Psi}^{DR}(\,\cdot\,;\eta)\|_{2,P}=o_p(1), \]

then conditional mean-zero and variance arguments give \(\mathcal E_n=o_p(n^{-1/2})\). Cross-fitting removes the need for a Donsker-class restriction in this fixed-\(K\) argument, but it does not remove score convergence, overlap, moment, or remainder conditions. Without sample splitting, a standard alternative is to impose suitable Donsker or entropy conditions and stochastic equicontinuity. Leave-one-out or growing-\(K\) schemes require separate analysis.

The product bound gives \(\bar{\mathcal B}_n=O_p(r_{e,n}r_{m,n})\). Hence a useful sufficient condition for the nuisance drift to be negligible on the root-\(n\) scale is

\[ \boxed{r_{e,n}r_{m,n}=o_p(n^{-1/2}).} \]

Together with square integrability of the EIF and the preceding conditions,

\[ \sqrt n\left\{\widehat\Psi_{\mathrm{CF}}^{DR}-\Psi(P)\right\}=\frac1{\sqrt n}\sum_{i=1}^nD_{\Psi}^*(O_i;P)+o_p(1)\rightsquigarrow N(0,\sigma^2),\qquad \sigma^2=E_P\{D_{\Psi}^*(O;P)^2\}. \]

For iid data, define \(\hat D_i=\phi_{\Psi}^{DR}\{O_i;\hat\eta^{(-k(i))}\}-\widehat\Psi_{\mathrm{CF}}^{DR}\) and \(\hat\sigma^2=n^{-1}\sum_{i=1}^n\hat D_i^2\). The canonical Wald interval is

\[ \widehat\Psi_{\mathrm{CF}}^{DR}\pm1.96\sqrt{\hat\sigma^2/n}. \]

This variance formula must be modified for clustering, repeated observations, survey sampling, or other dependence.

5.7 The \(n^{-1/4}\) Balanced Benchmark

Suppose

\[ r_{e,n}=O_p(n^{-\alpha_e}), \qquad r_{m,n}=O_p(n^{-\alpha_m}). \]

Then

\[ r_{e,n}r_{m,n} = O_p\{n^{-(\alpha_e+\alpha_m)}\}. \]

A useful sufficient exponent condition is

\[ \boxed{ \alpha_e+\alpha_m>\frac12. } \]

If the two rates are balanced, the familiar sufficient condition is

\[ r_{e,n}=o_p(n^{-1/4}), \qquad r_{m,n}=o_p(n^{-1/4}). \]

Exactly \(O_p(n^{-1/4})\) for both nuisances gives only an \(O_p(n^{-1/2})\) product. At that boundary, the scaled drift need not disappear.

Ordered-pair interpretation. A point \((\alpha_e,\alpha_m)=(0.30,0.30)\) means that the propensity error is of order \(n^{-0.30}\) and the combined outcome-regression error is of order \(n^{-0.30}\). Their product has order \(n^{-0.60}\), which is smaller than \(n^{-1/2}\). The rates need not be equal: \((0.50,0.10)\) also has exponent sum \(0.60\).

5.8 Connecting the Nonparametric Rate to the DR Benchmark

If both nuisance sides have the same generic \(s\)-smooth, \(d\)-dimensional rate,

\[ \alpha_e=\alpha_m=\frac{s}{2s+d}, \]

then the strict product-rate condition becomes

\[ \frac{2s}{2s+d}>\frac12 \quad\Longleftrightarrow\quad 2s>d. \]

Thus, the generic smoothness rate is faster than the balanced \(n^{-1/4}\) benchmark when \(s>d/2\). At \(s=d/2\), it lies exactly on the boundary.

Two different ideas must be separated.

  • A learner-specific optimal rate describes how quickly a nuisance function can be estimated over a specified function class under a specified loss. For example, the generic root-\(L_2\) benchmark for an \(s\)-smooth function of \(d\) continuous predictors is \(n^{-s/(2s+d)}\). A different class, dimension, smoothness assumption, or loss produces a different rate.
  • The DR inferential threshold comes from the remainder, not from minimax prediction theory. Root-\(n\) inference requires \[ r_{e,n}r_{m,n}=o_p(n^{-1/2}). \] If the two errors have the same order, \(r_{e,n}\asymp r_{m,n}\asymp r_n\), this reduces to \(r_n^2=o_p(n^{-1/2})\), or \[ r_n=o_p(n^{-1/4}). \]
  • The equality case is only a boundary. If both errors are merely \(O_p(n^{-1/4})\), their product is generally only \(O_p(n^{-1/2})\), not the required little-\(o_p(n^{-1/2})\). Additional structure would be needed to make the scaled remainder vanish.
  • Equal rates are not required. If \(r_{e,n}=O_p(n^{-\alpha_e})\) and \(r_{m,n}=O_p(n^{-\alpha_m})\), the relevant condition is \(\alpha_e+\alpha_m>1/2\). One nuisance may converge more slowly than \(n^{-1/4}\) when the other converges sufficiently faster.

Therefore, \(n^{-1/4}\) is a balanced sufficient-rate benchmark for this DR remainder, not a universal optimum and not a claim that every nuisance learner should attain exactly that rate.

5.9 Model DR, Rate DR, and Local Efficiency Are Different

Claim Sufficient conditions in this lecture What follows
Neyman orthogonality The pathwise derivative of \(M\{\Psi(P),\eta_r\}\) is zero at the truth First-order local nuisance perturbations cancel in the population moment
Model double robustness \(\bar e=e\), or both \(\bar m_1=m_1\) and \(\bar m_0=m_0\) The population drift is zero, so consistency can hold
Rate double robustness Fixed finite \(K\), foldwise score convergence, overlap and moment bounds, and \(r_{e,n}r_{m,n}=o_p(n^{-1/2})\) Root-\(n\) asymptotic linearity with the canonical first-order term
Local efficiency At the intersection model, the estimator is asymptotically linear with influence function \(D_{\Psi}^*\) The point estimator attains the semiparametric efficiency bound
Valid canonical Wald inference Local efficiency plus consistent estimation of \(E_P\{D_{\Psi}^*(O;P)^2\}\) The canonical influence-function confidence interval has asymptotically correct coverage

These are not interchangeable. If one nuisance is estimated at \(O_p(n^{-1/2})\) while the other is consistent at \(o_p(1)\), their product is \(o_p(n^{-1/2})\). If the second nuisance instead has a fixed nonzero misspecification error, the product is generally only \(O_p(n^{-1/2})\). The estimator may remain consistent by model double robustness, but the canonical efficient-influence-function approximation need not be the complete first-order expansion.

\[ O_p(n^{-1/2})\times O_p(1)=O_p(n^{-1/2}),\qquad\text{not necessarily }o_p(n^{-1/2}). \]

Under a union model with one misspecified parametric nuisance, a valid influence function can contain extra terms for estimating the nuisance probability limits; under the intersection model where both nuisances are correct, those terms simplify and local efficiency can be recovered. This distinction is explicit in Vermeulen (2016, pp. 9-17). Gao, Zhang, and Yang (2025, pp. 6-7) use the term rate double robustness for analogous sums of products of nuisance errors.

Sensitivity robustness is different again. It asks whether conclusions persist as an unidentified causal assumption is varied through scientifically plausible sensitivity parameters. Flexible nuisance learning, orthogonality, and cross-fitting do not identify such parameters or repair failed exchangeability. Gao, Zhang, and Yang (2025) carefully separate nuisance-rate robustness from sensitivity analysis in an externally controlled-trial setting.

6 Cross-Fitting

6.1 Honest Nuisance Prediction for Every Observation

With \(K\) folds, cross-fitting repeats two operations:

  1. fit the nuisance functions using observations outside fold \(k\); and
  2. evaluate the DR (AIPW) score on observations inside fold \(k\).

Cross-fitting makes the nuisance prediction for each observation out-of-fold. It helps control the empirical-process or own-observation term that can arise with adaptive learners.

Cross-fitting does not make a misspecified nuisance function correct, improve poor overlap, increase a learner’s intrinsic convergence rate, or repair violations of causal identification. Hyperparameter tuning and preprocessing must also occur inside each training fold.

7 Practical Workflow

7.1 What to Examine

A careful nuisance-model workflow includes:

  1. Prespecify the target and score. Decide which causal parameter and orthogonal score are being estimated.
  2. Use honest prediction. Tune and preprocess inside training folds; generate out-of-fold nuisance predictions.
  3. Check propensity behavior. Examine overlap, extreme fitted propensities, weight tails, and effective sample size.
  4. Check outcome prediction. Examine held-out loss and calibration overall and within treatment groups.
  5. Use multiple reasonable learners. Compare a transparent parametric baseline with prespecified flexible alternatives.
  6. Assess target stability. Compare estimates across learners, fold assignments, and defensible truncation choices.
  7. Report uncertainty honestly. Match the variance calculation to cross-fitting, clustering, survey design, and the estimator actually used.

7.2 What Prediction Metrics Cannot Prove

Good held-out prediction is useful but not equivalent to the nuisance conditions needed for a causal target.

  • Overall outcome prediction may be dominated by regions with abundant data rather than the overlap-sensitive regions that matter for treatment contrasts.
  • A high propensity-score AUC can indicate strong treatment separation and therefore poor overlap; it is not automatically desirable.
  • Calibration to observed outcomes does not reveal counterfactual calibration.
  • Similar causal estimates across learners are reassuring sensitivity evidence, not proof of identification.

7.3 Software Used in This Lecture

Task R package and function
Correct and misspecified parametric regression stats::lm() and stats::glm()
Histogram regression Base R findInterval() and tapply()
Cross-fitting Base R sample splitting plus stats::predict()
Tables knitr::kable()
Figures Base R graphics

No external machine-learning package is required. The purpose is to isolate the statistical ideas rather than compare software ecosystems.

8 Key Takeaways

  1. A nuisance function supports estimation of the target parameter but is not itself the final estimand.
  2. Total nuisance error contains both finite-sample estimation error and approximation error from the working class.
  3. A correctly specified fixed-dimensional parametric nuisance can achieve a root-\(n\) function rate; a misspecified fixed class can converge at a root-\(n\) rate to the wrong pseudo-truth.
  4. A generic \(s\)-smooth, \(d\)-dimensional nonparametric root-\(L_2\) rate is \(n^{-s/(2s+d)}\), illustrating the curse of dimensionality.
  5. The exact DR (AIPW) drift is a product of propensity and outcome nuisance errors.
  6. Model double robustness gives consistency if either the propensity nuisance or both outcome nuisances are consistently estimated.
  7. Root-\(n\) efficient-influence-function inference requires the stronger product-rate condition \(r_{e,n}r_{m,n}=o_p(n^{-1/2})\), plus positivity and empirical-process control.
  8. The symmetric \(n^{-1/4}\) rate is a critical benchmark, not a universal optimal nuisance rate.
  9. Cross-fitting controls own-observation dependence; it does not repair misspecification or causal-identification failures.

9 Technical References Used

  • Kennedy, Edward H. Semiparametric Doubly Robust Targeted Double Machine Learning: A Review. Chapter 10, especially pp. 207–236. DOI: 10.1201/9781003216223-10.
  • Ahrens, Achim, Victor Chernozhukov, Christian Hansen, Damian Kozbur, Mark Schaffer, and Thomas Wiemann. Robust Causal Inference Using Double/Debiased Machine Learning: A Guide for Empirical Research. 2024 guide and lecture materials.
  • Gao, Chenyin, Xiang Zhang, and Shu Yang. Doubly Robust Omnibus Sensitivity Analysis of Externally Controlled Trials with Intercurrent Events. Biometrics 81(2), 2025, ujaf047. DOI: 10.1093/biomtc/ujaf047.
  • Mai, Jiajun. Semiparametric Estimation in Targeted Learning by Numerical Construction of the Efficient Influence Function. MSc thesis, McGill University, 2019.
  • Vermeulen, Karel. A Double Robust Extension of the Mann–Whitney Test to Adjust for Confounding in Observational Studies. Master dissertation, Ghent University, 2016.

These sources play different roles. Kennedy supplies the general pathwise-differentiability, one-step, and cross-fitting theory; Ahrens et al. motivate honest learner comparison and practical DML; Gao, Zhang, and Yang illustrate rate double robustness and its separation from sensitivity analysis; Mai clarifies the requirements for a valid EIF; and Vermeulen distinguishes union-model consistency from intersection-model efficiency. The last paper studies a pairwise Mann–Whitney estimand and is not used as the proof of the risk-difference remainder.

10 Reproducibility Notes

The pseudo-truth Monte Carlo study and the learner-rate simulation each use 120 repetitions per sample size, controlled by the environment variable NUISANCE_RATE_SIM_REPS. Fixed seeds make the displayed teaching results reproducible.

For a fast smoke test, set NUISANCE_RATE_SIM_REPS to a small value such as 5. For smoother Monte Carlo summaries, increase it after checking runtime and Monte Carlo stability. The default value is a teaching default, not a universal research recommendation.

This lecture provides the nuisance-rate and misspecification foundation for the Doubly Robust Estimation lecture. The next application is to carry these conditions into the smoking–CVD example, compare nuisance learners, and examine how the DR correction and its influence-function uncertainty behave in that substantive setting.