# Using Joint Variable Importance Plots to Prioritize Variables in Assessing
the Impact of Glyburide on Adverse Birth Outcomes

Lauren D<EMAIL_ADDRESS>[ Yeyi<EMAIL_ADDRESS>[ Amanda L<EMAIL_ADDRESS>[ Rana F.
<EMAIL_ADDRESS>[ Samuel D.
<EMAIL_ADDRESS>[ Division of Biostatistics, University of
California, Berkeleypresep=, ]e1  Kaiser Permanente Northern California
Division of Researchpresep=, ]e2,e3,e4  Department of Statistics, University
of California, Berkeleypresep=, ]e5

###### Abstract

The only pharmacologic treatment for gestational diabetes (GDM) approved by
U.S. Food and Drug Administration is insulin. However, due to improved ease of
use and lower cost, oral antidiabetic medications, such as glyburide, have
been prescribed more commonly than insulin. We investigate glyburide’s impact
on two adverse perinatal outcomes compared to medical nutritional therapy, the
universal first-line therapy, in an observational study of a large,
population-based cohort. At the design stage, we employ matching to select
comparable treated (received glyburide) and controls (received medical
nutritional therapy). Multiple background variables were associated with GDM
treatment modality and perinatal outcomes; however, there is ambiguity about
which potential variables should be selected for matching. Standard selection
methods based on treatment imbalance neglect the variable’s relationship with
the outcome. Thus, we propose the joint variable importance plot (jointVIP) to
guide variable prioritization for this study. These plots enhance variable
comparisons via bias curves derived under the omitted variable bias framework
and can translate into recommended values for tuning parameters in existing
methods. After forming matched pairs, we conduct inference and perform
sensitivity analyses to assess the potential role of unmeasured confounding.
Understanding glyburide safety can better inform future pharmacologic
treatment strategies to manage GDM.

Observational Study,

Research Design,

Causality,

Gestational Diabetes,

Electronic Health Records,

Confounding Variables,

###### keywords:

, and

##  1 Introduction

###  1.1 The motivating problem: treatment for gestational diabetes

Gestational diabetes mellitus (GDM) is one of the most common pregnancy
metabolic complications in the U.S., with increasing prevalence above 35% over
the past decades, fueling the type 2 diabetes epidemic [Ferrara et al. (2004);
Zhu and Zhang (2016)]. Individuals with GDM and their infants are at higher
risk of adverse health outcomes than their non-GDM exposed counterparts. The
American Diabetes Association (2021) and American College of Obstetricians and
Gynecologists (2018) recommend insulin therapy as the standard pharmacologic
therapy for GDM. Although insulin is the only approved pharmacologic treatment
for GDM according to the U.S. Food and Drug Administration, use of oral
antidiabetic medications, such as glyburide and metformin, have significantly
increased from 7.4% in 2000 to 64.5% in 2011 in the U.S. due to lower cost and
better ease of use [Castillo et al. (2014); Goetzl and Wilkins (2002)]. The
safety of glyburide, however, remains a contention due to potential transfer
to the fetus through the placenta [American College of Obstetricians and
Gynecologists (2018)]. While some studies indicated an increased risk of
perinatal complications associated with glyburide use, others found no
difference [Rochon et al. (2006); Behrashi et al. (2016)]. The question
remains: does glyburide increase the risk of adverse perinatal outcomes in
real-world settings?

Since this is a question about a treatment effect, we use the formal framework
of causal inference to investigate it. We compare glyburide against medical
nutritional therapy (MNT), the universal first-line therapy for GDM. We assess
possible increased risk of two selected outcomes, Cesarean section delivery
(C-section) and admission of the newborn to the neonatal intensive care unit
(NICU) which may not only result in remarkable patient and health care level
burden but also increase risk of adverse maternal and neonatal outcomes in
later life, respectively.

###  1.2 Review: observational study design

Consider a hypothetical randomized control trial (RCT) in which patients
diagnosed with GDM are randomly assigned to receive either glyburide and MNT
together or MNT alone. The two comparison groups in this RCT would tend (in
large samples) to be balanced, or similar, on observed and unobserved baseline
attributes. The remaining differences in adverse perinatal outcomes could thus
be safely attributed to the impact of glyburide. However, RCTs are expensive
and potentially risky for pregnant individuals due to the unknown effect of
treatment on the individual and their offspring(s). In contrast, observational
studies based on abundant pre-existing electronic health records face fewer
practical barriers but greater statistical ones. There is no guarantee that
patients receiving glyburide in real-world clinical settings are comparable to
those receiving MNT alone. In these settings, one approach is to create new
comparison groups to enable treatment effect estimation amongst similar
individuals. Matching [Stuart (2010)] and weighting [Ben-Michael et al.
(2021)] methods are commonly used to achieve this by increasing similarity in
observed baseline variables through taking appropriate subsets or transforming
the raw groups. By analogy to the ideal RCT, in which a study design would be
developed using baseline attributes before any outcomes are collected, the
matching or weighting step is known as the design stage of the observational
study. The ”design principle” dictates that study outcomes may not be used
during the design stage; this protects researchers against biases that arise
when the same outcome data are used both to select a design (or model) and to
measure a causal effect [Rubin (2007)]. It is worth noting though that
matching and weighting adjust for observed differences between groups only, in
contrast to the ideal RCT design that addresses both observed and unobserved
differences. Sensitivity analyses must be conducted in observational studies
after outcome analysis to address possible unobserved group differences.

Under the aforementioned design principle, researchers have access only to
background variables and the treatment variable when conducting matching or
weighting. Thus, traditionally, researchers focus on eliminating treatment-
control imbalances in the background variables [Hainmueller (2012); Diamond
and Sekhon (2013)]. To avoid post-treatment adjustments during the design
stage, these background variables are measured prior to the treatment and
outcome. It is common to summarize each variable’s imbalance using
standardized mean differences, presented in a balance table or a graphical
form, such as a Love plot [Ahmed et al. (2006); Stuart et al. (2011); Greifer
(2021); Hansen and Bowers (2008); Rosenbaum and Rubin (1985)]. Matching or
weighting methods are then tailored to ensure that the largest imbalances are
eliminated. While balancing all background variables perfectly across groups
is sufficient to remove observed sources of confounding, this is often
impossible in datasets with a large number of measured variables. Prioritizing
variables for matching or weighting according to their imbalance ignores
important information about the role of each variable in the outcome model. In
particular, variables that are strongly related to treatment, but unrelated to
outcomes, are not confounders, while variables with only moderate imbalance
may be strong confounders if they are strongly associated with the outcome.
When choosing which background variables to prioritize for adjustment, solely
focusing on the treatment imbalance can risk ignoring important variables that
should take precedence due to their outcome importance. To avoid violating the
design principle, recent proposals discuss utilizing a separate sample (or
pilot sample), often only consisting of controls, to fit a prognostic score
(i.e., fitting an outcome model using this separate sample) during the design
stage [Hansen (2008); Antonelli et al. (2018); Aikens, Greaves and Baiocchi
(2020)]. However, these approaches are not yet widely used.

###  1.3 Proposed solution: joint outcome-treatment variable importance plots

We propose selecting high-priority variables for adjustment during the design
stage using the joint outcome-treatment variable importance plot (jointVIP).
The jointVIP represents each variable in two dimensions: one dimension
describing treatment model importance as measured by the standardized mean
difference (SMD), and one dimension describing outcome model importance,
measured by outcome correlation. To respect the design principle, we use a
pilot sample consisting of external controls. In addition, under a set of
simple working models, the bias incurred by ignoring each variable can be
derived and represented on the plot using bias curves, enhancing opportunities
for variable comparisons. In the study of glyburide safety, the jointVIP
provides clear recommendations for specifying important parameters in a
balance-constrained matching procedure. Resulting matched samples constructed
for multiple outcomes exhibit a high degree of treatment-outcome balance in
important variables and permit intuitive randomization inference and
sensitivity analysis procedures.

The paper is structured as follows. Section 2 introduces the study population
of pregnant individuals with GDM. Section 3 proposes the underlying framework
and construction of the jointVIP. Section 4 describes our design choices based
on jointVIPs to form matched pairs. Section 5 uses randomization inference to
assess the causal impact of glyburide on adverse perinatal outcomes. Section 6
discusses implications of our study results for GDM treatment and broader
potential uses of the jointVIP idea.

##  2 Data

###  2.1 Study population

The study population consists of Kaiser Permanente Northern California (KPNC)
members. KPNC is an integrated healthcare delivery system serving over 4.5
million members. KPNC members are representative of the socio-demographic
status of individuals residing in the geographic areas served [Gordon and Lin
(2016)]. Individuals with GDM are identified through the KPNC Pregnancy
Glucose Tolerance and GDM Registry, which is an active surveillance registry
including laboratory data that are used to inform GDM screening and diagnosis.
Pre-existing type 1 or type 2 diabetes is automatically excluded from this
registry. Our data-only study is approved by the KPNC Institutional Review
Board that waived the requirement for informed consent from participants.
Individuals who are diagnosed with GDM receive MNT as the universal first line
of therapy. Pharmacologic treatment, including oral antidiabetic medications
(glyburide, metformin, or other) and/or insulin, is prescribed in addition to
MNT if glycemic control goals are not met. The decision of which pharmacologic
therapy to recommend is made based on patient-physician discussions weighing
the benefits and risks of each option.

Figure 1: Comparison between the Love plot and the joint outcome-treatment
variable importance plot

BMI: body mass index, C-C: Carpenter-Coustan, C-section: Cesarean section,
GDM: gestational diabetes, OGTT: oral glucose tolerance test, OVB: omitted
variable bias, PI: Pacific Islander, SD: standard deviation, SMD: standardized
mean difference
The Love plot (A) is a common practice in visualizing treatment imbalance
using standard SMD, which calculates using a pooled standard deviation from
both the treated and control groups. The joint variable importance plot (B) is
our proposed method to visualize both treatment imbalance and outcome
correlation marginally building on the OVB framework using a slightly
different OVB-based SMD. The OVB-based SMD uses an external (pilot sample)
control standard deviation.

###  2.2 Descriptive analysis

We conduct our analysis on a dataset covering the years 2011-2021 based on the
delivery date of the index pregnancy complicated by GDM, and we focus on two
selected outcomes: one maternal, C-section delivery, and one neonatal, NICU
admission. These adverse perinatal outcomes are chosen due to significant and
potentially lasting impact on patients and their offspring. We use a separate
dataset from 2007-2010 to inform our observational study design, specifically
to evaluate each variable’s outcome relationship. While both 2007-2010 and
2011-2021 datasets contain GDM individuals, we choose the more recent
2011-2021 data to reflect a timely and current analysis. Individuals with GDM
who received MNT alone constituted our control group while those who received
glyburide as the only pharmacologic therapy in addition to MNT constituted our
treatment group. There are 54 common variables between the 2007-2010 and
2011-2021 datasets that we extracted (including indicators of missing data as
variables). Table 1 summarizes selected background variables (see Appendix
Table LABEL:long for the full data summary). As a pre-processing step to treat
missingness, we impute the missing variables, separately for each year, using
random forests; this algorithm was chosen because of its ability to handle
mixed-data types and potential non-parametric, complex relationship that
exists in this data [Stekhoven and Bühlmann (2012)]. More information about
the pattern of missing values and the imputation procedure are reported in
Appendix Table 6.

We visualize variable comparisons using a Love plot (Figure 1.A) and a
jointVIP for the C-section delivery outcome (Figure 1.B). The Love plot shows
absolute standardized mean differences sorted in order of size. In the Love
plot, gestational age at GDM diagnosis, oral glucose tolerance test for
fasting blood glucose, and GDM diagnosed by the Carpenter-Coustan criteria are
variables with the largest imbalance between the treated and control groups.
In addition to those variables, the jointVIP identifies pre-pregnancy body
mass index (BMI) in the obese category as an important variable primarily due
to its role in the outcome model. Essentially, the jointVIP extends the Love
plot by adding the outcome correlation on the y-axis with slightly modified
definition of SMD based on the omitted variable bias framework. We discuss the
construction and interpretation of the jointVIP in Section 3. For versions of
the Love plot and jointVIP with signed (rather than absolute) quantities
plotted on the two axes, see Appendix Figure 4.

Table 1: Summary of selected baseline variables for pregnant individuals with
gestational diabetes.

2007-2010 2011-2021 2011-2021 Control Control Treated n = 7,526 n = 19,183 n =
10,786 Family history of diabetes = yes (%) 198 (2.6) 1,202 (6.3) 822 (7.6)
OGTT for fasting blood glucose = abnormal (%) 2,165 (28.8) 3,762 (19.6) 4,512
(41.8) GDM severity = severe (%) 775 (10.3) 1,744 (9.1) 1,215 (11.3) GDM
diagnosed by the C-C criteria = yes (%) 7,323 (97.3) 17,303 (90.2) 8,419
(78.1) Age (%) Under 25 552 (7.3) 1,029 (5.4) 349 (3.2) Between 25-29 1,665
(22.1) 3,762 (19.6) 1,866 (17.3) Between 30-34 2,584 (34.3) 7,101 (37.0) 4,165
(38.6) Over 35 2,725 (36.2) 7,291 (38.0) 4,406 (40.8) Gestational age at GDM
diagnosis (mean (SD)) 26.06 (6.10) 26.86 (6.02) 23.53 (7.12) History of
macrosomia = yes (%) 69 (0.9) 145 (0.8) 147 (1.4) History of GDM = yes (%) 856
(11.4) 3,401 (17.7) 2,608 (24.2) Parity (%) 0 3,121 (41.5) 7,611 (39.7) 3,873
(35.9) 1 2,347 (31.2) 6,566 (34.2) 3,994 (37.0) more than 2 2,058 (27.3) 5,006
(26.1) 2,919 (27.1) Pre-pregnancy BMI (%) Underweight 100 (1.3) 391 (2.0) 76
(0.7) Normal 1,921 (25.5) 5,107 (26.6) 1,706 (15.8) Overweight 2,847 (37.8)
6,369 (33.2) 3,361 (31.2) Obese 2,658 (35.3) 7,316 (38.1) 5,643 (52.3)
Race/ethnicity (%) Asian or Pacific Islander 2,919 (38.8) 8,553 (44.6) 4,560
(42.3) Hispanic 2,322 (30.9) 5,013 (26.1) 2,923 (27.1) White 1,602 (21.3)
4,090 (21.3) 2,381 (22.1) Black or African American 315 (4.2) 736 (3.8) 410
(3.8) Other or unknown 368 (4.9) 791 (4.1) 512 (4.7) Pre-pregnancy pre-
diabetes = yes (%) 479 (6.4) 1,802 (9.4) 1,915 (17.8) Glucose challenge test
value (mean (SD)) 169.43 (22.38) 169.71 (22.14) 173.22 (24.32)

BMI: body mass index, C-C: Carpenter-Coustan, GDM: gestational diabetes, OGTT:
oral glucose tolerance test, SD: standard deviation
A normal OGTT fasting blood glucose level is lower than 95 mg/dL. Abnormal
indicates a result higher than normal.

##  3 Joint variable importance plots to prioritize potential confounders

The jointVIP is a novel tool to help design observational studies by selecting
key potential confounders for adjustment (e.g., by matching or weighting). To
construct and motivate jointVIP, we draw on concepts from causal inference.
Let $Z$ be a binary indicator of treatment, with $Z=1$ indicating glyburide
and $Z=0$ indicating MNT. We follow the potential outcomes framework proposed
by Neyman (1923) and Rubin (1974), letting $Y(1)$ represent the counterfactual
outcome that would have been observed if the individual in question had
received glyburide, and letting $Y(0)$ represent the counterfactual outcome
that would have been observed if the same individual received MNT alone (so
that the observed outcome $Y=ZY(1)+(1-Z)Y(0)$). We assume one individual’s
treatment does not affect another’s outcome, and there are no other versions
of treatment. See Imbens and Rubin (2015), Chapter 1.6, for further discussion
on the stable unit treatment value assumption. In addition, we operate under a
working assumption that no unmeasured confounding is present (i.e., potential
outcomes for treated and controls are independent of treatment given measured
background variables), an assumption we later relax in Section 5.

###  3.1 Joint variable importance plot construction

The high-level purpose of the jointVIP is to illustrate two different
dimensions of a variable’s possible role as a confounder – its imbalance, or
association with the treatment variable, and its association with the outcome
– on two numerical axes, with each variable plotted as an individual point. We
now discuss the specific measures of variable importance represented on each
axis.

For treatment model importance, described by the x-axis, we use SMDs, or
differences between the treated mean and the control mean divided by external
control (from the pilot sample) standard deviation for the background variable
in question. SMDs allow intuitive comparisons across variables with very
different scales, including both binary and continuous variables, and are
widely used to assess imbalance. In fact, a variation of the SMD forms the
basis for traditional visualizations such as the Love plot. Thus, using SMD
allows jointVIP to preserve all information typically contained in the Love
plot while also adding new insights.

For outcome model importance, represented on the y-axis, we compute the
Pearson correlation between each variable and the outcome in a pilot sample
consisting only of controls. Pearson correlation is a familiar, bounded
quantity and makes sense not only for relationships between two continuous
variables but also for binary and binary variables (phi coefficient), and
binary and continuous variables (point biserial correlation) [Pearson (1895)].
The outcome relationship is only calculated on the control sample to avoid
having to model treatment effects.

Two quantities that often come up in discussions of confounding are the
propensity score [Rosenbaum and Rubin (1983)] and the prognostic score [Hansen
(2008)], which summarize the joint impact of background variables on the
treatment probability and the baseline outcome, respectively. Although
estimates for these scores are produced by fitting models to other background
variables, the resulting fitted values are also variables and may be
represented in the jointVIP (or indeed in the Love plot or balance table). By
construction, the propensity score tends to have a high treatment imbalance,
and the prognostic score tends to have a high outcome correlation, but the
jointVIP may be helpful in illuminating the other dimension for each type of
score and in helping researchers to decide how useful it is in design. For
example, researchers may not want to adjust for the propensity score if it has
a near-zero relationship with the outcome variable as discussed by Brooks and
Ohsfeldt (2013) and Pimentel, Small and Rosenbaum (2016).

Using the controls from the analysis sample for computing outcome
correlations, as opposed to our recommended pilot sample approach, would
violate the design principle and potentially introduce bias into our treatment
effect estimates. As an example, suppose treated and control samples exhibit
an imbalance on several continuous background variables (with treated
individuals taking on larger values) but the study outcome is independent of
all these variables in the population. If we compute sample outcome
correlations in the analysis control sample and select the variable with the
largest such (positive) sample correlation as the sole criterion for matching,
we are essentially matching on the variable with the largest spurious
correlation (with the random noise in the current sample). Because of the
imbalance, the matching algorithm will systematically select controls with
large values for the spuriously correlated variable. Hence, the result will
have large positive outcome errors that introduce positive bias into the
average outcome for the matched controls. For an almost-identical example in
the context of prognostic score estimation with additional discussion and
simulation results, see Hansen (2008), and for a discussion of similar
pitfalls arising in adjustments for randomized trials see Abadie, Chingos and
West (2018).

The bottom line is that the pilot sample must not overlap with the main
analysis sample, although it must preserve the outcome-variable relationships
present in the analysis sample to be useful for jointVIP construction. Such a
pilot sample may be created by selecting a small (10-20%) portion of the
control sample at random from the full control group. To ensure the pilot
sample is drawn from the portion of the control space most relevant for the
observational study, Aikens, Greaves and Baiocchi (2020) suggest an
alternative strategy: conduct an initial round of matching on a standard
Mahalanobis distance to pair each treated subject to two controls, then select
one control from each set at random to construct the pilot sample. Our study
instead uses a separate, historical dataset (collected from 2007-2010) as the
pilot sample, allowing us to preserve the entire analysis sample (collected
from 2011-2021).

###  3.2 Addition of bias curves for variable comparison

Plotting points by SMD and outcome correlation gives a raw version of the
jointVIP that presents much valuable information. However, it can be difficult
in this version of the plot to determine the relative importance of two
variables represented by two distant points, one with a high outcome
correlation and low SMD, and the other with a high SMD and low outcome
correlation. A good way to approach such comparisons is to think about the
relative sizes of the biases contributed by ignoring each variable, since our
ultimate goal is to avoid biases in treatment effect estimation. Inspired by
Cinelli and Hazlett (2020) and Soriano et al. (2021), we plot quantities which
correspond to bias estimates for the treatment effect under a simple set of
statistical models and which take on the convenient form of hyperbolic curves
on the jointVIP. These curves, which arise from the classical omitted variable
bias (OVB) framework, provide valuable context about the relative amounts of
potential bias contributed by distant variables, making comparisons easier.

For any background variable $X_{j}$ with $j\in{1,...,J}$), consider the sample
least-squares fit of $Y$ on $X_{j}$ and $Z$:

$\displaystyle Y$ $\displaystyle=Z\tau_{0}+X_{j}\beta_{j}+\hat{\epsilon}$ (1)

Here $\hat{\epsilon}$ is a residual. In addiition, consider two related sample
regressions:

$\displaystyle Y$ $\displaystyle=Z\tau+\hat{e}$ (2) $\displaystyle X_{j}$
$\displaystyle=Z\Delta_{j}+\hat{u}$ (3)

Following Cochran’s formula [Cox (2007)], we may use (3) to rewrite (1) and
obtain a new representation for (2):

$\displaystyle\begin{split}Y{}&=X_{j}\beta_{j}+Z\tau_{0}+\hat{\epsilon}\\\
&=(Z\Delta_{j}+\hat{u})\beta_{j}+Z\tau_{0}+\hat{\epsilon}\\\
&=Z\Delta_{j}\beta_{j}+\hat{u}\beta_{j}+Z\tau_{0}+\hat{\epsilon}\\\
&=Z(\Delta_{j}\beta_{j}+\tau_{0})+(\hat{u}\beta_{j}+\hat{\epsilon})\end{split}$
(4)

Since the new error term $(\hat{u}\beta_{j}+\hat{\epsilon})$ is orthogonal to
$Z$ by the construction of residuals $\hat{u}$ and $\hat{\epsilon}$, we have
$\tau=(\Delta_{j}\beta_{j}+\tau_{0})$ and
$\hat{e}=\hat{u}\beta_{j}+\hat{\epsilon}$. Note that until now we have made no
model assumptions, merely fit regressions using sample quantities; however, if
we add a working assumption that model (1) is correct (i.e. that
$E(Y|X_{j},Z)=\beta^{pop}_{j}X_{j}+\tau^{pop}Z$ for some parameters
$\beta_{j}^{pop}$ and $\tau^{pop}$), then the difference

$\displaystyle\tau-\tau_{0}$ $\displaystyle=\Delta_{j}\beta_{j}$ (5)

is an estimate of the large-sample bias incurred by estimating treatment
effects via regression on $Z$ alone, ignoring $X_{j}$. We use the abbreviation
OVB to refer to quantity (5) in what follows.

Importantly for our purposes, the OVB can be rewritten in terms of Pearson
correlation between $X_{j}$ and $Y$ and a standardized mean difference with
normalization by the control standard deviation. The key is that when equation
(1) is fit on controls alone (as it will be in our pilot-sample approach),
both (1) and (3) are simple regressions. We now rewrite the corresponding
simple regression equations using familiar simple regression formulae.
$S_{Y_{pilot}}$ and $S_{X_{j,pilot}}$ denote the standard deviation of the
pilot sample for outcome and the standard deviation of the confounder in
question respectively. We include the pilot and analysis notations for
clarity.

$\displaystyle\beta_{j}$
$\displaystyle=r_{X_{j,pilot},Y_{pilot}}\frac{S_{Y_{pilot}}}{S_{X_{j,pilot}}}$
(6) $\displaystyle\Delta_{j}$
$\displaystyle=\bar{X}_{j1,analysis}-\bar{X}_{j0,analysis}$ (7)

Using (3), we obtain (7), where $\bar{X}_{j1,analysis}$ and
$\bar{X}_{j0,analysis}$ denote variable $j$’s sample means among treated
subjects and controls, respectively, in the analysis sample. Substituting into
expression (5) and rearranging, we obtain:

$\displaystyle\frac{\Delta_{j}\beta_{j}}{S_{Y_{pilot}}}$
$\displaystyle=r_{X_{j,pilot},Y_{pilot}}\frac{(\bar{X}_{j1,analysis}-\bar{X}_{j0,analysis})}{S_{X_{j,pilot}}}$
(8)

The left-hand side is a conveniently normalized version of the OVB that is
invariant to multiplicative rescalings of the outcome, and the right-hand side
is a product between a Pearson correlation computed in the pilot sample and a
standardized difference defined as follows:

$\displaystyle\text{OVB-based
SMD}=\frac{\bar{X}_{j1,analysis}-\bar{X}_{j0,analysis}}{S_{X_{j,pilot}}}$ (9)

Since the standardized outcome bias is a product of two terms, level sets for
bias take the form of hyperbolic curves on the jointVIP, which can be seen
plotted in Figures 2-3 to demarcate ”equivalent” levels of confounding under
the crude one-confounder model. In addition, an OVB number may be computed for
any individual variable using its respective SMD and outcome correlation, and
color-coding based on these quantities is used for plotting points in Figures
2-3.

Note that the definition for OVB-based SMD differs from the standard SMD
definition [Rosenbaum and Rubin (1985)], which uses pooled standard deviation
with equal weights for the treated and control group variances in the
denominator. However, calculating the SMD with a slightly different
standardization is not an uncommon approach. For example, Austin (2009) uses a
slightly different calculation of the pooled standard deviation for binary
variables. Hansen and Bowers (2008) use a pooled standard deviation with
weights that vary with the relative sizes of the treated and control groups.
Another variant, introduced by Stuart (2010) divides the numerator by the
standard deviation from only the treated group. This variant is similar to an
effect size estimator from Glass (1976), which standardizes the mean
difference by dividing the standard deviation from the control group [Hedges
(1981)].

A jointVIP implementation and code examples can be found in open-source
software R package jointVIP, currently hosted at GitHub. This package was used
to produce the figures given here.

##  4 Application: selecting a study design using jointVIP

###  4.1 The need for prioritization

Understanding how background variables contribute differentially to
confounding bias is especially important in settings where many variables are
present but not all can be perfectly controlled. In a regression model, this
might occur if the number of measured variables exceeds the number of
observations. Even when the problem is not high-dimensional in this particular
sense, adjusting for all background variables may not be practical. In our
study of pregnant individuals with GDM, we have 54 distinct background
variables, including missingness indicators. Even if all the variables were
binary, there would be $2^{54}\approx 1.8\times 10^{16}$ possible distinct
combinations, and it would be impossible to match all or even most treated
subjects to exactly identical controls despite our relatively large sample
($n$ = 29,969). Given that disagreement on some variables is inevitable, it
becomes vital to select the appropriate variables to match on versus the ones
to ignore or prioritize less highly. Many modern causal inference methods are
designed specifically with such prioritization in mind, offering the
researcher tuning parameters to incorporate substantive or data-driven
knowledge about which variables are likely to matter most [Sales, Hansen and
Rowan (2018); Athey, Imbens and Wager (2018); Wang and Zubizarreta (2020);
Ben-Michael, Feller and Rothstein (2021)]. The jointVIP can be very helpful in
selecting values for these tuning parameters. We illustrate this process in
the context of our study evaluating the impact of glyburide on GDM-complicated
pregnancies, in which we form matched pairs of treated and control individuals
using optimal matching with refined covariate balance, as implemented in R
package rcbalance [Pimentel (2016)].

###  4.2 Optimal matching with refined covariate balance: review

Despite the impossibility of matching each treated subject to an exactly
identical control, it may be possible to create matched pairs in such a way
that many individual background variables are balanced overall among the
treated subjects and the control subjects selected for the match. For example,
suppose that we form a matched pair in which both subjects have high BMI but
the treated subject is 23 years old while the control is 31 years old. While
this pair alone is poorly matched on age, in another pair we could match a 31
years old treated subject with moderate BMI to a 23 years old control with
moderate BMI; taking the two pairs together and comparing the _group_ of
matched controls to the _group_ of matched treated subjects, we have identical
marginal distributions for both age and BMI. In contrast to traditional
optimal matching, which focuses entirely on similarity within individual
pairs, optimal matching with balance constraints restricts attention to
matched-pair configurations that achieve this form of marginal balance,
sometimes called fine balance, on one or more specified variables [Rosenbaum,
Ross and Silber (2007); Yang et al. (2012)].

Refined covariate balance adapts the idea of marginal balance to datasets with
many different background variables of differing levels of importance by
allowing researchers to specify several tiers of nominal variables to be
balanced, given in descending order of priority [Pimentel et al. (2015)]. The
top-tier variables are balanced as closely as possible (as though they were
the only variables present), the second-tier variables are then balanced as
closely as remains possible, and so on, with optimally-close pairs selected
only among configurations that ensure all these constraints are met. More
precisely, to exhibit refined covariate balance given a particular sequence of
balance tiers (or sets of variables) $\mathcal{X}_{1},\ldots,\mathcal{X}_{L}$,
a candidate match $m$ (consisting of a specific subset of the overall treated
group and a corresponding matched subset of the overall control group) chosen
from a family of allowed matches $\mathcal{M}$ must satisfy the following
conditions, where $TV(m,\mathcal{S})$ is the total variation distance between
the empirical joint distributions of all background variables in set
$\mathcal{S}$ in the matched treated group and in the matched control group
specified by $m$:

  1. 1.

Joint imbalance on the variables in set $\mathcal{X}_{1}$ is minimized, i.e.

$TV(m,\mathcal{X}_{1})=\min_{m^{\prime}\in\mathcal{M}}TV(m^{\prime},\mathcal{X}_{1})$

  2. 2.

Joint imbalance on the variables in set $\mathcal{X}_{1}\cup\mathcal{X}_{2}$
is minimized subject to 1., i.e.

$TV(m,\mathcal{X}_{1}\cup\mathcal{X}_{2})=\min_{m^{\prime}\in\mathcal{M},m^{\prime}\text{
satisfies 1.}}TV(m^{\prime},\mathcal{X}_{1}\cup\mathcal{X}_{2})$

$\vdots$

  3. L.

Joint imbalance on the variables in set
$\mathcal{X}_{1}\cup\ldots\cup\mathcal{X}_{L}$ is minimized subject to 1., 2.,
…, L-1, i.e.

$TV(m,\mathcal{X}_{1}\cup\ldots\cup\mathcal{X}_{L})=\min_{m^{\prime}\in\mathcal{M},m^{\prime}\text{
satisfies 1., 2., \ldots,
L-1.}}TV(m^{\prime},\mathcal{X}_{1}\cup\ldots\cup\mathcal{X}_{L})$

For a full technical characterization and an efficient algorithm for computing
matches under such constraints see Pimentel et al. (2015). At a high level,
this constraint allows top-priority variables and their interactions to be
balanced as though they were the only variables in the study, with lower-
priority variables receiving further attention as possible. Typically a
covariate distance such as the robust Mahalanobis distance, defined using all
background variables is minimized under the refined covariate balance so that
even variables not included in any tier receive some attention. See Rosenbaum
(2010), Chapter 8.3, for details on robust Mahalanobis distance.

Refined balance may be combined with other restrictions on matching, notably
exact matching, in which subjects in any pair are required to share identical
values for some chosen nominal covariate, and caliper matching, in which
subjects in any pair are required to differ by no more than some prespecified
amount (the caliper) on some chosen continuous variable [Rosenbaum (2020)].
Exact matching and caliper matching (with a sufficiently narrow caliper) are
more stringent requirements than requiring marginal balance on a given
variable, in the sense that both produce perfect or near-perfect marginal
balance, while enforcing marginal balance need to not produce similar matches
in any given pair. When it is not possible to match all treated subjects to
distinct controls while respecting the exact match constraint or the caliper,
the largest possible match that does respect those constraints is formed
optimally, and the remaining treated subjects are excluded [Rosenbaum (2012);
Pimentel and Kelz (2020)].

While refined covariate balance offers substantial flexibility to the
researcher, it also relies on strong substantive knowledge to specify the
balance tiers in a reasonable manner. Frequently it is not immediately clear
how to organize a group of background variables into balance tiers in a
principled way. The jointVIP offers a data-driven approach in settings where
ambiguity remains even after accounting for substantive knowledge.

Figure 2: Pre-post match results for Cesarean section delivery

BMI: body mass index, C-section: Cesarean section, GDM: gestational diabetes,
OGTT: oral glucose tolerance test, SMD: standardized mean difference

Figure 3: Pre-post match results for neonatal intensive care unit admission

GDM: gestational diabetes, jointVIP: joint [treatment-outcome] variable
importance plot, NICU: neonatal intensive care unit, OGTT: oral glucose
tolerance test, SMD: standardized mean difference

###  4.3 The jointVIP as a guide for balance constraint selection

In what follows we use jointVIP, in concert with clinical knowledge, to select
balance tiers for refined covariate balance and calipers for optimal matching.
Since we have two different outcomes, C-section delivery and NICU admission,
we create separate jointVIPs and conduct separate matches based on different
balance constraints for each outcome.

Firstly, the variable with the largest potential bias in both plots is the
propensity score (calculated by LASSO regression). To ensure particularly
stringent control of the propensity score, we impose a caliper equal to 0.2
standard deviations of the fitted propensity score values in the entire
sample. Using a caliper on the propensity score is also a natural choice
because our approach to inference relies on similar propensity scores within
matched pairs [Pimentel (2022)], as discussed in Section 5.1. Secondly, we
note that one variable, year, is omitted entirely from the jointVIP due to a
perfect lack of overlap between the pilot sample (collected 2007-2010) and the
analysis sample (collected 2011-2021). To address substantive concerns about
potential for temporal shifts in the standard of care in the absence of
reliable outcome correlations for this variable, we match exactly on year for
both outcomes.

We address potential bias from additional variables by imposing a series of
refined balance constraints tailored to each outcome. Since the prognostic
score (fit in the pilot sample using LASSO regression) ranks among the
variables contributing the largest potential bias for both outcomes, we
include quintiles of the prognostic score in balance tiers for both matched
designs (first tier for C-section delivery and second tier for NICU
admission). We go on to include all variables contributing potential bias
(rounded to the third decimal place) greater than or equal to 0.010 in some
balance tier for each design, with variables that contributed larger amounts
of bias in higher tiers. Table 2 summarizes the chosen balance tiers for each
design. Specific potential bias values can be found in Appendix Table 8 column
Pre-matched bias. In addition, gestational age at GDM diagnosis is a
continuous variable and had to be discretized for compatibility with refined
covariate balance. Three categories were used — before 24 weeks, between 24 to
32 weeks, and over 32 weeks — with cutoffs chosen according to natural
thresholds in the empirical distribution. These cutoffs differ slightly from
standard clinical thresholds; for more discussion of the choice of categories,
see Appendix Figures 5.A-5.B.

Balance tier | C-section delivery | NICU admission
---|---|---
1 | Prognostic score quintile | GDM diagnosed by Carpenter-Coustan criteria
2 | OGTT for fasting blood glucose | OGTT for fasting blood glucose
| Obese pre-pregnancy BMI | Gestational age category at GDM diagnosis
|  | Prognostic score quintile
3 | GDM diagnosed by Carpenter-Coustan criteria |
| Gestational age category at GDM diagnosis |
| Pre-pregnancy pre-diabetes |
| Normal pre-pregnancy BMI |
Table 2: Balance tiers for refined covariate balance for each outcome, chosen
using jointVIP plots.

Matches were formed for each design using rcbalance package in R. Post-matched
jointVIP results, reflecting new levels of balance after matching, are plotted
in Figure 2.B and Figure 3.B. For variables that were specified, post-matched
biases are compared to pre-matched biases in Appendix Table 8, which shows all
background variables and summary measures to have small biases (around 0.005
or less) post-matching. Note in particular that variables with high outcome
correlation are balanced as well, a feature of the design that traditional
methods based on Love plots are not equipped to guarantee.

For both matched designs, 2,093 treated subjects are excluded due to caliper
and exact matching constraints, and 8,693 pairs are matched; although slightly
different sets of treated subjects are selected for exclusion in the two
matches, the two groups of excluded treated subjects have similar summary
statistics (see Appendix Table LABEL:excluded). The excluded treated people
differ from the included treated in several ways, for example, there are more
excluded treated individuals who have a pre-pregnancy BMI in the obese
category and earlier GDM diagnoses (Appendix Table LABEL:excluded). Those who
are excluded may have signs of severe GDM and have a high probability of
treatment without comparable controls for matching. As such, we note that the
average risk differences for C-section and NICU admission are best understood
not as estimates of an average treatment effect on the treated [Stuart
(2010)], but as an average effect on a ”marginal” population consisting of
individuals for whom treatment by either arm is reasonably likely [Rosenbaum
(2012); Li and Greene (2013); Li, Thomas and Li (2019); Greifer and Stuart
(2021)]. This estimand, while less common in theoretical discussions of causal
inference, adheres more closely to the substantive quantity of interest for
physicians who are typically more interested in guidance for patients with
equipoise, and less interested in effects on patients who would clearly be
assigned either MNT or glyburide in the large majority of cases.

In summary, using jointVIP and clinical expertise to select reasonable tuning
parameters produce well-balanced groups that accounted for each variable’s
relationship against the outcome and its treatment-control imbalance. These
large, well-balanced samples mimic the kinds of datasets that might be
produced in RCTs assessing the effect of glyburide versus MNT on risk of
adverse perinatal outcomes and provide a strong basis for interrogating causal
effects of glyburide among our study subjects with GDM.

##  5 Does glyburide have adverse impacts on perinatal outcome?

###  5.1 Inference and sensitivity analysis for matched pairs: review

We index matched pair by $i=1,\cdots,I$, and individuals in the matched-pair
by $k=1,2$. For a matched pair $i$, one person is treated with glyburide,
$Z_{ik}=1$, and the other with MNT, $Z_{ik}=0$, hence we have
$Z_{i1}+Z_{i2}=1$. Let $\mathcal{Z}$ denote the event that $Z_{i1}+Z_{i2}=1$
for each matched pair $i$. Each subject $ik$ has corresponding potential
outcomes $Y_{{ik}}(1)$ and $Y_{{ik}}(0)$ for treatment with and without
glyburide respectively. We write
$\mathcal{F}=\\{(Y_{{ik}}(1),Y_{{ik}}(0),\boldsymbol{x}_{ik}),i=1,\cdots,I,k=1,2\\}$.
For our study, both outcomes of interest are binary.

We test the Fisher (1936) sharp null hypothesis,
$H_{0}:Y_{{ik}}(1)=Y_{{ik}}(0)$ for all $i,k$, for C-section delivery and
again for NICU admission. Assuming that paired subjects are equally likely to
receive glyburide or MNT, we can test this hypothesis by repeatedly permuting
treatment indicators within pairs (independently across pairs) with
probability $\frac{1}{2}$. Since under the sharp null the outcomes remain
identical regardless of treatment assignment, we can compute a test statistic
under each permutation using observed outcomes and compare the actual observed
value of the test statistic to this reference distribution to conduct
inference. For binary outcomes, in particular, we may apply McNemar’s test
[McNemar (1947)]. Note that while such tests are sometimes described as
randomization tests due to their origins in the analysis of randomized trials,
in an observational study like this one where no randomization is present they
are better understood as quasi-randomization tests in the sense of Zhang and
Zhao (2022).

Note that the above procedure relies on the assumption
$Pr(Z_{ik}=1|\mathcal{F},\mathcal{Z})=\frac{1}{2}$ with independent assignment
for each pair, which is true when unobserved confounding is absent and
propensity scores are matched exactly. In real observational studies this
assumption may fail, and sensitivity analysis is needed to probe the
robustness of the initial findings to such failures [Rosenbaum (2010) §3]. A
convenient model for sensitivity analysis allows probabilities to differ
within matched pairs but only by a fixed amount:

$\displaystyle\frac{1}{1+\Gamma}\leq
Pr(Z_{ik}=1|\mathcal{F},\mathcal{Z})\leq\frac{\Gamma}{1+\Gamma}$ (10)

where the sensitivity parameter $\Gamma\geq 1$ controls the degree or strength
of residual confounding. When model (10) holds for a given $\Gamma$, we no
longer know the exact permutation distribution that should be used to test
$H_{0}$, but under some mild conditions on the test statistic (satisfied by
McNemar’s test) it is possible to conduct ”worst-case” inference, producing an
upper bound on the p-value associated with the true distribution by maximizing
over distributions allowed by the model. In practice we may not know the true
$\Gamma$, but for an initially significant result one can conduct repeated
sensitivity analyses for successively larger $\Gamma$-values until the
sensitivity analysis no longer rejects the null, reporting the ”threshold”
value of $\Gamma$ needed to overturn the causal conclusion.

###  5.2 Results

Table 3: Matched analysis for Cesarean-section delivery

| Treated with glyburide
---|---
|  | C-section | not C-section
Control (MNT only) | C-section | $1078$ | $1844$
| not C-section | $2023$ | $3748$

C-section: Cesarean-section, MNT: medical nutritional therapy

For the C-section outcome, of the $2\times 8,693$ individuals who are matched
in pairs, 6,023 (34.64%) individuals delivered by C-section. Matched results
are shown in Table 3. For control (MNT only) individuals, 33.61% delivered by
C-section, and for treated (glyburide and MNT) individuals, 35.67% delivered
by C-section (raw treatment-control difference of 2.06%). McNemar’s test
yields a one-sided p-value of 0.0020. Evaluating at significance level 0.05,
there is evidence to reject the null hypothesis under no unmeasured
confounding assumption. We find the threshold $\Gamma$ to be 1.041, which
indicates that a relatively small degree unmeasured confounding can explain
away the causal effect detected.

Table 4: Matched analysis for neonatal intensive care unit admission

| Treated with glyburide
---|---
|  | NICU | not NICU
Control (MNT only) | NICU | $163$ | $996$
| not NICU | $1004$ | $6530$

MNT: medical nutritional therapy, NICU: neonatal intensive care unit

For the NICU outcome, of the $2\times 8,693$ individuals, 2,326 (13.38%)
individuals’ infants were admitted to the NICU. Matched results are shown in
Table 4. For control individuals, 13.33% of infants were admitted, and for
treated individuals, 13.42% of whose infants were admitted (raw treatment
control difference 0.09%). McNemar’s test produces a one-sided p-value of
0.429, so the sharp null hypothesis is plausible even without accounting for
unobserved confounding.

##  6 Discussion

To understand whether glyburide increases the risks of perinatal complications
among individuals with GDM, we examined its impact on Cesarean section
(C-section) delivery as a maternal outcome of interest and neonatal intensive
care unit (NICU) admission as a neonatal outcome of interest. We introduced a
novel tool, the joint treatment-outcome variable importance plot (jointVIP),
to help identify variables for prioritization at the design stage. We
demonstrated the utility of jointVIP to inform our match designs in a
population-based pregnancy cohort with GDM, who delivered between 2011 and
2021. We found no significant matched risk difference using randomization
inference in NICU admission between treated (glyburide and MNT) and control
(MNT only) pairs. Although there is a slight matched risk difference in the
C-section delivery, the significant effect can be entirely explained away by a
mild unmeasured confounder $\Gamma=1.041$. For example, under the patient-
physician shared joint decision model, patient-physician pairs with a high
level of confidence in medications may be more inclined both to choose
glyburide and more likely to receive C-section, even in the absence of a
causal link. As such we find no substantial evidence that glyburide is causing
either of the adverse perinatal outcomes in this study.

In the absence of large-scale RCTs on this topic, observational studies are
essential to understand the safety of glyburide and provide better care for
individuals with GDM. Our work helps fill this important gap by providing
evidence to inform policies around glyburide use. We did not find strong
evidence for glyburide causing an increased incidence of C-section deliveries
and/or NICU admissions amongst almost 8.7 thousand matched pairs based on our
study design. However, further evidence indicating the safety of glyburide in
the GDM-complicated population is needed before official recommendations can
be made. Specifically, other medication safety-related outcomes, for example,
placental abruption, small for gestational age, and preterm birth, need to be
evaluated carefully.

JointVIP is also useful as a tool for applied statisticians to select
variables to balance during the observational study design phase. One notable
advantage of theover traditional methods is the visual ease of comparison
between the marginal relationships for each variable against both the outcome
and treatment. Since the outcome relationship is considered, methods
leveraging jointVIP offer better bias reduction and the potential for
increased robustness against unmeasured confounders [Rosenbaum (2005)].

The uses of jointVIP extend beyond the specific cases considered in our study.
In our application we illustrated how jointVIP can guide the choice of
parameters for matching with refined covariate balance; however, jointVIP can
help provide a natural choice for tuning parameters in many other existing
methods. For example, inverse values of the outcome correlations plotted on
the y-axis of the jointVIP could be used as balance tolerances in stable
balancing weights and cardinality matching [Zubizarreta (2015); Zubizarreta,
Paredes and Rosenbaum (2014)] to ensure that variables important in the
outcome model are balanced more stringently.

Another important potential use of jointVIP is to identify variables with high
treatment model importance but negligible outcome model importance. A
substantial body of the literature cautions against focusing adjustments on
such variables, which are sometimes referred to as instrumental variables or
prods. In particular, when unobserved bias is present adjusting for these
variables can increase bias [Brooks and Ohsfeldt (2013); Pimentel, Small and
Rosenbaum (2016); Ding, VanderWeele and Robins (2017)]. In studies with low
degrees of confounding in observed variables, propensity scores themselves may
take on a similar pathological role [King and Nielsen (2019)]. JointVIP
highlights any such variables clearly, placing them in the upper reaches of
the x-axis but near zero on the y-axis. In addition, using jointVIP to guide
design provides a natural opportunity either to exclude such variables from
matching or weighting techniques, or (in cases where it is not entirely clear
whether a variable should be excluded) to construct multiple control groups
that adjust for these variables differently as detailed in Pimentel, Small and
Rosenbaum (2016).

Several other authors have discussed ideas closely related to the joint
variable-importance plot. Meng (2018) presents a formula for bias due to
selective sampling that is very similar to the formula we use to compute
estimates of omitted variable bias, involving a product between an outcome-
selection covariance and a measure of the degree of selectivity. Zhao and Yang
(2022) propose variable selection for fitting generalized propensity scores
using measures of outcome importance and provide supporting theory suggesting
the optimality of this approach. Aikens, Greaves and Baiocchi (2020) and
Aikens and Baiocchi (2022) construct an alternative design-stage visualization
based partially on a pilot sample incorporating outcomes, the assignment-
control (AC) plot. In contrast to jointVIP however, the AC plot represents
subjects rather than variables on the plot, using the estimated prognostic
score and propensity score values on the axes. AC plots and jointVIP thus
provide valuable complementary representations of observational data at the
design stage, and we recommend that investigators use both.

Although jointVIP can inform design decisions for observational studies, it
should not be taken as a substitute for consulting scientific collaborators
about substantive priorities. Clinical expertise is essential in the design
and interpretation of observational studies using data from health records
[Rosenbaum and Silber (2001)]. In our study, clinical collaborators helped
identify several variables for which the definition in the electronic health
database had shifted subtly between pilot sample collection and analysis
sample collection, and these were then excluded or recoded appropriately to
ensure congruence between the variables in the two samples. More generally,
when there is a lack of overlap in background variable values between the
pilot sample and the analysis sample, as with our year variable or as may
occur with very rare variables or categories, jointVIP will be unable to
provide guidance and researchers must rely entirely on substantive
considerations.

One area for further practical consideration and development is the
appropriate choice of the pilot sample. Ideally, a pilot sample is drawn from
the same joint variable-outcome population distribution that produced the
analysis sample and is large enough to allow precise measurement of outcome
correlations without restricting the size of the analysis sample. In practice,
there are tradeoffs among these goals. In our study, we used a historical
sample predating the collection of the analysis sample, which ensured large
sample sizes for both pilot and analysis samples; however, we must rely on an
assumption that variable-outcome relationships do not change dramatically
between the collection of the two samples. We do observe some substantial
differences in the variable distributions themselves, most notably in the
patterns of missingness between the two samples (shown in Appendix Table 6).
To allay these concerns, one could instead extract a subset of the analysis
sample to use as the pilot sample as described in Section 3.1 at the cost of
reduced sample size for both samples. In ongoing work, we aim to study
additional outcomes in the Kaiser Permanente GDM cohort using this strategy
for jointVIP construction.

Another area for future study is the coordination of a common subpopulation or
match across multiple outcomes. In this study we conducted two different
matches for two different outcomes, producing inferences over two distinct
subsets of the original study population. Fortuitously the two subpopulations
appeared very similar in aggregate, and our goal of testing a sharp null
hypothesis of zero effect has a clear interpretation even when conducted in
slightly distinct populations. However, it is not guaranteed that different
matching specifications driven by jointVIPs for distinct outcomes will select
similar subpopulations, and in estimation-focused studies, it may be difficult
to explain a set of outcome estimates for estimands defined over such distinct
subpopulations. An alternate approach might be to integrate information across
multiple jointVIPs for multiple outcomes, either into a single matching or
weighting specification or by conducting several matches or weighting
procedures that all target a common template distribution as described in
Bennett, Vielma and Zubizarreta (2020).

##

Table 5: Summary of all baseline variables of pregnant individuals with gestational diabetes. |  | 2007-2010 | 2011-2021 | 2011-2021
---|---|---|---|---
|  | Control | Control | Treated
|  | n = 7,526 | n = 19,183 | n = 10,786
Age (%) | Under 25 | 552 (7.3) | 1,029 (5.4) | 349 (3.2)
| Between 25-29 | 1,665 (22.1) | 3,762 (19.6) | 1,866 (17.3)
| Between 30-34 | 2,584 (34.3) | 7,101 (37.0) | 4,165 (38.6)
| Over 35 | 2,725 (36.2) | 7,291 (38.0) | 4,406 (40.8)
KP member 6 months prior to pregnancy = yes (%) |  | 5,891 (78.3) | 15,752 (82.1) | 9,004 (83.5)
Median housing income (%) | Less than $ 40,000 | 892 (11.9) | 784 (4.1) | 355 (3.3)
| $ 40,000 - $ 59,999 | 1,800 (23.9) | 2,807 (14.6) | 1,325 (12.3)
| $ 60,000 - $ 79,999 | 1,903 (25.3) | 3,769 (19.6) | 2,132 (19.8)
| $ 80,000 and above | 2,931 (38.9) | 11,823 (61.6) | 6,974 (64.7)
Parity (%) | 0 | 3,121 (41.5) | 7,611 (39.7) | 3,873 (35.9)
| 1 | 2,347 (31.2) | 6,566 (34.2) | 3,994 (37.0)
| 2 or more | 2,058 (27.3) | 5,006 (26.1) | 2,919 (27.1)
Pre-pregnancy BMI (%) | Underweight | 100 (1.3) | 391 (2.0) | 76 (0.7)
| Normal | 1,921 (25.5) | 5,107 (26.6) | 1,706 (15.8)
| Overweight | 2,847 (37.8) | 6,369 (33.2) | 3,361 (31.2)
| Obese | 2,658 (35.3) | 7,316 (38.1) | 5,643 (52.3)
Race/ethnicity (%) |  Asian or Pacific Islander | 2,919 (38.8) | 8,553 (44.6) | 4,560 (42.3)
| Hispanic (%) | 2,322 (30.9) | 5,013 (26.1) | 2,923 (27.1)
| White (%) | 1,602 (21.3) | 4,090 (21.3) | 2,381 (22.1)
|  Black or African American | 315 (4.2) | 736 (3.8) | 410 (3.8)
| Other/Unknown (%) | 368 (4.9) | 791 (4.1) | 512 (4.7)
Singleton pregnancy = yes (%) |  | 7,286 (96.8) | 18,579 (96.9) | 10,567 (98.0)
Alcohol consumption prior to pregnancy = yes (%) |  | 2,915 (38.7) | 7,715 (40.2) | 4,954 (45.9)
Alcohol consumption during pregnancy = yes (%) |  | 437 (5.8) | 1,893 (9.9) | 946 (8.8)
Smoking prior to pregnancy = yes (%) |  | 611 (8.1) | 919 (4.8) | 690 (6.4)
Smoking during pregnancy = yes (%) |  | 207 (2.8) | 314 (1.6) | 222 (2.1)
Arrythmia diagnosis = yes (%) |  | 52 (0.7) | 194 (1.0) | 107 (1.0)
Asthma diagnosis = yes (%) |  | 737 (9.8) | 2,279 (11.9) | 1,507 (14.0)
Chronic hypertension = yes (%) |  | 435 (5.8) | 846 (4.4) | 591 (5.5)
Depression diagnosis prior to pregnancy = yes (%) |  | 624 (8.3) | 1,538 (8.0) | 1,020 (9.5)
Depression diagnosis during pregnancy = yes (%) |  | 455 (6.0) | 2,639 (13.8) | 1,626 (15.1)
Dyslipidemia diagnosis = yes (%) |  | 305 (4.1) | 1173 (6.1) | 861 (8.0)
Family history of diabetes = yes (%) |  | 198 (2.6) | 1202 (6.3) | 822 (7.6)
History of abortive outcome = yes (%) |  | 694 (9.2) | 2,253 (11.7) | 1,412 (13.1)
History of GDM = yes (%) |  | 856 (11.4) | 3,401 (17.7) | 2,608 (24.2)
History of macrosomia = yes (%) |  | 69 (0.9) | 145 (0.8) | 147 (1.4)
Polycystic ovary syndrome by diagnosis = yes (%) |  | 263 (3.5) | 831 (4.3) | 653 (6.1)
Pre-pregnancy pre-diabetes = yes (%) |  | 479 (6.4) | 1,802 (9.4) | 1,915 (17.8)
Count of blood pressure measurements taken (mean (SD)) |  | 6.21 (9.35) | 4.65 (3.49) | 5.11 (3.83)
Diastolic blood pressure prior to pregnancy (mean (SD)) |  | 72.89 (8.32) | 71.44 (10.15) | 72.47 (10.08)
Systolic blood pressure prior to pregnancy (mean (SD)) |  | 116.69 (11.53) | 117.26 (13.38) | 118.86 (13.36)
Diastolic blood pressure prior to GDM diagnosis (mean (SD)) |  | 68.63 (8.76) | 66.82 (9.70) | 67.83 (9.59)
Systolic blood pressure prior to GDM diagnosis (mean (SD)) |  | 114.45 (12.17) | 114.43 (13.14) | 116.37 (12.94)
Average diastolic blood pressure prior to pregnancy (mean (SD)) |  | 72.96 (7.38) | 71.26 (8.05) | 72.24 (8.06)
Average systolic blood pressure prior to pregnancy (mean (SD)) |  | 117.33 (10.42) | 117.14 (11.00) | 118.70 (10.98)
Median diastolic blood pressure prior to pregnancy (mean (SD)) |  | 72.98 (7.52) | 71.22 (8.26) | 72.21 (8.29)
Median systolic blood pressure prior to pregnancy (mean (SD)) |  | 117.16 (10.56) | 116.86 (11.23) | 118.41 (11.25)
Average diastolic blood pressure prior to GDM (mean (SD)) |  | 69.26 (7.37) | 67.68 (8.08) | 68.89 (8.18)
Average systolic blood pressure prior to GDM (mean (SD)) |  | 115.11 (10.55) | 115.38 (11.15) | 117.40 (11.17)
Median diastolic blood pressure prior to GDM (mean (SD)) |  | 69.20 (7.55) | 67.56 (8.29) | 68.80 (8.38)
Median systolic blood pressure prior to GDM (mean (SD)) |  | 114.91 (10.73) | 115.13 (11.38) | 117.16 (11.40)
Infant sex (%) | Female | 3,577 (47.5) | 9,250 (48.2) | 5,145 (47.7)
| Male | 3,855 (51.2) | 9,862 (51.4) | 5,605 (52.0)
| Unknown | 94 (1.2) | 71 (0.4) | 36 (0.3)
Glucose challenge test value (mean (SD)) |  | 169.43 (22.38) | 169.71 (22.14) | 173.22 (24.32)
Gestational age at GDM diagnosis (mean (SD)) |  | 26.06 (6.10) | 26.86 (6.02) | 23.53 (7.12)
Gestational weight gain up to GDM diagnosis (mean (SD)) |  | 15.28 (11.04) | 14.82 (10.99) | 12.86 (11.35)
Gestational hypertension = yes (%) |  | 385 (5.1) | 1,546 (8.1) | 884 (8.2)
GDM diagnosed by the C-C criteria = yes (%) |  | 7,323 (97.3) | 17,303 (90.2) | 8,419 (78.1)
One-hour OGTT = abnormal (%) |  | 6,536 (86.8) | 15,534 (81.0) | 8,244 (76.4)
Two-hour OGTT = abnormal (%) |  | 6,662 (88.5) | 15,829 (82.5) | 7,082 (65.7)
Three-hour OGTT = abnormal (%) |  | 2,875 (38.2) | 7,194 (37.5) | 2,966 (27.5)
OGTT for fasting blood glucose = abnormal (%) |  | 2,165 (28.8) | 3,762 (19.6) | 4,512 (41.8)
GDM severity = severe (%) |  | 775 (10.3) | 1,744 (9.1) | 1,215 (11.3)
BMI: body mass index, C-C: Carpenter-Coustan, GDM: gestational diabetes, KP:
Kaiser Permanente,
OGTT: oral glucose tolerance test, SD: standard deviation.
A normal OGTT fasting blood glucose level is lower than 95 mg/dL.
A normal one-hour OGTT blood glucose level is lower than 180 mg/dL.
A normal two-hour OGTT blood glucose level is lower than 155 mg/dL.
A normal three-hour OGTT blood glucose level is lower than 140 mg/dL.

Table 6: Missingness summary and out-of-bag imputation error estimate.

Variables with missingness indicators
---
| 2007-2010 | 2011-2021 | 2011-2021
| Control | Control | Treated
| n = 7,526 | n = 19,183 | n = 10,786
Gestational weight gain up to GDM diagnosis (%) | 1,391 (18.5) | 2,504 (13.1) | 1,307 (12.1)
Blood pressure measured prior to pregnancy (%) | 2,988 (39.7) | 2,914 (15.2) | 1,353 (12.5)
Blood pressure measured prior to GDM (%) | 658 (8.7) | 35 (0.2) | 46 (0.4)
Pre-pregnancy BMI (%) | 1,326 (17.6) | 866 (4.5) | 350 (3.2)
Glucose challenge test value (%) | 204 (2.7) | 152 (0.8) | 104 (1.0)
Median housing income | 4 (0.1) | 7 (0.0) | 0 (0.0)
Parity1 | 0 (0.0) | 26 (0.0) | 13 (0.0)
Out-of-bag imputation error
Year NRMSE PFC
2007 2.14*10-6 7.96*10-2
2008 2.27*10-6 0.00
2009 2.35*10-6 8.02*10-2
2010 2.35*10-6 7.83*10-2
2011 3.83*10-5 1.01*10-1
2012 4.17*10-5 0.00
2013 4.09*10-5 9.73*10-2
2014 4.08*10-5 0.00
2015 4.13*10-5 0.00
2016 4.13*10-5 1.69*10-1
2017 4.04*10-5 8.37*10-2
2018 3.98*10-5 1.76*10-1
2019 3.79*10-5 0.00
2020 3.43*10-5 1.82*10-1
2021 3.28*10-5 1.82*10-1

BMI: body mass index, GDM: gestational diabetes, NRMSE: Root mean squared
error, PFC: proportion of falsely classified entries
1Parity is only missing in analysis 2011-2021 dataset. Since the missingness
is quite small compared to data available, this indicator is dropped after
imputation.
2007-2010 is the pilot data, and 2011-2021 is the analysis data.
The out-of-bag error imputation is calculated separately for continuous
(NRMSE) and categorical variables (PFC).

Figure 4: Comparison between the Love plot and the joint outcome-treatment
variable importance plot with raw measures

BMI: body mass index, C-C: Carpenter-Coustan, C-section: Cesarean section,
GDM: gestational diabetes, OGTT: oral glucose tolerance test, OVB: omitted
variable bias, PI: Pacific Islander, SD: standard deviation, SMD: standardized
mean difference
The Love plot (A) is a common practice in visualizing treatment imbalance
using standard SMD. The joint variable importance plot (B) is our proposed
method to visualize both treatment imbalance and outcome correlation
marginally building on the OVB framework using a slightly different OVB-based
SMD.

Figure 5: Categorization of gestational age at GDM diagnosis

GDM: gestational diabetes
The categorization of gestational age at GDM diagnosis is used only for the
refined covariate balance matching method. The choice is based on the
empirical distribution (Fig. 5.A) over the typical clinical diagnosis time
(Fig. 5.B). The majority of diagnosis (over 70%) was performed in the 24-32
weeks gestational age time span.

Table 7: Summary of all baseline variables of post-matched treated pregnant individuals. |  | C-section delivery | C-section delivery | NICU admission
---|---|---|---|---
|  | included | excluded | excluded
|  | n = 8,693 | n = 2,093 | n = 2,093
Age (%) | Under 25 | 322 (3.7) | 27 (1.3) | 33 (1.6)
| Between 25-29 | 1,522 (17.5) | 344 (16.4) | 353 (16.9)
| Between 30-34 | 3,405 (39.2) | 760 (36.3) | 789 (37.7)
| Over 35 | 3,444 (39.6) | 962 (46.0) | 918 (43.9)
KP member 6 months prior to pregnancy = yes (%) |  | 7,229 (83.2) | 1,775 (84.8) | 1,762 (84.2)
Median housing income (%) | Less than $ 40,000 | 286 (3.3) | 69 (3.3) | 71 (3.4)
| $ 40,000 - $ 59,999 | 1,051 (12.1) | 274 (13.1) | 264 (12.6)
| $ 60,000 - $ 79,999 | 1,720 (19.8) | 412 (19.7) | 415 (19.8)
| $ 80,000 and above | 5,636 (64.8) | 1,338 (63.9) | 1,343 (64.2)
Parity (%) | 0 | 3,214 (37.0) | 659 (31.5) | 676 (32.3)
| 1 | 3,203 (36.8) | 791 (37.8) | 769 (36.7)
| 2 or more | 2,276 (26.2) | 643 (30.7) | 648 (31.0)
Pre-pregnancy BMI (%) | Underweight | 69 (0.8) | 7 (0.3) | 5 (0.2)
| Normal | 1,612 (18.5) | 94 (4.5) | 133 (6.4)
| Overweight | 2,862 (32.9) | 499 (23.8) | 538 (25.7)
| Obese | 4,150 (47.7) | 1,493 (71.3) | 1,417 (67.7)
Race/ethnicity (%) |  Asian or Pacific Islander | 3,796 (43.7) | 764 (36.5) | 763 (36.5)
| Hispanic (%) | 2,316 (26.6) | 607 (29.0) | 618 (29.5)
| White (%) | 1,880 (21.6) ) | 501 (23.9) | 502 (24.0)
|  Black or African American | 315 (3.6) | 95 (4.5) | 89 (4.3)
| Other/Unknown (%) | 386 (4.4) | 126 (6.0) | 121 (5.8)
Singleton pregnancy = yes (%) |  | 8,508 (97.9) | 2,059 (98.4) | 2,062 (98.5)
Alcohol consumption prior to pregnancy = yes (%) |  | 3,731 (42.9) | 1,223 (58.4) | 1,207 (57.7)
Alcohol consumption during pregnancy = yes (%) |  | 750 (8.6) | 196 (9.4) | 184 (8.8)
Smoking prior to pregnancy = yes (%) |  | 467 (5.4) | 223 (10.7) | 216 (10.3)
Smoking during pregnancy = yes (%) |  | 156 (1.8) | 66 (3.2) | 60 (2.9)
Arrythmia diagnosis = yes (%) |  | 81 (0.9) | 26 (1.2) | 27 (1.3)
Asthma diagnosis = yes (%) |  | 1,138 (13.1) | 369 (17.6) | 375 (17.9)
Chronic hypertension = yes (%) |  | 397 (4.6) | 194 (9.3) | 177 (8.5)
Depression diagnosis prior to pregnancy = yes (%) |  | 786 (9.0) | 234 (11.2) | 237 (11.3)
Depression diagnosis during pregnancy = yes (%) |  | 1,230 (14.1) | 396 (18.9) | 383 (18.3)
Dyslipidemia diagnosis = yes (%) |  | 619 (7.1) | 242 (11.6) | 243 (11.6)
Family history of diabetes = yes (%) |  | 590 (6.8) | 232 (11.1) | 231 (11.0)
History of abortive outcome = yes (%) |  | 1,102 (12.7) | 310 (14.8) | 314 (15.0)
History of GDM = yes (%) |  | 1,930 (22.2) | 678 (32.4) | 674 (32.2)
History of macrosomia = yes (%) |  | 105 (1.2) | 42 (2.0) | 47 (2.2)
Polycystic ovary syndrome by diagnosis = yes (%) |  | 468 (5.4) | 185 (8.8) | 176 (8.4)
Pre-pregnancy pre-diabetes = yes (%) |  | 1,213 (14.0) | 702 (33.5) | 677 (32.3)
Count of blood pressure measurements taken (mean (SD)) |  | 4.95 (3.58) | 5.77 (4.67) | 5.68 (4.48)
Diastolic blood pressure prior to pregnancy (mean (SD)) |  | 72.05 (9.99) | 74.21 (10.27) | 73.99 (10.36)
Systolic blood pressure prior to pregnancy (mean (SD)) |  | 118.10 (13.26) | 121.97 (13.30) | 121.70 (13.43)
Diastolic blood pressure prior to GDM diagnosis (mean (SD)) |  | 67.29 (9.40) | 70.04 (10.04) | 69.79 (10.02)
Systolic blood pressure prior to GDM diagnosis (mean (SD)) |  | 115.50 (12.67) | 119.96 (13.46) | 119.76 (13.50)
Average diastolic blood pressure prior to pregnancy (mean (SD)) |  | 71.77 (7.93) | 74.15 (8.29) | 73.89 (8.25)
Average systolic blood pressure prior to pregnancy (mean (SD)) |  | 117.91 (10.81) | 121.98 (11.07) | 121.58 (11.08)
Median diastolic blood pressure prior to pregnancy (mean (SD)) |  | 71.75 (8.16) | 74.15 (8.54)) | 73.89 (8.52)
Median systolic blood pressure prior to pregnancy (mean (SD)) |  | 117.62 (11.06) | 121.71 (11.42) | 121.29 (11.43)
Average diastolic blood pressure prior to GDM (mean (SD)) |  | 68.36 (7.96) | 71.06 (8.71) | 70.94 (8.71)
Average systolic blood pressure prior to GDM (mean (SD)) |  | 116.54 (10.92) | 120.99 (11.49) | 120.79 (11.54)
Median diastolic blood pressure prior to GDM (mean (SD)) |  | 68.27 (8.16) | 70.98 (8.89) | 70.88 (8.89)
Median systolic blood pressure prior to GDM (mean (SD)) |  | 116.27 (11.15) | 120.85 (11.67) | 120.60 (11.72)
Infant sex (%) | Female | 4,134 (47.6) | 1,011 (48.3) | 982 (46.9)
| Male | 4528 (52.1) | 1077 (51.5) | 1103 (52.7)
| Unknown | 31 (0.4) | 5 (0.2) | 8 (0.4)
Glucose challenge test value (mean (SD)) |  | 172.47 (23.80) | 176.32 (26.16) | 175.70 (25.98)
Gestational age at GDM diagnosis (mean (SD)) |  | 24.52 (6.57) | 19.39 (7.80) | 19.26 (7.81)
Gestational weight gain up to GDM diagnosis (mean (SD)) |  | 13.63 (11.16) | 9.64 (11.56) | 9.57 (11.33)
Gestational hypertension = yes (%) |  | 634 (7.3) | 250 (11.9) | 231 (11.0)
GDM diagnosed by the C-C criteria = yes (%) |  | 7,107 (81.8) | 1,312 (62.7) | 1,339 (64.0)
One-hour OGTT = abnormal (%) |  | 6,825 (78.5) | 1,419 (67.8) | 1,427 (68.2)
Two-hour OGTT = abnormal (%) |  | 6,146 (70.7) | 936 (44.7) | 953 (45.5)
Three-hour OGTT = abnormal (%) |  | 2,575 (29.6) | 391 (18.7) | 402 (19.2)
OGTT for fasting blood glucose = abnormal (%) |  | 3,095 (35.6) | 1,417 (67.7) | 1,450 (69.3)
GDM severity = severe (%) |  | 936 (10.8) | 279 (13.3) | 276 (13.2)
BMI: body mass index, C-C: Carpenter-Coustan, GDM: gestational diabetes, KP:
Kaiser Permanente,
OGTT: oral glucose tolerance test, SD: standard deviation.
A normal OGTT fasting blood glucose level is lower than 95 mg/dL.
A normal one-hour OGTT blood glucose level is lower than 180 mg/dL.
A normal two-hour OGTT blood glucose level is lower than 155 mg/dL.
A normal three-hour OGTT blood glucose level is lower than 140 mg/dL.

