#  Manipulating a Continuous Instrumental Variable in an Observational Study
of Premature Babies:
Algorithm, Partial Identification Bounds, and Inference under Randomization
and Biased Randomization Assumptions

Zhe Chen † Department of Statistics, University of Illinois Urbana-Champaign
Min Haeng Cho † Department of Biostatistics, University of Washington  Bo
Zhang Assistant Professor of Biostatistics, Vaccine and Infectious Disease
Division, Fred Hutchinson Cancer Center. Email<EMAIL_ADDRESS>† The
first two authors contributed equally. Vaccine and Infectious Disease
Division, Fred Hutchinson Cancer Center

Abstract: Regionalization of intensive care for premature babies refers to a
triage system of mothers with high-risk pregnancies to hospitals of varied
capabilities based on risks faced by infants. Due to the limited capacity of
high-level hospitals, which are equipped with advanced expertise to provide
critical care, understanding the effect of delivering premature babies at such
hospitals on infant mortality for different subgroups of high-risk mothers
could facilitate the design of an efficient perinatal regionalization system.
Towards answering this question, Baiocchi et al. (2010) proposed to
_strengthen_ an excess-travel-time-based, continuous instrumental variable
(IV) in an IV-based, matched-pair design by switching focus to a smaller
cohort amenable to being paired with a larger separation in the IV dose. Three
elements changed with the strengthened IV: the study cohort, compliance rate
and latent complier subgroup. Here, we introduce a non-bipartite, template
matching algorithm that embeds data into a target, pair-randomized
encouragement trial which maintains fidelity to the original study cohort
while strengthening the IV. We then study randomization-based and IV-
dependent, biased-randomization-based inference of partial identification
bounds for the sample average treatment effect (SATE) in an IV-based matched
pair design, which deviates from the usual effect ratio estimand in that the
SATE is agnostic to the IV and who is matched to whom, although a strengthened
IV design could narrow the partial identification bounds. Based on our
proposed strengthened-IV design, we found that delivering at a high-level NICU
reduced preterm babies’ mortality rate compared to a low-level NICU for
$81,766\times 2=163,532$ mothers and their preterm babies and the effect
appeared to be minimal among non-black, low-risk mothers.

Keywords: Biased randomization; Heterogeneous treatment effect; Instrumental
variable; Partial identification bound; Target trial emulation

##  1 Introduction

###  1.1 Does delivery at a high-level NICU reduce infant mortality?

Regionalization of intensive care for premature babies refers to a triage
system of mothers with high-risk pregnancies to hospitals of varied
capabilities based on risks faced by infants. In the analyses of premature
births in Pennsylvania, Baiocchi et al., (2010) and Lorch et al., (2012) found
that delivering premature babies at a high-level neonatal intensive care unit
(NICU), broadly defined as a medical unit equipped with advanced technical
expertise and capacity to provide critical care for newborns and premature
infants, was associated with lower infant mortality rate for a latent complier
subgroup. Yet, due to the limited capacity of high-level NICUs, it is not
practical to triage all mothers with high-risk pregnancies to high-level
NICUs. Some previous works have further explored the treatment effect
heterogeneity of delivering at a high-level versus low-level NICU. For
instance, Yang et al., (2014) found that delivering at a high-level NICU
significantly reduced deaths for babies of small gestational age but had a
negligible effect for almost mature babies. More recently, Chen and Zhang,
(2023) estimated individual treatment rules subject to different levels of
resource constraints and found that mothers’ age and race/ethnicity may also
be potential effect modifiers. These findings appear to align well with the
clinical literature. According to a $2015$ National Vital Statistics Reports
(Mathews et al.,, 2015), black infants in the US had a $2.2$-fold greater
mortality rate than white infants. On the other hand, it has long been
established that the risk of miscarriage increases with age; in particular,
Hansen, (1986) found that pregnant women aged $35$ or older experience an
increased risk of intrauterine fetal death, pregnancy-induced hypertension,
gestational diabetes, and delivery by cesarean. Motivated by these findings,
we utilize observational data from the Commonwealth of Pennsylvania to
investigate whether delivering preterm babies at a high-level NICU reduces
infant mortality for all mothers in Pennsylvania and certain high-risk
subgroups like black mothers and mother who had reached or were past the
advanced maternal age of $35$.

###  1.2 Excess travel time as an instrumental variable

To estimate the treatment effect from retrospective observational data, a
naïve study design would compare mothers who received care at a high-level
NICU with those at a low-level NICU. However, as pointed out in numerous
previous works (Baiocchi et al.,, 2010; Lorch et al.,, 2012; Yang et al.,,
2014; Pu and Zhang,, 2021; Chen and Zhang,, 2023), treatment effect estimates
derived from such a naïve comparison would be biased because receipt of
neonatal intensive care at a high-level NICU as opposed to a low-level one
could be easily confounded by unmeasured variables not captured in
administrative databases.

When faced with unmeasured confounding bias, researchers often resort to
quasi-experimental devices (see, e.g., Cook et al.,, 2002; Rosenbaum,, 2010).
Instrumental variable (IV) methods are among the most popular (Haavelmo,,
1943; Angrist and Imbens,, 1995; Angrist et al.,, 1996). Roughly speaking, a
valid IV is a variable independent of unmeasured treatment-outcome
confounders, possibly within strata defined by observed covariates, and whose
sole effect on the outcome is through its effect on the treatment (see Section
4 for a formal definition in the context of a continuous IV). In a seminal
paper, Angrist et al., (1996) showed that under an additional “monotonicity”
or “no defiers” assumption, a valid IV could be used to nonparametrically
identify the treatment effect among a well-defined, albeit latent, subgroup
referred to as “compliers.”

In their original study design, Baiocchi et al., (2010) used excess travel
time, or the difference in travel times (in minutes) from a mother’s zip code
to the nearest high-level and low-level NICUs, as an IV for delivering at a
high-level versus low-level NICU. According to this definition, a small excess
travel time indicates proximity to a high-level NICU relative to a low-level
NICU and corresponds to a strong encouragement to deliver at a high-level
NICU. Similar proximity-based IVs have been extensively used in health
services research (Newhouse and McClellan,, 1998; Baiocchi et al.,, 2014).
Using an IV-based method, Baiocchi et al., (2010) concluded that for every
$100$ additional mothers encouraged by the IV to deliver at a high-level NICU,
$0.9$ additional infant deaths could be avoided.

###  1.3 Compiler average treatment effect; criticism; strengthening a
continuous IV

The extensive use of instrumental variables in empirical research is not
without criticism. For instance, Deaton, (2009) argued that IV estimation
methods are analogous to letting the light fall where it may when looking for
an object and then proclaiming that whatever the light illuminates is what has
been sought after all along. Heckman and Urzua, (2010) expressed similar
sentiments. The crux of their arguments is concerning the tension between
internal and external validities of IV estimates: it is, after all, not clear
if the proclaimed treatment effect among compliers could be readily
generalized to a target subgroup of interest. The inconvenient fact that the
complier subgroup is not identified from the population and could possibly be
altered by different incentive structures only makes the generalizability of
IV estimates more challenging (Joffe,, 2011).

Responding to these criticisms, Imbens, (2010) argued that the internal
validity of an IV estimate is often superior to that of other estimands for
observational data, and that “a credible estimate of the average effect for a
subpopulation is preferred to an estimate of the average for the overall
population with little credibility” (Imbens,, 2010). As discussed by Imbens,
(2010), two approaches could complement an IV estimate. First, although a
valid IV cannot point identify the average treatment effect (ATE), it can
partially identify it; see Swanson et al., (2018) for a recent review and
references therein. Alternatively, some researchers propose additional
identification assumptions that allow the complier average treatment effect to
be generalized to the entire population. Some examples include the principal
ignorability assumption (Jo and Stuart,, 2009; Angrist and Fernandez-Val,,
2010; Ding and Lu,, 2017) and the homogeneity/no-interaction assumption
(Hernán and Robins,, 2006; Wang and Tchetgen Tchetgen,, 2018), among others.

In their original analysis, Baiocchi et al., (2010) used a study design
technique called non-bipartite matching (Lu et al.,, 2001, 2011; Rigdon et
al.,, 2018; Zhang et al.,, 2023) to construct matched pairs comprising
comparable mothers with different excess travel times. Their key innovation is
to consider two study designs: one that utilized all data and created $99,174$
pairs (Design I) and the other that “strengthened” the IV by only pairing two
similar mothers whose IVs were far apart, yielding $49,587$ pairs (Design II).

Three things changed between Design I and Design II. First, the population
under investigation changed; for instance, the proportion of white mothers
increased from $70\%$ in Design I to $85\%$ in Design II (Baiocchi et al.,,
2010, Table 1). Second, the compliance rate changed. Baiocchi et al.,’s (2010)
key insight is that, by forcing the continuous IV to be further apart within
each matched pair, the compliance rate induced by this strengthened IV
increased, improving statistical inference of the effect ratio estimand, an
analogue of the Wald estimate in a matched-pair design. Third, the latent
complier subgroup changed. The original IV in Design I and its strengthened
version in Design II targeted a sample average treatment effect over different
complier subgroups. Baiocchi et al., (2010, Section 5) argued that the effect
ratio estimate obtained from the strengthened design was more “typical” in the
sense that a typical mother was more likely to respond to a strong incentive
in Design II as opposed to a much weaker incentive in Design I.

###  1.4 Manipulating a continuous IV with two purposes; caveats

In his perspective piece, Deaton, (2009) argued that standard statistical
practice would first define a parameter of interest and then construct an
estimator targeting the parameter. In the causal inference literature, a
causal estimand of interest is often defined first, followed by identification
assumptions and estimation procedures operated under these identification
assumptions. On the contrary, IV-based procedures implicitly let the chosen
IV, which may have little to do with the scientific question at hand,
determine the causal estimand.

In this article, we propose a study design framework that builds upon Baiocchi
et al.,’s (2010) key insight and show that an ideal paradigm—stating a well-
defined causal estimand for an identifiable target subgroup and _irrespective
of_ the IV, followed by identification and estimation using an IV—could be
achieved by manipulating a continuous IV in a design-based framework
(Rosenbaum, 2002a, ; Imbens and Rosenbaum,, 2005). Our proposed study design
starts by building matched pairs that satisfy the following three criteria:
(i) the covariates of matched mothers are deemed well balanced according to
formal diagnostic tests; (ii) the covariate distribution of the matched sample
closely resembles that of a target study population; and (iii) the within-pair
difference in the IV is maximized to the extent that the data can afford, with
the goal of maximizing the IV-induced compliance rate. We then propose a
randomization inference (RI) approach to infer partial identification bounds
of the sample average treatment effect and extend the approach to accommodate
an IV-dose-dependent, biased randomization scheme. Our framework deviates from
previous works that conduct RI-based inference for the structural parameter in
a constant, proportional treatment effect model (Imbens and Rosenbaum,, 2005;
Small and Rosenbaum,, 2008) or the effect ratio estimand (Baiocchi et al.,,
2010; Zhang et al.,, 2022). By strengthening the IV and maximizing the
compliance rate, a careful design could help derive more informative bounds
for the desired effect. In Deaton,’s (2009) analogy, researchers clearly state
what is desired to be sought after and then adjust the brightness of the light
(i.e., strengthen the continuous IV) to illuminate the object.

A valid, continuous IV is at the heart of our proposed design. A binary IV
cannot be strengthened; it appears that researchers have to accept the IV as
it is and its associated complier subgroup. Fortunately, a lot of commonly
used IVs, such as variables related to geography and nature, are often
continuous. Alternatively, one may combine multiple binary IVs to construct a
many-category IV that is amenable to being strengthened. An important caveat
is regarding the validity of the IV: if the putative IV is not valid, then
strengthening it could increase the bias if the unmeasured IV-outcome
confounder is correlated with the continuous IV (Heng et al.,, 2023).

##  2 Algorithm

###  2.1 Design considerations

Three elements are key to our proposed study design. First, mothers with
similar observed covariates but different IV-defined exposures, e.g., excess
travel times, need be paired together. One approach would be to dichotomize
the continuous exposure and adopt a conventional matching algorithm designed
for a binary exposure (Rosenbaum,, 2020); however, this practice leads to a
potential violation of the stable unit treatment value assumption (SUTVA)
(Rubin,, 1980, 1986). As an alternative, non-bipartite matching creates
matched pairs without dichotomizing the continuous exposure (Lu et al.,, 2001,
2011; Baiocchi et al.,, 2010; Zhang et al.,, 2023). Second, when the data can
afford, it is useful to maximize the within-matched-pair difference in the IV
dose so that the compliance rate can be maximized. To achieve this, we will
resort to a design device called “dose caliper” (see, e.g., Zhang et al.,,
2023, Section 5). Third, we would like the final matched sample to closely
resemble the template, i.e., a random sample from a target subgroup, to
enhance the generalizability of the clinical conclusions. This design
consideration is referred to as template matching in the literature and has
been studied extensively for a binary exposure (Bennett et al.,, 2020; Zhang,,
2023). Yet to our knowledge, statistical matching methods that build
representative matched pairs with a many-category or continuous exposure are
currently lacking.

###  2.2 A non-bipartite matching algorithm that constructs representative
samples

To accommodate the three design considerations, we propose a modification of
the optimal non-bipartite matching algorithm in Baiocchi et al., (2010). Our
key idea is to add e “phantom units,” or “sinks,” that represent units in the
template and design a discrepancy matrix in such a way that study participants
who do not resemble the template are matched to the sinks and thus eliminated
from the final matched sample.

The discrepancy matrix is built using the following steps. Suppose there are
$n$ observational units and $e<n$ sinks representing a target subgroup. Define
the following $(n+e)\times(n+e)$ distance matrix:

