Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
63,620 characters · 13 sections · 66 citation commands
Demystifying and avoiding the OLS “weighting problem”: Unmodeled heterogeneity and straightforward solutions
\doparttoc \faketableofcontents
\affil[1]{Department of Statistics and Data Science, UCLA} \affil[2]{Department of Political Science, UCLA}
\doublespacing
When estimating the effect of a treatment on an outcome of interest by adjusting for covariates ($X$), researchers typically hope to interpret their result as a well-defined causal quantity, such as the average effect over strata such as the average treatment effect (ATE). Despite the introduction of many more flexible estimation procedures, regression adjustment remains the “workhorse approach” when adjusting for covariates to make causal claims in many disciplines aronow_does_2016.\footnote{In political science, for example, keele_adjusting_2010 found that regression adjustment was used in analyzing 95% of experiments in American Political Science Review, 95% in Journal of Politics, and 74% in American Journal of Political Science. The ubiquity of regression adjustment approaches in practice is echoed by authors in economics and other disciplines as well, such as angrist_mostly_2009,humphreys_bounds_2009, chattopadhyay2023implied. This is perhaps due to familiarity, ease of use, the efficiency of OLS and its suitably good performance in many contexts hoffmann_2023, green_analyzing_2011, kang_demystifying_2007, and well-established uncertainty estimation considerations and tools.} For a binary treatment ($D$) and outcome ($Y$), a simple bivariate regression of $Y$ on $D$ gives the difference in means estimate (the mean of $Y$ when $D=1$ minus the mean of $Y$ when $D=0$). However, this equivalence breaks down when covariates $X$ are included in the regression (as in $Y \sim D + X$), if the treatment effect may vary across values of $X$. Instead, the coefficient produced by this regression represents a weighted average of strata-specific difference in means estimates, where the weights are larger for strata with probability of treatment closer to 50% angrist_estimating_1998, angrist_mostly_2009. Such a weighted average is not typically of direct interest, and under severe enough heterogeneity in treatment effects, this can lead to widely incorrect substantive conclusions, as demonstrated below.
An ongoing literature regards these weights as a nuisance to be coped with or incorporated into analysis. Accordingly, authors have proposed diagnostics that index the potential for bias due to these weights or tools to aid interpretation given these weights aronow_does_2016, sloczynski_interpreting_2022, chattopadhyay_lmw_2023, chattopadhyay2023implied, kline2011, and have analyzed the behavior of proposed regression estimators in different settings in light of these weights (e.g. chernozhukov2013average).
In this article, we seek to demystify these weights, offering a simple derivation and conceptual clarification of why they arise, what they represent, and how to avoid them. We offer two main theoretical clarifications. First, under heterogeneous treatment effects (and varying probability of treatment), the linear regression of $Y$ on $D$ and $X$ will be misspecified (e.g. rubin1979, imbens_recent_2009, imbens2015). Re-expressing the regression coefficient as a weighted average of strata-wise effects simply offers a way to see how this misspecification causes the regression coefficient to diverge from the ATE. Such weights emerge, as we show, simply by relying on the Frisch-Waugh-Lowell (FWL) theorem frisch-waugh-1933,lovell-1963 to first construct the unit-level (not strata-level) weights. We can then group these weights to form a new expression of the strata-wise weights. Our expression does not rely on the usual assumption (as in angrist_estimating_1998 and later sources) that $D$ is linear in $X$, but reduces to that expression in the special case where $D$ is linear in $X$. This result corroborates the recent independent contribution of hahn2023properties.
Second, we clarify the assumption under which approaches can avoid these weights, which determines what (existing) estimators avoid this problem. Rather than attempting to diagnose, interpret, or otherwise cope with the undesired affect of these weights on interpretation, they can be avoided by any estimator that works under the separate linearity assumption, meaning each potential outcome $Y(d)$ is assumed to be linear in $X$, as opposed to the single linearity assumption that required $Y$ itself to be linear in $X$ and $D$. Equivalent assumptions have been stated in analyses by hahn2023properties, sloczynski_interpreting_2022, kline2011 and imbens_recent_2009. Fortunately, many well-known estimation approaches can be justified under the separate linearity framework including (i) regression imputation/g-computation/T-learner/multi-regression peters1941method, belson1956technique, robins_new_1986, kunzel_metalearners_2019, chattopadhyay2023implied; (ii) regression of $Y$ on $D$ and $X$ and their interaction lin_agnostic_2013, and (iii) balancing/calibration weights that achieve mean balance on $X$ in the treated and control groups (e.g. hainmueller_entropy_2012). Approaches (i) and (ii) are shown to be identical in the context of OLS models without regularization. Mean balancing weights (iii), when constructed to target the ATE, can be justified by the same assumption of separate linearity, but uses a different estimation approach and produces different point estimates. All three of these strategies produce unbiased ATE estimates under separate linearity, without the “weighting problem” suffered under the single linearity regression. These benefits do come at an efficiency cost, though it is low in settings with few covariates relative to sample size.
These considerations apply broadly to observational research in which an assumption of no unobserved confounding is relied upon. However, they can also apply to the analysis of even randomized experiments, concurring with lin_agnostic_2013. In particular, we illustrate these implications and recommendations for block-randomized experiments.
We consider settings where we are interested in estimating the average effect of binary treatment $D$ on outcome $Y$, while accounting for confounder $X$. The average treatment effect is defined as
where $Y_i(1)$ and $Y_i(0)$ denote the potential outcomes under treatment and control, respectively neyman1923applications, rubin1974estimating. In order to estimate this average treatment effect using observed data, we require an absence of unobserved confounders, concisely expressed as the conditional ignorability assumption,
For some discrete variables $X$ that satisfy conditional ignorability, and given consistency of the potential outcomes, the ATE can be identified according to:
where we suppress the $i$ subscripts for legibility and $DIM_x$ is the estimand for the difference in means in the stratum where $X=x$. In a given sample, this is approximated using the analog estimator,\footnote{This quantity could also properly be labeled as the (estimate of) the sample average treatment effect (SATE) rather than the ATE. However, we maintain the ATE notation for simplicity as its expectation over-samples is still the ATE, presuming the sample in question is a probability sample from the population of interest.}
The term $\hat{P}(X=x)$ gives the empirical proportion of units in each stratum in the sample, and can be thought of as the “natural” strata-wise weights, as they marginalize over the stratum according to the fraction of units falling in that stratum. It would be natural to form an analog estimator for this through sub-classification/stratification, simply computing the difference in means in each stratum of $X$ and compiling them per expression (ref). However, suppose we instead attempt to estimate the treatment effect by fitting a regression according to the model
While regression-based estimation of treatment effects has been widely used in practice across disciplines for decades, as angrist_estimating_1998 and many scholars since then have emphasized, $\hat{\tau}_{reg}$ do not in general yield $\hat{\tau}_{ATE}$, even when conditional ignorability holds. Rather, it can be understood as a version of Expression (ref) but with weights on each $\widehat{DIM}_x$ that differ from $\hat{P}(X=x)$. We now turn to a ground-up analysis of these strata-wise weights intended to demystify them and to recognize them as the natural consequence of misspecification of the linear model under heterogeneous effects.
While our plan is to (re-)derive strata-wise weights on the $\widehat{DIM}_x$ components, we start with “unit-level weights”, which arise simply because the regression coefficient is a weighted sum of $Y_i$ values. Specifically, let $\hat{d}(X_i) = \hat{p}(D_i = 1 | X_i) = \mathbf{X}_i \hat{\beta}$ where $\mathbf{X}_i =
$, and $\hat{\beta}$ is the estimated coefficients from a linear regression of $D$ on $X$. The unit-wise weighting formula that corresponds to the OLS estimate is defined according to
for some weights $w_i$. By the FWL theorem frisch-waugh-1933, lovell-1963,
which already takes the form of Equation (ref), with weights given by
Connection to propensity score. Though these individual-level weights are an intermediate step in our analysis, we give some consideration to their form here. First, as the denominator is a scalar, the behavior of these weights is revealed through the numerator, $D_i - \hat{d}(X_i)$. This is simply $1 - \hat{d}(X_i)$ for treated units and $\hat{d}(X_i)$ for control units. Compare this to inverse-propensity score weights (IPW), which are proportional to $1/\hat{d}(X_i)$ for treated units and $1/(1-\hat{d}(X_i))$ for control units. Like the IPW weights, these unit-level weights in Expression (ref) place greater weight on units when their treatment status is more “surprising” given the covariate values. Unlike IPW, the linear regression weights do not require constructing a ratio that has a denominator that can become close to zero, which can create explosive weights under IPW.\footnote{Equivalent unit-level weights are explored in depth by chattopadhyay2023implied and in Chattopadhyay2024Causation, where the authors usefully employ these weights to analogize how OLS implicitly compares (weighted) mean outcomes for the treated to (weighted) mean outcomes for the controls, sharpening the analogy to the difference-in-means estimator one might apply in a randomized experiment and discussing implications for qualities such as effective balance and sample size.} This similarity of the unit-level weights to the propensity score approach is explored in detail in kline2011.
Negative individual-level weights. The weights $w_i$ are constructed to be both positive and negative because Equation (ref) is structured as a single weighted sum/average. An isomorphic but useful way of signing the weights is to instead seek the weights that would appear in a “weighted difference in means” estimator,
where $\tilde{w}_i = D_i w_i + (1-D_i)(-w_i)$, which simply changes the sign for control units in order to accommodate the subtraction rather than the summation in the form of Expression (ref).
The weights $\tilde{w}$ signed as in Expression (ref) will often be positive, and any positive weight may be naturally interpreted as the relative contribution of each observation. However, they can be negative as well. Chattopadhyay2024Causation note that negative weights “translate into forming effect estimates that may lie outside the support or range of the observed data”, or that they “[leave] room for biases due to extrapolation from an incorrectly specified model.” borusyak2024negative also discuss negative weights, emphasizing that under randomization, these weights will all be non-negative, and so do not pose a problem for the resulting estimand.
We note that the meaning of negative individual weights is closely tied to estimates one would obtain from a linear probability model regressing the treatment indicators on $X$, even though that model is not run. Specifically, weights will be negative for a treated unit when $\hat{d}(X_i)>1$, and for control units when $\hat{d}(X_i)<0$. If a treated unit lies at an $X$ position “so unique to treated units” that a linearly-fitted model produces a predicted value greater than 1 there, it will have a negative weight. This can be thought of as indicating that a treated unit is “more like treated units than any control units”, so that it cannot be well compared to control units. The symmetric argument holds for control units with a negative value.
In this way, negative individual level weights do indicate extreme non-overlap, but in a very peculiar, model dependent sense and not in any non-parametric way that corresponds to broader notions of common support. These weights are also not a sure diagnostic for poor overlap or extreme counterfactuals: one could easily imagine a situation where all units above a certain value of $X$ are treated and all those below another value of $X$ are controls, indicating poor overlap. However, there can be many (or all) units outside the area of common support but with positive weights, because the fitted linear probability model for $D$ given $X$ does not exceed $1$ (for treated units) or drop below $0$ (for control units). Thus, we emphasize that negative weights can provide a warning that “some units look too much like treated units to have valid comparisons” (and likewise for controls), but a positive weight is not a guarantee of meaningful comparability in the sense of finding units of the opposite treatment value nearby in the covariate space.
With the individual-level weights in hand, we can now construct the strata-level weights of interest to our analysis. These weights are interesting because they allow us to see how the regression coefficient compares to ATE as a combination of average treatment effects by stratum as in Equation (ref). To do so we consider the case of discrete $X$ and the strata-wise weights, meaning those that can be expressed as,
where $w_x$ is the weight for the strata in which $X=x$ and $\hat{\tau}_x = \hat{\mathbb{E}}[Y_i(1) - Y_i(0) | X_i = x]$ is the conditional average treatment effect for subgroup $X_i = x$. The resulting $\hat{\tau}_{reg}$ correspond to the ATE only when $w_x = P(X=x)$. The best known expression for such strata-wise weights of OLS arises from angrist_estimating_1998,
This form shows that regression weights strata-wise $\widehat{DIM}_x$ components not according to $\hat{P}(X_i=x)$ as required to represent the ATE, but proportionally to $\hat{d}(X_i) (1 - \hat{d}(X_i))\hat{P}(X_i=x)$. Such a regression puts the most weight on strata in which the conditional variance of treatment status is largest, i.e. when $\hat{d}(X_i)$ is nearest to 50%. An intuition for this notes that OLS seeks to minimize squared error, and the opportunity to learn the most from strata with middling probabilities of treatment angrist_mostly_2009. Absent treatment effect heterogeneity, this is unproblematic. However, if there are high levels of effect heterogeneity in our data, then depending on how they correspond to strata of $X$ and the probability of treatment in those strata, these weights could move the regression coefficient far from the average treatment effect aronow_does_2016, humphreys_bounds_2009, angrist_mostly_2009,sloczynski_interpreting_2022. goldsmith2022contamination also discusses similar weights in the context of estimating effects for multiple treatments at once, where this leads to “contamination” of effect estimates.
These weights, however, are obtained under the additional assumption that the probability of treatment is linear in $X$, and not just modeled that way. angrist_mostly_2009 satisfy this by assuming the corresponding outcome regression can be saturated in $X$. This is often a reasonable assumption with discrete, low-dimensional $X$. Nevertheless we may wish to have a more general expression for the weights that does not rely on such an assumption, either for completeness or in service of generalizing to the case where $X$ is not discrete or is otherwise infeasible to include in a saturating form (i.e. $X$ with numerous multiple dimensions and levels). The issues that could arise when this type of assumption is not satistifed are highlighted by NBERw29709, who similarly analyze the weighted average representation of the two stage least squares estimate for the local average treatment effect presented by angrist_mostly_2009, and show that the equivalent linearity assumption that is invoked is often not satisfied in practice.
To obtain a more general expression for the weights, we simply organize the individual-level weights above into strata,
Rearranging these terms in the form $\hat{\tau}_{reg} = \sum_x w_x \tau_x$ is not in general possible, as the expected outcome under treatment and the expected outcome under control are weighted differently:
where $\pi_x = \hat{P}(D_i=1 | X_i = x)$ is the true, or correctly specified, probability of treatment in the sample for a given stratum. Thus the strata-wise weighting generally imposed by OLS for a given sample involves a combination of the true probability of treatment in the sample given $X$, and the probability of treatment given $X$ estimated using a linear model. This result is equivalent to the weighted representation independently developed by hahn2023properties, who also considers a special case for when the outcome model for the control group is linear. We can also formulate our result in terms of the discrepancy, $a_x$, between the true and linearly-approximated probability of treatment given $X$, so that $\hat{d}(X_i) = \pi_x + a_x$. The general strata-wise weights corresponding to OLS are then
Appendix (ref) gives additional details. Notice that this expression is infeasible to compute when $D$ is not linear in $X$, as $\pi_x$ is unknown. However, if $a_x = 0$ (i.e. $D$ is truly linear in $X$), this representation reduces to the variance-weighted representation of the regression coefficient given by angrist_estimating_1998. For many purposes then the angrist_estimating_1998 representation provides a clear and intuitive conception of the weights. Nevertheless, we found it necessary to have this more complete formulation in order to obtain the correct answer, as our simulations below show.
An OLS model regressing $Y$ on $D$ and $X$ alone would be correctly specified if the true conditional expectation function, $\mathbb{E}[Y | X, D]$ is linear in $D$ and $X$ as in $\mathbb{E}[Y | X, D] = \beta_0 + \beta_1 X + \beta_2D$. We call such an assumption the single linearity assumption, because it requires a single assumption about $Y$ being linear in some terms. Note that this allows for random variation in treatment effects not correlated with $X$, but it does not accommodate treatment effect variation that is correlated with $X$.
Accordingly, if the ATE rather than the regression coefficient is the target of inference, a natural solution is not to characterize or diagnose or bound the difference between the coefficient and the ATE, but rather to avoid the underlying misspecification problem. The strata-wise weights above merely describe the regression coefficient in other terms, and as such, elucidate the impact of misspecification for the difference between the coefficient and the ATE. Resolving this misspecification bias is a natural solution.
Here we consider the minimal change to regression practice that an investigator otherwise comfortable with a linear model could employ to avoid this misspecification concern and the resulting weighting behavior. We first define the separate linearity assumption, requiring that each of the potential outcomes is separately linear in $X$,
This assumption is slightly weaker than single linearity, which effectively forces $\alpha=\gamma$. While we label this assumption for easy reference and to distinguish it from single linearity, it is not intended to be novel, and further can be understood as the (often implicit) motivation for a number of longstanding approaches described next imbens2015, kline2011, imbens_recent_2009.
We turn next to a variety of known and straightforward estimation approaches that are in keeping with this assumption and that wholly avoid the awkwardness of “regression's weighting problem” rather than employing diagnostic or interpretational aids.
\paragraph{Interactions, imputation, g-computation, and stratification.}
Fortunately, a number of existing approaches are suggested by such an assumption. One alternative estimation approach, proposed by lin_agnostic_2013 initially to mitigate bias in covariate-adjusted estimates from randomized experiments, is to use a regression model that includes an interaction term for the treatment and confounder. goldsmith2022contamination similarly propose interacted models in the context of multi-valued treatments.
In the binary treatment setting we consider, the treatment effect estimate, $\hat{\tau}_{interact}$, is the estimated coefficient on the treatment variable $D_i$ in the regression of $Y_i$ on $D_i$, $X_i - \bar{X}$, and $D_i(X_i - \bar{X})$. The centering of $X$ in these terms is useful for interpretation: while the fitted models with and without centering $X$ are isomorphic, centering $X$ allows the coefficient on $D$ to represent the average effect “when $X$ is at its mean”. By linearity, this is also exactly the average marginal effect of $D$ on $Y$ taken across observations at their observed values of $X$.\footnote{This centering procedure and the resulting interpretation is complicated when one level of $X$ (or the intercept) must be dropped, as when $X$ represents block indicators or more generally group fixed effect. We describe this concern and propose a solution for this in Section (ref).}
Another straightforward approach to estimate the ATE under separate linearity is to run separate regressions of $Y$ on $X$ for the treatment group and the control group, and use the results from each regression to predict the unobserved outcomes in the other group. Then, using these predicted outcomes, we can estimate the individual treatment effect for each subject, and take the average of the estimated individual treatment effects to get the ATE. This is also known as g-computation, meaning that an estimate from g-computation will also be unbiased for the ATE in settings with high levels of heterogeneity in treatment effect and probability of treatment assignment robins_new_1986, snowden_implementation_2011. It has also been referred to even more recently as “multi-regression” Chattopadhyay2024Causation, who note earlier uses of this approach as far back as peters1941method and belson1956technique. It is also equivalent to the Oaxaca-Blinder decomposition oaxaca1973, blinder1973.
Specifically, let $\hat{\mu}_0(X_i)$ be the model fit to the control units and $\hat{\mu}_1(X_i)$ be the model fit to the treated units,
Because this explicitly puts weights of $1/N$ on every unit, the estimate does not suffer from the “weighting problem".\footnote{We note the equality of Expression (ref) to Expression (ref) above implies that one may either (i) compare each unit's observed outcome under the realized treatment status to the modeled outcome for that unit under the opposite, or (ii) for each unit, compare the modeled outcome under treatment to the modeled outcome under control, without using the observed outcome. This is a result of relying on OLS for each outcome model, since the average fitted value from a given model will be precisely equal to the average observed outcome over the same group. Such a property does not hold with estimators that cannot guarantee $\overline{\hat{\epsilon}}=0$.} Further, the Lin estimate and the regression imputation estimate are easily shown to be identical in this context. Specifically, the Lin estimate models $\mathbb{E}[Y|D,X]$ as $\beta_0 + \beta_1 D + \beta_2 X + \beta_3 DX$, and thus implies
Meanwhile, the imputation estimate is based on two regression models, one for the untreated group and one for the treated group:
Comparing this to the above equations, we can prove $\alpha_0=\beta_0$, $\alpha_1 =\beta_2$, $\gamma_0 = \beta_0 + \beta_1$, and $\gamma_1 = \beta_2 + \beta_3$ by showing the equivalency of the minimization problems in question (see Appendix (ref) for details). Using these equivalencies, we can show that the ATE estimate from regression imputation is equivalent to the estimate from the interacted regression,
Since these two estimation methods are identical under OLS, they can be used interchangeably. One can directly show the unbiasedness of either under assumptions of consistency and conditional ignorability,
Both are also identical to the stratification estimate when $X$ is discrete,
\paragraph{Mean balancing weights.} We can also connect the separate linearity assumption to a justification for calibration/balancing weights. Consider mean balancing weights for both the treated and control units, so that weighted average of covariates for each is equal to the overall unweighted average of covariates,
While many choices of weights can satisfy these constraints (when feasible), it is desirable to minimize some measure of their variation. In the case of maximum entropy weights hainmueller_entropy_2012, this is done by maximizing the entropy, $\sum_i w_i log(w_i/q_i)$. The key idea is that by achieving equal means on $X$ in the treated and control group (both equal to the full samples mean on $X$), then any linear function of $X$---which include $Y(1)$ and $Y(0)$ if separate linearity holds--- will also have equal means in these groups. A difference in means estimator using these weights is then unbiased for the ATE, without requiring any appeal to the relationship between this weighting procedure and propensity score modeling.\footnote{See Appendix (ref) for proof, with similar results targeting the ATT in hazlett2020kernel. See zhao2017entropy for a detailed analysis of identification and double-robustness of mean balancing weights.} kline2011 discusses how the regression imputation approach can be formulated as a weighting estimator that balances the means of covariates. An interesting feature of the mean balancing approach is that while it is justified by the same assumptions as regression imputation or the interactive regression above, their estimation strategies are different. Balancing weights do not require estimating the (nuisance) coefficients of any model. This may improve tolerance to misspecification (i.e. non-linear conditional expectations for each potential outcome) at the cost of variance, as borne out in simulations below.
To verify that these estimators behave as our analysis claims, we explore the performance of different estimation approaches under three different data generating processes, first with a discrete covariate $X$ and then with a continuous one.
The first three DGPs involve a binary treatment variable $D$, a discrete covariate $X$ in the range $[-3, 3]$, and an outcome $Y$ which depends on $D$ and $X$. For each simulation setting, the tables in Figure (ref) show the possible values of $X$, the corresponding probability of treatment, probability that $X=x$, and the average treatment effect $\tau_x$ for the subgroup where $X=x$. In all simulations, noise is added to the outcome to achieve an $R^2$ of 0.33 between the systematic (noiseless) portion of $Y$ and the final $Y$ with noise.
As expected, when there is a high level of heterogeneity in treatment effect and in probability of treatment, the regression adjustment estimate will have a substantial amount of bias. Figure (ref) shows the effect estimates when there is heterogeneity in both treatment effect and probability of treatment between subgroups. In the first two settings, the potential outcomes are linear in $X$. Here we see that the OLS estimate (reg) is heavily biased, and would lead investigators to conclude there is a statistically significant positive effect, though the true ATE is zero. Notably, slightly more extreme simulations would even make it possible for OLS to produce the incorrect sign for the treatment effect estimate. Meanwhile, regression imputation (impute), the Lin estimator (interact), mean balancing (meanbal), and matching (match) all successfully address this concern. We recommend the use of regression imputation or equivalently the Lin/interacted adjustment as a simple way to improve on conventional OLS estimation.
Naturally, in the case where the potential outcomes are not linear in the treatment and covariates, the estimates from the interacted model, imputation, and mean balancing are all biased for the ATE. This is expected, as the assumptions of both single and separate linearity are violated. Addressing such non-linearity requires non-linear estimators, such as matching. We also note that the variance-weighted estimate (using Expression (ref)) does not always reproduce the actual OLS estimate, as in the first simulation setting, where $P(D=1|X)$ is not linear in $X$. This is easily avoided by fully saturating the model in $X$, although such a manuever would not work for the continuous cases below.
The more general (but infeasible) weighting representation (Expression (ref)) reproduces the OLS estimate exactly regardless of the functional form of $P(D_i|X_i)$.
While these results show the expected variability in estimates under resampling from the given DGP, for an investigator working with one observed dataset, some form of estimated standard error is vitally important to inference. Table (ref) reports the average analytically estimated standard errors for DGP1 above, still with discrete $X$. Results are similar in other settings.
For stratification, we take a weighted sum of strata-specific Neyman variances. Here and below we write expressions in the plug-in/analog sample estimator form.
For the simple and interacted regression, we calculate the HC2 standard error of the coefficient on the treatment variable. For regression imputation, we calculate the standard error again in the Neyman style,
where $\hat{\Sigma}_1$ and $\hat{\Sigma}_0$ are the estimated variance-covariance matrices for the treatment model and the control model, respectively.
For mean balancing we show two types of analytical standard errors. First, we use the (HC2) standard error from the weighted regression of just $Y$ on $D$. Second, meanbal-adj uses the HC2 standard errors from a regression of $Y$ on $D$ and $X$, again with the estimated weights. Both methods produce identical point estimates (when perfect mean balance is achieved by the weights), but the analytical standard errors of the meanbal-adj approach benefit from partialing out the $X$. This is akin to how conventional OLS standard errors, under a fixed design, partial $X$ out of $Y$ so that the estimates are built on the conditional/residual variance of $Y$ rather than the total variance. For matching we use the Abadie-Imbens standard error abadie_large_2006.
We find, first, that the empirical standard errors are extremely similar across the methods relying on single or double linearity (reg, impute, interact, meanbal, meanbal-adj). The approaches not relying on single or separate linearity (stratify and \texttt{match}) show somewhat larger standard errors, as expected given their greater flexibility. Second, each method's average analytical standard is within 5% of the empirical standard error, with the exception of \texttt{meanbal}, with an average analytical standard error almost 40% larger than the empirical value. This is repaired, however, by using \texttt{meanbal-adj}.
Finally, some investigators may be concerned with efficiency in the sense of sampling variability in the estimate and its consequences for inference. This can be understood as a bias-variance tradeoff. The comparison between OLS (reg) and the Lin/interaction approach (interact) offers a simple starting point, since the later will add one additional parameter to the regression for every covariate dimension. Because the number of observations is large relative to the number of covariates, this has very little impact on the estimate's variability across resamples. For example in DGP 1, Table (ref) shows that the empirical SE is only about 5% larger for interact than for reg. The behavior of impute is of course identical. The empirical SE from meanbal and \texttt{meanbal\-adj} are similarly only 6% larger than from \texttt{reg}. If an investigator is primarily concerned with root mean square error (RMSE) of the estimates around the true value, RMSE values fall by nearly half for each of these methods relative to \texttt{reg}. That is, the small loss of efficiency in these settings (increasing variance) is more than made up for by reductions in bias as they factor into the RMSE.
It is important to recognize, however, the favorable nature of our simulation setting in this regard. If the number of covariates was large enough relative to the sample size, if treatment probability varied little by stratum of $X$, and/or if efficiency was a greater concern than bias or RMSE, then investigators might have cause to prefer OLS and adopting its weighted-ATE as the target estimand for their inferential purposes.
While we have considered discrete $X$ thus far, we do so for the sake of intuition regarding strata, but the lessons apply to settings with continuous $X$ as well. In the simulations shown in Figures (ref) and (ref), $X$ is randomly sampled from uniformly from $[-3, 3]$, and $Y(1) = Y(0) + 3X$. However, $P(D|X)$ takes on a different form in each setting. For the first continuous specification, $P(D \vert X)$ takes a logistic form. In the second, $P(D \vert X)$ increase linearly with $X$. For the third specification, $P(D \vert X)$ increases linearly in $X$, but $Y(1)$ is nonlinear in $X$.
In each of these settings, we see that the OLS estimate is biased for the true ATE. As before, the interaction, imputation, and mean balancing estimators perform well, except in Figure (ref) where even separate linearity fails.
These lessons may at first seem to be directed to researchers working with observational data such that conditioning on $X$ is a requirement. However, they also apply in settings where conditioning on $X$ is used to improve precision/ minimize finite sample differences from the expected result. Our advice essentially echoes that of lin_agnostic_2013 and a broad subsequent literature calling for a model that interacts treatment with the (centered) covariates. As noted above, this is identical to analogous imputation/g-computation/T-learner approaches when using OLS for the underlying models.
Consider first the cases where investigators have a randomized experiment, but adjust for covariates in the analysis to improve precision. When $D$ is fully randomized, the probability of treatment across values of $X$ will typically not vary greatly, except by chance. This implies that the bias due to regression's weighting behavior will typically be small. Nevertheless, any discrepancy between the ATE and coefficient on account of these weights can be avoided entirely.
Second, block randomization is a powerful tool often employed by experimentalists to reduce finite sample deviations from expectation. For the estimates to reflect the reduced variance this affords, it is standard to regress the outcome ($Y$) on the treatment ($D$) and block indicators ($B$) as fixed effect covariates. The weights implied by regression of $Y$ on $D$ and block indicators $B$ would be given by
where it is innocuous to assume that $P(D=1|B)$ is linear in the block indicators, as this regression would be fully saturated in $B$. If all blocks have the same probability of treatment, the weights are equal, so the regression coefficient will represent the difference in means per block, averaged over the size of each block, i.e. the ATE. This will often be the case. However, if the treatment probability varies by block, either by design or otherwise, then the resulting coefficient estimate would not in general be the ATE. Rather, regression will put higher weight on blocks with probabilities of treatment nearer to 50%.
As above, this is simply a consequence of misspecification generated by heterogeneous treatment effect estimates by block. Including interactions between the treatment and the block fixed effects would address this, under the separate linearity assumption,
One complication when applying the interacted approach here is that care must be taken regarding the interpretation due to the interaction. In typical usage, one block indicator will be omitted to avoid co-linearity with the intercept. If no centering/de-meaning is done on the block indicators, the coefficient estimate would represent the estimated effect (difference in means) in whichever block had its indicator omitted. The solution of centering covariates lin_agnostic_2013 as utilized above is now complicated by this omission. It is possible to omit the intercept, rather than the indicator for one block, and utilize the centering. However, a more general solution is to avoid any consideration of centering and omitting one level/the intercept, and works when dealing with one or multiple categorical variables. This is to simply compute the marginal effect estimate for each observation, in whatever block it is in, $\frac{\partial Y}{\partial D}\big|_{B=b}$. Averaging these across observations (giving equal weight to each observations) produces the average marginal effect (AME),
where $\widehat{\frac{\partial Y}{\partial D}}\big|_{B=b}$ is the estimated marginal effect in block $b$ where this individual unit is found. For example, in the regression
the estimated marginal effect in block $B=b$ is $\hat{\tau} + \alpha_p$, and so the average marginal effect of interest over the whole sample is simply
Using this approach to interpret the fitted regression with interactions always produces an estimate with the desired interpretation of “the estimated marginal effect, which can differ by block, averaged over blocks according to how many observations fall in each.” In the case of block indicators or other categorical variables, this approach avoids errors in relation to choices about what levels to omit in the regression, whether to omit the intercept, or having to center these variables. The variance for this estimator is
For comparison, we also consider two other approaches. One practice is to weight the blocked fixed effects regression by the (stabilized) inverse probability of treatment assignment for each block,
These weights are constructed so that within each block, the treated and control units make up equal weighted proportions. They therefore neutralize any differences in the $P(D=1|B)$ across blocks. If we perform a weighted regression of $Y$ on $D$ and the block dummy variables using these weights, the coefficient estimate for $D$ will be unbiased for the ATE.\footnote{A weighted difference in means (rather than regression including $B$) with these weights would produce the same point estimate, but the improvement in efficiency obtained under the block randomization design will not be fully reflected in the estimated standard error. This occurs because the block indicators will be orthogonal to treatment (once weighted) and so do not affect the coefficient estimate on treatment, but the omission of $B$ prevents the model from reducing the residuals that enter the standard error.}
Finally we also the stratification estimator (Expression (ref)) but with blocks as strata,
Figure (ref) compares estimates from the “plain” block fixed effects (block FE) regression to the interacted regression with the average marginal effect interpretation (block FE interact), the IPW-weighted regression (block FE IPW), and the block DIM average (mean block DIM) in a simulation setting where treatment probability varies by block and there is heterogeneity in treatment effect. As expected, the block fixed effects OLS model suffers from the weighting problem. All three alternative approaches produce unbiased and identical estimates.
We recommend the use of one of these alternative estimators to improve upon estimation of the ATE under block randomization.\footnote{These results are consistent with those demonstrated in \href{https://declaredesign.org/blog/posts/biased-fixed-effects.html}{https://declaredesign.org/blog/posts/biased-fixed-effects.html}, which use the DeclareDesign simulation approach blair2019declaring and conclude that the mean blockwise DiM/ stratification approach is suitable.}
Given the wide use of OLS regression of $Y$ on $D$ and $X$ in practice, and the potential gap between the coefficient and target quantities like the ATE, prior work has provided several forms of guidance to practitioners regarding this weighting problem. humphreys_bounds_2009 notes the conditions under which the OLS result will be nearer to ATT or ATC, and observes monotonicity conditions under which the coefficient estimate will fall between these two endpoints. aronow_does_2016 suggest characterizing the “effective sample” for which the regression coefficient would be an unbiased estimate of the ATE by applying the variance weights to the covariate distribution. Other papers have provided diagnostics for quantifying the severity of this bias. chattopadhyay2023implied characterize the implied unit-level weights of regression adjustment (as we do as an intermediate result), and propose diagnostics that use these weights to analyze covariate balance, extrapolation outside the support, variance of the estimator, and effective sample size. These ideas relate closely to those of aronow_does_2016 in that they compare the covariate distribution targeted by the estimator to that of the target population. sloczynski_interpreting_2022 derives a diagnostic related to variation in treatment probability,
where $\rho$ is the unconditional probability of treatment. The bias due to heterogeneity will be $\delta (\tau_{ATC} - \tau_{ATT})$ where $\tau_{ATC}$ and $\tau_{ATT}$ are the average treatment effect among treated and control, respectively. Thus $\delta$ captures one aspect of the potential for bias---variation in treatment probability---but reasoning about the severity of bias requires knowledge of the heterogeneity in treatment effects, which is unknown to the investigator.
However, investigators could instead consider different modeling choices that avoid this “weighting problem” altogether. We have sought to clarify the theoretical considerations required to do so, first by understanding the weighting problem as a symptom of misspecification when the treatment effect (and probability of treatment) vary with $X$. We provide a more general expression for the “weights of regression” by algebraically manipulating a representation of the OLS coefficient. The key fact is that the apparent “weights of regression” simply account for the way misspecification alters the coefficient away from the ATE. Second, viewed this way, a natural proposal for research practice is to rely instead on specifications that will not necessarily be violated by effect heterogeneity. To focus on minimal changes in practice, the conventional single linearity assumption (that $Y$ is linear in $D$ and $X$) can be relaxed at least to the separate linearity assumptions (each potential outcome is separately linear in $X$). This recapitulates equivalent assumptions employed in analyses by hahn2023properties, sloczynski_interpreting_2022, kline2011 and imbens_recent_2009. It also clarifies which estimators are suitable to avoid these weights, and whose properties we can consider by comparison. Fortunately, a number of existing, straightforward estimation approaches produce unbiased ATE estimates under this assumption. First, regression imputation (including g-estimation and the T-learner) with OLS as well as including the interaction of $X$ and $D$ (as in lin_agnostic_2013) are all suitable, and identical to each other in this setting. This longstanding approach dates back to at least peters1941method and belson1956technique, as noted by Chattopadhyay2024Causation who label it as “multi-regression”. In addition, “mean balancing” estimators applied to the ATE---those that choose weights to achieve the same mean of $X$ for the treated and control group as in the full sample---are also unbiased under separate linearity. These are not identical to the above as they avoid fitting any model, and our simulations suggest this property can partially mitigate bias when even the separate linearity assumption fails. These approaches can be useful not only in observational studies (so far as investigators can claim conditional ignorability), but also when using covariate adjustment after randomization, including the special case of analyzing block randomized experiments.
Avoiding these weights through these alternatives may be desirable in many cases, but is not without cost of limitation. The natural cost to first consider (relative to regression under single linearity) is the potential loss in efficiency. The interact/imputation/g-computation/T-learner approach adds an additional nuisance parameter per covariate. The circumstances investigated here---with just one covariate, a high degree of treatment effect heterogeneity, and large variation in treatment probability---clearly favor these approaches, where they eliminate bias and reduce RMSE by roughly half, while increasing the standard error by only about 5% in our settings. However, there may be settings in which the improved bias does not so assuredly dominate the efficiency cost. For example, where there are many covariates relative to the sample size, little difference in the probability of treatment by stratum, or little theoretical reason to expect treatment effect heterogeneity, investigators may have cause to prefer OLS if efficiency concerns are of paramount interest. In those cases, diagnostics or interpretational aids that gauge the severity of the misspecification bias (weights) may be useful, depending on how investigators prefer to tradeoff bias and efficiency in their inferential or decision-making goals.
Finally, when the linearity in potential outcomes assumption is not satisfied, none of these methods can guarantee unbiasedness. Our work here regards only the relaxation from single linearity to separate linearity. Investigators not satisfied with the linear approximation to treatment effects in this way could reasonably consider more flexible approaches---though doing so may incur additional uncertainty costs as illustrated in Table (ref). Nevertheless, we agree with assessments such as keele_adjusting_2010 and aronow_does_2016 that current practice largely remains reliant on linear approximations to adjust for covariates while hoping to interpret the result as a meaningful causal effect. Our analysis suggests that, especially where there are enough observations relative to covariates to support imputation/g-computation/T-learner, interactive regression, or mean balancing, this may be preferable to suffering regression's “weighting problem”, as it avoids this form of bias under slightly weaker assumptions and requires only a small change to current practice.