C-section delivery
---
Background variable | Pre-matched bias | Post-matched bias
Propensity score1 | 0.074 | 0.0056
Prognostic score | 0.043 | 0.0025
OGTT for fasting blood glucose | 0.023 | 0.0000
Obese pre-pregnancy BMI | 0.023 | 0.0000
GDM diagnosed by the C-C criteria | 0.017 | 0.0001
Gestational age at GDM diagnosis | 0.013 | 0.0018
Pre-pregnancy pre-diabetes | 0.012 | 0.0005
Normal pre-pregnancy BMI | 0.010 | 0.0001
NICU admission
Background variable | Pre-matched bias | Post-matched bias
Propensity score | 0.035 | 0.0029
GDM diagnosed by the C-C criteria | 0.018 | 0.0000
OGTT for fasting blood glucose | 0.012 | 0.0001
Gestational age at GDM diagnosis | 0.010 | 0.0012
Prognostic score1 | 0.010 | 0.0049

Table 8: Pre-and-post matching comparison of background variables with high
omitted variable bias.

BMI: body mass index, C-C: Carpenter-Coustan, C-section: Cesarean section,
GDM: gestational diabetes, NICU: neonatal intensive care unit, OGTT: oral
glucose tolerance test, OVB: omitted variable bias
1Denotes the maximum post-matched OVB-based bias for that outcome.
The OVB-based bias is the standardized OVB, which is divided by outcome from
the pilot sample.

