To fix the concepts I've introduced in the previous sections we're going to take a long, hard look at how linear regression might be used to estimate the effect of a treatment in an interventional study (either observational or randomized). First we'll motivate the definition of a treatment effect as an estimand. We'll then describe the usual maximum-likelihood version of linear regression as an estimator, which behaves well under the assumptions of the linear model. We'll then take this estimator outside of its "natural habitat" and show how it can be badly biased in generic observational studies. Lastly, we will show that this estimator is actually asymptotically unbiased in randomized trials, despite being motivated by a generally incorrect parametric model!

The point of this long exercise is to give you an idea of the danger of using "simple" methods without thinking critically about the relevant statistical model and estimand. It should also motivate our search for estimators that behave predictably and optimally across larger classes of problems.

Throughout this section (and throughout most of this book) the estimand we'll be interested in is the (statistical) average treatment effect (ATE). Given the joint distribution of variables $Y, A, X$ with $A \in \{0,1\}$ the average treatment effect is defined as $\psi = E[E[Y|X,A=1] - E[Y|X,A=0]]$. Note that the outer expectation here is over the full marginal distribution of $X$, not $X|A$.

We usually call $Y$ the outcome, $A$ the treatment, and $X$ the covariates. You can imagine, for example, a situation where $X_i$ represents a person's age, sex, and family history of disease prior to intervention, $A_i$ indicates whether that person received a new experimental treatment or a placebo treatment, and $Y_i$ indicates their condition after the treatment. Throughout this book I'll refer to this kind of data structure (where we have covariates, treatment, and an outcome) as interventional data.

The conditional expectation function $\mu_a(x)=E[Y|X=x,A=a]$ can thus be interpreted as the average outcome among those people with covariates $X=x$ who were treated with $A=a$. The difference $\mu_1(x) - \mu_0(x)$ describes how the much average outcome increases for people with $X=x$ if they are given the treatment $A=1$ instead of $A=0$. The ATE is just the average of this difference across the entire population (i.e. the average over the distribution of $X$).

You can see why this quantity is of immense interest. If we can assure ourselves that the treatment can be changed independently of any other unobserved variables that would also affect the outcome, then the ATE tells us exactly how much better off we'd be prescribing the new treatment relative to the placebo in this population. It directly answers the question: "how well does this treatment work?" in terms of the population and outcome in question.

Strictly speaking, what we're talking about here is called the statistical ATE. That distinguishes it from something called the causal ATE which will be part of our investigations in the next chapter. We'll see that the two are the same as long as we can satisfy certain assumptions.

Remember that, as with all estimands, the ATE is a function of the joint distribution of the data. You can see that if you write things more explicitly as $\psi(P) = E_{P_X}[ E_{P_{Y|A,X}}[Y|X,A=1] - E_{P_{Y|A,X}}[Y|X,A=0] ]$. Different distribution, possibly different ATE.

Now we've defined our general data structure and the estimand- all without having collected any real data at all. But let's say now we go and collect $n$ observations of $(Y_i, A_i, X_i)$. What do we do with that data in order to get an estimate of the ATE?

One very popular approach (but often unjustified, as we will see) is to use a linear regression estimator. The linear regression estimator is defined as the solution

Where we'll take $\hat\psi = \hat\gamma$ to be our estimate of the true, unknown ATE $\psi$. If we abbreviate $Z=[1,X^\top,A]^\top$ and $\theta = [\alpha, \beta^\top, \gamma]^\top$ we can abbreviate this as $\hat\theta = \argmin_\theta \sum (Y-Z^\top\theta)^2$. Taking a derivative with respect to $\theta$ and setting equal to zero gives the well-known analytical solution $\hat\theta = [\sum(Z_iZ_i^\top)]^{-1}\sum Z_iY_i$. The law of large numbers and limit arithmetic show that $\hat\theta \overset{p}{\to} \theta^* = E[ZZ^\top]^{-1}E[ZY]$ so as the number of samples we collect gets big, our estimates converge to this point $\theta^*$ that depends on what the data-generating distribution is, whatever it may be.

Unless you're already familiar with linear regression, this may look completely arbitrary. So why is this a common approach? Well, let's imagine that our statistical model consists only of parametric distributions that have the following form:

We'll call this statistical model $\mathcal M_{lin-norm}$ (linear-normal). In other words, we assume that $(A,X)$ has some unspecified joint distribution, but that $Y$ is linearly related to $A$ and $X$ plus some independent, normally distributed "error". Assuming linearity or gaussian errors is generally ridiculous in the real world, but let's go with it for now so we can motivate our estimator.

