A robust approach to quantifying uncertainty in matching problems of causal inference

TL;DR

This paper introduces a new methodology to quantify and account for the uncertainty in match assignment in observational causal inference, ensuring more robust hypothesis testing regardless of the specific matching used.

Contribution

It develops a discrete optimization-based approach to create robust tests that consider all plausible match assignments, a previously overlooked source of uncertainty.

Findings

The proposed method can produce more reliable causal conclusions.

It effectively accounts for match assignment variability in both simulated and real data.

The approach is computationally feasible using standard integer linear programming techniques.

Abstract

Unquantified sources of uncertainty in observational causal analyses can break the integrity of the results. One would never want another analyst to repeat a calculation with the same dataset, using a seemingly identical procedure, only to find a different conclusion. However, as we show in this work, there is a typical source of uncertainty that is essentially never considered in observational causal studies: the choice of match assignment for matched groups, that is, which unit is matched to which other unit before a hypothesis test is conducted. The choice of match assignment is anything but innocuous, and can have a surprisingly large influence on the causal conclusions. Given that a vast number of causal inference studies test hypotheses on treatment effects after treatment cases are matched with similar control cases, we should find a way to quantify how much this extra source of uncertainty impacts results. What we would really like to be able to report is that \emph{no matter} which match assignment is made, as long as the match is sufficiently good, then the hypothesis test result still holds. In this paper, we provide methodology based on discrete optimization to create robust tests that explicitly account for this possibility. We formulate robust tests for binary and continuous data based on common test statistics as integer linear programs solvable with common methodologies. We study the finite-sample behavior of our test statistic in the discrete-data case. We apply our methods to simulated and real-world datasets and show that they can produce useful results in practical applied settings.

Figures

Tables

Table 1. Table 1 : Agreement between rejection decisions on the same data between different matching methods. In each a 𝑎 a / b 𝑏 b entry, the a 𝑎 a is the agreement on rejecting the null between methods, and b 𝑏 b measures how often covariate balance agrees between the two methods. The bottom rows show some summary statistics of the replication experiment. Results in this table are a summary of results from Morucci and Rudin ( 2018 ) . This table shows that a seemingly arbitrary choice of matching method can heavily impact the results.

Rejection agreement / Balance agreement
		CEM	Genetic	Optimal	Nearest	Subclass	Optimal
	CEM	–	0.15 / 0.82	0.40 / 0.68	0.60 / 0.71	0.41 / 0.74	0.23 / 0.62
	Genetic		–	0.27 / 0.63	0.51 / 0.68	0.68 / 0.63	0.21 / 0.83
	Optimal			–	0.42 / 0.99	0.44 / 0.94	0.40 / 0.75
	Nearest				–	0.64 / 0.95	0.31 / 0.78
	Subclass					–	0.62 / 0.83
	Optimal						–
Number of hypotheses tested: 68 times 6 methods to test each hypothesis
% of tests with at least one rejection and one non-rejection: 64.7
% of tests with at least one positive rejection and one negative rejection: 5.9
% agreement between two methods, on average: 40.5

Table 2. Table 2 : Methods used in the simulation studies.

Name	Description	Citation
robust	Robust test	This paper
true	True potential outcomes	Benchmark
naive	No matching at all	Benchmark
pscore	Propensity score matching	(Rosenbaum 1984)
l2	Matching on the $L_2$ distance of the covariates	(e.g., Abadie and Imbens 2006)
optimal	Optimal matching via Mixed Integer Programming	(Zubizarreta 2012)

Table 3. Table 3 : What pairs can be created by exchanging matches. Cells are the resulting pairs when a pair in the left margin is exchanged with a pair in the top row.

	A	B	C	D
A	A,A	A,B	A,C	B,C
B	A,B	B,B	A,D	B,D
C	A,C	A,D	C,C	C,D
D	B,C	B,D	C,D	D,D

Table 4. Table 4 : P-values for McNemar and Z-test conducted on datasets for Case Studies 1 and 2 after matching with different traditional methods

	No Matching	Pscore	L2	Optimal
GLOW	0.00	1.00	1.00	0.42
Bike Share	0.00	0.20	0.32	0.00

Equations

$$\tau =\frac{1}{N}\sum _{i=1}^N{\mathbb{E}}_{Y_i|X_i=x_i}[Y_i(1)-Y_i(0)|X_i=x_i].$$

$${\mathbb{H}}_0:\tau =0|X_1=x_1,\dots ,X_N=x_N.$$

$${\mathbb{H}}_0|\mathbf{A}=\mathbf{a}:\frac{1}{M}\sum _{i=1}^N\mathbb{E}[Y(1)-Y(0)|X=x_i]\mathbb{I}(\mathbf{a}(i)\ne \varnothing )=0.$$

$$\sum _{\mathbf{a}\in {\mathcal{A}}_{\text{good}}}\mathrm{Pr}({\Psi }_{\mathcal{D}}(\mathbf{a})\le {\psi }_{\mathcal{D}}(\mathbf{a})|{\mathbb{H}}_0,\mathbf{A}=\mathbf{a})\mathrm{Pr}(\mathbf{A}=\mathbf{a}|{\mathbb{H}}_0)$$

$$={\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}[F({\psi }_{\mathcal{D}}(\mathbf{A}))].$$

$${\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}[F({\psi }_{\mathcal{D}}(\mathbf{A}))]$$

$${\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}[F({\psi }_{\mathcal{D}}(\mathbf{A}))]\le F({\psi }^{+}).$$

$${\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}[1-F({\psi }_{\mathcal{D}}(\mathbf{A}))]$$

$$2\min \left\{{\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}\left[F({\psi }_{\mathcal{D}}(\mathbf{A}))\right],{\mathbb{E}}_{\mathbf{A}|{\mathbb{H}}_0}\left[1-F({\psi }_{\mathcal{D}}(\mathbf{A}))\right]\right\}\le 2\min \{F({\psi }^{+}),1-F({\psi }^{-})\}.$$

$$\chi =\frac{B(\mathbf{a})-C(\mathbf{a})-1}{\sqrt{B(\mathbf{a})+C(\mathbf{a})}}=\frac{{\sum }_{i=1}^{N^t}{\sum }_{j=1}^{N^c}{a}_{ij}(y_i^t-y_j^c)-1}{\sqrt{{\sum }_{i=1}^{N^t}{\sum }_{j=1}^{N^c}{a}_{ij}(y_i^t+y_j^c-2y_i^t{y}_j^c)}}.$$

$${\chi }^{+}$$

$${\text{Maximize/Minimize}}_{\mathbf{a}}\chi (\mathbf{a})=\left[\frac{B(\mathbf{a})-C(\mathbf{a})-1}{\sqrt{m}}\right]\text{ subject to:}$$

$$\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}{a}_{ij}{y}_i^t(1-y_j^c)$$

$$\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}{a}_{ij}{y}_j^c(1-y_i^t)$$

$$B(\mathbf{a})+C(\mathbf{a})$$

$$\sum _{i=1}^{N^t}{a}_{ij}$$

$$\sum _{j=1}^{N^c}{a}_{ij}$$

$$a_{ij}$$

$$z_{\mathbf{y}}(\mathbf{a})$$

$$z^{+}:=\underset{\mathbf{a}\in {\mathcal{A}}_{\text{good}}}{\max }\frac{{\bar{d}}_{\mathbf{a}}\sqrt{M}}{{\hat{\sigma }}_{\mathbf{a}}},z^{-}:=\underset{\mathbf{a}\in {\mathcal{A}}_{\text{good}}}{\min }\frac{{\bar{d}}_{\mathbf{a}}\sqrt{M}}{{\hat{\sigma }}_{\mathbf{a}}}.$$

$${\text{Maximize/Minimize}}_{\mathbf{a}}\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c)a_{ij},$$

$$\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c{)}^2{a}_{ij}\le b_l.$$

$${\text{Maximize/Minimize}}_{\mathbf{a}}\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c){\text{a}}_{ij}\text{(Treatment effect), subject to}$$

$$\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c{)}^2{a}_{ij}$$

$$\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}{a}_{ij}$$

$$\sum _{i=1}^{N^t}{a}_{ij}$$

$$\sum _{j=1}^{N^c}{a}_{ij}$$

$$a_{ij}$$

$$\underset{l}{\max }\frac{{\bar{d}}_{{\mathbf{a}}_l}\sqrt{M}}{\sqrt{\frac{1}{M}{b}_l-({\bar{d}}_{{\mathbf{a}}_l}{)}^2}}\le \underset{\mathbf{a}}{\max }z(\mathbf{a})\le \underset{l}{\max }\frac{{\bar{d}}_{{\mathbf{a}}_l}\sqrt{M}}{\sqrt{\frac{1}{M}{b}_{l-1}-({\bar{d}}_{{\mathbf{a}}_l}{)}^2}}.$$

$$\frac{{\bar{d}}_{{\mathbf{a}}_l}\sqrt{M}}{\sqrt{\frac{1}{M}{b}_{l-1}-({\bar{d}}_{{\mathbf{a}}_l}{)}^2}}<\underset{\ell }{\max }\frac{{\bar{d}}_{{\mathbf{a}}_{\ell }}\sqrt{M}}{\sqrt{\frac{1}{M}{b}_{\ell }-({\bar{d}}_{{\mathbf{a}}_{\ell }}{)}^2}}.$$

$$\underset{l}{\max }F(f_1(x_l),b_l)\le \underset{x}{\max }F(f_1(x),f_2(x))\le F(f_1(x_{l^{\ast }}),b_{l^{\ast }-1})\le \underset{l}{\max }F(f_1(x_l),b_{l-1}).$$

$$\underset{l}{\min }F(f_1(x_l),b_l)\ge \underset{x}{\min }F(f_1(x),f_2(x))\ge F(f_1(x_{l^{\ast }}),b_{l^{\ast }-1})\ge \underset{l}{\min }F(f_1(x_l),b_{l-1}).$$

$$f_1(\mathbf{a})=\frac{1}{M}\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c)a_{ij},f_2(\mathbf{a})=\sum _{i=1}^{N^t}\sum _{j=1}^{N^c}(y_i^t-y_j^c{)}^2{a}_{ij},\text{ and}$$

$$F(f_1(\mathbf{a}),f_2(\mathbf{a}))=\frac{f_1(\mathbf{a})\sqrt{M}}{\sqrt{\frac{1}{M}{f}_2(\mathbf{a})-(f_1(\mathbf{a}){)}^2}}.$$

$${\mathbb{H}}_0^{sharp}:=Y_i(1)=Y_i(0).$$

Full text

1 Introduction

We have a reproducibility crisis in science. Part of the reason for the crisis is sources of uncertainty in the analysis pipeline that are not accounted for. Observational causal studies have serious problems with reproducibility, despite the fact that these studies underlie important policy decisions. In our paper, we focus on a source of uncertainty that can have a large impact on the result, but is almost never mentioned: match assignments of treatment units to control units.

Classically, assignments of treatment and control units to matches are constructed using a single method aimed at constructing pairs that achieve a good level of some measure of quality predefined by the analyst (e.g., Optimal Matching, Rosenbaum 1989, or Genetic Matching, Diamond and Sekhon 2013). Choosing a single fixed match ignores a major source of uncertainty, which is the design itself, or in other words, the uncertainty related to the choice of experimenter. What if there were two possible equally good match assignments, one where the treatment effect estimate is very strong and one where it is nonexistent? When we report a result on a particular matching assignment, we thus ignore the possibility of the opposite result occurring on an equally good assignment. It is entirely possible that two separate researchers studying the same effect on the same data, using two different equally good sets of matched groups, would get results that disagree. When researchers follow the classic pipeline of match-then-estimate, their hypothesis tests are conditioned on the match assignment; in other words, only uncertainty after matching is considered, and not the uncertainty stemming from the match assignment itself.

Our goal is to create robust matched pairs hypothesis tests for observational data. These tests implicitly consider all possible reasonably good assignments and consider the range of possible outcomes for tests on these data. This is a more computationally demanding approach to hypothesis testing than the standard approach where one considers just a single assignment, but the result is then more robust to the choice of experimenter who chooses the match assignment. It is computationally infeasible (and perhaps not very enlightening) to explicitly compute all possible assignments, but it is possible to look at the range of outcomes associated with them. In particular, our algorithms compute the maximum and minimum of quantities like the test statistic values and their associated $p$-values, in accordance with the principles of robust optimization.

After motivation and formalization for our framework in Sections 1.1 and 2, we offer a robust formulation for McNemar’s statistic for binary outcomes in Section 3. In Section 4 we formulate an integer linear program (ILP) to compute the robust version of the canonical $z$-test in a general case. We then study the finite-sample distribution of our robust version of McNemar’s test that includes the uncertainty stemming from the match assignment. Finally, we present evidence of the performance of our methods in Section 6 using simulated data where the ground truth is known, and two real-data case studies are discussed in Section 7. We know of no previous approaches to matching in observational causal inference that use robust optimization to handle uncertainty in the match assignment.

1.1 Matching Analyses and Their Sensitivity to Arbitrary Choices

Existing studies that use matching for hypothesis testing all roughly follow the following template: (Step 1) Choose a test statistic, (Step 2) Define criteria that acceptable matches should satisfy (e.g., covariate balance), (Step 3) Use an algorithm to find matched groups (and corresponding subsamples) that satisfy the criteria, (Step 4) Implicitly choose the subsample of interest to be defined by the matched data, (Step 5) Compute the test statistic and its p-value on the matched data (see, e.g., Rosenbaum 2010). This procedure does not explicitly incorporate uncertainty arising from Steps 3 and 4: values of the test statistic computed under different but equally good matched subsamples could differ. These matches are indistinguishable in terms of quality. Which of the matched sets should be taken as representative of the whole sample? There is no clear answer to this question.

Even a trivial source of randomness, such as the order in which the data appear in the dataset, could produce two equivalently good match assignments under a variety of popular matching methods, such as nearest-neighbors. If we apply the same test statistic to each of the two matched samples, we would be implicitly testing two different hypotheses: the first one is defined on the subset of the data selected by the first match assignment, and the second hypothesis on the second. Under the assumption that both matches are equally as good, common statistical tests will yield valid results for each of these two hypotheses, yet the hypothesis test conducted on the first matched sample could reject while the test conducted on the second could not. Since both assignments achieve similar quality, there is no clear way to choose one result over the other, and arbitrarily reporting one of the two would hide the fact that the other exists. Such uncertainty would likely be of interest to policy makers and researchers who would want their results to be reproducible under all arbitrary choices of analyst.

There is no easy fix. Assigning multiple control units to each treatment unit, or using deterministic matching methods that yield a unique result for each dataset, ignore the possibility of a slightly different – but still good – match assignment giving opposite results. Asymptotics (e.g., Abadie and Imbens 2011) provide no remedy, since asymptotic results are usually intended to apply to all good matching algorithms asymptotically; the problem we encounter in practice occurs because we always work with finite samples. The problem extends not only to treatment effects but also to variance estimates and beyond.

1.2 Empirical Evidence of the Problem

The problem we identified is serious in practice. Morucci and Rudin (2018) replicate several social science studies – all published in top journals after 2010 – that use matching. They perform the same hypothesis tests after applying several popular matching algorithms. Table 1 reports agreement between any two of six different matching methods on the same hypothesis test, i.e., the percentage of the tests replicated that have the same result (reject/fail to reject) under each pair of methods. If our hypothesis were false, and matching at the same level of quality always produced the same results, then we would see high agreement in the table. This is clearly not the case. Additionally, Table 1 also reports a measure of similarity of the balance between the datasets produced by each pair of matching methods. We measure balance as the proportion of covariates whose difference in means between the matched treated and control group is statistically insignificant at the 5% level. The numbers in $b$ of each $a/b$ entry of Table 1 represent the correlation between indicator vectors that keep track of which covariates were balance by each method according to the criterion just presented. Clearly, while imperfectly, methods tend to “agree” on balance more than they do on rejection: agreement on balance is always greater than agreement on rejection, sometimes by a large margin. This implies that differences in balance produced are not, at least in full, to blame for the problem of different matching methods leading to different conclusions: even when methods produce similar balance, they still disagree.

On average, two methods alone agree on rejecting vs. not rejecting the same null hypothesis only 40.5% of times but achieve similar balance 70.5% of times. These numbers drop dramatically if any group of three methods is considered, and are virtually zero if agreement among all six methods at once is considered. As stated before, we believe that this happens because different matching methods choose among matches of similar quality arbitrarily and in different ways. According to these results, seemingly arbitrary choices of match assignment affect conclusions. This is true in practical scenarios that involve real-world data, and that addressing it is an important step towards fully robust data analysis.

2 Proposed Approach

Throughout this paper we adopt the potential outcomes framework (see Holland 1986, Rubin 1974). For each unit $i=1,\dots,N$, we have potential outcomes $Y_{i}(t)\in\mathbb{R}$, where $t\in\{0,1\}$. As is standard in causal inference settings, there are $N^{t}$ units that receive the treatment and $N^{c}$ units that receive the control condition, with $N^{c}+N^{t}=N$. We denote the condition administered to each unit with $T_{i}\in\{0,1\}$. We never observe realizations of $Y_{i}(1)$ and $Y_{i}(0)$ for each unit at the same time, but only of $Y_{i}=Y_{i}(1)T_{i}+Y_{i}(0)(1-T_{i})$. We also observe a $P$-dimensional vector of covariates $X_{i}$ taking value in some set $\mathcal{C}$ for every unit. We will denote the complete set of random variables representing the data with $\mathcal{D}=\{X_{i},T_{i},Y_{i}\}_{i=1}^{N}$, with each element being an independent draw from some probability distribution. Analogously, we will denote the set of observed values for the data with the lowercase equivalents of the notation above: $\{x_{i},t_{i},y_{i}\}_{i=1}^{N}$. We assume that potential outcomes are related to the covariates as follows: $Y_{i}(t)=\mu_{t}(x_{i})+\nu_{i}$, with $\mathbb{E}[\nu_{i}|X_{i}=x_{i}]=0$. To conduct our hypothesis test we will also notate our observations as: $(x_{1}^{t},y_{1}^{t}),...,(x_{N^{t}}^{t},y_{N^{t}}^{t})$, and $(x_{1}^{c},y_{1}^{c}),...,(x_{N^{c}}^{c},y_{N^{c}}^{c})$, where $x_{i}^{t},y_{i}^{t}$ and $x_{i}^{c},y_{i}^{c}$ represents observed values of $X$ and $Y$ for the $i^{th}$ treated and control units respectively. We make the classical assumptions of Conditional Ignorability of treatment assignment and Stable Unit Treatment Value (SUTVA):

Assumption 1**.**

(Strong Conditional Ignorability) For any unit in the sample $i=1,\dots,N$ treatment allocation is independent of potential outcomes conditional on the observed covariates, that is: T_{i}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}(Y_{i}(1),Y_{i}(0))|X_{i}.

Assumption 2**.**

(SUTVA) A unit’s potential outcomes depend only on that unit’s treatment level, that is: for all units, $i$: $Y_{i}(t_{1},\dots,t_{N})=Y_{i}(t_{i})$, where $Y_{i}(t_{1},\dots,t_{N})$ is the potential outcome for unit $i$ under all units’ treatment assignments.

