Residual covariance and nuisance-estimation error

Population nuisance bias, held-out evaluation, and conditions for first-order inference

Before reading: Conditional expectation and L2 geometry, Model restrictions and efficiency, and From an influence function to an estimator. This application reuses the capstone’s nuisance-bias and held-out empirical-term arguments for a covariance target.

1 The problem this page solves

The Model restrictions and efficiency page introduced an asymptotic efficiency benchmark and showed how its influence function can depend on unknown nuisance quantities. Here the nuisance is a regression function: how a gene’s mean, or variance, depends on covariates such as cell composition, age, batch, or an interferon score. The estimand here is the average conditional covariance between two genes after removing their conditional means given Z, rather than a causal or conditional-independence claim. We do not want to assume the regressions are linear, so we estimate them flexibly, and then we need to know how the error in those flexible estimates propagates into the estimand. For this target, an orthogonal residual-product equation removes the first-order population nuisance error. Its exact remaining bias and its empirical fluctuation still need control.

This transfer application reuses Neyman orthogonality and cross-fitting for a different target. The task is to determine which derivatives and remainder terms change for covariance.

2 The conditional covariance, two ways

Let Z be a covariate (say, a cell-composition score), and let m_X(z) = E[X \mid Z = z] and m_Y(z) = E[Y \mid Z = z]. The estimand is the average conditional covariance

\psi = E\big[\operatorname{Cov}(X, Y \mid Z)\big].

This is the covariance between two genes after removing what each shares with Z, averaged over the population. There are at least two correct formulas for it:

\psi = E\big[(X-m_X(Z))(Y-m_Y(Z))\big]

\psi = E[XY]-E[m_X(Z)m_Y(Z)]

The first is the residual form and the second the difference form. Both are exactly right if m_X and m_Y are the true regression functions. Both can be turned into estimators by plugging in estimates \hat m_X and \hat m_Y and averaging. They behave completely differently when the estimates are imperfect.

To verify their equality, expand the residual product. The tower property gives E[Xm_Y(Z)]=E[m_X(Z)m_Y(Z)], and the same identity holds with X and Y reversed. Two subtractions and one addition of that mean-product term leave the difference form. Algebraic equivalence at the truth does not imply equal sensitivity to estimated inputs.

Write W=(X,Y,Z) for one person’s record. Throughout, P f=E[f(W)] denotes a population average, P_n f=n^{-1}\sum_i f(W_i) a sample average, and \|f\|_2=\{E[f(W)^2]\}^{1/2} the population root-mean-square size. Existence of the covariance needs second moments of X and Y. A CLT for their residual product needs the product itself to have finite variance; suitable fourth moments are a convenient stronger condition.

3 Neyman orthogonality

We want a formula whose population value changes little when a fitted conditional mean is slightly wrong. Write the trained errors as \delta_X=\hat m_X-m_X and \delta_Y=\hat m_Y-m_Y. A fitted residual is U_X-\delta_X, so the exact expansion is

\begin{aligned} E[(U_X-\delta_X)(U_Y-\delta_Y)]&=\psi-E[U_X\delta_Y]-E[U_Y\delta_X] +E[\delta_X\delta_Y]\\ &=\psi+E[\delta_X\delta_Y]. \end{aligned}

The two single-error terms vanish because each error is a function of Z and each true residual has conditional mean zero given Z. Only their product remains. In contrast, plugging into the difference form leaves -E[m_X\delta_Y+m_Y\delta_X+\delta_X\delta_Y], including single-error terms. This exact comparison motivates the derivative definition that follows.

Hold the population law P fixed and write the population estimating map as M(m), with nuisance arguments m and true pair m_0=(m_X,m_Y) satisfying M(m_0)=\psi. The target is fixed; we are varying the nuisance arguments supplied to the formula. Perturb them along an admissible direction h. The equation is Neyman orthogonal when

\frac{\partial}{\partial\varepsilon}M(m_0+\varepsilon h) \Big|_{\varepsilon=0}=0.

The derivative is taken at the true nuisance for the permitted directions. A quadratic remainder additionally needs smoothness and norm controls; a zero derivative alone is not a universal rate theorem. For the residual product below, the second-order term follows from exact algebra.

Residual form. Write U_X=X-m_X(Z) and U_Y=Y-m_Y(Z) for the true conditional-mean residuals. Perturb m_X by \varepsilon h_X(Z):

\begin{aligned} \frac{\partial}{\partial\varepsilon}E[(U_X-\varepsilon h_X)U_Y]\Big|_0&=-E[h_X(Z)U_Y]\\ &=-E[h_X(Z)E(U_Y\mid Z)]\\ &=0, \end{aligned}

because E[Y - m_Y(Z) \mid Z] = 0 by definition of m_Y. The same holds for m_Y. The residual form is orthogonal. Its second-order term is E[h_X(Z) h_Y(Z)]\,\varepsilon^2: the product of the two nuisance errors.

Difference form. Perturb m_X:

\frac{\partial}{\partial\varepsilon} \{E[XY]-E[(m_X+\varepsilon h_X)m_Y]\}\Big|_0=-E[h_X(Z)m_Y(Z)],

This derivative is generally nonzero, so the difference-form plug-in is not generally orthogonal. Particular directions can still give zero derivative; orthogonality requires the condition for every relevant direction.

Here it is numerically. Inject a known error into the nuisance (inflate the true regression function by a factor 1 + \varepsilon) and watch the two formulas.

Code

set.seed(1)
n <- 20000
z <- runif(n); e1 <- rnorm(n); e2 <- rnorm(n)
x <- z + e1                 # m_X(z) = z
y <- z + 0.5 * e1 + e2      # m_Y(z) = z;  Cov(X, Y | Z) = 0.5 exactly
m_x <- z; m_y <- z          # the true nuisances

Residual   <- function(mx, my) mean((x - mx) * (y - my))
Difference <- function(mx, my) mean(x * y) - mean(mx * my)

eps <- c(0, 0.05, 0.1, 0.2, 0.4)
out <- sapply(eps, function(e) c(residual_form = Residual((1 + e) * m_x, (1 + e) * m_y),
                                 difference_form = Difference((1 + e) * m_x, (1 + e) * m_y)))
colnames(out) <- paste0("eps=", eps)
round(out, 3)

                eps=0 eps=0.05 eps=0.1 eps=0.2 eps=0.4
residual_form   0.494    0.495   0.498   0.508   0.548
difference_form 0.493    0.459   0.423   0.346   0.173

With m_X = m_Y = Z and Z uniform, the residual form’s error is \varepsilon^2 E[Z^2] = \varepsilon^2/3 (about 0.003 at \varepsilon = 0.1) and the difference form’s is -2\varepsilon/3 - \varepsilon^2/3 (about -0.07 at \varepsilon = 0.1). Those are population biases for this perturbation. The displayed sample averages also contain Monte Carlo fluctuation; these numerical biases are specific to the perturbation.