$\mathcal{M}=\left(\begin{array}[]{cccc|ccc}\infty&\delta_{1,2}&\cdots&\delta_{1,n}&\Delta_{1,n+1}&\cdots&\Delta_{1,n+e}\\\
\delta_{2,1}&\infty&\cdots&\delta_{2,n}&\Delta_{2,n+1}&\cdots&\Delta_{2,n+e}\\\
\vdots&&\ddots&\vdots&\vdots&\ddots&\vdots\\\
\delta_{n,1}&\delta_{n,2}&\cdots&\infty&\Delta_{n,n+1}&\cdots&\Delta_{n,n+e}\\\
\hline\cr\Delta_{n+1,1}&\Delta_{n+1,2}&\cdots&\Delta_{n+1,n}&\infty&\cdots&\infty\\\
\vdots&&\ddots&\vdots&\infty&\ddots&\infty\\\
\Delta_{n+e,1}&\Delta_{n+e,2}&\cdots&\Delta_{n+e,n}&\infty&\cdots&\infty\\\
\end{array}\right),$

where $\mathcal{M}_{i,j}$ denotes the $(i,j)$-th entry of $\mathcal{M}.$ Here,
$\mathcal{M}_{i,j}$ represents a distance between two observational units for
$i,j\in[n],$ between two sinks for $i,j\in\\{n+1,n+2,\dots,n+e\\},$ and
between an observational unit and a sink for the remaining entries.

For $i,j\in[n]$, let $\delta_{i,j}$ be a measure of the distance between the
observed covariates $x_{i}$ and $x_{j}$ of units $i$ and $j,$ respectively.
Intuitively, $\delta_{i,j}$ measures homogeneity between two observational
units and can take different specifications. For example, $\delta_{i,j}$ can
be the absolute difference in the estimated propensity scores of units $i$ and
$j$ or (rank-based) Mahalanobis distance between $x_{i}$ and $x_{j}$
(Rosenbaum and Rubin,, 1985). To strengthen the continuous IV, we can
additionally incorporate the dose caliper into $\delta_{i,j}$ (Zhang et al.,,
2023). That is, we let
$\delta_{i,j}=\widetilde{\delta}_{i,j}+C\times\mathbb{I}\\{|\widetilde{Z}_{i}-\widetilde{Z}_{j}|\leq\tau\\},$
where $\widetilde{\delta}_{i,j}$ is a measure of similarity between $x_{i}$
and $x_{j},$ $\widetilde{Z}_{i}$ and $\widetilde{Z}_{j}$ are the continuous IV
doses of units $i$ and $j$, respectively, $\tau$ is the dose caliper and $C$
is a large penalty applied when $\widetilde{Z}_{i}$ and $\widetilde{Z}_{j}$
differ by less than or equal to the caliper size. By setting $C=\infty$ or a
large number, the dose caliper can be implemented as a hard or soft
constraint, respectively (Zhang et al.,, 2023). In the NICU study, such a dose
caliper would force excess travel times within matched pairs to be far apart
by a pre-specified amount $\tau$, thus strengthening the IV and improving the
IV-induced compliance rate. Note that we always set $\delta_{i,i}=\infty$ for
$i\in[n]$ to prevent observational units from being matched to themselves. By
symmetry, $\delta_{i,j}=\delta_{j,i}.$

For $i\in[n]$ and $j\in\\{n+1,n+2,\dots,n+e\\},$ let $\Delta_{i,j}$ be a
measure of the distance between the observed covariates $x_{i}$ and $x_{j}$ of
observational unit $i$ and sink $j$, which represents a unit from the
template. Intuitively, if unit $i$ is dissimilar to sink $j$, we would want
the algorithm to pair them together in order to eliminate the $i$-th unit, in
which case $\Delta_{i,j}$ would need to be small. Hence, one possible
specification is to let $\Delta_{i,j}$ equal an arbitrary large number minus
the (rank-based) Mahalanobis distance between $x_{i}$ and $x_{j}$. On top of
that, we can further incorporate a “reverse” generalizability caliper defined
by a generalizability score such as the “probability of participation,” i.e.,
the conditional probability of being selected into the template given observed
covariates, thereby controlling the degree of resemblance between unit $i$ and
sink $j$ (Cole and Stuart,, 2010; Stuart et al.,, 2011). To this end, let
$\Delta_{i,j}=\widetilde{C}-\widetilde{\Delta}_{i,j}+D\times\mathbb{I}\\{|\hat{S}_{i}-\hat{S}_{j}|\leq\gamma\\},$
where $\widetilde{C}$ is an arbitrary number added to ensure that
$\Delta_{i,j}$ is non-negative; $\widetilde{\Delta}_{i,j}$ is rank-based
Mahalanobis distance between $x_{i}$ and $x_{j}$; $\hat{S}_{i}$ and
$\hat{S}_{j}$ are the estimated probability of participation for units $i$ and
$j$, respectively; $\gamma$ is the generalizability score caliper; and
$D=\infty$ or is a large penalty applied when $\hat{S}_{i}$ and $\hat{S}_{j}$
differ by less than or equal to the caliper size $\gamma$. If the
generalizability scores of unit $i$ and sink $j$ are similar, this implies
that unit $i$ closely resembles sink $j$. By imposing a large penalty $D$, we
prevent them from being matched and thus retain unit $i$ in our final matched
sample. By symmetry, $\Delta_{i,j}=\Delta_{j,i}$.

Lastly, $\mathcal{M}_{i,j}=\infty$ for $i,j\in\\{n+1,n+2,\dots,n+e\\}.$ That
is, we impose an arbitrarily large penalty for any entry corresponding to two
sinks, or to a sink and itself, to prevent them from being matched into pairs.
In practice, $\mathcal{M}_{i,j}$ can be set to a large finite number.

An optimal non-bipartite matching algorithm takes as an input the distance
matrix $\mathcal{M}$ and $n+e$ units and yields $(n+e)/2$ matched pairs that
minimize the sum of the discrepancies. The final matched sample comprises
$(n-e)/2$ pairs of observational units and the remaining $e$ pairs consisting
of an observational unit and a phantom unit. The algorithm runs in polynomial
time and is available in the R package nbpmatching (Derigs,, 1988; Lu et al.,,
2011).

###  2.3 Informal and formal balance diagnostics

An IV-based matched cohort design seeks to embed the retrospective
observational data into an “approximate” randomized encouragement experiment.
A popular downstream analysis approach is to conduct randomization inference
conditional on all the potential outcomes, where the only source of randomness
in statistical inference comes from the probabilistic assignment mechanism of
the IV doses (Imbens and Rosenbaum,, 2005; Baiocchi et al.,, 2010; Heng et
al.,, 2023). Informal and formal balance diagnostics are often used to justify
the _randomization assumption_ , which states that the two IV doses within
each matched pair are randomly assigned (Silber et al.,, 2001; Gagnon-Bartsch
and Shem-Tov,, 2019; Yu,, 2020). For instance, Gagnon-Bartsch and Shem-Tov,’s
(2019) _Classification Permutation Test_ (CPT) provides a powerful machine-
learning-based test for the randomization assumption. If the randomization
assumption is rejected, then some relaxed version of it, such as the biased
randomization assumption, may be tested and a minimal degree of residual
imbalance due to inexact matching on the observed covariates may be quantified
(Chen et al.,, 2023). If there is strong evidence against the randomization
assumption, then one could opt to conduct the primary analysis under a biased
randomization scheme. Regardless of whether the randomization or biased
randomization assumption is used in the primary analysis, a sensitivity
analysis that further examines the impact of unmeasured IV-outcome confounding
is recommended.

##  3 Data and study design

###  3.1 NICU Data; target trial; details on statistical matching

We considered data on $181,762$ preterm babies in the Commonwealth of
Pennsylvania between 1995 and 2004. As stipulated by the American Academy of
Pediatrics, there are six levels of NICUs of increasing technical expertise
and capability: 1, 2, 3A, 3B, 3C and 3D, as well as regional centers 4
(Baiocchi et al.,, 2010). We followed Baiocchi et al., (2010) and defined a
hospital to be high level if it delivered at least 50 preterm babies per year
on average and its NICU was of level 3A-3D or 4; we defined a hospital as low
level if it delivered fewer than 50 preterm babies annually or its NICU was
below level 3A.

Our goal is to embed the retrospective, observational data into a hypothetical
randomized encouragement trial in such a way that the instrument is
strengthened to the extent that the data can afford and the final matched
sample closely resembles the template. We constructed and compared two matched
samples, one constructed using a usual technique (M0) and the other based on
the algorithm proposed in Section 2 (M1). For both designs, we matched mothers
on birth weight, parity, insurance type (fee for service or others), below
poverty, mother’s age, mother’s education, and a single birth indicator, and
matched exactly on the year the data was collected. To prepare for the
subgroup analyses based on M1, we further matched on a categorical, high-risk
variable defined by a combination of a mom’s race, age, and gestational age,
motivated by clinical literature discussed in Section 1.1: the high risk
category equals $0$ if a mother is black, $1$ if a mother is at least $35$
years old and her gestational age does not exceed $36$, $2$ if a mother is
younger than $35$ and her gestational age does not exceed $36$, and $3$ if a
mother’s gestational age exceeds $36$.

The first matched sample M0 retained and matched all $181,722$ mothers into
pairs using a conventional non-bipartite matching algorithm (Lu et al.,, 2001,
2011). For the second matched sample M1, we first created a template subgroup
of mothers by randomly sampling without replacement $10$ percent of mothers
from the entire dataset spanning $1995$ to $2004$. We then implemented the
statistical matching algorithm introduced in Section 2.2, which includes
several design devices such as a generalizability caliper and an IV dose
caliper. We implemented the “reverse” generalizability caliper by first
estimating the “probability of participation” as the generalizability score
(Cole and Stuart,, 2010). Specifically, we pooled the corresponding template
and the entire observational data, created a binary indicator variable $Y_{i}$
($1$ if participant $i$ belongs to the template and $0$ otherwise), and
estimated the “probability of participation” by fitting a logistic regression
model. Following the notation in Section 2.2, we let
$D=500,\,\tilde{C}=1+\max_{i,j}{\tilde{\Delta}}_{i,j},$ and $\gamma$ equal to
the standard deviation of $|S_{i}-S_{j}|$ over all observational units $i$ and
template units $j.$ Then, to implement the dose caliper, we let $C=250$ and
used rank-based Mahalanobis distance between two observational units $i,j$ for
$\tilde{\delta}_{i,j}$. Finally, we considered a dose caliper $\tau$ equal to
the standard deviation of the excess travel times from all mothers. Design M1
yielded $81,766$ matched pairs of mothers whose within-matched-pair rank-based
Mahalanobis distance was minimized subject to the constraints enforced by the
above generalizability and IV dose calipers.

###  3.2 Matched samples

Table 1 contrasts the average excess travel times and covariate balance
between mothers encouraged to deliver at high-level and low-level NICUs in the
matched samples M0 and M1. In both matches, the absolute SMD of almost every
covariate between matched pairs of mothers does not exceed $0.1$, or one-tenth
of one standard deviation (Silber et al.,, 2001). We will formally test the
randomization and relaxed randomization assumption before conducting inference
in Section 7.

Table 1: Degree of encouragement and covariate balance between mothers
encouraged to deliver at high-level (near) versus low-level (far) NICUs in two
matched comparisons using the $1995-2004$ data. Matches M0 and M1 are created
without and with a dose caliper, respectively. For each variable, we report
the mean and the absolute standardized difference. GA: gestational age.

| Match 0 (M0) | Match 1 (M1)
---|---|---
|  $90,861$ matched pairs |  $81,766$ matched pairs
| Near | Far | Absolute | Near | Far | Absolute
| Mean | Mean | SMD | Mean | Mean | SMD
Instrumental variable |  |  |  |  |  |
Excess travel time (min) | 6.07 | 20.97 | 0.93 | 2.07 | 27.10 | 1.88
Covariates |  |  |  |  |  |
Birthweight (g) | 2583 | 2584 | 0.00 | 2597 | 2597 | 0.00
Gestational age (weeks) | 35.12 | 35.13 | 0.00 | 35.21 | 35.19 | 0.01
Gestational diabetes (0/1) | 0.05 | 0.05 | 0.01 | 0.05 | 0.05 | 0.01
Single birth (0/1) | 0.83 | 0.83 | 0.00 | 0.85 | 0.84 | 0.02
Parity | 2.11 | 2.12 | 0.01 | 2.06 | 2.11 | 0.03
Mother’s age (years) | 28.05 | 28.04 | 0.00 | 28.03 | 27.80 | 0.04
Mother’s education (scale) | 3.69 | 3.68 | 0.00 | 3.73 | 3.63 | 0.08
Black (0/1) | 0.16 | 0.16 | 0.00 | 0.16 | 0.16 | 0.00
Below poverty (proportion) | 0.13 | 0.13 | 0.01 | 0.12 | 0.12 | 0.06
Fee for service (0/1) | 0.21 | 0.21 | 0.00 | 0.20 | 0.21 | 0.05
High risk category |  |  |  |  |  |
Black | 0.16 | 0.16 | 0.00 | 0.16 | 0.16 | 0
Non-Black, age $\geq 35$, GA $\leq 36$ | 0.08 | 0.08 | 0.00 | 0.08 | 0.08 | 0
Non-Black, age $<35$, GA $\leq 36$ | 0.40 | 0.39 | 0.00 | 0.40 | 0.40 | 0
Non-Black, GA $>36$ | 0.36 | 0.36 | 0.00 | 0.36 | 0.36 | 0

