# Efficient algorithms for building representative matched pairs with enhanced
generalizability

Bo Zhang Bo Zhang is an assistant professor of biostatistics at the Fred
Hutchinson Cancer Center. Email<EMAIL_ADDRESS>
Abstract: Many recent efforts center on assessing the ability of real-world
evidence (RWE) generated from non-randomized, observational data to produce
results compatible with those from randomized controlled trials (RCTs). One
noticeable endeavor is the RCT DUPLICATE initiative (Franklin et al., 2020,
2021). To better reconcile findings from an observational study and an RCT, or
two observational studies based on different databases, it is desirable to
eliminate differences between study populations. We outline an efficient,
network-flow-based statistical matching algorithm that designs well-matched
pairs from observational data that resemble the covariate distributions of a
target population, for instance, the target-RCT-eligible population in the RCT
DUPLICATE initiative studies or a generic population of scientific interest.
We demonstrate the usefulness of the method by revisiting the inconsistency
regarding a cardioprotective effect of the hormone replacement therapy (HRT)
in the Women’s Health Initiative (WHI) clinical trial and corresponding
observational study. We found that the discrepancy between the trial and
observational study persisted in a design that adjusted for study populations’
cardiovascular risk profile, but seemed to disappear in a study design that
further adjusted for the HRT initiation age and previous estrogen-plus-
progestin use. The proposed method is integrated into the R package match2C.

Keywords: generalizability, matching, RCT DUPLICATE Initiative; trial
emulation; Women’s Health Initiative

##  1 Introduction

###  1.1 RCT DUPLICATE Initiative: Comparing Observational Studies to
Randomized Controlled Trials (RCTs)

In a recent high-profile study published in Circulation, the RCT DUPLICATE
initiative (Franklin et al.,, 2020, 2021) designed $10$ observational studies
using retrospective, non-randomized claims data, and compared their real-
world-evidence-based (RWE-based) treatment effect estimates to those based on
$10$ randomized controlled trials (RCTs) investigating very similar clinical
questions. The RCT DUPLICATE initiative aims to build an empirical evidence
base for real world data through large-scale replication of RCTs and
understand to which clinical questions and by what analytic tools researchers
could draw credible causal conclusions from retrospective, non-randomized data
(e.g., electronic health records, administrative claims databases, and
diseases registries).

To enable a better comparison of effect estimates obtained from observational
and RCT data, it is desirable to design an observational study whose treated
and matched control groups are comparable to the RCT population in baseline
covariates. In their design stage, Franklin et al., (2021) carefully emulated
the RCT study population by applying the same inclusion and exclusion criteria
to the observational data prior to statistical matching; however, tangible and
potentially meaningful differences persist. Franklin et al., (2021) concluded:

> [I]nclusion and exclusion criteria from the trials could only be partially
> emulated, and even where fully emulated, the resulting distributions were at
> times meaningfully different between the RCT and RWE populations, possibly
> because of nonrepresentative participation in RCTs.

Discrepancies in age, race/ethnicity, and important preexisting comorbid
conditions between RCT and observational study populations are common in
Franklin et al., (2021)’s emulation studies. For instance, the Saxagliptin and
Cardiovascular Outcomes in Patients with Type 2 Diabetes Mellitus (SAVOR-TIMI
53) study (Scirica et al.,, 2013), one of the ten trials Franklin et al.,
(2021) emulated, enrolled $33.1\%$ female and $37.8\%$ with history of
myocardial infarction (MI); the observational study emulating the SAVOR-TIMI
53 trial, however, consisted of $46.8\%$ female and only $11.2\%$ with history
of MI (Franklin et al.,, 2021, Table 1). These differences persisted in the
final matched samples constructed by Franklin et al., (2021), and could
partially explain the disagreement in effect estimates derived from RCTs and
observational studies.

The ongoing RCT emulation study led by the RCT DUPLICATE initiative and many
similar endeavors to better reconcile RCT and observational study findings,
for instance, the Women’s Health Initiative study (Prentice et al.,, 2005;
Hernán et al.,, 2008), or findings derived from different observational
databases, motivate us to develop a transparent, efficient and easy-to-use
algorithm that constructs homogeneous matched pairs that resemble a template
(that is, a random sample from the target population) in key covariates.

There are two important applications for such an algorithm. First, in the
context of target trial emulation using observational data (Franklin et al.,,
2021), the target population would be the RCT-eligible population, the
template would consist of trial participants, and the algorithm could be used
to construct treated and matched control groups that resemble the target RCT
population. Second, the algorithm may be used to help assess inconsistencies
often seen among observational studies based on different databases. For
instance, in a recent study of the impact of intraoperative transesophageal
echocardiography (TEE) on patients’ clinical outcomes, Metkus et al., (2021)
conducted a matched cohort study and found an effect smaller than that derived
by MacKay et al., (2021). Metkus et al.,’s (2021) analysis was based on the
Society of Thoracic Surgeons Adult Cardiac Surgery Database (STS ACSD) while
MacKay et al.,’s (2021) analysis was based on Medicare beneficiaries who are
at higher peri-operative risk compared to the adult population in the STS
ACSD. In this example, Medicare beneficiaries would be the target population,
the template would consist of a random sample from Medicare claims data, and
our algorithm could be used to construct matched cohorts from the STS ACSD
that mimic Medicare beneficiaries in risk factors. This could enable
researchers to estimate the effect of TEE for an important population (elderly
Americans enrolled in Medicare) using STS ACSD, the world’s premier clinical
outcomes registry for adult cardiac surgery, and facilitate a more informed
comparison between results derived from different databases.

###  1.2 A Naive Method, An Existing Method, and A New Approach

One naive strategy to create homogeneous matched pairs resembling a template
would consist of two steps. Take the RCT DUPLICATE initiative as an example.
In the first step, researchers could select a subset of treated participants
from the observational data via matching on covariates collected by both the
RCT and the observational study. In the second step, controls in the
observational database are then matched to the treated participants selected
by the first step. Each step involves only two groups and could be done by
using standard statistical matching algorithms, for instance, network-flow-
based algorithms built upon bipartite networks, like Silber et al.,’s (2014)
template matching method, or mixed-integer-programming-based (MIP-based)
approaches (Zubizarreta,, 2012; Bennett et al.,, 2020).

$X_{1}$$X_{2}$ Template $X_{1}$$X_{2}$$X_{3}$T1T2 OBS Treated
$X_{1}$$X_{2}$$X_{3}$C1C2C3C4 OBS Control Figure 1: A toy example illustrating
the limitation of a naive two-step approach. Suppose our goal is to create one
treated-to-control pair. In the first step, treated units are matched to the
template to obtain a candidate treated group. In this toy example, both T1 and
T2 are good match to the template judging from two RCT covariates $X_{1}$ and
$X_{2}$. However, in view of the reservoir of control units, T1 is a better
pick as no control unit has an $X_{3}$ value similar to that of T2.

This two-step approach suffers from a major drawback. Franklin et al., (2021)
matched on more than $120$ covariates collected by the claims data to guard
against unmeasured confounding in their emulation of the target SAVOR-TIMI 53
trial, although the trial reported only around $20$ baseline covariates. In
the first step, there are many ways to design a smaller treated group similar
to the RCT-eligible population in $20$ RCT covariates; however, it is
difficult to determine which treated group designed from the first step should
be used to form the final match in the second step. A selected treated group
similar to the RCT-eligible population in $20$ RCT covariates could have poor
overlap with controls in the observational database in the other more than
$100$ covariates; see Figure 1 for a toy example illustrating this phenomenon.
In many applications, resemblance of the matched cohort data to some target
population is useful but not necessary; priority should be given to creating
well-matched treated and control groups to first maximize a study’s validity.
The naive two-step approach precludes a principled trade-off between a matched
cohort study’s internal validity and its generalizability to a target
population.

In principle, Bennett et al.,’s (2020) MIP-based method could be adapted to
designing a treated group and a matched control group that are well balanced
in many covariates while resembling a template in key covariates by modifying
constraints in their mathematical program. Although recent advancements in
computing power have made it more practical to solve large-scale MIP problems
with many complex constraints, the general MIP problems are still
theoretically intractable or NP-hard. From a very practical perspective, MIP-
based methods require installing a powerful commercial optimization routine
(e.g., Gurobi or IBM CPLEX), which is proprietary and could be an obstacle to
researchers. Other approaches to adjusting for sample selection bias include
weighting and doubly-robust methods (see, e.g., Stuart et al.,, 2011; Dahabreh
et al.,, 2019).

