\section{Introduction}

\subsection{RCT DUPLICATE Initiative: Comparing Observational Studies to Randomized Controlled Trials (RCTs)}
\label{subsec: intro RCT duplicate}
In a recent high-profile study published in {\it Circulation}, the RCT DUPLICATE initiative (\citealp{franklin2020nonrandomized, franklin2021emulating}) 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, \citet{franklin2021emulating} 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. \citet{franklin2021emulating} concluded:
\begin{quote}\small
    [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.
\end{quote}
Discrepancies in age, race/ethnicity, and important preexisting comorbid conditions between RCT and observational study populations are common in \citet{franklin2021emulating}'s emulation studies. For instance, the Saxagliptin and Cardiovascular Outcomes in Patients with Type 2 Diabetes Mellitus (SAVOR-TIMI 53) study (\citealp{scirica2013saxagliptin}), one of the ten trials \citet{franklin2021emulating} 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 (\citealp[Table 1]{franklin2021emulating}). These differences persisted in the final matched samples constructed by \citet{franklin2021emulating}, 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 (e.g., the Women’s Health Initiative Study; see, e.g., \citet{prentice2005combined} and \citet{hernan2008observational}) motivate us to develop a transparent, easy-to-use algorithm that designs matched samples homogeneous in many covariates while mimicking a target population (e.g., the RCT population) in some key covariates (e.g., potential effect modifiers), so that the scientific conclusions drawn from the observational data could be better generalized to other context.


\subsection{A Naive Method, An Existing Method, and A New Approach}
One naive approach to creating homogeneous matched pairs resembling a target template would first select a subset of treated units from the observational data via matching on covariates collected by both the RCT and the observational study. In the second step, control units in the observational data are then matched to the treated units selected by the first step. This strategy suffers from a major drawback. Take the RCT DUPLICATE initiative as an example. Although the benchmark SAVOR-TIMI 53 trial reported only around $20$ baseline covariates, \citet{franklin2021emulating} matched on more than $120$ to guard against unmeasured confounding in their emulation with observational data. In the first step, there are many ways to design a smaller treated group similar to the RCT population in $20$ RCT covariates; however, it is difficult to determine which treated group designed from the first step should be utilized to form the eventual match in the second step. A selected treated group similar to the RCT population in $20$ RCT covariates could have poor overlap with the control units in the observational database in the other more than $100$ covariates; see Figure \ref{fig: illustrate naive method} for a toy example that illustrates this phenomenon. In many cases, we prefer well-matched treated and control groups to first ensure internal validity of the matched cohort study; similarity to some target population should be taken as a perk. The naive approach precludes a principled trade-off between matched cohort study's internal validity and its generalizability to a target population.

\begin{figure}[ht]
    \centering
   \begin{tikzpicture}

\node at (0.25, 1.25) {$X_1$};
\node at (1, 1.25) {$X_2$};
\draw[blue, very thick] (0,0) rectangle (0.5,0.5);
\draw[blue, very thick] (1,0.25) circle (0.25);
\node[text width=2cm, align = center] at (0.6, -4) {\small\textsf{Template}};

\node at (4.25, 1.25) {$X_1$};
\node at (5, 1.25) {$X_2$};
\node at (5.75, 1.25) {$X_3$};

\node at (3.25,0.25) {$\textsf{T1}$};
\node at (3.25,-0.5) {$\textsf{T2}$};
\draw[blue, very thick] (4,0) rectangle (4.5,0.5);
\draw[blue, very thick] (5,0.25) circle (0.25);
\draw[blue, very thick] (5.75,0.25) circle (0.25);

\draw[blue, very thick] (4,-0.75) rectangle (4.5,-0.25);
\draw[blue, very thick] (5,-0.5) circle (0.25);
\draw[blue, very thick] (5.5,-0.75) -- (5.75,-0.25) -- (6,-0.75) -- cycle;

\node[text width=2cm, align = center] at (5, -4) {\small\textsf{OBS \\ Treated}};



\node at (9, 1.25) {$X_1$};
\node at (9.75, 1.25) {$X_2$};
\node at (10.5, 1.25) {$X_3$};


\node at (8,0.25) {$\textsf{C1}$};
\node at (8,-0.5) {$\textsf{C2}$};
\node at (8,-1.25) {$\textsf{C3}$};
\node at (8,-2) {$\textsf{C4}$};

\draw[blue, very thick] (9,0.25) circle (0.25);
\draw[blue, very thick] (9.75,0.25) circle (0.25);
\draw[blue, very thick] (10.5,0.25) circle (0.25);

\draw[blue, very thick] (8.75,-0.75) rectangle (9.25,-0.25);
\draw[blue, very thick] (9.75,-0.5) circle (0.25);
\draw[blue, very thick] (10.5,-0.5) circle (0.25);

\draw[blue, very thick] (8.75,-1.5) rectangle (9.25,-1);
\draw[blue, very thick] (9.75,-1.25) circle (0.25);
\draw[blue, very thick] (10.25,-1.5) rectangle (10.75, -1);

\draw[blue, very thick] (8.75,-2.25) rectangle (9.25,-1.75);
\draw[blue, very thick] (9.5,-2.25) rectangle (10, -1.75);
\draw[blue, very thick] (10.25,-2.25) rectangle (10.75, -1.75);

\node[text width=2cm, align = center] at (9.75, -4) {\small\textsf{OBS \\ Control}};


\end{tikzpicture}
\caption{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 $\textsf{T1}$ and $\textsf{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, $\textsf{T1}$ is a better pick as no control unit has an $X_3$ value similar to that of $\textsf{T2}$.}
\label{fig: illustrate naive method}
\end{figure}

Recently, \citet{bennett2020building} developed a mixed-integer-programming-based (MIP-based) approach that builds representative matched samples. MIP-based matching algorithms (\citealp{li2001balanced,zubizarreta2012using,bennett2020building}) encode requirements for different versions of covariate balance as constraints in a 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 some researchers. This being said, if one finds no difficulty installing the solver and that the match can be readily delivered by MIP-based methods, then this is great and one should do that. Other approaches to adjusting for sample selection bias include weighting and doubly-robust methods (see, e.g., \citealp{cole2010generalizing,stuart2011use,Dahabreh:2019aa}); see \citet{colnet2020causal,degtiar2021review} for a more comprehensive review. 

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 (\citealp{rosenbaum2002observational,rosenbaum2010design, stuart2010matching, austin2011introduction,rassen2012one}). 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 \citet{hernan2008observational}.

\section{Methodology}
\label{sec: basic structure}
\subsection{Basic Network Structure: Vertices and Edges}
\label{subsec: basic structure v and e}
We describe and illustrate the proposed method using the RCT emulation study, although the idea can be immediately generalized to accommodate an arbitrary target population other than an RCT population. We refer to study units from the randomized controlled trial as RCT units and those from the observational database as OBS units for short. 

Figure \ref{fig: tripartite new scheme} 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:
\begin{equation}
    \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\}.
\end{equation}

An ordered pair of vertices is referred to as an edge in the network. The basic structure in Figure \ref{fig: tripartite new scheme} consists of the following edges:
\begin{equation}
    \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\}.
\end{equation}
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$. 

\begin{figure}[ht]
\centering
\begin{tikzpicture}[thick, color = black,
 
  fsnode/.style={circle, fill=black, inner sep = 0pt, minimum size = 5pt},
  ssnode/.style={circle, fill=black, inner sep = 0pt, minimum size = 5pt},
  shorten >= 3pt,shorten <= 3pt
]


\begin{scope}[start chain=going below,node distance=7mm]
\foreach \i in {1,2,3}
  \node[fsnode,on chain] (r\i) [label=above left: {\small$\kappa_\i$} ] {};
\end{scope}

\begin{scope}[xshift=3cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4}
  \node[ssnode,on chain] (t\i) [label=above right: {\small$\tau_\i$}] {};
\end{scope}

\begin{scope}[xshift=6cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4}
  \node[ssnode,on chain] (tt\i) [label=above left: {\small$\overline\tau_\i$}] {};
\end{scope}

\begin{scope}[xshift=9cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4,5,6}
  \node[ssnode,on chain] (c\i) [label=above right: {\small$\gamma_\i$}] {};
\end{scope}

\node [circle, fill = black, inner sep = 0pt, minimum size = 5pt, label=left: $\xi$] at (-2, -0.9) (source) {};

\node [circle, fill = black, inner sep = 0pt, minimum size = 5pt, label=right: $\overline\xi$ ] at (12, -2.2) (sink) {};



\foreach \i in {1,2,3} {
   \draw[color=gray] (source) -- (r\i);
   }

\foreach \i in {1,2,3} {
   \foreach \j in {1,2,3,4} {
   \draw[color=gray] (r\i) -- (t\j);
   }
 } 
 
\foreach \i in {1,2,3,4} {
   \draw[color=gray] (t\i) -- (tt\i);
}  
 
\foreach \i in {1,2,3,4} {
   \foreach \j in {1,2,3,4,5,6} {
   \draw[color=gray] (tt\i) -- (c\j);
   }
} 

\foreach \i in {1,2,3,4,5,6} {
   \draw[color=gray] (c\i) -- (sink);
}


\node[text width=2cm, align = center] at (0,1.5) {\small\textsf{RCT \\ Treated}};
\node[text width=2cm, align = center] at (3,1.5) {\small\textsf{OBS \\ Treated}};
\node[text width=2cm, align = center] at (6,1.5) {\small\textsf{OBS \\ Treated}};
\node[text width=2cm, align = center] at (9,1.5) {\small\textsf{OBS \\ Control}};

\end{tikzpicture}
\caption{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.}
\label{fig: tripartite new scheme}
\end{figure}


\subsection{Basic Network Structure: Capacity and Matched Samples}
\label{subsec: basic structure: c and c}
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 \citet{franklin2021emulating}'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 \ref{fig: tripartite new scheme}, $\text{cap}\{(\xi, \kappa_r)\} = k$ for $r = 1, \dots, R$, and all other edges have capacity $1$. In Figure \ref{fig: tripartite new scheme}, 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 (\citealp{ahuja1988network,bertsekas1991linear,rosenbaum1989optimal}). 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 \ref{fig: tripartite new scheme} 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 \ref{fig: tripartite new scheme with feasible flow} 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 
\begin{equation}
    \mathcal{M} = \left\{(\tau_t, \gamma_c) ~\text{such that}~f\{(\tau_t, \overline\tau_t)\} = f\{(\overline\tau_t, \gamma_c)\} = 1\right\}.
\end{equation}
For instance, the matched samples returned by the feasible flow in Figure \ref{fig: tripartite new scheme with feasible flow} consist of $\mathcal{M} = \{(\tau_1, \gamma_1), (\tau_2, \gamma_5), (\tau_3, \gamma_4)\}$.


\begin{figure}[ht]
\centering
\begin{tikzpicture}[thick, color = black,
 
  fsnode/.style={circle, fill=black, inner sep = 0pt, minimum size = 5pt},
  ssnode/.style={circle, fill=black, inner sep = 0pt, minimum size = 5pt},
  shorten >= 3pt,shorten <= 3pt
]


\begin{scope}[start chain=going below,node distance=7mm]
\foreach \i in {1,2,3}
  \node[fsnode,on chain] (r\i) [label=above left: {\small$\kappa_\i$} ] {};
\end{scope}

\begin{scope}[xshift=3cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4}
  \node[ssnode,on chain] (t\i) [label=above right: {\small$\tau_\i$}] {};
\end{scope}

\begin{scope}[xshift=6cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4}
  \node[ssnode,on chain] (tt\i) [label=above left: {\small$\overline\tau_\i$}] {};
\end{scope}

\begin{scope}[xshift=9cm,yshift=0cm,start chain=going below,node distance=7mm]
\foreach \i in {1,2,3,4,5,6}
  \node[ssnode,on chain] (c\i) [label=above right: {\small$\gamma_\i$}] {};
\end{scope}

\node [circle, fill = black, inner sep = 0pt, minimum size = 5pt, label=left: $\xi$] at (-2, -0.9) (source) {};

