2.1 Data generation process

Theories of statistical inference require an axiom about the assumed data generation process in hypothetical repeated samples. In the matching literature, existing theories of inference for matching assume (usually implicitly) simple random sampling, which we define formally as follows:

In this article, we offer a new theory of inference for matching that replaces Axiom A0’ with an axiom based on stratified random sampling. Stratification is a well-known technique in statistics that has had a role in matching since at least Cochran (Reference Cochran1968) (see also Rubin Reference Rubin1977). To ease the exposition below, we denote by ${\mathcal{X}}$ the space of pretreatment covariates and offer a formal definition of stratification as follows.

For example, suppose that ${\mathcal{X}}$ consists of the variables $\mathtt{age}$ , $\mathtt{gender}$ , and earnings, that is, ${\mathcal{X}}=\{\mathtt{age},\mathtt{gender},\mathtt{earnings}\}$ . Then $\unicode[STIX]{x1D6F1}({\mathcal{X}})$ can be interpreted as the product (space) of variables $\mathtt{age}\times \mathtt{gender}\times \mathtt{earnings}=\unicode[STIX]{x1D6F1}({\mathcal{X}})$ . Therefore, in the example, one of the sets $A_{k}$ might be the subset of “young adult males making greater than $25,000,” that is, . When no ambiguity is introduced, we drop the subscript $k$ from $A_{k}$ . Stratified random sampling involves random sampling from within strata $A$ with given quotas proportional to the relative weight of the strata $\{W^{A},A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})\}$ .

Finally, we offer our alternative data generating process axiom.

In this alternative Axiom A0, the strata and the total number of observations for each hypothetical repeated sample and the observed sample are fixed. Then, the data set within each stratum is drawn according to simple random sampling from Axiom A0’.Footnote ¹

The axioms described in Axioms A0 and A0’ cannot be proven true or false on the basis of comparisons to a single observed data set, arguments about plausibility, or information about how matching methods are used. Because the repeated samples are strictly hypothetical, A0 and A0’ are not even statements that could be true or false in principle. Instead, the choice of an axiom merely defines how to interpret one’s causal inferences and uncertainty estimates, the specific type of repeated hypothetical samples, and the ultimate inferential target. As all matching methods use some kind of stratification of the covariates ${\mathcal{X}}$ , Axiom A0 highlights this fact and clarifies the theoretical assumptions necessary for valid inferences, rather than, as under Axiom A0’, keeping it hidden and leaving to applied researchers to deal with outside of the process of statistical inference.

2.2 Treatment assignment

Consider now the data generated in Axiom A0, where subject $i$ ( $i=1,\ldots ,n$ ) has been exposed to treatment $T_{i}=t$ , for $t\in {\mathcal{T}}$ , where ${\mathcal{T}}$ is either a subset of $\mathbb{R}$ or a set of (ordered or unordered) categories, $T_{i}$ is a random variable, and $t$ one possible value of it. Then ${\mathcal{Y}}=\{Y_{i}(t):t\in {\mathcal{T}},i=1,\ldots ,n\}$ is the set of potential outcomes, the possible values of the outcome variable when $T$ takes on different values. For each observation, we observe one and only one of the set of potential outcomes, that for which the treatment was actually assigned: $Y_{i}\equiv Y_{i}(T_{i})$ . In this setup, $T_{i}$ is a random variable, the potential outcomes are fixed constants for each value of $T_{i}$ , and $Y_{i}(T_{i})$ is a random variable, with randomness stemming solely from the data generation process for $T$ determining which of the potential outcomes is observed for each $i$ . Let $X_{i}$ be the $p\times 1$ vector ( $X\in {\mathcal{X}}$ ) of pretreatment covariates for subject $i$ .Footnote ²

2.3 Treatment effect

Let $t_{1}$ and $t_{2}$ be distinct values of $T$ that happen to be of interest, regardless of whether $T$ is binary, multicategory, or continuous (and which, for convenience, we refer to as the treated and control conditions, respectively). Assume $T$ is observed without error (an assumption we relax in Appendix B). Define the treatment effect for each observation as the difference between the corresponding two potential outcomes, $\text{TE}_{i}=Y_{i}(t_{1})-Y_{i}(t_{2})$ , of which at most only one is observed (this is known as the “Fundamental Problem of Causal Inference”; Holland Reference Holland1986). (Problems with multiple or continuous values of treatment variables have multiple treatment effects for each observation, but the same issues apply.)

The object of statistical inference is usually an average of treatment effects over a given subset of observations. Researchers then usually estimate one of two types of quantities. The first is the sample average treatment effect on the treated (SATT), for which the potential outcomes and thus TE $_{i}$ are considered fixed, and inference is for all treated units in the sample at hand: $\text{SATT}=\frac{1}{|M_{1}|}\sum _{i\in M_{1}}\text{TE}_{i}$ , with the control units used to help estimate this quantity (Imbens Reference Imbens2004, p. 6). Other causal quantities of this first type are averaged over different subsets of units, such as from the population, the subset of the population similar to $X$ , or all units in the sample or population regardless of the value of $T_{i}$ . Since a good estimate of one of these quantities will usually be a good estimate of the others, usually little attention is paid to the differences for point estimation, although there may be differences with respect to uncertainty estimates under some theories of inference (Imai Reference Imai2008; Imbens and Wooldridge Reference Imbens and Wooldridge2009).

The second type of causal quantity is when some treated units have no acceptable matches among a given control group and so are pruned along with unmatched controls, a common situation, which gives rise to “feasible” versions of SATT (which we label FSATT) or of the other quantities discussed above. This formalizes the common practice in many types of observational studies by focusing on quantities that can be estimated well (perhaps, in addition to estimating a more model dependent estimate of one of the original quantities) (see Crump et al. Reference Crump, Hotz, Imbens and Mitnik2009; Rubin Reference Rubin2010; Iacus, King, and Porro Reference Iacus, King and Porro2011), an issue we return to in Section 3.2. (In multilevel treatment applications, the researcher must choose whether to keep the feasible set the same across different treated units so that direct comparison of causal effects is possible, or to let the sets vary to make it easier to find matches.)