Compared to MIP-based statistical matching methods, network-flow-based methods
(see Supplemental Material A for a literature review) only require solving a
polynomial solvable problem that is tractable both in theory and practice, and
have proven successful in empirical comparative effectiveness research for
decades (Rosenbaum,, 2002, 2010; Stuart,, 2010; Austin,, 2011; Rassen et al.,,
2012). The primary goal of this article is to outline an efficient, network-
flow-based algorithm that designs matched pairs from observational data with
close resemblance to a target population. We demonstrate the usefulness of the
proposed method by revisiting the Women’s Health Initiative (WHI) study and
exploring how our method facilitates different study designs and yields
insight into the inconsistency between the WHI observational study and trial
findings. We found that the discrepancy regarding a cardioprotective effect of
the hormone replacement therapy persisted in a design that adjusted for the
cardiovascular risk profile differences between the observational and trial
data, but seemed to disappear in a design that further adjusted for the HRT
initiation age and previous estrogen-plus-progestin use, resonating with
similar findings in the seminal work by Hernán et al., (2008).

We conclude this section with an important caveat. Many reasons may contribute
to the disagreement between effect estimates derived from an observational
study and an RCT or from two observational studies. The synergy between effect
heterogeneity and difference in covariate distributions is only one of the
many possible reasons. Other important reasons include unmeasured confounding,
noncompliance in the RCT, difference in the definition of treatment and/or
clinical endpoints, among others. A matched cohort constructed from
observational data, no matter how balanced its treated and matched control
groups are, and how similar these two groups are to the target population, is
_not_ necessarily free from unmeasured confounding bias and in our opinion,
cannot replace an RCT. Nevertheless, addressing one critical and often
conspicuous issue, that is, difference in covariate distributions among
different studies, helps researchers focus on other issues and is a meaningful
step towards reconciling different study outcomes and understanding the
underlying mechanisms.

##  2 Methodology

###  2.1 Basic Network Structure: Vertices and Edges

We describe our proposed method using the RCT emulation study described in
Section 1.1 for illustration purpose; the method can accommodate an arbitrary
template from a target population other than the RCT-eligible population
(e.g., the population of Medicare beneficiaries as discussed in the TEE and
cardiac surgery example in Section 1.1).

Figure 2 depicts a basic version of the proposed network structure. There are
$R$ treated units from the target RCT. These $R$ RCT units
$\mathcal{K}=\\{\kappa_{1},\kappa_{2},\dots,\kappa_{R}\\}$ are represented by
nodes labeled $\kappa_{r}$, $r=1,2,\dots,R$. There are $T\geq R$ OBS treated
units and $C\geq T$ OBS control units from some administrative database. OBS
treated units $\mathcal{T}=\\{\tau_{1},\tau_{2},\dots,\tau_{T}\\}$ are
represented twice in the network, by nodes labeled $\tau_{t}$ and
$\overline{\tau}_{t}$, $t=1,2,\dots,T$, and OBS control units
$\mathcal{C}=\\{\gamma_{1},\gamma_{2},\dots,\gamma_{C}\\}$ are represented by
nodes labeled $\gamma_{c}$, $c=1,2,\dots,C$. In addition to $R+2\times T+C$
nodes representing RCT and OBS study units, there is a source node $\xi$ and a
sink node $\overline{\xi}$, so that the network consists of
$|\mathcal{V}|=R+2T+C+2=O(C)$ nodes in total:

$\mathcal{V}=\left\\{\xi,\kappa_{1},\dots,\kappa_{R},\tau_{1},\dots,\tau_{T},\overline{\tau}_{1},\dots,\overline{\tau}_{T},\gamma_{1},\dots,\gamma_{C},\overline{\xi}\right\\}.$
(1)

An ordered pair of vertices is referred to as an edge in the network. The
basic structure in Figure 2 consists of the following edges:

$\mathcal{E}=\left\\{(\xi,\kappa_{r}),(\kappa_{r},\tau_{t}),(\tau_{t},\overline{\tau}_{t}),(\overline{\tau}_{t},\gamma_{c}),(\gamma_{c},\overline{\xi}),~{}r=1,\dots,R,~{}t=1,\dots,T,~{}c=1,\dots,C\right\\}.$
(2)

There are a total of $|\mathcal{E}|=R+R\times T+T+T\times C+C=O(C^{2})$ edges
assuming $C$ is a constant multiple of $T$.

$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\tau_{1}$$\tau_{2}$$\tau_{3}$$\tau_{4}$$\overline{\tau}_{1}$$\overline{\tau}_{2}$$\overline{\tau}_{3}$$\overline{\tau}_{4}$$\gamma_{1}$$\gamma_{2}$$\gamma_{3}$$\gamma_{4}$$\gamma_{5}$$\gamma_{6}$$\xi$$\overline{\xi}$
RCT Treated OBS Treated OBS Treated OBS Control Figure 2: A dense network for
pair matching while emulating a target RCT population. Treated units from the
target RCT population appear on the far left as node $\kappa_{s}$,
$s=1,2,\dots,R$. Treated units from the observational database appear twice as
$\tau_{t}$ and $\overline{\tau}_{t}$, $t=1,2,\dots,T$. Control units from the
observational database appear on the far right as $\gamma_{c}$,
$c=1,2,\dots,C$. There is a source $\xi$ and a sink $\overline{\xi}$. We have
$R=3$, $T=4$, and $C=6$ in this simple example.

###  2.2 Basic Network Structure: Capacity and Matched Samples

To transform a statistical matching problem into an appropriate network-flow
optimization problem, one needs to carefully design the cost and capacity of
each edge. Fix an integer $k$ such that $1\leq k\leq\lfloor T/R\rfloor$, and
consider constructing $N=kR$ treated-to-control matched pairs from the
observational database. For instance, the SAVOR-TIMI 53 trial consists of
$R=8,280$ units assigned saxaglipitin, the intervention under evaluation, and
three retrospective databases available for emulation (Optum, MarketScan, and
Medicare) consist of $T=91,082$ subjects exposed to saxaglipitin according to
Franklin et al., (2021)’s study protocol registered at ClinicalTrials.gov
(identifier NCT03936023). Researcher could in principle choose any integer $k$
between $1$ and $\lfloor 91,082/8,280\rfloor=11$ in this example.

Let $\text{cap}(e)\geq 0$ denote the capacity of an edge $e\in\mathcal{E}$. In
the basic network structure depicted in Figure 2,
$\text{cap}\\{(\xi,\kappa_{r})\\}=k$ for $r=1,\dots,R$, and all other edges
have capacity $1$. In Figure 2, the source $\xi$ supplies $kR$ units of flow,
the sink $\overline{\xi}$ absorbs $kR$ units of flow, while all other nodes
preserve the flow by simply passing them along (Ahuja et al.,, 1988;
Bertsekas,, 1991; Rosenbaum,, 1989). A feasible flow $f(\cdot)$ of the
proposed network is formally defined as a mapping from the set of edges
$\mathcal{E}$ to $\\{0,1,2,\dots,k\\}$ such that (i) all capacity constraints
are respected, i.e., $0\leq f(e)\leq\text{cap}(e),~{}e\in\mathcal{E}$, (ii)
$kR$ units of flow are supplied at $\xi$ and absorbed at $\overline{\xi}$,
i.e., $\sum_{r=1}^{R}f\\{(\xi,\kappa_{r})\\}=kR$ and
$\sum_{c=1}^{C}f\\{(\gamma_{c},\overline{\xi})\\}=kR$, and (iii) the flow is
preserved at all nodes other than $\xi$ and $\overline{\xi}$, i.e.,
$\sum_{(a,b)\in\mathcal{E}\backslash\\{\xi,\overline{\xi}\\}}f\\{(a,b)\\}=\sum_{(b,c)\in\mathcal{E}\backslash\\{\xi,\overline{\xi}\\}}f\\{(b,c)\\}$
for all $b\in\mathcal{V}\backslash\\{\xi,\overline{\xi}\\}$.

It is beneficial to consider a concrete example. The toy example in Figure 2
has $R=3$, $T=4$, and $C=6$. Consider setting $k=1$ so that all edges in the
network have capacity $1$. Thick, black lines in Figure 3 correspond to one
(out of ${4\choose 3}\times 6\times 5\times 4=480$) feasible flows in this
network. The left part of the network helps select OBS treated units using RCT
units as a template, while the right part of the network performs the actual
statistical matching and outputs matched samples. Formally, the matched
samples $\mathcal{M}$ is defined by

$\mathcal{M}=\left\\{(\tau_{t},\gamma_{c})~{}\text{such
that}~{}f\\{(\tau_{t},\overline{\tau}_{t})\\}=f\\{(\overline{\tau}_{t},\gamma_{c})\\}=1\right\\}.$
(3)

For instance, the matched samples returned by the feasible flow in Figure 3
consist of
$\mathcal{M}=\\{(\tau_{1},\gamma_{1}),(\tau_{2},\gamma_{5}),(\tau_{3},\gamma_{4})\\}$.

$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\tau_{1}$$\tau_{2}$$\tau_{3}$$\tau_{4}$$\overline{\tau}_{1}$$\overline{\tau}_{2}$$\overline{\tau}_{3}$$\overline{\tau}_{4}$$\gamma_{1}$$\gamma_{2}$$\gamma_{3}$$\gamma_{4}$$\gamma_{5}$$\gamma_{6}$$\xi$$\overline{\xi}$
RCT Treated OBS Treated OBS Treated OBS Control Right Left Figure 3: A
feasible flow
$f(\cdot)\mathrel{\mathop{\ordinarycolon}}\mathcal{E}\mapsto\\{0,1\\}$ in a
dense network. Thick, black lines correspond to edges with $f(e)=1$ while
light gray lines correspond to edges with $f(e)=0$. The feasible flow yields a
matched sample of three pairs:
$\mathcal{M}=\\{(\tau_{1},\gamma_{1}),(\tau_{2},\gamma_{5}),(\tau_{3},\gamma_{4})\\}$.