Under these two assumption, the Sample Average Treatment Effect (SATE) is our quantity of interest and is defined as:

[TABLE]

Under SUTVA and Conditional Ignorability, $\tau$ is identified – it can be consistently estimated with the observed data – and hypothesis tests based on it can be conducted in a variety of ways. In general, we will consider testing a null hypothesis of 0 average treatment effect in the observed sample, which is canonically defined as follows:

[TABLE]

Note that the defining feature of the SATE is that the contextual covariates, $\mathbf{X}$ are only considered at the observed values, as implied by the definition of $\mathbb{H}_{0}$. As it is done in much of the causal inference literature, we would like to take advantage of matching in order to test $\mathbb{H}_{0}$. A matching operator determines which control is assigned to which treatment. In this paper we focus on one-to-one matching (though our results could be generalized to many-to-one or many-to-many matching), therefore we define the matching operator as follows.

Definition 1**.**

(Matching Operator) A matching operator $\mathbf{a}:\{1,...,N^{t}\}\rightarrow\{1,...,N^{c},\emptyset\}$ obeys the following: if $i\neq k$ and $\mathbf{a}(i)\neq\emptyset$ then $\mathbf{a}(i)\neq\mathbf{a}(k)$. That is, no two treatment units $i$ and $k$ are assigned to the same control unit.

We define the size of the matching, i.e.: the number of matched pairs, to be $M=\sum_{i=1}^{N_{t}}\mathbb{I}(\mathbf{a}(i)\neq\emptyset)$, with $\mathbb{I}(E)$ representing the indicator function for the event $E$ throughout the paper. The set of all matching operators is $\mathcal{A}$: the set of all possible assignments $\{\mathbf{a}\}$. Throughout the rest of the paper we also use $\mathbf{a}$ to represent a $N^{t}\times N^{c}$ matrix of matches such that $a_{ij}=1$ if $i$ and $j$ are matched together.

How should we make our matches? Choosing a matched set that leads to the most reliable estimates of $\tau$ is a canonical problem in observational causal inference, and indeed many good methods to find such an assignment have been proposed in the existing literature. One of the contributions of this paper is to not restrict consideration to a single good assignment, but, instead, to consider a set of potentially many equally good match assignments. We define such a set as $\mathcal{A}_{\textrm{good}}\subset\mathcal{A}$: the set of match assignments that satisfy some user-defined criterion of quality. Many criteria for defining a good match assignment already exist, e.g.: moment balance in the covariate distributions of treated and control samples after matching, low aggregate distance between units in covariate space, or similarity in propensity scores among matched units. Ultimately, $\mathcal{A}_{\textrm{good}}$ should be viewed as the set of match assignments that produce the most reliable estimates of $\tau$ with the given data. To represent the fact that one assignment is chosen among the many in $\mathcal{A}_{\textrm{good}}$, we will use the random variable $\mathbf{A}$, having domain equal to $\mathcal{A}_{\textrm{good}}$. The probability distribution of $\mathbf{A}$ is determined purely by human analysts’ choice of matching algorithm.

Consider, for now, testing $\mathbb{H}_{0}$ against a left-tailed alternative: $\mathbb{H}_{1}:\tau<0$ with a generic test statistic $\Psi_{\mathcal{D}}(\mathbf{a})$, dependent on both the data, $\mathcal{D}$, and the chosen match assignment, $\mathbf{a}$, and denote the observed value of $\Psi_{\mathcal{D}}(\mathbf{a})$ with $\psi_{\mathcal{D}}(\mathbf{a})$. Virtually all hypothesis tests done on data matched by any of the existing methods require the following assumption, that we also maintain:

Assumption 3**.**

(Test statistic distribution) Under Assumption 1 and Assumption 2, we assume that for any assignment $\mathbf{a}\in\mathcal{A}_{\textrm{good}}$: $\Pr(\Psi_{\mathcal{D}}(\mathbf{a})\leq\psi_{\mathcal{D}}(\mathbf{a})|\mathbb{H}_{0},\mathbf{A}=\mathbf{a})\approx F(\psi_{\mathcal{D}}(\mathbf{a}))$, for some known cumulative density function $F$ that does not depend on $\mathbf{a}$.

This assumption means that, if the matches are good enough, then the CDF of $\Psi_{\mathcal{D}}(\mathbf{a})$ under the null can be at least approximated well by some known function, $F$. This assumption is implicitly made in all studies that perform a hypothesis test on matched data. It is often justified by appealing to asymptotic arguments (Abadie and Imbens 2011), or by making parametric assumptions about the data (Rubin 2007). For example, if a researcher performs a z-test on matched data, they implicitly assume that $F$ is the CDF of the standard normal distribution.

The key problem we wish to highlight in this paper can be stated as follows: existing matching procedures wishing to test $\mathbb{H}_{0}$ under Assumption 3 use $\Pr(\Psi_{\mathcal{D}}(\mathbf{a})\leq\psi_{\mathcal{D}}(\mathbf{a})|\mathbb{H}_{0},\mathbf{A}=\mathbf{a})$ to compute the p-value for $\Psi_{\mathcal{D}}(\mathbf{a})$ under $\mathbb{H}_{0}$, however this p-value is both conditioned on $\mathbb{H}_{0}$ and the event $\mathbf{A}=\mathbf{a}$, i.e., conditional on assignment $\mathbf{a}$ being the one chosen among the many in $\mathcal{A}_{\textrm{good}}$. The hypothesis being tested here is not $\mathbb{H}_{0}$, but instead: $\mathbb{H}_{0}|\mathbf{A}=\mathbf{a}$:

[TABLE]

This is essentially a version of $\mathbb{H}_{0}$ restricted to consider only those units that do receive a match under assignment $\mathbf{a}$, and, most importantly, there is no guarantee that we will fail to reject $\mathbb{H}_{0}|\mathbf{A}=\mathbf{a}$ if $\mathbb{H}_{0}$ is true. Thus, while testing $\mathbb{H}_{0}|\mathbf{A}=\mathbf{a}$ rather than $\mathbb{H}_{0}$ seems not to be what the user ultimately desires, it is the calculation they typically perform. By overlooking the difference between these two quantities, they would fail to ask questions such as: how should we proceed if we find two assignments $\mathbf{a}_{1}\in\mathcal{A}_{\textrm{good}}$ and $\mathbf{a}_{2}\in\mathcal{A}_{\textrm{good}}$, that successfully reject $\mathbb{H}_{0}|\mathbf{A}=\mathbf{a}_{1}$ but fail to reject $\mathbb{H}_{0}|\mathbf{A}=\mathbf{a}_{2}$ at the same significance level? We now articulate our proposed solution to this problem.

Since we have modeled which assignment is selected among those in $\mathcal{A}_{\textrm{good}}$ as a random variable, it would be natural, but ultimately impossible, to propose computing p-values in expectation over the entirety of $\mathcal{A}_{\textrm{good}}$ to handle the issue just introduced. This would be possible under Assumption 3 by computing:

[TABLE]

This quantity cannot be computed due to the presence of $\Pr(\mathbf{A}=\mathbf{a}|\mathbb{H}_{0})$ in its definition: the probability of observing any match assignment among the ones in $\mathcal{A}_{\textrm{good}}$ cannot be modeled statistically. This is because this probability depends on factors such as which matching method is chosen by the analyst, which hyperparameter values the method is using, and how the method itself may determine how to choose one among many potentially equivalent matched sets. All these are variables that cannot be analytically studied in a general context, and will depend on the specific analyst, data, and method.

Since directly modeling $\Pr(\mathbf{A}=\mathbf{a}|\mathbb{H}_{0})$ is nonsensical, we propose to instead bound this quantity using the assumption that all analysts choose reasonably good match assignments. This assumption leads to the following simple relationship:

[TABLE]

In words, we bound the expected p-value over all good assignments with the largest p-value found among all good assignments. Bounding approaches to handle situations without clear probability distributions are common in robust optimization and robust statistics.

To further reduce the problem, we note that the CDF of the test statistic, $F(\psi_{\mathcal{D}}(\mathbf{a}))$, is monotonically increasing in $\psi_{\mathcal{D}}(\mathbf{a})$ (e.g., the larger the z-statistic, the larger the CDF of the normal), which implies that $\max_{\mathbf{a}\in\mathcal{A}_{\textrm{good}}}F(\psi_{\mathcal{D}}(\mathbf{a}))=F(\max_{\mathbf{a}\in\mathcal{A}_{\textrm{good}}}\psi_{\mathcal{D}}(\mathbf{a}))$. Letting $\psi^{+}=\max_{\mathbf{a}\in\mathcal{A}_{\textrm{good}}}\psi_{\mathcal{D}}(\mathbf{a})$ be the largest observed test statistic obtainable with a good match assignment, we have our final proposed bound for the robust p-value:

[TABLE]

Clearly, if $\mathbb{H}_{0}$ is rejected at a chosen level of significance by $F(\psi^{+})$, then it will also be rejected under $\mathbb{E}_{\mathbf{A}|\mathbb{H}_{0}}[F(\psi_{\mathcal{D}}({\mathbf{A}}))]$. It is in this sense that our proposed tests are robust, as they target minimization of the probability of Type-I error: incorrect rejection of $\mathbb{H}_{0}$. That is, our tests are less likely to reject the null hypothesis when it is true.

While the bound in Eq. (3) is a robust p-value against a left-sided alternative, the same set of arguments shown above can be applied to derive the following upper bound on the p-value for a right-sided alternative: $\mathbb{H}_{1}:\,\tau>0$:

[TABLE]

Finally, it follows from the bounds above that a robust p-value for testing $\mathbb{H}_{0}$ against a two-sided alternative: $\mathbb{H}_{1}:\,\tau\neq 0$ can be found by applying the canonical definition of p-value to the bounds introduced above:

[TABLE]

The problem we face is now one of finding, for a given dataset, the values of $\psi^{+}$ and $\psi^{-}$. Fortunately, this issue can be solved by adapting modern Mixed-Integer-Programming tools to finding the maximal and minimal matches that satisfy the set of constraints that define $\mathcal{A}_{\textrm{good}}$. In the rest of the paper, we propose MIP formulations and algorithms for computing these values for two of the most widely used test statistics.

2.1 Motivation for Using a Robust Approach

This robust approach is justified as follows:

Why modeling $\Pr(\mathbf{A}=\mathbf{a}|\mathbb{H}_{0})$ is not possible.

First, we do not want to consider how likely a certain assignment $\mathbf{a}$ is to appear. This would involve modeling human behavior of the experimenter, or arbitrary and complex choices of different matching algorithms. We do not want to place a distribution over the choice of algorithms or matches that an experimenter would choose. As Morgan and Winship (2007) note, there is no clear guidance on the choice of matching procedure. We do not presuppose a distribution over these procedures. In reality most researchers tend to use popular and widely-cited matching software packages, and they make this choice independently of the specific problem being studied.

Why the bounding approach is not too extreme.

It is inaccurate to assume that our tests might lead to results that are “too” extreme. This is because there is no reason why maxima or minima over $\mathcal{A}_{\textrm{good}}$ should be considered outlying or extreme, even if most of the matches in $\mathcal{A}_{\textrm{good}}$ are different than those at the extrema of the set. Note also that considering a result output by one of the extrema as an “outlier” in the sense that it is believed that most other results with good matches should be concentrated elsewhere in the space of $\psi_{\mathcal{D}}(\mathbf{a})$ is circular logic: if a set of assignments is included in $\mathcal{A}_{\textrm{good}}$, then it must be thought to be good enough according to the analyst’s criteria. Excluding this assignment because it produces values that are too extreme is tantamount to excluding observations from the dataset because they are not producing the desired result.

It is unlikely that there is only one good assignment.

It would be tempting to say that matches chosen are often the unique maximizers of some metric of quality, and that, because of this, there never is uncertainty over choice of matches in practice. This is not true in practice, as researchers often start with a pre-defined level of quality that matches must meet to be considered acceptable, and only report results from the assignment that maximizes that quality criterion if the maximum achieved is above that pre-defined quality threshold. Clearly, all assignments above that threshold should be considered, and that is what we propose in this work.

Asymptotic theory does not solve the problem.

Our robust method is also necessary because existing asymptotic arguments used to derive valid p-values and confidence intervals for matched data (see, e.g., Abadie and Imbens 2006, 2011). As argued previously, the defining feature of a “good” match assignment is that it leads to at least well approximated p-values for the desired test statistic. In the case of asymptotic approximation, it is assumed that $\Pr(\Psi_{\mathcal{D}}(\mathbf{a})\leq\psi_{\mathcal{D}}(\mathbf{a})|\mathbb{H}_{0},\mathbf{A}=\mathbf{a})\rightarrow F(\psi_{\mathcal{D}}(\mathbf{a}))$ as $N\rightarrow\infty$. This is assumed true for all assignments $\mathbf{a}\in\mathcal{A}_{\textrm{good}}$: while it might be the case that the p-values produced under all good assignments will converge to the same as the sample size grows, this is demonstrably not necessarily true in finite samples, hence the need for a robust procedure.

3 Robust McNemar’s Test

We consider the problem of hypothesis testing for binary outcomes, that is $y(t)\in\{0,1\}$. After matching, treated and control units will be paired together, and since the outcome is binary there can exist four types of pairs: $A$ is the number of pairs such that $y_{t}=0,\,y_{c}=0$, $B$ is the number of pairs such that $y_{t}=1,\,y_{c}=0$, $C$ is the number of pairs such that $y_{t}=0,\,y_{c}=1$, $D$ is the number of pairs such that $y_{t}=1,\,y_{c}=1$. More formally, let $B(\mathbf{a})=\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}a_{ij}y_{i}^{t}(1-y_{j}^{c})$ be the count of matched pairs in which the treated unit has outcome 1 and the control unit has outcome 0 under assignment $\mathbf{a}$, and let $C(\mathbf{a})=\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}a_{ij}y_{j}^{c}(1-y_{i}^{t})$ be the number of matched pairs in which the treated unit has outcome 0 and the control unit has outcome 1. We use the following test statistic for all null hypotheses:

[TABLE]

This is the classical formulation of McNemar’s statistic with a continuity correction at the numerator (see Tamhane and Dunlop 2000). Our proposed robust test statistic can be defined as the pair:

[TABLE]

In what follows, we outline a strategy to compute $(\chi^{+},\chi^{-})$ subject to general constraints that define $\mathcal{A}_{\textrm{good}}$. Theorem 1 in Section 3.2 shows that, in a special case where the constraints on $\mathcal{A}_{\textrm{good}}$ explicitly define strata for the data, the optimal statistic can be computed in linear time. In Section 5, we provide theoretical statements of the finite-sample behavior of the optimized statistics.

3.1 Optimizing McNemar’s statistic with general constraints

In this section we give an ILP formulation that optimizes $\chi$ with a pre-defined number of matches.

One can show that the number of pairs where the same outcome is realized for treatment and control (the number of $A+D$ pairs) is irrelevant, so we allow it to be chosen arbitrarily, with no constraints or variables defining it. The total number of pairs is thus also not relevant for this test. Therefore, we are permitted to choose only the total number of untied responses ($B$ and $C$ pairs), denoted $m$, which must always be greater than 0 for the statistic to exist. The problem can be solved either for a specific $m$ or globally for all $m>0$ by looping over all possible values of $m$ until the problem becomes infeasible. In most practical scenarios, the largest $m$ for which a solution is feasible is chosen.

Formulation 1: General ILP Formulation for McNemar’s Test

[TABLE]

[TABLE]

Equations (7) and (8) are used to define variables $B$ and $C$. To control the total number of untied responses, we incorporate Equation (9). Equations (10) and (11) confirm that only one treated/control unit will be assigned in a single pair. As it is common practice in matching, it is possible to test a null hypothesis for the SATT by adding the constraint $\sum_{i=1}^{N^{t}}\sum_{C}{a}_{ij}=N^{t}$, only in the case where $N^{t}\leq N^{c}$, so that all treatment units are matched.

3.2 Fast Optimization of McNemar’s Test Under Special Constraints

The general optimization problem stated in Formulation 1 can also be solved quickly when constraints have the following property: there exists a stratification of the data such that an assignment is a feasible solution if and only if it matches units exclusively within the same stratum. This is stated formally below:

Definition 2**.**

(Exclusively Binning Constraints) The constraints on $\mathcal{A}_{\textrm{good}}$ are exclusively binning if there exists a partition $\mathcal{S}$ of $\{1,\dots,N\}$ such that for all $\mathbf{a}\in\mathcal{A}$ we have $\mathbf{a}\in\mathcal{A}_{\textrm{good}}$ if and only if: $\forall\;i\in\{1\dots,N\}$, $i\in S$, and $\mathbf{a}(i)\in S$ for some subset $S\in\mathcal{S}$.

This type of grouping is commonly referred to as blocking (Imai et al. 2008), and there already exist matching methods that openly adopt it as a strategy to construct good quality matches (Iacus et al. 2012). Blocking can also occur naturally in data; there need only exist natural subcommunities. In practice, several types of constraints on the quality of matches, particularly balance constraints, as defined in Equation (31), can be implemented by making coarsening choices on the covariates (Iacus et al. 2011), and as such, these coarsening choices are exclusively binning constraints. Under exclusively binning constraints, solving the optimization problem to find $(\chi^{+},\chi^{-})$ becomes much simpler, as shown in the following theorem:

Theorem 1**.**

Let the constraints on $\mathcal{A}_{\textrm{good}}$ be Exclusively Binning Constraints. Then either the max or the min optimization problem in Formulation 1 can be solved in linear time, in the number of units $N$.

Proof of the theorem can be found in the appendix, where we explicitly state these linear-time algorithms, together with a proof of their correctness. These algorithms are not difficult, though stating them and providing the proof requires several pages. The algorithms work by first determining the sign of the optimal test statistic, and then by matching units in each stratum separately to optimize the test statistic once the sign is known. Maximizing (minimizing) the test statistic in each stratum is then accomplished by matching as many treated (control) units as possible with outcome 1 to control (treated) units with outcome 0. This procedure is intuitively linear in $N$, the total number of units. Figure 1 gives a simple example of how our linear-time algorithms optimize McNemar’s test in each stratum.

4 Robust Z-Test

In this section we consider the canonical $z$-test for estimating whether the difference in mean of the treatment and control populations is sufficiently greater than 0, when outcomes are real-valued, that is: $y(t)\in\mathbb{R}$. Again we will compute extreme values of $z$ that could occur from the set of feasible match assignments.

Again, after computing the extreme values of the z-statistics, p-values for them can be obtained with the canonical Normal asymptotic approximation. At that point, we can determine whether the hypothesis test result is robust to the choice of match assignment.

In what follows, $M$ is the total number of pairs, and $\hat{\sigma}$ is the sample standard deviation of the differences, $y^{t}_{i}-y_{\mathbf{a}(i)}^{c}$. The $z$-score is:

[TABLE]

Our robust statistic in this case is defined as the pair:

[TABLE]