A key difference between M0 and M1 is the level of separation in the IV
enforced by the study design. The average difference (larger minus smaller
within each matched pair) in excess travel times between mothers encouraged to
deliver at a low-level versus a high-level NICU for M0 and M1 are
$20.97-6.07=14.90$ and $27.10-2.07=25.03$ minutes, respectively. This is in
line with our study design: a large IV dose caliper enforces the excess travel
times for mothers in each matched pair to be further apart, thereby
strengthening the IV and improving the IV-induced compliance rate from
$65.4\%-45.6\%=19.8\%$ in M0 to $71.9\%-35.5\%=36.4\%$ in M1. According to the
analytical results in Heng et al., (2023), in order for M0 to attain the same
power as M1 when testing Fisher’s sharp null hypothesis, M0 needs to have a
sample size that is $(36.4\%/19.8\%)^{2}=3.38$ times that of M1. Design M0’s
sample size is only $1.1$ times that of M1; therefore, from the perspective of
testing Fisher’s sharp null hypothesis, the strengthened design M1 is more
advantageous.

Figure 1: Boxplots of the rank-based Mahalanobis distances of two types of
matched pairs in M1: (i) two observational units (Obs-Obs) or (ii) an
observational unit and a template unit (Obs-Template).

Furthermore, by using the template matching technique and enforcing the
generalizability caliper, the matched sample M1 also exhibits good
generalizability to the template. For example, M1 largely retains the original
composition of black mothers, who comprise $16\%$ of all the mothers in the
observational data, and of mothers with fee-for-service health insurance. This
contrasts the original strengthened IV design in Baiocchi et al., (2010)
which, in transitioning from the unstrengthened to strengthened IV matching
design, shifted the study cohort from $15\%$ to $5\%$ of black mothers and
from $21\%$ to $25\%$ of mothers with fee-for-service health insurance. Figure
1, which illustrates the distribution of rank-based Mahalanobis distances
between matched pairs retained in or eliminated from M1, further demonstrates
the effectiveness of our matching algorithm in maintaining good
generalizability of the final matched sample to the entire study cohort. For
the matched pairs that each consist of two mothers from the observational
database and thus comprise M1, the median rank-based Mahalanobis distance is
$12$ (interquartile range: [5, 31]). For the matched pairs that are eliminated
from M1, the median rank-based Mahalanobis distance is $279$ (interquartile
range: [226, 329]). This is in line with our expectation, as mothers who
resembled the template the least in terms of the covariates had large rank-
based Mahalanobis distances when matched to the template and were thus
eliminated from the final matched sample M1.

##  4 Notation: potential outcomes and estimands

###  4.1 Potential outcomes

Our proposed design consists of $I$ matched pairs of two study participants.
The $j$-th participant in the $i$-th pair is indexed by $ij$. Each participant
$ij$ is associated with a continuous IV
$\widetilde{Z}_{ij}=\widetilde{Z}_{ij}^{\text{obs}}$ e.g., excess travel time
in the NICU study. The design embeds the observational data into a
hypothetical randomized encouragement experiment as follows. Within each
matched pair $i$, fix the two doses of the continuous IV at
$\widetilde{Z}_{i1}^{\text{obs}}$ and $\widetilde{Z}_{i2}^{\text{obs}}$ and
flip a fair coin. If the coin lands heads, assign the IV dose
$\widetilde{Z}_{i1}^{\text{obs}}$ to $j=1$ and
$\widetilde{Z}_{i2}^{\text{obs}}$ to $j=2$; if the coin lands tails, assign
the IV dose $\widetilde{Z}_{i1}^{\text{obs}}$ to $j=2$ and
$\widetilde{Z}_{i2}^{\text{obs}}$ to $j=1$. Let
$\widetilde{\mathbf{Z}}=(\widetilde{Z}_{11},\dots,\widetilde{Z}_{I2})$ denote
the vector of $I\times 2=2I$ IV dose assignments. In a randomized
encouragement experiment, it is the encouragement (i.e., the continuous IV
$\widetilde{Z}$), not the treatment itself, that is randomized either by
nature or an experimenter.

Let $D_{ij}(\widetilde{\mathbf{Z}})$ denote the indicator of whether
participant $ij$ receives the binary treatment under the IV dose vector
$\widetilde{\mathbf{Z}}$. We make the SUTVA assumption (Rubin,, 1980, 1986) so
that $D_{ij}(\widetilde{\mathbf{Z}})$ depends on $\widetilde{\mathbf{Z}}$ only
via $\widetilde{Z}_{ij}$ such that $D_{ij}(\widetilde{\mathbf{Z}})$ can be
simplified to $D_{ij}(\widetilde{Z}_{ij})$. We adopt the notation in Zhang et
al., (2022) and Heng et al., (2023) and define the following potential
outcomes:

$\begin{split}&d_{Cij}\overset{\Delta}{=}D_{ij}(\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}})\leavevmode\nobreak\
\quad\text{and}\quad\leavevmode\nobreak\
d_{Tij}\overset{\Delta}{=}D_{ij}(\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}),\end{split}$
(1)

where $\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}}$
denotes the minimum of two IV doses and
$\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}$ the
maximum of two IV doses in the $i$-th matched pair. Importantly, participant
$ij$’s compliance status is determined by two IV doses in matched pair $i$ and
the compliance status is fixed once the study design is fixed. Participant
$ij$ is said to be a complier with respect to the IV dose pair
$(\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}},\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}})$
if $(d_{Tij},d_{Cij})=(1,0)$, an always-taker if $(d_{Tij},d_{Cij})=(1,1)$, a
never-taker if $(d_{Tij},d_{Cij})=(0,0)$, and a defier if
$(d_{Tij},d_{Cij})=(0,1)$, where $d_{Tij}$ and $d_{Cij}$ are defined as in
(1). We make the so-called monotonicity assumption (Angrist et al.,, 1996) and
assume there are no defiers. This is often a mild assumption as it merely
states that if unit $ij$ accepts the treatment under dose $\widetilde{Z}$,
then it will also accept the treatment under any dose $\widetilde{Z}^{\prime}$
such that $\widetilde{Z}^{\prime}>\widetilde{Z}$. Note that even under the
monotonicity assumption, researchers can still only observe one of the two
potential treatments received, so statistical inference is needed. Under the
monotonicity assumption, the compliance rate in the matched design with
respect to a continuous IV is then defined as

$\iota_{C}=\frac{1}{2I}\sum_{i=1}^{I}\sum_{j=1}^{2}(d_{Tij}-d_{Cij}),$ (2)

and a strengthened-IV design improves the compliance rate by maximizing the IV
dose difference within each matched pair (Baiocchi et al.,, 2010).

Let $R_{ij}(\widetilde{\mathbf{Z}},\mathbf{D})$ denote the outcome of $ij$
under $\widetilde{\mathbf{Z}}$ and $\mathbf{D}=(D_{11},\dots,D_{I2})$. We
assume SUTVA again so that $R_{ij}(\widetilde{\mathbf{Z}},\mathbf{D})$ depends
on $\widetilde{\mathbf{Z}}$ and $\mathbf{D}$ only via $\widetilde{Z}_{ij}$ and
$D_{ij}$. Lastly, we make the exclusion restriction assumption so that the
effect of $\widetilde{Z}_{ij}$ on $R_{ij}$ is only via $D_{ij}$; in other
words, $R_{ij}(\widetilde{Z}_{ij},D_{ij})$ can be written as
$R_{ij}(D_{ij}(\widetilde{Z}_{ij}))$ or $R_{ij}(\widetilde{Z}_{ij})$ as a
shorthand. Define the following potential outcomes:

$\begin{split}&r_{Cij}\overset{\Delta}{=}R_{ij}(\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}})\leavevmode\nobreak\
\quad\text{and}\quad\leavevmode\nobreak\
r_{Tij}\overset{\Delta}{=}R_{ij}(\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}),\end{split}$
(3)

in parallel to $(d_{Tij},d_{Cij})$ in (1). Analogous to $(d_{Tij},d_{Cij})$,
$(r_{Tij},r_{Cij})$ are defined with respect to the continuous IV dose pair
$(\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}},\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}})$.
In the NICU study, $r_{Cij}$ describes the potential infant mortality status
had the mother lived
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}}$
(as opposed to
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}$)
minutes farther away from a high-level NICU compared to a low-level NICU and
hence encouraged to attend a high-level NICU, while $r_{Tij}$ describes the
potential infant mortality status had the mother lived
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}$
(as opposed to
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}}$)
minutes farther away from a high-level NICU compared to a low-level NICU and
hence encouraged to attend a low-level NICU. Importantly, in either case, no
mother is forced to attend any hospital. Potential outcomes
$(r_{Tij},r_{Cij})$ are fixed once the design and hence the two IV doses
$(\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}},\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}})$
within each matched pair are fixed.

###  4.2 Effect ratio estimand; sample average treatment effect

Equipped with these potential outcomes, we now turn to causal estimands.
Baiocchi et al., (2010) consider the following effect ratio estimand:

$\lambda=\frac{\sum_{i=1}^{I}\sum_{j=1}^{2}(r_{Tij}-r_{Cij})}{\sum_{i=1}^{I}\sum_{j=1}^{2}(d_{Tij}-d_{Cij})},$

which is the ratio of the sample average treatment effect of the IV (possibly
after strengthening) on the outcome and that of the IV on the treatment. The
effect ratio estimand $\lambda$ is always well defined provided that the
denominator $\sum_{i=1}^{j}\sum_{j=1}^{2}(d_{Tij}-d_{Cij})\neq 0$, and it
coincides with the Wald estimator under an additional monotonicity assumption.
In the NICU study, the effect ratio can be interpreted as follows: for every
hundred mothers discouraged by the IV from delivering at a high-level NICU,
there were $100\times\lambda$ additional infant deaths (Baiocchi et al.,,
2010). As researchers switch to a different IV-based matched design, such as
one that further separates the two IV doses within each matched pair, the
effect ratio estimand will in general be different. In other words, _the
effect ratio estimand depends on who is matched to whom, and the design
dictates the estimand._

Lastly, we define the following potential outcomes:

$\begin{split}&r_{d=0,ij}\overset{\Delta}{=}R_{ij}(D_{ij}=0)\leavevmode\nobreak\
\quad\text{and}\quad\leavevmode\nobreak\
r_{d=1,ij}\overset{\Delta}{=}R_{ij}(D_{ij}=1),\end{split}$

which correspond to potential binary infant mortality outcomes had mother $ij$
attended a low-level NICU ($D_{ij}=1$) or a high-level NICU ($D_{ij}=0$).
Compared to $r_{Tij}-r_{Cij}$ that describes the unit-level intention-to-treat
effect (with respect to the designed IV), $r_{d=1,ij}-r_{d=0,ij}$ describes
the unit-level treatment effect and arguably is the ultimate estimand of
interest in the NICU study. Below, we will focus on the sample average
treatment effect of $D$ on $R$:

$\kappa=\frac{1}{2I}\sum_{i=1}^{I}\sum_{j=1}^{2}(r_{d=1,ij}-r_{d=0,ij}).$

Compared to the effect ratio estimand $\lambda$, the estimand $\kappa$ depends
on the study cohort but not on the strength of the IV or who is matched to
whom in the design phase.

###  4.3 Partial identification bounds

The estimand $\kappa$ cannot in general be point identified, as data provides
no treatment effect information for non-compliant participants. Partial
identification bounds provide an assumption-lean alternative. Define the
following index sets:
$I_{\text{COM}}=\\{(i,j)\mathrel{\mathop{\ordinarycolon}}(d_{Tij},d_{Cij})=(1,0)\\}$,
$I_{\text{AT}}=\\{(i,j)\mathrel{\mathop{\ordinarycolon}}(d_{Tij},d_{Cij})=(1,1)\\}$,
and
$I_{\text{NT}}=\\{(i,j)\mathrel{\mathop{\ordinarycolon}}(d_{Tij},d_{Cij})=(0,0)\\}$.
Note that

$\begin{split}\kappa=\frac{1}{N}\sum_{ij}(r_{d=1,ij}-r_{d=0,ij})=&\frac{1}{N}\left\\{\sum_{(i,j)\in
I_{\text{COM}}}(r_{d=1,ij}-r_{d=0,ij})+\sum_{(i,j)\in I_{\text{AT}}\cup
I_{\text{NT}}}(r_{d=1,ij}-r_{d=0,ij})\right\\},\end{split}$

where the first term can be further simplified to

$\sum_{(i,j)\in I_{\text{COM}}}(r_{d=1,ij}-r_{d=0,ij})=\sum_{(i,j)\in
I_{\text{COM}}}(r_{Tij}-r_{Cij})=\sum_{ij}(r_{Tij}-r_{Cij}).$

Let $r_{d=1,ij}-r_{d=0,ij}\in[K_{0},K_{1}]$ for all $ij$ and for some
$K_{0}\leq K_{1}<\infty.$ Without loss of generality, assume $K_{0}\geq 0$.
The second term is then bounded as follows:

$K_{0}\times(N-|I_{\text{COM}}|)\leq\sum_{(i,j)\in I_{\text{AT}}\cup
I_{\text{NT}}}(r_{d=1,ij}-r_{d=0,ij})\leq K_{1}\times(N-|I_{\text{COM}}|).$

Because $\sum_{ij}(d_{Tij}-d_{Cij})=\sum_{(i,j)\in
I_{\text{COM}}}(d_{Tij}-d_{Cij})=|I_{\text{COM}}|$, we then have

$K_{0}\times\left(N-\sum_{ij}(d_{Tij}-d_{Cij})\right)\leq\sum_{(i,j)\in
I_{\text{AT}}\cup I_{\text{NT}}}(r_{d=1,ij}-r_{d=0,ij})\leq
K_{1}\times\left(N-\sum_{ij}(d_{Tij}-d_{Cij})\right).$ (4)