\node [circle, fill = black, inner sep = 0pt, minimum size = 5pt, label=right: $\overline\xi$ ] at (12, -2.2) (sink) {};



\foreach \i in {1,2,3} {
   \draw[color = black, line width = 0.5mm] (source) -- (r\i);
   }

\foreach \i in {1,2,3} {
   \foreach \j in {1,2,3,4} {
   \draw[color=gray!50] (r\i) -- (t\j);
   }
 } 
 
\draw[color = black, line width = 0.5mm] (r1) -- (t2);
\draw[color = black, line width = 0.5mm] (r2) -- (t3);
\draw[color = black, line width = 0.5mm] (r3) -- (t1);

 
\foreach \i in {1,2,3,4} {
   \draw[color=gray!50] (t\i) -- (tt\i);
}  

\draw[color = black, line width = 0.5mm] (tt2) -- (t2);
\draw[color = black, line width = 0.5mm] (tt3) -- (t3);
\draw[color = black, line width = 0.5mm] (tt1) -- (t1);

 
\foreach \i in {1,2,3,4} {
   \foreach \j in {1,2,3,4,5,6} {
   \draw[color=gray!50] (tt\i) -- (c\j);
   }
} 

\foreach \i in {1,2,3,4,5,6} {
   \draw[color=gray!50] (c\i) -- (sink);
}

\draw[color = black, line width = 0.5mm] (tt2) -- (c5);
\draw[color = black, line width = 0.5mm] (tt3) -- (c4);
\draw[color = black, line width = 0.5mm] (tt1) -- (c1);

\draw[color = black, line width = 0.5mm] (c1) -- (sink);
\draw[color = black, line width = 0.5mm] (c4) -- (sink);
\draw[color = black, line width = 0.5mm] (c5) -- (sink);




\node[text width=2cm, align = center] at (0,1.5) {\small\textsf{RCT \\ Treated}};
\node[text width=2cm, align = center] at (3,1.5) {\small\textsf{OBS \\ Treated}};
\node[text width=2cm, align = center] at (6,1.5) {\small\textsf{OBS \\ Treated}};
\node[text width=2cm, align = center] at (9,1.5) {\small\textsf{OBS \\ Control}};

\node[text width=2cm, align = center] at (7.5,-4.5) {\Large\textsf{Right}};
\node[text width=2cm, align = center] at (1.5,-4.5) {\Large\textsf{Left}};

\end{tikzpicture}
\caption{A feasible flow $f(\cdot): \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)\}$.}
\label{fig: tripartite new scheme with feasible flow}
\end{figure}

\subsection{Basic Network Structure: Cost, Probability of Participation, Propensity Score}
\label{subsec: basic structure: cost}
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 \ref{fig: tripartite new scheme}, 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.,
\begin{equation}
    \text{cost}\{(\kappa_r, \tau_t)\} = \delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}}).
\end{equation}
 One intriguing strategy equals $\delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}})$ to a scalar ``balancing score" of $\widetilde{\boldsymbol{x}}$ (\citealp{rosenbaum1983central}), 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, \citet{stuart2011use} 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 (\citealp{stuart2011use}).

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.,
\begin{equation}
    \text{cost}\{(\overline\tau_t, \gamma_c)\} = \Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}.
\end{equation}
Different specifications of $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}$ have been extensively studied in the literature (\citealp{rosenbaum2002observational,rosenbaum2010design, stuart2010matching}). Perhaps the most widely used strategy is to equal $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}$ to the absolute difference in \citet{rosenbaum1983central}'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})$ (\citealp{cochran1973controlling,rubin1979using}), or the Mahalanobis distance within a propensity score caliper (\citealp{rosenbaum1985constructing}). 

\subsection{Minimum Cost Flow, Complexity, Trade-Off Between Internal and External Validity}
\label{subsec: tradeoff}
The cost of a feasible flow $f(\cdot)$ in the proposed network is equal to
\begin{equation}
\label{eq: cost of a feasible flow}
    \text{cost}(f) = \sum_{e \in \mathcal{E}} f(e)\cdot \text{cost}(e) = \underbrace{\sum_{\substack{(\kappa_r, \tau_t) \in \mathcal{E}: \\f\{(\kappa_r, \tau_t)\} = 1}} \delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}})}_{S1} ~+ \underbrace{\sum_{\substack{(\overline\tau_t, \gamma_c) \in \mathcal{E}: \\f\{(\overline\tau_t, \gamma_c)\} = 1}} \Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}}_{S2},
\end{equation}
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 (\citealp{korte2011combinatorial}) 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 $\textsf{R}$, the minimum cost flow can be found via \citet{bertsekas1991linear}'s auction algorithm implemented by \citet{bertsekas1988relax} and made available by \citet{hansen2007optmatch} and \citet{pimentel2015large}.

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, \citet{franklin2021emulating} 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_{\substack{(\overline\tau_t, \gamma_c) \in \mathcal{E}: \\f\{(\overline\tau_t, \gamma_c)\} = 1}} \Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}$ in expression \eqref{eq: cost of a feasible flow} 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_{\substack{(\kappa_r, \tau_t) \in \mathcal{E}: \\f\{(\kappa_r, \tau_t)\} = 1}} \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 enable researchers explore the trade-off between matched samples' internal and external validity, 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 a weighted average of $S1$ and $S2$:
\begin{equation}
\label{eq: cost of a feasible flow weighted}
    \text{cost}(f) = \sum_{e \in \mathcal{E}} f(e)\cdot \text{cost}(e) = \lambda\cdot\bigg\{\sum_{\substack{(\kappa_r, \tau_t) \in \mathcal{E}: \\f\{(\kappa_r, \tau_t)\} = 1}} \delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}})\bigg\} ~+ \lambda\cdot\bigg\{\sum_{\substack{(\overline\tau_t, \gamma_c) \in \mathcal{E}: \\f\{(\overline\tau_t, \gamma_c)\} = 1}} \Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}\bigg\}.
\end{equation}
A small $\lambda$ value gives priority to matched samples' internal validity, while a large $\lambda$ value generalizability to the target population. A similar weighting scheme is also used in \citet{zhang2021matching} but for a different purpose.



\section{Simulation Study and Software Availability}
\label{sec: simulation and software}

We illustrated the proposed method with an example in Supplemental Material C, and described a comprehensive simulation study in Supplemental Material D. In the simulation study, we examined how the matched samples delivered by the proposed method improves the performance of the downstream outcome analysis. In particular, we considered settings where the estimand of interest is the sample average treatment effect of a target template, the treatment effect is heterogeneous, and the effect modifiers' distributions are different in the target and observational data. Our proposed algorithm is able to achieve significant bias reduction compared to statistical matching algorithms that pays no attention to the target template. The proposed matching algorithm is integrated in the package \textsf{match2C} available via the statistical computing software \textsf{R} (\citealp{R_software}) with a detailed tutorial.


\section{Revisiting the Women’s Health Initiative (WHI) Study}
\subsection{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 \citet{study1998design}. 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 (\citealp{study1998design}). The WHI estrogen-plus-progestin (E + P for short) trial found a rather surprising elevation in coronary heart disease risk (\citealp{writing2002risks}) and sparked a lot of discussion concerning 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 (\citealp{humphrey2002postmenopausal, rutter2007identifying,Yu2021_WHI}); (ii) biological differences between trial participants and those in the observational study (\citealp{Michels2003editorial}); (iii) differences in time since menopause at hormone therapy initiation (\citealp{prentice2005combined,willett2006re, rossouw2007postmenopausal,hernan2008observational, prentice2009benefits}). Table \ref{tb: WHI baseline characteristics} 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, have slightly smaller BMI and 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.

\begin{table}[ht]
\caption{Important baseline characteristics of the WHI observational study and WHI trial subjects. Mean (SE) are reported for continuous variables and count (\%) for categorical variables.}
\label{tb: WHI baseline characteristics}
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{lcccccc}
  \hline \\ [-0.8em]
  &\multicolumn{2}{c}{WHI Observational Data} &&\multicolumn{2}{c}{WHI Trial}\\ \\ [-0.8em] \cline{2-3}\cline{5-6}
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{Never/Past} \\ \textbf{Users} \\($\mathbf{n = 75303}$)\end{tabular}}
  & \multirow{4}{*}{\begin{tabular}{c}\textbf{Current} \\ \textbf{Users} \\($\mathbf{n = 18340}$)\end{tabular}}&
  & \multirow{4}{*}{\begin{tabular}{c} \textbf{Control} \\($\mathbf{n = 8102}$)\end{tabular}}
  & \multirow{4}{*}{\begin{tabular}{c}\textbf{E+P} \\ \textbf{Intervention} \\($\mathbf{n = 8506}$)\end{tabular}}\\ \\ \\ \\
  
  \textbf{Age at screening} &  64.30 (7.37) &  60.84 (6.69) &&  63.33 (7.11) &  63.23 (7.13) \\ 
  \textbf{Race/Ethnicity}  \\ 
  \hspace{0.5cm}White &  61704 (81.9)  &  16285 (88.8)  &&   6805 (84.0)  &   7141 (84.0) \\ 
  \hspace{0.5cm}Black/Hispanic &  10166 (13.5)  &   1074 ( 5.9)  &&    989 (12.2)  &   1019 (12.0) \\ 
  \hspace{0.5cm}Other &   3433 ( 4.6)  &    981 ( 5.3)  &&    308 ( 3.8)  &    346 ( 4.1)   \\ 
  \textbf{Education} \\ 
  \hspace{0.5cm}College or above &  30418 (40.4)  &  10044 (54.8)  &&   3011 (37.2)  &   3111 (36.6)  \\ 
  \hspace{0.5cm}Some college &  28074 (37.3)  &   5853 (31.9)  &&   3060 (37.8)  &   3357 (39.5) \\ 
  \hspace{0.5cm}High school diploma/GED &  13068 (17.4)  &   2050 (11.2)  &&   1609 (19.9)  &   1615 (19.0) \\ 
  \hspace{0.5cm}Other &    3743 ( 5.0)  &    393 (2.1)  &&     422 (5.2)  &  423 (5.0) \\ 

  \textbf{Blood Pressure} \\
  \hspace{0.5cm}Systolic & 128 (18.12) & 123 (16.92) && 128 (17.53) & 128 (17.63)  \\ 
  \hspace{0.5cm}Diastolic &  75 (9.42) &  74 (8.99) &&  76 (9.09) &  76 (9.12)   \\ 
  \textbf{BMI} &  27.62 (5.98) &  25.84 (5.15) &&  28.50 (5.91) &  28.46 (5.82)   \\ 
  \textbf{Smoking} \\ 
  \hspace{0.5cm}NA &   1004 ( 1.3)  &    220 ( 1.2)  &&     98 ( 1.2)  &     83 ( 1.0) \\ 
  \hspace{0.5cm}Current Smoker &   4850 ( 6.4)  &    940 ( 5.1)  &&    838 (10.4)  &    880 (10.3)   \\ 
  \hspace{0.5cm}Never Smoked &  38296 (50.9)  &   8710 (47.5)  &&   3999 (49.4)  &   4178 (49.1)   \\ 
  \hspace{0.5cm}Past Smoker &  31048 (41.3)  &   8448 (46.1)  &&   3157 (39.0)  &   3362 (39.5)   \\ 
  \textbf{No. of PA episodes} \\
  \hspace{0.5cm}Total &   5.28 (4.12) &   5.80 (4.12) &&   4.77 (4.06) &   4.74 (4.10)  \\ 
  \hspace{0.5cm}Medium to strenuous &   2.94 (3.37) &   3.59 (3.57) &&   2.58 (3.24) &   2.50 (3.21)   \\ 
   \textbf{HRT Initiation}  \\
  \hspace{0.5cm}NA &  67259 (89.3)  & 0 &&   6706 (82.8)  &      0   \\ 
  \hspace{0.5cm}Age at initiation   &  52.04 (7.35) &  53.87 (6.84) &&  53.57 (6.53) &  61.93 (8.15) \\ 
  \textbf{Previous E+P use in years} &   0.47 (2.12) &   7.02 (5.64) &&   0.64 (2.23) &   0.70 (2.34)   \\ 
   \hline