Below, we provide general ILP formulations for computing the robust statistic under (balance and other types of) constraints on $\mathcal{A}_{\textrm{good}}$ by devising a linearized formulation of the $z$-statistic optimization problem that allows it to be solved with any ILP solver. ILP formulations that are slightly different from each other (that we will discuss) can handle testing of $\mathbb{H}_{0}$ for the ATE and ATT.

4.1 Computing $(z^{+},z^{-})$ Under General Constraints

The $z$ statistic is clearly not linear in the decision variables (the match assignments). If one were to optimize it directly, a solution could be approximated using a MINLP (mixed-integer nonlinear programming solver) but guarantees on the optimality of the solution might take an incredibly long time. In what follows, we show how this problem can be simplified to be solved by an algorithm that solves several linear integer programming problems instead. This algorithm benefits from the computational speed of ILP solvers, compared to MINLP solvers, and has a guarantee on the optimality of the solution.

To create the ILP formulation, we note that the objective is increasing in the average of the differences (this term appears both in the numerator and denominator), and it is decreasing in the sum of the squared differences (this term is the first term of $\hat{\sigma}$). We then replace the nonlinear objective in (15), expanded in (14), as follows:

[TABLE]

which is now linear in the decision variables. The quantity in (16) is the estimated treatment effect. Simultaneously, we will limit the sum of squared differences term in the denominator by $b_{l}$, which is now a parameter rather than a decision variable. Thus, we will optimize treatment effect subject to a bound on the variance. We accomplish this by introducing a new constraint:

[TABLE]

Putting this together, the new formulation is an ILP. We simply need to solve it for many values of $b_{l}$.

Formulation 3: ILP formulation for $z$-test

[TABLE]

[TABLE]

This formulation optimizes treatment effect, subject to the variance of the treatment effect being small. This formulation can be used by itself to find the range of reasonable treatment effects, given a fixed bound $b_{l}$ on the variance. The problem of testing $\mathbb{H}_{0}$ for the ATT and under full matching can be formulated by setting $M=N^{t}$ (where $N^{t}$ is the number of treatment points) in Formulation 3.

Modification 1: Formulation 3 can also be modified to solve the problem of choosing both treatment and control populations simultaneously. This is also the setting of Rosenbaum (2012). In that case, the mean is taken over the region of overlap between the control and treatment populations, removing extreme regions. This setting can be handled by looping over increasing values of $M$ until the program becomes infeasible. We would choose the solutions corresponding to the largest values of $M$.

Modification 2: For testing a sharp null hypothesis of 0 treatment effect, a simplified formulation is possible, because the sample variance is fixed and known in this case. We would use a special case of Formulation 3, where the variance constraint in (19) is replaced with an equality constraint, requiring the variance of the solution to be equal to the known sample variance. (The formula for the variance computed under the sharp null is available in textbooks, e.g., Imbens and Rubin 2015.)

Algorithm for Optimizing z-score: Let us get back to optimizing the z-score. Our algorithm will solve this formulation for many different values of $b_{l}$ to find the optimal z-scores and p-values. Let us denote the solution of the maximization problem as $\mathbf{a}_{l}$. Here, $\mathbf{a}_{l}$ is an optimal match assignment for a specific value of $b_{l}$. The indices of the match assignment $\mathbf{a}_{l}$ are, as usual, $ij$, which are pairs of treatment and control units. Using $\mathbf{a}_{l}$, we will then be able to bound the value of $z$. Shortly we will use Theorem 2 to prove bounds on the z-score as follows:

[TABLE]

That is, when $b_{l-1}$ is close to $b_{l}$, we have little uncertainty about $\max_{\mathbf{a}}z(\mathbf{a})$.

Using this bound, we can now formulate an algorithm (Algorithm 6 in Appendix Appendix D:) to choose progressively finer meshes for $b_{l}$ to maintain the guarantee on the quality of the solution for $\max_{\mathbf{a}}z(\mathbf{a})$, repeatedly solving the ILP Formulation. An analogous algorithm (with some signs flipped) will allow us to compute $\min_{\mathbf{a}}z(\mathbf{a})$. Algorithm 6 in Appendix Appendix D: works as follows:

• •

It first solves the ILP Formulation by relaxing (removing) the first constraint (upper bound on sample standard deviation).

• •

It then uses the resulting matches to compute the first upper bound on the standard deviation, $b^{(0)}_{L}$ (Line 2 in Algorithm 6).

• •

The algorithm then creates a coarse mesh $b_{1}^{(iter)},\dots,b_{L}^{(iter)}$ where $b_{1}^{(iter)}<b_{l}^{(iter)}<b_{L}^{(iter)}$ (Line 3) new refined mesh will be created at each iteration, and we denote iterations by $iter$. We want the interval $[b_{1}^{(iter)},b_{L}^{(iter)}]$ to be wide enough to contain the true value of $f_{2}(\mathbf{a}^{*}):={\displaystyle\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}(y_{i}^{t}-y_{j}^{c})^{2}a^{*}_{ij}},$ where $\mathbf{a}^{*}\in\operatorname*{arg\,max}z(\mathbf{a})$, which we do not know and are trying to obtain. Determining which procedure to use to create this mesh is left up to the user, the $b_{l}^{(iter)}$ could be chosen evenly spaced, though they do not need to be. Note that the choice of $(b_{1}^{(iter)},\dots,b_{L}^{(iter)})$ at each iteration does not affect the optimality of the solution, only the speed at which it is obtained.

• •

The algorithm then determines the sign of the maximal $z$-statistic, if negative (Line 17), it will try to maximize the denominator (variance) of the $z$ statistic, and if positive it will try to minimize it (Line 5).

• •

The algorithm then computes the solution to the ILP Formulation as well as upper and lower bounds for the solution using (25) (Lines 7-8 and 21-22).

• •

We then do some re-indexing. We take the union of all grid points $\bigcup_{l,(iter)}b_{l}^{(iter)}$ created over all iterations, order them and create a ordered single vector $\mathbf{b}=[b_{1},...,b_{L}]$. This is a single 1–D grid. On each point of this grid, we have match assignment $\mathbf{a}_{l}$ that maximizes $z$, subject to a constraint $b_{l}$ on the variance.

• •

For each $l$, the algorithm determines whether the interval $[b_{l-1},b_{l}]$ can be excluded because it provably does not contain a $f_{2}(\mathbf{a})$ value corresponding to the maximum value of $z(\mathbf{a})$ (Lines 9-12 or 22-25). In particular, we know from the bounds in (25) and from Theorem 2 that if the upper bound on the objective for a particular $b_{l^{{}^{\prime}}}$ is lower than all lower bounds for the optimal solution $\mathbf{a}^{*}$ then $l^{{}^{\prime}}$ cannot equal $l^{*}$ and the interval $[b_{l^{{}^{\prime}}-1},b_{l^{{}^{\prime}}}]$ can be excluded from further exploration. Specifically, we check for each $l$ whether

[TABLE]

If this holds for some $l$, it means $l$ cannot equal $l^{*}$ and the interval $[b_{l-1},b_{l}]$ can be excluded from further exploration.

• •

The intervals that remain included after this procedure are then refined again at each iteration, thus creating finer and finer meshes, and the process repeats on these finer meshes until most intervals are excluded and the desired tolerance ($\epsilon$) is achieved.

• •

The output is a binary matrix of match assignments, $\mathbf{a}$, and the optimal value of $z(\mathbf{a})$.

Correctness of the algorithm follows directly from optimality of the solutions at the bounds mentioned earlier and through use of the following theorem:

Theorem 2**.**

*(Optimal Solution of ILP with Upper Bound on Variance) Let the functions $f_{1}$ and $f_{2}$ be real-valued functions of $x\in X$, and let $F(f_{1}(x),f_{2}(x))$ be monotonically increasing in $f_{1}$ and monotonically decreasing in $f_{2}$. Consider the optimization problem of finding $x^{*}\in\operatorname*{arg\,max}_{x}F(f_{1}(x),f_{2}(x)),$ and assume we are given $[b_{1},b_{2},...,b_{l},..,b_{L}]$ that span a wide enough range so that $x^{*}$ obeys: $b_{l^{*}-1}\leq f_{2}(x^{*})\leq b_{l^{*}}\textrm{ for some }l^{*}\in\{1,...,L\}$.

**(1)*For $x_{l}\in\operatorname*{arg\,max}\limits_{x:f_{2}(x)\leq b_{l}}F(f_{1}(x),b_{l})=\operatorname*{arg\,max}\limits_{x:f_{2}(x)\leq b_{l}}f_{1}(x),$ where the equality follows because $F$ monotonically increases in $f_{1}$, we have:

[TABLE]

(2)For $x_{l}\in\operatorname*{arg\,min}\limits_{x:f_{2}(x)\leq b_{l}}F(f_{1}(x),b_{l})=\operatorname*{arg\,min}\limits_{x:f_{2}(x)\leq b_{l}}f_{1}(x),$ we have:

[TABLE]

Proof.

Proof See Appendix Appendix B:. ∎

This theorem bounds the optimal value of $F$ along the whole regime of $x$ in terms of the values computed at the $L$ grid points. Note that the objective function of Formulations 1 and 2 is exactly of the form of Theorem 2, where

[TABLE]

[TABLE]

Note that monotonicity of $F$ in $f_{2}$ is ensured in the algorithm by the conditions on the sign of the treatment effect at lines 5 and 17 of the algorithm. The extra constraints on $\mathbf{a}$ in the ILP Formulation are compatible with Theorem 2. Thus, the bounds in (25) are direct results of Theorem 2 applied to the $z$ statistic as an objective function. Finally, minimization of $z$ can be achieved with Algorithm 6 by flipping the treatment indicator, i.e., calling control units as treated and treated units as control in the input to the algorithm, then running the procedure as is.

5 Approximate and Exact Distributions of Robust Test Statistics

As introduced in Section 2, our method aims to find an upper bound on the p-value of a desired test statistic over the set of good matches. While not directly used, a quantity relevant to our robust tests is the distribution of the test statistics under the maximal and minimal match assignments themselves, that is the distribution of $\Psi^{+}_{\mathcal{D}}=\max_{\mathbf{a}\in\mathcal{A}_{\textrm{good}}}\Psi_{\mathcal{D}}(\mathbf{a})$, and $\Psi^{-}_{\mathcal{D}}=\min_{\mathbf{a}\in\mathcal{A}_{\textrm{good}}}\Psi_{\mathcal{D}}(\mathbf{a})$ which are random variables denoting a generic test statistic ($\Psi$) under the match assignment that maximizes and minimizes it respectively. While we do not employ these distributions to derive p-values for our tests, we still choose to study them here to better understand the statistical behavior of our robust tests.

Specifically, we will study the distribution of the robust version of McNemar’s test taking into account the randomness due to the choice of match assignment, in order to provide analysts with exact distributions for the robust statistics. Without Exclusively Binning Constraints as defined in Section 3, the joint distributions can be extremely complex, but with these constraints, the calculation becomes much clearer. As a reminder, these are constraints that uniquely define a division of the units into strata, such that matches between two units in the same stratum always satisfy the quality constraints that define $\mathcal{A}_{\textrm{good}}$. To maintain exposition, we refrain from stating our results for this section formally (as the statements of these distributions are analytical but cumbersome), and instead give informal statements of our theorems. We refer to the appendix for rigorous versions of all our results and their proofs. This computation can be used to derive exact p-values for our robust McNemar’s test. Specifically, in the following sections, we show that analytical formulas exist for both null distributions we seek under Exclusively Binning Constraints, and that these formulas can be used to generate a lookup table for the distributions of $(\chi^{+},\chi^{-})$ in polynomial time.

5.1 Exact Randomization Distribution of $(\chi^{+},\chi^{-})$

Let us compute the finite sample joint distribution of $(\chi^{+},\chi^{-})$ under Fisher’s sharp null hypothesis of no treatment effect, which is defined as, for all $i$:

[TABLE]

We compute this distribution for the case in which constraints are exclusively binning and $\mathcal{S}$ can be constructed prior to matching. In this case, potential outcomes for all of the units in our sample are fixed and nonrandom: the only randomness in our calculations will stem from the treatment assignment distribution.

Our approach to constructing this randomization distribution without conditioning on the match assignment is as follows: we assume that any units that could be matched together under a good assignment have the same probability of receiving the treatment, and that, after treatment is assigned, matches are made to compute $(\chi^{+},\chi^{-})$ with the ILP in Formulation 1. This ILP is modified to include Exclusively Binning Constraints in addition to those in the formulation, and to make as many matches as possible within each of the strata defined by the constraints. In this setting, we assume that treatment assignment is random within the strata defined by the Exclusively Binning Constraint. This implies that Assumption 1 becomes:

Assumption 4**.**

(Stratified Conditional Ignorability) All units $i$ in stratum $l$ receive the treatment with fixed probability $e_{l}$. That is: for all $l=1,\dots,L$: $\Pr(T_{i}=1|X=x_{i},i\in S_{l})=\Pr(T_{i}=1|i\in S_{l})=e_{l}$, with $0<e_{l}<1$.

Let us introduce some notation to represent outcome counts within strata. We use $U_{l}$ to denote the number of treated units in stratum $l$ with outcome equal to 1, $V_{l}$ to denote the number of treated units with outcome 0, $\eta_{l}$, the number of control units with outcome 1, and $\nu_{l}$ the number of control units with outcome 0. Note that, by $\mathbb{H}_{0}^{sharp}$, we have the number of units with outcome 1 and 0 fixed in each stratum, thus, for stratum $l$, we use $N_{l}^{1}$ and $N_{l}^{0}$ to denote the number of units with outcome 1 and 0 in stratum $l$ respectively. Under Assumptions 2-4, these count variables follow the following generating process:

[TABLE]

One can view this data generation process as a results of our assumptions, including those of $\mathbb{H}_{0}^{sharp}$. The distribution of $(\chi^{+},\chi^{-})$ under $\mathbb{H}_{0}^{sharp}$ has the following properties:

Theorem 3**.**

(Randomization Distribution of $(\chi^{+},\chi^{-})$) Under Assumptions 1, 2, 4, and conditionally on $\mathbb{H}_{0}^{sharp}$:

• •

For a given stratum, $l$, the joint distribution of $B_{l}$ and $C_{l}$ (the counts of optimal matched pairs) can be expressed as a function of the distribution of the count variables $U_{l}$, $V_{l}$, $\nu_{l}$, $N_{l}^{t}$, and $N_{l}^{c}$ defined in (26)-(28).

• •

Given the signs of $(\chi^{+},\chi^{-})$, for any two strata, $l$ and $l^{\prime}$, $B_{l}$ and $C_{l}$ are independent of $B_{l^{\prime}}$ and $C_{l^{\prime}}$.

• •

The joint distribution of $(\chi^{+},\chi^{-})$ is the convolution of the joint distributions over $B_{l},C_{l}$ with each other.

A precise formulation for the distribution of test statistics in the theorem above is given in Theorem 7 in Appendix Appendix E:.

5.2 Exact Conditional Distribution of $(\chi^{+},\chi^{-})$

In this section, we seek a formulation for $\Pr(\chi^{-}=s,\chi^{+}=r|X)$, the distribution of robust McNemar’s statistics under exclusively binning constraints and a fixed number of treated units, $N_{l}^{t}$ in each stratum, but without conditioning on a specific match assignment. Potential outcomes are treated as random quantities in this case, and treated and control outcomes are assumed to be equal only on average. This distribution permits testing of hypotheses other than $\mathbb{H}_{0}^{sharp}$, and is more general than the randomization distribution derived in the previous subsection. Now denote with $N_{l}^{t}$ the number of treated units in the stratum $l$, and with $N_{l}^{c}$ the number of control units in that same stratum. In addition, let $p_{l}^{t}=\Pr(Y_{i}=1|i\in S_{l},T_{i}=1)$, and $p_{l}^{c}=\Pr(Y_{i}=1|i\in S_{l},T_{i}=0)$. For each stratum, $l=1,\dots,L$, the data are generated as follows:

[TABLE]

As we did before, we can use the distributions of these count variables to arrive at a formulation for the joint distribution of $(\chi^{+},\chi^{-})$ under the DGP just outlined.

Theorem 4**.**

(Conditional Distribution of $(\chi^{+},\chi^{-})$) Under Assumptions 1, 2, 4, and, for all strata, $l$, conditionally on $N_{l}^{t},N_{l}^{c}$:

• •

For a given stratum, $l$, the joint distribution of $B_{l}$ and $C_{l}$ (the counts of optimal matched pairs) can be expressed as a function of the distribution of the count variables $U_{l}$, $V_{l}$, $\nu_{l}$, $N_{l}^{t}$, and $N_{l}^{c}$ defined in (29)–(30).

• •

Given the signs of $(\chi^{+},\chi^{-})$, for any two strata, $l$ and $l^{\prime}$, $B_{l}$ and $C_{l}$ are independent of $B_{l^{\prime}}$ and $C_{l^{\prime}}$.

• •

The joint distribution of $(\chi^{+},\chi^{-})$ is the convolution of the joint distributions over $B_{l},C_{l}$ with each other.

A precise formulation for the distribution in the theorem above is in Theorem 8 in Appendix Appendix F:. Note that this distribution can be used to test the null hypothesis of no effect in each stratum by requiring $p_{l}^{t}=p_{l}^{c}$ for all $l$ in the general formulation we provide in the theorem.

5.3 Fast Computation of Lookup Tables for $(\chi^{+},\chi^{-})$

Another important advantage of the theorems given in the last two sections is that they permit us to generate lookup tables for the finite-sample null distributions of interest in polynomial time. Let $M=\sum_{l=1}^{L}M_{l}$ be the total number of matches made. A naïve procedure for computing a lookup table for the distribution in Theorems 3 and 4 for a given dataset based on brute force enumeration has computational complexity $O(2^{N^{M}+(N^{M})^{2}})$, because naively computing the domain of these distributions is a version of the subset-sum problem.

Fortunately, the complexity of computing the distribution of $(\chi^{+},\chi^{-})$ under both $\mathbb{H}_{0}$ and $\mathbb{H}_{0}^{sharp}$ can be significantly reduced because by Thms. 4 and 3, these distributions are convolutions of simpler distributions. We take advantage of this fact and of existing fast convolution algorithms to establish the following result:

Theorem 5**.**

There exists an algorithm that creates a probability table for the null distribution of $(\chi^{+},\chi^{-})$ given in Theorems 3 and 4 in $O(4MN^{4M}\log N)$.

We give proof of this theorem and outline the algorithm in question in the appendix. Note that the theorem directly implies that the worst-case time of computing the exact distributions of $(\chi^{+},\chi^{-})$ is polynomial in $N$. Figure 10 in Appendix Appendix H: shows marginal distributions of $(\chi^{+},\chi^{-})$ computed with the algorithm in Theorem 5, and shows that there is a location and scale difference between the distribution of $(\chi^{+},\chi^{-})$ when the null is true and when it is not. This demonstrates our tests’ capacity to detect full-sample treatment effects.

6 Simulations

We present a set of experiments on simulated datasets that demonstrated the usefulness and performance of our method in a series of settings. We compare the performance of both our robust tests against several popular matching methods and find that our proposed tests tend to display better performance (lower Type-I error) in several cases, without sacrificing statistical power.

