Causal Inference by Minimizing the Dual Norm of Bias: Kernel Matching & Weighting Estimators for Causal Effects Nathan Kallus 0 School of Operations Research and Information Engineering Cornell University and Cornell Tech New York , NY 10011 , USA

-

We consider the problem of estimating causal effects from observational data and propose a novel framework for matching- and weighting-based causal estimators. The framework is based on expressing the bias of a causal estimator as an operator on the unknown conditional expectation function of outcomes and formulating the dual norm of the bias as the norm of this operator with respect to a function space that represents the potential structure for outcomes. We give the term worst-case bias minimizing (WCBM) to estimators that minimize this quantity for some function space and show that a great variety of existing causal estimators belong to this family, including one-to-one matching (with or without replacement), coarsened exact matching, and mean-matched sampling. We propose a range of new, kernel-based matching and weighting estimators that arise when one minimizes the dual norm of the bias with respect to a reproducing kernel Hilbert space. Depending on the case, these estimators can be solved either in closed form, using quadratic optimization, or using integer optimization. In numerical experiments, the new, kernel-based estimators outperform all standard causal estimators in estimation error, providing a successful balance between generality and efficiency. they become comparable. Comparable for the purpose of causal inference means as similar as possible in some observed covariates. The covariates constitute the relevant information known about each observational subject and, as long as these covariates account for any confounding between the effects of treatment and the effects of selfselection, making the groups comparable with respect to these makes the groups comparable for the purpose of causal inference.

Matching and weighting have been some of the most popular ways to achieve this comparability [ 4, 29, 40 ]. In matching, we sample a (multi-)subset from the groups to get samples that are more similar to one another than the original samples. For example, in one-to-one matching [ 28 ], one composes a matched sample out of pairs of treated and control subjects so that the total pairwise distance between covariate vectors is small or even minimal, mimicking a randomized matched-pair experiment [ 12 ]. If we allow subjects to be paired with replacement, we can have a sample with duplicates. Weighting is a generalization where we can assign weights that are not integer multiples. For example, in coarsened exact matching (CEM) [ 17 ], one coarsens the covariates to create strata and re-weights the samples so that they have equal frequency in each stratum, mimicking a randomized block experiment [ 7 ].

Matching and weighting is employed for two purposes: (a) to reduce error due to confounding and (b) to reduce error due to imbalance. In a controlled experiment, (a) is achieved by randomization and (b) is achieved by, e.g., blocking (see [ 19 ] for more about balance in controlled experiments). For example, in an experiment on the effect of a drug on mortality, a randomized block design that ensures balance in important covariates such as age is generally more powerful than a completely randomized design, while both are unbiased (unconfounded). In observational studies, a very popular matching method to achieve (a) is propensity score matching (PSM) [ 30 ]. However, because it does nothing toward purpose (b) and because it depends on strong modelling assumptions to fit a correct propensity model, it is sometimes advocated that one match on the covariates themselves rather than estimated propensity scores [ 20 ]. A related weighting method is propensity score weighting (PSW) [ 13 ]. However, it too does nothing toward (b) and relies on strong modelling assumptions. Even if the model is correct, the estimated weights can be very unstable leading to a practice of ad hoc trimming, which may re-introduce confounding bias [ 6, 5 ]. Here, partly for these reasons, we focus on matching and weighting that address both (a) and (b) by balancing the covariates themselves rather than imputed estimates of propensities. In this paper, we develop a novel and encompassing framework for estimators that balance the covariates via weighting and matching. There are many different such estimators and each addresses imbalance differently. Our framework teases out how a particular notion of imbalance corresponds to a notion of structure. By decomposing the error of matching and weighting estimators, we formulate the bias of the estimator as an operator on the conditional expectation of outcomes given covariates. This conditional expectation function is unknown (or else there would be no need to experiment) and when one considers what the worst-case bias may be over a space of possible such functions one recovers the dual norm of the bias if the space is a Banach space. The dual norm of the bias is an observable quantity, expressed only in terms of the given data. We term any estimator that chooses matched subsamples or weights by minimizing the worst-case bias as worst-case bias minimizing (WCBM). A surprising result is that a great variety of standard methods used in the practice of causal inference are all WCBM. This observation leads us to consider new methods that are WCBM. Using reproducing kernel Hilbert spaces (RKHS) to express structure we obtain a new class of kernel-based matching and weighting causal estimators. 2

Set up We begin by describing the set up. We consider an observational study with n subjects. We index the subjects by i = 1, . . . , n. We let this order be arbitrary so that the subjects are exchangeable. Of these, n1 received a treatment whose effect is of interest (denoted by Ti = 1) and n0 received a control treatment against which we want to compare (denoted by Ti = 0). Let T0 = {i : Ti = 0} and T1 = {i : Ti = 1} be the sets of subjects that received treatment and control, respectively. We let T = (T1, . . . , Tn) denote the collection of treatment assignments, which constitutes part of the observed data.

Using Neyman-Rubin potential outcome notation [ 36 ], we let Yi(0), Yi(1) be the (real-valued) potential outcomes for subject i. We observe the outcome for the treatment to which subject i was exposed, Yi = Yi(Ti). And, Y (1 Ti) represents the unobserved, counterfactual outcome we would have observed if subject i were exposed to the opposite treatment. Y (1 Ti) is missing data. Throughout the paper, for these to be well defined, we assume that the stable unit treatment value assumption (SUTVA) holds [ 32 ], which requires that which treatment one of subjects experiences not affect the outcomes of another subject and that potential outcomes are fixed as which treatment is experienced changes (so only which one we observe, and hence Yi, is affected).