\end{tabular}}
\end{table}

Our proposed matching algorithm is suited for addressing concerns (i), (ii), and (iii). The first concern of the internal validity of the observational study can be greatly alleviated if matched OBS samples are balanced for a large number of baseline covariates including detailed demographics, preexisting comorbid conditions, personal habits, etc. Rich covariate information collected by the WHI study offered us an opportunity to do this. The second and third concerns about the comparability between OBS and RCT samples in their cardiovascular risk profile and HRT initiation time can be at least partially addressed by treating the RCT samples as a template and forcing matched OBS samples to resemble the RCT cohort in these aspects. 

Reconciling RCT and OBS results goes beyond addressing these three concerns. One 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 (\citealp{willett2006re}) while most trial participants were never-users of the therapy by the time of randomization. The lack of overlap between RCT and OBS samples in previous estrogen-plus-progestin usage (due to the study design) makes it virtually impossible to balance this aspect between RCT and OBS samples using any matching tool; nevertheless, in our opinion, addressing or just alleviating some major concerns (e.g., concerns (i), (ii), and (iii)) helps researchers critically examine the remaining explanations of the discrepancy and is a meaningful step towards ultimately reconciling the RCT and OBS findings.

\subsection{Template and Three Matched Samples}
We aim to construct well-matched samples from the WHI observational study that resemble the WHI trial units in the following sense: (i) matched treated and control units resemble the WHI trial units in cardiovascular risk profile, and (ii) matched treated units resemble the WHI trial intervention group in the HRT initiation time.

To this end, our template consists of a random sample of $1000$ WHI trial units in the intervention group with the following covariates: risk factors listed in Table \ref{tb: WHI baseline characteristics} plus the HRT initiation time. In many applications, researchers may not have individual-level data but only the joint distribution of important baseline covariates; see, e.g., \citet{cole2010generalizing} when they try to generalize the AIDS Clinical Trial Group 320 Study (ACTG 320 Trial) to the HIV-infected individuals in the United States. If this is the case, researchers could readily simulate individual-level data subject to the marginal joint distribution constraints, as suggested by \citet{bennett2020building}.

We applied the proposed matching algorithm to constructing matched samples of different sizes corresponding to choosing different $k$'s in Section \ref{subsec: basic structure: c and c}. Our desired matched OBS samples would have similar cardiovascular risk profile as the template, and 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 particular, in addition to risk factors listed in Table \ref{tb: WHI baseline characteristics}, we further matched on subjects' region, partner's education level, income level, marital status, reproductive history, and eight important preexisting conditions.

Table \ref{tbl: matched samples I} summarizes results from two matched samples constructed using the proposed template matching algorithm with different parameters. Match $\textsf{M1}$ constructed $10,000$ matched pairs of two units: one OBS treated unit and the other OBS control unit. Match $\textsf{M1}$ used a Mahalanobis distance with an estimated generalizability score caliper (caliper size = $0.05$) on the left and a Mahalanobis distance with an estimated propensity score caliper (caliper size = $0.15$) on the right; see Figure \ref{fig: tripartite new scheme with feasible flow}. Because the large sample size ($\approx 100,000$ in the observational study), we applied a ``hard" caliper in the sense that edges connecting OBS treated and OBS control units are removed whenever the two units differ in their estimated propensity score by more than the caliper size; in this way, the network structure is largely sparsified and computation is boosted (\citealp{yu2020matching}). We set $\lambda = 10$ to give some priority to mimicking the target template. Match $\textsf{M2}$ is similar to $\textsf{M1}$ except that we only formed $3,000$ matched pairs. Lastly, match $\textsf{M0}$ formed $18,340$ matched pairs exhausting each and every treated OBS unit. The balance table and distributions of the estimated propensity score in the treated and matched control groups of $\textsf{M0}$ can be found in the Supplemental Material E.

Judging from internal validity, all three matched samples $\textsf{M0}$, $\textsf{M1}$, and $\textsf{M2}$ are acceptable; in fact, the absolute standardized mean differences (SMDs) of all cardiovascular risk factors and additional OBS covariates are less than $0.1$, or one-tenth of one pooled standard deviation (\citealp{silber2001multivariate,austin2015moving}); in fact, most absolute SMDs are much less than $0.1$. However, the three matches differ, sometimes significantly, in their resemblance to the target template. In particular, matched samples $\textsf{M0}$ differ from the target template in its cardiovascular risk profile (e.g., the percentage of black/Hispanic is $6\%$ in $\textsf{M0}$ compared to $12\%$ in the template; the percentage with a college degree or above is $55\%$ in $\textsf{M0}$ compared to merely $38\%$ in the template), and HRT initiation age ($54$ in $\textsf{M0}$ compared to $62$ in the template, which translates to a $6.30$-year difference in the previous HRT usage). 

On the other hand, both $\textsf{M1}$ and $\textsf{M2}$ are similar to the trial population in the cardiovascular risk profile: the percentage of white (black/Hispanic) women was $85\%$ ($11\%$) in $\textsf{M1}$, $84\%$ ($12\%$) in $\textsf{M2}$, and $84\%$ ($12\%$) in the template, compared to $89\%$ ($6\%$) in $\textsf{M0}$; the percentage with a college degree or above is $40\%$ in $\textsf{M1}$, $36\%$ in $\textsf{M2}$, and $38\%$ in the template, compared to $55\%$ in $\textsf{M0}$. Moreover, the treated group of $\textsf{M1}$ now initiated their HRT at an average age of $57$ with an average $5.61$-year of HRT usage prior to enrolling in the WHI study, and the treated group of $\textsf{M2}$ initiated their HRT at an average age of $61$ with an average of $3.28$-year previous HRT usage. The HRT initiation age of treated units in \textsf{M2} is now very close to that of the trial population, with a difference less than one year, though there is still a potentially meaningful discrepancy in the previous HRT usage ($3.28$ years in \textsf{M2} vs $0.70$ years in the template).

\begin{table}[ht]
    \centering
    \begin{tabular}{ccccccc}\hline \\ [-0.6em]
         & & &\multicolumn{3}{c}{\textsf{Resemblance to the Template}}\\ \\ [-0.6em] \cline{4-6}
         &\multirow{3}{*}{\begin{tabular}{c}{\small \textsf{Study}} \\ {\small \textsf{Samples}}\end{tabular}} & \multirow{3}{*}{\begin{tabular}{c}{\small\textsf{Internal}} \\ {\small\textsf{Validity}} \end{tabular}} &\multirow{3}{*}{\begin{tabular}{c}{\small\textsf{Cardiovascular}} \\ {\small\textsf{Risk Profile}} \end{tabular}}
         &\multirow{3}{*}{\begin{tabular}{c}{\small\textsf{HRT}} \\{\small\textsf{Initiation}} \end{tabular}}
         &\multirow{3}{*}{\begin{tabular}{c}{\small\textsf{Previous}} \\ {\small\textsf{E+P Use}} \end{tabular}} 
         &\multirow{3}{*}{\begin{tabular}{c}{\small\textsf{Cardioprotective}} \\ {\small\textsf{Effect of HRT}} \end{tabular}}\\ \\ \\
         & \textsf{Before} & Very Poor & Poor & Very Poor & Very Poor & Very Large\\
         & \textsf{M0} & Good & Poor & Very Poor & Very Poor & Persisted\\
         & \textsf{M1} &Good & Good & Poor & Poor & Persisted\\
         & \textsf{M2} &Good & Good & Good & Fair & Disappeared\\ \hline
    \end{tabular}
    \caption{A qualitative summary of the internal and external validity of observational study data before matching and three matched samples.}
    \label{tbl: HRT summary three matches}
\end{table}


Table \ref{tbl: HRT summary three matches} gives a qualitative summary of three different study designs (\textsf{M0}, \textsf{M1}, and \textsf{M2}): in short, while all three designs have stellar internal validity, \textsf{M2} most resembles the target template, with \textsf{M1} coming in the second, and \textsf{M0} the last. 


\subsection{Comparing Survival Outcomes}
Figure \ref{fig: real data K-M curve} plots the Kaplan-Meier curves (\citealp{kaplan1958nonparametric}) in the treatment and control groups in each study design. Compared to the entire reservoir of control OBS units in the unadjusted analysis, matched control units in \textsf{M0} have much higher survival probabilities, reflecting the fact that these matched controls are more identical to the treated units and a lot healthier; however, the cardioprotective effect of HRT still exists ($P$-value = $0.02$ according to \citet{obrien1987paired}'s matched-pair Prentice-Wilcoxon test) in \textsf{M0} and is consistent with the previous analysis of the WHI observational data (\citealp{prentice2005combined,prentice2006combined}). Compared to their counterparts in \textsf{M0}, both treated and matched comparison groups in \textsf{M1} have lower survival probabilities; however, a qualitatively similar cardioprotective effect persists in the design \textsf{M1} ($P$-value = $0.09$). Recall that \textsf{M1} has largely mimicked the cardiovascular risk profile of the target template; therefore, it seems that the difference in cardiovascular risk profile between observational study and trial populations is \emph{not} sufficient in explaining the inconsistency in the trial and observational study results. Compared to \textsf{M0} and \textsf{M1}, \textsf{M2} best resembles the trial population in their HRT initiation age and previous HRT use; indeed, when we examine the survival outcomes of matched OBS units in \textsf{M2}, the cardioproctective effect seemed to disappear ($P$-value = $0.65$). Three comparisons facilitated by three different study designs \textsf{M0}, \textsf{M1}, and \textsf{M2} seemed to support \citet{willett2006re}'s assessment that HRT initiation age and previous HRT usage played a key role in the discrepancy between WHI observational study and trial findings, and resonates with similar findings in \citet{hernan2008observational}.


\begin{figure}[ht]
   \centering
     \subfloat[Unadjusted]{\includegraphics[width = 0.49\columnwidth]{unadj.pdf}}
     \subfloat[Design $\textsf{M0}$]{\includegraphics[width = 0.49\columnwidth]{all.pdf}}\\
     \subfloat[Design $\textsf{M1}$]{\includegraphics[width = 0.49\columnwidth]{k10.pdf}}
      \subfloat[Design $\textsf{M2}$]{\includegraphics[width = 0.49\columnwidth]{k3.pdf}}
     \caption{Survival outcomes for coronary heart disease (CHD) in the unadjusted observational study samples (top left), design \textsf{M0} (top right), design \textsf{M1} (bottom left), and design \textsf{M2} (bottom right).}
\label{fig: real data K-M curve}
\end{figure}


\begin{table}[ht]
\centering
\caption{\small Balance table of baseline covariates in WHI observational study before matching, two matched samples constructed using the proposed algorithm, and the matched template. Match $\textsf{M1}$ constructed $10000$ matched pairs and Match $\textsf{M2}$ $3000$ matched pairs.}
\label{tbl: matched samples I}
\resizebox{\textwidth}{!}{
\begin{tabular}{lccccccccc}
  \hline
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{OBS} \\ \textbf{Treated}\end{tabular}} 
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{Matched} \\ \textbf{Treated} \\ \textbf{M1}\end{tabular}}
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{Matched} \\ \textbf{Treated} \\\textbf{M2}\end{tabular}}
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{Template}\end{tabular}}
   &\multirow{4}{*}{\begin{tabular}{c}\textbf{Matched} \\ \textbf{Control} \\\textbf{M2}\end{tabular}}
   &\multirow{4}{*}{\begin{tabular}{c}\textbf{Matched} \\ \textbf{Control} \\\textbf{M1}\end{tabular}}
    &\multirow{4}{*}{\begin{tabular}{c}\textbf{OBS} \\ \textbf{Control}\end{tabular}}
   &\multirow{4}{*}{\begin{tabular}{c}\textbf{SMD} \\ \textbf{M1}\end{tabular}}
    &\multirow{4}{*}{\begin{tabular}{c}\textbf{SMD} \\ \textbf{M2}\end{tabular}}
   \\ \\ \\ \\ \\ [-0.6em]
   
      \multicolumn{10}{c}{\large\textbf{Sample Size}}\\  \\ [-0.6em]
   \hspace{0.2cm}\textbf{N} &18340 & 10000 &3000 &1000 &3000 &10000 &75303 & \\ \\ [-0.6em]
   
   \multicolumn{10}{c}{\large\textbf{HRT Prior Usage}}\\  \\ [-0.6em]
   \textbf{Age at initiation} \\
   \hspace{0.2cm}NA & 0 &0 &0 &0 &0.88 &0.87 &0.89\\
   \hspace{0.2cm}Age & 53.87 &57.02 &60.96 &61.93 &52.95 &52.21 &52.04\\
   \textbf{Previous E+P use, yrs} &7.02 &5.61 &3.28 &0.70 &0.50 &0.51  &0.47\\ \\ [-0.6em]
   
   \multicolumn{10}{c}{\large\textbf{Covariates Collected in RCT and OBS}}\\  \\ [-0.6em]