In all our simulations, we pre-specify the number of units $N$, number of treated units $N^{t}$, and number of covariates $P$. We will vary $N$ throughout, however we will keep $N^{t}$ at 20% of the whole sample and $P=20$ unless otherwise specified (simulations with varying $P$ and constant $N$ are available in the appendix). We then generate data as follows, for $i=1,\dots,N$:

[TABLE]

Propensity scores are then simulated according to:

[TABLE]

where $\beta_{0}$ is a $P$-dimensional vector of coefficients, which are pre-specified before data generation. Treatment $t_{i}$ is then set to 1 for exactly $N_{t}$ units, such that each unit has $t_{i}=1$ with probability $e_{i}$. We then introduce two different data generating processes for our outcome variable:

[TABLE]

To simulate binary data for McNemar’s test we then set $y_{i}=\mathbb{I}[y_{i}^{*}>0]$, and to simulate continuous data for the Z-test we set $y_{i}=y_{i}^{*}$. The coefficient vectors $\beta_{0},\dots,\beta_{4}$, as well as $\tau$, are calibrated at each simulation so that the population ATT has a Cohen’s d-statistic (defined as: $\frac{\mathbb{E}[Y(1)-Y(0)|T=1]}{Var[Y(1)-Y(0)|T=1]}$) of a pre-specified value. In most settings, we will simulate data for $d=0,0.2,0.5,0.8,1.2,2$, which are commonly understood to correspond to null, tiny, weak, moderate, strong, and very strong treatment effects. For each setting, we simulate 1000 datasets and compute the proportion in which each method rejects the null hypothesis of 0 treatment effect at the 5% level as our estimate of the rejection rate of each method.

We compare against several commonly employed matching methods, as well as benchmark approaches such as a completely naïve test without any matching, and an idealized hypothesis test conducted on both potential outcomes, which is never possible in practice. We summarize the methods included in our simulation studies in Table 2.

To ensure the closest possible comparison, we constrain all methods to match all treated units to one control unit without replacement, that is $M=N^{t}$. Given that one of the main points of our paper is that there is not clear guidance on how to select hyperparameters for any given matching method, we try to use default or author-recommended values whenever possible for all the other methods tested. When this is not possible, we set hyperparameters to mimic those of the robust test as much as possible.

In all simulations, the set of good matches considered by our tests will be defined by all those matches that satisfy a mean balance constraint, as well as a caliper on the propensity score distance of the units matched. These constraints are implemented in our MIPs as follows:

• •

Mean Balance: $|\frac{1}{M}\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}a_{ij}(x^{t}_{ip}-x^{c}_{jp})|\leq\epsilon\times(\sigma(x_{p}^{t})/2+\sigma(x_{p}^{c})/2)$, where $\sigma_{p}(x_{p}^{t})$, $\sigma_{p}(x_{p}^{c})$ are the standard deviation of covariate $x_{p}$ in the treated and control set respectively. ($P$ constraints total).

• •

Caliper: $a_{ij}|e^{t}_{i}-e^{c}_{j}|\leq\epsilon$ for $i=1,\dots N^{t}$, $j=1,\dots,N^{c}$, where $e^{t}_{i}$ is the propensity score for treated unit $i$, and $e^{c}_{j}$ the analogue for control unit $j$ ($N^{t}\times N^{c}$ constraints total).

The value of $\epsilon$ is chosen to be the smallest feasible one on a grid between 0.01 and 1.

6.1 Results

We first present results from comparison of our robust McNemar’s test against the other methods. For all methods, p-values were computed after matching using the conventional formula for exact p-values for McNemar’s test (see Tamhane and Dunlop 2000).

Results are reported in Figure 2. We see that the robust test performs well as sample size grows: it displays low probability of rejection when the treatment effect is 0, and probability of rejection grows as the strength of the treatment also grows. When there is no ATT in the population, the robust test does display an error rate, but this is only when the sample size is low ($\leq$100), and the error rate decreases as the sample size grows, eventually converging with the one for the idealized McNemar’s test performed with the true potential outcomes. This is not true for methods we compare against, which seem to display an increase in error rate (incorrect rejection) as the sample size increases. When an ATT is present in the population, the robust test displays statistical power comparable to or better than that of other matching methods, indicating that our bounding approach is not too extreme, given that matches are chosen well.

Figure 3 reports results for the robust Z-test in the simulated data settings introduced previously. This version of our robust test performs even more strongly than the robust McNemar’s test, compared to other prominent matching methods. In this case, the robust z-test displays a good rejection rate even at low sample values ($N=50$) when the true ATT is 0, and is well powered as the ATT increases in strength, as well as when sample size grows.

7 Case Studies

In this section we apply our proposed methods to two real world datasets. We show that our test statistics can produce robust results in both datasets, which suggests that our methods can have wide practical applicability. Some of our results show that null hypotheses can be rejected robustly, while in some other cases they show that there is not enough evidence to reject the null hypothesis once the additional uncertainty from the choice of matching procedure is quantified via our robust tests.

7.1 Case Study 1: The Effect of Smoking on Osteoporosis in Women

In this case study we used GLOW (Global Longitudinal study of Osteoporosis in Women) data used in the study of Hosmer et al. (2013). Each data point represents a patient and the outcome is whether or not the person developed a bone fracture. The treatment is smoking. We match on several pre-treatment covariates: age, weight, height, and BMI. As there are several more control than treated units, we test $\mathbb{H}_{0}$ for the ATT with the general formulation of McNemar’s test introduced before, by matching almost all of the treated units.

Figure 4 shows results for testing $\mathbb{H}_{0}$ with the general program in Formulation 1, without binning constraints. We include a balance constraint in the formulation, by requiring that matched units $i,j$ respect $dist_{ij}\leq 0$, where $dist_{ij}$ is 0 if the sum of the differences of all the covariates is 6 or less and 1 otherwise. We choose the value 6 because this is the smallest caliper on the absolute distance between units that still permits all treated units to be matched. This figure lends evidence to the fact that $\mathbb{H}_{0}$ can be rejected robustly when the match assignment induces better balance in the data. We can conclude this because both min and max p-values are below 0.01 at each value of $M$, signifying that results would be statistically significant and positive using any good quality match assignment.

This case study shows the ability of our tests to detect effects when they are robust and present, even with a small number of observations.

7.2 Case Study 2: The Effect of Mist on Bike Sharing Usage

In our second case study, we used 2 years (2011-2012) of bike sharing data from the Capital Bike Sharing (CBS) system (see Fanaee-T and Gama 2014) from Washington DC. We study the effect of misty weather on the number of bikes rented. Control covariates include Season, Year, Workday, Temperature, Humidity and Wind Speed. Additional information on this case study is available in Appendix Appendix K:.

Figures 5 and 6 show the upper and lower bounds for the maximum objective function value for different $b_{l}$ with n=30 for the maximization problem in Figure 5 and minimization problem in Figure 6. These figures illustrate the meshes at different scales within the algorithm. We computed p-values for $Z^{+}$ and $Z^{-}$ under several counts of matched pairs: M=30, 50, 70, 90, 110. For $M=30$ through $M=90$, the p-value for $Z^{+}$ was 0 and the p-value for $Z^{-}$ was 1, while both p-values become 1 when the number of matches is 110. This is shown in Figure 7. The problem becomes infeasible for a larger number of matches. This illustrates that there is a lot of uncertainty associated with the choice of match assignment – that is, a reasonable experimenter choosing 90 matched pairs can find a p-value of $\sim$ 0 and declare a statistically significant difference while another experimenter can find a p-value of $\sim$ 1 and declare the opposite. In this case it is truly unclear whether or not mist has an effect on the total number of rental bikes. Figure 7 shows robust p-values for different amounts of matched units.

8 Conclusion

Believing hypothesis test results conducted from matched pairs studies on observational data can be perilous. These studies typically ignore the uncertainty associated with the choice of matching method, and in particular, how the experimenter or an algorithm chooses the matched groups: we have given both simulated and real-data evidence of this problem. We want to know that for any reasonable choice of experimenter who chooses the assignments, the result of the test would be the same. In this work, we have addressed the issue above by introducing robust test statistics that consider extrema over all possible good matches. This is justified because p-values obtained for the extrema must, by definition, include p-values obtained under any other possible good match. We have provided practical implementations of this principle for both discrete and continuous data, as well as theoretical and empirical analyses of the performance of our methods.

Acknowledgement The authors express their gratitude to the Natural Sciences and Engineering Research Council of Canada (NSERC) for partial financial support of this research.

Appendices

Appendix A: Constraining the Quality of the Matches

Determining the quality requirements that matches should satisfy to be usable is ultimately up to the experimenter, however there are several general types of constraints that matches in the uncertainty set $\mathcal{A}_{\textrm{good}}$ should most often obey. Let $dist(x,x^{\prime})$ be a metric on the space of $x$, some of these constraints are:

• •

(Calipers) When $\mathbf{a}(i)\neq\emptyset$ then $\textrm{dist}(x_{i}^{t},x_{\mathbf{a}(i)}^{c})\leq\epsilon$.

• •

(Covariate balance, mean of chosen treatment units similar to mean of control group) $\forall\textrm{ covariates }p$ we have:

[TABLE]

• •

(Maximizing the fitness of the matches) In general, one can optimize any measure of user-defined fitness for the assignment, and then constrain $\mathcal{A}_{\textrm{good}}$ to include all other feasible assignments at or near that fitness level, by including the following constraints:

[TABLE]

where Maxfit is precomputed as: $\textrm{Maxfit}=\max_{\mathbf{a}\in\mathcal{A}}\textrm{Fitness}(\mathbf{a},\{x_{i}^{t}\}_{i=1}^{N^{t}},\{x_{i}^{c}\}_{i=1}^{N^{c}}).$ If one desires the range of results for all maximally fit pairs and no other pairs, $\epsilon$ can be set to 0.

We note that, to constitute a ILP formulation, the chosen type of fitness function must be expressed as some linear function of the matching indicators. While it is technically true that not all fitness functions can be encoded in this way, many popular constraints for matching can be. For example, Zubizarreta [2012] gives linear constraint formulations for a number of popular match quality constraints. However, despite the fact that many match assignments may exist that satisfy these constraints, these methodologies produce only one, and it is not chosen uniformly at random from any known distribution.

In what follows, we provide special cases of the Robust Procedure for two specific hypothesis tests, McNemar’s test for binary outcomes and the $z$-test for real-valued outcomes. We outline strategies for computing these statistics as well as their distributions under the hypotheses of interest. We begin with McNemar’s test.

Appendix B: Proof of Theorem 2

Proof.

Proof We prove the result for statement (1), as the proof of statement (2) is exactly symmetrical. Because the $b_{l}$’s are defined on a pre-specified grid, we know the maximum value of $F$ may not occur at one of the grid points. Since by definition of $x_{l}$ we know that $f_{2}(x_{l})\leq b_{l}$, and since $F$ is decreasing in its second argument, we have $F(f_{1}(x_{l}),b_{l})\leq F(f_{1}(x_{l}),f_{2}(x_{l}))$ for each $l$, and taking a max over all $l$:

[TABLE]

This is the left inequality of the bound. The rest of the proof deals with the right inequalities. First, it is true that:

[TABLE]

If this were not true then either $\displaystyle f_{1}(x^{*})<\max_{x:f_{2}(x)=f_{2}(x^{*})}f_{1}(x)$ or $\displaystyle f_{1}(x^{*})>\max_{x:f_{2}(x)=f_{2}(x^{*})}f_{1}(x)$. If the first inequality were true then

[TABLE]

which contradicts the definition of $x^{*}$. The second option also cannot be true as we know there exists a solution $x^{*}$ so that the maximum is attained with $f_{1}$ and $f_{2}$ values $f_{1}(x^{*})$ and $f_{2}(x^{*})$. So we can say that (32) holds.

From (32) and using $l^{*}$ defined in the statement of the theorem, we can derive:

[TABLE]

where we used that the set $\{x:f_{2}(x)=f_{2}(x^{*})\}$ is smaller than the set $\{x:f_{2}(x)\leq f_{2}(x^{*})\}$ which is smaller than $\{x:f_{2}(x)\leq b_{l^{*}}\}$, since $f_{2}(x^{*})\leq b_{l^{*}}$ by definition of $l^{*}$. Thus, $f_{1}(x^{*})\leq f_{1}(x_{l^{*}})$. Now,

[TABLE]

Here the first inequality above follows from the definition of $l^{*}$, $b_{l^{*}-1}\leq f_{2}(x^{*})$, and the fact that $F$ decreases in the second argument. The second inequality comes from (33) and the fact that $F$ is increasing in its first argument. The third inequality follows from taking a maximum over all $l$ rather than using $l^{*}$. The proof is complete. ∎

Appendix C: Algorithms that maximize and minimize $\chi$ under Exclusively Binning Constraints

In the following sections, we offer exact formulas and computational algorithms for the distribution of the test statistic values obtained with the results in this section, and when matches are considered as part of the statistic.

If the constraints on $\mathcal{A}_{\textrm{good}}$ are exclusively binning, then $\chi$ can be optimized quickly and easily in both directions if the partition $\mathcal{S}$ is constructed first. Before showing this, we introduce an IP formulation for McNemar’s test with exclusively binning constraints.

As a reminder, we are in the following situation: we have $N^{t}$ treated units and $N^{c}$ control units measured on P features, $X$ that take value in a finite set, $\mathcal{X}$. All these units also have an outcome $Y\in\{0,1\}$ and treatment $T\in\{0,1\}$. Since the constraints on the optimization problem are exclusively binning, then based on the values of $X$ we group the units into $L$ strata such that each stratum $S_{1},\dots,S_{L}$ contains $N_{l}^{t}$ and $N_{l}^{c}$ units. In each stratum there will be $U_{l}$ treated units with outcome $Y_{i}=1$, $\eta_{l}$ control units with outcome $Y_{j}=1$, $V_{l}=N_{l}^{t}-U_{l}$ treatment units with outcome 0, and $\nu_{l}=N_{l}^{c}-\eta_{l}$ control units with outcome 1. We then would like to create pairs of units within each stratum, such that each pair contains exactly one treated unit ($T_{i}=1$) and one control unit ($T_{j}=0$). Once we have created the pairs we compute:

• •

$A_{l}$: the number of matched pairs in stratum $l$ such that the treated unit has $Y=1$ and the control unit has $Y=0$. We refer to these pairs as $A$-pairs.

• •

$B_{l}$: the number of matched pairs in stratum $l$ such that the treated unit has outcome 1 and the control unit 0. We refer to these pairs as $B$-pairs.

• •

$C_{l}$, the number of matched pairs in stratum $l$ such that the treated unit has outcome [math] and the control unit 1. We refer to these pairs as $C$-pairs.

• •

$D_{l}$, the number of matched pairs in stratum $l$ such that both the treated and control units have outcome [math], referred to as $D$-pairs.

The following tables summarize the data and the statistics we are interested in:

We then take the sums of $B_{l}$ and $C_{l}$ across all strata to obtain $B=\sum_{l=1}^{L}B_{l}$ and $C=\sum_{l=1}^{L}C_{l}$. Finally, we use these two quantities to compute:

[TABLE]

We would like to make the matches within each stratum such that $\chi$ is either maximized or minimized, assuming that we must match as many units as possible. Throughout the rest of this document we use a $+$ superscript to denote values of $A$, $B$, $C$, $D$ and $\chi$ output by maximizing $\chi$ and a $-$ superscript to denote the corresponding values obtained by minimizing it. We limit our analysis to the case in which the maximum number of feasible matches needs to be achieved, that is, exactly $M=\sum_{l=1}^{L}\min(N_{l}^{t},N_{l}^{c})$ must to be made:

Formulation 2: IP formulation for McNemar’s test with exclusively binning constraints

[TABLE]

subject to:

[TABLE]

Note that analysts can remove this part of the formulation and introduce a fixed number of matches as constraint in a way similar to Constraint (9) if desired.

This problem can be solved in linear time and without one of the canonical IP solution methods by using the fact that the strata defined by the exclusively binning constraints can each be optimized separately, once the direction of the resulting statistic is known. In stratum $S_{l}$ there will be $U_{l}$ treated units with outcome $Y_{i}=1$, $\eta_{l}$ control units with outcome $Y_{j}=1$, $V_{l}=N_{l}^{t}-U_{l}$ treatment units with outcome 0, and $\nu_{l}=N_{l}^{c}-\eta_{l}$ control units with outcome 1. This is summarized in Figure 8.

To ensure that as many units as possible are matched, within each stratum we make exactly $M_{l}=\min(N_{l}^{t},N_{l}^{c})$ matches. We would like to make the matches within each stratum such that $\chi$ is either maximized or minimized. Algorithm 1 maximizes $\chi$,

Appendix C: .1 Correctness of the Optimization Algorithms

As it is clear from their definitions, these algorithms do not require any of the conventional MIP solving techniques and as such are much faster: the running time of Algorithm 1 is clearly linear in $N$, the number of units. This constitutes a substantial speed up over solving the problem with a regular MIP solver. The following theorem states the correctness of the algorithm for solving Formulation 2:

Theorem 6**.**

(Correctness of Algorithm 1) Algorithm 1 globally solves Formulation 2 for the Maximal value of $\chi$.

What follows is a proof of the algorithms’ correctness. The proof is structured into 4 different claims and a theorem equivalent to Theorem 6 following directly from these claims. Before introducing these claims it is useful to summarize the various ways in which units can be matched and unmatched to lead to different pair types: the cells in Table 3 show what pairs can be created by unmaking two other pairs and matching their units across. For example, if we unmake an $A$-pair and a $D$-pair, we are left with a treated and a control unit with $Y=1$ and a treated and a control unit with $Y=0$: if we match them across we obtain a $B$ and a $C$-pair. In the rest of this document, we refer to this operation of unmaking two pairs and matching across their units as exchanging pair one with pair two.

Note also that, if the number of treated units to be matched equals the number of control units, and all units are matched in some way, the only way we can change those matches is by performing one or more of the operations detailed in the table. This is, of course, only possible if the required pairs are present among the existing matches for example, we cannot unmake a $A$-pair if there are none made already. For claims 1-3, assume that any two units can be matched together. Since we don’t explicitly consider the different strata in these claims, we omit the $l$ subscript from the notation.

Claim 1**.**

Suppose we want to maximize $\chi$ subject to the constraint that we must make as many matches as possible, that is we must make $M=\min(N^{t},N^{c})$ matches, and suppose that $N^{c}>N^{t}$ so that $M=N^{t}$. Then it is always optimal to leave unmatched $\min(N^{c}-N^{t},\eta)$ control units with outcome 1 and $\max(N^{c}-N^{t}-\eta,0)$ control units with outcome 0. Suppose instead that $N^{t}\geq N^{c}$, then it is always optimal to leave unmatched $\min(N^{t}-N^{c},V_{l})$ treated units with outcome 0 and $\max(U_{l}-N_{l}^{c},0)$ treated units with outcome 1.

Proof.