2.4 Assumptions

We now describe Assumptions A1–A3, which establish the theoretical background needed to justify valid causal inference under finite data with stratified random sampling as defined in Axiom A0; this set of assumptions can be seen as a natural stratum-wide extension of the pointwise theory by Rosenbaum and Rubin (Reference Rosenbaum and Rubin1983), which differs because it builds off Axiom A0’ instead.

The first assumption (which we generalize in Appendix B) helps to precisely define the variables used in the analysis.

SUTVA can be interpreted in at least three ways (see VanderWeele and Hernan Reference VanderWeele and Hernan2012). The first is “logical consistency,” which connects potential outcomes to the observed values and thus rules out a situation where say $Y_{i}(0)=5$ if $T_{i}=1$ but $Y_{i}(0)=12$ if $T_{i}=0$ (Robins Reference Robins1986). The second is “no interference,” which indicates that the observed value $T_{i}$ does not affect the values of $\{Y_{i}(t):t\in {\mathcal{T}}\}$ or $\{Y_{j}(t):t\in {\mathcal{T}},\forall ~j\neq i\}$ (Cox Reference Cox1958). And finally, SUTVA requires that the treatment assignment process produce one potential outcome value for any (true) treatment value (Neyman Reference Neyman1935).

To use our theory to justify a matching method requires that the information in these strata, and the variables that generate them, be taken into account. The theory does not require that our specific formalization of these strata be used in a matching method, only that the information is controlled for in some way. This can be done by directly matching on $A$ , using some function of $A$ in covariates to control for, or some type of weighting that takes account of $A$ . An example is given in Section 5.

We now introduce the second assumption, which ensures that the pretreatment covariates defining the strata are sufficient to adjust for any biases. (This assumption serves the same purpose as the “no omitted variable bias” assumption in classical econometrics, but without having to assume a particular functional form.) Thus, by conditioning on the values of $X$ encoded in the strata $A$ , we define the following.

For example, under A2, the distribution of potential outcomes under control $Y(0)$ is the same for the unobserved treated units and as the observed control units; below, this will enable us to estimate the causal effect by using the observed outcome variable in the control group.

Apart from the sampling framework, Assumption A2 can be thought of as a degenerate version of the Conditioning At Random (CAR) assumption in Heitjan and Rubin (Reference Heitjan and Rubin1991) with conditioning fixed. CAR was designed to draw inferences from coarsened data, when the original uncoarsened data are not observed. In the present framework, $\unicode[STIX]{x1D6F1}({\mathcal{X}})$ represents only a stratification of the reference population and each stratum $A$ in that definition is fixed in repeated sampling. A special case of Assumption A2, with sets $A$ fixed to singletons (i.e., taking $A=\{X=x\}$ ), is known as “weak unconfoundedness” used under exact matching theory (Imbens Reference Imbens2000; Lechner Reference Lechner, Lechner and Pfeiffer2001; Imai and van Dyk Reference Imai and van Dyk2004) and first articulated in Rosenbaum and Rubin (Reference Rosenbaum and Rubin1983).

Finally, any matching theory requires a version of the common support assumption, that is, for any unit with observed treatment condition $T_{i}=t_{1}$ and covariates $X_{i}\in A$ , it is also possible to observe a unit with the counterfactual treatment, $T_{i}=t_{2}$ , and the covariate values in the same set $A$ . This assumption rules out trying to estimate the causal effect of United Nations (UN) interventions in civil wars even though the UN intervenes only when it is likely to succeed (King and Zeng Reference King and Zeng2006). In less extreme cases, it is possible to narrow the quantity of interest to a portion of the sample space (and thus the data) where common support does exist. More formally, we introduce this version that works under the stratified random sampling Axiom A0.

Assumptions A2 and A3 make the search for counterfactuals easier since all observations in the vicinity of (i.e., in the same strata as) a unit, rather than only those with exactly the same covariate values, are now acceptable matches. (The combination of the pointwise versions of both A2 and A3 is often referred to as “strong ignorability” (Rosenbaum and Rubin Reference Rosenbaum and Rubin1983; Abadie and Imbens Reference Abadie and Imbens2002).) Assumption A3 also requires that at least one treated and one control unit (or one in each treatment regime) appear within every stratum, and so A3 imposes constraints on the weights.

2.5 Identification of the treatment effect

We show here that Assumptions A1–A3 enable point identification of the causal effect in the presence of approximate matching. Identification for the expected value of this quantity can be established under the new assumptions by noting, for each $A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})$ , that

$$\begin{eqnarray}E\{Y(t)|A\}\stackrel{\mathbf{A2}}{=}E\{Y(t)|T=t,A\}=E\{Y|T=t,A\},\end{eqnarray}$$

which means that within set $A_{k}$ , we can average over the observed $Y$ corresponding to the observed values of the treatment $T$ rather than unobserved potential outcomes for which the treatment was not assigned. The result is that the average causal effect within the set $A$ , which we denote by $\unicode[STIX]{x1D70F}^{A}$ , can be written as the difference in two means of observed variables, and so is easy to estimate:

(1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}^{A}=E\{Y(t_{1})-Y(t_{2})|A\}=E\{Y|T=t_{1},A\}-E\{Y|T=t_{2},A\},\end{eqnarray}$$

for any $t_{1}\neq t_{2}\in {\mathcal{T}}$ . That is, (1) simplifies the task of estimating the causal effect in approximate matching in that it allows one to consider the means of the treated and control groups separately, within each set $A$ , and to take the weighted average over all strata $A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})$ afterward. To take this weighted average, we use Assumption A3:

(2) $$\begin{eqnarray}E(Y(t))\stackrel{\mathbf{A3}}{=}E(E\{Y(t)|A\}),\end{eqnarray}$$

which is exactly what we need to calculate the average causal effect $\unicode[STIX]{x1D70F}=E(Y(t_{1}))-E(Y(t_{2}))$ . Assumption A3 is required because otherwise $E\{Y(t)|A\}$ may not exist for one of the two values of $t=t_{1}$ or $t=t_{2}$ for some stratum $A$ , in which case $E(Y(t))$ , would not exist and the overall causal effect would not be identified.

3.1 Difference in means estimator

To describe the property of the estimators, we adapt the notation of Abadie and Imbens (Reference Abadie and Imbens2011) (which operates under Axiom A0’) and rewrite the causal quantity of interest as the weighted sum computed within each stratum $A$ from (1):

(3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F} & = & \displaystyle \frac{1}{m_{1}}\mathop{\sum }_{i\in M_{1}}E\{\text{TE}_{i}\}=\frac{1}{m_{1}}\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\mathop{\sum }_{i\in M_{1}^{A}}E\{Y_{i}(t_{1})-Y_{i}(t_{2})|X_{i}\in A\}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{m_{1}}\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\mathop{\sum }_{i\in M_{1}^{A}}(\unicode[STIX]{x1D707}_{1}^{A}-\unicode[STIX]{x1D707}_{2}^{A})=\frac{1}{m_{1}}\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}(\unicode[STIX]{x1D707}_{1}^{A}-\unicode[STIX]{x1D707}_{2}^{A})m_{1}^{A}=\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\unicode[STIX]{x1D70F}^{A}W_{1}^{A},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}_{k}^{A}=E\{Y(t_{k})|X\in A\}$ ( $k=1,2$ ) and $\unicode[STIX]{x1D70F}^{A}$ is the treatment effect within set $A$ as in (1). Consider now an estimator $\hat{\unicode[STIX]{x1D70F}}$ for $\unicode[STIX]{x1D70F}$ based on this weighted average:

(4) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D70F}}=\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\hat{\unicode[STIX]{x1D70F}}^{A}W_{1}^{A}=\frac{1}{m_{1}}\mathop{\sum }_{i\in M_{1}^{A}}(Y_{i}(t_{1})-{\hat{Y}}_{i}(t_{2})), & & \displaystyle\end{eqnarray}$$

where $\hat{\unicode[STIX]{x1D70F}}^{A}$ is the simple difference in means within the set $A$ , that is:

(5) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D70F}}^{A} & = & \displaystyle \frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}\left(Y_{i}-{\hat{Y}}_{i}(t_{2})\right)=\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}\left(Y_{i}-\frac{1}{m_{2}^{A}}\mathop{\sum }_{j\in M_{2}^{A}}Y_{j}\right)\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}Y_{i}-\frac{1}{m_{2}^{A}}\mathop{\sum }_{j\in M_{2}^{A}}Y_{j}.\end{eqnarray}$$

Finally, we have the main result (see Appendix A for a proof).

Given that the sets of the partition $\unicode[STIX]{x1D6F1}({\mathcal{X}})$ are disjoint, it is straightforward to obtain the variance $\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}}^{2}=\text{Var}(\hat{\unicode[STIX]{x1D70F}})$ of the causal effect. If we denote by $\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}^{A}}^{2}$ the variance of the stratum-level estimates $\hat{\unicode[STIX]{x1D70F}}^{A}$ in (5), we have $\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}}^{2}=\sum _{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\left(\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}^{A}}W_{1}^{A}\right)^{2}$ .

3.2 Simplified inference through weighted least squares

The direct approach to estimating the treatment effect by strata and then aggregating is useful to define the matching estimator, but it is more convenient to rewrite the estimation problem in an equivalent way as a weighted least squares (WLS) problem. This approach provides a easy procedure for computing standard errors, even for multilevel treatment (see Section 3.3) or when one or more strata contain only one treated unit and one control unit (see Section 4). In this latter case, one cannot directly estimate the variance within the strata $\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}^{A}}^{2}$ but we can still obtain an estimate of it by applying whatever estimator one would have applied to the data set without matching.

We now introduce the weights we use to simplify the estimator in (4) and re-express them as the difference in weighted means. For all observations, define the weights $w_{i}$ as

