Abstract
Under current policy decision making paradigm, we make or evaluate a policy decision by intervening different socioeconomic parameters and analyzing the impact of those interventions. This process involves identifying the causal relation between interventions and outcomes. Matching method is one of the popular techniques to identify such causal relations. However, in onetoone matching, when a treatment or control unit has multiple pair assignment options with similar match quality, different matching algorithms often assign different pairs. Since, all the matching algorithms assign pair without considering the outcomes, it is possible that with same data and same hypothesis, different experimenters can make different conclusions. This problem becomes more prominent in case of largescale observational studies. Recently, a robust approach is proposed to tackle the uncertainty which uses discrete optimization techniques to explore all possible assignments. Though optimization techniques are very efficient in its own way, they are not scalable to big data. In this work, we consider causal inference testing with binary outcomes and propose computationally efficient algorithms that are scalable to largescale observational studies. By leveraging the structure of the optimization model, we propose a robustness condition which further reduces the computational burden. We validate the efficiency of the proposed algorithms by testing the causal relation between Hospital Readmission Reduction Program (HRRP) and readmission to different hospital (nonindex readmission) on the State of California Patient Discharge Database from 2010 to 2014. Our result shows that HRRP has a causal relation with the increase in nonindex readmission and the proposed algorithms proved to be highly scalable in testing causal relations from largescale observational studies.
Key words: Causal Inference, Policy Evaluation, Matching Method, Robustness, Public policy, Observational Study, Big Data.
Introduction
Effective and evidencebased public policy decisions aim to manipulate one or many socioeconomic variables and analyze their impact on the desired outcomes [1]. The impact assessment is not associational but causal [1, 2] which requires an understanding of the counterfactual—the difference in outcomes with or without the presence of the policy [3]. This is also true for any post policy evaluation [1]. A policy maker may design multiple policies and calculate the causal quantities including the effect of the proposed policies on different recipient groups, effect over time, possible tradeoffs between competing goals, and finally, choose the optimal policy [4]. The gold standard approach for calculating those causal quantities is conducting a randomized experiment [5, 6, 7, 8]
. In a randomized experiment, the experimenter can assign an observation to either a treatment or control group randomly; this randomness can avoid bias and eliminate confounding effects of covariates and thus can achieve unbiased estimation of treatment effects. In this case, a possible association between treatment and outcome will imply causation. However, many studies in health care, social science, economics, and epidemiology cannot be designed as a randomized experiment due to legal or ethical reasons; randomization can also be impractical, time consuming, or very expensive. Hence, those experiments are performed on data that is collected as a natural process. Such experiments are called observational studies (also refereed as natural experiments or quasiexperiments)
[9] and can be implemented in a prospective (collecting sample data as natural observation over time) or retrospective (experimenting on already collected data) way.Making causal inference from an observational study lacks the experimental elements of randomization on all possible background covariates (the observed and unobserved characteristics of a sample unit) [10, 11] and prone to bias, and systematic confounding on covariates. However, with proper understanding of the underlying process and careful control of nonrandomized data, it is possible to make reasonable estimation of the causal effect [5]. Researchers have been utilizing matching methods for identifying causality since the 1940s [10] and it is one of the most popular methods. It was used or referred in as many as 486,000 academic articles involving causal inference (see LABEL:s:1). Matching methods examine the possibility of restoring or replicating properties of randomization based on the observed covariates [10]. In fact, matching attempts to retrieve the latent randomization within the observational data [12]
. Being true to it’s name, matching methods aim to find a control group which is identical to the treatment group in terms of joint distribution of the observed covariates. As discussed by Stuart
[10], and Zubizarreta [13], matching the empirical distribution of the covariates has several significant advantages. For example, matching forces the experimenter to closely examine the data, checking the common support on the covariates, and make the experimenter aware on the quality of inference. Even though the matching process can be complex, the outcome analysis is often done with simple methods [14]. For instance, the Rubin Causal Model (also known as Potential Outcome Framework) estimates the causal effect as the difference of expected outcomes between the control group and the treatment group [15]. Due to it’s simplistic architecture and many other attractive properties (see [10, 13, 16]), matching has been used to make policy decision or policy evaluation in health care [17, 18, 19, 20], education [21, 22], economics [23], law [24], and politics [25].In this paper, we adopt a robust methodology recently proposed by Morucci et al. [26] and extend it to accommodate causal inference from big data observational studies. We show the efficiency of the proposed methods by evaluating the impact of Hospital Readmission Reduction Program (HRRP) [27] policy on Nonindex readmission—readmission to a hospital that is different from the hospital which discharged the patient.
Motivation and contribution
Motivation
The objective of current onetoone matching paradigm under potential outcome framework is to find pairs between samples from treatment group and from control group . A pair is assigned in such a way that and are exactly same or similar on a specific, predetermined set of covariates : . Over the years, researchers developed a wide array of algorithms to find such pairs, for example, Propensity Score matching [14], Mahalanobis Distance matching [14], Nearest Neighbour Greedy matching [28], Coarsed Exact Matching [29], and Genetic matching [30] are among the most popular algorithms. All these algorithms (including those not listed here) disregard the outcomes of corresponding pairs in the assignment process. Though the matching process reduces bias in treatment effect estimation, disregarding the outcomes in the assignment process introduces a new source of uncertainty. If a sample has multiple possible pair assignments and have similar covariate balance but different outcomes (i.e., ), by assigning pairs without considering the outcomes, an experimenter can estimate multiple degrees of causal effect (one for each possible assignment). Similarly, a sample from control group can have multiple possible assignment options . A possible scenario is presented in Fig 1 where within each circle we have multiple pair assignment options with almost similar match quality but different outcomes (outcomes are presented as the size of the data points). In that case, different experimenters using different matching algorithms can get different pairs, hence, their causal effect estimates and conclusions on the experiment can be different. It is possible that two researchers having the exact same hypothesis and using the exact same data but, two different matching algorithms can achieve completely opposite results due to this uncertainty. This problem is exacerbating for studies involving big data as we may have more pair assignment options. Therefore, making policy decisions, in health care or any other field, using matching method that disregards uncertainty due to pair assignment can make a harmful impact on the society.
On the other hand, in 2012, congress adopted HRRP under Patient Protection and Affordable Care Act (PPACA) [27] to increase quality of care and reduce hospital readmission rate. HRRP penalizes the hospital which discharged the patient (index hospital) if the patient with some specific conditions (i.e., Pneumonia, Acute Mayocardial Infraction (AMI), Congestive Heart Failure (CHF)) is readmitted within 30 days of discharge. The index hospital is always penalized even if the patient is readmitted to different hospital (nonindex hospital) [27]. Though readmission to the index hospital is following a decreasing trend over post HRRP periods, nonindex readmissions are increasing [31]. This increase in nonindex readmission rate—approximately one fifth of all readmissions for medicare patients [31, 32]
—creates suspicion that hospitals are possibly discouraging patients for readmission to avoid penalties introduced by HRRP. Moreover, a recent study identified that nonindex readmissions are associated with higher odds of inhospital mortality and longer length of stay
[33]. Therefore, we aim to identify whether HRRP has a causal relation to the increase in nonindex readmissions. Finding such causal relation involves analyzing large volume of health care data and matching method would be vulnerable to the uncertainty discussed above. The robust method proposed in [26] to handle such uncertainty requires solving multiple Integer Programming (IP) models (a minimization and a maximization problem) iteratively. Using stateoftheart integer programming solvers to solve those IP models for big data observational studies will be computationally expensive.Contribution
In this work, we extend the robust causal inference testing method proposed by Morucci et al. [26]
to handle large scale observational studies with binary outcomes. To handle big data, first, we propose a robustness condition which identifies when a robust solution is possible and combines the maximization and minimization problems into a single problem. Second, we propose an efficient algorithm to calculate the test statistics for the robust condition. In addition, we propose two algorithms—one to solve the minimization problem and one to solve the maximization problem—for any condition which will show the degree of uncertainty for a selected number of matched pair. Finally, we implement the algorithms by testing causal effect of HRRP to nonindex readmission using the State of California Patient Discharge Data and compare the computational efficiency with canonical IP solvers.
Remark 1.
Please note, by “Robust” we imply “Robust to the choice of matching method”: if represents a set of all possible matching algorithms, a researcher choosing any algorithm and testing a hypothesis of causal effect will get the same result if she has chosen algorithm . Also, we are considering matching as preprocessing and plan to achieve robust test result from a largescale observational study for a given set of good matches identified by any matching algorithm .
Causal inference with matching method and robust test
In the Rubin Causal Model, a sample unit from a set of observations can have two outcomes or responses. The response is called treatment response when the unit receives certain treatment () and control response when unit does not receive treatment (). It is assumed that the treatment assignment to any unit does not interfere with the outcome of other units [34]. This assumption is commonly known as Stable Unit Treatment Value Assumption (SUTVA). Under this assumption, the treatment effect on a sample unit is calculated as . However, it is impossible to observe the counterfactual scenario for the same sample [15]. Under a certain treatment regime and identical conditions, we can only observe or for sample : [5, 15]. Therefore, we cannot directly measure the treatment effect at an individual level. On the other hand, causal inference literature offers a statistical solution to this fundamental problem by taking expectation over the observation set , formally called Average Treatment Effect (ATE).
(1) 
The ATE as defined in Eq 1 provides the opportunity to divide into the treatment group when and control group when such that and work with their expectations. So, we can construct the ATE as but, this form of ATE implicitly assumes that the potential responses are independent of treatment assignment: . Though this independence assumption holds in randomized experiments, in general, it does not hold for observational studies as the experimenter rarely has control on the treatment assignment process. This problem is solved by making an assumption known as Strong Ignorability [7]. Let and be the set of pretreatment background variables (covariates) which characterizes the observations. The strong ignorability assumption states that the potential responses are independent of treatment assignment when conditioned on the covariates: and every unit
has a positive probability to receiving treatment:
. Another commonly used estimate of causal effect is Average Treatment Effect on Treated (ATT) which is defined under slightly relaxed assumption ().(2) 
Both of these estimates are prone to bias as the treatment assignment process is not random. In matching method, an unbiased estimate of causal inference can be achieved if treatment unit is exactly matched with a control unit in terms of their covariate set [7]. However, in most of the applications, it is impossible to achieve exact matching [13, 7, 35, 36]. A wide variety of matching methods are employed to make pairs as similar as possible [13, 37, 7] or finding a subset of control group samples which is similar to the treatment group samples in the joint distribution of the covariate set [35, 29]. In this work, we consider onetoone matching which aims to find a pair that is matched (either exactly or by some user defined balance function) on a set of covariates .
Before explaining the difference between classical method of causal inference [5, 14, 15] and robust causal inference testing approach [26], let us define the set of good match and the pair assignment variables .
Definition 1.
(A set of Good Match) A set of good match includes treatment group samples and control group samples that satisfies certain covariate balance criteria defined under matching algorithm .
Definition 2.
(Pair Assignment Operator) The Pair Assignment Operator is a binary assignment variable where if sample is paired with a sample and the pair ; otherwise.
For a given set of possible matches
, we can perform hypothesis test in the following form with the null hypothesis being no causal effect and alternative being the opposite.
(3) 
(4) 
Under the classical approach of matching method, we can test these hypotheses first by defining a test statistic , specifying an imbalance measure along with a tolerance limit on the imbalance. Then, we apply a matching algorithm to find the set of good match which satisfy the imbalance limit, otherwise we tune the allowable imbalance limit to generate . Robust approach differs from the classical approach moving forward from here (see Fig 2). The classical approach picks one (out of many) possible combination of pairs from and conducts the hypothesis test wherein, the robust approach calculate the maximum and minimum value of the test statistic and corresponding pvalues to explore all possible assignment combinations within which does not increase imbalance under Definition 1. The test will be robust if both and produce same conclusion on the hypothesis. We formally define the Robust Test in Definition 3.
Definition 3.
(Robust Test) Let be the level of significance set for the hypothesis and are the test statistics calculated from , then, testing is called robust if max(pvalue(), pvalue() or min(pvalue(), pvalue(). Testing is called absoluterobust when pvalue() pvalue().
Calculating the test statistic generates an integer programming model which is computationally expensive for large scale data. In the following section, we propose a robustness condition following the Robust Test definition which will allow us to calculate a for absoluterobust test and we can avoid solving two integer programming problems. In this work, we are interested in testing the hypothesis stated in Eqs (3,4) for binary outcomes: with the McNemar’s test [38] as proposed in [26]. From now on, we use the terms “robust” and “absoluterobust” interchangeably and when we use robust, we are implying to the absoluterobust test defined in Definition 3.
Robust McNemar’s test
McNemar’s test is the ideal candidate for testing hypothesis in Eqs (3,4) as it deals with ontoone matched pairs. It operates on a contingency table (see Table 1) and the test statistics under null hypothesis assumes that the marginal proportions are homogeneous. Among the four types of matched pairs, we are mainly interested in the discordant pairs and where is the pair assignment operator defined in Definition 2. Here, counts the number of pairs where treatment units has outcomes 0: and control units has outcomes 1: and counts the discordant pairs where and . Under the assumption of having at least 1 discordant pair: we will use the test statistic as defined in Eq (5) to test both hypotheses.
(5) 
Treatment  

Yes  No  
Yes  A  B  
Control  No  C  D 
Morucci et al. [26] proposed the following integer programming model that explores all possible assignment options and calculate maximum and minimum possible test statistics and , respectively.
Subject to:
(6)  
(7)  
(8)  
(9)  
(10) 
(11) 
As it is shown in Table 1, is the total number of untied responses when is matched with . Similarly, is total number of untied responses when is matched with . Therefore, both . Under this definition of and , we provide the following propositions on the objective function of the robust McNemar’s test and its optimal values.
Proposition 1.
The objective function , has the following properties:

For any , is strictly increasing in for

For any , is monotonically decreasing in for and strictly decreasing for
Proof.
Let , then for any , we have
(12) 
which implies is strictly increasing in for a fixed . Similarly, let , then for any we have,
(13) 
this proves the claims of Proposition 1. ∎
Before further discussion on and the optimality conditions, we introduce following notations and definitions of maximum untied responses, for both and .
is the number of treatment units in the matched set with positive outcome
is the number of treatment units in the matched set with negative outcome
is the number of control units in the matched set with positive outcome
is the number of control units in the matched set with negative outcome
Definition 4.
(Maximum type one discordant pair) in the maximum number of possible pairs between and where the treated observation has negative (“No”) outcome but the untreated (control) observation has positive (“Yes”) outcome, i.e.,
Definition 5.
(Maximum type two discordant pair) in the maximum number of possible pairs between and where the treated observation has positive (“Yes”) outcome but the untreated (control) observation has negative (“No”) outcome, i.e.,
For a fix value of , the McNemar’s test model becomes linear and the objective functions become,
(14) 
Using the property of explained in Proposition 1, we can find the optimal solution.
Proposition 2.
Let and denote as the total number of discordant pairs, then the optimal pair is given by:
Proof.
From Proposition 1, we know that is monotonically decreasing in C when . Therefore, in the minimization problem, assignment will be made to maximize until we are about to violate constraint . When the total number of discordant pairs is set to , will take the value of and will take the value of just to satisfy the total number of discordant pair constraints and the solution will be optimal. If , the new and the minimum value will be achieved at and .
Similarly, from Proposition 1, we know that is strictly increasing in for any . So, in the maximization problem, pair assignment will be made to maximize within the feasible region. When the total number of discordant pair is set to , at optimal solution, will take the value of and will take the value of just to stay in the feasible region. When , the will take the value and the maximum value will be achieved at and .
∎
Proposition 3.
For the linear model, an absoluterobust estimate will be achieved if and only if the total number of discordant pair .
Proof.
According to the proposed approach to the causal inference estimate, an absoluterobust estimate is achieved when and is equal. For the McNemar’s test model, the model becomes infeasible when is set to as we can only have number of total untied responses. So feasible range of is: .
To prove the Proposition 3, we first set to it’s maximum value . Using Proposition 2, in this case, the optimal solution for the problem is: , and the optimal solution for problem is: , . So, for case, we get and the solution is absoluterobust.
Conversely, can take any integer value in the range which can lead to the following six cases. For each of the cases, we will find the optimal solution using Proposition 2.

: The optimal solution for the minimization problem is and the maximization problem is .

: The optimal solution for the minimization problem is and the maximization problem is .

: The optimal solution for the minimization problem is and the maximization problem is .

: The optimal solution for the minimization problem is and the maximization problem is .

: The optimal solution for the minimization problem is and the maximization problem is .

: The optimal solution for the minimization problem is and the maximization problem is .
For all of the above six cases, , hence, the solution is not absoluterobust.
Therefore, the total number of discordant pairs have to be to get an absoluterobust estimate. ∎
As we can see from the Proposition 2, the optimization problem has become a counting problem and can be solved efficiently for big data. However, the optimal solution under Proposition 2 disregards the assignment constraints Eqs (9,10) and addition userdefined constraints. To find the optimal solution using the result from Proposition 2 that is feasible, we take a twostep approach. At the first step, we handle the userdefined constraints to find a good set of match as a preprocessing step. We can use any offtheshelf matching algorithm for that purpose or define a separate pair assignment model with different covariate balance measure to find . At the second step, we partition the set of good match into partitions such that within a partition , any treatment unit can be matched with any control unit . A formal definition of a partition is provided below.
Definition 6.
(Partition of ) is a partition if any treatment unit is a good match to any control unit and . Reverse is also true.
Construction of partitions under Definition 6 ensures that only good matches are considered for assignment. In addition, Definition 4 calculates by pairing negative outcomes of treatment units and positive outcomes of control units which inherently satisfies the pair assignment constraints Eqs (9,10). Similarly, we calculate by assigning a pair between samples with positive treatment outcomes and negative control outcomes. Therefore, none of the treatment or control unit is used more than once in the pair assignment process which satisfies the pair assignment constraints Eqs (9,10).
Now, using the above mentioned results, we propose Algorithm 1 which identifies the robustness condition and corresponding test statistics .
By the sketch of the Algorithm 1, it seems like we are only matching the discordant pairs and ignoring the other possible pair assignments in the data, which is not true. In Proposition 4, we show that we match maximum possible pairs.
Proposition 4.
Algorithm 1 ensures that the maximum possible pairs are matched in .
Proof.
To prove this Proposition, we only need to show that in any partition , Algorithm 1 matches maximum possible pairs. Then, we can sum the maximum pair assignments across the partitions to achieve maximum possible pairs assignment in .
Lets consider a partition where denotes the number of treatment samples and denotes the number of control samples. Hence, the maximum number of pairs we can assign in is . We will use to represent the number of samples with positive outcomes () and to represent the number of samples with negative outcomes (). After assigning the discordant pairs () and () as we did in Algorithm 1, we are left with treatment samples and control samples. Now, we can assign the remaining treatment and control samples into the other two types of pairs and to their limit:
It is trivial to show that the for partition ,
An example of maximum pair assignment is provided in Fig 3 where treatment outcomes () are sorted in descending order and control outcomes () are sorted in ascending order. In the left panel, we can have pairs at maximum. After assigning and according to Algorithm 1, we can assign only one pair to and . Therefore, we achieve the maximum number of pair assignments. We follow the similar procedure in the middle and right panels. ∎
In the Algorithm 1, we calculate the test statistics for robust condition which an experimenter can use to find the corresponding pvalue and compare with a predefined level of significance to make conclusion on the hypothesis of no causal relation. Decisions made in this process will be free of uncertainty and robust to the choice of matching algorithms. If different experimenters perform matching on same data using different matching algorithms but follow the above mentioned procedure, all of their conclusions will be exactly the same.
Regarding the computational complexity arises due to big data, our proposed algorithm only involves counting elements in vectors and few algebraic operations. The counting processing can be done with the summation of vectors as we are dealing with only binary outcomes: summation implies the total number of positive outcomes and we can calculate the negative outcomes by subtracting it from the size of the vector. In addition, we only need to solve the problem once—at robustness condition. Therefore, the proposed algorithm will be highly efficient for big data.
While the Algorithm 1 directly calculates test statistics at robustness condition, a researcher might be interested in exploring the degree of uncertainty in the causal inference test. She may want to see how the uncertainty changes towards the robust estimate with respect to the number of discordant pairs matched. For this purpose, we propose the following two algorithms (2, 3) following the result of Proposition 2.
Numerical experiment
In this section, we present the efficiency of the proposed algorithms with data from the State of California Patient Discharge Database and answer an interesting hypothesis on the effectiveness of HRRP implemented in October 1, 2012.
Readmission rate is considered an important measure of hospital or care quality. To increase the care quality and hold hospitals accountable, US Congress introduces HRRP under PPACA in 2012 [27]. The most important feature of this program is that the index hospital (the hospital that discharged the patient) is penalized if patients with Pneumonia, Congestive Heart Failure (CHF), and Acute Mayocardial Infraction (AMI) are readmitted (to the index hospital or any other hospital) within 30 days of discharge. During the post HRRP period, the overall rate of readmission is decreasing which the proponents of HRRP are attributing as the success of the policy. However, in this period, readmission to different hospitals (nonindex readmission) has been increasing [31, 32]. Nonindex readmissions are considered as a discontinuity of care and found to be associated with longer lengths of stay and higher inhospital mortality rates. There is a chance that hospitals are possibly discouraging patients for readmission to avoid penalties introduced by HRRP. To check the cause behind the increase in nonindex readmission post HRRP periods, we make the following hypothesis and test it with the proposed algorithms with the level of significance .
H_{0}: HRRP has no causal relation with the increase in nonindex readmission
H_{1}: HRRP has a positive causal relation with the increase in nonindex readmission
Data description and covariate balance
In this research, we primarily used the patient discharge data between 2010 to 2014 from California. We collected this nonpublic data from the California Office of Statewide Health Planning and Development (OSHPD) which collects inpatient data from California licensed hospitals. Each patient in this data set is presented with a unique identifier which is used to determine if a patient is readmitted in the future. In addition, the data set also contains patient level information such as, ICD9 codes for clinical diagnosis, comorbidity index, age, gender, discharge destination, patients’ zip code, and insurance information. The patients are tracked with unique identifiers to find out if they were readmitted within 30 days of discharge and when a readmission is found, we identify the destination hospital of that readmission. Then, a binary variable is created with 0 if the patient is readmitted to the same hospital or 1 if different hospital. To test the hypothesis, we will use this variable as our outcome:
if readmitted to a different hospital or if readmitted to the same hospital.Moreover, the OSHPD data set is merged with publicly available data from Centers for Medicare and Medicaid Services, American Association Annual Hospital Survey and the Area resource file. From these additional data sources, we get important hospital level information like, teaching status (membership status of the Council of Teaching Hospitals), ownership type (public, nonprofit, investor owned), hospital size based on number of beds (small: below 100 beds, medium: 101 to 399 beds, and large: 400 and above beds), hospital locations (rural, metro). We also include patient household incomes which we consider as the median income of patients’ residence zip codes. We divide the data into two sets: before and after October 1, 2012, the implementation date of HRRP. The treatment here is the implementation of HRRP, all the readmissions between February 1, 2010 and September 30, 2012 is considered as control group (treatment ) and all the readmission from October 1, 2012 to November 30, 2014 is considered as treatment group (treatment ). To capture any potential readmission within 30 days period of an index discharge, admissions before February 1, 2010 and beyond November 30, 2014 are excluded. A descriptive view of readmitted patients’ characteristics is presented in Table 2.
1.3 in0in
Variable  All Readmission  Index Hospital  Nonindex Hospital  Before HRRP  After HRRP 

Readmitted Patients  90553  67341  23212  53353  37200 
Demographic Characteristics  
Age  
020  635 (0.70)  505 (0.75)  130 (0.56)  427 (0.8)  208 (0.56) 
2130  1073 (1.18)  717 (1.06)  356 (1.53)  566 (1.06)  507 (1.36) 
3140  2186 (2.41)  1471 (2.18)  715 (3.08)  1269 (2.38)  917 (2.47) 
4150  6336 (7.00)  4196 (6.23)  2140 (9.22)  3714 (6.96)  2622 (7.05) 
5165  21018 (23.21)  14470 (21.49)  6548 (28.21)  11950 (22.4)  9068 (24.38) 
65 and above  59305 (65.49)  45982 (68.28)  13323 (57.4)  35427 (66.4)  23878 (64.19) 
Gender  
Female  45124 (49.80)  34240 (50.80)  10884 (46.9)  27049 (59.94)  18075 (40.06) 
Male  45429 (50.20)  33101 (49.20)  12328 (53.1)  26304 (57.9)  19125 (42.1) 
Household Income  
Quartile 1 
22428 (24.77)  15640 (23.23)  6788 (29.24)  13108 (24.57)  9320 (25.05) 
Quartile 2  22629 (24.99)  16633 (24.70)  5996 (25.83)  13331 (24.99)  9298 (24.99) 
Quartile 3  22450 (24.79)  17160 (25.48)  5290 (22.79)  13165 (24.68)  9285 (24.96) 
Quartile 4  23046 (25.45)  17908 (26.59)  5138 (22.14)  13749 (25.77)  9297 (24.99) 
Clinical Characteristics  
Primary Diagnosis  
CHF  50151 (55.40)  37404 (55.50)  12747 (54.9)  29351 (55.01)  20800 (55.91) 
AMI  11917 (13.20)  8148 (12.10)  3769 (16.2)  6865 (12.87)  5052 (13.58) 
Pneumonia  28485 (31.40)  21789 (32.40)  6696 (28.9)  17137 (32.12)  11348 (30.51) 
Charlson Comorbidity Index  
Low (02)  35394 (39.09)  25884 (38.44)  9510 (40.97)  21282 (39.89)  14112 (37.94) 
Medium (36)  51301 (56.65)  38454 (57.10)  12847 (55.35)  29850 (55.95)  21451 (57.66) 
Medium High (710)  3396 (3.75  2628 (3.90)  768 (3.31)  1944 (3.64)  1452 (3.9) 
High (10 and above)  462 (0.51)  375 (0.56)  87 (0.37)  277 (0.52)  185 (0.5) 
Hospital Characteristics  
Teaching Status  
Teaching Hospital  10261 (11.30)  7706 (11.40)  2555 (11)  5882 (11.02)  4379 (11.77) 
Nonteaching Hospital  80272 (88.70)  59635 (88.50)  20657 (89)  47471 (88.98)  32821 (88.23) 
Ownership Type  
Nonprofit Hospital  58592 (64.70)  45210 (67.10)  13382 (57.6)  34252 (64.2)  24340 (65.43) 
Investor Hospital  17902 (19.80)  11389 (16.90)  6513 (28.1)  10839 (20.32)  7063 (18.99) 
Public Hospital  14059 (15.50)  10742 (16.00)  3317 (14.3)  8262 (15.49)  5797 (15.58) 
Hospital Size  
Small (below 100 beds)  4982 (5.50)  3453 (5.10)  1529 (6.6)  2979 (5.58)  2003 (5.38) 
Medium (100399 beds)  61167 (67.60)  45149 (67.10)  16018 (69)  36314 (68.06)  24853 (66.81) 
Large (400 and above beds)  24404 (26.90)  18739 (27.80)  5665 (24.4)  14060 (26.35)  10344 (27.81) 
Hospital Location  
Rural  2426 (2.68)  1809 (2.69)  617 (2.66)  1348 (2.53)  1078 (2.90) 
Metro  88127 (97.32)  65532 (97.31)  22595 (97.34)  52005 (97.47)  36122 (97.10) 
The entries in each cell is presented as Number of patients “N (%)” form. From February 1, 2010 to September 30, 2012 is considered “Before HRRP” period. From October 1, 2012 to December 31, 2014 is considered “After HRRP” period. CHFCongestive Heart Failure, AMIAcute Mayocardial Infraction.
We match the patients based on the following covariates: Age, Gender, Primary diagnosis, Household Income, Charlson Comorbidity Index, Hospital Location, Hospitals’ Teaching Status, Hospitals’ Ownership Status, and Hospital Size. We divide the covariates into two groups: 1) discrete and 2) continuous. The discrete covariates (i.e., Gender, Primary Diagnosis, Hospital Location, Hospitals’ Teaching Status, Hospitals’ Ownership Status, and Hospital Size) are matched exactly. The continuous covariates (i.e., Age, Household Income, and Charlson Comorbidity Index) are first, divided into categories as shown in Table 2, then, the categories are matched exactly. This matching strategy resulted in 1822 partitions of data; within a partition, any treatment sample can be matched with any control sample. Though this matching approach seems ad hoc in nature, it is very similar to the well known method called Coarsed Exact matching (CEM) [29] with 1822 bins. Traditionally, CEM is implemented with much lesser number bins due to the lack of common support between treatment and control groups but a higher number of bins make finer covariate balance [39, 40] which is the objective of any matching method. However, to implement the proposed algorithms, an experimenter is not limited to CEM or the matching method we used. Given a good set of matchec created under Definition 1, we can always create the partitions under Definition 6.
Experiment and result
To test the hypothesis H_{0}, first, we perform the matching operations in R [41] to get matched sets. Then, using the matched sets of data, we calculate the test statistic and using the 1) optimization model with an Integer Programming solver and 2) with the proposed algorithms. The integer programming model is implemented in AMPL [42] and solved with the commercial solver CPLEX [43]. We implement the Algorithm 1, 2, and 3 in R [41]. All the experiments are performed in a Dell Precision workstation with 64 GB RAM, Intel(R) Xeon(R) CPU E52670 v3 processor running at 2.30 GHz.
Table 3 shows the comparison of solutions obtained using an optimization model with CPLEX iterating over different values of discordant pairs () and proposed algorithm at robustness condition. The range of pvalue achievable corresponding to the test statistics and is presented in Fig 4. The proposed Algorithm 1 directly identifies the robustness condition which is and and calculates the test statistics . The computation time required by the proposed algorithms are very insignificant compared to the time required by the optimization model. A major implication of the robust McNemar’s test is that if the same experiment is conducted with as many as 19,000 discordant pairs, we can achieve any pvalue between 0 to 0.23 (see Fig 4); some experimenter might rejected the hypothesis and some might fail to reject the null hypothesis. Both the experimenters, in this case, are right but the conclusion differed due to the fact that they choose different pairs. Any policy decision made using the matching method without considering this uncertainty has a possibility to fail.
In regards to the hypothesis we made at the beginning of this section, we can make a conclusion on the hypothesis by using the pvalue calculated at the robustness condition or the result from Table 3 and corresponding pvalues from Fig 4. We can see that, both the maximum and minimum pvalue when we match more than 20,000 discordant pairs. Therefore, we can reject the null hypothesis of no causal effect and conclude that HRRP is a cause for increase in the nonindex readmissions. This result points out that not only the readmission rate but also the nonindex readmission rate should be considered as a measure of health care quality and health care policy makers should take another look at the implication of HRRP policy.
Conclusion
Any policy decision or evaluation requires identifying the causal relation between policy alternatives and potential outcomes. Matching methods have become very popular in identifying such causal relations. However, in onetoone matching, when we have multiple pair assignment options, matching method is vulnerable to uncertainty as the pair construction process does not consider outcomes. In this paper, we consider the integer programming model for robust causal inference testing approach with binary outcomes proposed by Morucci et al. [26] and develop scalable algorithms that can be used for largescale observational studies. We identify a robustness condition which combines the maximization and minimization problem proposed in [26]. Instead of solving two computationally expensive integer programming models iteratively by increasing the number of discordant pairs until a robust estimate is achieved, we convert the problems into counting problems through a series of propositions. In addition, the proposed Algorithm 1 only solves one problem instead of two separate problems and it is computationally efficient. Numerical experiment conducted on State of California Patient Discharge Database shows that the proposed algorithms are highly scalable. The numerical experiment shows an interesting result on a well appreciated health care policy—HRRP—proposed under PPACA in 2012 to improve health care quality. We identify that HRRP is a cause to the increase of nonindex readmission which is associated to higher inhospital mortality rate and longer length of stay. Though the numerical experiment is performed with around 100,000 samples, the algorithms proposed in this paper can handle observational studies with millions of samples efficiently without further modification. In the future, we plan to develop similar robust causal inference testing algorithms with continuous outcomes for largescale observational studies.
Acknowledgments
We thank Mr. Tasnim Ibn Faiz, PhD candidate, MIE, Northeastern University for his help and advice on programming in AMPL. We also thank Mr. Md Mahmudul Hasan, PhD candidate, MIE, Northeastern University for sharing his knowledge on the OSHPD database and supporting us throughout the data analysis process.
References
 [1] Boniface Essama Nssah. Propensity score matching and policy impact analysis: A demonstration in eviews, volume 3877. World Bank Publications, 2006.
 [2] Judea Pearl. Causal inference in statistics: An overview. Statistics surveys, 3:96–146, 2009.
 [3] Jon Kleinberg, Jens Ludwig, Sendhil Mullainathan, and Ziad Obermeyer. Prediction policy problems. American Economic Review, 105(5):491–95, 2015.
 [4] Tristan Zajonc. Essays on causal inference for public policy. PhD thesis, 2012.
 [5] Donald B Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of educational Psychology, 66(5):688, 1974.
 [6] Peter C. Austin. An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behavioral Research, 46(3):399–424, 2011. PMID: 21818162.
 [7] Paul R. Rosenbaum and Donald B. Rubin. The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1):41–55, 1983.
 [8] Susan Athey and Guido W Imbens. The state of applied econometrics: Causality and policy evaluation. Journal of Economic Perspectives, 31(2):3–32, 2017.
 [9] Paul R Rosenbaum. Observational studies. In Observational studies, pages 1–17. Springer, 2002.
 [10] Elizabeth A. Stuart. Matching methods for causal inference: A review and a look forward. Statist. Sci., 25(1):1–21, 02 2010.
 [11] Donna F Stroup, Jesse A Berlin, Sally C Morton, Ingram Olkin, G David Williamson, Drummond Rennie, David Moher, Betsy J Becker, Theresa Ann Sipe, Stephen B Thacker, et al. Metaanalysis of observational studies in epidemiology: a proposal for reporting. Jama, 283(15):2008–2012, 2000.
 [12] Ben B Hansen. Full matching in an observational study of coaching for the sat. Journal of the American Statistical Association, 99(467):609–618, 2004.
 [13] José R Zubizarreta. Using mixed integer programming for matching in an observational study of kidney failure after surgery. Journal of the American Statistical Association, 107(500):1360–1371, 2012.
 [14] Paul R. Rosenbaum and Donald B. Rubin. Constructing a control group using multivariate matched sampling methods that incorporate the propensity score. The American Statistician, 39(1):33–38, 1985.
 [15] Paul W. Holland. Statistics and causal inference. Journal of the American Statistical Association, 81(396):945–960, 1986.
 [16] Stephen L Morgan and David J Harding. Matching estimators of causal effects: Prospects and pitfalls in theory and practice. Sociological methods & research, 35(1):3–60, 2006.
 [17] Nicholas A Christakis and Theodore J Iwashyna. The health impact of health care on families: a matched cohort study of hospice use by decedents and mortality outcomes in surviving, widowed spouses. Social science & medicine, 57(3):465–475, 2003.
 [18] Yuji Akematsu and Masatsugu Tsuji. Measuring the effect of telecare on medical expenditures without bias using the propensity score matching method. Telemedicine and eHealth, 18(10):743–747, 2012.
 [19] Astrid Kiil. Does employmentbased private health insurance increase the use of covered health care services? a matching estimator approach. International journal of health care finance and economics, 12(1):1–38, 2012.
 [20] Nazmi Sari and Meric Osman. The effects of patient education programs on medication use among asthma and copd patients: a propensity score matching with a differenceindifference regression approach. BMC health services research, 15(1):332, 2015.
 [21] José R Zubizarreta and Luke Keele. Optimal multilevel matching in clustered observational studies: A case study of the effectiveness of private schools under a largescale voucher system. Journal of the American Statistical Association, 112(518):547–560, 2017.
 [22] Guanglei Hong and Stephen W Raudenbush. Evaluating kindergarten retention policy: A case study of causal inference for multilevel observational data. Journal of the American Statistical Association, 101(475):901–910, 2006.
 [23] Rajeev H Dehejia and Sadek Wahba. Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs. Journal of the American statistical Association, 94(448):1053–1062, 1999.
 [24] Lee Epstein, Daniel E Ho, Gary King, and Jeffrey A Segal. The supreme court during crisis: How war affects only nonwar cases. NYUL rev., 80:1, 2005.
 [25] Michael C Herron and Jonathan Wand. Assessing partisan bias in voting technology: The case of the 2004 new hampshire recount. Electoral Studies, 26(2):247–261, 2007.
 [26] Marco Morucci, Md NoorEAlam, and Cynthia Rudin. Hypothesis tests that are robust to choice of matching method. arXiv preprint arXiv:1812.02227, 2018.
 [27] Colleen K McIlvennan, Zubin J Eapen, and Larry A Allen. Hospital readmissions reduction program. Circulation, 131(20):1796–1803, 2015.
 [28] Peter C Austin. Some methods of propensityscore matching had superior performance to others: results of an empirical investigation and monte carlo simulations. Biometrical Journal: Journal of Mathematical Methods in Biosciences, 51(1):171–184, 2009.
 [29] Stefano Iacus, Gary King, Giuseppe Porro, et al. Cem: software for coarsened exact matching. Journal of statistical Software, 30(13):1–27, 2009.
 [30] Alexis Diamond and Jasjeet S Sekhon. Genetic matching for estimating causal effects: A general multivariate matching method for achieving balance in observational studies. Review of Economics and Statistics, 95(3):932–945, 2013.
 [31] Min Chen. Reducing excess hospital readmissions: Does destination matter? International journal of health economics and management, 18(1):67–82, 2018.
 [32] Min Chen and David C Grabowski. Hospital readmissions reduction program: intended and unintended effects. Medical Care Research and Review, 2017.
 [33] Robert E Burke, Christine D Jones, Patrick Hosokawa, Thomas J Glorioso, Eric A Coleman, and Adit A Ginde. Influence of nonindex hospital readmission on length of stay and mortality. Medical care, 56(1):85–90, 2018.
 [34] Donald B. Rubin. Bayesian inference for causal effects: The role of randomization. Ann. Statist., 6(1):34–58, 01 1978.
 [35] Alexander G. Nikolaev, Sheldon H. Jacobson, Wendy K. Tam Cho, Jason J. Sauppe, and Edward C. Sewell. Balance optimization subset selection (boss): An alternative approach for causal inference with observational data. Operations Research, 61(2):398–412, 2013.
 [36] Gary King and Richard Nielsen. Why propensity scores should not be used for matching. https://gking.harvard.edu/publications/whyPropensityScoresShouldNotBeUsedFormatching (Under Review), 2016.
 [37] José R Zubizarreta. Stable weights that balance covariates for estimation with incomplete outcome data. Journal of the American Statistical Association, 110(511):910–922, 2015.
 [38] Quinn McNemar. Note on the sampling error of the difference between correlated proportions or percentages. Psychometrika, 12(2):153–157, 1947.
 [39] Stefano M Iacus, Gary King, and Giuseppe Porro. Causal inference without balance checking: Coarsened exact matching. Political analysis, 20(1):1–24, 2012.
 [40] Stefano M. Iacus, Gary King, and Giuseppe Porro. Multivariate matching methods that are monotonic imbalance bounding. Journal of the American Statistical Association, 106(493):345–361, 2011.
 [41] Ross Ihaka and Robert Gentleman. R: a language for data analysis and graphics. Journal of computational and graphical statistics, 5(3):299–314, 1996.
 [42] Robert Fourer, David M Gay, and Brian W Kernighan. AMPL: A mathematical programming language. AT & T Bell Laboratories Murray Hill, NJ 07974, 1987.
 [43] IBM ILOG CPLEX. V12. 1: User’s manual for cplex. International Business Machines Corporation, 46(53):157, 2009.
Comments
There are no comments yet.