Proof. We will show that, if $N^{c}>N^{t}$ and we must make $N^{t}$ matches, then it is optimal to first leave unmatched as many control units with outcome 1 as possible, that is $\min(N^{c}-N^{t},\eta)$, and, if after these units have been excluded, there still are more control than treated units, to leave unmatched the remaining control units with outcome 0. Suppose initially that $N^{c}=N^{t}+1$, fix $\chi=\frac{B-C-1}{\sqrt{B+C}}$ and assume that there are exactly $N^{t}-1$ already matched treatment units and at least two leftover control units. Since we must make exactly $N^{t}$ matches we can only match the leftover treatment unit with one of the two controls.

Assume that the two control units are $u_{1}$ with outcome $Y_{u_{1}}=1$ and $u_{0}$ with $Y_{u_{0}}=0$. There are two possible scenarios: first, the currently unmatched treatment unit has outcome 1: if we match it with $u_{0}$ we get a $B$-pair with corresponding value of $\chi$:

[TABLE]

if we match the treatment unit with $u_{1}$ we get a $A$-pair, and:

[TABLE]

With some algebra we can see that, $\chi_{u_{0}}\geq\chi_{u_{1}}$, which implies that we always gain more from matching a treated unit with outcome 1 to a control unit with outcome 0. In this case, leaving $u_{1}$ unmatched is the optimal choice.

Now suppose that the leftover treatment unit has outcome 0. Then, if we match it with $u_{0}$ we form a $D$-pair and get:

[TABLE]

no change from the initial $\chi$. If we match this treatment unit with $u_{1}$ we have formed a $C$-pair instead get:

[TABLE]

again some algebra reveals that $\chi_{u_{0}}\geq\chi_{u_{1}}$: the value of $\chi$ is maximized by choosing to match the treatment unit with $u_{0}$ in this case as well. This shows that, when there is a choice of multiple control units to match with a treatment unit the control with outcome 1 should always be left out if we wish to maximize $\chi$.

We now show that, if $N^{t}=N^{c}+1$ then it is always optimal to leave unmatched a treated unit with outcome 0 instead of one with outcome 1 if $\chi$ is to be maximized. Let there be one unmatched control unit and two candidate treatment units for it to be matched with, $u_{1}$ such that $Y_{u_{1}}^{t}=1$ and $u_{0}$ such that $Y_{u_{0}}=0$. Let the value of $\chi$ without those units be $\chi=\frac{B+C-1}{\sqrt{B+C+1}}$. If the unmatched control has outcome 1 and we match it to $u_{0}$ we end up with a $C$-pair and an updated value of $\chi$ equal to:

[TABLE]

If we instead match the control unit with outcome 1 to $u_{0}$ we produce an $A$-pair, and the following value of $\chi$:

[TABLE]

Then, with some algebra we can see that $\chi_{u_{0}}\leq\chi_{u_{1}}$, implying that, in order to maximize $\chi$, the optimal strategy is to leave $u_{0}$ unmatched. If the control unit to be matched has outcome 0 and we match it with $u_{0}$ we get a $D$-pair and $\chi_{u_{0}}=\chi.$ If we instead match the control unit with outcome 0 to $u_{1}$ we get a $B$-pair and:

[TABLE]

Again, with some algebra we see that $\chi_{u_{1}}\geq\chi_{u_{0}}$, implying that leaving $u_{0}$ unmatched is optimal in this case as well. This shows that leaving treated units with outcome 0 unmatched if $N^{t}=N^{c}+1$ is the optimal strategy to maximize $\chi$.

For the case in which $N^{t}=N^{c}+k$ we can proceed by induction on the number of matched control units: assume inductively that the optimal choice at $k-1$ is to leave unmatched a treated unit with outcome 0. At match $k$ we will have a value $\chi_{k}$ and one unmatched control unit and two candidate treatments. By the above, the value of $\chi_{k}$ that we would get by leaving the treated unit with outcome 0 unmatched is always larger than the one we would get by leaving the unit with outcome 1 unmatched. Because of this it is optimal to leave the treated unit with outcome 0 unmatched also for the $k^{th}$ match. This proves that, regardless of the difference between the number of treated and control units, it is optimal to leave as many treated units with outcome 0 as possible unmatched. If $N^{t}\geq N^{c}$, the largest possible number of treated units with outcome 0 that can be left unmatched is clearly $\min(N^{t}-N^{c},V)$, either we exhaust the difference between $N^{t}$ and $N^{c}$ by not matching treatment units with outcome 0, or we exhaust all $V$ treated units with outcome 0 and still have leftover treated units that don’t have a match in the control group. In this second case, we will have to make up this difference by leaving unmatched treated units with outcome 1 in excess: the precise amount of which is $N^{t}-N^{c}-V=N^{t}-N^{c}-N^{t}+U=U-N^{c}$. A symmetrical argument shows that, if $N^{c}=N^{t}+k$ it is always optimal to leave as many control units with outcome 1 as possible unmatched, before leaving units with outcome 0 unmatched, and that the largest possible amount of control units with outcome 1 that can be left unmatched in this case is $\min(N^{c}-N^{t},\eta)$, and the amount of control units with outcome 0 is $N^{c}-N^{t}-\eta$, in case all $\eta$ control units with outcome 1 are left unmatched without being able to exhaust the difference between unmatched units.

Finally, note that the proof for the case in which we wish to minimize $\chi$ is exactly symmetrical to this one.∎

Claim 2**.**

If $N^{t}=N^{c}=M$, and exactly $M$ matches must be made, then $B^{+}-C^{+}=B^{-}-C^{-}=U-\eta$ independently of how units are matched.

Proof.

Proof. Let $W$ be the matching in which all units are paired in a way such that $B=\min(U,\nu)$ and $C=\min(\eta,V)$. Since $M=U+V=\eta+\nu$, then $U\geq\nu\iff\eta\geq V$, so it must be that either: $B=U,C=\eta$ or $B=V,C=\nu$. By definition of $V$ and $\nu$ we have that in both cases: $B-C=U-\eta$. Now we will prove that this equality must hold for any other match in which all units are matched. Let $B$ be the number of $B$-pairs created with that matching and $C$ the number of $C$-pairs. Consider any other match $W^{\prime}\neq W$ also satisfying the fact that all $M$ treatment and control units are matched and let $B^{\prime}$ and $C^{\prime}$ be the counts of $B$ and $C$-pairs generated by $W^{\prime}$. Since exactly the same number of units are matched in $W$ and $W^{\prime}$ it must be that there exists some sequence of exchange operations that, if applied to $W$ generates $W^{\prime}$. We now proceed by induction on $k$, the number of operations applied to $W$ to get to $W^{\prime}$, starting with $k=1$. Note first, by Table 3 that the only operations that can alter the number of $B$ and $C$-pairs are: exchanging an $A$-pair with a $D$-pair, obtaining a $B$-pair and a $C$-pair, and exchanging a $B$-pair with a $C$-pair, obtaining an $A$-pair and a $D$-pair. In the first case, we unmake an $A$-pair and a $D$-pair in $W$ and use those units to make a $B$ and a $C$-pair in $W^{\prime}$ so the respective counts are now $B^{\prime}=B+1$ and $C^{\prime}=C+1$, which implies that $B^{\prime}-C^{\prime}=B-C$. In the other case we have $B^{\prime}=B-1$, $C^{\prime}=C-1$ and $B^{\prime}-C^{\prime}=B-C$. Finally, suppose inductively that $B^{(k)}-C^{(k)}=B-C$ after the $k$th exchange operation. Again, by Table 3 the only operations that can alter the counts of $B$ and $C$ are exactly the two discussed in the base case ; and by the same reasoning, we conclude that $B^{(k+1)}-C^{(k+1)}=B-C$. ∎

Claim 3**.**

If $N^{t}=N^{c}=M$ and exactly $M$ matches must be made, then Algorithms 2 and 3 respectively make the matches that maximize and minimize $B+C$.

Proof.

Proof. Consider Algorithm 3 first. We will only show the statement for this algorithm as the proof of the correctness of the other algorithm is exactly symmetrical to this. Let $W^{*}$ be the match output by this algorithm and let $A^{*},B^{*},C^{*},D^{*}$ be the respective pair counts under $W^{*}$. By lines 1-4 of Algorithm 3 we know that:

[TABLE]

By definition of $U,V,\eta,\nu$, $A^{*}$ and $D^{*}$ are the largest possible number of $A$ and $D$-pairs that can be made with the given units. It can be seen from the definitions above that one of $C^{*}$ or $B^{*}$ will always be 0: in case $A^{*}=U$ then $U-A^{*}=0$ and $B^{*}=\min(U-A^{*},\nu-D^{*})=0$. In case $A^{*}=\eta$ we have $\eta-A^{*}=0$ and $C^{*}=\min(\eta-A^{*},V-D^{*})=0$ by consequence. Since exactly $M$ units must be matched, and $N^{t}=N^{c}=M$ by assumption, then the only operations allowed to change the matches from $W^{*}$ are those in table 3, as there are no leftover units unmatched. Note that the only operation in the table that would allow for a decrease in $B$ and $C$ is exchanging$B$ with $C$, but this operation can never be performed on $W^{*}$ as either it contains no $C$-pairs or no $B$-pairs. Then it must be that $B^{*}$ and $C^{*}$ are the smallest number of $B$ and $C$-pairs that can be made with the existing units and thus that they minimize $B+C$. ∎

Appendix C: .1.1 Proof of Theorem 6

Proof.

Proof. Consider first the problem of maximizing $\chi$. The observations are divided into $l$ strata such that $N_{l}^{t}$ treated units and $N_{l}^{c}$ control units are in each stratum: we can only match units that are in the same stratum. If we denote the count of $B$-pairs and $C$-pairs in stratum $l$ with $B_{l}$ and $C_{l}$ then the objective function for the problem is

[TABLE]

As $TE$ is a separable function of $B_{l}$ and $C_{l}$, optimizing the former equals maximizing each of the latter individually, the same is true for $SD$. Note first, that, by the form of the objective function, we must make as many matches as possible. This number of matches is exactly $M=\sum_{i=1}^{L}M_{l}=\sum_{i=1}^{L}\min(N_{l}^{t},N_{l}^{c})$ because, by constraints (10) and (11), each treatment and control unit can only be matched once. By Constraint (38) we know that exactly $M_{l}=\min(N_{l}^{t},N_{l}^{c})$ matched pairs must be constructed in each stratum, therefore, units in excess of $M_{l}$ must be discarded in each stratum. By Claim 1 if there are more control units to be matched than treated units in one stratum it is always optimal to leave unmatched the amounts of units described in Claim 1 in each stratum, independently of how matches are made in other strata by separability of $TE$. The algorithm does this explicitly when defining $U_{l}^{+},V_{l}^{+},\eta_{l}^{+},\nu_{l}^{+}$ at lines 2-5. These are updated counts of units to be matched, such that, for all $l$, $U_{l}^{+}+V_{l}^{+}=\eta_{l}^{+}+\nu_{l}^{+}=M_{l}$ by definition of these quantities, therefore, the count of treated and control units to be matched is equal in all strata and equal to $M_{l}$. Note that definitions of the counts given in lines 2-5 of Algorithm 1 are precisely the initial number of each outcome-treatment pair minus the optimal amount of units to be left unmatched given in Claim 1. By Claim 2 we know that if the number of treated units to be matched is equal to the number of control units to be matched, as it is the case after line 5 of the algorithm, then $TE=\sum_{l=1}^{L}U_{l}^{+}-\eta_{l}^{+}$ independently of how matches are made. Because of this, $TE$ can be considered fixed at this point, and maximizing $\chi$ equates with maximizing $SD$ if $TE<0$ and minimizing it if $TE\geq 0$: this is checked explicitly by the algorithm at line 8. Lastly, if the algorithm calls $ComputeMinimizedSD$ if $TE\geq 0$ and $ComputeMinimizedSD$ if $TE<0$: these two procedures are shown to correctly maximize and minimize $B_{l}+C_{l}$ in each stratum by Claim 3. Note that $SD=\sum_{l=1}^{L}B_{l}+C_{l}$ is also separable in the strata, and therefore it can be optimized globally by separately optimizing $B_{l}+C_{l}$ in each stratum. Since $\sqrt{SD}$ is monotonic in $B_{l}+C_{l}$, we know that this quantity is also maximized or minimized in this way. This shows that Algorithm 4 globally maximizes $\chi$ over the set of allowed matches. Finally, note that constraints (7), and (8) are all not violated by definition of these quantities, and constraints (10) and (11) are not violated because no two units are matched more or than once. Finally, the additional exclusively binning constraints are obeyed by definition, as matches are made exclusively within each stratum. The proof of the correctness of Algorithm 5 is exactly symmetrical to this one: this symmetry is made apparent by the fact that running Algorithm 5 on the data is equivalent to flipping the treatment indicator and then running Algorithm 4 on the resulting data. ∎

Appendix D: Algorithm for Maximization of $Z$ statistic

Algorithm 6 is useful for maximization of $Z$. Minimization of $Z$ can be accomplished with a symmetric algorithm.

Appendix E: Randomization Distribution of $(\chi^{+},\chi^{-})$ Under Exclusively Binning Constraints

In this section we use the algorithm above to derive a randomization distribution for $(\chi^{+},\chi^{-})$ under the null hypothesis of no treatment effect. This allows us to test $\mathbb{H}_{0}^{sharp}$ under two types of uncertainty: in the data itself, and in the choices made by the analyst in a MaE model. As a reminder, under $\mathbb{H}_{0}^{sharp}$, potential outcomes are assumed to be fixed and invariant between treatment regimes. Under exclusively binning constraints, the units are divided into $S_{1},\dots S_{l}$ strata, and matches are allowed only in the same stratum.

Given $L$ different levels in which the covariates are grouped, and that there are $N_{l}^{t}$ treatment units and $N_{l}^{c}$ control units in each stratum $l$. From the two assumptions before it follows that the data are generated as follows:

[TABLE]

As stated in the paper, all these quantities have interpretations in our matching framework. Specifically, within one stratum, I.E., at one level of $X$, called $x_{l}$: $U_{l}$ is the number of treated units with outcome 1,$V_{l}$ is the number of units with $T=1$ and $Y=0$, $\eta_{l}$ is the number of units with $T=0$ and $Y=1$ and $\nu_{l}$ is the number of units with $T=0$ and $Y=0$. The null hypothesis of no treatment effect is encoded in the fact that the distributions of $U_{l}$ and $\eta_{l}$ differ only in the number of trials and not in the probability of a success. Note now that these pair counts are random variables: the algorithms make pairs to purposefully obtain optimal values of $\chi$, which in turn depends on the random data. Recall that $B$ and represents the number of pairs such that the treated unit has outcome 1 and the control unit has outcome 0 and $C$ the total number of pairs in which the control unit has outcome 1 and the treated unit 0.

Recall that the test statistic is defined as follows:

[TABLE]

where $B^{+}$ is the count of matched pairs produced by Algorithm 4 such that the treated unit has outcome 1 and the control unit in the pair has outcome 0, $C^{+}$ is the count of pairs where the opposite is true, and $B^{-}$ and $C^{-}$ are the analogues produced by the minimization algorithm. For convenience, we also introduce ”truncated” versions of the variables above, letting $G_{l}^{-}=\max(N_{l}^{t}-N_{l}^{c},0)$ and $G_{l}^{+}=\max(N_{l}^{c}-N_{l}^{t},0)$:

[TABLE]

These definitions correspond to those introduced at lines 2-5 of Algorithms 4 and 5. Figure 9 summarizes the variables in our framework as well as how the maximization and minimization algorithms operate.

Appendix E: .1 Simplified Representations for McNemar’s Statistic

The following statement gives us simplified forms for the values of $TE$ and $SD$ output by the maximization and minimization algorithms respectively, we state it as a claim as it requires a simple yet nontrivial proof:

Claim 4**.**

Let $TE^{+}$, $SD^{+}$ and $TE^{-}$ and $SD^{-}$ represent the values of the numerator and denominator of $\chi$ as defined in Eq (47), and output by algorithms 4 and  5 respectively. Let $M_{l}=\min(N_{l}^{t},N_{l}^{c})$ denote the number of matches that are made in each stratum. Then they can be written as follows:

[TABLE]

Proof.

Proof. Note first that the definition of $TE^{+}$ and $TE^{-}$ follows directly from lines 3-5 and 8 of Algorithms 4 and 5 respectively. As for the definition of $SD$ we will prove only the claim for $SD^{+}$ as the proof for $SD^{-}$ is exactly symmetrical. Proving the claim implies showing that, denoting with $B^{+}$ and $C^{+}$ the count of $B$ and $C$ pairs that maximize $\chi$:

[TABLE]

Now define:

[TABLE]

It is clear from this definition that $SD^{+}=\sum_{l=1}^{L}SD_{l}^{+}$, and since $SD_{l}^{+}$ has the same definition in all strata, it suffices to prove that $B_{l}^{+}+C_{l}^{+}$ can be written as in Eq. (60) for one stratum to prove the equality in Eq. (58). The rest of the proof is concerned with establishing this result.

Consider the computation for $SD^{+}$ occurring in Alg. 4: the algorithm checks explicitly if $TE\geq 1$ at line 8 and, if true it calls Alg. 3 on inputs $(U_{l}^{+},V_{l}^{+},\eta_{l}^{+},\nu_{l}^{+})$ to generate the matched pair counts $A_{l}^{+},B_{l}^{+},C_{l}^{+},D_{l}^{+}$, if false it generates the same quantities by calling Alg. 2 on the same inputs.

Consider now the case in which $TE^{+}\geq 1$, and Algorithm 3 is called with $(U_{l}^{+},V_{l}^{+},\eta_{l}^{+},\nu_{l}^{+})$ as inputs. We know that the algorithm returns the following counts for $B$ and $C$ pairs:

[TABLE]

Case 1 for $TE^{+}\geq 1$: $U_{l}^{+}\geq\eta_{l}^{+}$.

Note first that, in this case:

[TABLE]

Because of this we can see that:

[TABLE]

Case 2 for $TE^{+}\geq 1$: $U_{l}^{+}<\eta_{l}^{+}$.

Note that in this case Eq. (65), implies that: $\nu_{l}^{+}<V_{l}^{+}$. Because of this we have:

[TABLE]

This shows $SD_{l}^{+}=S_{l}^{+}$ in the case in which $TE\geq 1$, the first of Eq. (58).

The proof is similar for the case in which $TE<1$. Now we must show that, in this case $SD_{l}^{+}=R^{+}$. In this case Algorithm 2 is called with $U_{l}^{+},V_{l}^{+},\eta_{l}^{+},\nu_{l}^{+}$ as inputs at line 16 of Algorithm 4. On those inputs, Algorithm 2 will return the following counts of $B$ and $C$ pairs:

[TABLE]

We proceed as before with two separate cases.

Case 1 for $TE<1$: $U_{l}^{+}\geq\nu_{l}^{+}$.

First note that:

[TABLE]

Second, the assumption of this case also implies:

[TABLE]

Putting these results together we obtain:

[TABLE]

Case 2 for $TE<1$: $U_{l}^{+}<\nu_{l}^{+}$.

This assumption of this case together with (68) implies $\eta_{l}^{+}<V_{l}^{+}$. Note also that:

[TABLE]