Code

grid <- seq(-0.3, 0.3, length.out = 201)
matplot(grid, cbind(residual = grid^2/3, difference = -(2*grid + grid^2)/3),
        type = "l", lty = c(1, 2), col = c("#2F80C8", "#D55E00"),
        xlab = "Nuisance perturbation epsilon", ylab = "Population bias")
abline(h = 0, v = 0, col = "gray70")
legend("topright", c("Residual form", "Difference form"), lty = 1:2,
       col = c("#2F80C8", "#D55E00"), bty = "n")

Bias curves show a flat tangent at zero for the residual form and a downward linear component for the difference form. — Figure 1: Exact population bias in the stated synthetic perturbation. The residual formula changes quadratically; the difference formula has a nonzero first-order slope. These curves are not sampled estimates.

The remainder of the residual form is a product of the two nuisance errors, so it also vanishes when either regression function is exactly right. Keep m_X exact and make m_Y badly wrong:

Code

m_y_wrong <- 3 * z + 1                                   # a badly wrong nuisance function
round(c(residual_form   = Residual(m_x, m_y_wrong),      # still near 0.5 (Monte Carlo SE about 0.02)
        difference_form = Difference(m_x, m_y_wrong)),   # far from 0.5
      3)

  residual_form difference_form 
          0.488          -0.675

The important condition is on the product of nuisance errors, not on the name of the learning algorithm. For root-n inference centered at the target, a usual sufficient condition here is

\|\delta_X\|_2\,\|\delta_Y\|_2=o_p(n^{-1/2}),

along with moment, consistency, and sampling conditions. Errors of exactly order n^{-1/4} each give a product of order n^{-1/2}, which need not be negligible; both being o_p(n^{-1/4}) is sufficient, and faster convergence of one can compensate for slower convergence of the other. Orthogonality removes the first-order population term, while cross-fitting helps control the empirical term. Neither promises valid inference for arbitrary nuisance fits. Chernozhukov et al. develop this distinction.

4 The connection to the efficient influence function

Efficient influence functions provide a systematic route to orthogonal equations. Orthogonality to nuisance scores is a model-geometric statement; Neyman orthogonality concerns derivatives of the chosen estimating equation in its nuisance arguments. For a valid influence-function construction these properties are connected under regularity conditions. We can verify the nuisance derivative directly for this residual-product equation. For \psi above, the EIF is

\operatorname{EIF}(X,Y,Z)=U_XU_Y-\psi,

and setting its sample average to zero gives the residual form. The naive difference-form plug-in has first-order nuisance sensitivity. It is not automatically an estimator with a different valid influence function: its behavior also depends on how the regressions are estimated. So the recipe for an orthogonal estimator is: find the EIF, and estimate by solving its empirical mean equals zero. The one-step estimator adds an estimated influence-function correction to an initial estimate. Writing P_n for a sample average,

\hat\psi_{\rm step}=\hat\psi_{\rm initial} +P_n\phi(W;\hat\psi_{\rm initial},\hat m).

Here the centering parameter is included in the estimated influence function. For this residual-covariance equation the initial centering cancels algebraically, leaving the residual-product average. Equivalence to an estimating-equation solution must be checked for the particular target and construction, including how the initial target estimate enters the estimated EIF.

Why is this the influence function of a target that itself contains conditional means? Along a smooth path with score s, differentiating E_t[U_{X,t}U_{Y,t}] gives a law-change term E[U_XU_Ys] and two mean-change terms. Each mean-change term vanishes by conditional residual centering, just as in the calculation above. Since E[s]=0, the derivative is E[(U_XU_Y-\psi)s]. In the unrestricted model this identifies the canonical gradient. It is a derivative of the whole functional, including its changing conditional means, rather than a calculation that assumes those means are known.

Deriving EIFs is a skill; Hines et al. (2022) teach it with the Gateaux-derivative method from the Influence functions page, and Kennedy (2022) catalogs common examples. Residual correlation requires a further normalization step below.

5 Cross-fitting

Cross-fitting separates a held-out observation from its own training fit. It does not make nuisance predictions correct or all held-out contributions jointly independent: different folds share training observations.

One more problem. If \hat m_X is fit on the same observations it is then evaluated on, the fit absorbs some of each observation’s own noise. The residual X_i - \hat m_X(Z_i) is then too small, and because Y_i - \hat m_Y(Z_i) is shrunk in a correlated way, the residual product is biased. For a fixed common linear smoother S, an error-covariance contribution is scaled by \operatorname{tr}\{(I-S)^\top(I-S)\}/N; a projection fit loses its rank in this trace. A smoothing spline is not generally a projection, and adaptive fits need their own analysis. The direction and size of the total finite-sample bias also depend on approximation error.

Cross-fitting prevents a held-out outcome from being used to train its own prediction. Split the sample into K folds. For each fold, fit the nuisances on the other K - 1 folds and evaluate the residuals on the held-out fold. Average the residual products over all folds. Each observation’s residual now uses a nuisance fit that never saw it.

Code

set.seed(2)
# Nonlinear nuisance so that a smoother is needed. Cov(X, Y | Z) = 0.5 again.
OneRep <- function(n, df = 15, K = 5) {
  z <- runif(n); e1 <- rnorm(n); e2 <- rnorm(n)
  m <- sin(2 * pi * z)
  x <- m + e1; y <- m + 0.5 * e1 + e2
  Fit <- function(zt, xt, ze) predict(smooth.spline(zt, xt, df = df), ze)$y
  in_sample <- mean((x - Fit(z, x, z)) * (y - Fit(z, y, z)))
  fold <- sample(rep(seq_len(K), length.out = n))
  r <- numeric(n)
  for (k in seq_len(K)) {
    tr <- fold != k; te <- fold == k
    r[te] <- (x[te] - Fit(z[tr], x[tr], z[te])) * (y[te] - Fit(z[tr], y[tr], z[te]))
  }
  c(in_sample = in_sample, cross_fit = mean(r))
}

# Bias is measured against the truth (0.5) and compared with the SD of the estimator.
Summarize <- function(n, R = 150) {
  res <- replicate(R, OneRep(n))
  bias <- rowMeans(res) - 0.5
  sd_hat <- apply(res, 1, sd)
  data.frame(n = n, estimator = rownames(res), bias = round(bias, 3),
             bias_mcse = round(sd_hat / sqrt(R), 4),
             sd = round(sd_hat, 3), bias_over_sd = round(bias / sd_hat, 2), row.names = NULL)
}
do.call(rbind, lapply(c(200, 800, 3200), Summarize))

     n estimator   bias bias_mcse    sd bias_over_sd
