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.
36,631 characters · 12 sections · 31 citation commands
Matching on What Matters: A Pseudo-Metric Learning Approach to Matching Estimation in High Dimensions
While experimentation approaches to causal inference have strong theoretical guarantees, they are impractical in some settings. This leads researchers to rely on observational methods which suppose that conditioning on available covariates is sufficient to yield unbiased estimates. This can be especially tenuous (or prone to researcher manipulation) when the conditioning variables are selected in an ad hoc way. Recognizing this limitation, there has been growing interest in importing machine learning (ML) tools to automate and accelerate the model selection process athey2017beyond. Most examples of this work chernozhukov2016double,wager2018estimation rely on “black-box” ML algorithms that may obscure the underlying counterfactual reasoning. By contrast, matching methods yield a set of direct unit-to-unit comparisons that are fully transparent to policy makers, but are much less robust to the curse of dimensionality presented by high-dimensional data. Nonetheless, there has been relatively little work on the problem of model selection in high-dimensions for these methods.
Instead, propensity score matching (PSM) remains the most commonly used matching method used in applied economics and related fields. This approach side-steps the curse of dimensionality by collapsing all covariates to a single number reflecting probability of treatment McCaffrey2004,Wyss2014. However, PSM has been criticized because it does not guarantee covariate balance in the resulting sample, thus allowing significant model dependence and inefficiency in downstream estimations King2016. Prognostic score matching Hansen2008 also has the advantage of collapsing the covariates used for matching to a single dimension. While it does not explicitly balance covariates, it does balance on counterfactual outcomes, eliminating the potentially large efficiency losses associated with PSM.
Matching directly on the space of covariates--using methods like Mahalanobos Distance Matching (MDM) or Coarsened Exact Matching (CEM)--has gained popularity in fields such as Statistics, Epidemiology, Sociology, and Political Science ho2007matching. This has the advantage of not requiring correct specification of either a propensity or prognostic model. Additionally, these methods provide better covariate balance and, to the extent those covariates drive variance in outcomes, improve efficiency of downstream estimations when compared to PSM. Finally, these methods prevent the creation of unintuitive matches that happen when units with very different covariate values have similar estimated propensity (or prognostic) scores. However, MDM and CEM are particularly vulnerable to the curse of dimensionality discussed in abadie2006large and demonstrated in Gu1993. That is, even for the case of two continuous covariates, the non-parametric bias resulting from covariate differences between matched pairs can be large enough to make asymptotically valid inference impossible. Along with li2017matching, we explore a middle ground to recast the problem of matching in a metric learning framework that maps our covariates to a low-dimensional space that facilitates “closer matches” while still capturing important aspects of unit-level heterogeneity. These methods are a natural choice as they were initially designed to predict similarity in high-dimensional problems like facial recognition parkhi2015deep, person re-identification liao2015person, or image retrieval hoi2010semi.
In order to proceed, we first establish very high-level conditions on the behavior of a learned metric in order to satisfy the consistency of matching estimators. Then, we utilize feed-forward neural networks (NN) and siamese neural networks (SNN) trained (separately) on outcomes and treatment labels and extract the semi-final layer of these networks as low-dimensional representations that are, nonetheless, sufficient to control for the confounding effects of our covariates. The dimension of this embedding layer can be easily altered to trade off a richer representation for concerns about non-parametric bias. Through simulation and comparison to two canonical experimental datasets from the applied social sciences, we show our methods perform better than available competitors.
Following the potential outcome framework, we let $T \in \{0,1\}$ denote a unit's treatment status and $Y(1), \ Y(0)$ denote the outcomes that would be realized for each unit if it were treated (or not). Additionally, we adopt the following conventional assumptions.
Implicit in Assumption (ref) is the Stable Unit Treatment Value Assumption (SUTVA) that rules out any contaminating effects of one unit's treatment onto another unit's outcomes.
The average treatment effect on the treated (ATT) or untreated (ATUT) populations are then estimated under nearest neighbor matching (NNM) as
where $j^*_d(i) = \operatorname*{arg\,min}_{k:T_i \neq T_k} d(X_i,X_k)$ gives the nearest neighbor with opposing treatment status. We write the $d$ subscript to emphasize the estimator's dependence on some (potentially learned) distance metric in our setting. Now we establish high-level conditions that $d$ must obey in order to satisfy consistency of the resulting estimator.
Theorem (ref) guarantees consistency for metrics that suitably penalize either deviations in propensity score or a specific prognostic score (depending on the estimand). It is immediate to see that it can be applied to the already well-established results for MDM matching as well as propensity score and prognostic score matching for the cases when those models are correctly specified. Furthermore, it demonstrates the consistency of methods which (given guarantees on convergence) learn a matching metric that utilizes outcomes as labels, provided that only outcomes for control units are used to learn a metric for estimating the ATT, or conversely, that only outcomes for treated units are used to learn a metric for estimating the ATUT.\footnote{However be combined with additional results about convergence of the underlying metric learning methods in order to insure consistency in practice.} This closely parallels the results in Antonelli2016 that require two separate prognostic scores for average treatment effect estimation.
Critically, it does not support the consistency of matching estimators that utilize metrics learned on the pooled samples of outcome labels, including the procedure detailed in li2017matching. To demonstrate the potential failure of consistency in such settings, we pose a simple example with just a single covariate ($x\sim U[-1,1]$) and outcomes given by
A metric learned to minimize mean squared error on the pooled (treatment and control) outcome data would converge to
That is, for a treated unit with $x_i=\frac{1}{2}$ we would be indifferent between matching to a control unit with $x=\frac{1}{2}$, implying $\widehat{\tau}_i=\frac{1}{2}$, and matching to a control unit with $x=-\frac{1}{2}$, implying $\widehat{\tau}_i=\frac{3}{2}$. It is immediately evident that such an estimator would be inconsistent for both the ATT and ATUT parameters. The example in (ref) is highly stylized, but as long as one supposes heterogeneous treatment effects, it will be straightforward to construct counter-examples to the consistency of matching metrics learned on pooled outcome data.
In addition to using a control-only or treatment-only subsample of training units, an alternative but somewhat more parametric attempt to mitigate potential bias resulting from heterogenous treatment effects is detailed in JohanssonEtAl_icml16, wherein treatment status is concatenated with the learned data representation only at the stage immediately before prediction.
In this section, we provide a description of a general procedure to pre-process data for matching estimation. This procedure is then combined with NNM to produce a full matching method. We describe in additional detail how this procedure can be applied to the pseudo-metric learning paradigm specifically.
Our general framework, which is conformable to a variety of pre-processing techniques, involves training a predictive algorithm to predict an outcome $Y$, as well as training a separate instance of the algorithm to predict treatment status $D$. We then combine the learned representation of the data from these two tasks--where the method of combination depends on the predictive algorithm used--and carry out nearest-neighbor matching on the resulting object. Computationally, the complexity of the matching process is dominated by the component ML models. As these are standard model types their complexity is well established in the literature.
The following methods employ different neural net architectures in order to discover some low-dimensional representation of the original covariates on which a standard distance metric (e.g. Euclidean) can be applied to, together, constitute a learned pseudo-metric. In analogy to the alternatives of propensity and prognostic score matching, we learn two separate low-dimensional representations of the data, one in which units with similar potential outcomes are nearby each other,\footnote{When training the neural nets which learn similarity in outcomes, we use control-only or treatment-only subsamples. This decision is motivated by the previous discussion in Section (ref) regarding heterogeneous treatment effects.} and one in which units with similar propensity for treatment are nearby each other.
For a given unit $i$, these separate representations are characterized by two vectors of continuous-valued features, $\mathbf{m_{i,y}}$ and $\mathbf{m_{i,d}}$. We propose performing matching on the space defined by the union of these two representations, $M_{union} = [M_{y} \ \ M_{d}]$, where $dim(M_{d}) = n \times z_d$ and $dim(M_{y}) = n \times z_y$. It is possible that $z_d = z_y$ if the researcher so chooses, but this need not be true. Note that unlike methods which aim to select some subset of the original covariates, the combination of the learned features discussed here will include all learned features resulting from estimating the separate treatment and outcome models. Thus, the researcher, through specifying the dimension of the hidden layers to be extracted and used in matching, can explicitly control the dimension of the matching space. As a practical matter, we drop learned features which are perfectly correlated with another learned feature and learned features with near-zero variance. We recommend standard hyperparameter tuning for both the number of hidden layers as well as their size.
NNs are, theoretically, universal function approximators Cybenko1989. The success of NNs derives from the capacity of a network's inner layers to learn successive transformations of the data. Extracting these inner layers and using them as input into a second model is a known technique for dimensionality reduction Hinton504.
Specifically, we will refer to the final hidden layers of the network that predict treatment and outcome as $M_{d}$ and $M_{y}$, respectively. These are of dimension $z_d$ and $z_y$ respectively--parameters that the researcher can control when designing the network architecture that will determine the dimension of the matching space. Additionally, denote the vectors of final hidden layer weights as $\sigma_{l}^{d}$ and $\sigma_{l}^y$, respectively, and let $e_{n} = [1 \ 1 \ 1 \ . . . \ 1]$ be a $(1 \times n)$ vectors of ones. Then, letting $\otimes$ be the Kronecker product operator, define $A$, an $(n \times z_d)$ matrix, and $B$, an $(n \times z_y)$ matrix as
Using the above definitions, and letting the operator $\odot$ refer to element-wise matrix multiplication, we propose matching on
The matrix $X_{s}$ can be thought of simply as the matrix of scaled learned features. If a feature in the final hidden layer receives a very large weight as it enters the output layer, it is desirable to also increase the scale of that feature in the downstream matching estimation. Similarly, if a feature in the final hidden layer receives a near-zero weight as it is passed to the output layer, we would like to assign less importance to matching closely on that feature than on others.
The siamese neural network structure is not designed to predict a value, but rather to learn a mapping to a low-dimensional manifold such that alike observations are nearby each other and dissimilar observations are far from each other on the manifold. Each pass through the network involves two observations, $X_{i}$ and $X_{j}$. Each observation is passed through identical networks, terminating in a layer whose length the researcher controls.
The distance between the two layers is then input into one of two loss functions. The SNN which is designed to learn similarity in outcomes uses the loss function described in Equation (2.2). The SNN which is designed to learn similarity in treatment propensity uses a standard contrastive loss function as in Equation (2.3).
We recommend extracting the final hidden layer of each SNN, and matching on the combination of these layers. Denote the layer extracted from the treatment-target SNN as $M_{d}$ and the layer extracted from the outcome-target SNN as $M_{y}$. In the SNN architecture there is no set of weights connecting the extracted layer to some terminal node, instead the extracted layer enters directly into the loss function (2.2) or (2.3). Thus, unlike in the NN Method, there is no need to weight the extracted layer before its input into the matching estimator. We propose matching on
The rows of the matrix $X_{s}$ can be seen as representing the position of each observation $i$ on the learned manifold, where the mapping to that manifold has been explicitly constructed such that alike pairs are nearby each other in terms of Euclidean distance. This structure is naturally amenable to a standard matching procedure which assumes that observations that are nearby each other in terms of of Euclidean distance are alike. Indeed, this transformation of the raw covariate space renders that implicit assumption in the standard matching procedure far more plausible.
In order to benchmark the performance of the NN and SNN methods proposed here, we compare them to a number of other methods described in Table (ref)--including variable selection methods that choose a subset of the raw features---across three different simulated environments characterized by DGPs of increasing complexity.
For all simulated environments, the following are true:
Where $X_{i}$ denotes the vector of covariates for observation $i$, $K$ is the dimension of the covariate space, $N$ is the number of observations in a given simulation, and $\beta_0$ is the true treatment effect. We examine performance under three different DGPs: a sparse linear DGP, a sparse linear DGP with quadratic terms, and a NN DGP. Complete descriptions of the simulated environment can be found in Supplement A.
As additional benchmarks, we include two methods which generate a matching space through variable selection. First, in accordance with Section (ref), we estimate a pair of LASSO regressions to choose matching features. In the second alternative, we estimate a pair of Regularized Random Forests (RRF) Deng2012 to choose matching features on the basis of variable importance. Supplement B contains a more complete treatment of these alternatives.
Examining Table (ref) it is evident that, when the DGP is relatively simple, both PSM and PSMSQ do quite well in RMSE terms. Importantly, the two neural-net-based methods we propose, NN and SNN, perform nearly identically to PSM and PSMSQ in this setting.
At a first pass, it may seem counter-intuitive that oracle estimators--endowed with perfect knowledge of the relevant covariates--will underperform supervised algorithms which seek to recover some form of this knowledge. However, the non-parametric bias induced by matching to a nearest neighbor is $O(n^{-1/K})$ for problems with $K$ continuous covariates abadie2006large. This explains the poor performance of the oracle estimators and further motivates the pseudo-metric learning paradigm.
In Table (ref) we see that, as the DGP grows slightly non-linear, the performance of PSM and PSMSQ begins to deteriorate. It should be expected that PSM, now explicitly misspecified, should begin to falter. However, PSMSQ sees a non-trivial uptick in RMSE as well.
The simulation results from the NN DGP in Table (ref) indicate that, as the DGP grows increasingly complex, PSM and PSMSQ become wholly unusable. SNN and NN exhibit a very gentle uptick in RMSE, but are still nearly unbiased. If the researcher takes seriously the possibility that the DGP may be highly non-linear, then it is prudent to employ something akin to NN or SNN--it is harmless if the DGP turns out to be simple, but provides significant gains if is not.
We first examine the performance of the methods proposed in Section (ref) on the canonical lalonde1986evaluating \href{http://users.nber.org/ rdehejia/data/nswdata2.html}{data set}. Following imbens2014, we look specifically at the experimental and non-experimental versions of the original data set introduced in Dehejia1999.
In the LaLonde data, applicants to National Supported Work (NSW), a labor market training program, were selected at random for participation. Using the this data alone, we can learn the experimental ATT. Combining these data with a non-experimental comparison group derived from the Panel Study of Income Dynamics (PSID), it is possible to benchmark the performance of various causal inference methods designed to recover the ATT from observational data.
The variable of interest is labor market earnings in 1978, and available features include earnings in 1974 and 1975 as well as various demographic variables. In Table (ref), we compare non-experimental results with the experimentally estimated ATT of \$1,794. The NN and SNN estimates, standard errors, and confidence intervals are taken from the average over 100 iterations of each method. In the case of the non-experimental LaLonde data, the control and treated groups are significantly different in the observed features. For that reason, we find find it prudent to average over different iterations of the neural-net based methods in a situation with such extreme baseline control-treated feature imbalance.
As is standard practice in cases of significant class imbalance, we recommend oversampling in NN and SNN training if the treated units comprise less than $10\%$ of the entire sample. Other methods of protecting against overfit may also be used as well.
Standard errors are calculated according to abadie2006large for all estimators besides OLS. Among all methods, SNN comes the closest to recovering the experimental treatment effect, with NN outperforming all non-SNN methods by a considerable margin as well. While OLS provides the lowest standard error, it is considerably biased.\footnote{We attempted to benchmark our results on this data set against those reported in li2017matching, but were unable to recreate their exact environment.}
We additionally examine the performance of the methods proposed in Section (ref) on a modified version of the IHDP (Infant Health and Development Program) \href{https://www.icpsr.umich.edu/icpsrweb/HMCA/studies/9795}{data set}. The IHDP data set originates from a randomized longitudinal trial designed to discover the effect of a comprehensive set of interventions designed to curb health and developmental issues in low birth weight, premature infants.
In order to simulate an observational study, we follow hill2011, li2017matching, and others in removing all children with non-white mothers from the treatment group. The subset of removed infants is clearly non-random and therefore will induce the type of structural feature imbalance that we typically expect in non-experimental data.
Because the theoretical framework we present maintains the unconfoundedness assumption, we simulate outcomes using only pretreatment features and treatment assignment. The response surfaces we use follow li2017matching exactly. A full description can be found in Supplement C.
We simulate 50 response surfaces and report the results in Table (ref). Both NN and SNN result in small bias and RMSE, providing further support for the efficacy of these pseudo-metric learning methods across a variety of settings. While li2017matching only report the absolute value of the error in average treatment effect, they report an error of $0.16$, which is bested by both the NN and SNN methods. \newline
When estimating effects of conditionally exogenous treatment, pre-processing data with matching is a popular tool, and researchers with high-dimensional data may face a difficult decision about how best to define an appropriate metric. Following li2017matching, we recast the problem of matching in a metric learning framework that maps our covariates to a low-dimensional space. This facilitates “closer matches” while still capturing important aspects of unit-level heterogeneity.
We provide general conditions under which a learned distance metric is sure to lead to consistent estimation. This illuminates the importance of using control-only or treatment-only data when using supervised learning to discover a metric from outcome data. We also provide applied researchers with two new methods that leverage both MLPs and siamese neural nets and which compare favorably to state-of-the-art alternatives.
\nocite{abadie2016matching}