# Nonparametric Causal Methods
```{r}
#| include: false
library(npcausal)
library(boot)
library(SuperLearner)
library(tmle)
library(tidyverse)
library(ranger)
library(ggplot2)
library(DoubleML)
library(mlr3)
library(mlr3learners)
```
> **Related reading**: For a hands-on TMLE simulation comparing
> super-learner libraries and assessing model misspecification, see the
> [companion blog chapter on TMLE](https://xiangao.github.io/blog_book/tmle.html)
> in *Topics on Econometrics and Causal Inference*.
## Why Do We Need Nonparametric Methods?
- Parametric models are often mis-specified
- Nonparametric models are more flexible; less assumptions
- What kind of nonparametrics are good?
Considerations:
- Plug-in (g-computation) estimators are easy.
- The problem is not the decomposition into $E[Y(0)]$ and $E[Y(1)]$ itself: with correctly specified models and regularity, a plug-in estimator can be efficient. The problem appears once the nuisance models are flexible and regularised. Then the error in the nuisance estimate enters the plug-in target *linearly*, and the plug-in has no term that cancels it. An influence-function estimator adds exactly such a term, and the error that survives is a product of two nuisance errors rather than one. @sec-plug-in-bias does the algebra.
- So they can be consistent, but the convergence rate inherits the (often slow) nuisance rate, which is especially bad in high dimensions.
- Ideally we remove that first-order sensitivity using an orthogonal influence function, and get $\sqrt{n}$ efficiency.
- Influence-function-based (one-step/AIPW/TMLE) estimators are shown to be efficient and orthogonal to first-order nuisance error.
## Influence Function Based Estimators
The semiparametric efficiency theory behind this section is set out in Bickel, Klaassen, Ritov and Wellner [-@bickel-klaassen-ritov-wellner-1993]; Kennedy [-@kennedy-2022-overview] gives a modern review aimed at exactly these estimators. Estimators based on efficient influence function are often $\sqrt n$ consistent, meaning they converge to the truth at the $1/\sqrt n$ rate. This means the estimation error not only goes to 0 as $n$ goes to infinity, but does so as fast as $1/ \sqrt n$ goes to 0. Why does this matter? Because we have limited sample size, when we are based on asymptotics, the faster the bias goes to 0 the smaller sample size we need.
If we have an estimator such that $$ \sqrt n (\hat \psi - \psi) = \sqrt n P_n(\phi) + o_P(1) \rightarrow N(0, var(\phi))$$ {#eq-nonparametric-1}
Then this estimator is $\sqrt n$ consistent and asymptotically normal; and $\phi$ is the influence function. The efficient influence function is the one with variance at the lower bound (similar to parametric case for Cramer-Rao lower bound).
Unfortunately there is no universal formula for efficient influence function for every problem. We need to derive it for each problem. For example, ATE has an IF based estimator, but ATT needs a different formula; and LATE has a different formula too, etc. Deriving them is hard!
Take the mean potential outcome under treatment, $\psi_1 = E \{ E(Y | X, A=1) \} = E[Y(1)]$ (the ATE is the difference $\psi_1 - \psi_0$ of the two arm means, and its EIF is the difference of the two arm-specific EIFs below). The efficient influence function for $\psi_1$ is $$\phi (Z; P) = \frac{A}{\pi (X)} \{ Y - \mu (X) \} + \mu(X) - \psi_1$$ {#eq-nonparametric-2} where $\pi (X) = P(A=1 | X)$ the propensity score model, $\mu (X) = E(Y |X, A=1)$, the outcome model.
For example, an IF-based bias-corrected estimator $$\hat \psi = P_n [ \frac{A}{\hat \pi (X)} \{ Y - \hat \mu (X) \} + \hat \mu(X) ]$$ {#eq-nonparametric-3} where $P_n$ means sample mean. This is the familiar AIPW estimator. It is known to be doubly robust.
Since we know the influence function, we can estimate the asymptotic variance. The uncentered summand $\frac{A}{\hat\pi(X)}\{Y-\hat\mu(X)\}+\hat\mu(X)$ differs from $\phi$ only by the constant $\psi_1$, so the two have the same variance. We can use any estimators for the nuisance parameters; the sample mean of the summand is then $\hat\psi_1$, its sample variance estimates $var(\phi)$, and the standard error of $\hat\psi_1$ is $sd(\hat\phi)/\sqrt{n}$ (differencing the two arms gives the ATE).
In this estimator, both $\hat \mu$ and $\hat \pi$ can be done using any estimator, such as machine learning. One caveat: $\sqrt n$ inference with flexible ML nuisances requires either Donsker-type complexity restrictions or sample splitting (cross-fitting), where the nuisance models are fit on one fold and evaluated on another. The manual estimators below use in-sample predictions for simplicity; the `npcausal` calls cross-fit with `nsplits=2`, and the DoubleML section below returns to sample splitting as a key ingredient.
### Why the Plug-In Is Off {#sec-plug-in-bias}
The bullets above claimed that the plug-in inherits the nuisance error at first
order and that the augmentation term removes it. That is algebra, not a
simulation finding, so it is worth doing before we run anything.
Take $\psi_1 = E[\mu(X)]$ with $\mu(x) = E(Y \mid X = x, A = 1)$. The plug-in
estimator is $\hat\psi_1^{\,\text{pi}} = P_n[\hat\mu(X)]$. Split its error into
a sampling part and a nuisance part:
$$
\hat\psi_1^{\,\text{pi}} - \psi_1
= \underbrace{(P_n - P)[\mu]}_{O_P(n^{-1/2})}
\; + \; \underbrace{P[\hat\mu - \mu]}_{\text{nuisance error, linear}}
\; + \; \underbrace{(P_n - P)[\hat\mu - \mu]}_{\text{smaller still}} .
$$ {#eq-nonparametric-plugin}
The middle term is $E[\hat\mu(X) - \mu(X)]$, and it is linear in the nuisance
error. If $\hat\mu$ converges at rate $n^{-\alpha}$, that term is of order
$n^{-\alpha}$, and a regularised learner in even a moderate number of
dimensions gives $\alpha < 1/2$. The nuisance error then dominates the sampling
error, the estimator is not $\sqrt n$ consistent, and nothing in the plug-in
cancels it. That is what it means to say the regularisation bias passes straight
into the average.
Now do the same for the AIPW estimator of @eq-nonparametric-3. Evaluate it under
the population distribution, holding the estimated nuisances fixed. Because
$E[A\{Y - \hat\mu(X)\} \mid X] = \pi(X)\{\mu(X) - \hat\mu(X)\}$,
$$
P\Bigl[\frac{A}{\hat\pi(X)}\bigl\{Y - \hat\mu(X)\bigr\} + \hat\mu(X)\Bigr]
= \psi_1 + E\Bigl[\frac{\pi(X) - \hat\pi(X)}{\hat\pi(X)}\,
\bigl\{\mu(X) - \hat\mu(X)\bigr\}\Bigr].
$$ {#eq-nonparametric-drift}
The bias is a *product* of the two nuisance errors, and that is the entire
difference between the two estimators. The term linear in $\hat\mu - \mu$ has
cancelled: the $\hat\mu$ that the outcome model carries in is subtracted back
out by the $\hat\mu$ sitting inside the augmentation, and what survives is
multiplied by the propensity-score error. By Cauchy-Schwarz the bias is bounded
by $\lVert \pi - \hat\pi \rVert \, \lVert \mu - \hat\mu \rVert$ up to the
overlap constant, so rates of $o(n^{-1/4})$ on each nuisance suffice to make it
$o(n^{-1/2})$ --- slower than $\sqrt n$ on each piece, fast enough for $\sqrt n$
on the target.
Two familiar facts fall out of @eq-nonparametric-drift directly. Double
robustness is immediate: if either factor is zero the whole bias is zero, which
is why one correct model is enough. And the $n^{-1/4}$ requirement is the reason
cross-fitting matters, since reusing the same observations to fit and to
evaluate $\hat\mu$ makes the product term behave worse than the rate
calculation assumes.
### Simulation Setup
The aim is to run one estimand through five estimators on a single data set ---
a plug-in outcome regression, AIPW computed by hand, AIPW from `npcausal`,
TMLE, and DoubleML --- and compare them against a truth we can compute exactly.
We simulate $n = 1000$ people with four covariates: $w_1, w_2 \sim
\text{Bernoulli}(0.5)$, $w_3 \sim U(0,4)$ and $w_4 \sim U(0,5)$, the last two
rounded to three decimals. Treatment is
$A \sim \text{Bernoulli}(\Lambda(-0.4 + 0.2w_2 + 0.15w_3 + 0.2w_4 + 0.15w_2w_4))$.
Both potential outcomes are Bernoulli with the same logit index,
$Y(d) \sim \text{Bernoulli}(\Lambda(-1 + d - 0.1w_1 + 0.3w_2 + 0.25w_3 + 0.2w_4 + 0.15w_2w_4))$,
and we observe $Y = AY(1) + (1-A)Y(0)$.
Treatment shifts the logit index by exactly 1, so the individual effect
$\Lambda(\text{lin} + 1) - \Lambda(\text{lin})$ varies across people even though
the shift is constant. Note also that $w_1$ enters the outcome but not the
treatment, and that both the propensity and the outcome carry a $w_2 w_4$
interaction. Neither nuisance function is linear in the covariates, which is
why flexible estimators are worth using here.
```{r npcausal2, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
set.seed(1)
n=1000
w1 <- rbinom(n, size=1, prob=0.5)
w2 <- rbinom(n, size=1, prob=0.5)
w3 <- round(runif(n, min=0, max=4), digits=3)
w4 <- round(runif(n, min=0, max=5), digits=3)
A <- rbinom(n, size=1, prob= plogis(-0.4 + 0.2*w2 + 0.15*w3 + 0.2*w4 + 0.15*w2*w4))
Y.1 <- rbinom(n, size=1, prob= plogis(-1 + 1 -0.1*w1 + 0.3*w2 + 0.25*w3 + 0.2*w4 + 0.15*w2*w4))
Y.0 <- rbinom(n, size=1, prob= plogis(-1 + 0 -0.1*w1 + 0.3*w2 + 0.25*w3 + 0.2*w4 + 0.15*w2*w4))
Y <- Y.1*A + Y.0*(1 - A)
W<-data.frame(cbind(w1,w2,w3,w4))
data <- data.frame(w1, w2, w3, w4, A, Y, Y.1, Y.0)
True_Psi <- mean(data$Y.1-data$Y.0);
cat(" True_Psi:", True_Psi)
```
The line above reports `mean(Y.1 - Y.0)` for this particular draw. That is a
legitimate quantity --- the *sample* average treatment effect --- but it is not a
fixed benchmark, and it is easy to mistake for one. The potential outcomes are
binary, so at $n = 1000$ this realized average has a sampling standard deviation
of about 0.020: two runs of the same DGP under different seeds can differ by 0.04
without anything being wrong. (The R and Julia editions of this chapter print
0.21 and 0.184 for exactly this reason.) Estimator standard errors in this
chapter are around 0.03, so a noisy target of that size would swamp the
comparison we are trying to make.
The *population* ATE is a definite number, and since the DGP is fully specified
we can integrate it out rather than simulate it:
```{r popate, warning=FALSE, cache=TRUE, message=FALSE}
# w1, w2 ~ Bernoulli(0.5); w3 ~ U(0,4); w4 ~ U(0,5). Average the true individual
# effect plogis(lin0 + 1) - plogis(lin0) over that covariate distribution.
pop_ate <- mean(sapply(list(c(0,0), c(0,1), c(1,0), c(1,1)), function(w)
integrate(function(w4) sapply(w4, function(u4)
integrate(function(w3) {
lin0 <- -1 - 0.1*w[1] + 0.3*w[2] + 0.25*w3 + 0.2*u4 + 0.15*w[2]*u4
plogis(lin0 + 1) - plogis(lin0)
}, 0, 4)$value / 4), 0, 5)$value / 5))
cat(sprintf("Realized-sample ATE (this draw): %.4f\n", True_Psi))
cat(sprintf("Exact population ATE: %.4f\n", pop_ate))
```
The realized draw is 0.2100 and the population value is 0.2028. Compare the
estimators below against 0.2028.
All of them need estimates of the outcome regression $\mu$ and the propensity
score $\pi$. We estimate both with the same super-learner library: MARS
(`SL.earth`), a generalised additive model, the lasso, logistic regression with
pairwise interactions, the sample mean, a random forest, and gradient boosting.
The super learner takes a cross-validated convex combination of these, so no
single functional form has to be right.
```{r npcausal3, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
SL.library <- c("SL.earth","SL.gam","SL.glmnet","SL.glm.interaction", "SL.mean","SL.ranger", "SL.xgboost")
```
### Plug-in Estimator from Outcome Regression
**For illustration only -- do not use in-sample ML nuisance predictions in applied work.** The manual `mu1hat`/`mu0hat`/`pi1hat` estimators below reuse the same data to fit and evaluate the outcome-regression and propensity models (no cross-fitting). With flexible machine-learning nuisance estimators this causes overfitting bias that invalidates the estimator's asymptotic normality and the validity of standard-error calculations based on it (Chernozhukov et al., 2018). Always use sample splitting / cross-fitting with ML nuisances -- as the `npcausal` calls below do (`nsplits=2`), and as DoubleML does later in this chapter.
We fit the outcome regression on all four covariates and treatment, then
predict each unit's outcome twice, once with $A$ set to 1 and once with $A$ set
to 0, and average the difference. This is the plug-in, or g-computation,
estimator: no propensity score anywhere.
```{r npcausal4, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
set.seed(1)
mufit <- SuperLearner(Y=data$Y, X=data %>% dplyr::select(w1,w2,w3,w4,A),
SL.library=SL.library, family="binomial")
# predict the PO Y^1
mu1hat<- predict(mufit, newdata=data %>% dplyr::select(w1,w2,w3,w4,A) %>% mutate(A=1))$pred
mu0hat<- predict(mufit, newdata=data %>% dplyr::select(w1,w2,w3,w4,A) %>% mutate(A=0))$pred
mean(mu1hat-mu0hat)
```
The plug-in gives 0.1825 against the true 0.2028, and comes with no standard
error. The gap is the linear nuisance term $E[\hat\mu - \mu]$ of
@eq-nonparametric-plugin showing up in one sample: the super learner's
regularisation bias passes into the average with nothing to cancel it. The
number illustrates the algebra; it is not the evidence for it.
### AIPW ATE
We now add the augmentation term. We fit the propensity score with the same
library, form the two arm summands
$\frac{A}{\hat\pi}(Y - \hat\mu_1) + \hat\mu_1$ and
$\frac{1-A}{1-\hat\pi}(Y - \hat\mu_0) + \hat\mu_0$, and average their
difference. The standard error is the sample standard deviation of that
difference divided by $\sqrt{n}$, the influence-function standard error. The
printed vector is the estimate followed by the lower and upper 95% bounds.
```{r npcausal6, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
set.seed(123)
pifit <- SuperLearner(Y=data$A, X=data %>% dplyr::select(w1,w2,w3,w4),
SL.library=SL.library, family="binomial")
pi1hat <- pifit$SL.predict
pi0hat <- 1- pi1hat
IF.1<-((data$A/pi1hat)*(data$Y-mu1hat)+mu1hat)
IF.0<-(((1-data$A)/pi0hat)*(data$Y-mu0hat)+mu0hat)
IF<-IF.1-IF.0
aipw.1<-mean(IF.1);aipw.0<-mean(IF.0)
aipw.manual<-aipw.1-aipw.0
ci.lb<-mean(IF)-qnorm(.975)*sd(IF)/sqrt(length(IF))
ci.ub<-mean(IF)+qnorm(.975)*sd(IF)/sqrt(length(IF))
res.manual.aipw<-c(aipw.manual,ci.lb, ci.ub)
res.manual.aipw
```
AIPW gives 0.1894 with a 95% interval of $[0.1292, 0.2495]$, which covers the
true 0.2028. The augmentation term moved the estimate from 0.1825 to 0.1894 and
supplied an interval the plug-in could not.
### AIPW ATE Using npcausal
The `npcausal` implementation follows Kennedy [-@kennedy-2022-overview]; the estimator itself is that of Robins, Rotnitzky and Zhao [-@robins-rotnitzky-zhao-1994].
The same estimator from a package, with one change that matters: `nsplits=2`
cross-fits the nuisances, so each unit's $\hat\mu$ and $\hat\pi$ come from a
model fitted without it. The output reports both arm means and their difference.
```{r npcausal7, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
library(npcausal)
set.seed(1)
aipw<- ate(y=Y, a=A, x=W, nsplits=2, sl.lib=SL.library)
```
Cross-fitting gives $E[Y(0)] = 0.5625$, $E[Y(1)] = 0.7568$, and a difference of
0.1943 with a standard error of 0.0316 and an interval of $[0.1324, 0.2563]$.
That is closer to the true 0.2028 than the in-sample version, and the standard
error is the one to trust.
### AIPW ATT
The ATT has its own efficient influence function. With
$p = P(A=1)$ and $\mu_0(X)=E[Y|X,A=0]$,
$$\phi_{ATT}(Z) = \frac{1}{p}\left\{ A(Y-\mu_0(X)) - (1-A)\frac{\pi(X)}{1-\pi(X)}(Y-\mu_0(X)) \right\} - \frac{A}{p}\,\psi_{ATT}.$$ {#eq-nonparametric-4}
Treated units contribute their outcome-model residual; control units
enter with the odds weight $\pi(X)/(1-\pi(X))$, which reweights them to
the covariate distribution of the treated. Note the correction term
uses the CONTROL units — an estimator computed on the treated
subsample alone would silently reduce to regression adjustment.
So we fit $\mu_0$ on the controls only, predict it for everyone, and build the
summand above from the propensity score already estimated. The printed vector
is again the estimate and the two 95% bounds.
```{r npcausal8b, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
# npcausal (attached above) exports its own SL.ranger, whose fit object is
# incompatible with SuperLearner's predict method; pin SuperLearner's wrapper
SL.ranger <- SuperLearner::SL.ranger
data2 <- data %>%
dplyr::mutate(pi0hat=as.numeric(pi0hat), pi1hat=as.numeric(pi1hat))
set.seed(1)
data0 <- data2 %>% dplyr::filter(A==0)
mufit <- SuperLearner(Y=data0$Y, X=data0 %>% dplyr::select(w1,w2,w3,w4),
SL.library=SL.library, family="binomial")
# predict mu0 for ALL units, treated and control
mu0hat_all <- as.numeric(predict(mufit, newdata=data2 %>% dplyr::select(w1,w2,w3,w4))$pred)
p1 <- mean(data2$A)
resid0 <- data2$Y - mu0hat_all
# uncentered IF-based ATT estimator
phi <- (data2$A*resid0 -
(1-data2$A)*(data2$pi1hat/data2$pi0hat)*resid0)/p1
aipw.manual.att <- mean(phi)
# centered EIF for the variance
phi_c <- phi - (data2$A/p1)*aipw.manual.att
ci.lb<-aipw.manual.att-qnorm(.975)*sd(phi_c)/sqrt(length(phi_c))
ci.ub<-aipw.manual.att+qnorm(.975)*sd(phi_c)/sqrt(length(phi_c))
res.manual.aipw.att<-c(aipw.manual.att,ci.lb, ci.ub)
res.manual.aipw.att
```
The ATT comes out at 0.1818 with an interval of $[0.1204, 0.2431]$.
### AIPW ATT Using npcausal
```{r npcausal10, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
library(npcausal)
set.seed(1)
aipw<- att(y=Y, a=A, x=W, nsplits=2, sl.lib=SL.library)
```
Cross-fitted, the ATT is 0.1730 with a standard error of 0.0529 and an interval
of $[0.0693, 0.2767]$.
The ATT is a different estimand from the ATE, and here they genuinely differ:
the treatment probability rises with $w_2$, $w_3$ and $w_4$, and the individual
effect varies with those same covariates, so averaging over the treated is not
the same as averaging over everyone. The two estimates are 0.1730 and 0.1943,
and with standard errors of 0.0529 and 0.0316 the data cannot separate them.
Note that the ATT standard error is two-thirds larger, because only the
treated units contribute to the average.
## TMLE
Targeted maximum likelihood estimation is due to van der Laan and Rubin [-@vanderlaan-rubin-2006]; van der Laan and Rose [-@vanderlaan-rose-2011] is the book-length treatment, and the source of the Super Learner used for the nuisance fits.
TMLE is a variant of IF based estimator. It is also a doubly robust estimator. See Mark van der Laan and his coauthors for the recent development of TMLE.
TMLE targets the same ATE with the same nuisance library, so we should expect
it to agree with `npcausal`.
```{r tmle1, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
library(tmle)
TMLE<- tmle(Y=data$Y,A=data$A,W=W, family="binomial", Q.SL.library=SL.library, g.SL.library=SL.library)
TMLE$estimates$ATE
```
TMLE gives 0.1893 with an interval of $[0.1288, 0.2498]$, against AIPW's 0.1894
and `npcausal`'s 0.1943. The three agree to well within a standard error, which
is what the theory predicts: they share an efficient influence function and
differ only in how they solve its estimating equation.
## DoubleML
Double/debiased machine learning is Chernozhukov et al. [-@chernozhukov-2018-dml]. The cross-fitting it relies on is what keeps the nuisance error from contaminating the influence function.
The R package DoubleML implements the double/debiased machine learning framework of Chernozhukov et al. (2018). It provides functionalities to estimate parameters in causal models based on machine learning methods. The double machine learning framework consists of three key ingredients: Neyman orthogonality, High-quality machine learning estimation and Sample splitting.
They consider a partially linear model:
$$ y_i = \theta d_i + g_0(x_i) + \eta_i $$ {#eq-nonparametric-5}
$$ d_i = m_0(x_i) + v_i $$ {#eq-nonparametric-6}
This model is quite general, except it does not allow interaction of $d$ and $x$; therefore no heterogeneous treatment effect across $x$. But "DoubleML" implements more than partially linear model, it actually allows for heterogeneous treatment effects, in models such as interactive regression model.
The basic idea of doubleML is to use any machine learning algorithm to estimate outcome equation ($l_0(X) = E(Y | X)$) and treatment equation ($m_0(X) = E(D | X)$). Then get the residuals, namely $\tilde Y=Y-\hat l_0(X)$ and $\tilde D = D - \hat m_0(X)$.
Then regress $\tilde Y$ on $\tilde D$. Based on FWL theorem, you get $\hat \theta$.
An important component here is to specify a Neyman-orthogonal score function. In the case of PLR, the "partialling-out" score (the DoubleML default, matching the residual-on-residual regression above) is
$$\psi (Z; \theta, \eta) = \big(Y- l(X) - \theta (D-m(X))\big)\big(D-m(X)\big) $$ {#eq-nonparametric-7}
The estimator $\hat \theta$ solves the equation that the sample mean of this score function is 0.
The estimator's variance comes from the usual sandwich formula built on the score, $\hat\sigma^2 = \hat J^{-1}\, \widehat{E[\psi^2]}\, \hat J^{-1}$ with $\hat J = -\frac{1}{n}\sum_i (D_i-\hat m(X_i))^2$.
We set up the data object first, naming the outcome, the treatment and the four
covariates:
```{r doubleml1, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
library(DoubleML)
dml_data = DoubleMLData$new(data,
y_col = "Y",
d_cols = "A",
x_cols = c("w1","w2","w3","w4"))
print(dml_data)
```
Both nuisance functions are estimated by a random forest with 500 trees, depth
at most 5, and a minimum node size of 2. DoubleML cross-fits them and solves
the partialling-out score.
```{r doubleml2, warning=FALSE, cache=TRUE, message=FALSE, echo=TRUE}
library(mlr3)
library(mlr3learners)
# suppress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")
learner = lrn("regr.ranger", num.trees=500, max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()
learner = lrn("regr.glmnet")
ml_g_sim = learner$clone()
ml_m_sim = learner$clone()
set.seed(123)
obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
obj_dml_plr$summary()
```
DoubleML returns 0.19257 with a standard error of 0.03136, in line with AIPW's
0.1894 and TMLE's 0.1893.
One caveat on what this number is. The partially linear model has no
treatment-covariate interaction, so its $\theta$ is not the ATE when effects are
heterogeneous, and in this DGP they are: the individual effect
$\Lambda(\text{lin}+1) - \Lambda(\text{lin})$ depends on the covariates. By the
Frisch-Waugh-Lovell argument behind the partialling-out score, $\theta$ is the
average of the conditional effects weighted by the conditional variance of
treatment, $\text{Var}(A \mid X)$, so people whose treatment is hardest to
predict count most. It equals the ATE only when the effect is constant. Here
the two are close, 0.193 against 0.203, but that is a property of this DGP and
not a general result.