1  200 in_sample -0.042    0.0070 0.086        -0.48
2  200 cross_fit  0.044    0.0083 0.102         0.43
3  800 in_sample -0.008    0.0036 0.044        -0.19
4  800 cross_fit  0.013    0.0038 0.047         0.27
5 3200 in_sample -0.002    0.0017 0.021        -0.09
6 3200 cross_fit  0.004    0.0017 0.021         0.17

Read bias against its Monte Carlo SE and the estimator’s sampling SD. The common-smoother trace can explain downward in-sample noise absorption. Cross-fitting can leave positive bias because the two estimated regression errors are correlated through shared training noise.

This experiment does not prove root-n inference. Keeping the spline’s nominal degrees of freedom fixed does not generally make its approximation error vanish for every nonlinear mean. Three sample sizes do not establish an asymptotic rate, and changing the fold count need not decrease bias monotonically. Tuning, approximation, moments, and the full remainder still matter.

5.1 A counterexample: separated folds with a bad nuisance learner

Here both true conditional means are zero and the true residual covariance is 0.5. A one-nearest-neighbor fit interpolates its training outcomes, giving zero in-sample residual products. Out of fold, it predicts a different person’s noise. Because the same neighbor supplies both predictions, their errors have covariance 0.5: the cross-fitted residual product has expectation 1, not 0.5. The learner’s nuisance error does not vanish. A training-set mean is a consistent nuisance fit in this deliberately flat model.

Code

set.seed(22)
BadLearner <- function(n = 250, K = 5) {
  z <- runif(n); x <- rnorm(n); y <- 0.5*x + rnorm(n)
  fold <- sample(rep(seq_len(K), length.out = n))
  nearest <- constant <- numeric(n)
  for (k in seq_len(K)) {
    tr <- which(fold != k); te <- which(fold == k)
    neighbor <- tr[max.col(-outer(z[te], z[tr], "-")^2, ties.method = "first")]
    nearest[te] <- (x[te]-x[neighbor])*(y[te]-y[neighbor])
    constant[te] <- (x[te]-mean(x[tr]))*(y[te]-mean(y[tr]))
  }
  c(in_sample_1NN = 0, cross_fit_1NN = mean(nearest), cross_fit_mean = mean(constant))
}
bad_fit <- replicate(120, BadLearner())
boxplot(t(bad_fit), names = c("In-sample\n1NN", "Cross-fit\n1NN", "Cross-fit\nmean"),
        col = c("gray80", "#D55E0033", "#2F80C833"), ylab = "Estimated residual covariance")
abline(h = 0.5, lty = 2)

Boxplots compare an in-sample nearest-neighbor estimate near zero, a cross-fitted nearest-neighbor estimate near one, and a cross-fitted training-mean estimate near the true covariance one half. — Figure 2: An interpolation counterexample with 120 independent synthetic datasets. Separating folds removes own-outcome reuse, but a one-nearest-neighbor nuisance still leaves population bias. A training mean works here because the true means are zero.

A finite-sample spread in the plot does not change the exact expectation argument. To see the mechanism directly, the held-out person and their selected training neighbor have independent residual pairs, each with covariance 0.5. Expanding (X_i-X_{\rm neighbor})(Y_i-Y_{\rm neighbor}) leaves two covariance terms, totaling 1, because predicting with one neighbor retains its noise even as the training set grows. Predicting with a training mean instead adds only 0.5/n_{\rm train} in expectation. That extra term vanishes as the training set grows. Cross-fitting removes self-use in both cases, but only the second learner estimates the zero conditional means consistently in this example.

6 The fold includes the fitting pipeline

For an independent-subject analysis, put all observations from one subject in the same fold. Estimate outcome-dependent preprocessing, feature selection, tuning, and nuisance regressions within the training data when they belong to the fitting pipeline; then apply those fitted operations to held-out observations. A fixed transformation chosen independently of the analysis can be applied consistently across folds. Shared preprocessing learned from all observations needs its own justification and can otherwise leak evaluation information into training.

For this average covariance, the held-out contribution and estimator are

\begin{aligned} R_i&=\{X_i-\hat m_X^{(-k)}(Z_i)\} \times\{Y_i-\hat m_Y^{(-k)}(Z_i)\},\\ \hat\psi&=n^{-1}\sum_iR_i. \end{aligned}

The estimated influence contributions are R_i-\hat\psi. For a residual correlation, estimate the covariance and both residual variances jointly; averaging stratum-specific correlations is a different target. Put v_X=E[U_X^2], v_Y=E[U_Y^2], and \rho=\psi/\sqrt{v_Xv_Y}, with both variances positive. The delta method gives

\phi_\rho=\frac{U_XU_Y-\psi}{\sqrt{v_Xv_Y}} -\frac{\rho}{2}\left\{ \frac{U_X^2-v_X}{v_X}+\frac{U_Y^2-v_Y}{v_Y}\right\}.

The first term tracks covariance uncertainty; the second tracks uncertainty in both normalizers. Replace all components by held-out estimates, then take the SD of the resulting influence values divided by \sqrt n under the inference conditions below. Reusing only the covariance SE would omit denominator uncertainty. Also, an exact m_X alone does not repair the Y variance: its population bias is E[\delta_Y^2]. The covariance’s either-mean-correct property therefore does not automatically transfer to correlation.

7 Conditions connecting the formula to an SE

Let D_0=(X-m_X)(Y-m_Y), and for fold k let H_k=D_{\hat m_k}-D_0. With w_k=n_k/n,

\hat\psi-\psi=(P_n-P)(D_0-\psi) +\sum_k w_k(P_{n,k}-P)H_k +\sum_k w_kE[\Delta_{X,k}\Delta_{Y,k}].

For i.i.d. units and a fixed number of nonvanishing folds, require 0<\operatorname{Var}(D_0)<\infty, \|H_k\|_2=o_p(1), and the product bias o_p(n^{-1/2}) in every fold. Conditional variance within each held-out fold controls the middle term; the proof in From an influence function to an estimator shows that step. Do not assume joint conditional independence across all folds.

Under these controls, averaging the true contributions gives a CLT and the sample SD of the estimated contributions divided by \sqrt n consistently estimates the SE. L2 consistency of mean fits alone need not imply L2 consistency of their products. Bounded outcomes/fits plus L2 nuisance consistency are one sufficient setup; suitable L4 moments and L4 nuisance consistency are another route. See Estimating nuisance functions for inference for a concrete learning-rate calculation.

8 Putting it together: the recipe we use

For an estimand with nuisances, the group’s estimators are built in this order.

Write the estimand as a functional \psi(F) and derive its efficient influence function.
Estimate the nuisances flexibly, with cross-fitting.
Construct the orthogonal formula or EIF estimating equation. Verify any claimed equivalence with a one-step or targeted construction for this particular target.
Standard error: for independent units and valid remainder conditions, use the sample standard deviation of estimated EIF values divided by \sqrt n. With clustered observations, aggregate contributions and use a design-appropriate variance calculation. Orthogonality alone does not establish these conditions.