With the two results above we obtain:

[TABLE]

This proves that $SD_{l}^{+}=R_{l}^{+}$ in the case in which $TE^{+}<1$, the second case in Equation (58). ∎

Appendix E: .2 Lemma 1

Lemma 1**.**

(Randomization Distribution of Truncated Variables) For $l=1,\dots,L$ let $(N_{l}^{1},N_{l}^{0},e_{l})$ be fixed and known and let $\mathbb{D}=(U_{l},V_{l},\eta_{l},\nu_{l},N_{l}^{t},N_{l}^{c})$ be drawn from the data generating process of equations (41) – (46). Let $a,b,c,d,m$ be elements of $\{0,\dots,\min(N_{l}^{0},N_{l}^{2})\}$ and let $m\in\{0,\dots,N_{l}\}$. The variables ($U_{l}^{+}$, $U_{l}^{-}$, $\eta_{l}^{+}$, $\eta_{l}^{-}$, $M_{l}$) have the following joint distribution:

[TABLE]

where:

[TABLE]

and $\Pr(U_{l}=j|\mathbb{H}_{0}^{sharp})=Bin(j,N_{l}^{1},e_{l})$ and $\Pr(V_{l}=k|\mathbb{H}_{0}^{sharp})=Bin(k,N_{l}^{0},e_{l})$.

Note that even though the domain of the distribution above cannot be expressed by a simpler formula, it is simple and computationally fast to enumerate it for finite $N_{l}$.

Proof.

Proof Throughout the proof we maintain Assumption 4 and all of the probability statements to follow are conditional on it holding. The proof is simple and follows by inspecting the definitions of the truncated variables. Note first that the form for $\Pr(U_{l}^{+}=a,U_{l}^{-}=b,\eta_{l}^{+}=c,\eta_{l}^{-}=d,M_{l}=m|\mathbb{H}_{0}^{sharp})$ given in Eq. (84) follows from the law of total probability, independence of $U_{l}$ and $V_{l}$ and conditional independence of $U_{l}^{+},V_{l}^{+},\eta_{l}^{+},\eta_{l}^{-},M_{l}$ given $U_{l}$ and $V_{l}$. It remains to show that the forms for the indicator functions in Equations (85)–(89) are those in the theorem. This can be done by inspecting the definitions of the truncated variables and by expanding them out into conditions on $U_{l}$ and $V_{l}$. In what follows we refer to the following definitions of the quantities employed, listed here with reference to where they are introduced in the paper:

[TABLE]

We now derive conditions on $U_{l}$ and $V_{l}$ that lead to the realization of the event $U_{l}^{+}=a$, we do so by expanding the definition of $U_{l}^{+}$ and considering all the cases it entails separately. Starting with the definition of $U_{l}^{+}$ we have:

[TABLE]

There are now two cases for the event $U_{l}^{+}=a$, depending on how the max function in the definition of $U_{l}^{+}$ is resolved:

Case 1 for $U_{l}^{+}=a$:

[TABLE]

Because of this our event of interest can be written as:

[TABLE]

which is the first condition in the first case of Equation (85). Note that the condition of this case given in Eq. (91) implies that $2U_{l}+V_{l}-N_{l}\geq 0$, and we can use the result in (92) to expand this condition as follows:

[TABLE]

which is the second condition in the first case of Equation (85).

Case 2 for $U_{l}^{+}=a$:

[TABLE]

Using this case and the result in (90), we can write the event of interest as

[TABLE]

the first condition in the second case of Equation (85). Second, the fact that the $\max$ function equals 0 in Case 2, implies that $2U_{l}+V_{l}-N_{l}\leq 0$ is another condition for this case. Using the result in (94) this can be simplified as:

[TABLE]

which is the second condition in the second case of Equation (85).

We proceed in the same manner for the event $U_{l}^{-}=b$. In Equation (86) we start by expanding the definition of $U_{l}^{-}$:

[TABLE]

Because of the two $\max$ functions we have four different possibilities for the value of $U_{l}^{-}$, all dependent on how the $\max$ functions resolve. As we did before, we can derive conditions on $U_{l}$ and $V_{l}$ by studying these four cases separately. We start each of the four cases by listing the ways the inner max and the outer max resolve.

Case 1 for $U_{l}^{-}=b$:

[TABLE]

First, we can use both conditions to simplify the event $U_{l}^{-}=b$ as follows:

[TABLE]

which is the first condition in the first case of Equation (86). Note now that the condition in (97) implies that $U_{l}-\max(2(U_{l}+V_{l})-N_{l},0)\geq 0$, combining this with the above we obtain:

[TABLE]

so this condition is always satisfied because of how we restrict the domain of $b$. The condition in (96) implies that $2(U_{l}+V_{l})-N_{l}\geq 0$, again we use the result in (99) to expand this as follows:

[TABLE]

the second condition in the first case of Eq. (86).

Case 2 for $U_{l}^{-}=b$:

[TABLE]

The first condition in the second case of Eq. (86) follows from using these conditions with the event $U_{l}^{-}=b$:

[TABLE]

The condition in (101) implies that $U_{l}-\max(2(U_{l}+V_{l})-N_{l})\geq 0$, using the above we see that this condition is always satisfied in this case:

[TABLE]

Because of this, this condition is omitted from the formulation in Eq. (86). Finally, condition (100) implies:

[TABLE]

the second condition of the case.

Case 3 for $U_{l}^{-}=b$:

[TABLE]

First, we have the event $U_{l}^{-}=b$ taking form:

[TABLE]

which leads us to the first condition in case 3 of Equation (86). For the second condition, start with (103), then we have:

[TABLE]

which is the form of the second condition of Case 4 in (86). Finally, from Condition (104) we have:

[TABLE]

Clearly, (106) is the last condition in the fourth case of Equation (86).

Case 4 for $U_{l}^{-}=b$:

[TABLE]

First, the event $U_{l}^{-}=b$ takes the following form in this case:

[TABLE]

which is the first condition in the fourth case of Eq. (86). Second, we have:

[TABLE]

but since $U_{l}$ is a binomial random variable its value can never be less than 0. Because of this the second condition in the case becomes:

[TABLE]

Finally we have:

[TABLE]

the third condition in the fourth case of (86).

As we did before, we can now expand the definition of $\eta_{l}^{+}$ to derive conditions on $U_{l}$ and $V_{l}$ that lead to the realization of the event $\eta_{l}^{+}=c$. Starting with the definition of $\eta_{l}^{+}$ we have:

[TABLE]

Since there are two nested $\max$ functions in the definition of $\eta_{l}^{+}$, there will be 4 cases that correspond to how the max functions are resolved; each one of these cases is going to represent a different value of $\eta_{l}^{+}$. Below we study each case separately and show that they lead to the four cases in the indicator function of Eq. (87).

Case 1 for $\eta_{l}^{+}=c$:

[TABLE]

First, use both condition to derive a form for the event $\eta_{l}^{+}=c$:

[TABLE]

which is the first condition in the first case of Equation (87). Second, we use the result above to rewrite the first condition:

[TABLE]

which is the second condition in the first case in Equation (87).

Case 2 for $\eta_{l}^{+}=c$:

[TABLE]

First, use the second condition to find a form for the event $\eta_{l}^{+}=c$

[TABLE]

Now use the second condition together with the above to derive:

[TABLE]

This, and Equation (119) are the conditions in the second case of Equation (87).

Case 3 for $\eta_{l}^{+}=c$:

[TABLE]

The first condition of Case 3 in Eq. (87) is given by:

[TABLE]

Now use (121) to rewrite the (120):

[TABLE]

this is the second condition of the third case in Eq. (87). The final condition in the case is given by:

[TABLE]

Case 4 for $\eta_{l}^{+}=c$:

[TABLE]

The first condition in the fourth case of Eq. (87) is obtained as follows:

[TABLE]

The second condition in the same case can be obtained by starting from (123):

[TABLE]

Since $U_{l}$ is a binomial random variable with number of trials $N_{l}^{1}$, it can never be greater then $N_{l}^{1}$, so the condition above becomes:

[TABLE]

which is the second condition in the fourth case of Eq. (87). Now rearrange the terms in the second condition and use the result above to obtain the final condition in the case:

[TABLE]

This concludes the derivation of Eq. (87).

As before, we derive the conditions on $U_{l}$ and $V_{l}$ that lead to the realization of the event $\eta_{l}^{-}=d$ by expanding its definition:

[TABLE]

Using this expansion, we see that there are two possible definitions for the event $\eta_{l}^{-}=d$ in terms of $U_{l}$ and $V_{l}$, both depending on how the $\max$ function is evaluated.

Case 1 for $\eta_{l}^{-}=d$:

[TABLE]

First, we use the above to rewrite the event $\eta_{l}^{-}=d$ for this case:

[TABLE]

Second, we use the result just derived to rewrite the condition for the max in (127) in this case:

[TABLE]

This and (128) are, respectively, the first and second conditions on the first case in Equation (88).

Case 2 for $\eta_{l}^{-}=d$:

[TABLE]

First, we use the above to rewrite the event $d=\eta_{l}^{-}$ for this case:

[TABLE]

Second, we rewrite the condition on the max with the above result:

[TABLE]

These are the condition in the second case of Equation (88).

Finally we have the event $M_{l}=m$: conditions on $U_{l}$ and $V_{l}$ that lead to its realization can again be established by expanding its definition, this will lead to the form of $M_{l}$ in Equation (89). Starting with the definition of $M_{l}$:

[TABLE]

Because of this we see that the definition of $M_{l}$ is composed of two cases that depend on how the $\min$ function is resolved.

Case 1 for $M_{l}=m$:

[TABLE]

First, use the expanded definition of $M_{l}$ with the case above to rewrite the event $M_{l}=m$:

[TABLE]

Second, we use the result just introduced to rewrite the condition on the min for this case:

[TABLE]

leading us to both conditions in the first case of Equation (89).

Case 2 for $M_{l}=m$:

[TABLE]

Again, we use the case above to rewrite the definition of the event $M_{l}=m$:

[TABLE]

Second, we use the above to rewrite the condition for this case:

[TABLE]

These are the two conditions in the definition of $M_{l}=m$ in equation $\eqref{Eq:MSharp}$.

Finally, the forms in equation (84) are simply obtained by conditioning on the event $U_{l}=j,V_{l}=k$ for any of the definitions above. This concludes the proof of the lemma. ∎

Appendix E: .3 Theorem 7

The joint null distribution of $(\chi^{+},\chi^{-})$ is given in the following theorem:

Theorem 7**.**

(Randomization Distribution of $(\chi^{+},\chi^{-})$). For $l=1,\dots,L$ let $(N_{l}^{1},N_{l}^{0},e_{l})$ be fixed and known and let data, $\mathbb{D}=(U_{l},V_{l},\eta_{l},\nu_{l},N_{l}^{t},N_{l}^{c})$, be drawn from the data generating process of Equations (41) – (46). Let $\chi^{+}$ be the maximum of Formulation 2 on $\mathbb{D}$, and let $\chi^{-}$ be the minimum of the problem on the same variables. Let $N^{m}=\sum_{l=1}^{L}\min(N_{l}^{1},N_{l}^{0})$ and define: $\mathcal{X}_{N^{m}}=\left\{\frac{b-c-1}{\sqrt{b+c+1}}:b,c\in\{0,\dots,N^{m}\}\right\}$, Additionally define:

[TABLE]

Let $\mathcal{H}(x,y,r,s)=\mathcal{A}(y)\times\mathcal{B}(y,s)\times\mathcal{C}(x)\times\mathcal{D}(x,r)$ be the Cartesian product of the sets above, such that each element of $\mathcal{H}(x,y,r,s)$ is a 4-tuple of vectors: $(\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d})$. Let $N^{1}=\sum_{l=1}^{L}N_{l}^{1}$; the pmf of $(\chi^{+},\chi^{-})$, for two values $(s,r)\in\mathcal{X}_{N^{m}}\times\mathcal{X}_{N^{m}}$ is:

[TABLE]

where:

[TABLE]

Note that the exact form of the probabilities in $h_{1},h_{2},h_{3}$ is given by Lemma 1. The distribution is essentially a product of many binomial distributions truncated to be defined only on portions of their domains that include possible ranges allowed by the constraints on $\mathcal{A}_{\textrm{good}}$.

Proof.

Proof.All the probability statements throughout the proof are made conditionally on $X$ and $\mathbb{H}_{0}^{sharp}$, for this reason we omit the conditional notation from these statements. Recall also that all the quantities representing counts of units in each stratum are nonnegative integers.

Note first that, for any stratum $l$:

[TABLE]

and

[TABLE]

This implies that:

[TABLE]

where recall that $N^{1}=\sum_{l=1}^{L}N_{l}^{1}$, as defined in the statement of the theorem. By the same argument we can see that $-N^{1}\leq TE^{-}\leq N^{1}$. This fact is useful to bound the domain of $TE^{+}$ and $TE^{-}$: while it is likely that not all integers between $-N^{1}$ and $N^{1}$ have positive probability for $TE^{+}$ and $TE^{-}$ under the distribution of $(\chi^{+},\chi^{-})$, we know that no integers outside the range above will have positive mass under that distribution.

To bound the domain of $\chi=\frac{B-C-1}{\sqrt{B+C+1}}$ note that, by definition, $B_{l}$ represents the number of matched pairs in stratum $l$ that have outcome 1 for the treated unit and 0 for the control unit, independently of how matches are made. Because of this we know that $B_{l}$ can never be greater than either the number of treated units with outcome 1 or control units with outcome 0 in stratum $l$, and, therefore: $B_{l}\leq\min(N_{l}^{1},N_{l}^{0})$, independently of how many units are assigned to treatment. For the same reason we conclude that $C_{l}\leq\min(N_{l}^{1},N_{l}^{0})$ in every stratum, independently of how many units are treated or how the matches are made. Recall that $N^{m}=\sum_{l=1}^{L}\min(N_{l}^{1},N_{l}^{0})$, using these facts we can see that:

[TABLE]

Note finally that, since $\chi\in\mathcal{X}_{N^{m}}$ regardless of how matches are made, we know that $\chi^{+}\in\mathcal{X}_{N^{m}}$ and $\chi^{-}\in\mathcal{X}_{N^{m}}$ because $\chi^{+}$ and $\chi^{-}$ are special cases of $\chi$ in which matches are made with Algorithms 4 and 5 respectively. Because of this we conclude that $(\chi^{+},\chi^{-})\in\mathcal{X}_{N^{m}}\times\mathcal{X}_{N^{m}}$. This set is fast to enumerate computationally and summations over its elements can be performed efficiently. While the distribution of $(\chi^{+},\chi^{-})$ likely does not place positive probability over all values in $\mathcal{X}_{N^{m}}\times\mathcal{X}_{N^{m}}$, we shall consider values in this set and derive conditions under which they have positive probability as well as what their probability is under this distribution.

Begin now by writing the pmf of the range $(\chi^{+},\chi^{-})$. As stated in Theorem 7, assume that $N_{l}^{1},N_{l}^{0}$ are fixed and known for all strata $l=1,\dots,L$. For any two values $(r,s)\in\mathcal{X}_{N^{m}}\times\mathcal{X}_{N^{m}}$ we have:

[TABLE]

The equality in (150) follows by the representation of $SD^{+}$ and $SD^{-}$ in (58) and (59) of Claim 4 respectively. Now we work with each of the four parts separately.

Starting with $f_{3}$, we now show that it must always be that $f_{3}=0$. Consider the event set $\{TE^{-}\geq 1,TE^{+}<1\}$. By definition of $\chi^{+}$, we know that, in this case:

[TABLE]

Therefore, the event $\{TE^{-}\geq 1,TE^{+}<1\}$ implies the event $\{\chi^{+}<\chi^{-}\}$, but this is a contradiction, since $\chi^{+}$ and $\chi^{-}$ are, respectively, the maximum and minimum of the same optimization problem (Formulation 2). Because of this it will never be the case that $\chi^{+}<\chi^{-}$. It must be, then, that the probability of this event under the distribution of interest is 0, and that, therefore, the probability of the event $\{TE^{-}\geq 1,TE^{+}<1\}$ is also 0. Since $f_{3}$ is the probability of the event set $\{TE^{-}\geq 1,\frac{TE^{-}-1}{\sqrt{R^{-}+1}}=s,TE^{+}<1,\frac{TE^{+}-1}{\sqrt{R^{+}+1}}=r\}$, which is included in the 0-probability set above, then it must be that $f_{3}=0$.

This leaves us with $f_{1},f_{2},f_{4}$ to write out. For convenience we repeat the definitions of the sets introduced in the statement of the theorem:

[TABLE]

In addition, let $\mathcal{H}(x,y,r,s)=\mathcal{A}(y)\times\mathcal{B}(y,s)\times\mathcal{C}(x)\times\mathcal{D}(x,r)$ be the Cartesian product of the sets above, such that each element of $\mathcal{H}(x,y,r,s)$ is a 4-tuple of vectors each of length $L$: $(\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d})$. We are now ready to derive $f_{1}(r,s)$:

[TABLE]

Line (155) follows from the definition of $f_{1}$ given in (150). Equation (157) follows from rearranging the terms in the previous line. Equation (159) follows from the definitions in Claim 4. Equation (160) follows from the definition of $\mathcal{H}(x,y,r,s)$, and equation (161) follows from independence of $TE^{+}_{l},TE^{-}_{l},S^{+}_{l},S^{-}_{l},R^{+}_{l},R_{l}^{-}$ in stratum $l$ from the same quantities in any other stratum. This independence is evident from the definitions of these quantities given in Claim 4. The same exact derivation steps leads us to a similar definition for $f_{2}$:

[TABLE]

In the case of $f_{4}$ we can follow the same steps to derive:

[TABLE]

Note that $f_{1}$, $f_{2}$ and $f_{4}$ are nonzero on disjoint parts of the sets $\{-N^{1},\dots,N^{1}\}$, as well as the fact that the definition of $\mathcal{A}(y),\mathcal{B}(y,z),\mathcal{C}(x),\mathcal{D}(x,r)$ and $\mathcal{H}(x,y,r,s)$ doesn’t change for any of the functions. For these reasons we can write the pmf of the range as in Eq (141). It remains to show that the three functions $h_{1}$, $h_{2}$, $h_{3}$ have the form stated in the theorem. This can be accomplished by expanding the inner probabilities of $TE_{l}^{+},TE^{-}_{l},SD^{+}_{l},SD^{-}_{l}$ in each stratum. Before doing this, it will be useful to note that, for any stratum $l$:

[TABLE]

Because of this we consider values of $M_{l}$ in the integer range $0,\dots,N_{l}/2$. Beginning with $\Pr(TE_{l}^{-}=a_{l},S_{l}^{-}=b_{l},TE_{l}^{+}=c_{l},R_{l}^{+}=d_{l})$:

[TABLE]