(6) $$\begin{eqnarray}w_{i}=\left\{\begin{array}{@{}ll@{}}1,\quad & \text{if }T_{i}=t_{1},\\ 0,\quad & \text{if }T_{i}=t_{2}\text{ and }i\not \in M_{2}^{A}\text{ for all }A,\\ {\displaystyle \frac{m_{1}^{A}}{m_{2}^{A}}}{\displaystyle \frac{m_{2}}{m_{1}}},\quad & \text{if }T_{i}=t_{2}\text{ and }i\in M_{2}^{A}\text{ for one }A.\end{array}\right.\end{eqnarray}$$

Then, the estimator $\hat{\unicode[STIX]{x1D70F}}$ in (4) can be rewritten as

$$\begin{eqnarray}\hat{\unicode[STIX]{x1D70F}}=\frac{1}{m_{1}}\mathop{\sum }_{i\in M_{1}}Y_{i}w_{i}-\frac{1}{m_{2}}\mathop{\sum }_{j\in M_{2}}Y_{j}w_{j},\end{eqnarray}$$

where the variance is the sum of the variances of the two quantities. Therefore, the standard error of $\hat{\unicode[STIX]{x1D70F}}$ is the usual standard error of estimates for regression analysis with weights. For example, consider the linear regression model:

$$\begin{eqnarray}Y_{i}=\unicode[STIX]{x1D6FD}_{0}+\unicode[STIX]{x1D6FD}_{1}T_{i}+\unicode[STIX]{x1D716}_{i},\quad \unicode[STIX]{x1D716}_{i}\sim N(0,1)\text{ and i.i.d.},\end{eqnarray}$$

where $\hat{\unicode[STIX]{x1D70F}}\equiv \hat{\unicode[STIX]{x1D6FD}}_{1}$ if weights $w_{i}$ are used in estimation. So the standard error of $\hat{\unicode[STIX]{x1D6FD}}_{1}$ can be obtained as the output of this WLS model, and is the correct estimate of $\unicode[STIX]{x1D70E}_{\hat{\unicode[STIX]{x1D70F}}}^{2}$ . Other models, such as generalized linear models (GLM) with weights, can also be estimated in a similar fashion. The only change that needs to be made to the estimator without matching is to include these weights.

3.3 Estimation with multilevel treatments

For more than two treatments we define the multitreatment weights as

$$\begin{eqnarray}w_{i}(k)=\left\{\begin{array}{@{}ll@{}}1,\quad & \text{if }T_{i}=t_{1},\\ 0,\quad & \text{if }T_{i}=t_{k}\text{ and }i\not \in M_{k}^{A}\text{ for all }A,\\ {\displaystyle \frac{m^{A}}{m_{k}^{A}}}{\displaystyle \frac{m_{k}}{m_{1}}},\quad & \text{if }T_{i}=t_{k}\text{ and }i\in M_{k}^{A}\text{ for one }A.\end{array}\right.\end{eqnarray}$$

Then, for each $k=2,3,\ldots$ , the treatment effect $\unicode[STIX]{x1D70F}(k)$ can be estimated as $\hat{\unicode[STIX]{x1D6FD}}_{1}(k)$ in

$$\begin{eqnarray}Y_{i}=\unicode[STIX]{x1D6FD}_{0}+\unicode[STIX]{x1D6FD}_{1}(k)T_{i}+\cdots +\unicode[STIX]{x1D716}_{i}\end{eqnarray}$$

with weights $w_{i}(k)$ and, again, the usual standard errors are correct as is.

3.4 Additional regression adjustment for further covariates

If Assumption A2 holds, then adjusting for covariates is unnecessary to ensure unbiasedness. If Assumption A2 holds but the analyst is unsure if it does, and so adjusts for pretreatment covariates (with interactions), then the downside is trivial (Lin Reference Lin2013; Miratrix, Sekhon, and Yu Reference Miratrix, Sekhon and Yu2013). But sometimes, the researcher may need to adjust for covariates via a model, even if they were not used during matching. In this situation, it is sufficient to proceed as one would without the matching step by including all the variables in the regression model

$$\begin{eqnarray}Y_{i}=\unicode[STIX]{x1D6FD}_{0}+\unicode[STIX]{x1D6FD}_{1}T_{i}+\unicode[STIX]{x1D6FE}_{1}X_{i1}+\cdots \unicode[STIX]{x1D6FE}_{d}X_{id}+\unicode[STIX]{x1D716}_{i},\end{eqnarray}$$

and using the weights in (6) to run the WLS regression. The estimated coefficient $\hat{\unicode[STIX]{x1D6FD}}_{1}$ is then the estimator of the treatment effect $\hat{\unicode[STIX]{x1D70F}}$ and its standard error is an unbiased estimator of its standard deviation, under the model.

3.5 Defining strata in observational data

One question that may arise in this framework, as in any stratified sampling, is how to choose the strata $A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})$ in a given problem? The answer is by definition application-specific, which can be an advantage in that it relies on variables in the original investigator-defined units of measurement, reflecting knowledge the investigator must have.

To show the applicability of our approach in observational studies, we take advantage of the fact that in many data sets variables referred to as “continuous” in fact often have natural breakpoints that may be as or more important than the continuous values. These may include grade school, high school, and college degrees for the variable “years of education”; the official poverty level for the variable “income”; or puberty, official retirement age, and so on, for the variable “age.” This understanding of measurement recognizes that, for another example, $33\,^{\circ }\text{F}$ may be closer to $200^{\circ }$ than to $31^{\circ }$ , at least for certain purposes. Most data analysts not only know this distinction well but use it routinely to collapse variables in their ordinary data analyses. Indeed, in analyses of sample surveys, continuous variables with no natural breakpoints, or where authors never use breakpoints to collapse variables or categories, are uncommon.

For another example, consider estimating the causal effect of the treatment variable “taking one introductory statistics course” on the outcome variable “income after college,” and where we also observe one pretreatment covariate “years of education,” along with its natural breakpoints at high school and college degrees. Assumption A2 says that it is sufficient to control for the coarsened three-category education variable (no high school degree, high school degree and possibly some college courses but no college degree, and college degree) rather than the full “years of education” variable. In this application, A2 is plausible if, as seems common, employers at least at first primarily value degree completion in setting salaries. Then, poststratification and matching within the strata is appropriate. If, instead, a major difference in expected income exists between those who have, for example, one versus three years of college, then there can be some degree of bias induced.

4.1 How bias arises?

To understand where bias may arise under Axiom A0 when some strata $A$ need to be enlarged or changed, we study the following bias decomposition, by adapting ideas designed to work under Axiom A0’ from Abadie and Imbens (Reference Abadie and Imbens2006, Reference Abadie and Imbens2011, Reference Abadie and Imbens2012). Let $\unicode[STIX]{x1D707}_{t}(x)=E\{Y(t)|X=x\}$ and $\unicode[STIX]{x1D707}(t_{k},x)=E\{Y|X=x,T=t_{k}\}$ . Under Assumption A2 we know that $\unicode[STIX]{x1D707}_{t_{k}}(x)\stackrel{\mathbf{A2}}{=}\unicode[STIX]{x1D707}(t_{k},x)\equiv \unicode[STIX]{x1D707}_{k}^{A}$ for all $\{X=x\}\subseteq A$ . Then the bias is written as:

$$\begin{eqnarray}\hat{\unicode[STIX]{x1D70F}}^{A}-\unicode[STIX]{x1D70F}^{A}=\mathop{\sum }_{A\in \unicode[STIX]{x1D6F1}({\mathcal{X}})}\{(\bar{\unicode[STIX]{x1D70F}}^{A}-\unicode[STIX]{x1D70F}^{A})+E^{A}+B^{A}\}W_{1}^{A},\end{eqnarray}$$

where

$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D70F}}^{A}=\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}(\unicode[STIX]{x1D707}_{t_{1}}(X_{i})-\unicode[STIX]{x1D707}_{t_{2}}(X_{i})) & \displaystyle \nonumber\\ \displaystyle & \displaystyle E^{A}=\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}\left((Y_{i}-\unicode[STIX]{x1D707}_{t_{1}}(X_{i}))-\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}\frac{1}{m_{2}^{A}}\mathop{\sum }_{j\in M_{2}^{A}}(Y_{j}-\unicode[STIX]{x1D707}_{t_{2}}(X_{j}))\right) & \displaystyle \nonumber\end{eqnarray}$$

and

$$\begin{eqnarray}B_{A}=\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}\frac{1}{m_{2}^{A}}\mathop{\sum }_{j\in M_{2}^{A}}(\unicode[STIX]{x1D707}_{t_{2}}(X_{i})-\unicode[STIX]{x1D707}_{t_{2}}(X_{j})),\end{eqnarray}$$

where $\unicode[STIX]{x1D707}_{t_{k}}(X)=\unicode[STIX]{x1D707}_{k}^{A}$ for $X\in A$ . Therefore, both $(\bar{\unicode[STIX]{x1D70F}}^{A}-\unicode[STIX]{x1D70F}^{A})$ and $E^{A}$ have zero expectation inside each set $A$ and $B^{A}=0$ . But if some of the sets $A^{\prime }$ are different from the original partition $A$ , or combined or enlarged, then Assumption A2 may not apply any longer; in general, $\unicode[STIX]{x1D707}_{t_{k}}(X)\neq \unicode[STIX]{x1D707}_{k}^{A}$ for $X\in A^{\prime }\neq A$ .

4.2 Nonparametric regression adjustment

One way to proceed is with the following regression adjustment, as in Abadie and Imbens (Reference Abadie and Imbens2011), that compensates for the bias due to the difference between $A$ and $A^{\prime }$ . Let $\hat{\unicode[STIX]{x1D707}}_{t_{2}|A}(x)$ be a (local) consistent estimator of $\unicode[STIX]{x1D707}_{t_{2}}(x)$ for $x\in A$ . In this case, one possible estimator is the following:

(7) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D70F}}^{A}=\frac{1}{m_{1}^{A}}\mathop{\sum }_{i\in M_{1}^{A}}(Y_{i}-\hat{\unicode[STIX]{x1D707}}_{t_{2}|A}(X_{i}))-\frac{1}{m_{2}^{A}}\mathop{\sum }_{j\in M_{2}^{A}}(Y_{j}-\hat{\unicode[STIX]{x1D707}}_{t_{2}|A}(X_{j})).\end{eqnarray}$$

This estimator is asymptotically unbiased if the number of control units in each stratum grows at the usual rate. If instead of using a local estimator $\hat{\unicode[STIX]{x1D707}}_{t_{2}|A}(x)$ we use a global estimator $\hat{\unicode[STIX]{x1D707}}_{t_{2}}(x)$ , that is, using all control units in the sample as in Abadie and Imbens (Reference Abadie and Imbens2011), then the calculation of the variance of the estimator is no longer obtained by simple weighting and the validity of the approach requires a treatment similar to the asymptotic theory of exact matching. More technical assumptions and regularity on the unknown functions $\unicode[STIX]{x1D707}_{t}(x)$ are needed to prove that the regression type estimator in (7) can compensate for the bias asymptotically but, essentially, it is required that, for some $r\geqslant 1$ , we impose $m_{1}^{r}/m_{2}\rightarrow \unicode[STIX]{x1D705}$ , with $0<\unicode[STIX]{x1D705}<\infty$ . A simplified statement is that $m_{1}/m_{2}^{4/k}\rightarrow 0$ , where $k$ is the number of continuous covariates in the data and this condition is equivalent to $m_{1}^{k/4}/m_{2}=m_{1}^{r}/m_{2}\rightarrow \unicode[STIX]{x1D705}$ . The proof of these results can be found in Abadie and Imbens (Reference Abadie and Imbens2011).

4.3 Asymptotic filling of the strata

If Assumption A3 is apparently violated because there are not enough observations in one or more strata, but we still believe A1–A3 to be true and we happen to be able to continue to collect data, then it is worth knowing that it is theoretically possible to fill all the strata in $\unicode[STIX]{x1D6F1}({\mathcal{X}})$ and obtain unbiased estimates of the treatment effect. This is theoretically possible under the additional assumption that $m_{1}^{r}/m_{2}\leqslant \unicode[STIX]{x1D705}$ , with $0<\unicode[STIX]{x1D705}<\infty$ , $r>k$ , and $k$ the number of continuous covariates. By Proposition 1 in Abadie and Imbens (Reference Abadie and Imbens2012), all the strata $A$ will be filled with probability one. This result is enough to obtain asymptotically unbiased estimates of the causal effect under the original assumptions A2–A3, without changing the initial partition $\unicode[STIX]{x1D6F1}({\mathcal{X}})$ or other technical smoothness assumptions on the functions $\unicode[STIX]{x1D707}_{t}(x)$ and $\hat{\unicode[STIX]{x1D707}}_{t|A}(x)$ . As such, one could use an asymptotic approximation to obtain estimates and standard errors, but it is considerably safer to use these results as a guide to future data collection.