###  2.3 Basic Network Structure: Cost, Probability of Participation,
Propensity Score

While network infrastructure, i.e., vertices, edges, and capacities,
determines the collection of all feasible flows, costs associated with each
edge help select one best suited for empirical researchers’ specific purposes.
Let $\text{cost}(e)$ denote the cost associated with edge $e\in\mathcal{E}$.
In the basic network structure depicted in Figure 2, we let
$\text{cost}\\{(\xi,\kappa_{r})\\}=\text{cost}\\{(\tau_{t},\overline{\tau}_{t})\\}=\text{cost}\\{\gamma_{c},\overline{\xi}\\}=0$
for $r=1,\dots,R$, $t=1,\dots,T$, and $c=1,\dots,C$.

Costs associated with edges $(\kappa_{r},\tau_{t})$ play an important role in
forcing selected OBS samples to mimic RCT units in some covariates. Suppose
that each RCT unit $\kappa_{r}$ is associated with a vector of covariates
$\widetilde{\boldsymbol{x}}$ and each OBS unit
$(\widetilde{\boldsymbol{x}},\boldsymbol{x})$. As in the RCT DUPLICATE
example, $\widetilde{\boldsymbol{x}}$ contains roughly $20$ covariates that
both RCT and observational database collected while $\boldsymbol{x}$ consist
of more than $100$ additional covariates available only in the claims
database. Let $\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})$
denote a measure of covariate distance between $\kappa_{r}$ and $\tau_{t}$ in
$\widetilde{\boldsymbol{x}}$. The cost associated with an edge of the form
$(\kappa_{r},\tau_{t})$ is equal to
$\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})$, i.e.,

$\text{cost}\\{(\kappa_{r},\tau_{t})\\}=\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}}).$
(4)

One intriguing strategy equals
$\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})$ to a scalar
“balancing score” of $\widetilde{\boldsymbol{x}}$ (Rosenbaum and Rubin,,
1983), so that matching on this balancing score stochastically balances
$\widetilde{\boldsymbol{x}}$. In the context of generalizing RCT’s effect
estimates to a target population, Stuart et al., (2011) defined and studied
the conditional probability of selecting into the RCT group rather than the
OBS group, which is referred to as the “probability of participation.” In
practice, researchers could collect the RCT treated units $\mathcal{K}$ and
OBS treated units $\mathcal{T}$, and estimate the conditional probability of
being selected into the RCT given covariates $\widetilde{\boldsymbol{x}}$
using a logistic regression (Stuart et al.,, 2011).

Lastly, let
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
denote a measure of distance between OBS treated unit $\tau_{t}$ and OBS
control unit $\gamma_{c}$ in their observed covariates
$(\widetilde{\boldsymbol{x}},\boldsymbol{x})$. The cost associated with an
edge $(\overline{\tau}_{t},\gamma_{c})$ is set to
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$,
i.e.,

$\text{cost}\\{(\overline{\tau}_{t},\gamma_{c})\\}=\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}.$
(5)

Different specifications of
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
have been extensively studied in the literature (Rosenbaum,, 2002, 2010;
Stuart,, 2010). Perhaps the most widely used strategy is to equal
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
to the absolute difference in Rosenbaum and Rubin, (1983)’s propensity score.
Alternatively, one may let
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
represent the Mahalanobis distance in
$(\widetilde{\boldsymbol{x}},\boldsymbol{x})$ (Cochran and Rubin,, 1973;
Rubin,, 1979), or the Mahalanobis distance within a propensity score caliper
(Rosenbaum and Rubin,, 1985).

###  2.4 Minimum Cost Flow, Complexity, Trade-Off Between Internal and
External Validity

The cost of a feasible flow $f(\cdot)$ in the proposed network is equal to

$\text{cost}(f)=\sum_{e\in\mathcal{E}}f(e)\cdot\text{cost}(e)=\underbrace{\sum_{\begin{subarray}{c}(\kappa_{r},\tau_{t})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\kappa_{r},\tau_{t})\\}=1\end{subarray}}\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})}_{S1}~{}+\underbrace{\sum_{\begin{subarray}{c}(\overline{\tau}_{t},\gamma_{c})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\overline{\tau}_{t},\gamma_{c})\\}=1\end{subarray}}\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}}_{S2},$
(6)

and a flow $f(\cdot)$ is a minimum cost flow if it is feasible and every other
feasible flow has a cost at least as high as $\text{cost}(f)$. A minimum cost
flow can be found in
$O(|\mathcal{V}|\cdot|\mathcal{E}|+|\mathcal{V}|^{2}\log(|\mathcal{V}|))$
operations (Korte and Vygen,, 2011) and in the proposed network, we have
$|\mathcal{V}|=O(C)$, $|\mathcal{E}|=O(C^{2})$, so that the computation
complexity
$O(|\mathcal{V}|\cdot|\mathcal{E}|+|\mathcal{V}|^{2}\log(|\mathcal{V}|))$
simplifies to $O(C^{3})$. We discuss ways to sparsify the network and speed up
computation in Supplemental Material B. In the statistical computing software
R, the minimum cost flow can be found via Bertsekas, (1991)’s auction
algorithm implemented by Bertsekas and Tseng, (1988) and made available by
Hansen, (2007) and Pimentel et al., (2015).

There is a tension between internal and external validity in a matched
observational study. Take the RCT DUPLICATE initiative as an example. In their
emulation of the SAVOR-TIMI 53 trial using observational data, Franklin et
al., (2021) matched closely on and balanced almost $120$ observed covariates;
however, the effect estimate obtained from observational data (HR: $0.81$;
$95\%$ CI: $(0.76,0.86)$) was significantly different from the RCT effect
estimate (HR: $1.00$; $95\%$ CI: $(0.89,1.12)$). There are many possible
explanations in this large discrepancy, one of which is a quite large
difference in baseline characteristics of the RCT versus OBS populations,
including a disparity in some important preexisting conditions like MI. It is
one possibility that both the observational study and RCT effect estimates are
internally valid; they are just reporting effects for different populations.

According to our formulation, the cost
$S1=\sum_{\begin{subarray}{c}(\overline{\tau}_{t},\gamma_{c})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\overline{\tau}_{t},\gamma_{c})\\}=1\end{subarray}}\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
in expression (6) measures the homogeneity between the matched treated and
control groups. Well-matched samples have a small $S1$ value and better
internal validity. On the other hand, the cost
$S2=\sum_{\begin{subarray}{c}(\kappa_{r},\tau_{t})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\kappa_{r},\tau_{t})\\}=1\end{subarray}}\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})$
measures how much the matched samples mimic a target population, and a small
$S2$ corresponds to improved generalizability to the target population.

To facilitate exploring the trade-off between a matched cohort study’s
internal validity and its generalizability to a target population, we replace
$\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
with
$\lambda\cdot\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}$
for some $\lambda>0$ so that $\text{cost}(f)$ becomes:

$\text{cost}(f)=\sum_{e\in\mathcal{E}}f(e)\cdot\text{cost}(e)=\sum_{\begin{subarray}{c}(\kappa_{r},\tau_{t})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\kappa_{r},\tau_{t})\\}=1\end{subarray}}\delta_{\kappa_{r},\tau_{t}}(\widetilde{\boldsymbol{x}})~{}+\lambda\cdot\bigg{\\{}\sum_{\begin{subarray}{c}(\overline{\tau}_{t},\gamma_{c})\in\mathcal{E}\mathrel{\mathop{\ordinarycolon}}\\\
f\\{(\overline{\tau}_{t},\gamma_{c})\\}=1\end{subarray}}\Delta_{\tau_{t},\gamma_{c}}\\{(\widetilde{\boldsymbol{x}},\boldsymbol{x})\\}\bigg{\\}}.$
(7)

According to this formulation, a large $\lambda$ value gives priority to a
matched cohort study’s internal validity, while a small $\lambda$ value
prioritizes its generalizability to the target population. A similar weighting
scheme is also used in Zhang et al., (2021) but for a different purpose.

In Web Appendix B, we further discuss how to speed up computation and
incorporate additional design features, including exact and near-exact
matching on effect modifiers, Rosenbaum,’s fine and near-fine balance (1989)
and how to force including or excluding certain treated units, by modifying
various aspects of the proposed network structure.

###  2.5 Practical considerations and software availability

The proposed matching algorithm is integrated in the package match2C available
via the statistical computing software R (R Core Team,, 2021) with a detailed
tutorial.