Let Xi, taking values in some X , be the side covariates that we observe for subject i. Let X = (X1, . . . , Xn) denote the collection of all baseline covariates of all n subjects, which constitues part of the observed data. The space X is general; assumptions about it will be specified as necessary. As an example, it can be composed real-valued vectors X ✓ Rd that include both discrete (dummy) and continuous variables.

We denote by TEi = Yi(1) Yi(0) the unobservable causal treatment effect for subject i. The primary quantity of interest for estimation is the sample average (causal) treatment effect on the treated sample (SATT): SATT = n11 Pi2T 1 TEi = n11 Pin=1 Ti(Yi(1) Yi(0)).

We consider estimators for SATT based on weighting. We restrict to honest weights that only depend on the observed X, T and not on any observed outcome data. (If we used outcome data one might complain that we are mining for an effect that is not there.) In particular, we will consider the choice of a weighting function W = W (X, T ) that produces a weight Wi 2 R for each subject i, leading to the estimator ⌧ˆW = Pn

i=1( 1)Ti+1WiYi.

Because we are estimating SATT and we in fact observe Yi(1) for each i 2 T 1, we always set Wi = 1/n1 for i 2 T 1, leading to estimators of the form ⌧ˆW = n11 Pi2T 1 Yi

Pi2T 0 WiYi.

We also always assume Pi2T 0 Wi = 1.

The bias of the estimator resulting from weights W is the difference between it and SATT, conditioned on all the observable data upon which the weights are based: bias = E [⌧ˆW

SATT | X, T ] .

We let W = W0 ⇥ W 1 denote the space of allowable weights, where W0 and W1 are the space of weights for the control and treated sample, respectively. We required that W0 ✓ { WT0 2 RT0 : Pi2T 0 Wi = 1} and that W1 = {(1/n1, . . . , 1/n1)}. If all weights in W0 are rational with a fixed denominator, we call ⌧ˆW a matching estimator because it is equivalent to constructing a (multi-)set from the control subjects to match the treated sample. We note some special cases of W0 that correspond to a variety of existing classes of estimators for SATT: • General weights:

W0general = {WT0 2 RT0 : Pi2T 0 Wi = 1}. • Nonnegative (probability) weighting:

W0nonnegative = {WT0 2 RT+0 : Pi2T 0 Wi = 1}. • Matching with fixed size n00 and without replacement:

W0w/o rep. = {WT0 2 { 0, 1/n00}T0 : Pi2T 0 Wi = 1}. • Matching with fixed size n00 and with replacement:

W0w/ rep. = {WT0 2 { 0, 1/n00, . . . }T0 : Pi2T 0 Wi = 1}. Note that, as estimators for a population effect (if we are to assume random sampling of subjects from a population), both SATT and ⌧ˆ preclude regression adjustments [ 9 ]. Nonetheless, without parametric assumptions necessary for such adjustments, assuming random sampling, SATT is the uniform minimum variance unbiased estimator for the treatment effect on the treated population [ 22 ]. In fact, matching and weighting are often regarded as means of reducing model dependence [ 14 ]. Specifically, if one matches very closely, then simple difference estimators and complicated regression estimators are all very close, so the point of which estimator to use after matching is largely moot.

A standing assumption in this paper, essential for causal inference from observational data, is that of weak ignorability in expectation.

Assumption 1. For each t = 0, 1 and i = 1, . . . , n, conditioned on Xi, Yi(t) is mean-independent of Ti and each value of Ti is possible. That is, for each t = 0, 1 and i = 1, . . . , n,

E [Yi(t) | Ti, Xi] = E [Yi(t) | Xi] , and

P Ti = t Xi > 0.

Ignorability, also known as unconfoundedness, means that we have the right covariates needed to separate the effect of the treatment itself from the effect of self-selection. For example, in an observational study, self-selection for “treatment” might imply affluence, which might imply good outcomes regardless of treatment, so if we control for income in Xi we can isolate this effect. The form of ignorability we use is termed “weak” because it need only apply for each t = 0, 1 separately, and it is termed “in expectation” because only mean-independence, rather than full stochastic independence, is assumed.

Aside: alternative frameworks for causal inference

The Neyman-Rubin potential outcome framework is not the only framework used to describe causal relationships. Other frameworks for causality include, most notably, Pearl’s framework of causal Bayesian networks and docalculus [ 23 ] as well as structural equation models (SEM) [ 10 ]. We do not consider the SEM framework because of its need for a priori models, the common restriction to linear relationships, incompatible notation, and the less clear question of model-free identifiability.

Pearl’s framework generalizes both potential outcomes and SEM [ 26, 25 ]. Inference in this framework depends on directed acyclic graph (DAG) models to describe a priori causal relationships. The standard practice in applications of the Neyman-Rubin framework is generally to condition on all observed covariates X that are potentially relevant [ 33 ], but one can easily come up with DAG constructions where the inclusion of a covariate in such conditioning can (asymptotically) bias causal estimates, the simplest of which is the M -graph [ 37, 24 ]. In effect, a causal DAG, correctly specified, can specify the correct subset of the covariates X that should be included in order to achieve Assumption 1. The estimation or validation of a causal DAG from data is an active field of research, e.g. [ 15, 38, 27 ]. 4

