Graph-Based Tests for Multivariate Covariate Balance Under Multi-Valued Treatments

In studies of causal effects, balancing the covariate distributions in each treatment group is essential for unbiased estimates. For example, in a study of the effects of smoking on medical expenditures, multiple smoking behaviors are possible: there may be never smokers, former smokers, and current smokers. To estimate the relative effects of the different exposures, we would like the groups to be balanced by various demographic characteristics such as age, race, and income. Distributional balance is expected by design in randomized control trials (RCT) for both measured and unmeasured confounders, but conducting a RCT is often not possible because of either cost or ethical restrictions, as would be the case in such a study of smoking. Instead, researchers are left to rely on observational studies to answer their causal questions.

If the number of covariates is few and the covariates have a small number of discrete levels, obtaining exact matches between treatment groups may be possible, thereby achieving exact distributional balance. However, such simple settings are rarely found in practice. This problem can be further compounded in the multi-valued treatment setting when there are more groups to match. Achieving distributional balance in such multi-valued treatment settings then relies on adjustments via methods such as propensity score matching [27] or cardinality matching [35]. Intuitively, these methods work by matching units that are close in terms of some covariate distance, which will make the distributions more similar between the treatment groups. But how do we know when the distributions between treatment groups are similar enough to proceed?

To answer this question, investigators typically consider univariate assessments of covariate balance. Such assessments can include assessing the mean of each covariate by treatment group before and after matching and a corresponding test for the difference of these means [30]. Other authors run a Kolmogorov-Smirnov test [19, 29] for each covariate [33, 12]. However, these methods will not generally diagnose differences in the joint covariate distribution between treatment groups.

Instead, we seek a way to identify differences in the joint distributions between the multiple groups after matching is completed—in other words, whether the observational study fails to approximate a randomized experiment in terms of observed covariate balance. Several common methods of assessing distributional balance such as likelihood ratio tests or examining the distribution of the propensity scores [28] require having correct probability models for the data. Unfortunately, these data models are rarely known in practice.

Fortunately, there are several non-parametric tests applicable to the multi-valued treatment setting and the ones we consider involve forming a graph-like structure on the data. Some work involves ranking the data in some fashion such that the observations are ordered in a one-dimensional space and hence form a graph [5]. [25], for example, consider a ranks-based test that uses the data ranks for each individual covariate in a Kruskal-Wallis test [21] while [10] use data-depth—measuring how far each observation is from the overall sample mean in terms of Mahalanobis distance—to calculate another ranks based test. [22] also uses a rank-based test to simultaneously examine location and scale changes. Rather than ranks, work by [24] instead forms a non-bipartite graph and extends the crossmatch test of [26] to multi-valued treatments.

In this vein, we consider several graphs structures on the data that range from Hamiltonian paths connecting all of the data into a single graph to kNN graphs. Hamiltonian paths are graphs that connects all of the observations in a dataset but with the constraint that each observation is only visited once. If we imagine the path as a highway, we will see all of the observations along the highway, but there is no way to see the old points we have visited without turning around—in graph terminology, there are no cycles. There are many such paths through the data under this broad definition, but we will define the shortest of such possible paths as the optimal Hamiltonian path. Under this graph, if the distributions vary between treatment groups we expect units from the same distribution to be closer to each other along the path.

For nearest neighbor graphs, the goal is not to find a path but to find matches—i.e., neighbors—that are close to each other in terms of a distance metric. This can be done in a way that allows for overlap between matched sets or in a way that forces those matched sets to be disjoint. In both types of nearest neighbor graphs, if the distributions vary between treatment groups, we expect observations from the same treatment group to be connected more frequently to each other than to observations from different treatmentes.

Several methods use these graphs for two-sample testing. In Hamiltonian paths, [7] use the Wald-Wolfowitz runs test [32] along the path to test for distributional similarity in the two-sample setting. For nearest neighbor graphs, the crossmatch test of [26] creates non-overlapping kNN graphs for $k=1$ and counts the number of matches that occur across samples. Recent work by [9] develops a similar test on greedy kNN graphs that instead counts the number of times nearest neighbors come from the same sample.

Unfortunately, there has not been much work extending these particular graph-based tests to the multisample setting. As previously mentioned, [24] extend Rosenbaum’s crossmatch test to the multisample setting, but similar examples exist neither for Hamiltonian paths nor for greedy nearest neighbor graphs. In addition, the empirical performance of these methods is unclear in practice, and to our knowledge, there is no readily available computational implementation of these methods.

To bridge these gaps with matched samples, we offer four contributions. First, we develop similar extensions for tests along the Hamiltonian path. Second, we provide a multisample extension to the nearest neighbor test of [9]. Third, we conduct an extensive simulation study of these tests in a variety of settings, finding that the kNN test outperforms others. Fourth, we implement all of these tests into the new R package graphTest.

Next, we briefly describe our setting in Section 2, give a motivating example in Section 3, and describe how we construct and perform tests on the graphs in Section 4. Then we provide results from simulations and a case study in Sections 5 and 6, respectively. Finally, Section 7 concludes.