Equation (163) follows from the representations of $TE$ and $SD$ given in Claim 4, Equation (164) from the definition of conditional probability, Equation (166) from the fact that $a$ and $b$ are constants and therefore they are independent from the other quantities in the equation. Equation (167) follows from the definition of conditional probability. Equations (168) follows from the fact that the event $\{M_{l}-|U_{l}^{+}+\eta_{l}^{+}-M_{l}|=d_{l}\}$ is equal to the event $\{M_{l}-U_{l}^{+}-\eta_{l}^{+}+M_{l}=d_{l}\}\cup\{M_{l}+U_{l}^{+}+\eta_{l}^{+}-M_{l}=d_{l},M_{l}\neq d_{l}\}$, and therefore the probability of its occurrence is equal to the sum of the probability of these two events. Equation (169) follow by rearranging the terms in the previous line, and Equation (170) from summing over values of $\eta_{l}^{-}$. The final line of the derivation is from the definition of $h_{1}$ given in the statement of the theorem. The following derivations for $h_{2}$ and $h_{3}$ follow exactly the same steps, starting with $h_{2}$:

[TABLE]

Finally, $h_{3}$ can be derived from:

[TABLE]

This concludes the derivation of all the quantities in Theorem 7. ∎

Appendix F: Distribution of $(\chi^{+},\chi^{-})$ Under Exclusively Binning Constraints

To compute the distribution we seek, we define the quantities $p_{l}^{t}=\Pr(Y=1|T=1,X=x_{i},i\in S_{l})$ and $\Pr(Y=1|T=0,X=x_{j},j\in S_{l})$. We can encode $\mathbb{H}_{0}$ in this distribution by requiring: $p_{l}^{t}=p_{l}^{c}=p_{l}$ for all $l$. In order for the distribution of interest to correctly be estimated, Assumption 4 must hold. Note that this assumption implies directly that $\Pr(Y(t)=1|X=x_{i},i\in S_{l})=\Pr(Y=1|T=t,X=x_{i},i\in S_{l})$. With it, and by requiring $p_{l}^{t}=p_{l}^{c}$ we obtain $\mathbb{H}_{0}$: $\Pr(Y=1|T=1,X=x_{i},i\in S_{l})=\Pr(Y(1)=1|X=x_{i},i\in S_{l})=\mathbb{E}(Y(1)|X=x_{i},i\in S_{l})=\mathbb{E}(Y(0)|X=x_{i},i\in S_{l})=\Pr(Y(0)=1|X=x_{i},i\in S_{l})=\Pr(Y=1|T=0,X=x_{i},i\in S_{l})$. As stated before, while Assumption 4 is necessary to obtain the distribution of $(\chi^{+},\chi^{-})$.

An estimate of $p_{l}^{t}$ and $p_{l}^{c}$ can be produced in an unbiased and consistent way with $\hat{p}_{l}^{t}=\frac{\sum_{i\in S_{l}}{y_{i}^{t}}}{N_{l}^{t}}$, or any other estimator of choice. When one wants to test the null hypothesis that $p_{l}^{t}=p_{l}^{c}$, one can use $\hat{p}_{l}^{t}=\hat{p}_{l}^{c}=\frac{\sum_{i\in S_{l}}y_{i}^{t}+\sum_{j\in S_{l}}y_{j}^{c}}{N_{l}}$.

Recall that exactly $M_{l}=\min(N_{l}^{t},N_{l}^{c})$ matches are made, with $M=\sum_{l=1}^{L}M_{l}$. Now denote with $N_{l}^{t}$ the number of treated units in stratum $S_{l}$ and with $N_{l}^{c}$ the number of control units in that same stratum. For each stratum, $l=1,\dots,L$, the data are generated as follows:

[TABLE]

The joint distribution of the truncated variables introduced in Equations (48)-(51) under this DGP is given in the following lemma:

Lemma 2**.**

(Distribution of Truncated Variables) For all $l=1,\dots,L$, let $(N_{l}^{t},N_{l}^{c},p_{l}^{t},p_{l}^{c})$ be fixed and known, and let $G_{l}^{+}=\max(N_{l}^{c}-N_{l}^{t},0)$, $G_{l}^{-}=\max(N_{l}^{t}-N_{l}^{c},0)$. Let data $\mathbb{D}=\{(U_{l},V_{l},\eta_{l},\nu_{l})\}_{l=1}^{L}$ be drawn according to Equations (172)-(175). The truncated variables have the following joint distributions:

[TABLE]

where $\Pr(U_{l}=x)=Bin(x,p_{l}^{t},N_{l}^{t})$, $\Pr(x\leq U_{l}\leq y)=\sum_{z=x}^{y}Bin(z,p_{l}^{t},N_{l}^{t})$, $\Pr(\eta_{l}=x)=Bin(x,p_{l}^{c},N_{l}^{c})$, and $\Pr(x\leq\eta_{l}\leq y)=\sum_{z=x}^{y}Bin(z,p_{l}^{c},N_{l}^{c})$.

Proof.

Proof We will prove the result for $\Pr(U_{l}^{-}=a,U_{l}^{+}=b)$ only, as the proof for $\Pr(\eta_{l}^{-}=a,\eta_{l}^{+}=b)$ is symmetrical. Recall first that we are now in a case in which the number of treated and control units in each stratum, $N_{l}^{t}$ and $N_{l}^{c}$ respectively, are fixed. We use this together with Assumption 4 to write the marginal distributions of $U_{l}^{+}$ and $U_{l}^{-}$:

[TABLE]

In the above, $G_{l}^{-}=\max(N_{l}^{t}-N_{l}^{c},0)$, a fixed quantity. The first equality in both formulas follows from the definition of $U_{l}^{+}$ and $U_{l}^{-}$ given in Equations (48) and (49) respectively. The two cases of each definition follow directly from these definitions. Given the marginal distributions above, we can easily find distributions for $U_{l}^{+}$ and $U_{l}^{-}$ conditional on $U_{l}=x$:

[TABLE]

This is because, if the number of treated and control units in each stratum is fixed, randomness in $U_{l}^{+}$ and $U_{l}^{-}$ comes only from $U_{l}$, as is evident from the definitions of $U_{l}^{+}$ and $U_{l}^{-}$ given in Eqns. (85) and (86). Therefore once that is fixed the two quantities become constant and, as such, independent of each other. Given this fact, it is evident that $\Pr(U_{l}^{-}=a,U_{l}^{+}=b|U_{l}=x)=\mathbb{I}(U_{l}^{-}=a,U_{l}^{+}=b|U_{l}=x)=\mathbb{I}(U_{l}^{-}=a|U_{l}=x)\mathbb{I}(U_{l}^{+}=b|U_{l}=x)$. This product of indicator functions can also be written as four conditions, each one representing the intersection of the event sets on which both distributions place nonzero probability:

[TABLE]

This is evident from Equations (178) and (179). The event sets in the equation can be simplified further as some of their elements are redundant, we show this in the following.

Case 1: $b<N_{l}^{c},x=b$ and $a=0,x\leq G_{l}^{-}$.

First, using the second and fourth condition we can write $b\leq G_{l}^{-}$. Then, since $b$ has to be less than both $G_{l}^{-}$ and $N_{l}^{c}$ we can rewrite the condition as: $b<\min(N_{l}^{c},G_{l}^{-}+1)$, where we replace the equality on $G_{l}^{-}$ with $b$ being strictly less than $G_{l}^{-}+1$. This is possible because $G_{l}^{+}$ is, by definition, a nonnegative integer. The final event here becomes: $x=b,\;b<\min(N_{l}^{c},G_{l}^{-}+1),\;a=0$.

Case 2: $b<N_{l}^{c},x=b$ and $a>0,x=a+G_{l}^{-}$.

First, using the first, second and fourth condition we have: $N_{l}^{c}>b=x=a+G_{l}^{-}$, which we rewrite as the event: $x=b,b=a+G_{l}^{-},a<N_{l}^{c}-G_{l}^{-}$. With this the final event becomes: $x=b,\;b=a+G_{l}^{-},\;0<a<N_{l}^{c}-G_{l}^{-}$.

Case 3: $b=N_{l}^{c},x\geq N_{l}^{c}$ and $a=0,x\leq G_{l}^{-}$.

Here we can combine the second and fourth conditions to obtain: $N_{l}^{c}\leq x\leq G_{l}^{-}$. This leads to the event: $b=N_{l}^{c},\;N_{l}^{c}\leq x\leq G_{l}^{-},\;a=0$.

Case 4: $b=N_{l}^{c},x\geq N_{l}^{c}$ and $a>0,x=a+G_{l}^{-}$.

First from the second and fourth events we have: $x=a+G_{l}^{-}\geq N_{l}^{c}$, which can be used to obtain the event: $a\geq N_{l}^{c}-G_{l}^{-}$. We also rewrite the third condition as $a\geq 1$, leading to the following final representation for the event of this case: $b=N_{l}^{c},\;a\geq\max(N_{l}^{c}-G_{l}^{-},1),\;x=a+G_{l}^{-}$.

Finally, we can put all these events together to obtain a simplified formulation for the joint conditional distribution of $U_{l}^{+},U_{l}^{-}$:

[TABLE]

The marginal distribution of $U_{l}^{+},U_{l}^{-}$ can now be found by using the law of total probability to sum over $U_{l}^{+}$:

[TABLE]

To obtain the final line the condition on $x$ in the event is substituted in the probability. This concludes the derivation of the pmf of $U_{l}^{+},U_{l}^{-}$ and the proof of this lemma.

∎

Using Lemma 2, these distributions are easy to enumerate and to construct a lookup table for. They are also clearly symmetrical and their form implies that $\Pr(U_{l}^{-}=a,U_{l}^{+}=b,\eta_{l}^{-}=c,\eta_{l}^{+}=d)=\Pr(U_{l}^{-}=a,U_{l}^{+}=b)\Pr(\eta_{l}^{-}=c,\eta_{l}^{+}=d)$. This fact is going to allow us to derive an expression for the distribution of $(\chi^{+},\chi^{-})$, given in the following theorem.

Theorem 8**.**

(Distribution of $(\chi^{+},\chi^{-})$) For all $l=1,\dots,L$, let $(N_{l}^{t},N_{l}^{c},p_{l}^{t},p_{l}^{c})$ be fixed and known, and let $M_{l}=\min(N_{l}^{t},N_{l}^{c})$, $M=\sum_{l=1}^{L}M_{l}$. Let data $\mathbb{D}=\{(U_{l},V_{l},\eta_{l},\nu_{l})\}_{l=1}^{L}$ be drawn according to Equations (172)-(175). Let $\chi^{+}$ be the maximum of Formulation 2 on $\mathbb{D}$ and let $\chi^{-}$ be the minimum of Formulation 2 also on $\mathbb{D}$. Let $\mathcal{X}_{M}:=\left\{\frac{b-c-1}{\sqrt{b+c+1}}:b,c\in\{0,\dots,M\}\right\}$. Additionally, let $\mathcal{A}(y),\mathcal{B}(y,s),\mathcal{C}(x),\mathcal{D}(x,r),\mathcal{H}(x,y,r,s)$ be defined as in Theorem 7. The pmf of $(\chi^{+},\chi^{-})$, for two values $s,r\in\mathcal{X}_{M}$ is given by:

[TABLE]

where:

[TABLE]

The probabilities $\Pr(U_{l}^{-}=a,U_{l}^{+}=b)$ and $\Pr(\eta_{l}^{-}=a,U_{l}^{+}=b)$ are given in Lemma 2 and depend on $N_{l}^{t},N_{l}^{c},p_{l}^{t},p_{l}^{c}$.

Proof.

Proof The proof follows closely the template of the proof of Theorem 7. Note first that in this case $N_{l}^{t}$ and $N_{l}^{c}$ are fixed and known for all strata $l$. We first bound the values of $U_{l}^{+}$ and $\eta_{l}^{+}$:

[TABLE]

From this it follows that $-M_{l}\leq TE_{l}^{+}=U_{l}^{+}-\eta_{l}^{+}\leq M_{l}$. Recall that $M=\sum_{l=1}^{M}M_{l}$ and that $TE^{+}=\sum_{l=1}^{L}TE^{+}_{l}$: it must be that $-M\leq TE^{+}\leq M$. An exactly symmetric argument shows that $-M\leq TE^{-}\leq M$. Because of the above the values that we will consider for $TE^{+}$ and $TE^{-}$ in this proof are the integers ranging between $-M$ and $M$. Finally, as detailed in the proof of Theorem 7, it must be that $B$ and $C$ are always less than $M$, the total number of matches made, as they represent counts of matched pairs. This allows us to bound the set of values that have positive probability under the distribution of $(\chi^{+},\chi^{-})$:

[TABLE]

Given the set above, we will consider values of $(\chi^{+},\chi^{-})$ in the set $\mathcal{X}_{M}$, even though not all of the set’s elements will have positive probability under the distribution of interest. As done in Theorem 7 we derive a set of conditions on values within $\mathcal{X}_{M}$ that determine whether or not those values have positive probability and, if yes, what is their probability under the distribution of $(\chi^{+},\chi^{-})$ under Assumption 4.

The steps leading to the form of (180) in the main statement of the theorem are exactly the same as those in the proof of Theorem 7, with the sole difference that the possible values of $TE^{+}$ and $TE^{-}$ now are between $-M$ and $M$. Because of this, the summations over $x$ and $y$ both in the definition of the joint probability of interest and in the functions, $f_{1},f_{2},f_{3},f_{4}$, introduced in the proof of Theorem 7, range from $-M$ to $M$ in this case. In particular, the summation over $x$ and $y$ in $f_{1}$ ranges from $-M$ to 0 in this case, the summation over $x$ in $f_{2}$ ranges from $1$ to $M$, the summation over $y$ in $f_{2}$ now ranges from $-M$ to 0, and the summations over $x$ and $y$ in $f_{4}$ range from 1 tom $M$ in this case. The definitions of $f_{1},f_{2},f_{4}$ are also exactly the same as those in the proof of Theorem 7, including the fact that $f_{3}=0$. The argument made for this in the proof of Theorem 7 applies exactly in the same way for this case.

The only other difference is in the inner probabilities for the joint distribution of the statistics in each stratum. These inner probabilities result in equations (181)-(183), and can be derived as follows:

[TABLE]

In the above, Equation (184) follows from the definitions of $TE_{l}^{-},S_{l}^{-},TE_{l}^{+},R_{l}^{+}$ given in Claim 4. Equation (185) follows from rearranging the terms in the last equality in the parentheses, and from plugging in the equality for $U_{l}^{-}-\eta_{l}^{-}$, (186) follows from the fact that $a_{l}$ and $b_{l}$ are both constants and so the event $|a_{l}|=b_{l}$ is independent from all the others inside the parentheses. Equation 187 follows from the fact that the event $\{|U_{l}^{+}+\eta_{l}^{+}-M_{l}|=M_{l}-d_{l}\}=\{U_{l}^{+}+\eta_{l}^{+}-M_{l}=M_{l}-d_{l}\}\cup\{U_{l}^{+}+\eta_{l}^{+}-M_{l}=d_{l}-M_{l},M_{l}\neq d_{l}\}$, This follows from the definition of the absolute value function applied to positive integers. Note that the two sets are disjoint, so the probability of their events can be added and that $M_{l}$ and $d_{l}$ are fixed quantities, so they are independent of all the other events in the statement. Equation (188) follows from rearranging the terms in the equalities for $\eta_{l}^{+}$. Equation (189) follows from summing over all values of $\eta_{l}^{-}$ and from plugging in equalities for $\eta_{l}^{+}$ into the equation for $U_{l}^{+}$. Finally, (190) follows from the definition of $g_{1}(a_{l},b_{l},c_{l},d_{l})$ given in the statement of the theorem. The derivations for $g_{2}$ and $g_{3}$ follow exactly the same steps. For $g_{2}$ we have:

[TABLE]

Finally, $g_{3}$ can be derived with the same steps as $g_{1}$:

[TABLE]

This concludes the proof of Theorem 8. ∎

Appendix G: Proof of Theorem 5

Recall that, by Theorems 8 and 7, the distributions of matched pairs in each stratum are independent of each other because of SUTVA, and because optimization of $\chi$ is performed almost independently in each stratum. Since the calculation of the full joint distribution of McNemar’s statistic is a convolution of independent discrete probabilities, it can be sped up with multidimensional Fast Fourier Transforms. The basic row-column mDFT algorithm reduces the complexity of having to generate $\mathcal{H}(x,y,r,z)$ to $O(4mN^{4m}\log N)$ [see, e.g., Dudgeon, 1983]. Here $N$ is the total number of units, $m$ is the largest number of matches made in any of the strata. This last quantity is repeated four times, as there are four variables for which probabilities must be computed. These are the densities being convolved together in the summation. This, coupled with simple enumeration to create lookup tables for distributions of $g_{1}$, $g_{2}$, $g_{3}$ in each stratum (See Theorems 7 and 8 for the definition of these quantities), as well as the final range distribution yields an algorithm that creates a probability table for the null distribution of $(\chi^{+},\chi^{-})$ in $O(4mN^{4m}\log N)$. This implies that the worst-case time of computing the distribution in Theorem 7 is polynomial.

Appendix H: Plots of exact distributions of $(\chi^{+},\chi^{-})$ on a simulated dataset

Figure 10 shows marginal distributions of $\chi^{+}$ and $\chi^{-}$ derived from the joint distributions in Theorems 7 and 8, respectively. We simulated datasets from the DGPs described in those sections and estimated the distributions with the formulas in Theorems 3 and 4, and the polynomial time algorithm suggested in Theorem 5. Panel (a) shows a location and scale difference in the distribution of $(\chi^{+},\chi^{-})$ between when the null is true and when it is not. This demonstrates our tests’ capacity to detect full-sample treatment effects.

Appendix I: Additional Simulations

Appendix I: .1 Performance of robust tests with increasing number of covariates

We present results of simulation studies under the same settings as those in Section 6, but we instead keep the number of units fixed at 200, and let the number of covariates grow.

Results for McNemar’s test are reported in Figure 11, and results for the z-test are reported in Figure 12. We see that, largely, our tests tend to perform well even in the increasing $P$ regime: both tests achieve the overall lowest error rate when the ATT is 0 (closest to the error rate of the idealized test that includes both true potential outcomes), and display good statistical power as the ATT increases in relative strength. The z-test is somewhat conservative when the true ATT is weak (Cohen’s d=0.2), and when there are many covariates: this is a desirable pattern since the larger number of covariates makes the matching problem harder. Notably, the benchmark test performed on the true potential outcomes displays similar behavior to our robust z-test in Figure 12 when $d=0.2$.

Appendix I: .2 Performance of robust tests with Exclusively Binning Constraints

We present performance results of our robust tests when constraints are exclusively binning, i.e., when data is perfectly stratified. In this setting we partition the $N$ units into strata such that each stratum contains exactly $25$ units, for example, when $N=100$ there will be 4 strata. We denote the number of strata for a given $N$ with $L_{N}$, and use the variable $s_{i}=1,\dots,L_{N}$ to denote which stratum unit $i$ is assigned to. For a set of fixed coefficients $\beta_{1},\dots,\beta_{L_{N}}$ that are calibrated as described prior, we generate:

[TABLE]

with treatments $t_{i}$ assigned to exactly $N^{t}$ units but each with probability $e_{i}$. For binary data, we set $y_{i}=\mathbb{I}[y_{i}^{*}>0]$, and for continuous data $y_{i}=y_{i}^{*}$, as done in all the other simulations. All other settings are as in the simulations presented in Section 6 of the main paper.