To successfully implement the method, users need to specify two parameters.
Parameter $k$ controls the size of the matched cohort and is determined by the
size of the template, size of the observational database, overlap of the
template and observational database in covariate distributions and power
considerations. In general, a small $k$ value corresponds to constructing a
small matched cohort and this small matched cohort by design will be more
closely matched and better resemble the template compared to a larger matched
cohort corresponding to a larger $k$ value; see the small illustrative example
in Web Appendix C.2 and simulation studies in Section 3. If researchers find
the matched comparison not adequately powered, then larger $k$ values should
be explored. Parameter $\lambda$ controls the trade-off between the validity
of the matched comparison and its generalizability to the target population.
By default, the parameter $\lambda$ is set to a large number to prioritize the
validity of the matched cohort study. When the treated and control groups in
the observational database are well-overlapped so that there are many possible
internally-valid matched cohort studies, researchers may then reduce $\lambda$
to further improve the matched cohort study’s generalizability.

Matching is part of the design of an observational study and should be carried
out without looking at the outcome data. Good practice includes keeping time-
stamped analysis logs for review and posting a detailed pre-analysis protocol;
see, e.g., Franklin et al., (2020, 2021). Provided that no outcome data are
viewed, researchers typically perform statistical matching multiple times and
select the design based on covariate balance. Recently, many formal diagnostic
tests of covariate balance have been proposed (Gagnon-Bartsch and Shem-Tov,,
2019). Researchers could perform a formal diagnostic test, e.g., Gagnon-
Bartsch and Shem-Tov,’s (2019) classification permutation test (CPT) to see if
there is any residual imbalance in observed covariates
$(\widetilde{\boldsymbol{x}},\boldsymbol{x})$ between two groups in the
matched cohort; analogously, a formal test could be carried out to examine
residual imbalance in common covariates $\widetilde{\boldsymbol{x}}$ between
the matched group and the template. We will discuss these aspects more
concretely when examining the WHI study in Section 4.

##  3 Simulation study

###  3.1 Goal; structure; measurement of success

Our primary goal in this section is to examine how the study design delivered
by different matching algorithms affects the bias. In particular, we are
interested in the case where the treatment effect is heterogeneous and the
distributions of effect modifiers vary in the observational database and the
target population. We considered the following target population:
$\mathbf{X}\sim\text{Multivariate
Normal}\left(\mathbf{\mu},\mathbf{\Sigma}\right)$ with
$\mathbf{\mu}=(0.25,0,0,0,0)^{\text{T}}$ and
$\mathbf{\Sigma}=\mathbf{I}_{5\times 5}$, and generated a template
$\mathcal{K}$ consisting of a random sample of size $300$ from the target
population. We consider an observational database with $|\mathcal{T}|=1,000$
treated and $|\mathcal{C}|=3,000$ control units. The data-generating process
for units in the observational database and statistical matching procedures to
be investigated are specified via the following factorial design:

Factor 1: Dimension of covariates in the observational database, $d$: $10$,
$30$, and $50$.

Factor 2: Overlap, $\theta$: $\mathbf{X}\sim\text{Multivariate
Normal}\left(\mathbf{\mu},\mathbf{\Sigma}\right)$, with $\mathbf{\mu}=(\theta
Z,0,\dots,0)^{\text{T}}$ and $\mathbf{\Sigma}=\mathbf{I}_{d\times d}$. We
consider $\theta=0.50$ and $0.75$.

Factor 3: Matching algorithms to be investigated, $\mathcal{M}$:

  1. 1.

$\mathcal{M}_{\textsf{opt}}$: matching according to two criteria (Zhang et
al.,, 2021): (i) minimizing the earthmover’s distance, a measure of distance
between two probability distributions (Levina and Bickel,, 2001), between the
distributions of the estimated propensity scores in the treated and matched
control groups, and (ii) minimizing the within-matched-pair robust Mahalanobis
distances. Algorithm $\mathcal{M}_{\textsf{opt}}$ produces $1,000$ matched
pairs and does not make use of the template.

  2. 2.

$\mathcal{M}_{\textsf{template},~{}k=1}$: matching according to the proposed
network structure (Figure 3) with $k=1$. Algorithm
$\mathcal{M}_{\textsf{template},~{}k=1}$ produces $300$ matched pairs.

  3. 3.

$\mathcal{M}_{\textsf{template},~{}k=2}$: similar to
$\mathcal{M}_{\textsf{template},~{}k=1}$ but with $k=2$. Algorithm
$\mathcal{M}_{\textsf{template},~{}k=2}$ produces $600$ matched pairs.

Factor 4: Tuning parameter in $\mathcal{M}_{\textsf{template},~{}k}$,
$\lambda$: $100$, $1$, $0.01$.

Factor $1$ and $2$ define the data-generating process for units in the
observational database. Factor $3$ and $4$ define a total of $1+2\times 3=7$
matching algorithms to be investigated. As discussed in Section 2.4, the
tuning parameter $\lambda$ controls the trade-off between a matched cohort
study’s validity and its generalizability to the target population.

For each unit, we further generate two potential outcomes:

$Y(0)\sim\text{Normal}(\mathbf{X}^{T}\boldsymbol{\nu},1),\qquad
Y(1)=Y(0)+\beta(X_{1}),$ (8)

and the observed outcome satisfies $Y=Z\cdot Y(1)+(1-Z)\cdot Y(0)$. Factor $5$
specifies the mean of the potential outcome $Y(0)$ and Factor $6$ specifies
the treatment effect:

Factor 5: Mean of $Y(0)$, $\mathbf{X}^{T}\boldsymbol{\nu}$: a constant
$\boldsymbol{\nu}$ vector with all entries equal to $0$, $0.05$ or $0.1$.

Factor 6: Treatment effect, $\beta(X_{1})$: a constant effect
$\beta(X_{1})=2$, a mildly heterogeneous effect $\beta(X_{1})=2-0.2X_{1}$, and
a strongly heterogeneous effect $\beta(X_{1})=2-X_{1}$.

The average treatment effect of the target population satisfies
$\textsf{ATE}_{\text{target}}=2$ when $\beta(X_{1})=2$,
$\textsf{ATE}_{\text{target}}=1.95$ when $\beta(X_{1})=2-0.2X_{1}$ and
$\textsf{ATE}_{\text{target}}=1.75$ when $\beta(X_{1})=2-X_{1}$.

There are multiple ways to analyze matched-pair data. Examples include a
parametric t-test, randomization inference (Rosenbaum,, 2002, 2010) and
regression adjustment (Rubin,, 1979; Ho et al.,, 2007). In this simulation
study, we calculated a difference-in-means estimator for matched data produced
by each of the $7$ algorithms and compared these $7$ effect estimates to the
target parameter $\textsf{ATE}_{\text{target}}$ in data-generating process.

###  3.2 Simulation results

Table 1 summarizes the percentage of bias of each difference-in-means
estimator when the overlap parameter $\theta=0.50$. Web Appendix D.1 reports
similar results for $\theta=0.75$.

We have observed a few consistent trends. First, when the treatment effect is
constant, the bias is relatively small under all data-generating processes and
statistical matching algorithms under consideration. Second, when the
treatment effect is heterogeneous and the effect modifier $X_{1}$ has a
different distribution in the template and the group of treated units, the
effect estimate obtained from $1000$ matched pairs constructed using algorithm
$\mathcal{M}_{\textsf{opt}}$ is clearly biased from
$\textsf{ATE}_{\text{target}}$, and the percentage of bias increases (i) as
the distributions of $X_{1}$ in the target population and in the group of OBS
treated units become increasingly dissimilar, that is, as $\theta$ increases,
and (ii) as effect modification becomes more pronounced, that is, from
$\beta(X_{1})=2$ (constant) to $\beta(X_{1})=2-0.2X_{1}$ (mild) to
$\beta(X_{1})=2-X_{1}$ (strong). In the most adversarial scenario with
$\beta(X_{1})=2-X_{1}$, the percentage bias of
$\widehat{\theta}_{\mathcal{M}_{\textsf{opt}}}$ can be as large as $25\%$. We
need to stress that although $\widehat{\theta}_{\mathcal{M}_{\textsf{opt}}}$
may not be generalized to the target population, it is an internally-valid
estimator for the average treatment effect on the treated.

Our proposed algorithm outperforms $\mathcal{M}_{\textsf{opt}}$ in bias
reduction against $\textsf{ATE}_{\text{target}}$ in all $6$ different
implementations under all data-generating processes considered in this
simulation study, although the gain in bias reduction differs from
implementation to implementation. In particular, we observe that the gain is
most pronounced when (i) $k$ is small so that a smaller treated group bearing
more resemblance to the target template is constructed, and (ii) $\lambda$ is
small so that the matching algorithm gives priority to resemblance to the
template.

Table 1: Percentage of bias with respect to $\textsf{ATE}_{\text{target}}$ of
$7$ difference-in-means estimator constructed from matched samples obtained
from each of the $7$ matching algorithms under consideration. The overlap
parameter $\theta$ in Factor 2 is equal to $0.50$. Each cell is averaged over
$1000$ simulations.