Since we have a parametric model, each distribution $P(Y,A,X) \in \mathcal M_{lin-norm}$ can be represented by the parameters $(\gamma, \beta, \alpha, \sigma)$ of $P(Y|A,X)$ in addition to the nonparametrically unspecified $P(A,X)$. We can therefore construct the (log-)likelihood function of the parameters and maximize it to find the parameters representing the conditional distribution most likely to have generated the data (leaving $U_X, U_A$ nonparametrically unspecified). When you work out the algebra, the optimization problem comes out looking exactly like what we have above. Computationally, it can be solved analytically or by gradient descent, but that's not our concern here. This is where the linear regression estimator comes from. Moreover, the theory of maximum likelihood estimation ensures that $\hat\gamma$ has a sampling distribution that itself converges to a normal centered at the true value $\gamma$.

However, this argument only justifies why $\hat\gamma$ is a reasonable estimator of the parameter $\gamma$ that exists in $\mathcal M_{lin-norm}$, but not why it is a reasonable estimator of the ATE $\psi$! The coefficient $\gamma$ is just some number that describes the distributions within the linear-normal model. Without further justification, it has no inherent meaning.

Within the linear-normal model, however, it turns out that $\gamma$ is in fact the same as the ATE. The proof is straightforward: taking the conditional expectation of $Y$ shows $\mu_a(x) = E[Y|A=a,X=x] = a\gamma + x\beta + \alpha$, so $\mu_1(X) - \mu_0(X) = \gamma$. Taking the expectation over $X$ changes nothing. Thus, for the linear-normal model, we have that $\psi = \gamma$.

That justifies linear regression when we know that the true probability distribution is in $\mathcal M_{lin-norm}$. However, in reality, that is never the case. If $P\notin \mathcal M_{lin-norm}$, then $\gamma$ doesn't even exist! The estimator $\hat\gamma$ can still be computed using arbitrary interventional data, but what does it mean? What can it mean? These are precisely the questions we will now turn to.

Let's keep the linear regression estimator, but change our statistical model to represent a generic observational study. In an observational study we follow a group of subjects who are assigned to their treatments by an unknown process and then we observe the outcomes later on. Based on that information alone, we can make no assumptions at all about the structure of the data-generating distribution:

The variables $U_X,U_A, U_Y$ could be anything with the right domains (e.g. $U_A \in \{0,1\}$) and could have any arbitrary joint distribution. The only exception is that here we've written $Y$ as the sum of its conditional mean $\mu_A(X) =E[Y|A,X]$ and the exogenous term $U_Y$, which forces $U_Y$ to be mean-zero. This is just a notational convenience, though: we could have just as easily defined some $\tilde U_Y$ to be arbitrary and then set $U_Y = \tilde U_Y - E[Y|A,X]$. We call this statistical model $\mathcal M_{obs}$ (observational). Note that $\mathcal M_{obs}$ is much "bigger" than $\mathcal M_{lin-norm}$ because it places fewer assumptions on what the joint distribution of $Y,A,X$ is.

Our question now is: does the linear regression estimator still work to estimate the ATE in this larger statistical model? To answer this we need to look at what $\hat\theta$ converges toward as we gather near-infinite amounts of data. We already claimed above that $\hat\theta \overset{P}{\to} \theta^* = E[ZZ^\top]^{-1}E[ZY]$. This result was arrived at without any mention of a specific statistical model, so it is true generically. Therefore the question comes down to establishing

Where the dots are stuff we don't care about (in the linear model, they correspond to estimates of $\alpha$ and $\beta$). If the central equality does not hold, then the estimate we get out of our linear regression estimator is likely to be near some value that isn't the ATE. If that were the case it would make no sense to use linear regression to try and estimate the ATE.

Unfortunately, that equality does not hold for $\mathcal M_{obs}$ in general. Obviously it holds inside $\mathcal M_{lin-norm}$, but the problem is that there are distributions outside of $\mathcal M_{lin-norm}$ where things go wrong. We can give one such example. Presume the true distribution is specified by this probability mass function:

Here $X$ has three categories and $A$ is binary with a joint distribution such that they are more likely to be equal in absolute value (i.e. $(X,A) =(1,1)$ or $(0,0)$ or $(-1,1)$) than not (i.e. $(1,0)$). The outcome $Y$ is 1 if any only if $X$ and $A$ are equal in absolute value. You might imagine this is a scenario where $X_i$ is some kind of U-shaped indicator of the type of treatment someone would benefit from. If $X_i=1$ or $-1$ then the best treatment is $A_i=1$ which results in the good outcome $Y_i=1$, but $A_i=0$ results in the bad outcome $Y_i=0$. On the other hand if $X=0$ then it is $A=0$ that gives the good outcome.

You can crank through the calculations (outlined in the toggle below) to find that in this scenario the ATE is $\psi = \frac{1}{6}$. However, you can also calculate that for this distribution the limit of the linear regression coefficients is $\theta^*= [\dots, \frac{101}{37}]^\top$. Clearly $\frac{1}{6} \ne \frac{101}{37}$ so we conclude that linear regression is not a consistent estimator of the ATE in a general observational model! Even if we had infinite data, we would still get the wrong answer.