Results for the robust McNemar’s test are presented in Figure 13.

Results for the robust z-test are presented in Figure 14

We see that both test perform well under these types of constraints: error rate for the z-test is almost 0 when the ATT is 0, and the test achieves good statistical power as the ATT grows in strength. A similar patterns is visible for McNemar’s test, though this test requires a slightly larger sample size ($N\geq 100$) in order to achieve a good false rejection rate.

Appendix I: .3 Performance of robust z-test with different types of constraints

We study whether the type of match quality constraints used in our robust test affects both the error rate of the robust tests as well as their runtime. All simulation settings are kept same as before, except that our robust test is ran 3 times per iteration, each time with a different type of constraint. The types of constraints we compare are as follows:

• •

Moment: We require balance between the first two moments of the empirical distributions of each covariate in the treated/control sets after matching. This is formulated as: $|\frac{1}{M}\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}a_{ij}(x^{t}_{ip}-x^{c}_{jp})|\leq\epsilon_{p}^{mean}\times(\sigma(x_{p}^{t})/2+\sigma(x_{p}^{c})/2)$, for $p=1,\dots,P$, and $|\frac{1}{M}\sum_{i=1}^{N^{t}}\sum_{j=1}^{N^{c}}a_{ij}(x^{t}_{ip}-x^{c}_{jp})^{2}|\leq\epsilon_{p}^{var}\times(\sigma^{2}(x_{p}^{t})/2+\sigma^{2}(x_{p}^{c})/2)$, for $p=1,\dots,P$ ($2P$ constraints total).

• •

Caliper: We only allow units to be matched if their propensity scores are close enough. This is formulated as: $a_{ij}|e^{t}_{i}-e^{c}_{j}|\leq\epsilon^{cal}$ for $i=1,\dots N^{t}$, $j=1,\dots,N^{c}$ ($N^{t}\times N^{c}$ constraints total).

• •

Quantile: We require the empirical distributions of each covariate in the treated/control sets after matching to be similar in terms of their quantiles. For each covariate $p=1,\dots,P$, and a pre-specified sequence of quantiles $\gamma_{p1},\dots,\gamma_{pL}$, such that $\gamma_{p1}\leq\gamma_{p2}\leq\dots\gamma_{pL}$, we require: $|\frac{1}{M}\sum_{i=1}^{N^{t}}a_{ij}\mathbb{I}[x^{t}_{ip}\leq\gamma_{p}]-\sum_{j=1}^{N^{c}}\sum_{i=1}^{N^{t}}a_{ij}\mathbb{I}[x^{c}_{jp}\leq\gamma_{p}]|\leq\epsilon^{qua}_{p}$ ($P\times L$ constraints total).

as previously done, tolerance values ($\epsilon^{mom}_{p}$, $\epsilon^{cal}$, $\epsilon_{p}^{qua}$) are chosen as the smallest feasible values on an evenly spaced grid between 0.1 and 1, and $\sigma_{p}^{2}(x_{p}^{t})$, $\sigma_{p}(x_{p}^{t})$, $\sigma_{p}^{2}(x_{p}^{c})$, $\sigma_{p}(x_{p}^{c})$ are the variance and standard deviation of covariate $x_{p}$ in the treated and control set respectively.

Results are reported in Figure 15. In both fixed $P$ and fixed $N$ regimes we see that the pscore caliper constraint seems to display the lowest rejection rate when the true ATT is 0, indicating that this type of constraint is particularly suitable for the settings we studied. Interestingly, the caliper constraint seems to also be the one able to achieve better statistical power as the true ATT increases in size in the fixed $P$ regime, while it displays worse statistical power than the other two constraints in the fixed-$N$ regime. In terms of runtime, there appears to not be virtually any difference between constraint types in the fixed-$P$ regime, however the caliper constraints is much faster and unaffected by the number of covariates in the fixed-$N$ regime, as the number of covariate grows, which is expected due to the propensity score being a 1-dimensional quantity for any amount of covariates. Also expectedly, the quantile constraint requires longer execution time as the number of covariates grows compared to the moment constraint.

The main takeaway from this excercise is that it is a good idea for analyst to include several constraints at once in their application of the robust tests, requiring several different metrics of balance to all be below a certain tolerance to admit a match assignment as good.

Appendix I: .4 Performance of robust z-test with different MIP approximation methods

Since MIP solution algorithms have a runtime that is exponential in the number of decision variables (units to match, in our case), we have implemented and tested two different approximation methods to speed up computation for our robust tests. Both methods are based on first solving a linear relaxation of the original MIP, and then solving an integer program that has a solution closest to that of the relaxation. The two methods are described as follows:

Approximation 1

: This is the same approximation method used in Zubizarreta [2012]. We first solve a linear relaxation of all the MIPs used for the robust z-test, (Formulation 2) and once an optimal solution is found via our proposed algorithm, we approximate a MIP solution by solving the following IP. Let $b_{11},\dots,b_{N^{t}N^{c}}$ be the continuous matching indicators output by solving the linear relaxation of Formulation 2. The approximation can be obtained by solving:

[TABLE]

The approximation tries match as many units entirely as the linear relaxation matches fractionally, while not matching any unit more than once. This approximation provides a substantial speedup over the MIP, since it removes all the additional constraints, however the solution could be sub-optimal and some constraints may be violated in practice. Results from comparisons on simulated data in Figure 16 show that, in practice, this approximation is almost as reliable as solving the exact MIP, making it a viable option in practice.

Approximation 2

: Approximation 2 solves the same problem as Approximation 1, except that all the additional constraints (including user-defined balance constraints) in Formulation 2 are kept. This ensures that the solution found satisfies the constraints specified, but is more expensive in terms of computation time. As Figure 16 shows, this approximation performs even more closely to the MIP: gains in computation time are offset by almost no loss in performance.

Appendix I: .5 Comparison To Matching Bounds [Morucci and Rudin, 2018]

We perform a comparison of our robust tests to Matching Bounds [Morucci and Rudin, 2018], a method for a different facet of the issue of multiple match assignment disagreeing. There are several substantial differences between the Matching Bounds approach of Morucci and Rudin [2018] and the robust testing approach proposed in this work. Namely that Matching Bounds is a tool for finding ranges of treatment effects, our robust tests are a tool for finding ranges of test statistics (i.e., we consider not just treatment effects but the uncertainty in their estimates). An analyst interested in knowing what possible treatment effects could be obtained by changing matched samples should employ Matching Bounds, while an analyst interested in obtaining the range of test statistics (and p-values) that could be obtained by matching in different ways should use the robust tests. Specifically [Morucci and Rudin, 2018] finds bounds on the raw value of the treatment effect of interest, our robust testing approach targets a test statistic for the treatment effect: this latter quantity depends both on TE value and on TE variance. This latter dependence implies that the optimum Z-value for a robust test will likely not be the same Z-value obtained by, for example, performing a Z-test via normal approximation on the set of units that maximizes/minimizes the TE as selected by solving the optimization problem in [Morucci and Rudin, 2018]. Specifically, a Z-test conducted on matches selected as in [Morucci and Rudin, 2018] will find a sub-optimal solution to the Z-statistic optimization problem, potentially leading to increased chances of type-I error (incorrect rejection of the null.).

To further demonstrate the differences between the two methods and their purposes, we have performed a set of simulations whose settings are similar to those in Section 6 of our paper and results are reported in Figure LABEL:fig:mbcomparison. One can see the benefit of finding the range of test statistics, as compared to just the range of treatment effects: P-values output by MBs are not as conservative as those from the robust test when the ATT is 0, and the bounds on the ATT output by MBs are better than those output by the robust test when the ATT is 0.

Appendix J: Agreement of traditional matching methods in case studies 1 and 2

As a supplement to our case studies, we present p-values obtained after applying traditional matching algorithms to the datasets of Case Studies 1 and 2. We do so to gauge whether there is disagreement in rejection decisions between matching methods on these datasets.

We have constrained all the matching methods used to respect the caliper constraint we employ in constructing our robust test, e.g., units whose absolute distance in the covariates exceeds the thresholds defined in Section 7 are disallowed from being matched. This ensures that all three methods compared produce at the very least balance on the covariates up to that threshold. In this sense, all these methods produce similar balance.

Results of this comparison are presented in Table 4. We see that, generally methods tend to fail to reject on the GLOW dataset: this is because all methods fail to maximize the effective sample size of McNemar’s test. They match units with the same outcome, and therefore all produce low-powered tests. Instead, we constrain our tests to output only matches that have a good effective sample size (almost all the treated units) – in this case we do indeed demonstrate that an effect is present.

For the Bike Sharing case study, we see that propensity score matching and L2 distance matching fail to reject the null hypothesis, while optimal constrained matching does in fact reject it. The disagreement here is likely due to the fact that several matched sets that meet the balance requirements are present, as we argue in our paper. Our application of the robust z-test to this dataset further corroborates this hypothesis by showing that, indeed, it cannot reject the null hypothesis over all possible good matches.

Appendix K: Additional Information about Case Study 2

In our second case study, we used 2 years (2011-2012) of bike sharing data from the Capital Bike Sharing (CBS) system [see Fanaee-T and Gama, 2014] from Washington DC comprised of 3,807,587 records over 731 days, from which we chose 247 treatment days and 463 control days according to the weather as follows: The control group consists of days with Weather 1: Clear, Few clouds, Partly cloudy, and Partly cloudy. The treatment group consists of days with Weather 2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist. The covariates for matching are as follows: Season (Spring; Summer; Fall; Winter), Year (2011; 2012), Workday (No; Yes), Temperature (maximum 41 degree Celsius), Humidity (maximum 100 percent), Wind speed (maximum 67). The outcome is the total number of rental bikes. We computed distance between days as follows: $\textrm{dist}_{ij}$ =1 if covariates season, year and workday were the same, and the differences in temperature, humidity and wind speed are less or equal to 2, 6 and 6, respectively for treated unit i and control unit j, 0 otherwise. Since the outcome variable is continuous, we focus on testing $\mathbb{H}_{0}^{SATT}$ by producing ranges of p-values with the method introduced in Section LABEL:Sec:approx. Optimization models for this case study have been implemented in AMPL [Fourer et al., 2002], and solved with the solver CPLEX [ILOG, 2007]. The bike sharing and GLOW datasets are publicly available. The reported solution time with a X64-based PC with Intel(R) Core(TM) i7-4790 CPU running at 3.60 GHz with 16 GB memory and mip gap of 0.001 to solve a single instance is less than 1 second for all the tests.

Appendix L: Additonal Case Studies

Appendix L: .1 Case Study 3: Training Program Evaluation

In this case study, we used training program evaluation data described in Dehejia and Wahba [1999], and Dehejia and Wahba [2002], which were drawn from Lalonde [1986]. This data set contains 15,992 control units and 297 treatment units. The covariates for matching are as follows: age, education, Black (1 if black, 0 otherwise), Hispanic (1 if Hispanic, 0 otherwise), married (1 if married, 0 otherwise), nodegree (1 if no degree, 0 otherwise), and earnings in 1975. The outcome is the earnings in 1978. The treatment variable is whether an individual receives job training. We computed distance between units as follows: $\textrm{dist}_{ij}$ =1 if covariates Black, Hispanic, married and nondegree were the same, and the differences in age, education and earnings in 1975 were less or equal to 5, 4 and 4000, respectively, for treated unit i and control unit j, 0 otherwise.

In the Figure 18, the P-value upper bounds are 1 and the P-value lower bounds are 0, illustrating that there is a lot of uncertainty associated with the choice of experimenter – one experimenter choosing 287 matched pairs can find a P-value of $\sim$ 0 and declare a statistically significant difference while another experimenter can find a P-value of $\sim$ 1 and declare the opposite. In this case it is truly unclear whether or not training has an effect on the earnings. To sanity check whether a reasonably sized effect would have been detected had one been present, we injected synthetic random noise (with normal distribution of mean $\simeq$ \$$10,000 and standard deviation \simeq$$$$100) on the treatment outcome, and the$z$-test robustly detects the treatment effect before the solutions become infeasible.

Appendix L: .2 Case Study 4: Crime and Transition into Adulthood

In this case study we have used the data from a U.S. Department of Justice study regarding crime during the transition to adulthood, for youths raised in out-of-home care [see Courtney and Cusick, 2010]. Each observation represents a youth, and the outcome is whether he or she committed a violent crime over the 3 waves of the study.

The (binary) covariates are as follows: hispanic, white, black, other race, alcohol or drug dependency, mental health diagnosis, history of neglect, entry over age of 12, entry age under 12, in school or employed, prior violent crime, any prior delinquency. The “treatment” variable is whether or not the individual is female; in particular we want to determine whether being female (controlling for race, criminal history, school/employment and relationship with parents) influences the probability of committing a violent crime. Here $\textrm{dist}_{ij}$ =1 whenever all covariates of treated unit i are the same as those of the control unit j, 0 otherwise.

Figure 7 is constructed in an analogous way to Figure 19 (using McNemar’s test rather the $z$-test) showing the total number of discordant pairs along the x axis. Here, any matched pairs assignment would show a significant difference for the risks of violence between males and females. This difference becomes more pronounced as the number of pairs increases. Thus, the outcome is robust to the choice of experimenter.

Appendix M: Robust Wilcoxon Signed Rank Test

In this section, we present an extension of our robust framework to the Wilcoxon Signed Rank Test. This extension is intended to show that our framework is general enough to be extendable to other test statistics, as well as provide one more robust test for practitioners.

The Wilcoxon Signed Rank Test is a more powerful test on the median than the sign test, but with a stronger assumption: it assumes that the population distribution is symmetric. We define the following parameters and decision variables to formulate the model:

• •

$T_{i}$  is the the outcome of a treated observation $i$ in the treatment group

• •

$C_{j}$  is the the outcome of a control observation $j$ in the control group

• •

$s_{ij}$  is the $ij$th element of a matrix. It takes value 1 if the covariates of treated observation $i$

  and control observation $j$ are similar enough to be a possible matched pair, otherwise 0

• •

$\delta_{ij}$  is a binary precomputed parameter, which is equal to 1 if $T_{i}$ is greater than $C_{j}$, 0 otherwise

• •

$g_{ijkl}$  is a binary precomputed parameter, which is equal to 1 if $\eta_{ij}$ is greater than $\eta_{kl}$, where $\eta_{ij}=\left|T_{i}-C_{j}\right|$, 0 otherwise

There are three matrices of decision variables, namely:

• •

$a_{ij}$  is a binary variable that is 1 if i and j are in the same pair, otherwise 0

• •

$z_{ijkl}$  is a binary variable that is 1 if $g_{ijkl}$, $a_{ij}$ and $a_{kl}$ all are equal to 1, 0 otherwise. Intuitively $z_{ijkl}$ is 1 only when $ij$ and $kl$ are both being used as pairs and when $ij$’s absolute difference is larger than $kl$’s absolute difference.

• •

$h_{ij}$  is an integer variable whose value is the rank of pair $ij$. It is the count of $kl$ pairs ranked beneath pair $ij$ (plus one, so that the lowest rank is 1 rather than 0).

To compute the test statistic in the traditional way, one would compute absolute differences $\eta_{ij}$ and rank order them. The test statistic is the sum of ranks of the positive differences.

The formulation is as follows :

[TABLE]

subject to:

[TABLE]

Equations (191) to (194) are used to ensure that $z_{ijkl}$=1 when $g_{ijkl}$, $a_{ij}$ and $a_{kl}$ all are one. Equation (195) is used to maintain covariates constraint such that only when the precomputed parameter $s_{ij}$ is 1 then $a_{ij}$ is allowed to take a value of 1. Constraints (196)-(198) are the same constraints as in the sign test, to make sure we have $n$ pairs with one treatment and control observation in each pair. Equation (199) defines binary variables $a_{ij}$. Equation (200) is used to calculate rank $h_{ij}$ for each pair $ij$. If the pair $ij$ is not being used then $a_{ij}$ will be 0, which means $z_{ijkl}$ will be 0 by Constraint (192). This means the only $ij$ pairs that have positive heights are those that are being used as matched pairs. The objective will only add up heights of pairs $ij$ for which the absolute differences are positive, which was precomputed as $\delta_{ij}$.

The z-score is computed from the optimal value of $w_{+}$ through the following formula, and the pvalue is computed as usual.

[TABLE]

Appendix N: Integer Linear Programming Basics

ILP techniques have become practical for many large-scale problems over the past decade, due to a combination of increased processor speeds and better ILP solvers. Any type of logical condition can be encoded as linear constraints in an ILP formulation with binary or integer variables. Consider two binary variables $x\in\{0,1\}$ and $y\in\{0,1\}$. The logical condition “if $y=0$ then $x=0$” can be simply encoded as

[TABLE]

Note that this condition imposes no condition on $x$ when $y=1$. Translating if-then constraints into linear constraints can sometimes be more complicated; suppose, we would like to encode the logical condition that if a function $f(w)$ is greater than 0, then another function $g(w)$ is greater or equal to 0. We can use the following two linear equations to do this, where $\theta$ is a binary variable and $M$ is a positive number that is larger than the maximum values of both $f$ and $g$:

[TABLE]

In order for $f(w)$ to be positive, then $\theta$ must be 0, in which case, $g(w)$ is then restricted to be positive. If $f(w)$ is negative, $\theta$ must be 0, in which case no restriction is placed on the sign of $g(w)$. (See for instance the textbook of Winston and Venkataramanan [2003], for more examples of if-then constraints).

ILP can capture other types of logical statements as well. Suppose we would like to incorporate a restriction such that the integer variable $S_{i}$ takes a value of $K$ only if $i=t$, and 0 otherwise. The following four if-then constraints can be used to express this statement, where $\lambda_{1}$ and $\lambda_{2}$ are binary variables:

[TABLE]

Each of these if-then constraints (4)–(7) can be converted to a set of equivalent linear equations, similar to what we described above. (See also Noor-E-Alam et al. [2012] and Winston and Venkataramanan [2003]).

There is no known polynomial-time algorithm for solving ILP problems as they are generally NP-hard, but they can be solved in practice by a number of well-known techniques (Wolsey [1998]). The LP relaxation of an ILP provides bounds on the optimal solution, where the LP relaxation of an ILP is where the integer constraints are relaxed and the variables are allowed to take non-integer (real) values. For instance, if we are solving a maximization problem, the solution of the LP relaxation can serve as an upper bound, since it solves a problem with a larger feasible region, and thus attains a value at least as high as that of the more restricted integer program. ILP solvers use branch-and-bound or cutting plane algorithms combined with other heuristics, and are useful for cases where the optimal integer optimal solution is not attained by the LP relaxation. The branch-and-bound algorithms often use LP relaxation and semi-relaxed problems as subroutines to obtain upper bounds and lower bounds, in order to determine how to traverse the branch-and-bound search tree [Chen et al., 2011, Wolsey, 1998]. The most popular ILP solvers such as CPLEX, Gurobi and MINTO each have different versions of branch-and-bound techniques with cutting plane algorithms and problem-specific heuristics.