Why Residualize Both Outcome and Treatment?
In a high-dimensional setting, the observed outcome and treatment can both depend strongly on many covariates. Estimating the treatment coefficient while controlling flexibly for those covariates is difficult: errors in an estimated adjustment function can otherwise feed directly into the coefficient estimate. Double machine learning (DML) addresses this problem by estimating two nuisance regressions, removing their fitted contributions, and using a score designed to be insensitive to small first-order errors in those regressions.
The target here is the coefficient on a scalar treatment in a partially linear model. The treatment may be an exposure, policy, or intervention indicator, but the model alone does not establish a causal interpretation. That interpretation requires an appropriate causal design and assumptions about confounding in addition to the regression conditions below.
The conditional-mean restriction on \(U\) implies \(E[U\mid X]=0\). The treatment residual \(V=D-m_0(X)\) is the part of treatment not predicted by \(X\). The condition \(E[V^2]>0\) is an identification condition: after adjusting for \(X\), some variation in treatment must remain. If \(V\) is always zero, the treatment is completely determined by the covariates and its separate coefficient cannot be identified from this model.
The Orthogonal Score and Its Estimator
For candidate nuisance functions \(l\) and \(m\), define the score $$ \psi(W;\theta,l,m)=(D-m(X))\{Y-l(X)-\theta(D-m(X))\}. $$ At the true nuisance functions and \(\theta=\theta_0\), the expression in braces is \(U\). Thus the true score is \(VU\), whose expectation is zero. The score uses the outcome residual \(Y-l(X)\) and the treatment residual \(D-m(X)\), rather than fitting the target coefficient to the unadjusted outcome and treatment.
The population moment identifies \(\theta_0\). In fact, at the true nuisances, $$ E[\psi(W;\theta,l_0,m_0)] =E[V\{U-(\theta-\theta_0)V\}] =-(\theta-\theta_0)E[V^2]. $$ The first term vanishes because \(E[U\mid D,X]=0\), and the remaining coefficient is nonzero by \(E[V^2]>0\). Therefore the unique zero of this moment is \(\theta_0\).
This formula follows by setting the average score to zero. Writing \(\tilde D_i=D_i-\hat m_i\) and \(\tilde Y_i=Y_i-\hat l_i\), the equation is \(\sum_i\tilde D_i(\tilde Y_i-\theta\tilde D_i)=0\). Solving the linear equation in \(\theta\) gives the ratio above. It is a regression of the cross-fitted outcome residual on the cross-fitted treatment residual with no additional intercept.
Worked Example: Solving the Cross-Fitted Score
Suppose four hypothetical observations have the following observed values and cross-fitted nuisance predictions. These small data are for arithmetic illustration, not for making a reliable large-sample inference.
| Observation | \(D_i\) | \(\hat m_i\) | \(Y_i\) | \(\hat l_i\) | \(\tilde D_i\) | \(\tilde Y_i\) |
|---|---|---|---|---|---|---|
| 1 | 3 | 2 | 10 | 7 | 1 | 3 |
| 2 | 2 | 3 | 8 | 9 | -1 | -1 |
| 3 | 5 | 3 | 13 | 9 | 2 | 4 |
| 4 | 1 | 3 | 7 | 9 | -2 | -2 |
The numerator is $$ \sum_i\tilde D_i\tilde Y_i =(1)(3)+(-1)(-1)+(2)(4)+(-2)(-2)=16, $$ and the denominator is $$ \sum_i\tilde D_i^2=1^2+(-1)^2+2^2+(-2)^2=10. $$ Consequently, \(\hat\theta=16/10=1.600\). In the units of this hypothetical outcome and treatment, the fitted residualized outcome changes by 1.600 outcome units for a one-unit increase in residualized treatment. This calculation by itself does not establish a causal effect.
Why the Score Is Orthogonal
The defining advantage of this score is Neyman orthogonality: at the true nuisance functions, the expected score has zero first derivative in nuisance-function directions. Informally, a small nuisance estimation error does not create a first-order change in the target moment. That property helps prevent moderate estimation error in flexible regressions from dominating inference on \(\theta_0\).
To verify the claim and fix the signs, take arbitrary perturbation functions \(h_l(X)\) and \(h_m(X)\), and use the paths \(l_r=l_0+r h_l\) and \(m_r=m_0-r h_m\). At the true coefficient, \(D-m_r(X)=V+r h_m(X)\). Also, since \(Y-l_0(X)-\theta_0V=U\), the second factor of the score along this path is \(U-rh_l(X)-\theta_0r h_m(X)\). Thus $$ \psi(W;\theta_0,l_r,m_r) =(V+r h_m(X))\{U-rh_l(X)-\theta_0r h_m(X)\}. $$ Differentiating at \(r=0\) gives $$ \left.\frac{\partial}{\partial r}\psi(W;\theta_0,l_r,m_r)\right|_{r=0} =-Vh_l(X)+h_m(X)U-\theta_0Vh_m(X). $$ Each term has expectation zero: \(E[V\mid X]=0\), while \(E[U\mid D,X]=0\) implies \(E[U\mid X]=0\). Therefore the derivative of the expected score is zero in every such direction. The sign of the final term is negative for the stated path; its expectation is zero, but the sign matters for a correct derivative calculation.
Worked Example: Checking a Nuisance Perturbation
Consider the derivative above conditional on \(X\). The first term has conditional expectation $$ E[-Vh_l(X)\mid X]=-h_l(X)E[V\mid X]=0. $$ The second has conditional expectation \(h_m(X)E[U\mid X]=0\), and the third has conditional expectation \(-\theta_0h_m(X)E[V\mid X]=0\). Taking expectations over \(X\) therefore gives zero for the whole derivative. This check uses the model's conditional-mean assumptions, not a claim that \(U\) and \(V\) are independent.
Orthogonality is a first-order property, not a guarantee that nuisance errors have no effect at all. To see the remaining error, let \(\delta_l(X)=\hat l(X)-l_0(X)\) and \(\delta_m(X)=\hat m(X)-m_0(X)\) denote fixed nuisance errors. At \(\theta_0\), the score factors are \(V-\delta_m(X)\) and \(U-\delta_l(X)+\theta_0\delta_m(X)\). Conditional-mean assumptions eliminate the first-order terms, leaving $$ E[\psi(W;\theta_0,\hat l,\hat m)] =E[\delta_m(X)\delta_l(X)]-\theta_0E[\delta_m(X)^2]. $$ The leading population error is consequently second order in the nuisance errors. This also shows why both nuisance regressions matter and why treatment-regression error can contribute a squared-error term.
Cross-Fitting and Conditions for Inference
Cross-fitting separates nuisance fitting from score evaluation. Divide the observations into folds. For each fold, estimate \(l_0\) and \(m_0\) using the other folds, then predict both nuisance functions for observations in the held-out fold. Combine the held-out predictions across folds to form \(\hat l_i\) and \(\hat m_i\), and calculate the score using all \(n\) observations. This reduces dependence between an observation's score contribution and the nuisance fit used to construct it. The fold-specific regressions can use flexible prediction methods, but the inferential result still depends on their accuracy and on the model conditions.
A common sufficient rate requirement is that the nuisance prediction errors, measured in suitable \(L_2\) norms, are each \(o_p(n^{-1/4})\), together with appropriate moment and stability conditions. Then their product is \(o_p(n^{-1/2})\), and the squared treatment-regression error is also \(o_p(n^{-1/2})\). More generally, the relevant requirement is that the second-order remainder be negligible relative to \(n^{-1/2}\), alongside conditions ensuring a stable positive residual-treatment variance and a central limit theorem for the score. These rates are sufficient conditions, not a claim that every method or data set automatically satisfies them.
Under independent, identically distributed sampling, the stated conditional-mean assumptions, \(E[V^2]>0\), suitable finite moments, valid cross-fitting, and nuisance rates that make the remainder negligible, the estimator is asymptotically normal. Its influence function is \(V U/E[V^2]\): the score \(VU\) is divided by the derivative magnitude \(E[V^2]\) of the population moment with respect to \(\theta\).
Worked Example: Computing a Score-Based Standard Error
Use the four residuals in the first example, for which \(\hat\theta=1.600\). The estimated score contributions are $$ \hat\psi_i=\tilde D_i(\tilde Y_i-1.600\tilde D_i), $$ so the four values are \(1(3-1.6)=1.4\), \((-1)(-1+1.6)=-0.6\), \(2(4-3.2)=1.6\), and \((-2)(-2+3.2)=-2.4\). Their average is zero, as required by the fitted score equation.
The average squared score is $$ \hat\Omega=\frac{1.4^2+(-0.6)^2+1.6^2+(-2.4)^2}{4} =\frac{10.64}{4}=2.660, $$ and \(\hat Q=(1+1+4+4)/4=2.500\). Thus $$ \widehat{SE}(\hat\theta) =\sqrt{\frac{2.660}{4(2.500)^2}} =\sqrt{0.1064}=0.3262. $$ For illustration only, a normal-approximation 95% interval is \(1.600\pm1.96(0.3262)=[0.9606,2.2394]\). In a setting where the asymptotic approximation is justified, this procedure has 95% repeated-sampling coverage; it does not mean there is a 95% probability that the fixed coefficient lies in this particular interval. With only four observations, the interval should not be treated as reliable inference.
Common Pitfalls and Interpretation
- Using in-sample nuisance predictions without justification. Flexible fits can adapt to the same observations used in the score and make standard asymptotic arguments unreliable. Cross-fitting uses held-out predictions to reduce this problem.
- Assuming orthogonality removes all nuisance error. It removes the first-order effect at the truth. The product and squared-error remainder must still be small enough for the intended inference.
- Ignoring residual treatment variation. If \(\hat Q\) is near zero, the treatment is poorly distinguished from the covariates, and both the estimate and its standard error can be unstable.
- Calling the coefficient causal automatically. The partially linear regression and its score do not by themselves rule out unobserved confounding, reverse causation, or other identification failures.
- Using the wrong sign in a perturbation calculation. For \(m_r=m_0-rh_m\), the treatment residual is \(V+rh_m\), and the score's second factor is \(U-rh_l-\theta_0rh_m\) when \(l_r=l_0+rh_l\). The derivative is \(-Vh_l+h_mU-\theta_0Vh_m\).
The practical value of DML is that it combines flexible nuisance prediction with a target score whose first-order sensitivity to nuisance errors vanishes. Its coefficient estimate and standard error are meaningful only when the residualized treatment identifies the target, the nuisance estimates meet appropriate accuracy conditions, and the model and sampling assumptions support the desired interpretation.
Check Your Understanding
Answer each question using the score and assumptions developed above.
- Why is \(E[V^2]>0\) needed to identify \(\theta_0\) from the population score?
- Write the sample score equation and solve it for \(\hat\theta\) in terms of residualized outcome and treatment.
- For the paths \(l_r=l_0+rh_l\) and \(m_r=m_0-rh_m\), what is the derivative of the score at \(r=0\), and why does its expectation vanish?
- What second-order terms remain in the expected score when the nuisance estimates have errors \(\delta_l\) and \(\delta_m\)?
- Why does a robust score-based standard error not make a four-observation example reliable for inference?