5.1 Data

For data, we consider the National Supported Work (NSW) training program used in the seminal paper by Lalonde (Reference Lalonde1986). In these data, the outcome variable is the real earnings of workers in 1978 ( $\mathtt{re78}$ ) and the treatment is the participation in the program. Pretreatment control variables include $\mathtt{age}$ ; years of $\mathtt{education}$ ; indicators for $\mathtt{black}$ and $\mathtt{hispanic}$ ; an indicator for marital status, $\mathtt{married}$ ; an indicator for not possessing a high school degree ( $\mathtt{nodegree}$ ); earnings in 1975, $\mathtt{re75}$ , and 1974, $\mathtt{re74}$ ; and unemployment status in both years, $\mathtt{u74}$ and $\mathtt{u75}$ . The data set contains 297 individuals exposed to treatment and 425 control units. This is in fact an experimental sample, although Lalonde (Reference Lalonde1986) analyzed it as an observational study to provide insights about matching methods. The quantity of interest is the sample average treatment effect on the treated (ATT), the increase in earnings in 1978 due to treatment.

5.2 Applying simple random sampling-based theory

We now show how three ways of satisfying the existing simple random sampling-based theory of inference all fail in these data. We begin with the simple random sampling Axiom A0’, and assume SUTVA, along with pointwise unconfoundedness and common support.

First, we apply exact matching theory, the hardest to satisfy but best case scenario, requiring exact matching on $X$ . In most applications with rich sets of covariates, few matches are available. In our application, we do happen to have 55 treated units that exactly match 74 control observations (this occurs because the otherwise continuous variables, $\mathtt{re74}$ and $\mathtt{re75}$ , have a point mass at zero, and other variables are discrete). Unfortunately, this small sample would leave confidence intervals on the ATT too wide to make useful inferences.

More generally, this exact matching approach may work for very large samples, when there is a high probability that match occurs without replacement for one-to-one nearest neighbor matching, and imbalance is zero (or very small according to some distance). In this situation, a simple regression model, with pretreatment covariates and a treatment indicator, will normally be able to take into account the remaining bias in either simple or stratified random sampling.

Second, we consider exact matching on the propensity score. If successful, this approach would yield less efficient estimates than exact matching on $X$ , and would introduce a variety of other serious problems, but causal estimates would at least be ex ante unbiased (King and Nielsen Reference King and Nielsen2017). To try this, we use the propensity score specification in Dehejia and Wahba (Reference Dehejia and Wahba1999), a logistic model for all the indicator variables, as well as $\mathtt{age}$ , $\mathtt{education}$ , $\mathtt{re74}$ , and $\mathtt{re75}$ and their squares. Unfortunately, as is typical in data sets with continuous covariates, lowering the bar for what constitutes a match in this way buys us zero additional matched observations. This is not a surprise, since propensity score matching requires exact matching on the propensity score, which does not happen with any higher probability than exact matching on $X$ as long as we have some continuous variables.

A final option to follow existing theory would be to have a very large data set. Although we do not have a large data set, in observational data analysis, the data set is whatever one chooses to include. In this case, we could add new data by gathering contemporaneous surveys on the same subject, of similar people, and treat them as part of the pool of potential control units (see Dehejia and Wahba Reference Dehejia and Wahba1999). In fact, adding external data has been tried in this application but turns out not to help because it greatly increases heterogeneity, does not markedly increase the information in the data as $n$ increases about the ATT, and so does not satisfy the conditions for the theory of inference to apply (see Smith and Todd Reference Smith and Todd2005; King, Lucas, and Nielsen Reference King, Lucas and Nielsen2017).

5.3 Applying approximate simple random sampling-based theory

At this point, we can see that applying an existing theory of inference based on simple random sampling, to generate valid causal inferences in these data without approximations, is not possible. Researchers in this situation typically try to come up with an approximate matching solution, but this leads to two problems.

First, approximate matching is not justified by the simple random sampling-based theory of inference, as the formal properties of the resulting statistical estimators do not hold. Second, one might think that small deviations from the theoretical requirement would be approximately unbiased, but this is untrue. No known theorem supported this claim and even exactly matched propensity scores imply only approximate matching on $X$ .

By looking at how imbalanced a data set is, we can get a feel for at least the potential bias due to failing to exactly matching on $X$ or on the propensity score. To illustrate, we use one-to-one nearest neighbor propensity score matching (NN-PSM) with a caliper of 0.001. This results in 100 treated units and 100 control units. The closest (inexact) match allowed by this procedure has a difference in propensity scores of only 0.000003, but yet still has substantial imbalance:

As can be seen, $\mathtt{education}$ (at 9 years for treated and 8 for control) is not far off, at least for a job training program. Also apparently close is $\mathtt{age}$ but, at 20 and 23, the impact could be determinative if the legal age of adulthood (21) impacts prospective employers’ hiring decisions. More serious is that the treated person in this pair is $\mathtt{hispanic}$ and employed in both 1974 and 1975, whereas the control person is $\mathtt{black}$ and unemployed in the same years. In this typical example, a practitioner would have to implicitly admit to not controlling for several of the variables they designated as pretreatment confounders, thus violating ignorability or to hope that the biases of one confounder miraculously cancel out the biases for another.

We also repeat here the same analysis for a larger caliper to further increase the number of matched units, and show its cost in terms of increasing imbalance further. Here we choose a larger caliper of 0.01, and find 219 treated units matched with 219 control units. The next best pair of matched units this brings in has a “small” propensity score difference of 0.000622, but with an obviously large imbalance on $X$ :