Heterogeneity Level $\beta(X_{1})$ |  DGP of $Y(0)$ $\nu$ |  |  $\mathcal{M}_{\textsf{opt}}$ |  | $\mathcal{M}_{\textsf{template}}$
---|---|---|---|---|---
|  |  |  |  |  $k=1$ $\lambda=0.01$ |  $k=1$ $\lambda=1$ |  $k=1$ $\lambda=100$ |  $k=2$ $\lambda=0.01$ |  $k=2$ $\lambda=1$ |  $k=2$ $\lambda=100$
$d=10$
Constant | 0 |  | 0.07% |  | 0.21% | -0.06% | -0.06% | -0.01% | 0.05% | -0.06%
0.05 |  | 0.00% |  | 0.33% | 0.18% | 0.14% | 0.23% | 0.25% | 0.44%
0.10 |  | 0.06% |  | 0.59% | 0.30% | 0.44% | 0.73% | 0.66% | 0.61%
Mild | 0 |  | -5.17% |  | -1.23% | -1.19% | -2.04% | -1.52% | -1.56% | -2.56%
0.05 |  | -4.80% |  | -0.11% | -0.15% | -1.25% | -0.53% | -0.65% | -1.77%
0.10 |  | -5.00% |  | 0.04% | -0.20% | -1.22% | -0.58% | -0.48% | -1.44%
Strong | 0 |  | -24.99% |  | -4.01% | -3.75% | -9.05% | -6.45% | -5.96% | -11.68%
0.05 |  | -25.18% |  | -3.10% | -2.92% | -8.40% | -5.50% | -5.60% | -11.43%
0.10 |  | -25.23% |  | -3.68% | -3.39% | -8.07% | -6.14% | -5.91% | -11.25%
$d=30$
Constant | 0 |  | 0.23% |  | 0.09% | 0.45% | 0.40% | 0.31% | 0.33% | 0.28%
0.05 |  | 0.13% |  | 1.04% | 1.27% | 1.37% | 1.13% | 1.19% | 1.29%
0.10 |  | 0.12% |  | 1.55% | 1.43% | 1.62% | 1.70% | 1.61% | 1.83%
Mild | 0 |  | -4.77% |  | -0.57% | -1.32% | -2.06% | -1.07% | -1.76% | -2.50%
0.05 |  | -5.07% |  | -0.55% | -0.83% | -1.56% | -0.67% | -1.27% | -1.88%
0.10 |  | -4.63% |  | 0.99% | 0.80% | 0.48% | 0.62% | 0.06% | -0.41%
Strong | 0 |  | -24.71% |  | -3.57% | -5.90% | -9.28% | -5.65% | -9.50% | -12.83%
0.05 |  | -24.89% |  | -3.13% | -5.44% | -8.80% | -5.30% | -8.96% | -12.37%
0.10 |  | -24.78% |  | -2.49% | -4.88% | -8.01% | -4.51% | -8.21% | -11.35%
$d=50$
Constant | 0 |  | 0.06% |  | 0.31% | -0.21% | -0.21% | 0.19% | -0.12% | -0.06%
0.05 |  | 0.07% |  | 0.77% | 1.15% | 1.12% | 0.99% | 1.29% | 1.29%
0.10 |  | -0.05% |  | 1.33% | 1.82% | 2.07% | 1.66% | 2.02% | 2.18%
Mild | 0 |  | -5.12% |  | -1.02% | -1.59% | -2.20% | -1.60% | -2.43% | -2.94%
0.05 |  | -5.15% |  | 0.37% | -0.69% | -1.26% | -0.37% | -1.16% | -1.63%
0.10 |  | -5.03% |  | 0.78% | 0.20% | -0.24% | 0.47% | 0.01% | -0.40%
Strong | 0 |  | -24.64% |  | -3.48% | -6.82% | -9.48% | -5.84% | -10.26% | -12.72%
0.05 |  | -25.00% |  | -3.22% | -6.76% | -9.06% | -5.39% | -9.88% | -12.21%
0.10 |  | -24.78% |  | -1.80% | -5.01% | -7.37% | -3.97% | -8.34% | -10.63%

###  3.3 Additional simulation studies

In Web Appendix D.2, we designed a simulation study to investigate the role of
unmeasured confounding. In the presence of unmeasured confounding in the
observational database, the bias consisted of two parts, a generalization bias
and an unmeasured confounding bias. By removing most of the generalization
bias using our proposed algorithm, researchers could better focus on the
unmeasured confounding bias. In Web Appendix D.3, we further compared the
computation cost of a network-flow-based algorithm and an MIP-based algorithm
implemented in the R package designmatch based on an open source optimization
routine GLPK. We found that the network-flow-based method largely outperforms
the MIP-based method in this comparison, although we note that the performance
of MIP-based methods would improve when they are implemented using a more
powerful commercial optimization routine like Gurobi.

##  4 Revisiting the Women’s Health Initiative (WHI) study

###  4.1 Background and our goal

The Women’s Health Initiative (WHI) is a combined clinical trial and
observational study. Postmenopausal women were screened for clinical trial
eligibility; those who were ineligible or unwilling to participate in the
trial were enrolled in the observational study. The design of the WHI study is
described in The Women’s Health Initiative Study Group, (1998). One important
goal of the WHI clinical trial is to evaluate the hypothesized
cardioprotective effect of postmenopausal hormone therapy, following a
substantial body of evidence from observational studies (The Women’s Health
Initiative Study Group,, 1998). The WHI estrogen-plus-progestin (E + P for
short) trial found a rather surprising elevation in coronary heart disease
risk (Writing Group for the Women’s Health Initiative Investigators,, 2002)
and sparked a lot of discussion regarding the discrepancy between clinical
trial and observational study results.

Many authors have speculated on why results would differ dramatically between
the trial and observational study. Some major concerns include: (i) potential
bias due to unmeasured confounding in observational studies (Humphrey et al.,,
2002); (ii) biological differences between trial participants and those in the
observational study (Michels and Manson,, 2003) and (iii) differences in time
since menopause at hormone therapy initiation (Prentice et al.,, 2005; Willett
et al.,, 2006; Hernán et al.,, 2008; Prentice et al.,, 2009), among others.
Table 2 summarizes some important baseline covariates in the WHI trial and
associated observational study and illustrates some of these concerns.
Compared to past and never users, current users of estrogen-plus-progestin in
the WHI observational study are younger, less likely to be black or Hispanic,
more educated, and have more physical activity episodes per week. On the other
hand, baseline covariates of participants in the control and intervention
groups are similar in the WHI trial by virtue of randomization. There is also
a tangible difference in socioeconomic status and smoking status between trial
and observational study participants. Moreover, trial participants initiated
their HRT at a much older age compared to the current users in the
observational study.

Table 2: Important baseline characteristics of the WHI observational study and
WHI trial subjects. Mean (SE) are reported for continuous variables and count
(%) for categorical variables.

| WHI Observational Data |  | WHI Trial
---|---|---|---
|  Never/Past Users ($\mathbf{n=75303}$) |  Current Users ($\mathbf{n=18340}$) |  |  Control ($\mathbf{n=8102}$) |  E+P Intervention ($\mathbf{n=8506}$)
Age at screening | 64.30 (7.37) | 60.84 (6.69) |  | 63.33 (7.11) | 63.23 (7.13)
Race/Ethnicity |  |  |  |  |
White | 61704 (81.9) | 16285 (88.8) |  | 6805 (84.0) | 7141 (84.0)
Black/Hispanic | 10166 (13.5) | 1074 ( 5.9) |  | 989 (12.2) | 1019 (12.0)
Other | 3433 ( 4.6) | 981 ( 5.3) |  | 308 ( 3.8) | 346 ( 4.1)
Education |  |  |  |  |
College or above | 30418 (40.4) | 10044 (54.8) |  | 3011 (37.2) | 3111 (36.6)
Some college | 28074 (37.3) | 5853 (31.9) |  | 3060 (37.8) | 3357 (39.5)
High school diploma/GED | 13068 (17.4) | 2050 (11.2) |  | 1609 (19.9) | 1615 (19.0)
Other | 3743 ( 5.0) | 393 (2.1) |  | 422 (5.2) | 423 (5.0)
Blood Pressure |  |  |  |  |
Systolic | 128 (18.12) | 123 (16.92) |  | 128 (17.53) | 128 (17.63)
Diastolic | 75 (9.42) | 74 (8.99) |  | 76 (9.09) | 76 (9.12)
BMI | 27.62 (5.98) | 25.84 (5.15) |  | 28.50 (5.91) | 28.46 (5.82)
Smoking |  |  |  |  |
NA | 1004 ( 1.3) | 220 ( 1.2) |  | 98 ( 1.2) | 83 ( 1.0)
Current Smoker | 4850 ( 6.4) | 940 ( 5.1) |  | 838 (10.4) | 880 (10.3)
Never Smoked | 38296 (50.9) | 8710 (47.5) |  | 3999 (49.4) | 4178 (49.1)
Past Smoker | 31048 (41.3) | 8448 (46.1) |  | 3157 (39.0) | 3362 (39.5)
No. of PA episodes |  |  |  |  |
Total | 5.28 (4.12) | 5.80 (4.12) |  | 4.77 (4.06) | 4.74 (4.10)
Medium to strenuous | 2.94 (3.37) | 3.59 (3.57) |  | 2.58 (3.24) | 2.50 (3.21)
HRT Initiation |  |  |  |  |
NA | 67259 (89.3) | 0 |  | 6706 (82.8) | 0
Age at initiation | 52.04 (7.35) | 53.87 (6.84) |  | 53.57 (6.53) | 61.93 (8.15)
Previous E+P use in years | 0.47 (2.12) | 7.02 (5.64) |  | 0.64 (2.23) | 0.70 (2.34)
Unopposed estrogen use ever |  |  |  |  |
Yes | 32705 (43.4) | 2170 (11.8) |  | 865 (10.7) | 903 (10.6)
No | 42598 (56.6) | 16170 (88.2) |  | 7237 (89.3) | 7603 (89.4)
Reproductive history |  |  |  |  |
No ovary removed | 48303 (64.1) | 16906 (92.2) |  | 7705 (95.1) | 8083 (95.0)