\textbf{Age at screening} & 60.84 &62.49 & 63.93 & 63.10 & 63.71 &62.40 & 64.30 & 0.01 & 0.03 \\ 
\textbf{Race/Ethnicity}\\
 \hspace{0.2cm}White & 0.89 &0.85 & 0.84 & 0.84 & 0.84 &0.85 &0.82 &-0.00 & -0.00 \\ 
   \hspace{0.2cm}Black/Hispanic & 0.06 &0.11 & 0.12 & 0.12 & 0.12 &0.11 & 0.14 &0.00 & 0.00 \\ 
   \textbf{Education}\\
 \hspace{0.2cm}College or above & 0.55 &0.40 & 0.36 & 0.38 & 0.36 &0.40 & 0.40 &0.00 & 0.00 \\ 
  \hspace{0.2cm} Some college & 0.32 &0.41 & 0.43 & 0.40 & 0.41 &0.41 & 0.37 &0.01 & 0.04 \\ 
  \hspace{0.2cm} High school diploma/GED & 0.11 &0.16 & 0.17 & 0.17 & 0.18 &0.16 & 0.17 &-0.00 & -0.04 \\ 
  \textbf{Blood pressure}\\
  \hspace{0.2cm}Systolic & 123.27 & 125.65 & 127.23 & 126.99 & 125.43 &127.71 & 127.86 &0.01 & -0.03 \\
  \hspace{0.2cm}Diastolic & 74.03 & 74.47 & 75.00 & 75.66 & 74.89 & 74.57 & 74.91 &-0.01  &0.01 \\ 
  \textbf{BMI} & 25.85 & 27.10 & 28.19 & 28.46 & 27.70 & 26.78 & 27.61 &0.06 &0.09 \\
 \textbf{Smoking}\\
 \hspace{0.2cm}Current smoker & 0.05 &0.09 & 0.11 & 0.10 & 0.11 &0.09 &0.06 &0.01 &0.00 \\ 
  \hspace{0.2cm}Never smoked & 0.47 &0.47 & 0.48 & 0.51 & 0.47 &0.47 &0.51 &0.00 & 0.03 \\ 
  \hspace{0.2cm}Past smoker & 0.46 &0.42 & 0.40 & 0.37 & 0.41 &0.43 & 0.41 &-0.01 &-0.03 \\ 
  \textbf{No. of PA episodes}\\
  \hspace{0.2cm}Total & 5.80 &5.25 & 4.77 & 4.71 & 4.79 &5.33 &5.29 &-0.02 &-0.00 \\ 
  \hspace{0.2cm}Medium to strenuous & 3.58 &3.05 & 2.61 & 2.56 & 2.57 &3.07 & 2.94 &-0.00 &0.01 \\ \\
   \multicolumn{10}{c}{\large\textbf{Additional OBS Covariates}}\\  \\ [-0.6em]
  \textbf{Region}\\
  \hspace{0.2cm}Midwest & 0.22 &0.22 & 0.23 &  & 0.23 &0.23 &0.22 &-0.01 &0.00 \\ 
  \hspace{0.2cm}Northeast & 0.18 &0.19 &0.20 &  & 0.23 &0.20 &0.24 &-0.01 &-0.06 \\ 
  \hspace{0.2cm}South & 0.25 &0.25 & 0.25 &  & 0.25 &0.26 & 0.26 &-0.02 &-0.00 \\ 
     \textbf{Partner's education}\\
  \hspace{0.2cm}College or above & 0.43 &0.36 & 0.31 &  & 0.30 &0.34 &0.30 &0.03 & 0.03 \\ 
  \hspace{0.2cm}Some college & 0.16 &0.19 & 0.18 &  & 0.18 &0.18 & 0.17 &0.00 &0.00 \\ 
   \hspace{0.2cm}High school diploma/GED & 0.07 &0.08 & 0.09 &  & 0.09 &0.08 & 0.09 &-0.00 &0.01 \\ 
   \textbf{Income}\\
   \hspace{0.2cm}Below 35K & 0.23 &0.29 & 0.35 &  & 0.37 &0.31 &0.40 &-0.03 & -0.04 \\ 
   \hspace{0.2cm} 35K - 75K & 0.42 &0.43 & 0.42 &  & 0.39 &0.42 & 0.36 &0.01 &0.06 \\ 
   \hspace{0.2cm}Above 75K & 0.29 &0.21 &0.17 &  & 0.17 &0.20 & 0.16 &0.03 &-0.00 \\ 
   \textbf{Marital status}\\
   \hspace{0.2cm}Married & 0.69 &0.66 &0.62 &  & 0.61 &0.65 &0.60 &0.03 & 0.02 \\ 
   \hspace{0.2cm}Divorced/Widowed & 0.26 &0.29 & 0.33 &  & 0.35 &0.31 & 0.35 &-0.03 &-0.04 \\ 
 \textbf{Employment status}\\
  \hspace{0.2cm}Yes & 0.45 &0.39 & 0.34 &  & 0.34 &0.39 & 0.32 &0.00 & 0.00 \\ 
  \hspace{0.2cm}No & 0.53 &0.59 & 0.64 &  & 0.64 &0.59 & 0.65 &0.00 & 0.00 \\ 
\textbf{Reproductive history}\\  
  \hspace{0.2cm}No ovary removed & 0.92 &0.91 &0.92 &  & 0.92 &0.91 & 0.64 &0.00 & 0.01 \\ 
\hspace{0.2cm}Oral contraceptive use ever & 0.53 &0.46 & 0.41 &  & 0.41 &0.46 & 0.37 &0.00 & 0.01 \\ 
  \hspace{0.2cm}OC duration in years & 5.59 &5.52 & 5.43 &  & 5.15 &5.27 &5.16 &0.07 & 0.08 \\ 
  \textbf{Preexisting Conditions}\\
  \hspace{0.2cm}Stroke & 0.01 &0.01 & 0.01 &  & 0.01 &0.01 & 0.02 &0.01 & 0.00 \\ 
  \hspace{0.2cm}MI & 0.01 &0.02 &0.02 &  & 0.02 &0.02 &0.03 &0.02 & 0.05 \\ 
  \hspace{0.2cm}CHF & 0.00 &0.01 & 0.01 &  & 0.01 &0.01 &0.01 &0.01 & 0.01 \\ 
  \hspace{0.2cm}Liver diseases & 0.02 &0.02 & 0.02 &  & 0.03 &0.02 & 0.02 &-0.03 &-0.02 \\ 
  \hspace{0.2cm}Hypertension & 0.25 &0.29 &0.34 &  & 0.33 &0.28 & 0.35 &0.03 & 0.02 \\ 
  \hspace{0.2cm}Fracture & 0.10 &0.13 & 0.15 &  & 0.14 &0.13 &0.14 &-0.00 & 0.04 \\ 
  \hspace{0.2cm}CABG/PTCA & 0.01 &0.02 & 0.02 &  & 0.02 &0.02 & 0.02 &0.01 & 0.02 \\ 
  \hspace{0.2cm}BRCA & 0.01 &0.01 & 0.01 &  & 0.01 &0.01 & 0.07 & 0.00 &0.00 \\ 
   \hline
\end{tabular}}
\end{table}



\section{Summary}
Motivated by the ongoing work undertaken by the RCT DUPLICATE Initiative, we proposed a statistical matching algorithm that constructs well-matched samples from large observational databases while mimicking a target population like the RCT population. By designing matched samples that resemble a target population, empirical researchers could potentially (i) better reconcile the sometimes conflicting RCT and observational 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 mimic the WHI trial population in the cardiovascular risk profile. However, after we further designed matched samples resembling the WHI trial intervention group in the HRT initiation age, the cardioprotctive effect seemed to disappear. Our findings seemed to provide some evidence for the argument that HRT initiation age might have played a major role in explaining the observational study and trial discrepancy (\citealp{willett2006re, hernan2008observational}).


\subsection{Statistical Matching, Network Flow, Recent Advancement}
\textcolor{red}{Move this section to the Appendix.}

Matching and subclassification are widely used in empirical research to embed noisy, non-randomized observational data into an approximate randomized controlled experiment, and facilitate analyzing data as such (\citealp{ho2007matching, rubin2007design,rubin2008objective,rosenbaum2002observational,rosenbaum2010design,rosenbaum2020modern,stuart2010matching,rassen2012one,imbens2015matching}). In a seminal paper, \citet{rosenbaum1989optimal} first bridged statistics literature on constructing matched samples and operations research and computer science literature on matching in graphs and networks. Prior to \citet{rosenbaum1989optimal}, statisticians and practitioners leveraged a ``greedy heuristic" in pair matching that starts with finding a matched pair with a minimum pre-specified distance, removes the pair from further consideration, and iterates the process until finding one control unit for each treated one in the study cohort. Some often-used distance metrics include the Mahalanobis distance in metric-based matching (\citealp{cochran1973controlling,rubin1980bias}), absolute difference in the estimated propensity score in propensity-score matching (\citealp{rosenbaum1983central}), or a combination of both (\citealp{rosenbaum1985constructing}).

\citet{rosenbaum1989optimal}'s key insight is to represent units to be matched, treated and control, as vertices in a bipartite graph and recast the statistical problem of constructing a matched control group as finding a minimum-cost network flow, a standard combinatorial optimization problem for which very efficient algorithms exist. Figure \ref{fig: bipartite schematic plot} illustrates this paradigm using a small example of three treated and five control units. Readers should imagine that the source node on the far left emits three units of flow, one unit arriving at the node $T_1$, one at $T_2$, and the other at $T_3$. Edges of the form $\textsf{e}(T_i, C_j)$ connecting $T_i$ and $C_j$, $i = 1, 2, 3$ and $j = 1, 2,..., 5$, each have capacity $1$ and are associated with a cost $\textsf{cost}\{\textsf{e}(T_i, C_j)\}$ equal to some covariate distance between $T_i$ and $C_j$. These three units of flow arriving at nodes $T_1$, $T_2$, and $T_3$ would choose to flow in such a way that minimizes the total cost; for instance, Figure \ref{fig: bipartite schematic plot} displays a minimum-cost flow (bold, black lines) consisting of $\{\textsf{e}(T_1, C_3), \textsf{e}(T_2, C_1), \textsf{e}(T_3, C_4)\}$. In this way, three matched pairs $\{(T_1, C_3), (T_2, C_1), (T_3, C_4)\}$ are constructed. The network structure displayed in Figure \ref{fig: bipartite schematic plot} can be readily modified to accommodate many useful design aspects, including exact matching, 1-to-$k$ matching, full matching (\citealp{hansen2004full,hansen2006optimal}), subset matching (\citealp{rosenbaum2007minimum}), among others, and is the cornerstone of some most influential statistical matching packages like \textsf{optmatch} (\citealp{hansen2007optmatch}) and \textsf{MatchIt} (\citealp{stuart2011matchit}) in the statistical computing software \textsf{R} (\citealp{team2013r}). Researchers have also proposed network-sparsification techniques to further facilitate statistical matching on large administrative datasets consisting of hundreds of thousands of observations (\citealp{pimentel2015large, yu2020matching}). 