Put together, without additional assumptions, the target estimand $\kappa$ is
lower bounded by

$\text{LB}\mathrel{\mathop{\ordinarycolon}}=\frac{1}{N}\sum_{ij}(r_{Tij}-r_{Cij})-K_{0}\times\frac{1}{N}\sum_{ij}(d_{Tij}-d_{Cij})+K_{0},$

and upper bounded by

$\text{UB}\mathrel{\mathop{\ordinarycolon}}=\frac{1}{N}\sum_{ij}(r_{Tij}-r_{Cij})-K_{1}\times\frac{1}{N}\sum_{ij}(d_{Tij}-d_{Cij})+K_{1}.$

Three remarks follow. First, the lower and upper bounds can be achieved when
the treatment effect among all the non-compliers attain the minimum and
maximum ($K_{0}$ and $K_{1}$, respectively). Second, the length of the partial
identification interval, $\text{UB}-\text{LB}$, equals
$(K_{1}-K_{0})\times(1-\iota_{C})$, where $\iota_{C}$ is the compliance rate
defined in (2). Thus, by maximizing the compliance rate, a strengthened-IV
design could help derive a narrower partial identification interval for
$\kappa$. In an ideal situation where $\iota_{C}=1$, we have
$\text{LB}=\text{UB}$ and the identification of $\kappa$ is achieved by
design. Third, by dividing $N-\sum_{ij}(d_{Tij}-d_{Cij})$ in equation (4),
$K_{0}$ and $K_{1}$ can also be interpreted as the lower and upper bounds of
the SATE among non-compliers.

The target parameters LB and UB involve both potential treatments received
$(d_{Tij},d_{Cij})$ and both potential clinical outcomes $(r_{Tij},r_{Cij})$;
for each study unit $ij$, either $(d_{Tij},r_{Tij})$ or $(d_{Cij},r_{Cij})$ is
observed, so inference is needed for the LB and UB.

##  5 Inference for partial identification bounds

###  5.1 Inference under a randomization assumption

One widely-adopted downstream analysis strategy for matched cohort data is
randomization inference treating the matched cohort data as from a finely
stratified experiment (Rosenbaum, 2002b, ; Rosenbaum,, 2010; Fogarty, 2018a, ;
Fogarty, 2018b, ). We will focus on developing a valid statistical test and
confidence statement for the partial identification bounds under the
randomization assumption in this section; methods that further take into
account within-pair residual bias that persists after matching (Gagnon-Bartsch
and Shem-Tov,, 2019; Chen et al.,, 2023) will be discussed next.

Write
$\mathcal{F}=\\{(\mathbf{x}_{ij},d_{Tij},d_{Cij},r_{Tij},r_{Cij})\mathrel{\mathop{\ordinarycolon}}i=1,\dots,I,j=1,2\\}$,
where $d_{Tij}$, $d_{Cij}$, $r_{Tij}$ and $r_{Cij}$ are defined as in (1) and
(3). Write $\mathcal{Z}$ for the set containing the $2^{I}$ possible values
$\widetilde{\mathbf{z}}$ of IV dose assignments $\widetilde{\mathbf{Z}}$, so
that $\widetilde{\mathbf{z}}\in\mathcal{Z}$ if each $\widetilde{{z}}_{ij}$
equals $\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}$ or
$\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}$. Let
$\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee}=(\widetilde{Z}^{\text{obs}}_{11}\vee\widetilde{Z}^{\text{obs}}_{12},\dots,\widetilde{Z}^{\text{obs}}_{I1}\vee\widetilde{Z}^{\text{obs}}_{I2})$
and
$\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge}=(\widetilde{Z}^{\text{obs}}_{11}\wedge\widetilde{Z}^{\text{obs}}_{12},\dots,\widetilde{Z}^{\text{obs}}_{I1}\wedge\widetilde{Z}^{\text{obs}}_{I2})$.
In randomization inference with a continuous IV $\widetilde{\mathbf{Z}}$, the
only probability distribution that enters statistical inference is the
conditional probability
$\pi_{i}=P(\widetilde{Z}_{i1}=\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2},\widetilde{Z}_{i2}=\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\mid\mathcal{F},\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee},\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge})$
that characterizes the underlying IV dose assignment mechanism in each matched
pair $i$. Assuming the IV is valid and the randomization assumption holds, we
have $\pi_{i}=1/2$ for $i=1,\cdots,I$. Define a shorthand $Z_{ij}=1$ if
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\vee\widetilde{Z}_{i2}^{\text{obs}}$,
$Z_{ij}=0$ if
$\widetilde{Z}_{ij}=\widetilde{Z}_{i1}^{\text{obs}}\wedge\widetilde{Z}_{i2}^{\text{obs}}$
and $\mathbf{Z}=(Z_{11},Z_{12},\dots,Z_{I2})$. Under the matched pair design
conditional on $\mathcal{Z}$, $Z_{i1}+Z_{i2}=1$ for each $i$. Let
$R_{ij}=Z_{ij}r_{Tij}+(1-Z_{ij})r_{Cij}$ denote the observed outcome and
$D_{ij}=Z_{ij}d_{Tij}+(1-Z_{ij})d_{Cij}$ denote the observed treatment
received. We consider the null hypothesis
$H_{0}^{L}\mathrel{\mathop{\ordinarycolon}}\text{LB}=l$ and an ever-larger
experiment with $I\rightarrow\infty$. Proposition 1 derives an asymptotically
valid test for $H_{0}^{L}$.

######  Proposition 1.

Define the test statistic

$T\left(l;K_{0}\right)=\frac{1}{I}\sum_{i=1}^{I}\left\\{\sum_{j=1}^{2}Z_{ij}\left(R_{ij}-K_{0}D_{ij}\right)-\sum_{j=1}^{2}\left(1-Z_{ij}\right)\left(R_{ij}-K_{0}D_{ij}\right)\right\\}-(l-K_{0})$

and
$S^{2}\left(K_{0}\right)=\frac{1}{I(I-1)}\sum_{i=1}^{I}\left(\hat{\tau}_{i}-\bar{\tau}\right)^{2},$
where $\hat{\tau}_{i}=\sum_{j=1}^{2}(2Z_{ij}-1)(R_{ij}-K_{0}D_{ij})$ and
$\bar{\tau}=\frac{1}{I}\sum_{i=1}^{I}\hat{\tau}_{i}$. Under mild regularity
conditions and conditional on
$\mathcal{F},\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee},\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge}$,
the test that rejects $H_{0}^{L}\mathrel{\mathop{\ordinarycolon}}\text{LB}=l$
when $|T(l;K_{0})|\geq z_{1-\alpha/2}\sqrt{S^{2}(K_{0})}$ is an asymptotically
valid level-$\alpha$ test, where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-th
quantile of the standard normal distribution.

######  Remark 1.

By replacing $K_{0}$ with $K_{1}$ in Proposition 1, a valid test can be
derived analogously for testing the null hypothesis
$H_{0}^{U}\mathrel{\mathop{\ordinarycolon}}\text{UB}=u$. A one-sided
level-$\alpha$ confidence interval of LB is
$\left[T(l;K_{0})+l-z_{1-\alpha}\sqrt{S^{2}(K_{0})},\,\infty\right)$, where
$z_{1-\alpha}$ is the $(1-\alpha)$-th quantile of the standard normal
distribution. A valid confidence interval of $[\text{LB},\text{UB}]$ may then
be obtained by combining the two one-sided, level-$\alpha/2$ confidence
intervals of LB and UB.

######  Remark 2.

The variance estimator $S^{2}(K_{0})$ in general overestimates the true
variance (Imai,, 2008). One set of sufficient conditions for it to be an
unbiased estimator of the true variance is when the proportional treatment
effect model $r_{Tij}-r_{Cij}=\beta(d_{Tij}-d_{Cij})$ holds (Small and
Rosenbaum,, 2008), the monotonicity assumption holds, and one of the following
three conditions holds: (i) all units are compliers; (ii) each pair has one
complier and one always-taker or never-taker; (iii) no units are compliers.
Such a condition appears restrictive in practice, so the variance estimator
$S^{2}(l)$ would tend to overestimate the true variance and the test based on
it would tend to be conservative.

A standard approach to improving the asymptotic precision of randomization-
based inference is through regression adjustment (see, e.g., Lin,, 2013;
Fogarty, 2018a, ; Fogarty, 2018b, ; Li and Ding,, 2020; Su and Ding,, 2021;
Lei and Ding,, 2021; Zhao and Ding,, 2022, among others). Proposition 2 adapts
this idea to making inference of the proposed partial identification bounds by
proposing a regression-assisted variance estimator for $T(l;K_{0})$.

######  Proposition 2.

Let $Q$ denote an $I\times p$ matrix with $p<I$ and
$H_{Q}=Q(Q^{T}Q)^{-1}Q^{T}$ be the hat matrix for $Q$ with $h_{Qi}$ on the
$i$-th diagonal. Let $\hat{\tau}_{Q}$ be an $I\times 1$ column vector with the
$i$-th entry $\hat{\tau}_{i}/\sqrt{1-h_{Qi}}$, where
$\hat{\tau}_{i}=\sum_{j=1}^{2}(2Z_{ij}-1)(R_{ij}-K_{0}D_{ij})$ as defined in
Proposition 1. Let $\boldsymbol{\tau}$ be an $I\times 1$ column vector with
the $i$-th entry
$\tau_{i}=\frac{1}{2}\sum_{j=1}^{2}\left\\{r_{Tij}-r_{Cij}-K_{0}(d_{Tij}-d_{Cij})\right\\}$.
Define $S^{2}_{Q}(K_{0})=I^{-2}\hat{\tau}^{T}_{Q}(I-H_{Q})\hat{\tau}_{Q}$.
Under conditions S1-S2 in Supplemental Material A, as $I\rightarrow\infty$,

$IS^{2}_{Q}(K_{0})-var\\{\sqrt{I}T(l;K_{0})\mid\mathcal{F},\widetilde{\mathbf{Z}}_{\vee},\widetilde{\mathbf{Z}}_{\wedge}\\}\xrightarrow{p}\lim\limits_{I\to\infty}\frac{1}{I}\boldsymbol{\tau}^{T}(I-H_{Q})\boldsymbol{\tau}\geq
0.$

######  Remark 3.

Proposition 2 indicates that $IS^{2}_{Q}(K_{0})$ is a conservative estimator
for the variance of $\sqrt{I}T(l;K_{0})$ in probability, so testing
$H_{0}^{L}$ based on the statistic $T(l;K_{0})/{\sqrt{S^{2}_{Q}(K_{0})}}$ and
a Gaussian reference distribution are asymptotically valid. When
$Q=\boldsymbol{e}$, a column vector of 1’s, $S^{2}_{Q}(K_{0})$ reduces to the
classical variance estimator $S^{2}(K_{0})$ as defined in Proposition 1. In
general, for $Q=\left[\boldsymbol{e},\boldsymbol{V}\right]$, where
$\boldsymbol{V}$ is a matrix whose $i$-th row contains centered covariate
$(\boldsymbol{x}_{i1},\boldsymbol{x}_{i2})$,
$IS^{2}_{\boldsymbol{e}}(K_{0})-IS^{2}_{Q}(K_{0})\xrightarrow{p}\beta_{V}^{T}\Sigma_{V}^{-1}\beta_{V}\geq
0$, where
$\Sigma_{V}=\lim\limits_{I\to\infty}I^{-1}\boldsymbol{V}^{T}\boldsymbol{V}$
and
$\beta_{V}=\lim\limits_{I\to\infty}I^{-1}\boldsymbol{V}^{T}\boldsymbol{\tau}$.
Therefore, $IS^{2}_{Q}(K_{0})$ is asymptotically no more conservative than
$IS^{2}(K_{0})$ if $Q$ contains any covariate information, though the degree
of conservativeness reduction depends on how well the covariates predict
$\hat{\tau}_{i}$ (Fogarty, 2018b, ).

###  5.2 Inference under a IV-dose-dependent, biased randomization assumption

Although matching helps remove most overt bias, residual bias due to inexact
matching and from unmeasured confounders could bias the IV dose assignment
within each pair and induce a so-called biased randomization scheme
(Rosenbaum, 2002b, ; Rosenbaum,, 2010; Fogarty,, 2020; Fogarty et al.,, 2021;
Chen et al.,, 2023). Randomization inference under a biased randomization
scheme is useful in two scenarios. First, when the minimal degree of bias due
to residual confounding from observed covariates can be quantified by formal
statistical tests (Chen et al.,, 2023), primary analyses should adopt a biased
randomization scheme. Second, since the observed data contains no information
about unmeasured confounding bias, a sensitivity analysis that further relaxes
the randomization or biased randomization assumption in the primary analysis
is always warranted.

Recall that
$\pi_{i}=P(\widetilde{Z}_{i1}=\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2},\widetilde{Z}_{i2}=\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\mid\mathcal{F},\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee},\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge})$.
Following Rosenbaum and Rubin, (2023), we consider a principal unobserved
covariate, $u_{ij}\in[0,1]$, associated with each participant $ij$ and
supported on the unit interval. We consider the following biased IV dose
assignment model $\mathcal{M}_{\gamma}$ within each matched pair $i$ with an
observed IV dose pair
$(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2},\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})$
and some $\gamma\geq 0$:

$\begin{split}\pi_{i}=\frac{\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\cdot
u_{i1}\right\\}}{\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\cdot
u_{i1}\right\\}+\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\cdot
u_{i2}\right\\}},\end{split}$

for $i=1,\dots,n$, or equivalently,

$\frac{1}{\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\right\\}+1}\leq\pi_{i}\leq\frac{\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\right\\}}{\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\right\\}+1}.$