Our proposed matching algorithm could help alleviate these concerns. The first
concern regarding the validity of the observational study could be alleviated
if matched observational study participants are balanced for a large number of
baseline covariates including detailed demographics, preexisting comorbid
conditions and personal habits. The second and third concerns about the
comparability between OBS and RCT participants in their cardiovascular risk
profile and HRT initiation time may be mitigated by treating the RCT cohort as
a template and constructing matched OBS pairs that resemble the RCT cohort in
these aspects.

Our proposed algorithm cannot fully address the concern regarding unmeasured
confounding in the observational study. Another concern that cannot be
adequately addressed is the difference in estrogen-plus-progestin usage time
prior to the WHI study between RCT and OBS samples: current users in the
observational study had used HRT for more than $7$ years on average, while
participants of the trial had used HRT for less than a year on average. This
difference is largely due to the design of the WHI study: the observational
study enrolled a cross-section of post-menopausal women that included many who
had been using the therapy for years (Willett et al.,, 2006) while most trial
participants were never-users of the therapy by the time of randomization.
This lack of overlap between RCT and OBS samples in previous estrogen-plus-
progestin usage makes it virtually impossible to balance this aspect between
RCT and OBS samples using any matching tool; nevertheless, in our opinion,
alleviating some major concerns helps researchers critically examine the
remaining explanations and is a meaningful step towards reconciling the WHI
RCT and OBS findings.

###  4.2 Template and three matched samples

Our goal is to construct well-matched pairs from the WHI observational study
that resemble the WHI trial participants in the following sense: (i) matched
pairs resemble the WHI trial participants in cardiovascular risk profile, and
(ii) treated participants in the matched pairs resemble the WHI trial
intervention group in their HRT initiation time. Our template consists of a
random sample of $1,000$ WHI trial participants in the intervention group with
the following covariates: risk factors listed in Table 2 plus the HRT
initiation time.

We applied the proposed matching algorithm to constructing matched pairs of
different sizes corresponding to choosing different $k$’s in Section 2.2. Our
desired matched OBS cohort would have similar cardiovascular risk profile as
the template and are closely matched for many other baseline covariates to
maximally guard against the unmeasured confounding bias due to the non-
randomized nature of the WHI observational study. In addition to risk factors
listed in Table 2, we further matched on participants’ region, partners’
education level, income level, marital status, reproductive history and eight
important preexisting conditions.

Table 3 summarizes results from two matched samples constructed using the
proposed algorithm with different parameters. Match M1 constructed $10,000$
matched pairs of two observational study participants: one treated and the
other control. Match M1 used a Mahalanobis distance with an estimated
generalizability score caliper (caliper size = $0.05$) and a Mahalanobis
distance with an estimated propensity score caliper (caliper size = $0.15$) on
the left and right part of the network depicted in Figure 3, respectively.
Because of the large sample size ($\approx 100,000$ in the observational
study), we applied a “hard” caliper in the sense that edges connecting one OBS
treated and one OBS control units are removed whenever they differ in their
estimated propensity score by more than the caliper size; in this way, the
network is sparsified and computation is boosted; see Web Appendix B for
details. We set $\lambda=100$ to give priority to internal validity of the
matched comparison. Match M2 is similar to M1 except that we only formed
$3,000$ matched pairs. Lastly, match M0 formed $18,340$ matched pairs
exhausting every treated observational study participant. The balance table
and propensity score distributions of M0 can be found in Web Appendix E.

Judging from internal validity, all three matched cohorts are acceptable: the
absolute standardized mean differences (SMDs) of most cardiovascular risk
factors and additional OBS covariates are less than $0.1$, or one-tenth of one
pooled standard deviation (Silber et al.,, 2001), and many are below $0.01$.
However, the three designs differ, sometimes significantly, in their
resemblance to the trial population. In particular, design M0 differs from the
trial participants in their cardiovascular risk profile (e.g., the percentage
of black/Hispanic is $6\%$ in M0 compared to $12\%$ in the template) and HRT
initiation age ($54$ in M0 compared to $62$ in the trial population). On the
other hand, both M1 and M2 are similar to the trial population in the
cardiovascular risk profile: the percentage of white (black/Hispanic) women
was $85\%$ ($10\%$) in M1, $84\%$ ($12\%$) in M2 and $84\%$ ($12\%$) in the
template, compared to $89\%$ ($6\%$) in M0; the percentage with a college
degree or above is $41\%$ in M1, $36\%$ in M2, and $38\%$ in the template,
compared to $55\%$ in M0. In fact, the matched treated group in either M1 or
M2 is now indistinguishable from the template based on the cardiovascular risk
profile as judged by Gagnon-Bartsch and Shem-Tov,’s classification permutation
test (2019). Moreover, the treated group of M1 now initiated their HRT at an
average age of $57$ with an average $5.54$-year of HRT usage prior to
enrolling in the WHI study, and the treated group of M2 initiated their HRT at
an average age of $61$ with an average of $3.24$-year previous HRT usage. The
HRT initiation age in the design M2 is now closer to that of the trial
population, though a potentially meaningful discrepancy in the previous HRT
usage persists ($3.24$ vs $0.70$ years).

###  4.3 Comparing survival outcomes

Figure 4 plots the Kaplan-Meier curves (Kaplan and Meier,, 1958) in the
treatment and control groups in each study design. Compared to all controls in
the unadjusted analysis (top left panel of Figure 4), matched controls in M0
have higher survival probability, reflecting the fact that these matched
controls are now more identical to the treated OBS participants and healthier;
however, the cardioprotective effect of HRT persists ($P$-value = $0.020$
according to O’Brien and Fleming,’s (1987) matched-pair Prentice-Wilcoxon
test) in M0 and is consistent with the previous analysis of the WHI
observational data (Prentice et al.,, 2005). Compared to their counterparts in
M0, treated and matched comparison groups in M1 have lower survival
probabilities; however, a qualitatively similar cardioprotective effect
persists in the design M1 ($P$-value = $0.003$). Recall that the design M1 has
largely mimicked the cardiovascular risk profile of the trial population;
therefore, it seems that the difference in cardiovascular risk profile between
observational study and trial populations is _not_ sufficient in explaining
the inconsistency in the trial and observational study results. Compared to M0
and M1, the matched design M2 best resembles the trial population in their HRT
initiation age and previous HRT use; when we examine the survival outcomes in
M2, the cardioproctective effect seemed to disappear ($P$-value = $0.495$).

Our results suggested that once we maximally restricted ourselves to the
observational study cohort resembling the RCT population, the comparison
between the treated and matched control groups seemed to start aligning
reasonably well with that in the RCT. Unfortunately, this finding still cannot
lend full credibility to the WHI observational study because of unmeasured
confounding heterogeneity: it is conceivable that little unmeasured
confounding is present in the 3,000 matched pairs in design M2, but sizable
unmeasured confounding could persist in the rest of the observational study
and still nullify the observational study result. With this important caveat
in mind, three comparisons facilitated by three different study designs M0, M1
and M2 seemed to support Willett et al.,’s (2006) assessment that HRT
initiation age and previous HRT usage played a key role in the discrepancy
between WHI observational study and trial findings, and resonate with similar
findings in Hernán et al., (2008).

(a) Unadjusted (b) Design M0

(c) Design M1 (d) Design M2

Figure 4: Survival outcomes for coronary heart disease (CHD) in the unadjusted
observational study samples (top left), design M0 (top right), design M1
(bottom left), and design M2 (bottom right). Table 3: Balance table of
baseline covariates in WHI observational study before matching, two matched
samples constructed using the proposed algorithm, and the matched template.
Match M1 constructed $10000$ matched pairs and Match M2 $3000$ matched pairs.