The worst-case bias

We define the conditional expectation of the control potential outcome given the covariates x as follows:

f0(x) = E ⇥Yi(0) Xi = x⇤ .

The non-random function f0 does not depend on i due to exchangeability. By the law of iterated expectation, the residual ✏i = Yi(0) f0(Xi) has mean 0, is meanindependent of Xi, and is uncorrelated with any function of Xi.

By conditioning on Xi, we can decompose the error of the estimator into two terms: error that can be controlled by matching on Xi and the orthogonal residual error, which cannot be controlled by Xi but which disappears in expectation due to ignorability.

Theorem 1. Under Assumption 1, the bias of ⌧ˆW is

E [⌧ˆW where

SATT | X, T ] = B(W ; f0), B(W ; f ) := n11 Pi2T 1 f (Xi)

 Moreover, letting we have that

E(W ) = n11 Pi2T 1 ✏i Pi2T 0 Wi✏i, ⌧ˆW

SATT = B(W ; f0) + E(W ), E [E(W ) | X, T ] = 0.

Pi2T 0 Wif (Xi).

 Proof of Theorem 1. Let us write SATT as

SATT = n11 Pi2T 1 Yi n11 Pi2T 1 Yi(0).

It is then clear that SATT differs from ⌧ˆW only in the second term, that is, ⌧ˆ SATT = n11 Pi2T 1 Yi(0) = Pn i=1( 1)Ti+1WiYi(0)

Pi2T 0 WiYi(0) = Pn i=1( 1)Ti+1Wi✏i, i=1( 1)Ti+1Wif0(Xi) + Pn where we recognize the last term as E(W ). For each term of E(W ) we have

E ⇥( 1)Ti+1Wi✏i X, T ⇤ = ( 1)Ti+1Wi (E [Yi(0)|X, T ] = ( 1)Ti+1Wi (E [Yi(0)|X] f0(Xi)) f0(Xi)) = 0, where the first equality is by definition of ✏i and the fact that Wi = Wi(X, T ) and the second is by Assumption 1.

The target of weighting or matching for causal inference is to eliminate bias in comparing the treatment and control samples. Theorem 1 provides an explicit form of the bias in terms of the observed covariates X. However, it involves the unknown function f0 : X ! R. As alluded to in Sec. 1, we consider weighting schemes that guard against any possible such function by minimizing the worst-case bias over the unit ball of a Banach space. A normed vector space is a Banach space if the corresponding metric space is complete (see [ 21 ] and Ch. 10 of [ 31 ] for more on Banach spaces). Let V denote the vector space of all functions X ! R under usual pointwise addition and scaling. Let F ✓ V be a subspace of functions, against which we wish to guard. Endow this space with a semi-norm k·k : F ! R (a seminorm can assign zero magnitude to nonzero vectors). For f 2 /F , let us write kf k = 1 . Thus, the assumption that f0 2 F is encapsulated by kf0k < 1 .

Given only that kf0k < 1 , we will consider weighting or matching schemes that choose W to minimize the worstcase bias,

maxkfkk f0k |B(W ; f )| = kf0k maxkfk 1 B(W ; f ), where the equality holds because B(W ; ↵f ) = ↵B (W ; f ) is degree-1 homogeneous and k↵f k = |↵ | kf k is degree1 positively homogeneous and symmetric. Clearly, it only matters that kf0k < 1 and the particular finite value of it does not change which W minimizes the above. In light of this, we define the worst-case bias as

B(W ; F ) = maxkfk 1 B(W ; f ). c), where c 2 R represents a constant function x 7! c. To eliminate this irrelevant mode of F , we can just consider the quotient space F /R, which consists of the equivalence classes [f ] = {f + c : c 2 R} endowed with the norm k[f ]k = minc2 R kf + ck. Note that by construction, B(W ; [f ]) = B(W, f ) is well defined. For brevity, we will simply refer to F and k·k when we mean F /R and the corresponding norm.

We will consider spaces (F , k·k) that satisfy the following conditions: Assumption 2. The space F is a Banach space. Assumption 3. For each W 2 W , f 7! continuous mapping F ! R.

B(W ; f ) is a Since B(W, f ) is also linear in f , these assumptions imply that, for each W , the operator B(W, ·) is in the continuous dual space of F . Hence,

B(W ; F ) = kB(W ; ·)k⇤ is precisely the dual norm of the bias, where the dual norm of a continuous linear operator A on a Banach space with norm k·k is kAk⇤ = supkuk 1 A(u). This also guarantees that B(W ; F ) is finite and well-defined.

Definition 1. A weighting (or matching) method W (T, X) is said to be worst-case bias minimizing (WCBM) if for some W and (F , k·k) satisfying Assumptions 2 and 3 we have

W (T, X) 2 arg minW 2W

B(W ; F ) 6= W.

Let Bmin(F ) = minW 2W B(W ; F ) be the optimal value. Clearly, if a weighting method W (T, X) is WCBM with (F , k·k) and W then the bias of ⌧ˆW is bounded by |B(W ; f0)|  k f0k Bmin(F ). 5

 Existing methods as WCBM

