# 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 real issue is that a flexible/regularized nuisance estimate carries a first-order *regularization bias* into the plug-in target, and the plug-in lacks the augmentation (correction) term that would make the estimator first-order insensitive to nuisance error.
- 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
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))$$
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$$ 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) ]$$ 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.
### Simulation Setup
```{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)
```
```{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
::: {.callout-warning}
**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.
:::
```{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)
```
### AIPW ATE
```{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 ATE Using npcausal
```{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)
```
### 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}.$$
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.
```{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
```
### 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)
```
## TMLE
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.
We should expect npcausal and tmle give similar answers.
```{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
```
## DoubleML
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 $$
$$ d_i = m_0(x_i) + v_i $$
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) $$
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$.
```{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)
```
```{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()
```