For sample-specific network summaries, define the population quantity and derive the relevant correction before assigning an observation-level interpretation. The Jackknife pseudo-values page explains why simply recentering a cohort-derived score is not by itself an inferential argument.

9 Words you will meet

Double robustness. A separate robustness property of some estimating equations with two nuisances (an outcome regression and a propensity or censoring model) where the estimator is consistent if either nuisance is correct. The residual form of \psi has this property: its remainder is the product E[h_X h_Y], so it is zero when either m_X or m_Y is exactly right, whatever the other is (see the check above). This same either-correct guarantee does not generally transfer to residual correlation: its two residual-variance denominators each have a squared-error remainder for their own mean fit.
Targeted learning, TMLE, double machine learning (DML). Related frameworks that can use orthogonal estimating constructions and sample splitting. Targeted learning emphasizes distributional updates, while DML emphasizes orthogonal moments with learned nuisances; neither label alone specifies a complete algorithm or proof. Some valid procedures avoid splitting by imposing different empirical-process conditions. For the same target and efficiency conditions, particular estimators can share an efficient first-order expansion. The frameworks are not interchangeable in every model or finite sample.

10 Exercises

The checkpoint is this page’s completion target. Additional practice is optional. Use the missing-mean argument as a reference, but derive the covariance remainder rather than assuming it is unchanged.

Checkpoint · COV-01 — Separate nuisance bias from held-out fluctuation

Task type: Derivation and rate reasoning.

Let m_X(Z)=E[X\mid Z], m_Y(Z)=E[Y\mid Z], and \psi=E[(X-m_X(Z))(Y-m_Y(Z))]. Hold training fits fixed, with errors \delta_X=\hat m_X-m_X and \delta_Y=\hat m_Y-m_Y. Assume the variables and fitted means have finite second moments.

Derive the population bias of the residual-product formula and of the difference formula E[XY]-E[\hat m_X(Z)\hat m_Y(Z)]. Which formula removes the terms involving only one nuisance error?
Take both errors to be the positive constant n^{-1/5}. Calculate the residual-product bias and its size after multiplication by \sqrt n. Would evaluating on held-out people remove it?
State a sufficient product-of-L2-errors condition for negligible root-n bias. Identify the separate fitted-minus-true empirical term and variance conditions still needed before reporting an influence-function SE.

Ready to continue when: you can derive the covariance nuisance-bias remainder and distinguish its rate requirement from cross-fitting and variance consistency.

Transfer: If one conditional mean is exactly correct and the other is misspecified, what happens to the population bias? Does that alone establish the efficient-SE argument?

Hint

Write each fitted residual as its true residual minus \delta. Condition on Z to remove the single-error terms in the residual product. For the rate comparison, multiply the product bias by \sqrt n before judging whether it disappears.

Solution outline — open after trying

Writing U_X=X-m_X(Z) and U_Y=Y-m_Y(Z) gives

E[(U_X-\delta_X)(U_Y-\delta_Y)]-\psi =E[\delta_X\delta_Y],

because both true residuals have conditional mean zero given Z. The difference formula instead has bias -E[m_X\delta_Y+m_Y\delta_X+\delta_X\delta_Y], which generally includes single-error terms.

With both errors equal to n^{-1/5}, the residual-product bias is n^{-2/5} and its scaled size is n^{1/10}, which diverges. Cross-fitting does not change this population expectation.

Cauchy–Schwarz bounds the bias by \|\delta_X\|_2\|\delta_Y\|_2, so a product o_p(n^{-1/2}) is sufficient. The held-out empirical term is (P_{n,k}-P)H_k, where H_k is the fitted-minus-true residual product. Its L2 size must vanish under the stated fold conditions. A CLT for the true contribution and consistent estimation of its finite, positive variance are also required; L2 consistency of the mean fits alone need not control their product’s L2 size.

If either error is zero, the population bias is zero. The remaining fit can still change the leading contribution and its variance, so this alone does not establish the efficient-SE argument.

Additional practice (optional)

These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.

Show that the variance \sigma^2 = E[(X - \mu)^2], viewed as a formula with nuisance \mu, is Neyman orthogonal in \mu. Then show that the formula E[X^2] - \mu^2 is not. Check both numerically by inflating \mu by 1 + \varepsilon.
In the orthogonality simulation, replace the inflation error with an additive error m + \varepsilon(Z - 1/2) and predict, before running, the leading term of each formula’s bias.
Change df in the cross-fitting example to 5 and to 40 and re-run at N = 200. How do the two biases change, and which grows faster with df? Then replace the smoothing spline with an average of 20 deep regression trees, each grown on a bootstrap resample (rpart is a recommended R package; use rpart.control() with cp=0, minsplit=5, and minbucket=3). Compare the in-sample and cross-fit biases. Why does a deep tree fit its own noise?
The residual correlation is \psi / \sqrt{\psi_{XX}\,\psi_{YY}} with \psi_{XX} = E[\operatorname{Var}(X \mid Z)]. Write the orthogonal estimators of the three pieces and combine them. What does the delta method say about the standard error?

11 Reading

The DML paper and Kennedy review provide the formal conditions behind the constructions.

Hines, O., Dukes, O., Diaz-Ordaz, K., and Vansteelandt, S. (2022). Demystifying statistical learning based on efficient influence functions. The American Statistician, 76(3), 292 to 304. An introductory derivation-focused treatment.
Kennedy, E. H. (2022). Semiparametric doubly robust targeted double machine learning: a review. arXiv:2203.06469. Influence-function construction, remainder decomposition, and efficiency conditions.
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), C1 to C68. Section 2 for Neyman orthogonality and cross-fitting.
Vansteelandt, S. and Dukes, O. (2022). Assumption-lean inference for generalised linear model parameters. Journal of the Royal Statistical Society B, 84(3), 657 to 685. Asks what a regression coefficient estimates when the model is wrong, defines an estimand that does not depend on the model being right, and builds the orthogonal estimator for it. The target remains defined independently of the correctness of the working regression model.
Tsiatis, A. A. (2006). Semiparametric Theory and Missing Data. Springer. Chapters 3 and 4 for the geometry, when you want the full version.