Surprisingly, a great many methods for causal inference that are standard in practice are also in fact WCBM. On the one hand, this interpretation gets at the core of the structural motivations behind many of these methods (e.g., “if you believe the conditional expectation is Lipschitz and nothing more then you should pairwise match”) and allows one to choose a method appropriate to one’s beliefs about problem structure. On the other hand, these results provide motivation that WCBM is the right framework in which to think about weighting and matching for causal inference and this motivates us to consider new WCBM methods in Sec. 6. 5.1

One-to-one matching Since Pn

i=1( 1)Ti+1Wi = 0, we have that B(W ; f ) is invariant to constant shifts to f , i.e., B(W ; f ) = B(W ; f + One-to-one (pairwise) matching is by far the most common matching method. In one-to-one matching, each treated subject is paired with exactly one control subject so that the sum of pairwise distances is minimized as measured by some distance metric (x, x0) on X [ 28 ]. Usually, the Mahalanobis metric is used: (x, x0) = q (x x0)⌃ ˆ 1(x x0, where ⌃ ˆ is the pooled sample covariance matrix. Oneto-one matching can be done either without replacement (each control subject used at most once) or with replacement (each control subject could be reused and matched to two or more treated subjects). The estimate of SATT is the average pairwise differences of outcomes. This estimator is exactly ⌧ˆW where the weight on control subject i is 1/n1 times the number of times subject i was matched, i.e., the matched control sample is the (multi-)set of control subjects that got matched to treated subjects.

One-to-one matching is WCBM.

Theorem 2. One-to-one matching with pairwise distance metric (x, x0) with replacement and without replacement are both WCBM with • kf k = supx6=x0 f(x()x,fx0()x0) , the Lipschitz constant of f ; • F = {f : kf k < 1} ; • W0 is either W0nonnegative or W0w/ rep. (with n00 = n1) if with replacement; and • W0 is either n

WT0 2 W 0nonnegative : n1Wi  1 8 io or

W0w/o rep. (with n00 = n1) if without replacement. Remark 1. Note that even if the weights are not restricted to be multiples of 1/n1, the optimal unrestricted weights will end up to be multiples of 1/n1 regardless. That is, the optimal weighting is optimal matching for Lipschitz functions.

Remark 2. Note that (F , k·k) is not a Banach space. In particular, constant functions have zero Lipschitz constant. However, as required, F /R is a Banach space and evaluation differences are continuous because they are bounded by the magnitude.

Remark 3. Algorithmically, one-to-one matching with replacement amounts to finding the control subject of minimal distance to each treated subject in a greedy manner. One-ton-one matching without replacement amounts to minimum-sum-of-distances bipartite matching with unbalanced parts, which is easily solved by the Ford-Fulkerson algorithm [ 8 ].

Proof of Theorem 2. Let D be the distance matrix Dii0 = (Xi, Xi0 ). For this choice of (F , k·k), by linear optimizaWi0 2 RT+0 P Pi2T 0 Wi0 = n1 Pi02T 0 Sii0 = 1 8 i 2 T 1

i2T 1 Sii0 Wi0 = 0 8 i0 2 T 0.

This describes a min-cost netwrok flow problem with sources T1 with inputs 1; nodes T0 with 0 exogenous flow; one sink with output n1; edges from each i 2 T 1 to each i0 2 T 0 with flow variable Sii0 , cost Dii0 , and without capacity; and edges from each i 2 T 0 to the sink with flow variable Wi0 and without cost or capacity. Because all data is integer, the optimal solution of W 0 = n1W is integer (see [ 1 ]). Hence, since W0w/ rep. ✓ Z/n1, the solution is the same when we restrict to W0 = W0w/ rep.. This solution (in terms of W 0) is equal to sending the whole input 1 from each source in T1 to the node in T0 with smallest distance and from there routing this flow to the sink, which corresponds exactly to one-to-one matching with replacement. In the case of no replacement and for W0 = {W 2 W0nonnegative : n1Wi  1 8 i}, using the transformation tion duality we get

sup P vi B(W ; F ) = n11 vi vi0  Dii0 8 i,i0 i2T 1 P n1Wivi i2T 0 ! = n11 minS s.t. S 2

Pn

i0=1 (Sii0 Pn

i0=1 (Sii0 Pi,i0 Dii0 Sii0

R+⇥ n n

Si0i) = 1 Si0i) =

8 i 2 T 1 n1Wi 8 i 2 T 0.

This describes a min-cost network flow problem with sources T1 with inputs 1, sinks T0 with outputs Wi, edges between every two nodes with costs Dii0 and without capacities. Consider any source i 2 T 1 and any sink i0 2 T 0 and any path i, i1, . . . , im, i0. By the triangle inequality, Dii0  Dii1 + Di1i2 + · · · + Dimi0 . Therefore, as there are no capacities, it is always preferable to send the flow from the sources to the sinks along the direct edges from T1 to T0. That is, we can eliminate all other edges and write B(W ; F ) = n11 minS

P

i2T 1, i02T 0 Dii0 Sii0 s.t. S 2 RT+1⇥T 0

P Pi02T 0 Sii0 = 1

i2T 1 Sii0 = n1Wi In the case of with replacement and W0 = using the transformation Wi0 = n1Wi, we get

This describes the same min-cost netwrok flow problem except that the edges from each i 2 T 0 to the sink have a capacity of 1. Because all data is integer, the optimal solution of S and W 0 = n1W is integer (see [ 1 ]). Hence, since W0w/o rep. ✓ Z/n1, the solution is the same when we restrict to W0 = W0w/o rep.. The optimal Sii0 is integer and so, by Pi02T 0 Sii0 = 1, for each i 2 T 1 there is exactly one i0 2 T 0 with Sii0 = 1 and all others are zero. Sii0 = 1 denotes matching i with i0. The optimal Wi0 is integral and so, by Wi0  1, Wi0 2 { 0, 1}. Hence, for each i 2 T 0, Pi02T 1 Sii0 2 { 0, 1} so we only use node i at most once. The cost of S is exactly the sum of pairwise distances in the match. Hence, the optimal solution corresponds exactly to one-to-one matching without replacement. CEM [ 17 ] is a weighting method whereby one coarsens the covariates into a few (M ) strata via a coarsening function C : X ! { 1, . . . , M }, and then one matches exactly within each stratum. For example, if there are 5 treated subjects and 3 control subjects in a given stratum then each of the 3 control subjects is given weight proportional to 5/3, whereas if there were 0 treated subject the weights would be 0. The case of a stratum containing only treated subjects is not allowed (no extrapolation). ([ 18 ] suggests that in this case one estimates the “feasible average treatment effect on the treated,” meaning to modify the sample of interest from the treated sample to the subset that has good matches.) Under Assumption 1, lack of any overlap is rare for large n.

Theorem 3. CEM with coarsening function C : X ! {1, . . . , M } is WCBM with • F = f : f 1(C 1(j)) = 1 8 j = 1, . . . , M , i.e., piece-wise constant on the coarsening partitions; • kf k = supx2X |f (x)| for f 2 F , otherwise 1 ; and • W0 is either W0general or W0nonnegative, assuming that each partition that contains a treatment subject also contains a control subject (no extrapolation). Very often, practitioners will evaluate the quality of a matched control sample by measuring the Mahalanobis distance between the matched control sample and the treated sample: MV (W ) =

V 1/2 ⇣ 1 P n1 i2T 1 Xi

Pi2T 0 WiXi ⌘ 2 , where X ✓ Rd and V is some positive definite matrix usually taken to be V = ⌃ ˆ , the pooled sample covariance matrix of X. This distance is a rotated 2-norm between the sample means. Mean-matched sampling finds a matched control sample of a prescribed size to minimize this distance.

Theorem 4. Mean-matched sampling n00 subjects (with or

without replacement) from the control set is WCBM with • F = x 7! 0 +

T x : 2 Rd ; • kx 7! 0 + otherwise; and

T xk = p T V + 02 and kf k = 1 • W0 is either W0w/ rep or W0w/o rep, respectively. Remark 4. Since finite, the space (F , k·k) is always a Banach space and evaluations (and hence their differences) are always continuous. See Thms. 5.33 and 5.35 of [ 16 ]. Proof of Theorem 4. By duality of norms,

B(W ; F ) = sup

T V  1

T 1 P Xi n1 i2T 1

P WiXi i2T 0 ! = MV (W ).

The optimal W minimizes this discrepancy over subsamples from control with the allowable size. 6

 Kernel WCBM methods

In the previous section we saw that a variety of standard methods for causal inference are WCBM. Each was recovered using a different form of structure on the conditional expectations of outcomes. In this section we develop a range of new WCBM based on kernels and their corresponding reproducing kernel Hilbert spaces (RKHS). Kernels are standard in machine learning (ML) as ways to generalize the structure of learned conditional expectation functions, like classifiers or regressors [ 34 ]. Kernels have many applications in statistics [ 2, 11, 41 ]. The same way kernels are used to generalize the structure of learned functions in ML, we can use these to generalize the structure of f0.

A Hilbert space is an inner-product space such that the norm induced by the inner product, kf k2 = hf, f i, yields (a) `1 norm (b) Quadratic (c) Cubic (d) Sinusoidal a Banach space. An RKHS F is a Hilbert space of functions for which, for every x 2 X , the map f 7! f (x) is a continuous mapping [ 2 ]. Continuity and the Riesz representation theorem imply that for each x 2 X there is K(x, ·) 2 F such that hK(x, ·), f (·)i = f (x) for every f 2 F . The symmetric map K : X ⇥ X ! R is called the reproducing kernel of F . The name is motivated by the fact that F = closure (span {K(x, ·) : x 2 X } ). Thus K fully characterizes F . Prominent examples of kernels for X ⇢ Rd are: (i) (ii)

The polynomial kernel Ks(x, x0) = (1 + xT x0/s)s, whose RKHS spans the finite-dimensional space of all polynomials of degree up to s.

The exponential kernel K(x, x0) = exT x0 , the infinite-dimensional limit of the polynomial kernel. (iii) The Gaussian kernel K(x, x0) = e kx x0k2 . The corresponding RKHS is infinite-dimensional [ 39 ]. For X 2 X n and a kernel K, the Gram matrix is Kij = K(Xi, Xj ), which is always positive semi-definite (PSD). Generally, we normalize the covariate data before putting it in a kernel so that the sample has zero sample mean and identity pooled sample covariance Note that any RKHS F satisfies Assumptions 2 and 3. As such it gives rise to WCBM matching and weighting methods.

Theorem 5. Let F be an RKHS with kernel K. Let K be the Gram matrix on X. Then,

⇣ B(W ; F ) =

WTT0 KT0T0 WT0 2k0T WT0 + k0T k0 ⌘1/2 , where k0 = KT0T1 en1 /n1.

Remark 5. If WT0 2 { 0, 1/n00}T0 then B(W ; F ) is exactly the kernel maximum mean discrepancy (MMD) statistic between the treated sample and the matched control sample. Kernel MMD is a common test statistic in twosample goodness-of-fit testing [ 11, 35 ]. We can interpret minimizing this discrepancy as trying to make the two samples appear to come from the exact same distribution.

Proof of Theorem 5. We have

B2(W ; F ) = maxkfk 1 Pin=1( 1)Ti+1Wif (Xi) = ⌦Pn i=1( 1)Ti+1WiK(Xi, ·)↵ i=1( 1)Ti+1WiK(Xi, ·), Pn

2 = Pn

i,j=1( 1)Ti+Tj WiWj Kij , which when written in block form gives rise to the result.

For F an RKHS, we refer to minimizing B(W ; F ) over nonnegative weights as kernel weighting. This can be formulated as a linearly-constrainted convex-quadratic optimization problem: argminWT0 2W 0nonnegative B(W ; F ) = argminW 2 Rn+0 :eTn0 W =1 W KT0T0 W

2k0T W .

This problem can be solved in polynomial time with interior point methods [ 3 ] and is amenable to solution with off-the-shelf solvers like Gurobi.

For F an RKHS, we refer to minimizing B(W ; F ) subsets or multisubsets as kernel matching. Kernel matching with replacement can be formulated as a linear-integer-constrainted convex-quadratic optimization problem: argminWT0 2W 0w/ rep. B(W ; F ) = n100 W 02 Zanr0 g:emTn0iWn 0=n00 ⇣ n10 W 0KT0T0 W 0 0

2k0T W 0⌘ , where we used the change of variables W 0 = n00WT0 . Kernel matching without replacement is similarly formulated: argminWT0 2W 0w/o rep. B(W ; F ) = n100 W 02{ 0,1a}rng0 m:eiTnn0 W 0=n00 ⇣ n10 W 0KT0T0 W 0 0

2k0T W 0⌘ .

Both problems are NP-hard (reducible to number partitioning for rank(KT0T0 ) = 1), but amenable to solution by off-the-shelf integer programming solvers like Gurobi.

RMSE 150 (a) `1 norm 0.50 ▽◇□ ◆●▲▼■ ◆●▲▼■ ◆●▲▼■ ◆●▲▼ ●▲▼ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ●▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ○△ ◇ ▽□ ◇ ◇ ◇■ ◆◇■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◆◇▼■ ◇▼■ ◆◇▼■ ◆◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇■ ◆●◇■ ◆●◇■ 00..1005 ○△ ◆ ▼ ▼ ▼ 0.10 ○▽△□ ▽□ ▽□ ▽ 0.05 ○△ ○△ ○△□ ○▽△□ ○▽△□ ○▽△□ ○△▽□ ○△▽□ △▽□ ○△▽□ △▽ ▽△ △▽ ○

△▽ △▽ △▽ ▽△ △▽ ▽△ ▽△ △▽ ▽△ ▽△ ▽△ △▽ △▽ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ ○□ 150 (c) Cubic 50 100 200 250 300 n ▽ 300 n 300 n 50 100

150 (b) Quadratic

200 RMSE 1 ○△◇▲ ▽ ▼ □ ◆ ■ ○△◇ ◇▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲

▲ 0.50 ◆▽●□▼■ ◆○△▽●▼□■ ◆○△●◇▼■ ◆●◇▼■ ◆△●◇▼ ◆●◇▼ ●◇▼ ●◇▼ ●◇▼ ▼ ▽ ○△ ○■ △■ ◆■ ◆■ ◆■ ◆●◇■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ●◇▼■ ◆●◇▼■ ●◇▼■ ●◇▼■ ◆●◇▼■ ●◇▼■ ●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ ◆●◇▼■ □ ○ ▽□ ▽□ ▽ ○△ ○△ ○△ ○△ ○△ □ ▽□ ▽□ ▽□ ▽□ ▽□ ○△ ○△ △ △ ◆ ◆ ◆ ◆ ◆ ○ ○ △ △

○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ ▽ ▽ ▽□ ▽□ ▽ □ □ ▽□ ▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□ ○▽□

□ 50 100 200

250 150 (d) Sinusoidal 300 n ● No matching ◆ CEM ▼ PSM ○ Exp kernel weight △ Exp kernel match ■ One-to-one ▲ Mahal. means ◇ Quad kernel weight □ Gauss kernel weight ▽ Gauss kernel match

 7 Numerical experiments

In this section, we study the comparative efficiency of various causal estimators, including our new kernel estimators. (a) `1 norm: f0(x) = |x1| + |x2|; Consider the following fictitious observational study with (b) quadratic: f0(x) = (x1 + x2) + (x1 + x2)2; one treatment and control. Subjects are drawn at random (c) cubic: f0(x) = (x1 + x2)2 + (x1 + x2)3; from a population. For each subject we observe a twodimensional vector of covariates Xi 2 R2. In the popula- (d) sin: f0(x) = sin(⇡ (x1 + x2)) + cos(⇡ (x1 x2)). tion, these are distributed as uniform on [ 1, 1 ]2. Each subject has either received treatment or control and we observe These are shown in Figure 1.

Ti. In the population, Ti is distributed as Bernoulli with probability 0.8/ 1 + p 2 kXik2 , which ranges 0.27 ⇠ FFoorr eeaacchh enxp=er1im0,e2n0t,w. e..c,o3n0s0id,ewreapvraordieutcyeo1f0e0stirmepaltiocrast:es. 0.8.

The potential outcomes are distributed as

Yi(0) = f0(Xi) + ✏0i, Yi(1) = f1(Xi) + ✏1i, cwuhseroen✏t0hie,✏c1aise⇠Nof s(m0,a0ll.1r)esiisduinadlenpoeinsdee(nvtanrioainscee. nWoet efxo-- (b) Omnale-btoip-aornteit:ewmeamtchaticnhg non1 cthoentmroaltrsiuxbojefcptsaiurswinisgeoMptai-plained by Xi) so to tease out the comparative efficiency in halanobis distances between treated and control submatching X (if residual noise is big, any method that only jects; (a) No matching: we take the whole control sample to be

the matched sample (Wi = 1/n0); matches on X will do badly). We let f1 be any function whatsoever. We consider a variety of possible cases for f0: (c) CEM: we find the largest b 0 such that coarsening each of the covariates into even bins {[ 1, 1 + 2b 1), . . . , [1 2b 1, 1]} leaves no box (product of two bins) that contains only treated subjects, then we perform exact matching within each box; (d) Mahal. means: we match n1 control subjects with replacement to minimize the Mahalanobis distance between the means of the two samples; (e) PSM: we match n1 control subjects using propensity score matching by fitting a logistic regression to impute propensity scores and doing optimal bipartite matching on imputed scores; (f) Quad kernel weight: we use nonnegative kernel

weighting with the quadratic kernel; (g) Exp kernel weight: we use nonnegative kernel weight

ing with the exponential kernel; (h) Gauss kernel weight: we use nonnegative kernel

weighting with the Gaussian kernel; (i) Exp kernel match: we match n1 control subjects with replacement using kernel matching with the exponential kernel; and (j)

Gauss kernel match: we match n1 control subjects with replacement using kernel matching with the

Gaussian kernel.

We use Gurobi v6.5 (www.gurobi.com) to solve all quadratic and integer optimization problems. For each estimator, we compute ⌧ˆW SATT. Then, we measure the RMSE over the 100 replicates, RMSE = (Eˆ100 h(⌧ˆW SATT)2i)1/2. We plot the results in Figure

2. Note the log scale.

The results clearly show the power of our approach. In each case, every one of our exponential- or Gaussian-kernelbased estimators outperforms standard causal estimators by an order of magnitude (base 10). The advantage is particularly noticeable in smaller samples and for our kernel weighting methods. This can be explained by the fact that it can be difficult to find a good control pair for every treated subject in small samples, and similarly it can be difficult to have a fine enough coarsening of the data without creating a stratum that only has treated subjects. At the same time, by optimizing the mismatch as characterized by the dual norm of the bias one can achieve small mismatch with even small samples.

Another observation is that matching based on parametric models can be fragile. This can be seen here for PSM, which is based on a misspecified logistic model, and also for estimators that match on X itself. We also see that mean-matched sampling does very poorly in every example, even doing worse than no matching. Indeed, matching the means only makes sense if the effect is purely linear. A linear model assumption is very fragile and even small violations can trip up mean-matched sampling. Similarly, matching per the quadratic kernel depends on an assumption of quadratic effect. Indeed, the estimator based on the quadratic kernel does the best of all estimators when the effect is quadratic (panel b). However, unlike linear, a quadratic model is generally more robust as quadratics can better approximate a wider range of functions. Accordingly, we see that the estimator based on the quadratic kernel has reasonable performance even when the effect is not quadratic (panels a and c), while extreme violations trip it up (panel d).

Overall, the universal kernels (exponential and Gaussian) seem to do the best by far. They appear to provide a good balance between generality of model with efficiency of balancing. And, fully optimizing mismatch as measured by the dual norm of the bias in their RKHS can lead to small objective value even for moderate n. 8

 Conclusion

We presented a novel framework for matching and weighting estimators for causal inference from observational data. The framework is based on minimizing the dual norm of the bias operator with respect to a space of possible conditional expectation functions. Many existing methods common in practice appear to fit this framework. The framework also gives rise to kernel-based estimators . that prove exceedingly successful in a numerical experiment. [41] K. Zhang, J. Peters, D. Janzing, and B. Scho¨lkopf.

Kernel-based conditional independence test and application in causal discovery. 2012.

 [1] R. Ahuja , T. Magnanti , and J. Orlin . Network flows: theory, algorithms, and applications . Prentice Hall, Upper Saddle River, 1993 . [2] A. Berlinet and C. Thomas-Agnan . Reproducing kernel Hilbert spaces in probability and statistics . Kluwer Academic, 2004 . [3] S. P. Boyd and L. Vandenberghe. Convex Optimization . Cambridge University Press, Cambridge, 2004 . [4] W. G. Cochran . Planning and analysis of observational studies . John Wiley & Sons, New York, 1983 . [5] R. K. Crump , V. J. Hotz , G. W. Imbens , and O. A. Mitnik . Dealing with limited overlap in estimation of average treatment effects . Biometrika , 2009 . [6] M. R. Elliott . Model averaging methods for weight trimming . Journal of official statistics , 24 ( 4 ): 517 , 2008 . [7] R. A. Fisher . Statistical methods for research workers . Oliver and Boyd , Edinburgh, 1925 . [8] L. R. Ford and D. R. Fulkerson . Maximal flow through a network . Canadian journal of Mathematics , 8 ( 3 ): 399 - 404 , 1956 . [9] D. A. Freedman . On regression adjustments to experimental data . Advances in Applied Mathematics , 40 ( 2 ): 180 - 193 , 2008 . [10] A. S. Goldberger . Structural equation methods in the social sciences . Econometrica , 1972 . [11] A. Gretton , K. M. Borgwardt , M. Rasch , B. Scho¨lkopf, and A. J. Smola . A kernel method for the two-sample-problem . In Advances in neural information processing systems , 2006 . [12] X. S. Gu and P. R. Rosenbaum . Comparison of multivariate matching methods: Structures, distances, and algorithms . Journal of Computational and Graphical Statistics , 2 ( 4 ): 405 - 420 , 1993 . [13] K. Hirano and G. W. Imbens . Estimation of causal effects using propensity score weighting: An application to data on right heart catheterization . Health Services and Outcomes research methodology , 2 ( 3 - 4): 259 - 278 , 2001 . [14] D. E. Ho , K. Imai , G. King , and E. A. Stuart . Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference . Political analysis , 15 ( 3 ): 199 - 236 , 2007 . [15] P. Hoyer , D. Janzing , J. Mooij , J. Peters , and B. Scho ¨lkopf. Nonlinear causal discovery with additive noise models . 2009 . [16] J. K. Hunter and B. Nachtergaele. Applied Analysis . World Scientific , 2001 . [17] S. M. Iacus , G. King , and G. Porro . Causal inference without balance checking: Coarsened exact matching . Political analysis, page mpr013 , 2011 . [18] S. M. Iacus , G. King , and G. Porro . A theory of statistical inference for matching methods in applied causal research . 2015 . [19] N. Kallus . Optimal a priori balance in the design of controlled experiments . 2014 . [20] G. King and R. Nielsen . Why propensity scores should not be used for matching . 2015 . [21] M. Ledoux and M. Talagrand . Probability in Banach Spaces: isoperimetry and processes . Springer, 1991 . [22] E. L. Lehmann and G. Casella . Theory of Point Estimation . Springer-Verlag, New York, 1998 . [23] J. Pearl . Causality: models, reasoning and inference . Cambridge University Press, 2000 . [24] J. Pearl . Remarks on the method of propensity score . Statistics in Medicine, 28 ( 9 ): 1415 - 1416 , 2009 . [25] J. Pearl . The causal foundation of structural equation modeling . In R. Hoyle, editor, Handbook of structural equation modeling , pages 68 - 91 . Sage, Newbury Park, 2012 . [26] J. Pearl et al. Causal inference in statistics: An overview. Statistics Surveys , 3 : 96 - 146 , 2009 . [27] J. Peters , J. Mooij , D. Janzing , and B. Scho ¨lkopf. Identifiability of causal graphs using functional models . 2012 . [28] P. R. Rosenbaum . Optimal matching for observational studies . Journal of the American Statistical Association , 84 ( 408 ): 1024 - 1032 , 1989 . [29] P. R. Rosenbaum . Observational studies. Springer, New York, 2002 . [30] P. R. Rosenbaum and D. B. Rubin . The central role of the propensity score in observational studies for causal effects . Biometrika , 70 ( 1 ): 41 - 55 , 1983 . [31] H. L. Royden. Real Analysis . Prentice Hall , 1988 . [32] D. B. Rubin . Comments on “randomization analysis of experimental data” . Journal of the American Statistical Association , 75 ( 371 ): 591 - 593 , 1980 . [33] D. B. Rubin . Should observational studies be designed to allow lack of balance in covariate distributions across treatment groups? Statistics in Medicine, 28 ( 9 ): 1420 - 1423 , 2009 . [34] B. Scholkopf and A. J. Smola . Learning with kernels: support vector machines, regularization, optimization, and beyond . MIT press, 2001 . [35] D. Sejdinovic , B. Sriperumbudur , A. Gretton , and K. Fukumizu . Equivalence of distance-based and RKHS-based statistics in hypothesis testing . Ann. Statis., 41 ( 5 ): 2263 - 2291 , 2013 . [36] J. S. Sekhon . The neyman-rubin model of causal inference and estimation via matching methods , 2008 . [37] I. Shrier . Propensity scores. Statistics in Medicine, 28 ( 8 ): 1317 - 1318 , 2009 . [38] O. Stegle , D. Janzing , K. Zhang , J. M. Mooij , and B. Scho ¨lkopf. Probabilistic latent variable models for distinguishing between cause and effect . In Advances in Neural Information Processing Systems , pages 1687 - 1695 , 2010 . [39] I. Steinwart , D. Hush , and C. Scovel . An explicit description of the reproducing kernel hilbert spaces of Gaussian RBF kernels . IEEE Trans. Inform. Theory , 52 ( 10 ): 4635 - 4643 , 2006 . [40] E. A. Stuart . Matching methods for causal inference: A review and a look forward . Statistical science: a review journal of the Institute of Mathematical Statistics , 25 ( 1 ): 1 , 2010 .