Double Machine Learning (DML)

Part of Statistical Methods — MORIE’s statistical-methods reference.

MORIE implements the Partially Linear Regression (PLR) and Interactive Regression (IRM) models from Chernozhukov et al. (2018) natively in both languages: cross-fitting, the nuisance learners and the orthogonal scores run on the package’s own cores (morie.fn.plr, morie.fn.irm; the same code in R, including morie_dml_clustered() for cluster-robust inference). Nothing external is installed or imported; the estimates were cross-validated against the DoubleML reference implementation when the native engines were written.

Partially Linear Regression

The PLR model posits

\[Y = \theta_0 D + g_0(X) + \varepsilon, \quad \mathbb{E}[\varepsilon \mid D, X] = 0\]
\[D = m_0(X) + \nu, \quad \mathbb{E}[\nu \mid X] = 0\]

where \(D\) is the treatment, \(X\) are confounders, and \(g_0, m_0\) are unknown nuisance functions estimated nonparametrically. The parameter of interest is \(\theta_0\) (ATE under PLR).

Cross-fitting

To avoid regularization bias from using the same sample for nuisance and target parameter estimation, DML uses K-fold cross-fitting:

  1. Split data into \(K\) folds.

  2. For each fold \(k\): fit nuisance models on the complement \(\mathcal{I}^c_k\), predict on \(\mathcal{I}_k\).

  3. Form residuals: \(\tilde{Y}_i = Y_i - \hat{g}_0(X_i)\), \(\tilde{D}_i = D_i - \hat{m}_0(X_i)\).

  4. Regress \(\tilde{Y}\) on \(\tilde{D}\) to obtain \(\hat{\theta}_0\).

This satisfies Neyman orthogonality and achieves \(\sqrt{n}\)-consistency under mild conditions on the nuisance estimators.

Neyman orthogonality

The score function \(\psi(W; \theta, \eta)\) satisfies

\[\partial_\eta \mathbb{E}[\psi(W; \theta_0, \eta_0)][\eta - \eta_0] = 0\]

ensuring that first-order errors in \(\hat{\eta}\) do not bias \(\hat{\theta}\).

MORIE implementation

Python entry points: morie.causal.estimate_double_ml() and morie.effects.estimate_plr() (both call morie.fn.plr).

Nuisance learners: ridge regressions for both \(g_0\) and \(m_0\), cross-fitted; the partialling-out score gives \(\hat\theta\). Default n_folds=5.

import random
from morie.causal import estimate_double_ml

# a frame, a CSV path or a dict of columns; here a synthetic confounded sample
r = random.Random(1)
x = [r.gauss(0, 1) for _ in range(200)]
t = [1 if xi + r.gauss(0, 1) > 0 else 0 for xi in x]
y = [0.5 * ti + xi + r.gauss(0, 1) for ti, xi in zip(t, x)]
result = estimate_double_ml({"y": y, "t": t, "x": x}, outcome="y", treatment="t", covariates=["x"])
print({k: round(v, 3) for k, v in result.items() if isinstance(v, float)})
# {'ate': 0.434, 'se': 0.171, 'ci_lower': 0.099, 'ci_upper': 0.769, 'pval': 0.011}

# on the CPADS PUMF (morie pull ocp21): outcome="heavy_drinking_30d", treatment="cannabis_any_use",
# covariates=["age_group", "gender", "province_region", "mental_health"]

(morie.estimate_ate() is the IPW-weighted estimator, not DML; see Causal Inference.)

Interactive Regression Model (IRM)

The IRM extends PLR to allow heterogeneous treatment effects. The model posits:

\[Y = g_0(D, X) + \varepsilon, \quad \mathbb{E}[\varepsilon \mid D, X] = 0\]
\[D = m_0(X) + \nu, \quad \mathbb{E}[\nu \mid X] = 0\]

Unlike PLR, the outcome function \(g_0(D, X)\) interacts treatment \(D\) with covariates \(X\), making the model suitable when CATE varies across individuals. The target estimand is the ATE:

\[\theta_0^{\text{IRM}} = \mathbb{E}[g_0(1, X) - g_0(0, X)]\]

with the doubly-robust score:

\[\psi^{\text{IRM}}_i = \hat{g}_0(1, X_i) - \hat{g}_0(0, X_i) + \frac{D_i \bigl(Y_i - \hat{g}_0(1, X_i)\bigr)}{\hat{m}_0(X_i)} - \frac{(1-D_i)\bigl(Y_i - \hat{g}_0(0, X_i)\bigr)}{1 - \hat{m}_0(X_i)}\]

Python entry point: morie.causal.estimate_irm() (native, morie.fn.irm): cross-fitted random-forest nuisances for \(g_0(d, X)\) and \(m_0(X)\), Neyman-orthogonal ATE score.

Partially Linear IV (PLIV) — LATE estimation

When treatment \(D\) is endogenous, a valid instrument \(Z\) (correlated with \(D\) but independent of \(\varepsilon\) given \(X\)) identifies the Local Average Treatment Effect (LATE):

\[\text{LATE} = \frac{\text{Cov}(Y, Z \mid X)}{\text{Cov}(D, Z \mid X)}\]

The Partially Linear IV model is:

\[Y = \theta_0 D + g_0(X) + \varepsilon\]
\[D = m_0(Z, X) + \nu\]

Cross-fitting proceeds as in PLR, with the additional first stage estimating \(\hat{m}_0(Z, X)\).

Python entry point: morie.effects.estimate_pliv()

Nuisance learner defaults

  • PLR, \(g_0\) and \(m_0\): native ridge regressions, cross-fitted.

  • IRM, \(g_0(d, X)\) and \(m_0(X)\): native random forests, cross-fitted.

  • PLIV: native first stage; falls back to 2SLS when the first stage is linear.

References

  • Chernozhukov V, Chetverikov D, Demirer M, Duflo E, Hansen C, Newey W, Robins J (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1–C68. https://doi.org/10.1111/ectj.12097

  • Bach P, Chernozhukov V, Kurz MS, Spindler M (2022). DoubleML — An object-oriented implementation of double machine learning in Python. JMLR, 23(53):1–6.

  • Imbens GW, Angrist JD (1994). Identification and estimation of local average treatment effects. Econometrica, 62(2):467–475. https://doi.org/10.2307/2951620