The biased randomization model $\mathcal{M}_{\gamma}$ is analogous to that in
Rosenbaum, (1989) and reduces to the classical Rosenbaum bounds model
(Rosenbaum,, 1987; Rosenbaum, 2002b, ) when
$(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2},\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})=(1,0)$.
Under $\mathcal{M}_{\gamma}$, the degree of bias in each matched pair depends
critically on the IV dose difference. For instance, if a matched pair happens
to pair two participants with the same IV dose, that is,
$\widetilde{Z}^{\text{obs}}_{i1}=\widetilde{Z}^{\text{obs}}_{i2}$, then
$\mathcal{M}_{\gamma}$ would imply $\pi_{i}=1/2$, that is, IV dose assignment
becomes randomized for this pair. As the IV dose difference increases,
$\mathcal{M}_{\gamma}$ would stipulate a larger degree of potential bias
within each pair. Instead of allowing a dose-dependent maximal bias within
each pair, Rosenbaum’s classical $\Gamma$ sensitivity analysis model
(Rosenbaum,, 1987; Fogarty et al.,, 2021) places a uniform bound on the
maximal biased randomization probability across all pairs, which could be
unduly conservative for some matched pairs, e.g., those with
$\widetilde{Z}^{\text{obs}}_{i1}=\widetilde{Z}^{\text{obs}}_{i2}$.

Write
$\mathcal{F}^{\prime}=\\{(\mathbf{x}_{ij},u_{ij},d_{Tij},d_{Cij},r_{Tij},r_{Cij})\mathrel{\mathop{\ordinarycolon}}0\leq{u}_{ij}\leq
1,\leavevmode\nobreak\ i=1,\dots,I,\leavevmode\nobreak\ j=1,2\\}$. Theorem 1
derives an asymptotically valid test for the null hypothesis $H_{0}^{L}$ under
the biased randomization model $\mathcal{M}_{\gamma}$ for some $\gamma\geq 0$,
.

######  Theorem 1.

For a fixed $\gamma\geq 0$, define
$\Gamma_{i}=\exp\left\\{\gamma\cdot\left(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\right)\right\\}\geq
1$ to be the IV-dose-dependent odds $\pi_{i}/(1-\pi_{i})$ under
$\mathcal{M}_{\gamma}$ in matched pair $i$. For each matched pair $i$, define
the following scaled treated-minus-control difference in the observed outcome:

$\hat{\tau}_{i,\Gamma_{i}}=\frac{\Gamma_{i}+1}{4\Gamma_{i}}\left\\{\left(\Gamma_{i}+1\right)\hat{\tau}_{i}-\left(\Gamma_{i}-1\right)\left|\hat{\tau}_{i}\right|\right\\},$
(5)

where $\hat{\tau}_{i}$ is defined as in Proposition 1. Let
$\bar{\tau}_{\gamma}=\frac{1}{I}\sum_{i=1}^{I}\hat{\tau}_{i,\Gamma_{i}}$ be
the sample average across all pairs and
$S^{2}\left(\gamma;K_{0}\right)=\frac{1}{I(I-1)}\sum_{i=1}^{I}\left(\hat{\tau}_{i,\Gamma_{i}}-\bar{\tau}_{\gamma}\right)^{2}$
the usual variance estimator for the sample mean. Under the biased IV
assignment model $\mathcal{M}_{\gamma}$ and conditional on
$\mathcal{F}^{\prime},\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee},\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge}$,
the test that rejects $H^{L}_{0}\mathrel{\mathop{\ordinarycolon}}\text{LB}=l$
when

$\bar{\tau}_{\gamma}-(l-K_{0})\geq
z_{1-\alpha}\sqrt{S^{2}\left(\gamma;K_{0}\right)}$ (6)

is an asymptotically valid level-$\alpha$ test, where $z_{1-\alpha}$ is the
$(1-\alpha)$-th quantile of the standard normal distribution. A valid one-
sided, level-$\alpha$ confidence interval can then be constructed by inverting
the hypothesis testing.

######  Remark 4.

When $\gamma=0$, $\mathcal{M}_{\gamma}$ reduces to the the randomization
scheme studied in Section 5.1. In this case, $\Gamma_{i}=1$ for all $i$, and
$\hat{\tau}_{i,\Gamma_{i}}$ defined in (5) reduces to the usual treated-minus-
control difference in the observed outcome in each pair. The test (6) and the
associated confidence interval reduce to those studied in Section 5.1 under
the randomization assumption.

######  Remark 5.

By replacing $K_{0}$ with $K_{1}$ in Theorem 1, a valid test can be derived
analogously for testing the null hypothesis
$H_{0}^{U}\mathrel{\mathop{\ordinarycolon}}\text{UB}=u$. A valid
level-$\alpha$ confidence interval of the partial identification bound
$[\text{LB},\text{UB}]$ may then be obtained by combining two one-sided,
level-$\alpha/2$ confidence intervals of LB and UB.

The test derived in Theorem 1 also applies to testing the effect ratio
estimand under $\mathcal{M}_{\gamma}$, as formalized in Proposition 3.

######  Proposition 3.

Consider testing the null hypothesis
$H^{\lambda}_{0}\mathrel{\mathop{\ordinarycolon}}\lambda=\lambda_{0}$
regarding the effect ratio under the biased IV assignment model
$\mathcal{M}_{\gamma}$. Let
$\hat{\tau}_{i}=\sum_{j=1}^{2}(2Z_{ij}-1)(R_{ij}-\lambda_{0}D_{ij})$. Define
$\hat{\tau}_{i,\Gamma_{i}}$, $\bar{\tau}_{\gamma}$, and
$S^{2}\left(\gamma;\lambda_{0}\right)$ as in Theorem 1. Conditional on
$\mathcal{F}^{\prime}$, $\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee}$, and
$\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge}$, the test that rejects
$H^{\lambda}_{0}$ when $\bar{\tau}_{\gamma}\geq
z_{1-\alpha}\sqrt{S^{2}\left(\gamma;\lambda_{0}\right)}$ is an asymptotically
valid level-$\alpha$ test.

##  6 Simulation studies

We have two goals in the simulation studies. In Section 6.1, we illustrate how
strengthening an IV in the design stage helps improve the compliance rate and
thus narrow the partial identification bounds and deliver more informative
inferential targets. In Section 6.2, we generate synthetic data modeled after
the real data from the case study and verify the level of the proposed tests.

###  6.1 Simulation study I: study design and its effect on the causal
estimand

We generate potential outcomes data for each participant
$n\in\\{1,\dots,N\\}$. We consider a 5-dimensional covariate as follows:
$X_{n,1}\sim\mathcal{N}(0,1),\,X_{n,2}\sim\mathcal{N}(2,5),\,X_{n,3}\sim\text{Unif}[1,3],\,X_{n,4}\sim\text{Unif}[-2,0]$
and $X_{n,5}\sim\text{Bernoulli}(0.5)$. Each participant is associated with a
step function $D_{n}(Z=z)$ with a jump at a randomly sampled threshold
$T_{n}\sim\text{Unif}[20,30]$. The step function $D_{n}(Z=z)$ describes the
potential treatment received for each participant as a function of the IV dose
$Z$. According to this data-generating process, each participant has one’s own
incentive structure. Factor $1$ specifies the distribution of the observed IV
dose $Z^{obs}_{n}$ for participant $n$:

Factor 1:

We consider two scenarios for $\widetilde{Z}^{obs}_{n}$: (1)
$\widetilde{Z}^{obs}_{n}\sim\text{Unif}[5,50]$ is randomized; and (2)
$\widetilde{Z}^{obs}_{n}=4X_{n,2}+6X_{n,3}+\varphi_{n}$ is correlated with
observed covariates, where $\varphi_{n}\sim\text{Unif}[0,2]$.

The observed treatment received $D^{obs}_{n}$ is then determined by
$D^{obs}_{n}=D_{n}(\widetilde{Z}^{obs}_{n})=\mathbb{I}(\widetilde{Z}_{n}^{obs}>T_{n}).$
Using the covariates and the observed IV dose $\widetilde{Z}^{obs}_{n},$ we
perform non-bipartite matching. Factor 2 specifies three matching algorithms
under consideration:

Factor 2:

$\widetilde{\mathcal{M}}_{1},\,\widetilde{\mathcal{M}}_{2}$ and
$\widetilde{\mathcal{M}}_{3}$ refer to matching designs with no caliper, a
small caliper of $7$ minutes, and a large caliper of $15$ minutes on the
absolute difference between two within-matched-pair observed IV doses,
respectively.

Each matching design yields $I=N/2$ matched pairs of participants, where each
matched pair is indexed by $i\in[I]$ and each within-matched-pair participant
indexed by $j\in\\{1,2\\}.$ The two observed IV doses within each matched
pair, $\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2}$ and
$\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2},$ are fixed in each
matched design. We then determine participant $ij$’s compliance status
$S_{ij}$ as follows:

$\displaystyle\text{If
}D_{ij}(\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2})$
$\displaystyle=D_{ij}(\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2})=0,\text{
then }S_{ij}=(0,0);$ $\displaystyle\text{If
}D_{ij}(\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2})$
$\displaystyle=0\text{ and
}D_{ij}(\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2})=1,\text{ then
}S_{ij}=(1,0);$ $\displaystyle\text{If
}D_{ij}(\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2})$
$\displaystyle=D_{ij}(\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2})=1,\text{
then }S_{ij}=(1,1),$

where $S_{ij}\in\\{(0,0),(1,0),(1,1)\\}$ indicates that participant $ij$ is a
never-taker, complier or always-taker with respect to the two observed IV
doses in matched pair $i$, respectively.

Finally, we generate the potential outcomes $r_{d=0,ij}\sim\mathcal{N}(0,1)$.
Factor 3 specifies the unit-level treatment effect and hence the data-
generating process for $r_{d=1,ij}$:

Factor 3:

We consider two scenarios for $r_{d=1,ij}$: (1)
$r_{d=1,ij}=r_{d=0,ij}+\psi_{ij}$, where the unit-level treatment effect
$\psi_{ij}$ is independent of participant $ij$’s compliance status and follows
$\psi_{ij}\sim\text{Unif}[4,6]$; and (2) $r_{d=1,ij}=r_{d=0,ij}+\phi_{ij},$
where the unit-level treatment effect $\phi_{ij}$ depends on participant
$ij$’s compliance status as follows: $\phi_{ij}\sim\text{Unif}[2,5]$ for a
complier, $\phi_{ij}\sim\text{Unif}[4,6]$ for an always-taker, and
$\phi_{ij}\sim\text{Unif}[1,3]$ for a never-taker.

Participant $ij$’s potential outcomes $r_{T_{ij}}$ and $r_{C_{ij}}$ are then
determined by the compliance status $S_{ij}$ and the potential outcomes
$(r_{d=0,ij},r_{d=1,ij})$ as follows: $(r_{T_{ij}},r_{C_{ij}})$ equals
$(r_{d=0,ij},r_{d=0,ij})$ for never-takers, $(r_{d=1,ij},r_{d=0,ij})$ for
compliers, and $(r_{d=1,ij},r_{d=1,ij})$ for always-takers. For each matched
dataset, we calculate $LB$ and $UB$ as described in Section 4.3 with
$(K_{0},K_{1})\in\\{(4,6),(1,6)\\}$ for Scenario 1 or 2 in Factor 3,
respectively. We generate data for $N=1,000$ participants and repeat the
simulation $1,000$ times under each data-generating process.

Table 2: Average widths of the partial identification bounds across $1,000$ simulations. In the first column, $F_{ij}$ refers to the $j$-th scenario in the $i$-th factor. | $\widetilde{\mathcal{M}}_{1}$ | $\widetilde{\mathcal{M}}_{2}$ | $\widetilde{\mathcal{M}}_{3}$
---|---|---|---
$F_{11}\times F_{31}$ | 1.03 | 0.72 | 0.41
$F_{11}\times F_{32}$ | 2.57 | 1.80 | 1.02
$F_{12}\times F_{31}$ | 1.83 | 1.33 | 0.83
$F_{12}\times F_{32}$ | 4.58 | 3.34 | 2.07

Table 2 summarizes the empirical means of the widths of partial identification
bounds, $\Delta=UB-LB$, across $1,000$ simulations for each combination of the
three factors. In the first column of Table 2, $F_{ij}$ denotes the $j$-th
scenario in the $i$-th factor; for instance, $F_{11}\times F_{32}$ denotes the
simulation setting of randomized IVs $\widetilde{Z}_{k}^{obs}$ and potential
outcomes $r_{d=1,k}$ that depend on the compliance status. Across the four
data-generating processes specified by Factors $1$ and $3$, the more we
strengthen the continuous IV, the narrower the widths of the partial
identification bounds, yielding more useful information on the SATE. This is
in line with our expectation, as the width of a partial identification bound
is proportional to the compliance rate, which increases with the IV dose
caliper.

Table S3 in Supplemental Material D further summarizes within-matched-pair
covariate balance under each data-generating process using the average
absolute value of the standardized mean difference (SMD) across $1,000$
simulations. The average absolute value of the SMD increases with the IV dose
caliper, with $X_{2}$ and $X_{3}$ exhibiting the largest increase when they
are correlated with the IV dose as specified by Factor 1 Scenario 2.
Nevertheless, SMDs are almost always less than $0.1$, which usually indicates
a good covariate balance in empirical studies (Silber et al.,, 2001).

###  6.2 Simulation study II: validity of the proposed tests

In this simulation study, we verify the validity of the proposed statistical
tests under randomization inference (Proposition 1) and biased randomization
inference (Theorem 1). For every unit $j$ in matched pair $i$, we first
generate potential outcomes data. Specifically, we sample two continuous IV
doses $\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2}$ and
$\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2}$ with replacement from
the real NICU dataset. Unit $ij$’s compliance status
$S_{ij}(\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2},\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2})\in\\{(1,0),(1,1),(0,0)\\}$
is generated as follows:

$\begin{split}&P\\{S_{ij}=(1,0)\\}=C^{-1};\\\
&P\\{S_{ij}=(1,1)\\}=C^{-1}\cdot\exp\\{-0.2\times(\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2}-\widetilde{Z}^{obs}_{i1}\wedge\widetilde{Z}^{obs}_{i2})\\};\\\
&P\\{S_{ij}=(0,0)\\}=C^{-1}\cdot\exp\\{-\alpha_{i}(\widetilde{Z}^{obs}_{i1}\vee\widetilde{Z}^{obs}_{i2})\\},\end{split}$

where $C$ is a normalizing constant, and $ij$’s compliance status depends on
the two IV doses in each matched pair. The latent compliance category $S_{ij}$
then determines $(d_{Tij},d_{Cij})$ as defined in (1).

We consider two scenarios for each of the two types of potential outcomes
$(r_{d=0,ij},r_{d=1,ij})$:

Continuous case:

Scenario 1: $r_{d=0,ij}\sim\mathcal{N}(0,1)$ and $r_{d=1,ij}=r_{d=0,ij}$;
Scenario 2: $r_{d=0,ij}\sim\mathcal{N}(0,1)$ and
$r_{d=1,ij}=r_{d=0,ij}+\tau_{ij},$ where $\tau_{ij}\sim\text{Unif}[-1,1]$.

Binary case:

Scenario 1: $P(r_{d=0,ij}=1)=\text{expit}(0.5),$ where
$\text{expit}(x)=\exp(x)/(1+\exp(x))$ is the inverse of the logit function,
and $r_{d=1,ij}=r_{d=0,ij}$; Scenario 2: $P(r_{d=0,ij}=1)=\text{expit}(0.5)$
and $P(r_{d=1,ij}=1)=\text{expit}(0.5+\tau_{ij}),$ where
$\tau_{ij}\sim\text{Unif}[-1,1]$.

For both outcome types, Scenario 1 represents a sharp null hypothesis of no
unit-level treatment effect, while Scenario 2 allows a heterogeneous treatment
effect. Once we generate the potential outcomes data, for every matched pair
$i$, we let the participant with the larger potential outcome under control,
i.e., $r_{d=0,ij}$, be more likely to receive the treatment under a biased IV
dose assignment mechanism $\pi_{i},$ which is introduced below. In this way,
the IV doses are correlated with both the treatment assignment and the
potential outcomes $(r_{d=0,ij},r_{d=1,ij})$.

According to this data generating process, $ij$’s potential outcomes
$(r_{T_{ij}},r_{C_{ij}})$ are completely specified by the compliance category
$S_{ij}$ and potential outcomes $(r_{d=1,ij},r_{d=0,ij})$. Finally, we
generate the observed IV dose $\widetilde{Z}_{ij}^{obs}$ according to one of
the following IV dose assignment mechanisms:

$\begin{split}\pi_{i}&=P(\widetilde{Z}_{i1}=\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2},\widetilde{Z}_{i2}=\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2}\mid\mathcal{F}^{\prime},\widetilde{\mathbf{Z}}^{\text{obs}}_{\vee},\widetilde{\mathbf{Z}}^{\text{obs}}_{\wedge})\\\
&\sim\text{Unif}\left[\frac{1}{2},\frac{\exp\left\\{\gamma(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})\right\\}}{\exp\left\\{\gamma(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})\right\\}+1}\right]\leavevmode\nobreak\
(\textbf{Mechanism I}),\leavevmode\nobreak\ \text{or}\\\
&=\frac{\exp\left\\{\gamma(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})\right\\}}{\exp\left\\{\gamma(\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})\right\\}+1}\leavevmode\nobreak\
(\textbf{Mechanism II}).\end{split}$

The observed data comprises the observed IV dose $\widetilde{Z}_{ij}^{obs},$
treatment $D_{ij},$ and outcome $R_{ij}$, with the observed $D_{ij}$ and
$R_{ij}$ completely determined by $\widetilde{Z}_{ij}^{obs}$ and the potential
outcomes data. We vary the following factors in the simulation study:

Factor 1:

Number of matched pairs $I$: $100$, $500$, $1000$ and $2000$.

Factor 2:

IV dose assignment probability $\gamma$: $0$, $0.025$ and $0.05$. If
$\gamma=0$, the two IV doses are randomly assigned within each pair and the
randomization assumption holds in either mechanism. Otherwise, the IV dose
assignment could be biased within each pair and the degree of bias depends on
the difference between the two IV doses.

In Scenario 1, we have $(K_{0},K_{1})=(0,0)$ regardless of whether the outcome
is continuous or binary. In Scenario 2, $(K_{0},K_{1})=(-1,1)$ for both the
continuous and binary potential outcomes. In each setting, we repeat the
simulation $1,000$ times. For each simulated dataset, we calculate the ground
truth
$LB_{0}=\frac{1}{N}\sum_{ij}(r_{Tij}-r_{Cij})-K_{0}\times\frac{1}{N}\sum_{ij}(d_{Tij}-d_{Cij})+K_{0}$
and construct the 95% confidence interval for $LB$ according to Proposition 1
and Theorem 1. We apply Proposition 1 and Theorem 1 to data generated under
the biased IV dose assignment scheme, i.e., $\gamma\neq 0$, to illustrate that
the 95% confidence intervals constructed according to Theorem 1 have the
desired level while those constructed using Proposition 1 may not.

Table 3 considers the case where the potential outcomes
$(r_{d=0,ij},r_{d=1,ij})$ are continuous and summarizes the coverage of $95\%$
confidence intervals under Proposition 1 and Theorem 1 across $1,000$
simulations in Scenario 1 and Scenario 2. In the absence of a biased IV dose
assignment, i.e., $\gamma=0$, the two IV assignment mechanisms collapse to
one, and the test statistic almost or fully achieves nominal coverage in both
scenarios across a range of sample sizes.

As the degree of bias in the IV dose assignment increases in either scenario,
the test exhibits exceedingly small or zero coverage under Proposition 1. In
such cases, conducting randomization inference for $LB$ overestimates $LB$ and
leads to under-coverage, even though the test tends to be conservative in the
first place. By contrast, the test constructed according to Theorem 1 tends to
achieve perfect coverage under Mechanism I and close to nominal coverage under
Mechanism II, whereby the degree of coverage varies with unit-level treatment
effect heterogeneity.

Two sources of conservativeness account for this trend. First, the test
constructed according to Theorem 1 assumes a worst-case allocation of bias
within each pair; hence, the test will tend to be more conservative for
Mechanism 1 compared to Mechanism 2. Second, even under the worst-case bias
allocation in Mechanism 2, the variance estimator of the test statistic is a
conservative estimator of the true variance in Scenario 2; therefore, the test
still tends to be conservative with a coverage often exceeding its nominal
level. We observe similar trends for the case where the potential outcomes
$(r_{d=0,ij},r_{d=1,ij})$ are binary; see Supplemental Material D for details.

Table 3: Simulation results when the potential outcomes
$(r_{d=0,ij},r_{d=1,ij})$ are continuous under two scenarios: causal null
hypothesis (Scenario 1) and heterogeneous unit-level treatment effect
(Scenario 2). We report the average coverage of $95\%$ confidence intervals
under Proposition 1 and Theorem 1 across $1,000$ simulations.

|  |  | Scenario 1 | Scenario 2
---|---|---|---|---
|  |  | 95% CI Coverage | 95% CI Coverage | 95% CI Coverage | 95% CI Coverage
| $I$ | $\gamma$ | Proposition 1 | Proposition 3 | Proposition 1 | Proposition 3
Mechanism I | 100 | 0.000 | 0.937 | 0.946 | 0.953 | 0.944
| 500 | 0.000 | 0.949 | 0.951 | 0.961 | 0.958
| 1000 | 0.000 | 0.955 | 0.945 | 0.972 | 0.971
| 2000 | 0.000 | 0.954 | 0.941 | 0.965 | 0.959
| 100 | 0.025 | 0.797 | 0.997 | 0.826 | 0.998
| 500 | 0.025 | 0.349 | 1.000 | 0.416 | 1.000
| 1000 | 0.025 | 0.100 | 1.000 | 0.142 | 1.000
| 2000 | 0.025 | 0.002 | 1.000 | 0.005 | 1.000
| 100 | 0.050 | 0.537 | 1.000 | 0.611 | 1.000
| 500 | 0.050 | 0.020 | 1.000 | 0.030 | 1.000
| 1000 | 0.050 | 0.000 | 1.000 | 0.000 | 1.000
| 2000 | 0.050 | 0.000 | 1.000 | 0.000 | 1.000
Mechanism II | 100 | 0.000 | 0.953 | 0.951 | 0.963 | 0.966
| 500 | 0.000 | 0.953 | 0.945 | 0.960 | 0.954
| 1000 | 0.000 | 0.945 | 0.963 | 0.961 | 0.964
| 2000 | 0.000 | 0.947 | 0.951 | 0.963 | 0.957
| 100 | 0.025 | 0.405 | 0.932 | 0.522 | 0.973
| 500 | 0.025 | 0.002 | 0.943 | 0.007 | 0.996
| 1000 | 0.025 | 0.000 | 0.946 | 0.000 | 0.998
| 2000 | 0.025 | 0.000 | 0.940 | 0.000 | 1.000
| 100 | 0.050 | 0.028 | 0.893 | 0.069 | 0.980
| 500 | 0.050 | 0.000 | 0.922 | 0.000 | 0.998
| 1000 | 0.050 | 0.000 | 0.920 | 0.000 | 1.000
| 2000 | 0.050 | 0.000 | 0.944 | 0.000 | 1.000

##  7 Revisiting the NICU study

###  7.1 Assessing the randomization assumption

As a first step towards a randomization-based outcome analysis for the design
M1, we assessed the randomization assumption using the Classification
Permutation Test (CPT) (Gagnon-Bartsch and Shem-Tov,, 2019). Specifically, we
tested the null hypothesis that the excess travel times of mothers within each
matched pair are randomly assigned. When implementing the CPT, we fit a
logistic regression model to predict which mother within each matched pair had
a larger excess travel time based on the observed covariates and used in-
sample classification accuracy rate as our test statistic. We then obtained
the permutation-based null distribution of the test statistic based on $500$
permutations of the IV doses. Finally, we computed an exact p-value by
comparing the observed test statistic to the reference distribution. Figure 2a
shows this null distribution with the observed test statistic superimposed. A
large value of the test statistic indicates that the covariates, even after
matching, are still predictive of the IV dose magnitude within each matched
pair and serves as evidence against the randomization assumption. Thus, the
randomization assumption for the design M1 was rejected (p-value $<0.001$).

To what extent was the randomization assumption violated? Towards
understanding this question, we proceeded to test the following biased
randomization assumption:
$H_{0,\Gamma}\mathrel{\mathop{\ordinarycolon}}\frac{1}{1+\Gamma}\leq\pi_{i}\leq\frac{\Gamma}{1+\Gamma}$
for all $i$ (Chen et al.,, 2023). The null hypothesis $H_{0,\Gamma}$ implies
that the IV-dose-dependent odds ratio $\Gamma_{i}$ for every matched pair $i$
within M1 is uniformly bounded between $1/\Gamma$ and $\Gamma$ and hence the
IV assignment probability $\pi_{i}$ is uniformly bounded between
$1/(1+\Gamma)$ and $\Gamma/(1+\Gamma)$. Using Chen et al.,’s (2023) sample-
splitting CPT, we rejected the null hypothesis $H_{0,\Gamma}$ for $\Gamma$ as
large as $1.17$; therefore, we have evidence that at least one matched pair
exhibits bias in the IV dose assignment odds $\Gamma_{i}$ at a magnitude
larger than $1.17$. This implies that the smallest $\gamma$ in
$\mathcal{M}_{\gamma}$ not falsified by the observed data and the diagnostic
test equals
$\gamma=\log(1.17)/(\max_{i}\widetilde{Z}^{\text{obs}}_{i1}\vee\widetilde{Z}^{\text{obs}}_{i2}-\widetilde{Z}^{\text{obs}}_{i1}\wedge\widetilde{Z}^{\text{obs}}_{i2})=\log(1.17)/129=0.0012$
for M1. Figure 2b displays the distribution of $\Gamma_{i}$’s across all
matched pairs in M1 when $\gamma=0.0012$. This value of $\gamma$ and its
associated sensitivity analysis model $\mathcal{M}_{\gamma}$ served as the
basis for our biased-randomization-based primary outcome analysis.

(a)  (b)

Figure 2: Panel (a): Distribution of the in-sample classification accuracy
rate under the null hypothesis for M1. The vertical bar denotes the observed
test statistic. Panel (b): density plot of $\Gamma_{i}$’s with $\gamma=0.0012$
in the design M1.

###  7.2 Primary analysis of the matched sample

We conducted a biased-randomization-based primary analysis for the matched
sample M1 under $\mathcal{M}_{\gamma}$ with $\gamma=0.0012$. We considered two
outcomes of interest. The primary outcome of interest was infant mortality
rate. We also considered a composite outcome that combines infant mortality
rate and length of stay in the intensive care unit (ICU) as an overall
measurement of health care quality; such a composite outcome is referred to as
a “placement of death” approach (Lin et al.,, 2017) and used in many clinical
settings (see, e.g., Courtright et al., 2024). Specifically, this composite
outcome equals the length of stay in the ICU (LOSI) if the baby survived and
the top 1% quantile of survivor LOSI if the baby died. Table 4 summarizes the
sample size, compliance rate, mortality, LOSI among survivors, and the
composite outcome in the near group (mothers with the smaller excess travel
time in each pair) and far group (mothers with the larger excess travel time
in each pair).