[Acknowledgments] The authors would like to thank all participants, clinical
staff, Regional Perinatal Service Center of Kaiser Permanente Northern
California. SDP was supported by the Hellman Family Fellowship and by the
National Science Foundation (grant 2142146). LDL was supported by National
Science Foundation Graduate Research Fellowship (grant DGE 2146752). YZ was
supported by the National Institute of Diabetes and Digestive and Kidney
Diseases (grant K01DK120807), National Heart, Lung, and Blood Institute
(R01HL157666), and the Kaiser Permanente Northern California Community
Benefits Program (grant RNG209492). The funders had no role in study design,
data collection and analysis, decision to publish, or preparation of the
manuscript.

## References

  * Abadie, Chingos and West (2018) [author] Abadie, AlbertoA., Chingos, Matthew MM. M. and West, Martin RM. R. (2018). Endogenous stratification in randomized experiments. Review of Economics and Statistics 100 567–580.
  * Ahmed et al. (2006) [author] Ahmed, AliA., Husain, AhsanA., Love, Thomas ET. E., Gambassi, GiovanniG., Dell’Italia, Louis JL. J., Francis, Gary SG. S., Gheorghiade, MihaiM., Allman, Richard MR. M., Meleth, SreelathaS. and Bourge, Robert CR. C. (2006). Heart failure, chronic diuretic use, and increase in mortality and hospitalization: an observational study using propensity score methods. European heart journal 27 1431-1439.
  * Aikens and Baiocchi (2022) [author] Aikens, Rachael CR. C. and Baiocchi, MichaelM. (2022). Assignment-Control Plots: A Visual Companion for Causal Inference Study Design. The American Statistician 1-13.
  * Aikens, Greaves and Baiocchi (2020) [author] Aikens, Rachael CR. C., Greaves, DylanD. and Baiocchi, MichaelM. (2020). A pilot design for observational studies: Using abundant data thoughtfully. Statistics in Medicine 39 4821-4840.
  * Antonelli et al. (2018) [author] Antonelli, JosephJ., Cefalu, MatthewM., Palmer, NathanN. and Agniel, DenisD. (2018). Doubly robust matching estimators for high dimensional confounding adjustment. Biometrics 74 1171-1179.
  * American Diabetes Association (2021) [author] American Diabetes Association (2021). 14\. Management of diabetes in pregnancy: standards of medical care in diabetes-2021. Diabetes Care 44 S200-S210.
  * Athey, Imbens and Wager (2018) [author] Athey, SusanS., Imbens, Guido WG. W. and Wager, StefanS. (2018). Approximate residual balancing: debiased inference of average treatment effects in high dimensions. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80 597–623.
  * Austin (2009) [author] Austin, Peter CP. C. (2009). Balance diagnostics for comparing the distribution of baseline covariates between treatment groups in propensity-score matched samples. Statistics in medicine 28 3083-3107.
  * Behrashi et al. (2016) [author] Behrashi, MitraM., Samimi, MansoorehM., Ghasemi, TayyebehT., Saberi, FarzanehF. and Atoof, FatemehF. (2016). Comparison of glibenclamide and insulin on neonatal outcomes in pregnant women with gestational diabetes. International Journal of Preventive Medicine 7\.
  * Ben-Michael, Feller and Rothstein (2021) [author] Ben-Michael, EliE., Feller, AviA. and Rothstein, JesseJ. (2021). Varying impacts of letters of recommendation on college admissions: Approximate balancing weights for subgroup effects in observational studies. arXiv preprint arXiv:2008.04394.
  * Ben-Michael et al. (2021) [author] Ben-Michael, EliE., Feller, AviA., Hirshberg, David AD. A. and Zubizarreta, José RJ. R. (2021). The Balancing Act in Causal Inference. arXiv preprint arXiv:2110.14831.
  * Bennett, Vielma and Zubizarreta (2020) [author] Bennett, MagdalenaM., Vielma, Juan PabloJ. P. and Zubizarreta, José RJ. R. (2020). Building representative matched samples with multi-valued treatments in large observational studies. Journal of Computational and Graphical Statistics 29 744–757.
  * Brooks and Ohsfeldt (2013) [author] Brooks, John MJ. M. and Ohsfeldt, Robert LR. L. (2013). Squeezing the balloon: propensity scores and unmeasured covariate balance. Health services research 48 1487–1507.
  * Castillo et al. (2014) [author] Castillo, Wendy CameloW. C., Boggess, KimK., Stürmer, TilT., Brookhart, M AlanM. A., Benjamin Jr, Daniel KD. K. and Funk, Michele JonssonM. J. (2014). Trends in glyburide compared with insulin use for gestational diabetes treatment in the United States, 2000-2011. Obstetrics and gynecology 123 1177\.
  * Cinelli and Hazlett (2020) [author] Cinelli, CarlosC. and Hazlett, ChadC. (2020). Making sense of sensitivity: Extending omitted variable bias. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 82 39-67.
  * Cox (2007) [author] Cox, DRD. (2007). On a generalization of a result of WG Cochran. Biometrika 94 755–759.
  * Diamond and Sekhon (2013) [author] Diamond, AlexisA. and Sekhon, Jasjeet SJ. S. (2013). Genetic matching for estimating causal effects: A general multivariate matching method for achieving balance in observational studies. Review of Economics and Statistics 95 932–945.
  * Ding, VanderWeele and Robins (2017) [author] Ding, PengP., VanderWeele, TJT. and Robins, James MJ. M. (2017). Instrumental variables as bias amplifiers with general outcome and confounding. Biometrika 104 291–302.
  * Ferrara et al. (2004) [author] Ferrara, AssiamiraA., Kahn, Henry SH. S., Quesenberry, Charles PC. P., Riley, CandiceC. and Hedderson, Monique MM. M. (2004). An increase in the incidence of gestational diabetes mellitus: Northern California, 1991-2000. Obstetrics & Gynecology 103 526-533.
  * Fisher (1936) [author] Fisher, Ronald AylmerR. A. (1936). Design of experiments. British Medical Journal 1 554\.
  * Glass (1976) [author] Glass, Gene VG. V. (1976). Primary, secondary, and meta-analysis of research. Educational researcher 5 3–8.
  * Goetzl and Wilkins (2002) [author] Goetzl, LauraL. and Wilkins, IsabelleI. (2002). Glyburide compared to insulin for the treatment of gestational diabetes mellitus: a cost analysis. Journal of Perinatology 22 403-406.
  * Gordon and Lin (2016) [author] Gordon, NancyN. and Lin, TeresaT. (2016). The Kaiser Permanente Northern California adult member health survey. The Permanente Journal 20\.
  * Greifer (2021) [author] Greifer, NoahN. (2021). cobalt: Covariate Balance Tables and Plots R package version 4.3.1.
  * Greifer and Stuart (2021) [author] Greifer, NoahN. and Stuart, Elizabeth AE. A. (2021). Choosing the estimand when matching or weighting in observational studies. arXiv preprint arXiv:2106.10577.
  * Hainmueller (2012) [author] Hainmueller, JensJ. (2012). Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies. Political analysis 20 25–46.
  * Hansen (2008) [author] Hansen, Ben BB. B. (2008). The prognostic analogue of the propensity score. Biometrika 95 481-488.
  * Hansen and Bowers (2008) [author] Hansen, Ben BB. B. and Bowers, JakeJ. (2008). Covariate balance in simple, stratified and clustered comparative studies. Statistical Science 219-236.
  * Hedges (1981) [author] Hedges, Larry VL. V. (1981). Distribution theory for Glass’s estimator of effect size and related estimators. journal of Educational Statistics 6 107–128.
  * Imbens and Rubin (2015) [author] Imbens, Guido WG. W. and Rubin, Donald BD. B. (2015). Causal inference in statistics, social, and biomedical sciences. Cambridge University Press.
  * King and Nielsen (2019) [author] King, GaryG. and Nielsen, RichardR. (2019). Why propensity scores should not be used for matching. Political Analysis 27 435–454.
  * Li and Greene (2013) [author] Li, LiangL. and Greene, TomT. (2013). A weighting analogue to pair matching in propensity score analysis. The international journal of biostatistics 9 215–234.
  * Li, Thomas and Li (2019) [author] Li, FanF., Thomas, Laine EL. E. and Li, FanF. (2019). Addressing extreme propensity scores via the overlap weights. American journal of epidemiology 188 250–257.
  * McNemar (1947) [author] McNemar, QuinnQ. (1947). Note on the sampling error of the difference between correlated proportions or percentages. Psychometrika 12 153–157.
  * Meng (2018) [author] Meng, Xiao-LiX.-L. (2018). Statistical paradises and paradoxes in big data (I): law of large populations, big data paradox, and the 2016 US presidential election. The Annals of Applied Statistics 12 685–726.
  * Neyman (1923) [author] Neyman, JerzyJ. (1923). On the application of probability theory to agricultural experiments. Essay on principles. Section 9. Translated in Statistical Science 5 (1990) 465-480.
  * American College of Obstetricians and Gynecologists (2018) [author] American College of Obstetricians and Gynecologists (2018). ACOG practice bulletin no. 190: gestational diabetes mellitus. Obstet Gynecol 131 e49-e64.
  * Pearson (1895) [author] Pearson, KarlK. (1895). VII. Note on regression and inheritance in the case of two parents. proceedings of the royal society of London 58 240–242.
  * Pimentel (2016) [author] Pimentel, Samuel DS. D. (2016). Large, Sparse Optimal Matching with R package rcbalance. Observational Studies 2 4-23.
  * Pimentel (2022) [author] Pimentel, Samuel DS. D. (2022). Covariate-adaptive randomization inference in matched designs. arXiv preprint arXiv:2207.05019.
  * Pimentel and Kelz (2020) [author] Pimentel, Samuel DS. D. and Kelz, Rachel RR. R. (2020). Optimal tradeoffs in matched designs comparing US-trained and internationally trained surgeons. Journal of the American Statistical Association 115 1675–1688.
  * Pimentel, Small and Rosenbaum (2016) [author] Pimentel, Samuel DS. D., Small, Dylan SD. S. and Rosenbaum, Paul RP. R. (2016). Constructed second control groups and attenuation of unmeasured biases. Journal of the American Statistical Association 111 1157-1167.
  * Pimentel et al. (2015) [author] Pimentel, Samuel DS. D., Kelz, Rachel RR. R., Silber, Jeffrey HJ. H. and Rosenbaum, Paul RP. R. (2015). Large, sparse optimal matching with refined covariate balance in an observational study of the health outcomes produced by new surgeons. Journal of the American Statistical Association 110 515-527.
  * Rochon et al. (2006) [author] Rochon, MeredithM., Rand, LarryL., Roth, LisaL. and Gaddipati, SreedharS. (2006). Glyburide for the management of gestational diabetes: risk factors predictive of failure and associated pregnancy outcomes. American journal of obstetrics and gynecology 195 1090-1094.
  * Rosenbaum (2005) [author] Rosenbaum, Paul RP. R. (2005). Heterogeneity and causality: Unit heterogeneity and design sensitivity in observational studies. The American Statistician 59 147–152.
  * Rosenbaum (2010) [author] Rosenbaum, Paul RP. R. (2010). Design of observational studies 10\. Springer.
  * Rosenbaum (2012) [author] Rosenbaum, Paul RP. R. (2012). Optimal matching of an optimally chosen subset in observational studies. Journal of Computational and Graphical Statistics 21 57–71.
  * Rosenbaum (2020) [author] Rosenbaum, Paul RP. R. (2020). Modern algorithms for matching in observational studies. Annual Review of Statistics and Its Application 7 143–176.
  * Rosenbaum, Ross and Silber (2007) [author] Rosenbaum, Paul RP. R., Ross, Richard NR. N. and Silber, Jeffrey HJ. H. (2007). Minimum distance matched sampling with fine balance in an observational study of treatment for ovarian cancer. Journal of the American Statistical Association 102 75–83.
  * Rosenbaum and Rubin (1983) [author] Rosenbaum, Paul RP. R. and Rubin, Donald BD. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika 70 41-55.
  * Rosenbaum and Rubin (1985) [author] Rosenbaum, Paul RP. R. and Rubin, Donald BD. B. (1985). The bias due to incomplete matching. Biometrics 103-116.
  * Rosenbaum and Silber (2001) [author] Rosenbaum, Paul RP. R. and Silber, Jeffrey HJ. H. (2001). Matching and thick description in an observational study of mortality after surgery. Biostatistics 2 217–232.
  * Rubin (1974) [author] Rubin, Donald BD. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of educational Psychology 66 688\.
  * Rubin (2007) [author] Rubin, Donald BD. B. (2007). The design versus the analysis of observational studies for causal effects: parallels with the design of randomized trials. Statistics in medicine 26 20-36.
  * Sales, Hansen and Rowan (2018) [author] Sales, Adam CA. C., Hansen, Ben BB. B. and Rowan, BrianB. (2018). Rebar: Reinforcing a matching estimator with predictions from high-dimensional covariates. Journal of Educational and Behavioral Statistics 43 3–31.
  * Soriano et al. (2021) [author] Soriano, DanD., Ben-Michael, EliE., Bickel, Peter JP. J., Feller, AviA. and Pimentel, Samuel DS. D. (2021). Interpretable sensitivity analysis for balancing weights. arXiv preprint arXiv:2102.13218.
  * Stekhoven and Bühlmann (2012) [author] Stekhoven, Daniel JD. J. and Bühlmann, PeterP. (2012). MissForest—non-parametric missing value imputation for mixed-type data. Bioinformatics 28 112–118.
  * Stuart (2010) [author] Stuart, Elizabeth AE. A. (2010). Matching methods for causal inference: A review and a look forward. Statistical science: a review journal of the Institute of Mathematical Statistics 25 1\.
  * Stuart et al. (2011) [author] Stuart, Elizabeth AE. A., King, GaryG., Imai, KosukeK. and Ho, DanielD. (2011). MatchIt: nonparametric preprocessing for parametric causal inference. Journal of statistical software.
  * Wang and Zubizarreta (2020) [author] Wang, YixinY. and Zubizarreta, Jose RJ. R. (2020). Minimal dispersion approximately balancing weights: asymptotic properties and practical considerations. Biometrika 107 93–105.
  * Yang et al. (2012) [author] Yang, DanD., Small, Dylan SD. S., Silber, Jeffrey HJ. H. and Rosenbaum, Paul RP. R. (2012). Optimal matching with minimal deviation from fine balance in a study of obesity and surgical outcomes. Biometrics 68 628–636.
  * Zhang and Zhao (2022) [author] Zhang, YaoY. and Zhao, QingyuanQ. (2022). What is a randomization test? arXiv preprint arXiv:2203.10980.
  * Zhao and Yang (2022) [author] Zhao, HongheH. and Yang, ShuS. (2022). Outcome-adjusted balance measure for generalized propensity score model selection. Journal of Statistical Planning and Inference 221 188–200.
  * Zhu and Zhang (2016) [author] Zhu, YeyiY. and Zhang, CuilinC. (2016). Prevalence of gestational diabetes and risk of progression to type 2 diabetes: a global perspective. Current diabetes reports 16 1-11.
  * Zubizarreta (2015) [author] Zubizarreta, José RJ. R. (2015). Stable weights that balance covariates for estimation with incomplete outcome data. Journal of the American Statistical Association 110 910–922.
  * Zubizarreta, Paredes and Rosenbaum (2014) [author] Zubizarreta, José RJ. R., Paredes, Ricardo DR. D. and Rosenbaum, Paul RP. R. (2014). Matching for balance, pairing for heterogeneity in an observational study of the effectiveness of for-profit and not-for-profit high schools in Chile. The Annals of Applied Statistics 8 204–231.