The differences between the treated and control groups on $X$ here are even more substantial, even with an only slightly larger caliper. Here we match a treated African American who dropped out of junior high school with no income, to a control group Hispanic who graduated from high school with more than $27,000 of income.

The two matched pairs of units we describe here are each intuitive and the degree of approximation is easy to understand. However, to understand the full degree of approximation for the entire matching solution requires performing this identical comparison on every pair of matched observations (100, 219, or 274, depending on the choice of caliper).

Although we do not offer an example until later, we also note that running NN-PSM with a caliper of 0.1 matches 274 treated units to 274 control units (i.e., 548 units).

As is clear from this discussion, the size of the propensity score caliper alone provides little intuition about the quality of the match, the degree of approximation to the requirements of the theory of inference, the resulting level of imbalance on $X$ , or the degree of statistical bias in the ATT. Since the required Axiom A0’ is an axiom, being able to clearly convey what it means is the only real requirement in applications.

5.4 Applying stratified random sampling-based theory

We now illustrate three advantages of replacing Axiom A0’ with Axiom A0, and thus thinking of the data in terms of stratification rather than simple random sampling.

The first advantage of stratification-based matching is ease of interpretation, which is essential in matching. Understanding assumptions—which by definition are unverifiable—is important in any empirical analysis. However, the sampling axiom in matching defines the inferences being made and thus also the meaning of the sampling distribution and standard errors. As such, without some clear understanding of the sampling process, making any inferences at all makes little sense.

To convey how one would interpret an application under this stratification-based theory of inference, and why matching is easier to understand than under simple random sampling-based approaches, we now give an example involving analysis choices in a real data set. For stratification-based inference, the key choice is the partition of $X$ , which we have been referring to as $A$ . In principle, this choice must be made prior to examining the data, or else the weights will not be fixed in repeated samples. We discuss different ways of interpreting this requirement so that it may be used in practice when not generating the data oneself.

A reasonable way to define $A$ , before seeing the data, is to define it based on information in the data set’s codebook. In the case of these data, a natural choice is to match all binary variables exactly, $\mathtt{age}$ according to the official US Bureau of Labor Statistics stratification (i.e., 16–24, 25–54, 55 and over), and for the variable $\mathtt{education}$ to coarsen by formal degrees—elementary [0,6], high school [7–12], undergraduate [13–16], and graduate [17,). The covariate $\mathtt{u74}$ (and $\mathtt{u75}$ ) is an indicator variable, which is nonzero when $\mathtt{re74}$ (and $\mathtt{re75}$ ) is nonzero. As a result, this continuous counterpart of the unemployment status ( $\mathtt{re74}$ and $\mathtt{re75}$ ) can be in principle dropped from the matching stage and eventually included in the model specification for the ATT estimation.

If we use this definition of $A$ , which we could plausibly have arrived at before examining the data, then we can think of the data generation process as stratified random sampling within this given partition. Then all hypothetical repeated samples, the resulting sampling distribution, and associated standard errors, confidence intervals, and hypothesis tests are defined as a consequence. As it happens, when we can try this stratification with our one observed data set, we find that we have 221 treated units matched with 313 control units.

Now suppose we examine the data and prefer to drop $\mathtt{re74}$ and $\mathtt{re75}$ from the partition in order to prune fewer observations. To make this decision statistically justifiable, two conditions must hold. The first condition is that we must ensure we do not violate set-wide ignorability (i.e., Assumptions A2 and A3), which will be satisfied in one or more of three situations: the two variables are unrelated to the treatment variable, unrelated to the outcome given the treatment, or included in an appropriate model during the estimation stage. The second condition involves conceptualizing the resulting strata $A$ . If the choice of $A$ is determined from the data (not merely the codebook), then the weights are random and, as a result, more complicated methods must be used for uncertainty estimates (point estimates remain unchanged). However, we may still be able to conceptualize $A$ as fixed ex ante if we can argue that we would have interpreted the partitions the same way if we had thought of the same reasoning before seeing the data. That is, sometimes seeing the data causes one to surface ideas that could easily have been specified ex ante. Of course, we should try to avoid the lure of post hoc, just-so stories, but if we do, we would be justified in interpreting $A$ as fixed, and then all the familiar methods are available for computing uncertainty estimates, such as standard errors and confidence intervals. Either way, the advantage of stratification is that understanding the sampling axiom, and how to think about the resulting data generation process, is straightforward. In this case, dropping $\mathtt{re74}$ and $\mathtt{re75}$ results in matching 278 treated units matched to 394 control units.

Now suppose we go another step and try to interpret our analyses without much prior knowledge of $A$ . Here, we first generate a large set of matching solutions. This need not be done in practice, but we find it useful here for illustrative purposes. To do this, we create 500 random stratifications of the covariate space by dividing the support of each pretreatment covariate into a random number of subintervals (chosen uniformly on the integers $1,\ldots ,15$ ). For comparison, we also generate 500 matching solutions from NN-PSM models, by randomly selecting propensity score models and its caliper, and 500 nearest neighbor, Mahalanobis distance matching (NN-MDM) solutions, by randomly selecting input variables and its caliper (both selecting from the set of all main, polynomial, and interaction terms up to the second degree, with a logistic specification as usual for the propensity score model). In real applications, imbalance is best measured on the scale of the original variables but, to save space for our methodological purposes, we use the average of the standardized difference in means applied to each matching variable. (Other measures of imbalance do not materially change our conclusions.) Then, in Figure 1, the vertical axis is this measure of imbalance and the horizontal axis is the number of matched units (scaled according to $1/n$ ). Each point in the plot corresponds to one randomly selected matching solution (with stratification solutions in blue, NN-PSM in green, and NN-MDM in red). Stratification solutions here are all based on coarsened exact matching, but our stratified theory of inference applies to any member of the class of “monotonic imbalance bounding” matching methods (Iacus, King, and Porro Reference Iacus, King and Porro2011).