Suppose we have covariate data from $G$ treatment groups: $\{X_{1}^{(1)},\ldots,X_{n_{1}}^{(1)}\},$ $\ldots,$ $\{X_{1}^{(G)},$ $\ldots,$ $X_{n_{G}}^{(G)}\}$ with $n_{1}+n_{2}+\cdots+n_{G}=N$. These groups could arise from sampling from different populations or because observations from the same population have chosen different treatments. As such, our goal is to see if after matching or some other adjustment, we can successfully evaluate if the distribution in each group is the same—i.e. $H_{0}:F_{1}=F_{2}=\cdots=F_{G}$.

Constructing graph-based tests studied in this work requires a notion of closeness between units and a way to denote which units are connected. We define $\mathbf{D}$ as the pairwise distance matrix between all points in the sample with entries $\mathbf{D}_{ij}={\mathbf{d}}(X_{i},X_{j}),\,\forall i\neq j,\,i,j\in\{1,...,N\}$, for some distance $\mathbf{d}:{}\operatorname{\mathbb{R}}^{d}\times{}\operatorname{\mathbb{R}}^{d}\mapsto{}\operatorname{\mathbb{R}}_{+}$. Further, let $\mathbf{M}$ be a matrix whose entries indicate whether unit $i$ is connected to unit $j$. If unit $i$ is connected to unit $j$, then $\mathbf{M}_{ij}=1$. Otherwise, $\mathbf{M}_{ij}=0$.

To better understand this setting and the challenges it presents, we now consider a small example where individuals can choose one of three treatments. Our setting has only two confounders, $X_{1},X_{2}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny iid}}}{\sim}}\mathcal{N}(0,1)$ and we then generate a treatment indicator $Z\sim\operatorname{Categorical}(\mathbf{p})$, where $\mathbf{p}=(p_{1},p_{2},p_{3}),^{\top}\,p_{i}=\frac{\exp(\eta_{i})}{\sum_{j=1}^{3}\exp(\eta_{j})}$ and $\eta_{1}=0.1X_{1}-0.1X_{2}-X_{1}X_{2},\eta_{2}=-0.2X_{1}+0.2X_{2}+0.5X_{1}^{2},\eta_{3}=-0.1X_{1}+0.2X_{2}-2X_{1}X_{2}.$ We set the sample size to $N=150$.

Typically, researchers will not know the true treatment-assignment mechanism and we emulate this situation in our setup. We use cardinality matching [35] to balance first moments between each group and the overall sample mean but ignore the covariate functions $X_{1}^{2}$ and $X_{1}X_{2}$ necessary to capture the true propensity score model. Matching towards the overall sample means will target the average treatment effects for the whole sample [36, 2]. We use the R package designmatch [34] and set the parameters of the matching such that the group means should differ from the sample mean by no more than 0.1 standard deviations. The results of such a matching for one draw of the data is shown in Figure 1.

The covariate balance before and after matching for this draw in Table 1 would indicate that the matching was a success. Further, the absolute standardized mean difference, defined here as $\frac{1}{3}\sum_{g=1}^{3}\frac{|\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu_{g}-\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu|}{\sqrt{\widehat{\operatorname{Var}}(X)}}$, also indicates that all groups are within 0.1 standard deviations of the sample mean on average, and $F$ tests for differences in means fail to reject the null hypothesis that the group means are the same.

Unfortunately, the matching model misses two important features of treatment selection: $X_{1}^{2}$ and $X_{1}X_{2}$. Running this experiment 1000 times, we see that while standardized differences of means $X_{1}$ and $X_{2}$ are well balanced on average, the other covariate functions are not (Figure 2(a)). Similarly, univariate tests of covariate means are insufficient to diagnose the statistically significant deviations in the other covariate functions (Figure 2(b)). Clearly, we need an omnibus test to diagnose failure without having to check every possible covariate moment.

With this goal in mind, we now examine some graph-based tests that may offer a solution.

This section considers two graph structures: Hamiltonian paths and nearest neighbor graphs. For each, we describe the graph structures, the algorithms to construct them, and the tests that we will run on these structures. Throughout, we will use the motivating example of the previous section for exposition in our figures.

We now present simulation results in a variety of settings. Note that we set $k=\lfloor 0.1N\rfloor$ in each case where we use the kNN test.

In the first simulation study, we return to the motivating example from Section 3 to see if the new tests can correctly diagnose deviations from the null hypothesis that $F_{1}=\cdots=F_{G}$ missed by univariate tests. Then in the second simulation study, we utilize the settings of [24] to evaluate the performance of the graph-based tests. Overall, we find greedy paths and greedy nearest neighbor constructions perform better than their optimal counterparts while the kNN test performs best on average.

Our goal is to use matching to make the distributions in the respective treatment groups more similar to the overall sample distribution. We will then use our graph-based tests to see if univariate tests miss any distributional differences. To do so, we use data from the 1987 National Medical Expenditures Survey data as collected by the TriMatch package in R [8].

In this paper, we proposed multisample extensions to the Hamiltonian path testing of [7] and the nearest neighbor testing of [9] and we evaluated these tests in a variety of settings. We considered various graph structures, algorithms to form these structures, and tests of distributional balance.

The author would like to thank Claire Chaumont, Gang Liu, Aarón Sonabend, Lorenzo Trippa, and José Zubizarreta for helpful comments and feedback on an earlier version of this manuscript. This research was funded by generous support from NIH grant 5T32CA009337-40, the Department of Biostatistics at the Harvard T.H. Chan School of Public Health, and the David Geffen Scholarship from the David Geffen School of Medicine at UCLA.