# Double Machine Learning: What kind of "naive" estimator do the authors use to get such a bias?

I am reading the double-machine learning paper by Chernozhukov et al. (2018), in Example 1.1. the authors consider the partially linear model:

$$Y = D \theta_0 + g_0(X) + U, E[U|X, D] = 0\\ D = m_0(X) + V, E[V|X] = 0$$

and then they're saying:

"A naive approach to estimation of $$\theta_0$$ using ML methods would be, for example, to construct a sophisticated ML estimator $$D\hat{\theta}_0 + \hat{g}_0(X)$$ for learning the regression function $$D \theta_0 + g(X)$$."

Then they simulate this approach and provide the following picture of the empirical distribution of $$\hat{\theta}_0 - \theta_0$$ based on simulations, illustrating that $$\hat{\theta}_0$$ is biased.

Does anybody have an idea of what kind of "sophisticated ML estimator" they used to get such results?

More generally, how would you estimate $$\theta_0$$ and $$g(X)$$ using standard ML methods?

What first came to my mind was to start with some guess for $$\theta_0$$, call it $$\tilde{\theta}$$ then run LightGBM regressor of $$Y - D\tilde{\theta}$$ on $$X$$, then compute a prediction $$\hat{g}_1(X)$$, compute $$\hat{Y}:=Y - \hat{g}_1(X)$$ and run a simple OLS of $$\hat{Y}$$ on $$D$$ getting a new value for $$\tilde{\theta}$$, then iterate until convergence.

For simulations, I took the same parameters $$n=500, p=20$$, $$X \sim U[0, 1]^p$$. For functions I took $$g_0(X) = 0.2+\sin(2 \pi X_1) + \frac{1}{1+e^{10*X_2}} + 0.4\max\{X_3, 0.3\}$$, $$D \sim Bern(\min\{\max\{0.3*X_1-X_2^2+X_4/3, 0.3\}, 1\})$$ (just random non-linear/non-smooth functions that came to my mind). Then I ran the above algorithm with LGB with default parameters, and got the following results (this is also the distribution of $$\hat{\theta}_0 - \theta$$ based on 700 Monte Carlo simulations. This doesn't at all look as biased as in the paper. I am sure what I was doing must be not the smartest thing to do, and there should be better standard methods, but I am a bit confused by the fact that even this simple method that didn't require any involved cross-fitting/orthogonalization/etc., as well as any tuning, worked not as bad as the paper says it should.

I think I figured the answer to my own question. Apparently, what they did was estimating $$E[Y|D, X]$$ (i.e. fitting a model using $$(D, X)$$ as features), and then computing the estimator of $$\theta_0$$ as $$\hat{\mathbb{E}}[Y|X, D=1] - \hat{\mathbb{E}}[Y|X, D=0]$$ and average the above quantity over $$X$$. Simulations (below) gave very similar (in terms of biasedness) distribution to that in the paper.
While a much better method of estimating $$\theta_0$$ was proposed by Kunzel et al. (2019) that they called "X-learner". The idea is that you use a subsample $$D=1$$ to estimate $$\mu_1(X):=\mathbb{E}[Y|X, D=1]$$, a subsample $$D=0$$ to estimate $$\mu_0(X):=\mathbb{E}[Y|X, D=0]$$, then for units with $$D=1$$ you compute a predicted $$\hat{\mu}_0(X)$$, and symetrically for units with $$D=0$$ you compute a predicted $$\hat{\mu}_1(X)$$. Finally, $$\hat{\theta}_0 = \frac{1}{n}\sum_{i:D_i = 1} (Y_i - \hat{\mu}_0(X_i)) + \frac{1}{n}\sum_{i:D_i = 0} (\hat{\mu}_1(X_i) - Y_i)$$