More recently, \citet{zhang2021matching} introduced a minimum-cost network flow algorithm built around a ``tripartite graph," where treated units appear twice, on the far left and far right, with control units sandwiched between them. Efforts to balance high-dimensional covariates (e.g., via the stochastic balancing property of the propensity score or by directly minimizing the earthmover distance between marginal distributions) are represented on the right, while efforts to find close pairings (e.g., exact matching on potential effect modifiers or quantiles of \citet{hansen2008prognostic}'s prognostic score) are represented on the left. By separating two sometimes-conflicting objectives, pairing and balancing, the network built around a tripartite graph is shown to deliver matched samples that are homogeneous in a few key covariates while maintaining good balance in many other covariates (\citealp{zhang2021matching}).


\begin{figure}[ht]
\centering
\begin{tikzpicture}[thick, color = gray,
  every node/.style={draw,circle},
  fsnode/.style={fill=black, inner sep = 0pt, minimum size = 5pt},
  ssnode/.style={fill=gray, inner sep = 0pt, minimum size = 5pt},
  every fit/.style={ellipse,draw,inner sep=-2pt,text width=2cm},
  shorten >= 3pt,shorten <= 3pt
]


\begin{scope}[start chain=going below,node distance=8mm]
\foreach \i in {1,2,3}
  \node[fsnode,on chain] (t\i) [label=above: {\small$T_\i$} ] {};
\end{scope}

\begin{scope}[xshift=3.5cm,yshift=0cm,start chain=going below,node distance=8mm]
\foreach \i in {1,2,3,4,5}
  \node[ssnode,on chain] (c\i) [label=above: {\small$C_\i$}] {};
\end{scope}

\node [fill = black, inner sep = 0pt, minimum size = 5pt, label=left: \textsf{Source}] at (-2, -1) (source) {};

\node [fill = black, inner sep = 0pt, minimum size = 5pt, label=right: \textsf{Sink}] at (5.5, -2) (sink) {};


\draw (t1) -- (c1);
\draw (t1) -- (c2);
\draw [line width = 0.5mm, color = black] (t1) -- (c3);
\draw (t1) -- (c4);
\draw (t1) -- (c5);

\draw [line width = 0.5mm, color = black] (t2) -- (c1);
\draw (t2) -- (c2);
\draw (t2) -- (c3);
\draw (t2) -- (c4);
\draw (t2) -- (c5);

\draw (t3) -- (c1);
\draw (t3) -- (c2);
\draw (t3) -- (c3);
\draw [line width = 0.5mm, color = black] (t3) -- (c4);
\draw (t3) -- (c5);

\draw [line width = 0.5mm, color = black] (source) -- (t1);
\draw [line width = 0.5mm, color = black] (source) -- (t2);
\draw [line width = 0.5mm, color = black] (source) -- (t3);

\draw [line width = 0.5mm, color = black] (c1) -- (sink);
\draw (c2) -- (sink);
\draw [line width = 0.5mm, color = black] (c3) -- (sink);
\draw [line width = 0.5mm, color = black] (c4) -- (sink);
\draw (c5) -- (sink);
\end{tikzpicture}
\caption{Network-flow representation of a small matching problem. Five control units $C_1$, $C_2$, ..., $C_5$ are to be matched to three treated units $T_1$, $T_2$, and $T_3$. Each edge connecting $T_i$, $i = 1,2,3$, and $C_j$, $j = 1,2,3,4,5$, is associated with a flow capacity (equal to $1$ in pair matching) and a cost equal to a pre-specified covariate distance between unit $T_i$ and unit $C_j$.}
\label{fig: bipartite schematic plot}
\end{figure}



\subsection{A Naive Method, An Existing Method, and A New Approach}
One naive approach to create homogeneous matched pairs resembling a target template would first select a subset of treated units from the observational data via matching on covariates collected by both RCT and observational study. In the second step, control units in the observational data are then matched to the treated units selected by the first step. This strategy suffers from a major drawback. Take the RCT DUPLICATE initiative as an example. Although the benchmark SAVOR-TIMI 53 trial reported only around $20$ baseline covariates, \citet{franklin2021emulating} matched on more than $120$ to guard against unmeasured confounding in their emulation with observational data. In the first step, there are many ways to design a smaller treated group similar to the RCT population in $20$ RCT covariates; however, it is difficult to determine which treated group designed from the first step should be utilized to form the eventual match in the second step. A selected treated group similar to the RCT population in $20$ RCT covariates could have poor overlap with the control units in the observational database in the other more than $100$ covariates; see Figure \ref{fig: illustrate naive method} for a toy example that illustrates this phenomenon. Ideally, we prefer well-matched treated and control groups using observational data to first ensure internal validity of the matched cohort study; similarity to some target population should be taken as a perk. The naive approach precludes a principled trade-off between matched cohort study's internal validity and its generalizability to a target population.

\begin{figure}[ht]
    \centering
   \begin{tikzpicture}

\node at (0.25, 1.25) {$X_1$};
\node at (1, 1.25) {$X_2$};
\draw[blue, very thick] (0,0) rectangle (0.5,0.5);
\draw[blue, very thick] (1,0.25) circle (0.25);
\node[text width=2cm, align = center] at (0.6, -4) {\small\textsf{Template}};

\node at (4.25, 1.25) {$X_1$};
\node at (5, 1.25) {$X_2$};
\node at (5.75, 1.25) {$X_3$};

\node at (3.25,0.25) {$\textsf{T1}$};
\node at (3.25,-0.5) {$\textsf{T2}$};
\draw[blue, very thick] (4,0) rectangle (4.5,0.5);
\draw[blue, very thick] (5,0.25) circle (0.25);
\draw[blue, very thick] (5.75,0.25) circle (0.25);

\draw[blue, very thick] (4,-0.75) rectangle (4.5,-0.25);
\draw[blue, very thick] (5,-0.5) circle (0.25);
\draw[blue, very thick] (5.5,-0.75) -- (5.75,-0.25) -- (6,-0.75) -- cycle;

\node[text width=2cm, align = center] at (5, -4) {\small\textsf{OBS \\ Treated}};



\node at (9, 1.25) {$X_1$};
\node at (9.75, 1.25) {$X_2$};
\node at (10.5, 1.25) {$X_3$};


\node at (8,0.25) {$\textsf{C1}$};
\node at (8,-0.5) {$\textsf{C2}$};
\node at (8,-1.25) {$\textsf{C3}$};
\node at (8,-2) {$\textsf{C4}$};

\draw[blue, very thick] (9,0.25) circle (0.25);
\draw[blue, very thick] (9.75,0.25) circle (0.25);
\draw[blue, very thick] (10.5,0.25) circle (0.25);

\draw[blue, very thick] (8.75,-0.75) rectangle (9.25,-0.25);
\draw[blue, very thick] (9.75,-0.5) circle (0.25);
\draw[blue, very thick] (10.5,-0.5) circle (0.25);

\draw[blue, very thick] (8.75,-1.5) rectangle (9.25,-1);
\draw[blue, very thick] (9.75,-1.25) circle (0.25);
\draw[blue, very thick] (10.25,-1.5) rectangle (10.75, -1);

\draw[blue, very thick] (8.75,-2.25) rectangle (9.25,-1.75);
\draw[blue, very thick] (9.5,-2.25) rectangle (10, -1.75);
\draw[blue, very thick] (10.25,-2.25) rectangle (10.75, -1.75);

\node[text width=2cm, align = center] at (9.75, -4) {\small\textsf{OBS \\ Control}};


\end{tikzpicture}
\caption{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 $\textsf{T1}$ and $\textsf{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, $\textsf{T1}$ is a better pick as no control unit has an $X_3$ value similar to that of $\textsf{T2}$.}
\label{fig: illustrate naive method}
\end{figure}

Recently, \citet{bennett2020building} developed a mixed-integer-programming-based (MIP-based) approach that builds representative matched samples. MIP-based matching algorithms (\citealp{li2001balanced,zubizarreta2012using,bennett2020building}) are more flexible compared to network-flow-based ones as requirements for different versions of covariate balance could be readily represented by constraints in a 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), and this is itself an obstacle to some empirical researchers. This being said, if one finds no difficulty installing a powerful solver and that the matching task can be readily solved using MIP-based methods, then this is great and one should do that.

Compared to MIP-based methods, network-flow-based methods 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 (\citealp{rosenbaum2002observational,rosenbaum2010design, stuart2010matching, austin2011introduction,rassen2012one}). The primary goal of this article is to outline an efficient, easy-to-use, 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 explore how the method facilitates different designs of the WHI observational study data and yields some useful, new insight into the well-aged inconsistency between WHI trial and observational study findings. 



\subsection{Additional Design and Computation Considerations}
\label{subsec: additional design considerations}
Many additional design techniques can be used in conjunction with the proposed method by modifying aspects of the basic network structure in Figure \ref{fig: tripartite new scheme}. For instance, it is often helpful to match exactly on the potential effect modifiers in the design stage of an observational study, so that researchers could perform subgroup analysis by doing hypothesis testing in each stratum defined by the effect modifiers; see, e.g., \citet{lee2018discovering,lee2018powerful,lee2021discovering}. If such a design is desired, then researchers should set $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\} = \infty$, or equivalently remove the edge connecting $\overline\tau_t$ and $\gamma_c$, for all $\tau_t$ and $\gamma_c$ disagreeing on the effect modifiers. \citet[Section 3.2]{rosenbaum1989optimal}'s fine balance strategy that forces equi-distribution of one nominal variable in the treated and matched control groups can also be readily accommodated by including an additional ``category" layer in the network; see also \citet{rosenbaum2007minimum, yang2012optimal, pimentel2015large,zhang2021matching}. In some practical situations, certain treated OBS units $\tau_t \in \mathcal{T}_{\textsf{sub}}$ are expected to be included in the matched samples. To achieve this, keep $\text{cost}\{(\tau_t, \overline\tau_t)\} = 0$ for $\tau_t \in \mathcal{T}_{\textsf{sub}}$ and set $\text{cost}\{(\tau_t, \overline\tau_t)\}$ to a large penalty for $\tau_t \in \mathcal{T}\backslash\mathcal{T}_{\textsf{sub}}$.

While the RCT is typically of a smaller, fixed sample size, modern observational databases could easily go beyond $1$ million observations. Consider a scenario with fixed RCT sample size $R$ and growing observational data sample sizes such that $T \rightarrow \infty$, $C \rightarrow \infty$, and $T/C \rightarrow \epsilon \in (0, 1)$. Recall the computation complexity of finding a minnimum cost flow is $O(|\mathcal{V}|\cdot|\mathcal{E}| + |\mathcal{V}|^2\log(|\mathcal{V}|))$. If the network can be made sparser so that $|\mathcal{E}| = O(|\mathcal{V}|\log(|\mathcal{V}|)) = O(C\log(C))$, then the complexity can be reduced from $O(C^3)$ to $O(C^2\log(C))$. One simplest way to sparsify the network is to leverage the propensity score caliper (\citealp{rosenbaum1985constructing,austin2011optimal}) and for each $\tau_t$, keep only edges of the form $e(\overline\tau_t, \gamma_c)$ for $C_1$ control units closest to $\tau_t$ in the estimated propensity score value. In this way, $|\mathcal{E}| = O(C_1 \cdot C)$, and the network is sparsified and computation complexity reduced provided $C_1 = O(\log C)$. Such a strategy is discussed in more detail in \citet{yu2020matching}.

\subsection{A Small Illustrative Example}
\label{subsec: a small example}
\textcolor{red}{Perhaps do a subsection on software and implementation.}
To showcase the utility of the proposed method, we consider a simple example where we generated a target RCT of size $R = 100$, and an observational study with $T = 500$ and $C = 1500$. The RCT collects $d_1 = 5$ covariates while the observational study collects an additional $5$ covariates. The RCT samples $\mathcal{K}$, OBS treated samples $\mathcal{T}$, and OBS control samples $\mathcal{C}$ have the same $\text{Normal}(0, 1)$ distribution for all but the first covariate $X_1$: $X_1 \sim \text{Normal}(0.25, 1)$ in the RCT samples, $X_1 \sim \text{Normal}(1, 1)$ in the OBS treated samples, and $X_1 \sim \text{Normal}(0, 1)$ in the OBS control samples. 