Table 4: Level of encouragement, proportion attending a high-level NICU, and
mortality in the near (mothers with the smaller excess travel time in each
pair) group and far (mothers with the larger excess travel time in each pair)
group.

|  | M1 | Inferential Results |
---|---|---|---|---
|  | Near Group | Far Group | Pt. Est. | 95% CI |
Primary analysis: | Sample size | 81,766 | 81,766 |  |  |
Entire cohort | Excess travel time (min) | 2.07 | 27.1 |  |  |
| High-level NICU (%) | 71.9 | 35.5 |  |  |
| Death (%) | 2.1 | 2.6 |  |  |
|  Effect ratio, $\lambda$ |  |  | 1.28 | [0.69, 1.90] |
|  SATE, $\kappa$ |  |  |  | Fig 3(a) |
| LOSI among survivors (days) | 8.16 | 7.65 |  |  |
| Composite outcome | 9.64 | 9.50 |  |  |
|  Effect ratio, $\lambda$ |  |  | -0.38 | [-1.34, 0.59] |
|  SATE, $\kappa$ |  |  |  | Fig 3(b) |
Subgroup Analysis I: | Sample size | 13,454 | 13,454 |  |  |
Black | Excess travel time (min) | -0.39 | 9.25 |  |  |
| High-level NICU (%) | 82.96 | 64.64 |  |  |
| Death (%) | 3.54 | 3.84 |  |  |
|  Effect ratio, $\lambda$ |  |  | 1.62 | [-1.16, 4.50] |
|  SATE, $\kappa$ |  |  |  | Fig 4(a) |
| LOSI among survivors (days) | 9.95 | 9.58 |  |  |
| Composite outcome | 12.81 | 12.72 |  |  |
|  Effect ratio, $\lambda$ |  |  | -0.48 | [-1.34, 0.59] |
|  SATE, $\kappa$ |  |  |  | Fig S1 |
Subgroup Analysis II: | Sample size | 6,863 | 6,863 |  |  |
Non-Black, Age $\geq 35$, | Excess travel time (min) | 1.95 | 26.01 |  |  |
Gestational age $\leq 36$ | High-level NICU (%) | 77.78 | 45.32 |  |  |
| Death (%) | 2.77 | 3.77 |  |  |
|  Effect ratio, $\lambda$ |  |  | 3.10 | [0.92, 5.42] |
|  SATE, $\kappa$ |  |  |  | Fig 4(b) |
| LOSI among survivors (days) | 11.84 | 11.36 |  |  |
| Composite outcome | 14.15 | 14.58 |  |  |
|  Effect ratio, $\lambda$ |  |  | 1.32 | [-1.96, 4.74] |
|  SATE, $\kappa$ |  |  |  | Fig S1 |
Subgroup Analysis III: | Sample size | 32,294 | 32,294 |  |  |
Non-Black, Age $\leq 35$, | Excess travel time (min) | 2.54 | 30.86 |  |  |
Gestational age $\leq 36$ | High-level NICU (%) | 71.29 | 37.40 |  |  |
| Death (%) | 2.31 | 2.85 |  |  |
|  Effect ratio, $\lambda$ |  |  | 2.02 | [0.85, 3.26] |
|  SATE, $\kappa$ |  |  |  | Fig 4(b) |
| LOSI among survivors (days) | 11.58 | 10.66 |  |  |
| Composite outcome | 13.42 | 13.72 |  |  |
|  Effect ratio, $\lambda$ |  |  | -0.81 | [-2.78, 1.21] |
|  SATE, $\kappa$ |  |  |  | Fig S1 |
Subgroup Analysis IV: | Sample size | 29,155 | 29,155 |  |  |
Gestational age $>$ 36 | Excess travel time (min) | 2.70 | 31.43 |  |  |
| High-level NICU (%) | 65.87 | 17.73 |  |  |
| Death (%) | 0.23 | 0.41 |  |  |
|  Effect ratio, $\lambda$ |  |  | 0.35 | [0.10, 0.61] |
|  SATE, $\kappa$ |  |  |  | Fig 4(b) |
| LOSI among survivors (days) | 2.72 | 2.60 |  |  |
| Composite outcome | 2.73 | 2.62 |  |  |
|  Effect ratio, $\lambda$ |  |  | -0.24 | [-0.42, -0.07] |
|  SATE, $\kappa$ |  |  |  | Fig S1 |

We first made inference for the effect ratio estimand under
$\mathcal{M}_{\gamma}$ with $\gamma=0.0012$ using Proposition 3. The effect
ratio of the mortality outcome (low-level NICU vs. high-level NICU) was
estimated to be $1.28\%$ with a $95\%$ confidence interval of
$[0.69\%,1.90\%]$, providing evidence that delivering at a low-level NICU
increased infants’ mortality rate among the $71.9\%-35.5\%=36.4\%$ compliers
in M1. Such an estimated treatment effect on compliers is substantial, almost
half of the overall mortality rate in the entire study cohort. On the other
hand, the effect on the composite outcome was estimated to be $-0.38$ among
compliers with a $95\%$ confidence interval of $[-1.34,0.59]$, not providing
evidence for an effect on the composite outcome.

(a) Mortality rate  (b) Composite outcome

Figure 3:  $95\%$ confidence interval plot of $\kappa$ for (a) the mortality
rate with $K_{0}=0$ and $0\leq K_{1}\leq 0.03$ and (b) for the composite
outcome with $0\leq-K_{0}=K_{1}=K\leq 2$, under model $\mathcal{M}_{\gamma}$
with $\gamma=0.0012.$ The dotted lines represent the point estimates while the
solid lines represent the one-sided confidence limits of UB (in blue) and LB
(in red), respectively.

We next considered inferring the sample average treatment effect for the
entire matched cohort M1 under $\mathcal{M}_{\gamma}$ with $\gamma=0.0012$
based on Theorem 1. The left panel of Figure 3 plots the 95% confidence
intervals of the mortality outcome for $81,766\times 2=163,532$ participants
in M1 against different hypothesized upper bound $K_{1}$ on the unidentified
SATE among non-compliers. We assumed that compared to a low-level NICU,
delivering at a high-level NICU would not hurt non-compliers so $K_{0}=0$. For
instance, the 95% confidence interval for the SATE was estimated to be
$[0.25\%,1.32\%]$ when $K_{1}=1\%$ and $[0.25\%,1.97\%]$ when $K_{1}=2\%$.
Results in Figure 3 simultaneously account for the sampling variability, the
bias in the instrumental variable (controlled by $\mathcal{M}_{\gamma}$), and
the unidentifiability of the SATE among non-compliers (controlled by the
sensitivity parameters $K_{0}$ and $K_{1}$). Finally, the right panel of
Figure 3 summarizes the results for the composite outcome for some selected
$K_{0}$ and $K_{1}$ values.

###  7.3 Subgroup analyses: Does the treatment effect vary by risk groups?

To further examine treatment effect heterogeneity, we perform an additional
biased-randomization-based primary outcome analysis for the matched sample M1
under Proposition 3 and model $\mathcal{M}_{\gamma}$ with $\gamma=0.0012$ on
four subgroups of mothers: (i) black mothers; (ii) non-black mothers who had
reached 35 and whose gestational age were at most 30 weeks (non-black, high
risk); (iii) non-black mothers who had not reached 35 and whose gestational
age were at most 30 weeks (non-black, medium risk); and (iv) non-black mothers
whose gestational age exceeded 30 weeks (non-black, low risk). Table 4
summarizes the information of sample size, compliance rate, and outcomes of
interest in the near and far groups for these subgroups in the design M1.

For the subgroup of black mothers, the estimated effect ratio (low-level NICU
vs. high-level NICU) is $1.62\%$ with a $95\%$ confidence interval of
$[-1.16\%,4.50\%].$ Thus, we did not find evidence for a treatment effect
among the $82.96\%-64.64\%=18.32\%$ compliers in $\textsf{M1}.$ Figure 4a
further plots the 95% CIs for the sample average treatment effect among black
mothers against different levels of $K_{1}$. Because of the low compliance
rate among black mothers (even after strengthening) and a smaller sample size,
both the effect ratio estimate and the partial identification bounds estimates
are less precise compared to the analysis based on the entire study cohort in
the design M1. The results from our analysis acknowledge that excess travel
time yields rather limited information regarding a potential treatment effect
of delivering at a high-level NICU vs. low-level NICU for the black mothers in
Pennsylvania, mostly as a consequence of the geography of the state.

Among the non-black mothers, the effect ratio (low-level NICU vs. high-level
NICU) was estimated to be $3.10\%$ with a $95\%$ confidence interval of
$[0.92\%,5.42\%]$ among the non-black, high risk subgroup; $2.02\%$ with a
$95\%$ confidence interval of $[0.85\%,3.26\%]$ among the non-black, medium
risk subgroup; and $0.35\%$ with a $95\%$ confidence interval of
$[0.10\%,0.61\%]$ among the non-black, low risk subgroup. Our analyses suggest
that the clinical benefit of delivering at a high-level NICU appeared to be
most pronounced among compliers in the high-risk subgroup, with a risk
difference estimate almost 10-fold that among compliers in the low-risk
subgroup.

Figure 4 further plots the 95% CIs of the sample average treatment effect for
each non-black subgroup against different values of $K_{1}$. Comparing the
partial identification bounds of three non-black risk groups, we have an
impression that the low-risk subgroup appeared to minimally benefit from the
high-level NICU; for instance, a $95\%$ CI for the SATE equals
$[0.05\%,0.39\%]$ when $K_{1}=0.2\%$, $[0.05\%,0.50\%]$ when $K_{1}=0.4\%$,
and $[0.05\%,0.61\%]$ when $K_{1}=0.6\%$. Interestingly, among the four
subgroups, the partial identification bounds are the narrowest for the non-
black, low-risk subgroup, as a consequence of the large sample size
($n=29,155\times 2=58,310$) and high compliance rate ($65.9\%-17.7\%=48.2\%$)
in this subgroup.

Finally, we repeated the subgroup analysis for the composite outcome. We did
not find any difference in the composite outcome among compliers in the black,
high-risk non-black, and the medium-risk non-black subgroups; however, the
effect ratio was estimated to be $-0.24$ with a $95\%$ confidence interval of
$[-0.42,-0.07]$ for low-risk non-black mothers, providing some additional
evidence that the triage system may consider sending low-risk mothers to low-
level NICUs, especially under a resource constraint for high-level NICUs. For
a $95\%$ confidence interval plot of the sample average treatment effect for
each subgroup; see Supplemental Material C.

(a) Black mothers  (b) Non-black, high-risk mothers  (c) Non-black, medium-
risk mothers  (d) Non-black, low-risk mothers

Figure 4:  $95\%$ confidence interval plot of $\kappa$ for (a) the mortality
rate with $K_{0}=0$ and $0\leq K_{1}\leq 0.03$ and (b) for the composite
outcome with $0\leq-K_{0}=K_{1}=K\leq 2$, under model $\mathcal{M}_{\gamma}$
with $\gamma=0.0012.$ The dotted lines represent the point estimates while the
solid lines represent the one-sided confidence limits of UB (in blue) and LB
(in red), respectively.

##  8 Discussion

Strengthening an instrumental variable when designing an observational study
is a novel and useful tool currently underused in empirical comparative
effectiveness research. It also has profound consequences. In their
discussion, Baiocchi et al., (2010, Section 5) discussed aspects that would
change as one switch between an unstrengthened IV design and a strengthened IV
design, or between two IV designs that strengthen an IV to different extent.
This article builds upon the discussion in Baiocchi et al., (2010) and related
critiques in Deaton, (2009) and proposes strategies to better reconcile
different strengthened IV designs. One promising strategy is to complement an
effect ratio estimate with partial identification bounds on the sample average
treatment effect. Unlike the effect ratio estimate that targets the SATE among
compliers, a study-design-dependent estimand, the SATE on the entire matched
cohort is agnostic about how one forms the matched pairs or matched sets. Just
like strengthening an IV improves the precision of the effect ratio estimate,
strengthening an IV also helps narrow the partial identification bounds by
improving the compliance rate.

We studied the effect of delivering at a high-level NICU versus a low-level
NICU using a comprehensive administrative database in Pennsylvania and a suite
of tools developed in this article, including an improved non-bipartite
matching algorithm that emulates a randomized encouragement trial, inferential
methods for constructing a valid level-$\alpha$ confidence interval for
partial identification bounds, and inferential methods that conducts valid
inference under a more realistic and less conservative, IV-dependent biased
randomization scheme. Our data analyses suggest that delivering at a high-
level NICU decreased preterm infants’ mortality rate. The effect was found to
be heterogeneous, with a most pronounced effect for non-black, high-risk
mothers but marginal for the non-black, low-risk mothers, providing some
evidence that a triage system may not need to send mothers and their preterm
babies to high-level NICUs. Interestingly, although black mothers were at
elevated risk and could potentially benefit from delivering at a high-level
NICU, our analysis did not find evidence for this, most likely because the
excess travel time was a weak IV, even after strengthening, for black mothers
in Pennsylvania.

In this study, we pre-specified several subgroups based on previous clinical
findings (Hansen,, 1986; Yang et al.,, 2014; Mathews et al.,, 2015). With a
binary treatment and in the context of testing Fisher’s sharp null hypothesis,
data-dependent methods that identify promising subgroups have been developed
and deployed in matched observational studies (Hsu et al.,, 2013, 2015; Lee et
al.,, 2018). It is of great interest to develop methods that could identify
subgroups most amenable to being strengthened and with promise for a large
treatment effect.

