---
title: "Using machine learning for causal effect in observational study"
date: "2017-09-21"
---
## A simulation for an OLS model
In an observational study, we need to assume we have the functional form to get causal effect estimated correctly, in addition to the assumption of treatment being exogenous.
The simulation makes that concrete. We draw $n = 2000$ observations of three
correlated standard normals $(x, z, w)$, with covariances $0.8^2$ between $x$ and
$w$, $0.5^2$ between $x$ and $z$, and $0.6^2$ between $z$ and $w$. The DGP then
runs on the *logs of their squares*, not on the variables themselves: treatment
is $A \sim \text{Bernoulli}(1/(1+e^{g}))$ with
$g = \log(w^2 + 0.01) + N(0,1)$, the untreated outcome is
$y_0 = \log(x^2 + 0.01) + N(0,1)$, and the individual treatment effect is
$\tau_i = 2 + N(0,1)$, so the true ATE is exactly 2.
The point is what an analyst sees. Treatment depends on $\log w^2$ and the
outcome on $\log x^2$, but only $x$, $z$ and $w$ are observed. This is the
common situation: we have the right variables and the wrong functional form.
```{r}
#| cache: true
#| warning: false
library(MASS)
library(ggplot2)
library(dplyr)
library(tmle)
library(glmnet)
set.seed(366)
nobs <- 2000
xw <- .8
xz <- .5
zw <- .6
nrow <- 3
ncol <- 3
covarMat = matrix( c(1^2, xz^2, xw^2, xz^2, 1^2, zw^2, xw^2, zw^2, 1^2 ) , nrow=ncol , ncol=ncol )
mu <- rep(0,3)
rawvars <- mvrnorm(n=nobs, mu=mu, Sigma=covarMat)
df <- as_tibble(rawvars, .name_repair = "minimal")
names(df) <- c('x','z','w')
df <- df %>%
# A small constant inside each log() keeps log(v^2) from diverging to
# large negative values when v is close to 0 (e.g. log(w^2) alone ranges
# down to about -15 here), which otherwise pushes the propensity score
# P(A=1|W) toward 0/1 for those units and violates positivity.
mutate(log.x=log(x^2 + 0.01), log.z=log(z^2 + 0.01), log.w=log(w^2 + 0.01), z.sqr=z^2, w.sqr=w^2) %>%
mutate(g.var= log.w + rnorm(nobs)) %>%
mutate(A = rbinom(nobs, 1, 1/(1+exp((g.var))))) %>%
mutate(y0=rnorm(nobs) + log.x) %>%
mutate(tau.true = 2 + rnorm(nobs), y1=y0+tau.true, treat=A, y = treat*y1 + (1-treat)*y0)
lm1 <- lm(y ~ A + log.w + log.x , data=df)
summary(lm1)
lm2 <- lm(y ~ A , data=df)
summary(lm2)
lm3 <- lm(y ~ A + w, data=df)
summary(lm3)
lm4 <- lm(y ~ A + w + x, data=df)
summary(lm4)
```
In this example, treatment assignment process is determined by logged
w, and outcome is determined by logged x and treatment. However, what
we observe is w and x. In observational studies, this happens all the
time. In fact, this is an ideal situation, that we observe variables
that are determinants of outcome, although we are not sure about the
functional form that determines the outcome. However, this example
shows that unless we have observed exactly the factors themselves (in
this case logged x, w, which determines the DGP), we have biased
estimates of the true treatment effect.
The four estimates are 1.955, 1.470, 1.469 and 1.464 against a true 2.
Model 1 is the only model with reasonable estimate of treatment effect
(which is 2 in this case), at 1.955 with a standard error of 0.067 — and it is
also the only one given the correct functional form, `log.x` and `log.w`. Its
`log.x` coefficient of 1.001 recovers the DGP's coefficient of 1 almost exactly. Model 2 is a model with endogeneity: A is
correlated with the missing variable logged x. Model 3 and 4 we have
x and w, but not logged, therefore still biased.
The lesson here is the functional form does matter. However, we have
no way of knowing the functional form. What can we do here?
```{r}
#| cache: true
#| warning: false
# Algorithm set trimmed to 7 fast learners for render speed.
# A production analysis would also include SL.randomForest, SL.gbm, SL.gam.
Q.SL.library <- c("SL.glmnet","SL.glm","SL.glm.interaction", "SL.rpart","SL.bayesglm","SL.step","SL.mean")
g.SL.library <- c("SL.glmnet","SL.glm","SL.glm.interaction", "SL.rpart","SL.bayesglm","SL.step","SL.mean")
# tmle1 uses x and w; tmle2 additionally includes z. Note z is not in the true
# treatment or outcome DGP -- it is only correlated with x and w -- so adding it
# should not change the estimate much.
tmle1 <- tmle(Y = df$y, A = df$treat, W = df[,c('x','w')], g.SL.library = g.SL.library , Q.SL.library = Q.SL.library)
tmle1
tmle2 <- tmle(Y = df$y, A = df$treat, W = df[,c('x','w', 'z')], g.SL.library = g.SL.library , Q.SL.library = Q.SL.library)
tmle2
```
We use [Mark van der Laan's TMLE method](http://biostats.bepress.com/ucbbiostat/paper275/).
It uses [SuperLearner](http://biostats.bepress.com/ucbbiostat/paper222/) as
the initial estimator. It's an ensemble of multiple machine learning
algorithms. Therefore it does not need to assume the functional form
of the DGP. Even if we don't have the variables that determines the
DGP of outcome, if we observe some functions (even nonlinear
functions) of these variables, we can still get reasonable estimates
of the treatment effect.
In this example, we used multiple popular machine learning algorithms
in modeling both treatment assignment process and the outcome process.
The first TMLE model is with x and w (note not the logged x and w
which are in the true DGP), the second one with an additional variable
z.
TMLE gives an additive effect of 1.9637 with $x$ and $w$, and 1.9540 when $z$ is
added — against a true 2, and against 1.464 to 1.470 for the linear models given
the same untransformed variables. The 95% intervals, $[1.823, 2.105]$ and
$[1.813, 2.095]$, both cover 2. Adding $z$ changes almost nothing, which is what
should happen: $z$ is in neither the treatment nor the outcome equation and is
merely correlated with the variables that are.
It seems that TMLE results are less biased than the linear models with
x and w. That is putting it mildly: the linear models miss by about 0.53 and
TMLE by about 0.04, an order of magnitude better, purely from not having to
guess the functional form. It may not be better than the linear model with logged x and
w, but in empirical studies, we often cannot assume we have the
variables in the DGP, but only some proxy of the variables in the DGP.
I'll do more simulations to see whether TMLE does perform better in
the situation that we are not sure about the functional form. We
should expect that is the case.
The simple `tmle` package example used here is for a binary treatment. TMLE itself is not limited to binary treatments: there are TMLE estimators for continuous, multivalued, longitudinal, stochastic, and survival settings (e.g. the `ltmle` and `lmtp` packages), though which estimand and package you use differs by setting.
It's about time we embrace machine learning techniques into studies of
causal effect in observational studies.
---
<!-- see-also-footer -->
*Systematic treatment: [R](https://xiangao.github.io/causal_econometrics_guide/nonparametric.html) · [Julia](https://xiangao.github.io/causal_econometrics_julia/nonparametric.html).*