Figure 1. Randomly created matching solutions: imbalance (vertically) by matched sample size (horizontally). Each dot represents a matching solution, including the original unmatched data (raw); three solutions with lower imbalance and more matched treated units than any other (best random stratifications, at the lower left); NN-PSM solutions with calipers of 0.1, 0.01, and 0.001; 500 solutions based on random stratifications of the covariate space (dots); 500 random NN-PSM solutions (stars), and 500 random NN-MDM solutions (squares). The plot represents solutions with at least 200 matched units.

The $\mathtt{raw}$ dot (at the middle left) corresponds to the original, unmatched data. In the same figure, we also include and label the three different NN-PSM matching solutions discussed above with calipers set to 0.1, 0.01, and 0.001 (across the middle of the graph from left to right). The dots marked with “best stratifications” (at the bottom left) represent matching solutions based on stratifications with the lowest imbalance for a given number of matched units or the largest number of matched units for any given imbalance. These solutions do not necessarily represent the theoretical frontier of imbalance and matched sample size, since they were generated randomly, but they are the best solutions among those in this graph (King, Lucas, and Nielsen Reference King, Lucas and Nielsen2017) and are still the best among those randomly generated.

Then, to convey how easy it is to understand a stratified matching solution, consider only the central dot of this sequence of “best stratifications.” This matching solution was constructed (by chance, i.e., randomly) using the cross products of the following strata.

The advantage of presenting these strata is that they convey all information necessary about the entire matching solution in an easy-to-understand and compact display. To use our stratified theory of inference, we need to imagine that our data, and all the repeated hypothetical samples, were generated by a stratified sampling design, based on these strata. In fact, the data are observational, and the hypothetical distributions do not and will not exist. However, we can still conceptualize what this distribution means as if these strata are fixed. The argument should be recognized as more of a stretch, since we did arrive at this stratification directly from the data, but stating this axiom about hypothetical (stratified) sampling replications cannot be proven wrong and so it is reasonable to use it as a way to interpret a matching-based estimator. Our main point here is that the axiom itself is easy to understand: all we need to do is to understand the strata defined above.

In contrast, to convey all information in an NN-PSM matching solution, we would need to understand every individual matched set in a matching solution, as we did above, but for 100, 219, or 274 individual matched sets (corresponding to calipers of 0.001, 0.01, and 0.1). This of course would be infeasible to comprehend all at once. With stratification, a researcher can more easily, quickly, and concisely understand the approximation and what assumptions are necessary to believe that bias is being constrained, without hundreds of separate evaluations. These stratifications might still be implausible as ex ante definitions for $A$ , but researchers will be able to understand, if appropriate justify, the assumptions more easily. In this example, we can see that this particular matching solution does not control for $\mathtt{black}$ or $\mathtt{re75}$ , as was the case for a particular pair of NN-PSM matched sets above, but this time we can see all the compromises from the entire data set at once, so that one can judge whether this approximation is justifiable. The problem of course is that even if one can understand a hundred or more stratifications, the axiom of simple random sampling requires exact matching on $X$ or the propensity score, not a nearest neighbor solution, or one within some caliper.

A second advantage of stratification-based matching is that imbalance tends to be lower than under other matching methods for any given number of matched observations. This is not a general claim, but it is a typical pattern in many applications (King, Lucas, and Nielsen Reference King, Lucas and Nielsen2017). This can be seen in Figure 1 by the blue stratification-based matching points appearing to the lower left, indicating lower imbalance given higher numbers of observations, whereas the green NN-PSM and red NN-MDM matching solutions appearing above and to the right indicate more imbalance or fewer matched observations.

A final advantage of stratification-based matching is that estimated treatment effects are often less variable, and thus somewhat more robust, than under other matching methods. To see this common, but also not universal, pattern we compute, for each of the matching solutions in Figure 1, an estimate of the ATT and standard error (by regressing the outcome variable on the treatment indicator and all pretreatment variables included in the matching solution). We then present, in Figure 2, all the ATT estimates (vertically) by the matched sample size (horizontally). Because of the enormous variability of NN-MDM, the plot is zoomed in, excluding points outside of the range of the visible axes.

Figure 2. Estimated ATT (vertically) by number of matching units (horizontally), for matching solutions in Figure 1. The ellipses are convex hulls of the plotted points and represent roughly the dispersion of the ATT estimates for each matching method. The plot represents solutions with at least 200 matched units.

This figure shows that for any given matched sample size (i.e., any point on the horizontal axis), the vertical variability of the ATT estimates is much larger for NN-PSM and NN-MDM than for stratification. Because Figure 2 does not reveal all the data, we present Table 1, which summarizes aspects of these same estimates. Table 1 demonstrates that stratification, in addition to having lower overall imbalance (i.e., finding better matched subgroups of treated and control units), is also the method that on average produces more matched units, a less variable and more robust ATT estimate, with a smaller standard error. In addition, the one-sigma Monte Carlo confidence interval for average treatment effect contains zero for both PSM and MDM, but not under stratification.

Table 1. Distribution of estimated ATTs, their standard error, and number of matched units for the data in Figure 2.

As we exemplify with this analysis, the choice of a theory of inference defines the nature of the hypothetical repeated samples used for statistical inference. Whether these samples are based on simple or stratified random sampling is not an assumption vulnerable to being proven wrong, but rather than an axiom that defines how we interpret standard errors and confidence intervals. As such, the critical question is not which is more appropriate but whether we are able to clarify the meaning of one’s uncertainty calculations. As we show here, under stratified random sampling, the assumptions and inferences are considerably clearer and easier to understand, and do not require asymptotic results, which is quite unlike the situation with most methods for simple random sampling-based inference.