We considered three matched samples. In the first matched comparison, $500$ controls were matched to $500$ OBS treated units without heed of the target RCT covariate distributions. We matched this sample according to two criteria: (i) the earthmover's distance between the distributions of the estimated propensity score in the treated group and matched control group is minimized, and (ii) subject to (i), the total robust Mahalanobis distances between treated and control units are minimized (\citealp{zhang2021matching}). We refer to this match $\textsf{M1}$. In the second matched comparison, we leveraged the proposed new network structure in Figure \ref{fig: tripartite new scheme} with $k = 1$ and $\lambda = 1$ and formed $100$ treated-to-control matched pairs. In this second match, we specified $\delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}})$ as the absolute difference in the estimated generalizability score between $\kappa_r$ and $\tau_t$ in $5$ RCT covariates, and let $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}$ be the robust Mahalanobis distance within a $0.05$ estimated propensity score caliper between $\tau_t$ and $\gamma_c$ in all $10$ covariates. We refer to this match $\textsf{M2}$. The third match $\textsf{M3}$ used the same network structure and distance specifications as $\textsf{M2}$, but with $k = 2$ so that $200$ treated-to-control pairs were formed.

Figure \ref{fig: illustrative example} displays the covariate distributions of three matched samples. Before matching, the treated and control groups are vastly different in $X_1$ but otherwise similar. Match $\textsf{M1}$ does a good job balancing $X_1$ via balancing the estimated propensity score; see the third boxplot in the left panel of Figure \ref{fig: illustrative example}. However, despite within-matched-samples close resemblance, the treated group $\mathcal{T}$ and matched controls produced by $\textsf{M1}$ differ systematically from the target template in the distribution of $X_1$. The $100$ matched pairs returned by $\textsf{M2}$ not only have good within-matched-sample-homogeneity, but also are much more similar to the target template in covariates' distributions. The same pattern holds for the $200$ matched pairs returned by $\textsf{M3}$.

\begin{figure}[ht]
   \centering
     \subfloat{\includegraphics[width = 0.49\columnwidth]{illustrate_X1.pdf}}
     \subfloat{\includegraphics[width = 0.49\columnwidth]{illustrate_X2.pdf}}
     \caption{Comparing the distributions of covariates $X_1$ (left panel) and $X_3$ (right panel) in three matched comparisons. Match $\textsf{M1}$ }
\label{fig: illustrative example}
\end{figure}

\section*{Supplemental Material A: Statistical Matching, Network Flow, Recent Advancement}
Matching and subclassification are widely used in empirical research to embed noisy, non-randomized observational data into an approximate randomized controlled experiment, and facilitate analyzing data as such (\citealp{ho2007matching, rubin2007design,rubin2008objective,rosenbaum2002observational,rosenbaum2010design,rosenbaum2020modern,stuart2010matching,rassen2012one,imbens2015matching}). In a seminal paper, \citet{rosenbaum1989optimal} first bridged statistics literature on constructing matched samples and operations research and computer science literature on matching in graphs and networks. Prior to \citet{rosenbaum1989optimal}, statisticians and practitioners leveraged a ``greedy heuristic" in pair matching that starts with finding a matched pair with a minimum pre-specified distance, removes the pair from further consideration, and iterates the process until finding one control unit for each treated one in the study cohort. Some often-used distance metrics include the Mahalanobis distance in metric-based matching (\citealp{cochran1973controlling,rubin1980bias}), absolute difference in the estimated propensity score in propensity-score matching (\citealp{rosenbaum1983central}), or a combination of both (\citealp{rosenbaum1985constructing}).

\citet{rosenbaum1989optimal}'s key insight is to represent units to be matched, treated and control, as vertices in a bipartite graph and recast the statistical problem of constructing a matched control group as finding a minimum-cost network flow, a standard combinatorial optimization problem for which very efficient algorithms exist. Figure \ref{fig: bipartite schematic plot} illustrates this paradigm using a small example of three treated and five control units. Readers should imagine that the source node on the far left emits three units of flow, one unit arriving at the node $T_1$, one at $T_2$, and the other at $T_3$. Edges of the form $\textsf{e}(T_i, C_j)$ connecting $T_i$ and $C_j$, $i = 1, 2, 3$ and $j = 1, 2,..., 5$, each have capacity $1$ and are associated with a cost $\textsf{cost}\{\textsf{e}(T_i, C_j)\}$ equal to some covariate distance between $T_i$ and $C_j$. These three units of flow arriving at nodes $T_1$, $T_2$, and $T_3$ would choose to flow in such a way that minimizes the total cost; for instance, Figure \ref{fig: bipartite schematic plot} displays a minimum-cost flow (bold, black lines) consisting of $\{\textsf{e}(T_1, C_3), \textsf{e}(T_2, C_1), \textsf{e}(T_3, C_4)\}$. In this way, three matched pairs $\{(T_1, C_3), (T_2, C_1), (T_3, C_4)\}$ are constructed. The network structure displayed in Figure \ref{fig: bipartite schematic plot} can be readily modified to accommodate many useful design aspects, including exact matching, 1-to-$k$ matching, full matching (\citealp{hansen2004full,hansen2006optimal}), subset matching (\citealp{rosenbaum2007minimum}), among others, and is the cornerstone of some most influential statistical matching packages like \textsf{optmatch} (\citealp{hansen2007optmatch}) and \textsf{MatchIt} (\citealp{stuart2011matchit}) in the statistical computing software \textsf{R} (\citealp{team2013r}). Researchers have also proposed network-sparsification techniques to further facilitate statistical matching on large administrative datasets consisting of hundreds of thousands of observations (\citealp{pimentel2015large, yu2020matching}). 

More recently, \citet{zhang2021matching} introduced a minimum-cost network flow algorithm built around a ``tripartite graph," where treated units appear twice, on the far left and far right, with control units sandwiched between them. Efforts to balance high-dimensional covariates (e.g., via the stochastic balancing property of the propensity score or by directly minimizing the earthmover distance between marginal distributions) are represented on the right, while efforts to find close pairings (e.g., exact matching on potential effect modifiers or quantiles of \citet{hansen2008prognostic}'s prognostic score) are represented on the left. By separating two sometimes-conflicting objectives, pairing and balancing, the network built around a tripartite graph is shown to deliver matched samples that are homogeneous in a few key covariates while maintaining good balance in many other covariates (\citealp{zhang2021matching}).


\begin{figure}[ht]
\centering
\begin{tikzpicture}[thick, color = gray,
  every node/.style={draw,circle},
  fsnode/.style={fill=black, inner sep = 0pt, minimum size = 5pt},
  ssnode/.style={fill=gray, inner sep = 0pt, minimum size = 5pt},
  every fit/.style={ellipse,draw,inner sep=-2pt,text width=2cm},
  shorten >= 3pt,shorten <= 3pt
]


\begin{scope}[start chain=going below,node distance=8mm]
\foreach \i in {1,2,3}
  \node[fsnode,on chain] (t\i) [label=above: {\small$T_\i$} ] {};
\end{scope}

\begin{scope}[xshift=3.5cm,yshift=0cm,start chain=going below,node distance=8mm]
\foreach \i in {1,2,3,4,5}
  \node[ssnode,on chain] (c\i) [label=above: {\small$C_\i$}] {};
\end{scope}

\node [fill = black, inner sep = 0pt, minimum size = 5pt, label=left: \textsf{Source}] at (-2, -1) (source) {};

\node [fill = black, inner sep = 0pt, minimum size = 5pt, label=right: \textsf{Sink}] at (5.5, -2) (sink) {};


\draw (t1) -- (c1);
\draw (t1) -- (c2);
\draw [line width = 0.5mm, color = black] (t1) -- (c3);
\draw (t1) -- (c4);
\draw (t1) -- (c5);

\draw [line width = 0.5mm, color = black] (t2) -- (c1);
\draw (t2) -- (c2);
\draw (t2) -- (c3);
\draw (t2) -- (c4);
\draw (t2) -- (c5);

\draw (t3) -- (c1);
\draw (t3) -- (c2);
\draw (t3) -- (c3);
\draw [line width = 0.5mm, color = black] (t3) -- (c4);
\draw (t3) -- (c5);

\draw [line width = 0.5mm, color = black] (source) -- (t1);
\draw [line width = 0.5mm, color = black] (source) -- (t2);
\draw [line width = 0.5mm, color = black] (source) -- (t3);

\draw [line width = 0.5mm, color = black] (c1) -- (sink);
\draw (c2) -- (sink);
\draw [line width = 0.5mm, color = black] (c3) -- (sink);
\draw [line width = 0.5mm, color = black] (c4) -- (sink);
\draw (c5) -- (sink);
\end{tikzpicture}
\caption{Network-flow representation of a small matching problem. Five control units $C_1$, $C_2$, ..., $C_5$ are to be matched to three treated units $T_1$, $T_2$, and $T_3$. Each edge connecting $T_i$, $i = 1,2,3$, and $C_j$, $j = 1,2,3,4,5$, is associated with a flow capacity (equal to $1$ in pair matching) and a cost equal to a pre-specified covariate distance between unit $T_i$ and unit $C_j$.}
\label{fig: bipartite schematic plot}
\end{figure}

Representing statistical matching as a two-part network as depicted in Figure \ref{fig: bipartite schematic plot} is an attractive, general conceptual framework. Recent advancements in network-flow-based matching algorithms often involve adding an additional ``category" layer to the two-part network and using this additional layer to balance one or more nominal covariates in the treated and matched control groups (\citealp{rosenbaum1989optimal,yang2012optimal, pimentel2015large}). More recently, \citet{zhang2021matching} introduced a minimum-cost network flow algorithm built around a ``tripartite graph," where treated units appear twice, on the far left and far right, with control units sandwiched between them. Efforts to balance high-dimensional covariates (e.g., via the stochastic balancing property of the propensity score or by directly minimizing the earthmover's distance between marginal distributions) are represented on the right, while efforts to find close pairings (e.g., exact matching on potential effect modifiers or quantiles of \citet{hansen2008prognostic}'s prognostic score) are represented on the left. By separating two sometimes-conflicting objectives, pairing and balancing, the network built around a tripartite graph is shown to deliver matched samples that are homogeneous in a few key covariates while maintaining good balance in many other covariates (\citealp{zhang2021matching}).


\section*{Supplemental Material B: Additional Design and Computation Considerations}
Many additional design techniques can be used in conjunction with the proposed method by modifying aspects of the basic network structure in Figure \ref{fig: tripartite new scheme}. For instance, it is often helpful to match exactly on the potential effect modifiers in the design stage of an observational study, so that researchers could perform subgroup analysis by doing hypothesis testing in each stratum defined by the effect modifiers; see, e.g., \citet{lee2018discovering,lee2018powerful,lee2021discovering}. If such a design is desired, then researchers should set $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\} = \infty$, or equivalently remove the edge connecting $\overline\tau_t$ and $\gamma_c$, for all $\tau_t$ and $\gamma_c$ disagreeing on the effect modifiers. \citet[Section 3.2]{rosenbaum1989optimal}'s fine balance strategy that forces equi-distribution of one nominal variable in the treated and matched control groups can also be readily accommodated by including an additional ``category" layer in the network; see also \citet{rosenbaum2007minimum, yang2012optimal, pimentel2015large,zhang2021matching}. In some practical situations, certain treated OBS units $\tau_t \in \mathcal{T}_{\textsf{sub}}$ are expected to be included in the matched samples. To achieve this, keep $\text{cost}\{(\tau_t, \overline\tau_t)\} = 0$ for $\tau_t \in \mathcal{T}_{\textsf{sub}}$ and set $\text{cost}\{(\tau_t, \overline\tau_t)\}$ to a large penalty for $\tau_t \in \mathcal{T}\backslash\mathcal{T}_{\textsf{sub}}$.