Now let's work with a generic randomized trial. In a randomized trial we take a group of subjects, flip a coin to decide who gets what treatment, and then observe the outcomes later on. Based on that information alone, we can make no assumptions about the distribution of the covariates and how the affect the outcome (e.g. linearly), nor can we assume how the treatment factors into the outcome. What we can say for sure is that the treatment is independent of the measured covariates and of anything else that could possibly affect the outcome.

Here we are assured that $U_Y, A\perp X$ and $E[U_Y]=0$. We'll call this statistical model $\mathcal M_{RCT}$ (randomized controlled trial). Any distribution of the form above is a special case of our other model $\mathcal M_{obs}$, so we have that $\mathcal M_{RCT} \subset \mathcal M_{obs}$. The relationship with $\mathcal M_{lin-norm}$ is also clear: there is some overlap (when $U_Y$ is normal and $\mu_A(X)$ is linear in $A$ and $X$) but there are also distributions in $\mathcal M_{lin-norm}$ that aren't in $\mathcal M_{RCT}$ because the latter imposes a very specific form on the distribution of $A$. Therefore the two intersect but neither is a subset of the other.

We know the linear regression estimator is consistent in $\mathcal M_{lin-norm}$. We also know there are distributions in $\mathcal M_{obs}$ for which it is not. But are any of those distributions in $\mathcal M_{RCT}$?

To figure this out all we need to do is calculate the linear regression limit $\theta^* = E[ZZ^\top]^{-1}E[ZY]$ under the assumption that $U_Y, A\perp X$ and check that its last element is indeed the ATE $\psi = E[\mu_1(X) - \mu_0(X)]$.

We'll use one additional simplifying "assumption" here, which is that $E[X] = 0$. This isn't actually an assumption because we can always center the observed covariates (call them $\tilde X$) and in the limit the empirical centering is equivalent to de-meaning. This is useful because $E[XA]=0$ (one of the terms in $E[ZZ^\top]$) follows from the independence $A \perp X$ and $E[X]=0$. Moreover $E[XX^\top] = V[X]$. Thus

Now we turn our attention to $E[ZY]$. The terms in $E[ZY]$ are $[E[Y], E[X^\top Y], E[AY]]^\top$ but since we're only interested in the last element of the vector $E[ZZ^\top]^{-1}E[ZY]$ and we see there is a 0 in the middle position of the last row of $E[ZZ^\top]^{-1}$, it turns out that $E[X^\top Y]$ doesn't matter. Thus $\gamma^* = -2E[Y] +4E[AY]$. Now we use the fact that $E[U_Y]=0$ and also that $E[U_YA]=0$ by the independence $A \perp U_Y$.

And thus by multiplication and subtraction we arrive at $\gamma^*=E[\mu_1(X) - \mu_0(X)]$, which is exactly $\psi$, the ATE. We have therefore proven that the linear regression coefficient is a consistent estimator of ATE in a randomized trial.

The takeaway message from all of this is that you need to think critically about the estimand, the estimator and the statistical model together to be able to claim that your method "works". Linear regression is widely used in the observational study literature because introductory statistics classes fail to emphasize that the maximum likelihood results require totally unrealistic assumptions (linearity, normal errors). In fact, many curricula fail to even mention that generally what is of interest is a parameter that is defined outside of a parametric model (e.g. the ATE $\psi$ vs. some fictitious coefficient $\gamma$).

However, we also showed that, in some cases, linear regression still works to estimate the ATE despite the maximum likelihood assumptions not holding. This justifies the use of linear regression to estimate the ATE in randomized trials, which is also common practice. The point is, however, that this should be surprising rather than expected!

We haven't discussed this, but it turns out that the usual way of estimating the standard error (and thus confidence intervals and p-values) is not consistent outside of the linear-normal model, even in a randomized trial! That's why the sandwich estimator exists.

And what if we change the estimand? If $Y$ is binary, we might consider using logistic regression to estimate a marginal odds ratio instead of the ATE. It's common practice to exponentiate the coefficient from a logistic regression fit and call it an odds ratio representing some effect of treatment. But this, too, is a problem, because even in the relatively safe haven of a randomized trial, one can show that this estimator is consistent for a weird estimand called the conditional odds ratio and not the marginal odds ratio that is most often of interest! Things are even worse in survival settings: the hazard ratio is essentially meaningless without the assumption of proportional hazards!

The fact that inference is not taught from this perspective continues to create confusion in many fields. The big problem isn't about p-hacking, or frequentist vs. bayesian statistics, or this or that model. It's about not thinking clearly about the question and the assumptions.

Is it sometimes difficult to keep all of this in your head? Sure. But we promise that also gets much easier with practice! And it's certainly easier in the long run to understand the goal of your analysis and a framework in which to think than it is to try and memorize a flowchart of a hundred different methods that may all be wrong anyways. So think about the model, think about the estimand, and make sure your estimator makes sense in that context!