|  OBS Treated |  Matched Treated M1 |  Matched Treated M2 |  Template |  Matched Control M2 |  Matched Control M1 |  OBS Control |  SMD M1 |  SMD M2
---|---|---|---|---|---|---|---|---|---
Sample Size
N | 18340 | 10000 | 3000 | 1000 | 3000 | 10000 | 75303 |  |
HRT Prior Usage
Age at initiation |  |  |  |  |  |  |  |  |
NA | 0 | 0 | 0 | 0 | 0.88 | 0.86 | 0.89 |  |
Age | 53.87 | 56.83 | 60.67 | 61.93 | 53.14 | 52.44 | 52.04 |  |
Previous E+P use, yrs | 7.02 | 5.54 | 3.24 | 0.70 | 0.45 | 0.51 | 0.47 |  |
Covariates Collected in RCT and OBS
Age at screening | 60.84 | 62.21 | 63.57 | 63.10 | 63.35 | 62.15 | 64.30 | 0.01 | 0.03
Race/Ethnicity |  |  |  |  |  |  |  |  |
White | 0.89 | 0.85 | 0.84 | 0.84 | 0.84 | 0.85 | 0.82 | -0.00 | -0.01
Black/Hispanic | 0.06 | 0.10 | 0.12 | 0.12 | 0.12 | 0.10 | 0.14 | 0.00 | 0.01
Education |  |  |  |  |  |  |  |  |
College or above | 0.55 | 0.41 | 0.36 | 0.38 | 0.36 | 0.40 | 0.40 | 0.01 | 0.00
Some college | 0.32 | 0.41 | 0.43 | 0.40 | 0.43 | 0.41 | 0.37 | -0.00 | 0.01
High school diploma/GED | 0.11 | 0.15 | 0.17 | 0.17 | 0.17 | 0.16 | 0.17 | -0.01 | -0.02
Blood pressure |  |  |  |  |  |  |  |  |
Systolic | 123.27 | 125.26 | 126.69 | 126.99 | 126.28 | 125.11 | 127.86 | 0.01 | 0.02
Diastolic | 74.03 | 74.42 | 75.01 | 75.66 | 74.81 | 74.36 | 74.91 | 0.01 | 0.02
BMI | 25.85 | 27.01 | 28.22 | 28.46 | 27.61 | 26.78 | 27.61 | 0.04 | 0.11
Smoking |  |  |  |  |  |  |  |  |
Current smoker | 0.05 | 0.09 | 0.11 | 0.10 | 0.11 | 0.09 | 0.06 | 0.00 | 0.01
Never smoked | 0.47 | 0.47 | 0.49 | 0.51 | 0.47 | 0.47 | 0.51 | 0.00 | 0.04
Past smoker | 0.46 | 0.43 | 0.39 | 0.37 | 0.41 | 0.43 | 0.41 | -0.00 | -0.05
No. of PA episodes |  |  |  |  |  |  |  |  |
Total | 5.80 | 5.23 | 4.79 | 4.71 | 4.83 | 5.22 | 5.29 | 0.00 | -0.01
Medium to strenuous | 3.58 | 3.05 | 2.59 | 2.56 | 2.67 | 3.01 | 2.94 | 0.01 | -0.02
Unopposed estrogen use ever | 0.12 | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 | 0.43 | 0.00 | -0.00
No ovary removed | 0.92 | 0.95 | 0.95 | 0.95 | 0.95 | 0.95 | 0.64 | -0.00 | 0.00
Additional OBS Covariates
Region |  |  |  |  |  |  |  |  |
Midwest | 0.22 | 0.22 | 0.22 |  | 0.24 | 0.24 | 0.22 | -0.04 | -0.04
Northeast | 0.18 | 0.20 | 0.21 |  | 0.24 | 0.21 | 0.24 | -0.03 | -0.07
South | 0.25 | 0.25 | 0.25 |  | 0.24 | 0.24 | 0.26 | 0.01 | 0.02
Partner’s education |  |  |  |  |  |  |  |  |
College or above | 0.43 | 0.37 | 0.32 |  | 0.31 | 0.35 | 0.30 | 0.03 | 0.03
Some college | 0.16 | 0.18 | 0.18 |  | 0.18 | 0.18 | 0.17 | 0.00 | -0.01
High school diploma/GED | 0.07 | 0.08 | 0.09 |  | 0.09 | 0.09 | 0.09 | -0.02 | -0.01
Income |  |  |  |  |  |  |  |  |
Below 35K | 0.23 | 0.29 | 0.35 |  | 0.37 | 0.30 | 0.40 | -0.03 | -0.03
35K - 75K | 0.42 | 0.43 | 0.41 |  | 0.39 | 0.43 | 0.36 | 0.00 | 0.03
Above 75K | 0.29 | 0.22 | 0.18 |  | 0.16 | 0.20 | 0.16 | 0.04 | 0.04
Marital status |  |  |  |  |  |  |  |  |
Married | 0.69 | 0.67 | 0.62 |  | 0.62 | 0.66 | 0.60 | 0.02 | 0.02
Divorced/Widowed | 0.26 | 0.29 | 0.33 |  | 0.33 | 0.30 | 0.35 | -0.01 | 0.00
Employment status |  |  |  |  |  |  |  |  |
Yes | 0.45 | 0.40 | 0.36 |  | 0.35 | 0.40 | 0.32 | 0.00 | 0.01
No | 0.53 | 0.58 | 0.62 |  | 0.63 | 0.58 | 0.65 | -0.00 | -0.01
Reproductive history |  |  |  |  |  |  |  |  |
Oral contraceptive use ever | 0.53 | 0.48 | 0.43 |  | 0.42 | 0.47 | 0.37 | 0.00 | 0.01
OC duration in years | 5.59 | 5.55 | 5.47 |  | 5.18 | 5.32 | 5.16 | 0.06 | 0.08
Preexisting Conditions |  |  |  |  |  |  |  |  |
Stroke | 0.01 | 0.01 | 0.01 |  | 0.01 | 0.01 | 0.02 | 0.01 | -0.00
MI | 0.01 | 0.02 | 0.02 |  | 0.02 | 0.02 | 0.03 | -0.01 | 0.01
CHF | 0.00 | 0.01 | 0.01 |  | 0.01 | 0.00 | 0.01 | 0.02 | 0.01
Liver diseases | 0.02 | 0.02 | 0.02 |  | 0.02 | 0.02 | 0.02 | -0.01 | 0.03
Hypertension | 0.25 | 0.29 | 0.34 |  | 0.31 | 0.28 | 0.35 | 0.02 | 0.07
Fracture | 0.10 | 0.12 | 0.16 |  | 0.13 | 0.12 | 0.14 | 0.03 | 0.09
CABG/PTCA | 0.01 | 0.01 | 0.02 |  | 0.01 | 0.02 | 0.02 | -0.01 | 0.04
BRCA | 0.01 | 0.01 | 0.01 |  | 0.01 | 0.01 | 0.07 | 0.00 | 0.00

##  5 Summary

We proposed a statistical matching algorithm that constructs well-matched
pairs from large observational databases that resemble a target population. By
designing a representative matched sample, empirical researchers could
potentially (i) better reconcile the sometimes conflicting study findings, and
(ii) answer a clinical/epidemiological query for a scientifically meaningful
population. We applied the proposed method to investigate the discrepancy
between the Women’s Health Initiative (WHI) observational study and clinical
trial findings on HRT’s effect on coronary heart diseases. The method
facilitated some interesting findings. In particular, we found that a matched
cohort study constructed from the WHI observational data still supported a
cardioprotective effect of HRT, and this cardioprotective effect persisted
even after the designed matched samples were forced to resemble the WHI trial
population in their cardiovascular risk profile. However, in a matched design
that further resembles the WHI trial intervention group in the HRT initiation
age, the cardioprotctive effect seemed to disappear. Our findings provide some
evidence for the argument that HRT initiation age might have played an
important role in explaining the observational study and trial discrepancy
(Willett et al.,, 2006; Hernán et al.,, 2008).

## Data Availability Statement

The data that support the findings of this paper are available from the
National Heart, Lung, and Blood Institute. Restrictions apply to the
availability of these data, which were used under license for this paper. Data
are available from the author with the permission of the NHLBI.