While the RCT is typically of a smaller, fixed sample size, modern observational databases could easily go beyond $1$ million observations. Consider a scenario with fixed RCT sample size $R$ and growing observational data sample sizes such that $T \rightarrow \infty$, $C \rightarrow \infty$, and $T/C \rightarrow \epsilon \in (0, 1)$. Recall the computation complexity of finding a minnimum cost flow is $O(|\mathcal{V}|\cdot|\mathcal{E}| + |\mathcal{V}|^2\log(|\mathcal{V}|))$. If the network can be made sparser so that $|\mathcal{E}| = O(|\mathcal{V}|\log(|\mathcal{V}|)) = O(C\log(C))$, then the complexity can be reduced from $O(C^3)$ to $O(C^2\log(C))$. One simplest way to sparsify the network is to leverage the propensity score caliper (\citealp{rosenbaum1985constructing,austin2011optimal}) and for each $\tau_t$, keep only edges of the form $e(\overline\tau_t, \gamma_c)$ for $C_1$ control units closest to $\tau_t$ in the estimated propensity score value. In this way, $|\mathcal{E}| = O(C_1 \cdot C)$, and the network is sparsified and computation complexity reduced provided $C_1 = O(\log C)$. Such a strategy is discussed in more detail in \citet{yu2020matching}.

\section*{Supplemental Material C: A Small Illustrative Example}
To showcase the utility of the proposed method, we consider a simple example where we generated a target RCT of size $R = 100$, and an observational study with $T = 500$ and $C = 1500$. The RCT collects $d_1 = 5$ covariates while the observational study collects an additional $5$ covariates. The RCT samples $\mathcal{K}$, OBS treated samples $\mathcal{T}$, and OBS control samples $\mathcal{C}$ have the same $\text{Normal}(0, 1)$ distribution for all but the first covariate $X_1$: $X_1 \sim \text{Normal}(0.25, 1)$ in the RCT samples, $X_1 \sim \text{Normal}(1, 1)$ in the OBS treated samples, and $X_1 \sim \text{Normal}(0, 1)$ in the OBS control samples. 

We considered three matched samples. In the first matched comparison, $500$ controls were matched to $500$ OBS treated units without heed of the target RCT covariate distributions. We matched this sample according to two criteria: (i) the earthmover's distance between the distributions of the estimated propensity score in the treated group and matched control group is minimized, and (ii) subject to (i), the total robust Mahalanobis distances between treated and control units are minimized (\citealp{zhang2021matching}). We refer to this match $\textsf{M1}$. In the second matched comparison, we leveraged the proposed new network structure in Figure \ref{fig: tripartite new scheme} with $k = 1$ and $\lambda = 1$ and formed $100$ treated-to-control matched pairs. In this second match, we specified $\delta_{\kappa_r, \tau_t} (\widetilde{\boldsymbol{x}})$ as the absolute difference in the estimated generalizability score between $\kappa_r$ and $\tau_t$ in $5$ RCT covariates, and let $\Delta_{\tau_t, \gamma_c}\{(\widetilde{\boldsymbol{x}}, \boldsymbol{x})\}$ be the robust Mahalanobis distance within a $0.05$ estimated propensity score caliper between $\tau_t$ and $\gamma_c$ in all $10$ covariates. We refer to this match $\textsf{M2}$. The third match $\textsf{M3}$ used the same network structure and distance specifications as $\textsf{M2}$, but with $k = 2$ so that $200$ treated-to-control pairs were formed.

Figure \ref{fig: illustrative example} displays the covariate distributions of three matched samples. Before matching, the treated and control groups are vastly different in $X_1$ but otherwise similar. Match $\textsf{M1}$ does a good job balancing $X_1$ via balancing the estimated propensity score; see the third boxplot in the left panel of Figure \ref{fig: illustrative example}. However, despite within-matched-samples close resemblance, the treated group $\mathcal{T}$ and matched controls produced by $\textsf{M1}$ differ systematically from the target template in the distribution of $X_1$. The $100$ matched pairs returned by $\textsf{M2}$ not only have good within-matched-sample-homogeneity, but also are much more similar to the target template in covariates' distributions. The same pattern holds for the $200$ matched pairs returned by $\textsf{M3}$.

\begin{figure}[ht]
   \centering
     \subfloat{\includegraphics[width = 0.49\columnwidth]{illustrate_X1.pdf}}
     \subfloat{\includegraphics[width = 0.49\columnwidth]{illustrate_X2.pdf}}
     \caption{Comparing the distributions of covariates $X_1$ (left panel) and $X_3$ (right panel) in three matched comparisons.}
\label{fig: illustrative example}
\end{figure}



\section*{Supplemental Material D: Simulation Studies}
\subsection*{D.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 performance of downstream statistical inference. 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 template. We generated a template $\mathcal{K}$ with $|\mathcal{K}| = 300$ and $\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}$. We consider an observational database with $|\mathcal{T}| = 1000$ treated and $|\mathcal{C}| = 3000$ 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:
\begin{description}
\item \textbf{Factor 1:} Dimension of covariates in the observational database, $d$: $10$, $30$, and $50$.
\item \textbf{Factor 2:} Overlap between template, treated ($Z = 1$), and control ($Z = 0$) units, $\beta$: 
$\mathbf{X} \sim \text{Multivariate Normal}\left(\mathbf\mu, \mathbf \Sigma\right)$, with $\mathbf \mu = (\theta Z, 0, 0, 0, 0)^{\text{T}}$ and $\mathbf\Sigma = \mathbf I_{d\times d}$. We consider $\theta = 0.25$, $0.50$, and $1$.

\item \textbf{Factor 3:} Matching algorithms to be investigated, $\mathcal{M}$: 
\begin{enumerate}
    \item $\mathcal{M}_{\textsf{opt}}$: matching according to two criteria (\citealp{zhang2021matching}): (i) minimizing the earthmover's distance between the distributions of the estimated propensity score in the treated and matched control groups, and (ii) subject to (i), minimizing the within-matched-pair robust Mahalanobis distances. Algorithm $\mathcal{M}_{\textsf{opt}}$ produces $1000$ matched pairs, and pays no attention to the target template.
    
    \item $\mathcal{M}_{\textsf{template},~k = 1}$: matching according to the proposed network structure (Figure \ref{fig: tripartite new scheme}) with $k = 1$. Algorithm $\mathcal{M}_{\textsf{template},~k = 1}$ produces $300$ matched pairs.
   
    \item $\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.
\end{enumerate}

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

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 under consideration. As discussed in Section \ref{subsec: tradeoff}, the tuning parameter $\lambda$ controls the trade-off between internal and external validity of effect estimates obtained from matched samples. 

For each unit, we further generate two potential outcomes according to the following data-generating process:
\begin{equation}
    \label{eq: simu outcome}
    Y(0) \sim N(0, 1),\qquad Y(1) = Y(0) + \beta(X_1),
\end{equation}
and the observed outcome satisfies $Y = Z\cdot Y(1) + (1-Z)\cdot Y(0)$. The last factor specifies the treatment effect:
\begin{description}
\item \textbf{Factor 5:} Treatment effect, $\beta(X_1)$: a constant treatment effect $\beta(X_1) = 2$, a mildly heterogeneous treatment effect $\beta(X_1) = 2 - 0.2X_1$, and a strongly heterogeneous treatment effect $\beta(X_1) = 2 - X_1$. 
\end{description}
When the treatment effect is heterogeneous and the effect modifier $X_1$ has a different distribution in the RCT group and the OBS treated group, then the average treated effect on the treated (ATT) estimate obtained from the matched observational study cannot be automatically generalized to the RCT cohort (\citealp{stuart2011use}). 

There are multiple ways to analyze matched pair data. Examples include parametric t-test, randomization inference (\citealp{rosenbaum2002observational,rosenbaum2010design}), and regression adjustment (\citealp{rubin1979using, ho2007matching}). In this simulation study, we report a simple difference-in-means estimator for matched data produced by each of the $7$ algorithms, and compare these $7$ effect estimates to the treatment effect averaged over the target template, i.e., $\textsf{SATE}_{\text{target}}$. We report the average bias and mean squared error of each of the $7$ effect estimates against $\textsf{SATE}_{\text{target}}$.

\subsection*{D.2: Simulation Results}
Table \ref{tbl: simulation results b = 1} summarizes the percentage of bias with respect to $\textsf{SATE}_{\text{target}}$ of each difference-in-means estimator $\widehat{\theta}_{\mathcal{M}}$ obtained from matched samples constructed from algorithm $\mathcal{M}$. 


We have observed a few trends consistent with both theory and intuition. First, when the treatment is constant, there is no generalizability issue and the bias is minuscule 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 treated units, effect estimate obtained from $1000$ matched pairs constructed using algorithm $\mathcal{M}_{\textsf{opt}}$ is clearly biased from $\textsf{SATE}_{\text{target}}$, and the percentage of bias increases (i) as the distributions of $X_1$ become increasingly dissimilar in the target template and in the OBS treated units (i.e., as $b$ increases), and (ii) as effect modification becomes more dramatic (i.e., from constant to mild to strong). In the most adversarial setting considered in this simulation study, i.e., when $b = 1$ and $\beta(X_1) = 2 - X_1$, the percentage bias of $\widehat{\theta}_{\mathcal{M}_{\textsf{opt}}}$ can be as large as $45\%$. We need to stress that although $\widehat{\theta}_{\mathcal{M}_{\textsf{opt}}}$ may not be generalized to the target template, it is a perfectly internally-valid estimator for the treatment effect averaged over covariates' distributions of OBS treated units.

