# Treatment effects via posterior prediction

**URL:** <https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346>\
**Category:** General\
**Created:** [April 3, 2019, 3:47pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346 "2019-04-03T15:47:20Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [April 3, 2019, 3:47pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/1 "2019-04-03T15:47:20Z")

</div>

These questions came out of conversations I had today with @mariel and @ignacio.

In the stan\_lm [vignette](https://cran.r-project.org/web/packages/rstanarm/vignettes/lm.html), @jonah and @bgoodri suggest investigating the average treatment effect using `posterior_predict`.

```
clouds_cf <- clouds
clouds_cf$seeding[] <- "yes"
y1_rep <- posterior_predict(post, newdata = clouds_cf)
clouds_cf$seeding[] <- "no"
y0_rep <- posterior_predict(post, newdata = clouds_cf)
qplot(x = c(y1_rep - y0_rep), geom = "histogram", xlab = "Estimated ATE")

```

From what I understand (please correct me if this is wrong), `posterior_predict` is not just calculating \hat y\_i = X\_i\hat\beta (the expected value of y); it is drawing a new \tilde y\_i from the posterior for each MCMC iteration. This incorporates the residual variability into the prediction.

First off, is there a problem with the above calculation, in that it assumes the draws `y1_rep` and `y0_rep` are independent? I would think you want to make this calculation in a way that takes into account the fact that these are the same people (really the same observation) in each case. Another way to look at it, is that I would expect y\_1 - y\_0 to be smaller if you are comparing those observations on the same person, as compared to two different people that happen to have the same covariates. I don’t think your suggested procedure accounts for this.

Second, I just want to point out that this is not the “standard” way of calculating an ATE, which is done by most software packages that I’ve come across. Most of the time, \hat y\_i = X\_i\hat\beta is used as the prediction, and the residual variability is ignored. When this is done, the issue I raise above is moot. I’m not saying that your approach is wrong, just pointing out that it is not the standard.

I’m curious about your thoughts to these issues. Thanks!

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [April 3, 2019, 3:55pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/2 "2019-04-03T15:55:50Z")

</div>

A third issue with the code in the vignette that was just pointed out by @mariel:

You calculate the posterior of y1-y0 for all individuals observations in the dataset. But then you just plot all the individual treatment effects (by vectorizing the resulting matrix). For the AVERAGE treatment effect, you should be averaging across observations (by taking `rowMeans(y1_rep - y0_rep)`). That totally changes the interpretation of what is being plotted. Do you agree?

---

<div class="post-metadata">

**Author:** ![tjmahr](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tjmahr/32/357_2.png) [@tjmahr](https://discourse.mc-stan.org/u/tjmahr)\
**Post date:** [April 3, 2019, 7:02pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/3 "2019-04-03T19:02:07Z")

</div>

> [@jgellar](#):
>
> Most of the time, ^yi=Xi^β\hat y\_i = X\_i\hat\beta is used as the prediction, and the residual variability is ignored. When this is done, the issue I raise above is moot.

Hmmm, but this is Bayes so we get to incorporate our uncertainty about the residual error into the predictions. rstanarm does have posterior\_linpred() for this conventional kind of prediction but the help pages advise against it. Maybe the phrase “average” treatment effect is a misnomer? Maybe “posterior-predictive treatment effect” is more apt.

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [April 3, 2019, 7:08pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/4 "2019-04-03T19:08:27Z")

</div>

It seems to me that you’re talking about two very different things in that paragraph, which correspond to points 2 and 3 that I raised in my first two posts.

1. You’re right, we are able to incorporate the uncertainty about the residual into our predictions. However, this isn’t limited to Bayesian analysis. Frequentists can do this as well. See, for example, the difference between “confidence intervals” and “prediction intervals”. Keep in mind that I don’t necessarily disagree with including the residual in the prediction (I still have to think about this more), but I do disagree with predicting y\_1 and y\_2 independently (unless @jonah or @bgoodri can convince me otherwise).

2. What that code displays is the distribution of individual-specific treatment effects. What we are typically interested in (for causal inference) is the average treatment effect. So yeah, if you call what is plotted in that plot an “average treatment effect”, that’s a misnomer. But the average is the thing we’re interested in, so I’d say it’s a “miscalculation” of what we want.

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [April 7, 2019, 4:09pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/5 "2019-04-07T16:09:55Z")

</div>

After thinking about this a couple days, there are two gripes I have with the way treatment effects are calculated in the vignette above.

1. My first gripe is pretty straightforward. You are not calculating and plotting the average treatment effect. You are plotting the distribution of individual treatment effects. That inflates the variance by a factor of \sqrt n.
2. The second gripe is more conceptual. The central issue is how correlated an observation is with its counterfactual. Until we get a time machine, that’s not something that we can really measure or know for sure. However, the standard in causal inference is to assume that they have a correlation of one. In other words, we assume the exact same residual for \tilde y\_{i1} and \tilde y\_{i0}. Thus, when you calculate \tilde y\_{i1} - \tilde y\_{i0}, the residuals cancel out, so you get the same thing as you would if you just did \hat y\_{i1} - \hat y\_{i0} (using the expected values). Whether that is appropriate or not could be debated, but it’s probably better than assuming the correlation is 0, which is what your code is doing. And this really isn’t a small difference - this is a **big** deal for estimating causal effects. In many cases the residual variance will be much larger than the variance of the impact of interest.

@shira, I’m wondering if you have any thoughts?

---

<div class="post-metadata">

**Author:** ![jgellar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgellar/32/2917_2.png) [@jgellar](https://discourse.mc-stan.org/u/jgellar)\
**Post date:** [April 12, 2019, 1:21am UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/6 "2019-04-12T01:21:55Z")

</div>

Just bumping this back up the list. I would very much like to hear a response from @jonah or @bgoodri, because I am now convinced that their approach for estimating treatment effects is flawed.

RE: @tjmahr’s point that “average treatment effect” may be a misnomer… the average treatment effect is a causal inference term - this is what we are trying to estimate. ATE = E[Y(1) - Y(0)]. You can’t blame the estimand if the estimator isn’t correct.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [April 12, 2019, 1:30am UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/7 "2019-04-12T01:30:48Z")

</div>

What we have in that vignette is wrong if you assume that the observed outcomes have measurement error that makes them differ from the potential outcome functions and that the measurement error is the same irrespective of which potential outcome function is realized. Under those assumptions, which I think a lot of people are comfortable making, then subtracting one from the other would cancel out the measurement error, leaving the difference in the potential outcome functions. The `posterior_predict` function is basically assuming independent measurement errors on top of the two potential outcome functions, in which case the subtraction does not cancel them.

---

<div class="post-metadata">

**Author:** ![pauper42](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/pauper42/32/4593_2.png) [@pauper42](https://discourse.mc-stan.org/u/pauper42)\
**Post date:** [April 18, 2019, 2:53pm UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/8 "2019-04-18T14:53:40Z")

</div>

@jgellar and @bgoodri , I have posted a related question here: [Help replicating example from Imbens and Rubin's Causal Inference book](https://discourse.mc-stan.org/t/help-replicating-example-from-imbens-and-rubins-causal-inference-book/8521?u=pauper42). I would love to hear your comments.

This approach to causal inference was first proposed by Donald Rubin (1978) and it’s explained quite well in the book [Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction](https://www.amazon.com/Causal-Inference-Statistics-Biomedical-Sciences/dp/0521885884/ref=sr_1_1?keywords=causal+inference+rubin&qid=1555598619&s=gateway&sr=8-1) (Chapter 8) where they call it _model-based_ imputation.

They do talk about the correlation between Y\_i(0) and Y\_i(1) and suggest that a conservative approach would have them perfectly correlated (p.175)

\begin{pmatrix} Y\_i(0) \\ Y\_i(1) \end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix} \mu\_0 \\ \mu\_1\end{pmatrix}, \begin{pmatrix} \sigma & \sigma\\ \sigma & \sigma\end{pmatrix}\right).

As you point out that correlation is not identified without a time-machine since we never observe both potential outcomes for the same individual, but that’s not a big deal for a Bayesian and it ends up reflected in a wider posterior distribution of the average treatment effects, as you also point out.

Also, about using rowMeans, I think you’re correct.

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [April 25, 2019, 3:19am UTC](https://discourse.mc-stan.org/t/treatment-effects-via-posterior-prediction/8346/9 "2019-04-25T03:19:45Z")

</div>

I totally missed this conversation, but I will contribute late. I struggled with this a lot while trying to decide whether I agreed with the claim in [this](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5790647/pdf/nihms934348.pdf) paper that posterior predictive draws could sometimes be replaced with their expectations. (I ultimately chose to make the assumption of independent errors in [my own research](https://arxiv.org/abs/1902.10613).)

Gripe #1 is a plot labelling issue, and I agree with @jgellar that you need to take the average at some point to call it an average treatment effect. I also agree that any definition of ATE (on the difference scale) that is not mathematically equivalent to \mathbb{E}\left[Y(1) - Y(0) \right] is plain wrong.

But isn’t gripe #2 a general gripe with causal inference? The relationship between residuals for Y\_i(0) and Y\_i(1) is unknown and unknowable. I’m not sure we should naturally assume the residuals have a correlation of 1; it’s just a conceptually simple way to explain how people have used parametric models for causal inference in the past.

Basically, I don’t see how you can get around making some kind of assumption here. At one end of the spectrum, we might make the strong and untestable assumption of completely independent errors. (Pearl’s non-parametric structural equation models with independent errors makes the _same_ independence assumption as the `posterior_predict` approach here, so I wouldn’t say there is a uniformly accepted “standard in causal inference.”) Or you can make a different strong and untestable assumption of perfectly correlated errors, in which case I think you should advertise that you are almost certainly being too conservative.