## References

  * Ahuja et al., (1988) Ahuja, R. K., Magnanti, T. L., and Orlin, J. B. (1988).  Network Flows.  Cambridge, Mass.
  * Austin, (2011) Austin, P. C. (2011).  An introduction to propensity score methods for reducing the effects of confounding in observational studies.  Multivariate Behavioral Research, 46(3):399–424.
  * Bennett et al., (2020) Bennett, M., Vielma, J. P., and Zubizarreta, J. R. (2020).  Building representative matched samples with multi-valued treatments in large observational studies.  Journal of Computational and Graphical Statistics, 29(4):744–757.
  * Bertsekas, (1991) Bertsekas, D. P. (1991).  Linear Network Optimization: Algorithms and Codes.  MIT Press.
  * Bertsekas and Tseng, (1988) Bertsekas, D. P. and Tseng, P. (1988).  The relax codes for linear minimum cost network flow problems.  Annals of Operations Research, 13(1):125–190.
  * Cochran and Rubin, (1973) Cochran, W. G. and Rubin, D. B. (1973).  Controlling bias in observational studies: A review.  Sankhyā: The Indian Journal of Statistics, Series A, pages 417–446.
  * Dahabreh et al., (2019) Dahabreh, I. J., Robertson, S. E., Tchetgen, E. J., Stuart, E. A., and Hernán, M. A. (2019).  Generalizing causal inferences from individuals in randomized trials to all trial-eligible individuals.  Biometrics, 75(2):685–694.
  * Franklin et al., (2021) Franklin, J. M., Patorno, E., Desai, R. J., Glynn, R. J., Martin, D., Quinto, K., Pawar, A., Bessette, L. G., Lee, H., Garry, E. M., et al. (2021).  Emulating randomized clinical trials with nonrandomized real-world evidence studies: first results from the RCT DUPLICATE Initiative.  Circulation, 143(10):1002–1013.
  * Franklin et al., (2020) Franklin, J. M., Pawar, A., Martin, D., Glynn, R. J., Levenson, M., Temple, R., and Schneeweiss, S. (2020).  Nonrandomized real-world evidence to support regulatory decision making: process for a randomized trial replication project.  Clinical Pharmacology & Therapeutics, 107(4):817–826.
  * Gagnon-Bartsch and Shem-Tov, (2019) Gagnon-Bartsch, J. and Shem-Tov, Y. (2019).  The classification permutation test: A flexible approach to testing for covariate imbalance in observational studies.  Annals of Applied Statistics, 13(3):1464–1483.
  * Hansen, (2007) Hansen, B. B. (2007).  Optmatch: Flexible, optimal matching for observational studies.  R News, 7(2):18–24.
  * Hernán et al., (2008) Hernán, M. A., Alonso, A., Logan, R., Grodstein, F., Michels, K. B., Stampfer, M. J., Willett, W. C., Manson, J. E., and Robins, J. M. (2008).  Observational studies analyzed like randomized experiments: an application to postmenopausal hormone therapy and coronary heart disease.  Epidemiology (Cambridge, Mass.), 19(6):766.
  * Ho et al., (2007) Ho, D. E., Imai, K., King, G., and Stuart, E. A. (2007).  Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference.  Political Analysis, 15(3):199–236.
  * Humphrey et al., (2002) Humphrey, L. L., Chan, B. K., and Sox, H. C. (2002).  Postmenopausal hormone replacement therapy and the primary prevention of cardiovascular disease.  Annals of Internal Medicine, 137(4):273–284.
  * Kaplan and Meier, (1958) Kaplan, E. L. and Meier, P. (1958).  Nonparametric estimation from incomplete observations.  Journal of the American Statistical Association, 53(282):457–481.
  * Korte and Vygen, (2011) Korte, B. H. and Vygen, J. (2011).  Combinatorial Optimization, volume 1.  Springer.
  * Levina and Bickel, (2001) Levina, E. and Bickel, P. (2001).  The earth mover’s distance is the mallows distance: Some insights from statistics.  In Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 2, pages 251–256. IEEE.
  * MacKay et al., (2021) MacKay, E. J., Zhang, B., Heng, S., Ye, T., Neuman, M. D., Augoustides, J. G., Feinman, J. W., Desai, N. D., and Groeneveld, P. W. (2021).  Association between transesophageal echocardiography and clinical outcomes after coronary artery bypass graft surgery.  Journal of the American Society of Echocardiography.
  * Metkus et al., (2021) Metkus, T. S., Thibault, D., Grant, M. C., et al. (2021).  Transesophageal echocardiography in patients undergoing coronary artery bypass graft surgery.  Journal of the American College of Cardiology, 78(2):112–122.
  * Michels and Manson, (2003) Michels, K. B. and Manson, J. E. (2003).  Postmenopausal hormone therapy.  Circulation, 107(14):1830–1833.
  * O’Brien and Fleming, (1987) O’Brien, P. C. and Fleming, T. R. (1987).  A paired Prentice-Wilcoxon test for censored paired data.  Biometrics, pages 169–180.
  * Pimentel et al., (2015) Pimentel, S. D., Kelz, R. R., Silber, J. H., and Rosenbaum, P. R. (2015).  Large, sparse optimal matching with refined covariate balance in an observational study of the health outcomes produced by new surgeons.  Journal of the American Statistical Association, 110(510):515–527.
  * Prentice et al., (2005) Prentice, R. L., Langer, R., Stefanick, M. L., Howard, B. V., Pettinger, M., Anderson, G., Barad, D., Curb, J. D., Kotchen, J., Kuller, L., et al. (2005).  Combined postmenopausal hormone therapy and cardiovascular disease: toward resolving the discrepancy between observational studies and the women’s health initiative clinical trial.  American Journal of Epidemiology, 162(5):404–414.
  * Prentice et al., (2009) Prentice, R. L., Manson, J. E., Langer, R. D., Anderson, G. L., Pettinger, M., Jackson, R. D., Johnson, K. C., Kuller, L. H., Lane, D. S., Wactawski-Wende, J., et al. (2009).  Benefits and risks of postmenopausal hormone therapy when it is initiated soon after menopause.  American Journal of Epidemiology, 170(1):12–23.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rassen et al., (2012) Rassen, J. A., Shelat, A. A., Myers, J., Glynn, R. J., Rothman, K. J., and Schneeweiss, S. (2012).  One-to-many propensity score matching in cohort studies.  Pharmacoepidemiology and Drug Safety, 21:69–80.
  * Rosenbaum, (1989) Rosenbaum, P. R. (1989).  Optimal matching for observational studies.  Journal of the American Statistical Association, 84(408):1024–1032.
  * Rosenbaum, (2002) Rosenbaum, P. R. (2002).  Observational Studies.  Springer.
  * Rosenbaum, (2010) Rosenbaum, P. R. (2010).  Design of Observational Studies.  Springer.
  * Rosenbaum and Rubin, (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55.
  * Rosenbaum and Rubin, (1985) Rosenbaum, P. R. and Rubin, D. B. (1985).  Constructing a control group using multivariate matched sampling methods that incorporate the propensity score.  The American Statistician, 39(1):33–38.
  * Rubin, (1979) Rubin, D. B. (1979).  Using multivariate matched sampling and regression adjustment to control bias in observational studies.  Journal of the American Statistical Association, 74(366a):318–328.
  * Scirica et al., (2013) Scirica, B. M., Bhatt, D. L., Braunwald, E., Steg, P. G., Davidson, J., Hirshberg, B., Ohman, P., Frederich, R., Wiviott, S. D., Hoffman, E. B., et al. (2013).  Saxagliptin and cardiovascular outcomes in patients with type 2 diabetes mellitus.  New England Journal of Medicine, 369(14):1317–1326.
  * Silber et al., (2014) Silber, J. H., Rosenbaum, P. R., Ross, R. N., Ludwig, J. M., Wang, W., Niknam, B. A., Mukherjee, N., Saynisch, P. A., Even-Shoshan, O., Kelz, R. R., et al. (2014).  Template matching for auditing hospital cost and quality.  Health Services Research, 49(5):1446–1474.
  * Silber et al., (2001) Silber, J. H., Rosenbaum, P. R., Trudeau, M. E., Even-Shoshan, O., Chen, W., Zhang, X., and Mosher, R. E. (2001).  Multivariate matching and bias reduction in the surgical outcomes study.  Medical Care, pages 1048–1064.
  * Stuart, (2010) Stuart, E. A. (2010).  Matching methods for causal inference: A review and a look forward.  Statistical Science, 25(1):1–21.
  * Stuart et al., (2011) Stuart, E. A., Cole, S. R., Bradshaw, C. P., and Leaf, P. J. (2011).  The use of propensity scores to assess the generalizability of results from randomized trials.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 174(2):369–386.
  * The Women’s Health Initiative Study Group, (1998) The Women’s Health Initiative Study Group (1998).  Design of the women’s health initiative clinical trial and observational study.  Controlled Clinical Trials, 19(1):61–109.
  * Willett et al., (2006) Willett, W. C., Manson, J. E., Grodstein, F., Stampfer, M. J., and Colditz, G. A. (2006).  Re:“combined postmenopausal hormone therapy and cardiovascular disease: toward resolving the discrepancy between observational studies and the women’s health initiative clinical trial”.  American Journal of Epidemiology, 163(11):1067–1068.
  * Writing Group for the Women’s Health Initiative Investigators, (2002) Writing Group for the Women’s Health Initiative Investigators (2002).  Risks and benefits of estrogen plus progestin in healthy postmenopausal women: principal results from the women’s health initiative randomized controlled trial.  JAMA, 288(3):321–333.
  * Zhang et al., (2021) Zhang, B., Small, D., Lasater, K., McHugh, M., Silber, J., and Rosenbaum, P. (2021).  Matching one sample according to two criteria in observational studies.  Journal of the American Statistical Association, (just-accepted):1–34.
  * Zubizarreta, (2012) Zubizarreta, J. R. (2012).  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.