Our proposed template matching algorithms outperform $\mathcal{M}_{\textsf{opt}}$ in bias reduction against $\textsf{SATE}_{\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 large so that the matching algorithm gives more priority to the left part compared to the right part of the network depicted in Figure \ref{fig: tripartite new scheme}, i.e., resemblance to the target template in covariates' distributions is emphasized over resemblance between matched treated and matched control groups. Both parameters, $k$ and $\lambda$, effectively allow a trade-off between the internal validity of a matched observational study and its generalizability to a target population. 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., \cite{franklin2020nonrandomized, franklin2021emulating}. Provided that no outcome data are viewed, researchers can feel free to perform multiple statistical matching and select the one achieving the best trade-off between internal and external validity in an real world problem.

\begin{table}[ht]
\centering
\caption{Percentage of bias with respect to $\textsf{SATE}_{\text{target}}$ of $7$ difference-in-means estimator constructed from matched samples obtained from each of the $7$ matching algorithms. Each cell is averaged over $1000$ simulations.}
\label{tbl: simulation results b = 1}
\resizebox{0.9\textwidth}{!}{
\begin{tabular}{cccccccccccc}\hline \\ [-0.8em]\multirow{3}{*}{\begin{tabular}{c}Heterogeneity \\ Level\end{tabular}}
&\multirow{3}{*}{\begin{tabular}{c}$b$\end{tabular}}&& \multirow{3}{*}{\begin{tabular}{c}$\mathcal{M}_{\textsf{opt}}$\end{tabular}} 
&& \multicolumn{6}{c}{$\mathcal{M}_{\textsf{template}}$} \\ \\ [-0.8em] \cline{6-11} \\ [-0.8em]
&&&& & \multirow{2}{*}{\begin{tabular}{c}$k = 1$ \\ $\lambda = 0.01$ \end{tabular}}  
& \multirow{2}{*}{\begin{tabular}{c}$k = 1$ \\ $\lambda = 1$ \end{tabular}} 
& \multirow{2}{*}{\begin{tabular}{c}$k = 1$ \\ $\lambda = 100$ \end{tabular}} 
& \multirow{2}{*}{\begin{tabular}{c}$k = 2$ \\ $\lambda = 0.01$ \end{tabular}}  
& \multirow{2}{*}{\begin{tabular}{c}$k = 2$ \\ $\lambda = 1$ \end{tabular}} 
& \multirow{2}{*}{\begin{tabular}{c}$k = 2$ \\ $\lambda = 100$ \end{tabular}} \\ \\ \\ [-0.8em]
    \multicolumn{11}{c}{$d = 10$} \\ \\ [-0.8em]
    
   
    \multirow{3}{*}{Constant} &0.25
    && 0.001 && 0.000 & 0.001 & 0.001 & 0.001 & 0.000 & 0.001 \\ 
    &0.50 && 0.000 && 0.001 & 0.001 & 0.001 & 0.001 & 0.000 & 0.000 \\ 
    & 1.00 && 0.000 && 0.002 & 0.000 & 0.001 & 0.000 & 0.000 & 0.000 \\ \\ [-0.8em]

       \multirow{3}{*}{Mild} &0.25
   && 0.007 && 0.003 & 0.000 & 0.002 & 0.004 & 0.001 & 0.002 \\ 
  &0.50 && 0.033 && 0.009 & 0.008 & 0.018 & 0.011 & 0.012 & 0.020 \\ 
  &1.00 && 0.084 && 0.019 & 0.020 & 0.044 & 0.034 & 0.034 & 0.050 \\ \\ [-0.8em]


  \multirow{3}{*}{Strong} &0.25
  && 0.042 && 0.007 & 0.003 & 0.001 & 0.011 & 0.006 & 0.010 \\ 
 &0.50 && 0.177 && 0.041 & 0.035 & 0.086 & 0.055 & 0.050 & 0.103 \\ 
 &1.00 && 0.451 && 0.102 & 0.098 & 0.223 & 0.188 & 0.186 & 0.262 \\ \\ [-0.8em]
     \multicolumn{11}{c}{$d = 30$} \\  \\ [-0.8em]
        
          \multirow{3}{*}{Constant} &0.25
    && 0.000 && 0.003 & 0.000 & 0.001 & 0.000 & 0.000 & 0.000 \\ 
    &0.50 && 0.001 && 0.001 & 0.001 & 0.002 & 0.001 & 0.002 & 0.002 \\ 
    &1.00 && 0.001 && 0.001 & 0.000 & 0.001 & 0.000 & 0.001 & 0.002 \\ \\ [-0.8em]
   
     \multirow{3}{*}{Mild} &0.25
    && 0.008 && 0.003 & 0.003 & 0.004 & 0.003 & 0.003 & 0.005 \\ 
    &0.50 && 0.031 && 0.008 & 0.009 & 0.015 & 0.009 & 0.013 & 0.019 \\ 
    &1.00 && 0.084 && 0.019 & 0.025 & 0.037 & 0.036 & 0.039 & 0.048 \\ \\ [-0.8em]
 
   \multirow{3}{*}{Strong} &0.25
  && 0.042 && 0.007 & 0.007 & 0.012 & 0.010 & 0.012 & 0.020 \\ 
 &0.50 && 0.178 && 0.040 & 0.058 & 0.095 & 0.053 & 0.080 & 0.114 \\ 
 &1.00 && 0.452 && 0.105 & 0.137 & 0.207 & 0.189 & 0.214 & 0.259 \\ \\ [-0.8em]
  \multicolumn{11}{c}{$d = 50$} \\  \\ [-0.8em]
 
   \multirow{3}{*}{Constant} &0.25
  && 0.001 && 0.000 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 \\ 
 &0.50 && 0.001 && 0.002 & 0.001 & 0.001 & 0.002 & 0.000 & 0.000 \\ 
 &1.00 && 0.000 && 0.000 & 0.001 & 0.000 & 0.001 & 0.001 & 0.001 \\ \\ [-0.8em]
 
   \multirow{3}{*}{Mild} &0.25
  && 0.006 && 0.001 & 0.001 & 0.001 & 0.002 & 0.001 & 0.003 \\ 
 &0.50 && 0.034 && 0.008 & 0.011 & 0.017 & 0.010 & 0.016 & 0.021 \\ 
 &1.00 && 0.084 && 0.019 & 0.028 & 0.038 & 0.035 & 0.042 & 0.049 \\ \\ [-0.8em]
 
   \multirow{3}{*}{Strong} &0.25
  && 0.040 && 0.007 & 0.006 & 0.013 & 0.009 & 0.014 & 0.021 \\ 
 &0.50 && 0.178 && 0.038 & 0.065 & 0.095 & 0.053 & 0.088 & 0.113 \\ 
 &1.00 && 0.451 && 0.102 & 0.150 & 0.205 & 0.188 & 0.224 & 0.260 \\ \\ [-0.8em] \hline
\end{tabular}}
\end{table}


\clearpage
\section*{Supplemental Material E: More Details on the WHI Case Study}
\subsection*{E.1: Plot of the Estimated Propensity Score}

\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{all_pscore_plot.pdf}
    \caption{Distributions of the estimated propensity score in the entire treated OBS group (blue), the entire control OBS group (red), and the matched control group (green).}
\end{figure}


\subsection*{E.2: Balance table of $\textsf{M0}$}

\begin{table}[H]
\centering
\caption{\small Balance table before matching and of $18,340$ matched pairs in the match $\textsf{M0}$}
\label{tbl: matched samples I}
\resizebox{0.78\textwidth}{!}{
\begin{tabular}{lccccccccc}
  \hline
  &\multirow{4}{*}{\begin{tabular}{c}\textbf{OBS} \\ \textbf{Treated} \\ \textbf{(n = 18340)}\end{tabular}} 
    &\multirow{4}{*}{\begin{tabular}{c}\textbf{OBS} \\ \textbf{Control} \\ \textbf{(n = 75303)}\end{tabular}}
        &\multirow{4}{*}{\begin{tabular}{c}\textbf{SMD} \\ \textbf{Before}\end{tabular}}
     &\multirow{4}{*}{\begin{tabular}{c}\textbf{Matched} \\ \textbf{Control} \\\textbf{(n=18340)}\end{tabular}}
   &\multirow{4}{*}{\begin{tabular}{c}\textbf{SMD} \\ \textbf{M0}\end{tabular}}
   \\ \\ \\ \\ \\ [-0.6em]
   
   \multicolumn{6}{c}{\large\textbf{HRT Prior Usage}}\\  \\ [-0.6em]
   \textbf{Age at initiation} \\
   \hspace{0.2cm}NA & 0 &0.89 & &0.85\\
   \hspace{0.2cm}Age & 53.87 &52.04 & & 51.75\\
   \textbf{Previous E+P use, yrs} & 7.02 & 0.47 & &0.55 &  \\ \\ [-0.6em]
   
   \multicolumn{6}{c}{\large\textbf{Covariates Collected in RCT and OBS}}\\  \\ [-0.6em]
\textbf{Age at screening} & 60.84 & 64.30 & -0.35 & 61.10 & -0.03 \\ 
\textbf{Race/Ethnicity}\\
 \hspace{0.2cm}White & 0.89 & 0.82 & 0.14 & 0.89 & -0.00 \\ 
   \hspace{0.2cm}Black/Hispanic & 0.06 & 0.14 & -0.18 & 0.06 & 0.00 \\  
   \textbf{Education}\\
 \hspace{0.2cm}College or above & 0.55 & 0.40 & 0.21 & 0.55 & 0.00 \\ 
  \hspace{0.2cm} Some college & 0.32 & 0.37 & -0.08 & 0.32 & 0.00 \\ 
  \hspace{0.2cm} High school diploma/GED & 0.11 & 0.17 & -0.13 & 0.11 & -0.00 \\ 
  \textbf{Blood pressure}\\
  \hspace{0.2cm}Systolic & 123.27 & 127.86 & -0.19 & 123.54 & -0.01 \\ 
  \hspace{0.2cm}Diastolic & 74.03 & 74.91 & -0.07 & 74.23 & -0.02 \\
  \textbf{BMI} & 25.85 & 27.61 & -0.22 & 25.99 & -0.02 \\
 \textbf{Smoking}\\
 \hspace{0.2cm}Current smoker & 0.05 & 0.06 & -0.04 & 0.05 & 0.00 \\ 
  \hspace{0.2cm}Never smoked & 0.47 & 0.51 & -0.05 & 0.48 & -0.00 \\  
  \hspace{0.2cm}Past smoker & 0.46 & 0.41 & 0.07 & 0.46 & 0.00 \\ 
  \textbf{No. of PA episodes}\\
  \hspace{0.2cm}Total & 5.80 & 5.29 & 0.09 & 5.74 & 0.01 \\ 
  \hspace{0.2cm}Medium to strenuous & 3.58 & 2.94 & 0.13 & 3.50 & 0.02 \\ \\
   \multicolumn{6}{c}{\large\textbf{Additional OBS Covariates}}\\  \\ [-0.6em]
  \textbf{Region}\\
  \hspace{0.2cm}Midwest & 0.22 & 0.22 & 0.00 & 0.23 & -0.01 \\ 
  \hspace{0.2cm}Northeast & 0.18 & 0.24 & -0.09 & 0.19 & -0.02 \\
  \hspace{0.2cm}South & 0.25 & 0.26 & -0.02 & 0.25 & 0.00 \\ 
     \textbf{Partner's education}\\
  \hspace{0.2cm}College or above & 0.43 & 0.30 & 0.20 & 0.42 & 0.02 \\ 
  \hspace{0.2cm}Some college & 0.16 & 0.17 & -0.01 & 0.17 & -0.01 \\  
   \hspace{0.2cm}High school diploma/GED & 0.07 & 0.09 & -0.05 & 0.07 & -0.00 \\ 
   \textbf{Income}\\
   \hspace{0.2cm}Below 35K & 0.23 & 0.40 & -0.26 & 0.24 & -0.01 \\ 
   \hspace{0.2cm} 35K - 75K & 0.42 & 0.36 & 0.09 & 0.44 & -0.02 \\  
   \hspace{0.2cm}Above 75K & 0.29 & 0.16 & 0.22 & 0.27 & 0.04 \\ 
   \textbf{Marital status}\\
   \hspace{0.2cm}Married & 0.69 & 0.60 & 0.13 & 0.69 & 0.01 \\
   \hspace{0.2cm}Divorced/Widowed & 0.26 & 0.35 & -0.13 & 0.27 & -0.01 \\ 
 \textbf{Employment status}\\
  \hspace{0.2cm}Yes & 0.45 & 0.32 & 0.19 & 0.45 & -0.00 \\ 
  \hspace{0.2cm}No & 0.53 & 0.65 & -0.17 & 0.53 & 0.00 \\  
\textbf{Reproductive history}\\  
  \hspace{0.2cm}No ovary removed & 0.92 & 0.64 & 0.51 & 0.92 & 0.00 \\ 
\hspace{0.2cm}Oral contraceptive use ever & 0.53 & 0.37 & 0.24 & 0.53 & 0.00 \\ 
  \hspace{0.2cm}OC duration in years & 5.59 & 5.16 & 0.09 & 5.41 & 0.04 \\
  \textbf{Preexisting Conditions}\\
  \hspace{0.2cm}Stroke & 0.01 & 0.02 & -0.08 & 0.00 & 0.01 \\ 
  \hspace{0.2cm}MI  & 0.01 & 0.03 & -0.07 & 0.01 & 0.00 \\ 
  \hspace{0.2cm}CHF & 0.00 & 0.01 & -0.06 & 0.00 & 0.00 \\
  \hspace{0.2cm}Liver diseases & 0.02 & 0.02 & -0.02 & 0.02 & -0.01 \\
  \hspace{0.2cm}Hypertension & 0.25 & 0.35 & -0.16 & 0.25 & -0.00 \\ 
  \hspace{0.2cm}Fracture & 0.10 & 0.14 & -0.09 & 0.11 & -0.01 \\ 
  \hspace{0.2cm}CABG/PTCA & 0.01 & 0.02 & -0.05 & 0.01 & 0.00 \\
  \hspace{0.2cm}BRCA & 0.01 & 0.07 & -0.22 & 0.01 & 0.00 \\ 
   \hline
\end{tabular}}
\end{table}