## References

  * Angrist and Fernandez-Val, (2010) Angrist, J. and Fernandez-Val, I. (2010).  Extrapolate-ing: External validity and overidentification in the late framework.  Technical report, National Bureau of Economic Research.
  * Angrist and Imbens, (1995) Angrist, J. D. and Imbens, G. W. (1995).  Two-stage least squares estimation of average causal effects in models with variable treatment intensity.  Journal of the American Statistical Association, 90(430):431–442.
  * Angrist et al., (1996) Angrist, J. D., Imbens, G. W., and Rubin, D. B. (1996).  Identification of causal effects using instrumental variables.  Journal of the American Statistical Association, 91(434):444–455.
  * Baiocchi et al., (2014) Baiocchi, M., Cheng, J., and Small, D. S. (2014).  Instrumental variable methods for causal inference.  Statistics in Medicine, 33(13):2297–2340.
  * Baiocchi et al., (2010) Baiocchi, M., Small, D. S., Lorch, S., and Rosenbaum, P. R. (2010).  Building a stronger instrument in an observational study of perinatal care for premature infants.  Journal of the American Statistical Association, 105(492):1285–1296.
  * Bennett et al., (2020) Bennett, M., Vielma, J. P., and Zubizarreta, J. R. (2020).  Building representative matched samples with multi-valued treatments in large observational studies.  Journal of Computational and Graphical Statistics, 29(4):744–757.
  * Chen et al., (2023) Chen, K., Heng, S., Long, Q., and Zhang, B. (2023).  Testing biased randomization assumptions and quantifying imperfect matching and residual confounding in matched observational studies.  Journal of Computational and Graphical Statistics, 32(2):528–538.
  * Chen and Zhang, (2023) Chen, S. and Zhang, B. (2023).  Estimating and improving dynamic treatment regimes with a time-varying instrumental variable.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(2):427–453.
  * Cole and Stuart, (2010) Cole, S. R. and Stuart, E. A. (2010).  Generalizing evidence from randomized clinical trials to target populations: the ACTG 320 trial.  American Journal of Epidemiology, 172(1):107–115.
  * Cook et al., (2002) Cook, T. D., Campbell, D. T., and Shadish, W. (2002).  Experimental and quasi-experimental designs for generalized causal inference.  Houghton Mifflin Boston, MA.
  * Courtright et al., (2024) Courtright, K. R., Madden, V., Bayes, B., Chowdhury, M., Whitman, C., Small, D. S., Harhay, M. O., Parra, S., Cooney-Zingman, E., Ersek, M., et al. (2024).  Default palliative care consultation for seriously ill hospitalized patients: A pragmatic cluster randomized trial.  JAMA, 331(3):224–232.
  * Deaton, (2009) Deaton, A. S. (2009).  Instruments of development: Randomization in the tropics, and the search for the elusive keys to economic development.  Technical report, National Bureau of Economic Research.
  * Derigs, (1988) Derigs, U. (1988).  Solving non-bipartite matching problems via shortest path techniques.  Annals of Operations Research, 13(1):225–261.
  * Ding and Lu, (2017) Ding, P. and Lu, J. (2017).  Principal stratification analysis using principal scores.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), pages 757–777.
  * (15) Fogarty, C. B. (2018a).  On mitigating the analytical limitations of finely stratified experiments.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 80(5):1035–1056.
  * (16) Fogarty, C. B. (2018b).  Regression-assisted inference for the average treatment effect in paired experiments.  Biometrika, 105(4):994–1000.
  * Fogarty, (2020) Fogarty, C. B. (2020).  Studentized sensitivity analysis for the sample average treatment effect in paired observational studies.  Journal of the American Statistical Association, 115(531):1518–1530.
  * Fogarty et al., (2021) Fogarty, C. B., Lee, K., Kelz, R. R., and Keele, L. J. (2021).  Biased encouragements and heterogeneous effects in an instrumental variable study of emergency general surgical outcomes.  Journal of the American Statistical Association, 116(536):1625–1636.
  * Gagnon-Bartsch and Shem-Tov, (2019) Gagnon-Bartsch, J. and Shem-Tov, Y. (2019).  The classification permutation test: A flexible approach to testing for covariate imbalance in observational studies.  Annals of Applied Statistics, 13(3):1464–1483.
  * Haavelmo, (1943) Haavelmo, T. (1943).  The statistical implications of a system of simultaneous equations.  Econometrica, Journal of the Econometric Society, pages 1–12.
  * Hansen, (1986) Hansen, J. P. (1986).  Older maternal age and pregnancy outcome: a review of the literature.  Obstetrical & Gynecological Survey, 41(11):726–742.
  * Heckman and Urzua, (2010) Heckman, J. J. and Urzua, S. (2010).  Comparing iv with structural models: What simple iv can and cannot identify.  Journal of Econometrics, 156(1):27–37.
  * Heng et al., (2023) Heng, S., Zhang, B., Han, X., Lorch, S. A., and Small, D. S. (2023).  Instrumental variables: to strengthen or not to strengthen?  Journal of the Royal Statistical Society Series A: Statistics in Society, page qnad075.
  * Hernán and Robins, (2006) Hernán, M. A. and Robins, J. M. (2006).  Instruments for causal inference: an epidemiologist’s dream?  Epidemiology, pages 360–372.
  * Hsu et al., (2013) Hsu, J. Y., Small, D. S., and Rosenbaum, P. R. (2013).  Effect modification and design sensitivity in observational studies.  Journal of the American Statistical Association, 108(501):135–148.
  * Hsu et al., (2015) Hsu, J. Y., Zubizarreta, J. R., Small, D. S., and Rosenbaum, P. R. (2015).  Strong control of the familywise error rate in observational studies that discover effect modification by exploratory methods.  Biometrika, 102(4):767–782.
  * Imai, (2008) Imai, K. (2008).  Variance identification and efficiency analysis in randomized experiments under the matched-pair design.  Statistics in Medicine, 27(24):4857–4873.
  * Imbens, (2010) Imbens, G. W. (2010).  Better late than nothing: Some comments on deaton (2009) and heckman and urzua (2009).  Journal of Economic Literature, 48(2):399–423.
  * Imbens and Rosenbaum, (2005) Imbens, G. W. and Rosenbaum, P. R. (2005).  Robust, accurate confidence intervals with a weak instrument: quarter of birth and education.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 168(1):109–126.
  * Jo and Stuart, (2009) Jo, B. and Stuart, E. A. (2009).  On the use of propensity scores in principal causal effect estimation.  Statistics in Medicine, 28(23):2857–2875.
  * Joffe, (2011) Joffe, M. (2011).  Principal stratification and attribution prohibition: good ideas taken too far.  The International Journal of Biostatistics, 7(1).
  * Lee et al., (2018) Lee, K., Small, D. S., Hsu, J. Y., Silber, J. H., and Rosenbaum, P. R. (2018).  Discovering effect modification in an observational study of surgical mortality at hospitals with superior nursing.  Journal of the Royal Statistical Society Series A: Statistics in Society, 181(2):535–546.
  * Lei and Ding, (2021) Lei, L. and Ding, P. (2021).  Regression adjustment in completely randomized experiments with a diverging number of covariates.  Biometrika, 108(4):815–828.
  * Li and Ding, (2020) Li, X. and Ding, P. (2020).  Rerandomization and regression adjustment.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(1):241–268.
  * Lin, (2013) Lin, W. (2013).  Agnostic notes on regression adjustments to experimental data: Reexamining freedman’s critique.  Annals of Applied Statistics.
  * Lin et al., (2017) Lin, W., Halpern, S. D., Prasad Kerlin, M., and Small, D. S. (2017).  A “placement of death” approach for studies of treatment effects on icu length of stay.  Statistical Methods in Medical Research, 26(1):292–311.
  * Lorch et al., (2012) Lorch, S. A., Baiocchi, M., Ahlberg, C. E., and Small, D. S. (2012).  The differential impact of delivery hospital on the outcomes of premature infants.  Pediatrics, 130(2):270–278.
  * Lu et al., (2011) Lu, B., Greevy, R., Xu, X., and Beck, C. (2011).  Optimal nonbipartite matching and its statistical applications.  The American Statistician, 65(1):21–30.
  * Lu et al., (2001) Lu, B., Zanutto, E., Hornik, R., and Rosenbaum, P. R. (2001).  Matching with doses in an observational study of a media campaign against drug abuse.  Journal of the American Statistical Association, 96(456):1245–1253.
  * Mathews et al., (2015) Mathews, T., MacDorman, M. F., and Thoma, M. E. (2015).  Infant mortality statistics from the 2013 period linked birth/infant death data set.
  * Newhouse and McClellan, (1998) Newhouse, J. P. and McClellan, M. (1998).  Econometrics in outcomes research: the use of instrumental variables.  Annual Review of Public Health, 19(1):17–34.
  * Pu and Zhang, (2021) Pu, H. and Zhang, B. (2021).  Estimating optimal treatment rules with an instrumental variable: A partial identification learning approach.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(2):318–345.
  * Rigdon et al., (2018) Rigdon, J., Baiocchi, M., and Basu, S. (2018).  Near-far matching in r: the nearfar package.  Journal of Statistical Software, 86(CS-5).
  * Rosenbaum and Rubin, (2023) Rosenbaum, P. and Rubin, D. (2023).  Propensity scores in the design of observational studies for causal effects.  Biometrika, 110(1):1–13.
  * Rosenbaum, (1987) Rosenbaum, P. R. (1987).  Sensitivity analysis for certain permutation inferences in matched observational studies.  Biometrika, 74(1):13–26.
  * Rosenbaum, (1989) Rosenbaum, P. R. (1989).  Sensitivity analysis for matched observational studies with many ordered treatments.  Scandinavian Journal of Statistics, pages 227–236.
  * (47) Rosenbaum, P. R. (2002a).  Covariance adjustment in randomized experiments and observational studies.  Statistical Science, 17(3):286–327.
  * (48) Rosenbaum, P. R. (2002b).  Observational Studies.  Springer.
  * Rosenbaum, (2010) Rosenbaum, P. R. (2010).  Design of Observational Studies, volume 10.  Springer.
  * Rosenbaum, (2020) Rosenbaum, P. R. (2020).  Modern algorithms for matching in observational studies.  Annual Review of Statistics and Its Application, 7:143–176.
  * Rosenbaum and Rubin, (1985) Rosenbaum, P. R. and Rubin, D. B. (1985).  Constructing a control group using multivariate matched sampling methods that incorporate the propensity score.  The American Statistician, 39(1):33–38.
  * Rubin, (1980) Rubin, D. B. (1980).  Randomization analysis of experimental data: The fisher randomization test comment.  Journal of the American Statistical Association, 75(371):591–593.
  * Rubin, (1986) Rubin, D. B. (1986).  Statistics and causal inference: Comment: Which ifs have causal answers.  Journal of the American Statistical Association, 81(396):961–962.
  * Silber et al., (2001) Silber, J. H., Rosenbaum, P. R., Trudeau, M. E., Even-Shoshan, O., Chen, W., Zhang, X., and Mosher, R. E. (2001).  Multivariate matching and bias reduction in the surgical outcomes study.  Medical Care, 39(10):1048–1064.
  * Small and Rosenbaum, (2008) Small, D. S. and Rosenbaum, P. R. (2008).  War and wages: the strength of instrumental variables and their sensitivity to unobserved biases.  Journal of the American Statistical Association, 103(483):924–933.
  * Stuart et al., (2011) Stuart, E. A., Cole, S. R., Bradshaw, C. P., and Leaf, P. J. (2011).  The use of propensity scores to assess the generalizability of results from randomized trials.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 174(2):369–386.
  * Su and Ding, (2021) Su, F. and Ding, P. (2021).  Model-assisted analyses of cluster-randomized experiments.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(5):994–1015.
  * Swanson et al., (2018) Swanson, S. A., Hernán, M. A., Miller, M., Robins, J. M., and Richardson, T. S. (2018).  Partial identification of the average treatment effect using instrumental variables: review of methods for binary instruments, treatments, and outcomes.  Journal of the American Statistical Association, 113(522):933–947.
  * Wang and Tchetgen Tchetgen, (2018) Wang, L. and Tchetgen Tchetgen, E. (2018).  Bounded, efficient and multiply robust estimation of average treatment effects using instrumental variables.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 80(3):531–550.
  * Yang et al., (2014) Yang, F., Lorch, S. A., and Small, D. S. (2014).  Estimation of causal effects using instrumental variables with nonignorable missing covariates: application to effect of type of delivery nicu on premature infants.  Annals of Applied Statistics, 8(1):48–73.
  * Yu, (2020) Yu, R. (2020).  Evaluating and improving a matched comparison of antidepressants and bone density.  Biometrics.
  * Zhang, (2023) Zhang, B. (2023).  Efficient Algorithms for Building Representative Matched Pairs with Enhanced Generalizability.  Biometrics, 79(4):3981–3997.
  * Zhang et al., (2022) Zhang, B., Heng, S., MacKay, E. J., and Ye, T. (2022).  Bridging preference-based instrumental variable studies and cluster-randomized encouragement experiments: Study design, noncompliance, and average cluster effect ratio.  Biometrics, 78(4):1639–1650.
  * Zhang et al., (2023) Zhang, B., Mackay, E. J., and Baiocchi, M. (2023).  Statistical matching and subclassification with a continuous dose: characterization, algorithm, and application to a health outcomes study.  The Annals of Applied Statistics, 17(1):454–475.
  * Zhao and Ding, (2022) Zhao, A. and Ding, P. (2022).  Regression-based causal inference with factorial experiments: estimands, model specifications and design-based properties.  Biometrika, 109(3):799–815.