--- title: "Residual covariance and nuisance-estimation error" subtitle: "Population nuisance bias, held-out evaluation, and conditions for first-order inference" guide-status: "Available" --- **Before reading:** [Conditional expectation and L2 geometry](../statistical-models/conditional-expectation.qmd), [Model restrictions and efficiency](../statistical-models/models-and-efficiency.qmd), and [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd). This application reuses the capstone's nuisance-bias and held-out empirical-term arguments for a covariance target. ## The problem this page solves The [Model restrictions and efficiency](../statistical-models/models-and-efficiency.qmd) page introduced an asymptotic efficiency benchmark and showed how its influence function can depend on unknown nuisance quantities. Here the nuisance is a regression function: how a gene's mean, or variance, depends on covariates such as cell composition, age, batch, or an interferon score. The estimand here is the average conditional covariance between two genes after removing their conditional means given $Z$, rather than a causal or conditional-independence claim. We do not want to assume the regressions are linear, so we estimate them flexibly, and then we need to know how the error in those flexible estimates propagates into the estimand. For this target, an orthogonal residual-product equation removes the first-order population nuisance error. Its exact remaining bias and its empirical fluctuation still need control. This transfer application reuses Neyman orthogonality and cross-fitting for a different target. The task is to determine which derivatives and remainder terms change for covariance. <span id="the-efficient-influence-function-briefly"></span> ## The conditional covariance, two ways Let $Z$ be a covariate (say, a cell-composition score), and let $m_X(z) = E[X \mid Z = z]$ and $m_Y(z) = E[Y \mid Z = z]$. The estimand is the average conditional covariance $$ \psi = E\big[\operatorname{Cov}(X, Y \mid Z)\big]. $$ This is the covariance between two genes after removing what each shares with $Z$, averaged over the population. There are at least two correct formulas for it: $$ \psi = E\big[(X-m_X(Z))(Y-m_Y(Z))\big] $$ $$ \psi = E[XY]-E[m_X(Z)m_Y(Z)] $$ The first is the **residual form** and the second the **difference form**. Both are exactly right if $m_X$ and $m_Y$ are the true regression functions. Both can be turned into estimators by plugging in estimates $\hat m_X$ and $\hat m_Y$ and averaging. They behave completely differently when the estimates are imperfect. To verify their equality, expand the residual product. The tower property gives $E[Xm_Y(Z)]=E[m_X(Z)m_Y(Z)]$, and the same identity holds with X and Y reversed. Two subtractions and one addition of that mean-product term leave the difference form. Algebraic equivalence **at the truth** does not imply equal sensitivity to estimated inputs. Write $W=(X,Y,Z)$ for one person's record. Throughout, $P f=E[f(W)]$ denotes a population average, $P_n f=n^{-1}\sum_i f(W_i)$ a sample average, and $\|f\|_2=\{E[f(W)^2]\}^{1/2}$ the population root-mean-square size. Existence of the covariance needs second moments of X and Y. A CLT for their residual product needs the product itself to have finite variance; suitable fourth moments are a convenient stronger condition. ## Neyman orthogonality We want a formula whose population value changes little when a fitted conditional mean is slightly wrong. Write the trained errors as $\delta_X=\hat m_X-m_X$ and $\delta_Y=\hat m_Y-m_Y$. A fitted residual is $U_X-\delta_X$, so the exact expansion is $$ \begin{aligned} E[(U_X-\delta_X)(U_Y-\delta_Y)]&=\psi-E[U_X\delta_Y]-E[U_Y\delta_X] +E[\delta_X\delta_Y]\\ &=\psi+E[\delta_X\delta_Y]. \end{aligned} $$ The two single-error terms vanish because each error is a function of Z and each true residual has conditional mean zero given Z. Only their product remains. In contrast, plugging into the difference form leaves $-E[m_X\delta_Y+m_Y\delta_X+\delta_X\delta_Y]$, including single-error terms. This exact comparison motivates the derivative definition that follows. Hold the population law $P$ fixed and write the population estimating map as $M(m)$, with nuisance arguments $m$ and true pair $m_0=(m_X,m_Y)$ satisfying $M(m_0)=\psi$. The target is fixed; we are varying the nuisance arguments supplied to the formula. Perturb them along an admissible direction $h$. The equation is **Neyman orthogonal** when $$ \frac{\partial}{\partial\varepsilon}M(m_0+\varepsilon h) \Big|_{\varepsilon=0}=0. $$ The derivative is taken at the true nuisance for the permitted directions. A quadratic remainder additionally needs smoothness and norm controls; a zero derivative alone is not a universal rate theorem. For the residual product below, the second-order term follows from exact algebra. **Residual form.** Write $U_X=X-m_X(Z)$ and $U_Y=Y-m_Y(Z)$ for the true conditional-mean residuals. Perturb $m_X$ by $\varepsilon h_X(Z)$: $$ \begin{aligned} \frac{\partial}{\partial\varepsilon}E[(U_X-\varepsilon h_X)U_Y]\Big|_0&=-E[h_X(Z)U_Y]\\ &=-E[h_X(Z)E(U_Y\mid Z)]\\ &=0, \end{aligned} $$ because $E[Y - m_Y(Z) \mid Z] = 0$ by definition of $m_Y$. The same holds for $m_Y$. The residual form is orthogonal. Its second-order term is $E[h_X(Z) h_Y(Z)]\,\varepsilon^2$: the product of the two nuisance errors. **Difference form.** Perturb $m_X$: $$ \frac{\partial}{\partial\varepsilon} \{E[XY]-E[(m_X+\varepsilon h_X)m_Y]\}\Big|_0=-E[h_X(Z)m_Y(Z)], $$ This derivative is generally nonzero, so the difference-form plug-in is not generally orthogonal. Particular directions can still give zero derivative; orthogonality requires the condition for every relevant direction. Here it is numerically. Inject a known error into the nuisance (inflate the true regression function by a factor $1 + \varepsilon$) and watch the two formulas. ```{r} set.seed(1) n <- 20000 z <- runif(n); e1 <- rnorm(n); e2 <- rnorm(n) x <- z + e1 # m_X(z) = z y <- z + 0.5 * e1 + e2 # m_Y(z) = z; Cov(X, Y | Z) = 0.5 exactly m_x <- z; m_y <- z # the true nuisances Residual <- function(mx, my) mean((x - mx) * (y - my)) Difference <- function(mx, my) mean(x * y) - mean(mx * my) eps <- c(0, 0.05, 0.1, 0.2, 0.4) out <- sapply(eps, function(e) c(residual_form = Residual((1 + e) * m_x, (1 + e) * m_y), difference_form = Difference((1 + e) * m_x, (1 + e) * m_y))) colnames(out) <- paste0("eps=", eps) round(out, 3) ``` With $m_X = m_Y = Z$ and $Z$ uniform, the residual form's error is $\varepsilon^2 E[Z^2] = \varepsilon^2/3$ (about 0.003 at $\varepsilon = 0.1$) and the difference form's is $-2\varepsilon/3 - \varepsilon^2/3$ (about $-0.07$ at $\varepsilon = 0.1$). Those are population biases for this perturbation. The displayed sample averages also contain Monte Carlo fluctuation; these numerical biases are specific to the perturbation. ```{r} #| label: fig-orthogonal-sensitivity #| fig-cap: "Exact population bias in the stated synthetic perturbation. The residual formula changes quadratically; the difference formula has a nonzero first-order slope. These curves are not sampled estimates." #| fig-alt: "Bias curves show a flat tangent at zero for the residual form and a downward linear component for the difference form." grid <- seq(-0.3, 0.3, length.out = 201) matplot(grid, cbind(residual = grid^2/3, difference = -(2*grid + grid^2)/3), type = "l", lty = c(1, 2), col = c("#2F80C8", "#D55E00"), xlab = "Nuisance perturbation epsilon", ylab = "Population bias") abline(h = 0, v = 0, col = "gray70") legend("topright", c("Residual form", "Difference form"), lty = 1:2, col = c("#2F80C8", "#D55E00"), bty = "n") ``` The remainder of the residual form is a *product* of the two nuisance errors, so it also vanishes when either regression function is exactly right. Keep $m_X$ exact and make $m_Y$ badly wrong: ```{r} m_y_wrong <- 3 * z + 1 # a badly wrong nuisance function round(c(residual_form = Residual(m_x, m_y_wrong), # still near 0.5 (Monte Carlo SE about 0.02) difference_form = Difference(m_x, m_y_wrong)), # far from 0.5 3) ``` The important condition is on the product of nuisance errors, not on the name of the learning algorithm. For root-$n$ inference centered at the target, a usual sufficient condition here is $$ \|\delta_X\|_2\,\|\delta_Y\|_2=o_p(n^{-1/2}), $$ along with moment, consistency, and sampling conditions. Errors of exactly order $n^{-1/4}$ each give a product of order $n^{-1/2}$, which need not be negligible; both being $o_p(n^{-1/4})$ is sufficient, and faster convergence of one can compensate for slower convergence of the other. Orthogonality removes the first-order population term, while cross-fitting helps control the empirical term. Neither promises valid inference for arbitrary nuisance fits. [Chernozhukov et al.](https://arxiv.org/abs/1608.00060) develop this distinction. ## The connection to the efficient influence function Efficient influence functions provide a systematic route to orthogonal equations. Orthogonality to nuisance *scores* is a model-geometric statement; Neyman orthogonality concerns derivatives of the chosen estimating equation in its nuisance arguments. For a valid influence-function construction these properties are connected under regularity conditions. We can verify the nuisance derivative directly for this residual-product equation. For $\psi$ above, the EIF is $$ \operatorname{EIF}(X,Y,Z)=U_XU_Y-\psi, $$ and setting its sample average to zero gives the residual form. The naive difference-form plug-in has first-order nuisance sensitivity. It is not automatically an estimator with a different valid influence function: its behavior also depends on how the regressions are estimated. So the recipe for an orthogonal estimator is: find the EIF, and estimate by solving its empirical mean equals zero. The **one-step estimator** adds an estimated influence-function correction to an initial estimate. Writing $P_n$ for a sample average, $$ \hat\psi_{\rm step}=\hat\psi_{\rm initial} +P_n\phi(W;\hat\psi_{\rm initial},\hat m). $$ Here the centering parameter is included in the estimated influence function. For this residual-covariance equation the initial centering cancels algebraically, leaving the residual-product average. Equivalence to an estimating-equation solution must be checked for the particular target and construction, including how the initial target estimate enters the estimated EIF. Why is this the influence function of a target that itself contains conditional means? Along a smooth path with score $s$, differentiating $E_t[U_{X,t}U_{Y,t}]$ gives a law-change term $E[U_XU_Ys]$ and two mean-change terms. Each mean-change term vanishes by conditional residual centering, just as in the calculation above. Since $E[s]=0$, the derivative is $E[(U_XU_Y-\psi)s]$. In the unrestricted model this identifies the canonical gradient. It is a derivative of the **whole functional**, including its changing conditional means, rather than a calculation that assumes those means are known. Deriving EIFs is a skill; Hines et al. (2022) teach it with the Gateaux-derivative method from the [Influence functions](../statistical-models/influence-functions.qmd) page, and Kennedy (2022) catalogs common examples. Residual correlation requires a further normalization step below. ## Cross-fitting Cross-fitting separates a held-out observation from its own training fit. It does not make nuisance predictions correct or all held-out contributions jointly independent: different folds share training observations. One more problem. If $\hat m_X$ is fit on the same observations it is then evaluated on, the fit absorbs some of each observation's own noise. The residual $X_i - \hat m_X(Z_i)$ is then too small, and because $Y_i - \hat m_Y(Z_i)$ is shrunk in a correlated way, the residual product is biased. For a fixed common linear smoother $S$, an error-covariance contribution is scaled by $\operatorname{tr}\{(I-S)^\top(I-S)\}/N$; a projection fit loses its rank in this trace. A smoothing spline is not generally a projection, and adaptive fits need their own analysis. The direction and size of the total finite-sample bias also depend on approximation error. **Cross-fitting** prevents a held-out outcome from being used to train its own prediction. Split the sample into $K$ folds. For each fold, fit the nuisances on the other $K - 1$ folds and evaluate the residuals on the held-out fold. Average the residual products over all folds. Each observation's residual now uses a nuisance fit that never saw it. ```{r} set.seed(2) # Nonlinear nuisance so that a smoother is needed. Cov(X, Y | Z) = 0.5 again. OneRep <- function(n, df = 15, K = 5) { z <- runif(n); e1 <- rnorm(n); e2 <- rnorm(n) m <- sin(2 * pi * z) x <- m + e1; y <- m + 0.5 * e1 + e2 Fit <- function(zt, xt, ze) predict(smooth.spline(zt, xt, df = df), ze)$y in_sample <- mean((x - Fit(z, x, z)) * (y - Fit(z, y, z))) fold <- sample(rep(seq_len(K), length.out = n)) r <- numeric(n) for (k in seq_len(K)) { tr <- fold != k; te <- fold == k r[te] <- (x[te] - Fit(z[tr], x[tr], z[te])) * (y[te] - Fit(z[tr], y[tr], z[te])) } c(in_sample = in_sample, cross_fit = mean(r)) } # Bias is measured against the truth (0.5) and compared with the SD of the estimator. Summarize <- function(n, R = 150) { res <- replicate(R, OneRep(n)) bias <- rowMeans(res) - 0.5 sd_hat <- apply(res, 1, sd) data.frame(n = n, estimator = rownames(res), bias = round(bias, 3), bias_mcse = round(sd_hat / sqrt(R), 4), sd = round(sd_hat, 3), bias_over_sd = round(bias / sd_hat, 2), row.names = NULL) } do.call(rbind, lapply(c(200, 800, 3200), Summarize)) ``` Read bias against its Monte Carlo SE and the estimator's sampling SD. The common-smoother trace can explain downward in-sample noise absorption. Cross-fitting can leave positive bias because the two estimated regression errors are correlated through shared training noise. **This experiment does not prove root-$n$ inference.** Keeping the spline's nominal degrees of freedom fixed does not generally make its approximation error vanish for every nonlinear mean. Three sample sizes do not establish an asymptotic rate, and changing the fold count need not decrease bias monotonically. Tuning, approximation, moments, and the full remainder still matter. ### A counterexample: separated folds with a bad nuisance learner Here both true conditional means are zero and the true residual covariance is 0.5. A one-nearest-neighbor fit interpolates its training outcomes, giving zero in-sample residual products. Out of fold, it predicts a different person's noise. Because the same neighbor supplies both predictions, their errors have covariance 0.5: the cross-fitted residual product has expectation 1, not 0.5. The learner's nuisance error does not vanish. A training-set mean is a consistent nuisance fit in this deliberately flat model. ```{r} #| label: fig-crossfit-counterexample #| fig-cap: "An interpolation counterexample with 120 independent synthetic datasets. Separating folds removes own-outcome reuse, but a one-nearest-neighbor nuisance still leaves population bias. A training mean works here because the true means are zero." #| fig-alt: "Boxplots compare an in-sample nearest-neighbor estimate near zero, a cross-fitted nearest-neighbor estimate near one, and a cross-fitted training-mean estimate near the true covariance one half." set.seed(22) BadLearner <- function(n = 250, K = 5) { z <- runif(n); x <- rnorm(n); y <- 0.5*x + rnorm(n) fold <- sample(rep(seq_len(K), length.out = n)) nearest <- constant <- numeric(n) for (k in seq_len(K)) { tr <- which(fold != k); te <- which(fold == k) neighbor <- tr[max.col(-outer(z[te], z[tr], "-")^2, ties.method = "first")] nearest[te] <- (x[te]-x[neighbor])*(y[te]-y[neighbor]) constant[te] <- (x[te]-mean(x[tr]))*(y[te]-mean(y[tr])) } c(in_sample_1NN = 0, cross_fit_1NN = mean(nearest), cross_fit_mean = mean(constant)) } bad_fit <- replicate(120, BadLearner()) boxplot(t(bad_fit), names = c("In-sample\n1NN", "Cross-fit\n1NN", "Cross-fit\nmean"), col = c("gray80", "#D55E0033", "#2F80C833"), ylab = "Estimated residual covariance") abline(h = 0.5, lty = 2) ``` A finite-sample spread in the plot does not change the exact expectation argument. To see the mechanism directly, the held-out person and their selected training neighbor have independent residual pairs, each with covariance 0.5. Expanding $(X_i-X_{\rm neighbor})(Y_i-Y_{\rm neighbor})$ leaves **two** covariance terms, totaling 1, because predicting with one neighbor retains its noise even as the training set grows. Predicting with a training mean instead adds only $0.5/n_{\rm train}$ in expectation. That extra term vanishes as the training set grows. Cross-fitting removes self-use in both cases, but only the second learner estimates the zero conditional means consistently in this example. ## The fold includes the fitting pipeline For an independent-subject analysis, put all observations from one subject in the same fold. Estimate outcome-dependent preprocessing, feature selection, tuning, and nuisance regressions within the training data when they belong to the fitting pipeline; then apply those fitted operations to held-out observations. A fixed transformation chosen independently of the analysis can be applied consistently across folds. Shared preprocessing learned from all observations needs its own justification and can otherwise leak evaluation information into training. For this average covariance, the held-out contribution and estimator are $$ \begin{aligned} R_i&=\{X_i-\hat m_X^{(-k)}(Z_i)\} \times\{Y_i-\hat m_Y^{(-k)}(Z_i)\},\\ \hat\psi&=n^{-1}\sum_iR_i. \end{aligned} $$ The estimated influence contributions are $R_i-\hat\psi$. For a residual correlation, estimate the covariance and both residual variances jointly; averaging stratum-specific correlations is a different target. Put $v_X=E[U_X^2]$, $v_Y=E[U_Y^2]$, and $\rho=\psi/\sqrt{v_Xv_Y}$, with both variances positive. The delta method gives $$ \phi_\rho=\frac{U_XU_Y-\psi}{\sqrt{v_Xv_Y}} -\frac{\rho}{2}\left\{ \frac{U_X^2-v_X}{v_X}+\frac{U_Y^2-v_Y}{v_Y}\right\}. $$ The first term tracks covariance uncertainty; the second tracks uncertainty in both normalizers. Replace all components by held-out estimates, then take the SD of the resulting influence values divided by $\sqrt n$ under the inference conditions below. Reusing only the covariance SE would omit denominator uncertainty. Also, an exact $m_X$ alone does not repair the Y variance: its population bias is $E[\delta_Y^2]$. The covariance's either-mean-correct property therefore does not automatically transfer to correlation. ## Conditions connecting the formula to an SE Let $D_0=(X-m_X)(Y-m_Y)$, and for fold $k$ let $H_k=D_{\hat m_k}-D_0$. With $w_k=n_k/n$, $$ \hat\psi-\psi=(P_n-P)(D_0-\psi) +\sum_k w_k(P_{n,k}-P)H_k +\sum_k w_kE[\Delta_{X,k}\Delta_{Y,k}]. $$ For i.i.d. units and a fixed number of nonvanishing folds, require $0<\operatorname{Var}(D_0)<\infty$, $\|H_k\|_2=o_p(1)$, and the product bias $o_p(n^{-1/2})$ in every fold. Conditional variance within each held-out fold controls the middle term; the proof in [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd#checkpoint-sp-12) shows that step. Do not assume joint conditional independence across all folds. Under these controls, averaging the true contributions gives a CLT and the sample SD of the estimated contributions divided by $\sqrt n$ consistently estimates the SE. L2 consistency of mean fits alone need not imply L2 consistency of their products. Bounded outcomes/fits plus L2 nuisance consistency are one sufficient setup; suitable L4 moments and L4 nuisance consistency are another route. See [Estimating nuisance functions for inference](../semiparametric-methods/nuisance-estimation.qmd) for a concrete learning-rate calculation. ## Putting it together: the recipe we use For an estimand with nuisances, the group's estimators are built in this order. 1. Write the estimand as a functional $\psi(F)$ and derive its efficient influence function. 2. Estimate the nuisances flexibly, with cross-fitting. 3. Construct the orthogonal formula or EIF estimating equation. Verify any claimed equivalence with a one-step or targeted construction for this particular target. 4. Standard error: for independent units and valid remainder conditions, use the sample standard deviation of estimated EIF values divided by $\sqrt n$. With clustered observations, aggregate contributions and use a design-appropriate variance calculation. Orthogonality alone does not establish these conditions. For sample-specific network summaries, define the population quantity and derive the relevant correction before assigning an observation-level interpretation. The [Jackknife pseudo-values](jackknife-pseudo-values.qmd) page explains why simply recentering a cohort-derived score is not by itself an inferential argument. ## Words you will meet - **Double robustness.** A separate robustness property of some estimating equations with two nuisances (an outcome regression and a propensity or censoring model) where the estimator is consistent if *either* nuisance is correct. The residual form of $\psi$ has this property: its remainder is the product $E[h_X h_Y]$, so it is zero when either $m_X$ or $m_Y$ is exactly right, whatever the other is (see the check above). This same either-correct guarantee does not generally transfer to residual *correlation*: its two residual-variance denominators each have a squared-error remainder for their own mean fit. - **Targeted learning, TMLE, double machine learning (DML).** Related frameworks that can use orthogonal estimating constructions and sample splitting. Targeted learning emphasizes distributional updates, while DML emphasizes orthogonal moments with learned nuisances; neither label alone specifies a complete algorithm or proof. Some valid procedures avoid splitting by imposing different empirical-process conditions. For the same target and efficiency conditions, particular estimators can share an efficient first-order expansion. The frameworks are not interchangeable in every model or finite sample. ## Exercises **The checkpoint is this page's completion target. Additional practice is optional.** Use the missing-mean argument as a reference, but derive the covariance remainder rather than assuming it is unchanged. ::: {#checkpoint-cov-01 .checkpoint} ### Checkpoint · COV-01 — Separate nuisance bias from held-out fluctuation {.unnumbered} **Task type:** Derivation and rate reasoning. Let $m_X(Z)=E[X\mid Z]$, $m_Y(Z)=E[Y\mid Z]$, and $\psi=E[(X-m_X(Z))(Y-m_Y(Z))]$. Hold training fits fixed, with errors $\delta_X=\hat m_X-m_X$ and $\delta_Y=\hat m_Y-m_Y$. Assume the variables and fitted means have finite second moments. 1. Derive the population bias of the residual-product formula and of the difference formula $E[XY]-E[\hat m_X(Z)\hat m_Y(Z)]$. Which formula removes the terms involving only one nuisance error? 2. Take both errors to be the positive constant $n^{-1/5}$. Calculate the residual-product bias and its size after multiplication by $\sqrt n$. Would evaluating on held-out people remove it? 3. State a sufficient product-of-L2-errors condition for negligible root-$n$ bias. Identify the separate fitted-minus-true empirical term and variance conditions still needed before reporting an influence-function SE. **Ready to continue when:** you can derive the covariance nuisance-bias remainder and distinguish its rate requirement from cross-fitting and variance consistency. **Transfer:** If one conditional mean is exactly correct and the other is misspecified, what happens to the population bias? Does that alone establish the efficient-SE argument? <details class="checkpoint-help"> <summary>Hint</summary> Write each fitted residual as its true residual minus $\delta$. Condition on $Z$ to remove the single-error terms in the residual product. For the rate comparison, multiply the product bias by $\sqrt n$ before judging whether it disappears. </details> <details class="checkpoint-help"> <summary>Solution outline — open after trying</summary> Writing $U_X=X-m_X(Z)$ and $U_Y=Y-m_Y(Z)$ gives $$ E[(U_X-\delta_X)(U_Y-\delta_Y)]-\psi =E[\delta_X\delta_Y], $$ because both true residuals have conditional mean zero given $Z$. The difference formula instead has bias $-E[m_X\delta_Y+m_Y\delta_X+\delta_X\delta_Y]$, which generally includes single-error terms. With both errors equal to $n^{-1/5}$, the residual-product bias is $n^{-2/5}$ and its scaled size is $n^{1/10}$, which diverges. Cross-fitting does not change this population expectation. Cauchy–Schwarz bounds the bias by $\|\delta_X\|_2\|\delta_Y\|_2$, so a product $o_p(n^{-1/2})$ is sufficient. The held-out empirical term is $(P_{n,k}-P)H_k$, where $H_k$ is the fitted-minus-true residual product. Its L2 size must vanish under the stated fold conditions. A CLT for the true contribution and consistent estimation of its finite, positive variance are also required; L2 consistency of the mean fits alone need not control their product's L2 size. If either error is zero, the population bias is zero. The remaining fit can still change the leading contribution and its variance, so this alone does not establish the efficient-SE argument. </details> ::: ### Additional practice (optional) {.unnumbered} These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work. 1. Show that the variance $\sigma^2 = E[(X - \mu)^2]$, viewed as a formula with nuisance $\mu$, is Neyman orthogonal in $\mu$. Then show that the formula $E[X^2] - \mu^2$ is not. Check both numerically by inflating $\mu$ by $1 + \varepsilon$. 2. In the orthogonality simulation, replace the inflation error with an additive error $m + \varepsilon(Z - 1/2)$ and predict, before running, the leading term of each formula's bias. 3. Change `df` in the cross-fitting example to 5 and to 40 and re-run at $N = 200$. How do the two biases change, and which grows faster with `df`? Then replace the smoothing spline with an average of 20 deep regression trees, each grown on a bootstrap resample (`rpart` is a recommended R package; use `rpart.control()` with `cp=0`, `minsplit=5`, and `minbucket=3`). Compare the in-sample and cross-fit biases. Why does a deep tree fit its own noise? 4. The residual correlation is $\psi / \sqrt{\psi_{XX}\,\psi_{YY}}$ with $\psi_{XX} = E[\operatorname{Var}(X \mid Z)]$. Write the orthogonal estimators of the three pieces and combine them. What does the delta method say about the standard error? ## Reading The [DML paper](https://arxiv.org/abs/1608.00060) and [Kennedy review](https://arxiv.org/html/2203.06469v2) provide the formal conditions behind the constructions. - Hines, O., Dukes, O., Diaz-Ordaz, K., and Vansteelandt, S. (2022). Demystifying statistical learning based on efficient influence functions. *The American Statistician*, 76(3), 292 to 304. An introductory derivation-focused treatment. - Kennedy, E. H. (2022). Semiparametric doubly robust targeted double machine learning: a review. arXiv:2203.06469. Influence-function construction, remainder decomposition, and efficiency conditions. - Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. *The Econometrics Journal*, 21(1), C1 to C68. Section 2 for Neyman orthogonality and cross-fitting. - Vansteelandt, S. and Dukes, O. (2022). Assumption-lean inference for generalised linear model parameters. *Journal of the Royal Statistical Society B*, 84(3), 657 to 685. Asks what a regression coefficient estimates when the model is wrong, defines an estimand that does not depend on the model being right, and builds the orthogonal estimator for it. The target remains defined independently of the correctness of the working regression model. - Tsiatis, A. A. (2006). *Semiparametric Theory and Missing Data*. Springer. Chapters 3 and 4 for the geometry, when you want the full version. <span id="continue-through-the-collection"></span>