L0pt

# Scalable kernel balancing weights in a nationwide observational study of
hospital profit status and heart attack outcomes††thanks:  This work was
supported through a grant from the Alfred P. Sloan Foundation (G-2020-13946)
and an award from the Patient Centered Outcomes Research Initiative (PCORI,
ME-2022C1- 25648). The authors report there are no competing interests to
declare. This work was completed while Kwangho Kim was a research associate at
Harvard Medical School.

Kwangho Kim Department of Statistics, Korea University, 145 Anam-ro, Seongbuk-
gu, Seoul, 02841, South Korea; email<EMAIL_ADDRESS>Bijan A. Niknam
Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health,
615 N Wolfe St, Baltimore, MD 21205, United States; email<EMAIL_ADDRESS>José R. Zubizarreta Departments of Health Care Policy, Biostatistics, and
Statistics, Harvard University, 180 Longwood Avenue, Office 307-D, Boston, MA
02115, United States; email<EMAIL_ADDRESS>
###### Abstract

Weighting is a general and often-used method for statistical adjustment.
Weighting has two objectives: first, to balance covariate distributions, and
second, to ensure that the weights have minimal dispersion and thus produce a
more stable estimator. A recent, increasingly common approach directly
optimizes the weights toward these two objectives. However, this approach has
not yet been feasible in large-scale datasets when investigators wish to
flexibly balance general basis functions in an extended feature space. For
example, many balancing approaches cannot scale to national-level health
services research studies. To address this practical problem, we describe a
scalable and flexible approach to weighting that integrates a basis expansion
in a reproducing kernel Hilbert space with state-of-the-art convex
optimization techniques. Specifically, we use the rank-restricted Nyström
method to efficiently compute a kernel basis for balancing in nearly linear
time and space, and then use the specialized first-order alternating direction
method of multipliers to rapidly find the optimal weights. In an extensive
simulation study, we provide new insights into the performance of weighting
estimators in large datasets, showing that the proposed approach substantially
outperforms others in terms of accuracy and speed. Finally, we use this
weighting approach to conduct a national study of the relationship between
hospital profit status and heart attack outcomes in a comprehensive dataset of
1.27 million patients. We find that for-profit hospitals use interventional
cardiology to treat heart attacks at similar rates as other hospitals, but
have higher mortality and readmission rates.

Keywords: Causal Inference; Observational Studies; Weighting; Convex
Optimization; Propensity Score

##  1 Introduction

###  1.1 Hospital profit status and heart attack care

For-profit organizations, including hospitals, are widespread across the US
healthcare system. The impact of profit motive on treatment and outcomes is a
major focus of healthcare research. For-profit status may possibly incentivize
hospitals to provide more efficient care, leading to better patient outcomes.
For-profit hospitals may also invest more in technology and human capital to
improve competitiveness and potentially patient outcomes. Alternatively,
higher reimbursement for interventional treatment may possibly incentivize
them to use invasive treatments more frequently than medically necessary. For
example, Medicare pays hospitals more for heart attack patients who receive
percutaneous coronary intervention (PCI), a treatment that is essential in
some cases but of ambiguous benefit in others, potentially incentivizing
unnecessary use. Previous studies of subgroups of hospitals reached
conflicting conclusions on both PCI treatment rates and patient outcomes
(Sloan et al. 2003; Shah et al. 2007), emphasizing the need for a
representative study of a comprehensive national sample.

###  1.2 Weighting for covariate adjustment in large-scale datasets

Covariate adjustment is essential in observational studies to remove bias due
to differences in observed confounders. One method for adjustment is
weighting, which uses each unit’s probability of treatment assignment to form
weighted samples that are free from differences due to observed confounders
(Austin and Stuart 2015). Weighting attempts to remove bias by balancing
observed covariate distributions while also producing robust and stable
estimators. An appealing property of weighting is that it does not require
explicit modeling of the outcome, and hence is part of the design stage of the
study (Rubin 2006). Further, a single set of weights can in principle be used
to examine multiple outcomes (Little and Rubin 2019).

Procedurally, there are two broad approaches to weighting (Ben-Michael et al.
2021). The modeling approach emphasizes the fit of a model for the probability
of treatment assignment given observed covariates, or propensity score
(Rosenbaum and Rubin 1983). However, covariate balance may be inadequate due
to model misspecification or small samples; moreover, the weights may be
highly variable and produce unstable estimators (Kang and Schafer 2007;
Zubizarreta 2015). To address these challenges, the balancing approach
attempts to balance covariates directly, often by solving a convex
optimization problem to find the weights of minimum dispersion that balance
functions of the covariates (see, e.g., Hainmueller 2012, Zubizarreta 2015).
This affords control over the form and degree of balance while also targeting
estimator stability. Thus, the balancing approach can have better finite-
sample performance (Chattopadhyay et al. 2020).

Despite these appealing properties, the computational cost of many weighting
algorithms prevents their scaling to the massive dimensions of many modern
data sources in health sciences, which can reach millions of observations and
hundreds of covariates. Moreover, to improve estimation accuracy, researchers
increasingly recognize the importance of balancing richer covariate function
classes beyond simple means and other low-dimensional parametric summaries.
However, this additional complexity makes execution of weighting methods at
scale even more challenging.

###  1.3 Contribution and outline

We propose a weighting approach that addresses practical barriers to the use
of balancing methods in big data. Following Wong and Chan (2018) and Hazlett
(2020), we adopt a kernel balancing approach which assumes the outcome
regression functions are in a flexible function space represented by a kernel.
Further, we use the recently developed rank-restricted Nyström approximation
(Wang et al. 2019) that gives a provable approximate solution to the efficient
kernel basis calculation problem. This allows us to construct the kernel basis
functions in near-linear time and space, with strong relative-error
guarantees. We then incorporate these functions into a quadratic program with
linear constraints, which are efficiently solved via a specialized first-order
alternating direction method of multipliers (Stellato et al. 2020), yielding
the stable kernel balancing weights of minimum dispersion that approximately
balance the covariate functions (Zubizarreta 2015). We analyze the worst-case
bias bound of the resulting linear estimator and discuss its implications for
parameter selection. We benchmark the proposed weighting estimator’s
performance against common alternatives in extensive simulation studies,
providing new insights into how weighting estimators perform at scale, finding
that the proposed approach considerably outperforms others in terms of both
accuracy and speed. Finally, we use the new approach in a national study of
heart attack treatment and outcomes by hospital profit status.

The paper outline is as follows. Section 2 describes the setup, the balancing
approach to weighting, and kernel balancing. Section 3 details the rank-
restricted Nyström approximation, the operator splitting solver for quadratic
programs, and the proposed scalable kernel balancing approach integrating the
two concepts, along with an analysis of the worst-case bias bound for the
resulting weighting estimator. Section 4 outlines the simulation study
settings, examines the performance of the proposed estimator, and provides
general insights into the speed and accuracy of various weighting estimators.
Section 5 examines the relationship between hospital profit status and heart
attack treatment and outcomes in a national Medicare dataset. Section 6
concludes.

##  2 Estimation framework

###  2.1 Two approaches to weighting

Consider an observational study with $n$ triplets
$(X,A,Y)\overset{iid}{\sim}\mathbb{P}$, where
$X\in\mathcal{X}\subset\mathbb{R}^{d}$ is a vector of observed covariates, $A$
is a binary treatment assignment indicator with the number of treated units
$n_{t}=\sum_{i=1}^{n}A_{i}$ and control units $n_{c}=n-n_{t}$, and
$Y\in\mathbb{R}$ is a real-valued outcome. Under the potential outcomes
framework for causal inference (Rubin 1974), we write
$Y_{i}=A_{i}Y_{i}^{1}+(1-A_{i})Y_{i}^{0}$, where $Y_{i}^{a}$ is the potential
outcome for unit $i$ under treatment $A=a$; $a=0,1$. In our case study, $A=1$
if a patient was admitted to a for-profit hospital, and $A=0$ otherwise. Here,
we wish to estimate the average treatment effect on the treated (ATT)111While
we focus on estimating the ATT, our discussion also applies to the average
treatment effect (ATE), and more generally to the target average treatment
effect (TATE), $\mathbb{E}_{\mathbb{T}}(Y^{1}-Y^{0})$, where $\mathbb{T}$
represents the measure of the target population of interest (Kern et al.
2016). The ATE and ATT are special cases of the TATE where
$\mathbb{T}=\mathbb{P}$ and $\mathbb{T}=\mathbb{P}(A=1)$, respectively.,

$\mathbb{E}(Y^{1}-Y^{0}\mid A=1).$ (1)

In order to identify the ATT from observational data, we invoke the following
assumptions (Rosenbaum and Rubin 1983):

Assumption 2.1: (i) Consistency, $Y=Y^{a}$ if $A=a$; (ii) Exchangeability,
$A\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}Y^{0}\mid X$; (iii) Positivity,
$\mathbb{P}(A=0\mid X=x)>0$ a.s. $[\mathbb{P}]$.

Given these assumptions, the ATT in (1) is identified by

$\displaystyle\mathbb{E}(Y\mid A=1)-\mathbb{E}\\{\mathbb{E}(Y\mid X,A=0)\mid
A=1\\}\equiv\mathbb{E}(Y\mid A=1)-\psi.$

The first term can be estimated easily by taking a sample average. However,
estimating $\psi$, the mean conditional potential outcomes of the treated
under control, requires not only the aforementioned assumptions but also
careful covariate adjustment.

A simple yet general weighting estimator for $\psi$ is given by

$\widehat{\psi}=\sum_{\\{i\mid
A=0\\}}\hat{w}_{i}Y_{i}\quad\text{for}\quad\hat{w}\geq 0.$ (2)

where $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{n_{c}})\in\mathbb{R}^{n_{c}}$ is
the set of estimated weights that removes imbalances in covariate
distributions such that after weighting adjustment, the control sample
resembles the treated sample in aggregate.

In the modeling approach to weighting, we explicitly model the unknown
propensity score $\pi(X)=\mathbb{P}(A=1\mid X)$, and then construct the
weights $\hat{w}$ using the estimated $\widehat{\pi}$. One of the most widely-
used modeling-based weighting estimators is the _Hajek estimator_ defined by

$\displaystyle\widehat{\psi}_{\text{mod}}=\frac{1}{n}\sum_{i=1}^{n}\left\\{\frac{\widehat{\pi}(X_{i})(1-A_{i})}{1-\widehat{\pi}(X_{i})}Y_{i}\right\\}\Bigg{/}\frac{1}{n}\sum_{i=1}^{n}\left\\{\frac{\widehat{\pi}(X_{i})(1-A_{i})}{1-\widehat{\pi}(X_{i})}\right\\}.$
(3)

The performance of model-based weighting estimators hinges upon the accuracy
of estimation of $\widehat{\pi}$. Often, flexible nonparametric models are
employed for $\widehat{\pi}$, but these approaches encounter difficulties when
scaled to very large datasets (e.g., Zhou et al. 2017).

The balancing approach to weighting is a more recent alternative approach
which avoids explicit modeling of $\pi$ and rather directly optimizes for the
weights. This approach can directly target not only covariate balance, but
also the dispersion of the weights, thereby yielding a more stable estimator.
Consequently, the balancing approach has been found to provide better finite-
sample performance than the modeling approach (Chattopadhyay et al. 2020).
Here, we focus on the stable balancing weights (SBW) (Zubizarreta 2015), which
solve the following quadratic program

$\displaystyle\underset{w\in\mathbb{R}^{n_{c}}}{\text{minimize}}\quad\sum_{\\{i\mid
A=0\\}}(w_{i}-\bar{w})^{2}$ (4) $\displaystyle\text{subject to
}\quad\left|\sum_{\\{i\mid
A=0\\}}w_{i}B_{b}(X_{i})-\frac{1}{n_{t}}\sum_{\\{i\mid
A=1\\}}B_{b}(X_{i})\right|\leq\delta,\quad b=1,2,...,M,$
$\displaystyle\phantom{\text{subject to }}\quad\,\sum_{\\{i\mid
A=0\\}}w_{i}=1,\quad w\geq 0,$

where $\bar{w}=\frac{1}{n_{c}}\sum_{\\{i\mid A=0\\}}w_{i}$, $B_{1},...,B_{M}$
are $M$ real-valued basis functions of the covariates that span the model
space in which the unknown outcome regression function is assumed to lie, and
$\delta$ is a prespecified tolerance level for the maximum acceptable
imbalance. We use $\widehat{\psi}_{\text{bal}}$ to denote the weighting
estimator when the set of weights $\hat{w}$ represents the solution of (4). We
focus on the stable balancing weights for two reasons. First, the objective
function targets the variance of the weights and hence (4) is a convex
quadratic program (QP) to which we may apply a state-of-the-art algorithm to
improve scalability. Second, the optimal weights have minimum variance and
balance the covariates, thereby controlling the bias and variance of the
resulting estimator, and achieving the semiparametric efficiency bound under
certain regularity conditions (Wang and Zubizarreta 2020). Further,
approximate rather than exact balancing makes the approach’s performance less
sensitive to settings with limited covariate overlap. The parameter $\delta$
controls the trade-offs between covariate balance and the variance of the
weights, and therefore the bias-variance trade-off of the resulting estimator.
A smaller value of $\delta$ reduces imbalances and therefore bias due to
differences in observed covariates, while a larger $\delta$ produces a set of
weights with smaller variance and therefore a more stable estimator. In
practice, subject-matter experts are typically consulted to ensure the choice
of $\delta$ is suitable for the application at hand. However, some recent work
has explored how to choose $\delta$ in an empirical fashion. For example,
appealing to asymptotic theory, Hirshberg and Wager (2020) and Hirshberg et
al. (2021) recommend selecting $\delta$ in the context of the dual formulation
of the weighting problem and setting it to be roughly equivalent to the
conditional variance of the outcome. Other automatic selection approaches have
also been suggested (Kallus 2020; Wang and Zubizarreta 2020; Zhao 2019);
however, more theoretical and empirical work remains to fully characterize
their performance.

###  2.2 Kernel balancing

Let $\mu_{0}(X)=\mathbb{E}(Y^{0}\mid X)$ be the outcome regression function in
some model space $\mathcal{M}$. It can be shown that the bias of the weighting
estimator (2) with the weights $\hat{w}$ such that

$\displaystyle\max_{\mu_{0}\in\mathcal{M}}\left|\sum_{\\{i\mid
A=0\\}}\hat{w}_{i}\mu_{0}(X_{i})-\frac{1}{n_{t}}\sum_{\\{i\mid
A=1\\}}\mu_{0}(X_{i})\right|\leq\delta$ (5)

is on the order of $\delta$ (Ben-Michael et al. 2021). Hence, to maximize the
chance that the weighting adjustment eliminates the covariate imbalance in
$\mu_{0}$, it is desirable to choose a flexible set of basis functions
$\\{B_{1},...,B_{M}\\}$ that spans a general model space for $\mu_{0}$. One
such space is the Reproducing Kernel Hilbert Space (RKHS). For a Mercer kernel
$K:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ and a given sample of
size $n$, the RKHS $\mathcal{H}_{K}$ is defined by the completion of the
function space

$\left\\{f:f(x)=\sum_{i=1}^{n}\alpha_{i}K(X_{i},x),\,\alpha_{i}\in\mathbb{R}\right\\}$

with respect to the norm
$\|f\|_{\mathcal{H}_{K}}=\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_{i}\alpha_{j}K(X_{i},X_{j})}=\sqrt{\alpha^{\top}\bm{K}\alpha}$,
$\forall f\in\mathcal{H}_{K}$, where $\alpha=(\alpha_{1},\ldots,$
$\alpha_{n})^{\top}\in\mathbb{R}^{n}$ and $\bm{K}$ is the $n\times n$ kernel
matrix in our observational data with $\bm{K}_{ij}=K(X_{i},X_{j})$. For a
clearer conceptualization of the function space $\mathcal{H}_{K}$, define a
feature map
$\Phi(x):x\mapsto[\sqrt{\lambda_{1}}\varphi_{1}(x),\sqrt{\lambda_{2}}\varphi_{2}(x),\ldots]$,
where $\\{\lambda_{j}\\}$ and $\\{\varphi_{j}\\}$ are, respectively, the
eigenvalues and orthonormal eigenfunctions of the kernel operator such that
$\int K(x,y)\varphi_{j}(y)dy=\lambda_{j}\varphi_{j}(x)$. Then by Mercer’s
Theorem it can be shown that $\Phi(x)^{\top}\Phi(y)=K(x,y)$, and that for any
$f\in\mathcal{H}_{K}$, $f$ can be expanded either in terms of $K$ or in terms
of the bases $\\{\varphi_{j}(\cdot)\\}_{j=1}^{\infty}$, such that

$\displaystyle
f(x)=\sum_{i=1}^{n}\alpha_{i}K(X_{i},x)=\sum_{j=1}^{\infty}\beta_{j}\varphi_{j}(x)$
(6)

for some $\alpha_{i},\beta_{j}\in\mathbb{R}$. In the kernel balancing
approach, we assume that $\mu_{0}\in\mathcal{H}_{K}$ as follows.

Assumption 2.2: $\mathcal{M}=\mathcal{H}_{K}$ for a Mercer kernel $K$ and its
induced RKHS $\mathcal{H}_{K}$.

Letting the balancing constraints in (4) be $B_{b}(\cdot)=\varphi_{b}(\cdot)$,
$b=1,\ldots,\infty$, we pursue approximate mean balance on the basis set for
$\mathcal{H}_{K}$, $\\{\varphi_{j}(\cdot)\\}_{j=1}^{\infty}$, to control bias
according to (5). By (6), this is equivalent to using
$B_{b}(\cdot)=\bm{K}_{b}^{\top}$, $b=1,\ldots,n$, where $\bm{K}_{i}$ denotes
the $i$-th row of $\bm{K}$. In other words, one may use the $n$ columns of
$\bm{K}$ as bases for $\mu_{0}$ rather than
$\\{\varphi_{j}(\cdot)\\}_{j=1}^{\infty}$. Note that $\bm{K}$ can be
constructed without computing $\varphi_{b}$. Achieving mean balance on columns
of $\bm{K}$ can also guarantee covariate balance in $\widehat{\mu}_{0}$. If we
find $\widehat{\mu}_{0}$ as the solution to the regularized problem

$\displaystyle\underset{\mu_{0}\in\mathcal{H}_{K}}{\min}\sum_{i=1}^{n}\mathcal{L}(\mu_{0}(X_{i}),Y^{0}_{i})+\tau\Omega(\|\mu_{0}\|_{\mathcal{H}_{K}}),$
(7)

for some loss function $\mathcal{L}:\mathbb{R}^{2}\rightarrow\mathbb{R}$,
monotone increasing function $\Omega$, and $\tau>0$, then by the representer
theorem (Wahba 1990), $\widehat{\mu}_{0}$ always has the form
$\widehat{\mu}_{0}(x)=\sum_{i=1}^{n}\alpha^{*}_{i}K(X_{i},x)=\sum_{j=1}^{\infty}{\beta}^{*}_{j}\varphi_{j}(x)$
for some $\alpha^{*}_{i},{\beta}^{*}_{j}\in\mathbb{R}$. Hence, the same
argument applies as above.

However, incorporating $n$ balancing constraints into the optimization is
challenging for large $n$. Admitting the spectral decomposition, one may write
$\bm{K}=U\Lambda U^{\top}$ where $\Lambda$ is the diagonal matrix with a set
of eigenvalues in decreasing order. Then under Assumption (2.2), the bias of
the linear estimator (2) due to imbalance on $\bm{K}$ is given by

$\displaystyle\text{Bias}(\widehat{\psi};\bm{K})\coloneqq
C_{0}\left(\widehat{w}^{\top}U_{c}-\frac{1}{n_{t}}\mathbb{1}_{n_{t}}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha_{0},$ (8)

for some $C_{0}>0$ and $\alpha_{0}$ such that
$\|\mu_{0}\|_{\mathcal{H}_{K}}=\sqrt{\alpha_{0}^{\top}\bm{K}\alpha_{0}}$,
where $U_{c}$ ($U_{t}$) are rows of $U$ corresponding to the control (treated)
units and $\mathbb{1}_{n_{t}}$ is the vector of ones of length $n_{t}$. Based
on this finding, Hazlett (2020) proposed balancing the first $r\ll n$
eigenvectors of $\bm{K}$, i.e., the first $r$ columns of $U$, rather than
balancing all $n$ columns of $\bm{K}$. While easier than balancing all $n$
columns, this can still be prohibitive due to the high computational cost of
the eigenvalue decomposition; standard algorithms for the eigenvalue
decomposition of a dense $n\times n$ matrix require $O(n^{2})$-space and have
$O(n^{3})$-time complexity. In the next section, we discuss an effective
approximation method for $\bm{K}$ and a scalable strategy for the
corresponding basis function construction.

######  Remark 1 (Positivity under a model class).

Although positivity is not strictly required for identification under the
assumption that $\mathcal{M}=\mathcal{H}_{K}$, it may play a role in
estimation by influencing the feasibility of finding the weights in the
balancing problem. At the extreme, one may not feasibly find weights if no
control units have similar rows of $\bm{K}$ as a given treated unit. However,
the notion of positivity in this setting is more subtle than the standard
positivity assumption, as there are several factors simultaneously affecting
the feasibility of weighting solutions, including the type of kernel,
bandwidth, tolerance level, and the observed sample at hand.

######  Remark 2 (Considerations for Kernel selection).

For a given a kernel $K$, the corresponding RKHS is unique. One desirable
property for $K$ is _universal approximation_ , which allows us to approximate
an arbitrary continuous function uniformly on the feature space induced by
$\Phi$ (Micchelli et al. 2006). Following Hazlett (2020), we use the Gaussian
kernel, which satisfies the universal approximation property. However, other
kernels could also be considered depending on prior knowledge of the
properties of the outcome regression function (see Micchelli et al., 2006,
Section 3, for examples of universal kernels). When setting the kernel
bandwidth, we likewise follow Hazlett (2020) by setting it to the column rank
of $X$. Like the type of kernel, the chosen bandwidth may also affect both the
performance and feasibility of the kernel balancing estimator. The optimal
choices of the kernel and its bandwidth are important considerations in causal
inference that warrant further investigation.

##  3 Toward scalable and flexible weighting

###  3.1 Scalable kernel basis calculation using the rank-restricted Nyström
approximation

The kernel approach provides a flexible way to balance functions of the
covariates. However, its expensive computational cost presently prevents its
use in large datasets. This problem can be surmounted with low-rank matrix
approximations. Rather than compute the eigenvalue decomposition of the entire
kernel, one may instead perform a partial singular value decomposition (SVD)
(Lehoucq et al. 1998). However, the performance of this approach depends
heavily on the singular spectrum structure of $\bm{K}$, and only realizes
computational gains when $\bm{K}$ is sparse. Another alternative is to use a
class of randomized algorithms such as the randomized QR decomposition or
randomized SVD (Halko et al. 2011), but despite having good convergence, these
algorithms still generally require $O(n^{2})$ time.

A more general alternative is the Nyström method (e.g., Gittens and Mahoney
2016), an efficient technique for low-rank approximation of a full kernel
matrix based on a subset of its columns. The method can substantially reduce
the time complexity since only a portion of $\bm{K}$ is sampled for
computation. However, when $n$ is extremely large, even a small matrix may no
longer be small in practical terms, and the accuracy of the matrix
approximation drastically declines if the number of sampled columns is too
small. In this context, Pourkamali-Anaraki et al. (2018) note the limitations
of traditional Nyström methods for very large data sets, a problem we address
using a more recently introduced form of the Nyström method.

We adopt the rank-restricted Nyström method recently proposed by Wang et al.
(2019). Unlike traditional Nyström methods, the rank-restricted method allows
the target rank to be much smaller than the number of sampled columns while
still providing guarantees on the relative-error bound. For
$\bm{K}\in\mathbb{R}^{n\times n}$, we first sample a set of $m\ll n$ column
indices $\mathcal{I}_{m}$ (with the magnitude of sampled columns $m$ also
termed a _sketch size_) . While various sampling or _sketching_ schemes have
been proposed (see, e.g., Kumar et al. 2012, and Table 1), we use uniform
sampling here as its minimal computational cost per sampled column guarantees
$O(n)$ time and hence has advantages for scalability. After sampling columns
of $\bm{K}$, we form a matrix $C\in\mathbb{R}^{n\times c}$ with
$C_{ij}=K(X_{i},X_{j})$, $i\in\\{1,...,n\\}$, $j\in\mathcal{I}_{m}$, and a
matrix $W\in\mathbb{R}^{m\times m}$ with $W_{ij}=K(X_{i},X_{j})$,
$i,j\in\mathcal{I}_{m}$. Then, the standard Nyström approximation for $\bm{K}$
is

$\displaystyle\widetilde{\bm{K}}=CW^{+}C^{\top},$ (9)

where $W^{+}$ is the Moore-Penrose (pseudo) inverse of $W$.

Sketching method | Sketch size ($c$) | Time complexity
---|---|---
Uniform sampling | $O\left(\frac{s\mu}{\epsilon}+s\mu\log s\right)$ | $O(nc)$
Leverage sampling | $O\left(\frac{s}{\epsilon}+s\log s\right)$ | $O(n^{2}s)$
Gaussian projection | $O\left(\frac{s}{\epsilon}\right)$ | $O(n^{2}c)$
Subsampled randomized Hadamard transform | $O\left(\left(s+\log n\right)\left(\epsilon^{-1}+\log s\right)\right)$ | $O(n^{2}\log c)$
CountSketch | $O\left(\frac{s}{\epsilon}+s^{2}\right)$ | $O(n^{2})$
Table 1: Sufficient sketch sizes for popular sketching methods to construct
the proposed Nyström approximation $\left(\widetilde{\bm{K}}\right)_{s}$ (Wang
et al. 2019, Lemma 10). Here, $\epsilon\in(0,1)$ and
$\mu\in\left[1,\frac{n}{s}\right]$ denote the error parameter and the row
coherence of $V_{s}$, respectively, where ${V}_{s}{\Sigma}_{s}{V}_{s}^{\top}$
is the truncated SVD of $\bm{K}_{s}$.

Because small singular values in $W$ often result in numerical instability and
thus large errors in $W^{+}$, we mitigate this problem by adopting a commonly
used regularization heuristic. Specifically, for a _regularization parameter_
$l<m$, we compute the SVD of $W$ as $U_{W}\Lambda_{W}U_{W}^{\top}$ and let
$W_{l}=U_{W,l}\Lambda_{W,l}U_{W,l}^{\top}$ be the truncated SVD of $W$, where
$\operatorname{diag}(\sigma_{1},...,\sigma_{l})=\Lambda_{W,l}\in\mathbb{R}^{l\times
l}$. Then, we use

$\displaystyle\widetilde{\bm{K}}_{l}=CW_{l}^{-1}C=CU_{W,l}\Lambda_{W,l}^{-1}U_{W,l}^{\top}C^{\top}$
(10)

as a regularized version of $\widetilde{\bm{K}}$.

Finally, we compute the rank-restricted Nyström approximation with a _target
rank_ $s$ as

$\displaystyle(\widetilde{\bm{K}}_{l})_{s}=DD^{\top}\quad\text{where}\quad
D=CU_{W,l}\Lambda_{W,l}^{-1/2}\tilde{V}_{R,s}\equiv
R\tilde{V}_{R,s}\in\mathbb{R}^{n\times s},$ (11)

where $\tilde{V}_{R,s}$ is the dominant $s$ right singular vectors of
$R=CU_{W,l}\Lambda_{W,l}^{-1/2}$, $s<l$.

Because the rank-restricted Nyström approximation offers an effective way to
construct the $s$ most relevant balancing constraints, the approach computes
the weights faster while also yielding a more stable estimator. With uniform
sampling, the time and space complexities to compute
$(\widetilde{\bm{K}}_{l})_{s}$ are $O(m^{3}+nml)$ and $O(nm+m^{2})$,
respectively. Since $m$, a proxy for our computational resources, increases at
much slower rates than $n$ in typical settings, both are now near-linear in
$n$. Hence, (11) approximates the eigenvector-based kernel bases with major
computing gains. This method is also much more computationally efficient in
large sample sizes than the rank-$s$ randomized SVD, whose time complexity is
$O(n^{2}s+s^{3})$ (Li et al. 2014).

###  3.2 Scalable balancing weights using OSQP

Having described a practical way of computing the kernel representation, we
next address optimization of the weights at scale. For large-scale convex
optimization, the _alternating direction method of multipliers_ (ADMM) is
likely the most common solution method due to its computational efficiency
(Boyd et al. 2011). However, ADMM has some limitations, such as its inability
to detect infeasibility, the sensitivity of its performance to parameter
selection and data setting, and its diminished speed advantage on high-
accuracy settings. Recently, Stellato et al. (2020) developed the operator
splitting solver for quadratic programs (OSQP), a state-of-the-art ADMM-based
solver for general QPs that addresses these problems.

The OSQP algorithm is particularly well-suited to solving for the stable
balancing weights in (4). First, the most time-intensive step of OSQP is
solving the linear system of equations (Line 9 in Algorithm 1); however, since
the coefficient matrix of the linear system is symmetric quasi-definite and
sparse (with sparsity $\approx\frac{1}{n_{c}}$), a number of efficient
algorithms can be applied to this step, such as the QDLDL factorization (Davis
2005). Second, the coefficient matrix of the linear system is always non-
singular regardless of the step-size parameter values, preserving numerical
stability without sacrificing speed. Finally, when careful tuning of $\delta$
is desired (for example, by using the bootstrap (Chattopadhyay et al. 2020;
Wang and Zubizarreta 2020)), OSQP can execute this step much more efficiently
through factorization caching and warm starting. See Appendix B for more
details on ADMM and OSQP.

###  3.3 Implementation

This section describes an algorithm to formulate the proposed estimator
incorporating kernel-based stable balancing weights. Assume we have obtained
the rank-$s$ Nyström approximation $(\widetilde{\bm{K}}_{l})_{s}=DD^{\top}$ as
in (11), and for the matrix $D$, let $D_{t}\in\mathbb{R}^{n_{t}\times s}$ and
$D_{c}\in\mathbb{R}^{n_{c}\times s}$ be the rows of $D$ corresponding to the
treated and control units, respectively. Given these inputs, the proposed
kernel-based stable balancing weights are given by the optimal solution to the
following quadratic program:

$\displaystyle\underset{\begin{subarray}{c}w\in\mathbb{R}^{n_{c}}\\\ z\in
R^{1+n_{c}+l}\end{subarray}}{\text{minimize}}\quad
w^{\top}w-\frac{1}{n_{c}}w^{\top}\mathbbm{1}_{n_{c}}$ (12)
$\displaystyle\text{subject to}\quad Qw=z,z\in\mathcal{C}(\delta),$

for
$Q=\begin{bmatrix}\mathbbm{1}_{n_{c}}&\mathbb{I}_{n_{c}}&D_{c}\end{bmatrix}^{\top}$
and a set $\mathcal{C}(\delta)=\left\\{z\mid l(\delta)\leq z\leq
u(\delta)\right\\}$ where

$\displaystyle
l(\delta)=\begin{bmatrix}1&\mathbb{0}_{n_{c}}&\overline{D_{t}}-\delta\end{bmatrix}^{\top},\quad
u(\delta)=\begin{bmatrix}1&\mathbbm{1}_{n_{c}}&\overline{D_{t}}+\delta\end{bmatrix}^{\top},\quad\delta=\left(\delta_{1},...,\delta_{s}\right).$

Here, $\mathbb{I}_{r}\in\mathbb{R}^{r\times r}$ is the $r\times r$ identity
matrix, and $\mathbbm{1}_{r},\mathbb{0}_{r}\in\mathbb{R}^{1\times r}$, the
vectors of ones and zeros of length $r$, respectively; $\delta$ is a vector of
maximum allowable imbalances specified by the investigator, and
$\overline{D_{t}}\in\mathbb{R}^{1\times s}$ is a row vector containing the
mean of each column of $D_{t}$. Program (12) can be efficiently solved via
OSQP. Algorithm 1 describes computation of the proposed estimator in detail.
Note that the kernel balance constraints can also be used alongside
traditional balance requirements for low-dimensional summaries of the
covariates, such as their means.

Program (12) expresses the stable balancing weights problem (4) in the
standard nomenclature of ADMM, with the approximated kernel basis $D_{ib}$ as
$B_{b}(X_{i})$. Computation of $D$ is informed by the scaling parameters
$m,l,s$. The available computational capacity should be considered when
selecting the sketch size $m$. Choices for the regularization parameter $l$
and the target rank $s$ are related to the relative-error bound analysis in
the next section.

1 input: Sample $Z_{1},...,Z_{n}$, vector $\delta$, integers $s<l<m\ll n$, and
kernel $K(\cdot,\cdot)$

2 Sample a set $\mathcal{I}_{m}$ of $m\ll n$ indices in $\\{1,...,n\\}$

3 Form $C\in\mathbb{R}^{n\times m}$ with $C_{ij}=K(X_{i},X_{j})$,
$i\in\\{1,...,n\\}$, $j\in\mathcal{I}_{m}$, and $W\in\mathbb{R}^{m\times m}$
with $W_{ij}=K(X_{i},X_{j})$, $i,j\in\mathcal{I}_{m}$

4 Compute the rank-$l$ truncated SVD of $W$ as
$W_{l}=U_{W,l}\Lambda_{W,l}U_{W,l}^{\top}$

5 Compute the rank-$s$ truncated SVD of $R=CU_{W,l}\Lambda_{W,l}^{-1/2}$ as
$\tilde{U}_{R,s}\tilde{\Lambda}_{R,s}\tilde{V}_{R,s}^{\top}$, then compute
$D=R\tilde{V}_{R,s}$

6 Compute $Q,l(\delta),u(\delta)$ using (12)

7 Choose parameters $\rho>0$, $\sigma>0$, and $\upalpha\in(0,2)$; initialize
$w^{0}$, $y^{0}$, $z^{0}$, and $\nu^{0}$; and set $k=0$

8 repeat

9 $(\widetilde{w}^{k+1},{\nu}^{k+1})\leftarrow$ solve
$\begin{bmatrix}-(1+\sigma)\mathbb{I}_{n_{c}}&-Q^{\top}\\\
-Q&\rho^{-1}\mathbb{I}_{1+n_{c}+s}\end{bmatrix}\begin{bmatrix}\widetilde{w}^{k+1}\\\
\nu^{k+1}\end{bmatrix}=-\begin{bmatrix}\sigma
w^{k}-1/{n_{c}}\mathbbm{1}_{n_{c}}\\\ z^{k}-\rho^{-1}y^{k}\end{bmatrix}$

10 $\widetilde{z}^{k+1}\leftarrow z^{k}+\rho^{-1}(\nu^{k+1}-y^{k})$

11 $w^{k+1}\leftarrow\upalpha\widetilde{w}^{k+1}+(1-\upalpha)w^{k}$

12
$z^{k+1}\leftarrow\Pi_{\mathcal{C}(\delta)}\left(\upalpha\widetilde{z}^{k+1}+(1-\upalpha)z^{k}+\rho^{-1}y^{k}\right)$

13 $y^{k+1}\leftarrow
y^{k}+\rho\left(\upalpha\widetilde{z}^{k+1}+(1-\upalpha)z^{k}-z^{k+1}\right)$

14 $w\leftarrow w_{k}$, $k\leftarrow k+1$

15until  _termination criterion satisfied_ ;

output: $\sum_{i=1}^{n_{c}}w_{i}Y_{c,i}$, where $Y_{c}$ is an outcome vector
for the control units

Algorithm 1 Fast stable kernel balancing weights

###  3.4 Analysis of the worst-case bias bound

In the next theorem, we analyze a worst-case bound on the bias,
$\text{Bias}(\widehat{\psi};\bm{K})$, in (8):

$\displaystyle\underset{\alpha}{\sup}\left|\left(\widehat{w}^{\top}U_{c}-\frac{1}{n_{t}}\mathbbm{1}_{n_{t}}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha\right|,$

where $\widehat{w}$ is the proposed kernel-based estimator obtained via (12).

######  Theorem 3.1.

Suppose that Assumptions (2.1) and (2.2) hold. Further assume that
$\|\widetilde{\bm{K}}-\widetilde{\bm{K}}_{l}\|_{2}\leq\frac{\sigma_{l}-\sigma_{l+1}}{2}$
and $\sup\|\alpha\|_{2}<\infty$. Let $\epsilon\in(0,1)$ be an error parameter,
and $W_{l}$ and $\bm{K}_{s}$ be the best low-rank approximations for $W^{+}$
and $\bm{K}$ with ranks $l$ and $s$, respectively. Then, for all the column
sampling schemes in Table 1, with high probability at least $0.9$,

$\displaystyle\underset{\alpha}{\sup}\left|\left(\widehat{w}^{\top}U_{c}-\frac{1}{n_{t}}\mathbbm{1}_{n_{t}}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha\right|$
$\displaystyle=(1+\epsilon)\left\|\bm{K}-\bm{K}_{s}\right\|_{*}$ (13)
$\displaystyle\quad+O\left(\|W^{+}-W^{-1}_{l}\|_{2}+\|W^{+}-W^{-1}_{l}\|_{F}^{2}+\delta\right),$

where $\bm{K}_{s}$ is the best rank-$s$ approximation of $\bm{K}$.
$\|\cdot\|_{2}$, $\|\cdot\|_{*}$ and $\|\cdot\|_{F}$ are the 2-norm, the trace
norm, and the Frobenius norm, respectively.

See Appendix A in Supplementary Materials for the formal definitions of the
matrix norms and the proof. The above theorem provides an upper bound on the
bias and dissects it into three components. The first component,
$(1+\epsilon)\left\|\bm{K}-\bm{K}_{s}\right\|_{*}$, pertains to the bias due
to the degree of kernel matrix approximation. It provides guarantees on
relative-error trace norm approximation of the rank-$s$ Nyström approximation
without regularization. While the magnitude of this bias component may be
reduced by choosing a larger value for $s$, this decision must be informed by
choices for the other scaling parameters $l$ and $m$ as well as our error
tolerance (see Remark 3). The second component of the bias,
$\|W^{+}-W^{-1}_{l}\|_{2}+\|W^{+}-W^{-1}_{l}\|_{F}^{2}$, represents the
additive error due to the regularization step described in (10), and declines
to zero with no regularization. Typically, we expect the second component to
be much smaller than the first. The third and final component is the bias from
the residual covariate imbalance after weighting, which is directly controlled
by the tolerance vector $\delta$ governing the balancing constraints.

######  Remark 3.

[Scaling parameter selection] The upper bound in (3.13) depends on sketching
schemes, as well as the scaling parameters $m,l,s$, error parameter
$\epsilon$, and matrix coherence $\mu\in\left[1,\frac{n}{s}\right]$ of the
dominant $s$-dimensional singular space of $\bm{K}$. Although it is often not
straightforward to characterize this upper bound in real datasets, some
guidelines for selecting scaling parameters could be established. In Table 1,
we provide sufficient sketch sizes $m$ in terms of $s,\epsilon,\mu$. While a
larger $s$ better approximates the original kernel matrix, it must be some
fraction of $m$, and it is desirable to set $m$ as small as possible since it
is a proxy for available computational resources. Thus, balancing these two
sources of error requires choosing an $s$ value that is sufficiently large to
reduce error from matrix approximation, while also being a sufficiently small
fraction of $m$. It is known that when uniform sampling is used to form the
Nystrom approximation $m=\tilde{O}\left(\frac{s\mu}{\epsilon}\right)$ (Wang et
al. 2019, Lemma 10). Hence, if one sets $\log s\sim\frac{1}{\epsilon}$, then
the chosen $m$ should be on the order of $\frac{s\mu}{\epsilon}$. When the
budget of column samples $m$ must to be fixed, the target rank could be
$s\sim\frac{m}{\mu/\epsilon}$ where $m$ is at least $\mu/\epsilon>1$ times
larger than $s$. The coherence parameter $\mu$ can be estimated from the
truncated SVD of $\bm{K}$ (see, e.g., Mohri and Talwalkar 2011). The
regularization parameter $l$ should remain as close to $m$ as possible as long
as $\sigma_{l}$ is not extremely small, guaranteeing numerical stability of
$\Lambda_{W,l}^{-1}$. Unfortunately, there is no theory for the choice of $l$
as it is related to the finite precision error. Following Wang et al. (2019),
we apply a simple heuristic that $l=\lceil\frac{s+m}{2}\rceil$. Future
theoretical and empirical work is required for optimal scaling parameter
selection, as well as characterizing how large the upper bound could be in
real datasets.

##  4 Simulation study

We study the computational performance and estimation accuracy of several
modeling and balancing weighting methods. We extend the simulation design by
Hainmueller (2012) to encompass nonlinear treatment assignment processes and
sample sizes up to a million.

###  4.1 Study design

Consider an observational study with continuous and binary covariates
generated as follows

$\displaystyle\begin{bmatrix}X_{1}\\\ X_{2}\\\ X_{3}\end{bmatrix}\sim
N\begin{pmatrix}\begin{bmatrix}0\\\ 0\\\
0\end{bmatrix},\begin{bmatrix}2&1&-1\\\ 1&1&-0.5\\\
-1&-0.5&1\end{bmatrix}\end{pmatrix},\quad X_{4}\sim\text{Unif}[-3,3],\quad
X_{5}\sim\chi^{2}_{1},\quad X_{6}\sim\text{Bern}[0.5].$

Using these six covariates, we generate the outcome and treatment variables as

$\displaystyle Y=(X_{1}+X_{2}+X_{5})^{2}+\eta,\quad\eta\sim N(0,1),$ (14)
$\displaystyle
A=\mathbbm{1}\left\\{X_{1}^{2}+2X_{2}^{2}-2X_{3}^{2}-(X_{4}+1)^{3}-0.5\log(X_{5}+10)+X_{6}-1.5+\varepsilon>0\right\\},$
(15)

where $\varepsilon\sim N(0,\sigma^{2}_{\varepsilon})$. We consider two
settings with $\sigma^{2}_{\varepsilon}=30$ and
$\sigma^{2}_{\varepsilon}=100$, hereafter referred to as the weak overlap and
strong overlap settings, respectively. We define untransformed covariates
$X=(X_{1},X_{2},X_{3},X_{4},X_{5},X_{6})$ and transformed covariates
$X^{*}=(X_{1}X_{3},X_{2}^{2},X_{4},X_{5},X_{6})$. For each
$\sigma^{2}_{\varepsilon}$, we consider each interaction of treatment
assignment mechanism and accuracy of covariate specification across five
sample sizes: $n\in\\{2\cdot 10^{3},5\cdot 10^{3},10^{4},10^{5},10^{6}\\}$. In
this design, the ATT equals zero by construction. To separate design from
analysis and focus on scalability and speed, we based the weighting estimator
(2) solely on $\widehat{\pi}$; however, it can be augmented by an outcome
model to form a doubly-robust estimator as in Robins and Rotnitzky (1995).

Estimation of the ATT by weighting requires explicit modeling of $\pi$ for
$\widehat{\psi}_{\text{mod}}$ in the modeling approach, or solving the QP (12)
to obtain $\widehat{\psi}_{\text{bal}}$ in the balancing approach. We consider
various methods, including some considered useful for large-scale modeling or
optimization. For practicality, we focus on methods in readily available R
packages. For consistent comparisons of balancing estimators, we fix
$\delta=0.0005$ across all methods. However, we also evaluate performance when
selecting $\delta$ in a data-driven fashion using the algorithm in
Chattopadhyay et al. (2020).

We consider nine methods for $\widehat{\psi}_{\text{mod}}$ and seven methods
for $\widehat{\psi}_{\text{bal}}$, for a total of 16 methods classified into
six groups (see Table 2. For $\widehat{\psi}_{\text{mod}}$, Group M1 methods
use logistic regression (or GLM) to model the propensity score. Group M2
contains lasso algorithms with the tuning parameter selected via 10-fold cross
validation. In Groups M1 and M2, the regression models for $\pi$ add quadratic
terms and all pairwise products of observed covariates. Group M3 includes
three common nonparametric algorithms: random forests, kernel ridge
regression, and generalized additive models using splines. For
$\widehat{\psi}_{\text{bal}}$, we solve (12) where we use the Gaussian kernel
and set $m=300$, $l=200$, $s=100$ for the Nyström approximation procedure.
Group B1 includes open-source QP solvers in R. Group B2 studies two commercial
solvers. Finally, Group B3 includes three state-of-the-art ADMM-based solvers
aimed at reducing processing time and used with the kernel balancing approach.
For OSQP, we use the same parameters as Stellato et al. (2020).

We assess estimator performance using the root-mean-squared error (RMSE)
defined by
$\left\\{\frac{1}{S}\sum_{j=1}^{S}(\widehat{\psi}_{j}-\psi)^{2}\right\\}^{1/2}$
averaged across $S=500$ simulations, and quantify the computational cost in
average central processing unit (CPU) time across simulations, which we report
in seconds unless otherwise noted. All simulations were completed using the
Harvard Faculty of Arts and Sciences Research Computing clusters with
available memory fixed at 128GB. In the results tables, we write ‘F’ when the
optimization process failed in more than half the simulations, and ‘M’ when
the system ran out of memory in more than half the simulations.

Table 2: Description of modeling and solution methods used in the simulation study.  | Method | R package | Description
---|---|---|---
Modeling Approach |  M1. Logistic regression | glm |  | Commonly used method to fit logistic
---
regression models in R
sgdglm |  | Stochastic gradient descent methods for
---
estimation with big data (Tran et al. 2015)
M2. Regularized regression | glmnet |  | Lasso via coordinate descent
---
(Friedman et al. 2010)
biglasso |  | Scalable lasso for big data
---
(Zeng and Breheny 2017)
oem |  | Lasso for tall (large n) data
---
(Xiong et al. 2016)
penreg |  | ADMM-based lasso implementation
---
(Boyd et al. 2011) (admm.lasso in results)
M3. Non-parametric regression | ranger | Fast implementation for Random Forests
DRR |  | Fast implementation for Kernel Ridge
---
Regression (kernel.ridge in results)
| bam |  | Generalized additive model using splines
---
for big data (Wood et al. 2015)
Balancing Approach | B1. Open-source solvers | quadprog |  | Commonly-used QP solver based on the
---
dual method (Goldfarb and Idnani 1983)
qpoases |  | QP solver based on parametric
---
active-set algorithm (Ferreau et al. 2014)
B2. Commercial solvers | gurobi |  | Large scale commercial optimizer
---
by Gurobi Optimization
Rmosek |  | Large scale commercial optimizer
---
by Mosek ApS
B3. ADMM-based solvers | osqp |  | Operator splitting solver for QP (OSQP)
---
(Stellato et al. 2020)
pogs |  | ADMM-based graph form solver (POGS)
---
(Fougner and Boyd 2018)
scs |  | Operator splitting conic solver (SCS)
---
(O’donoghue et al. 2016)

###  4.2 Empirical results

Table 3: Root-mean square error (RMSE) under weak overlap. Each cell contains the RMSE for each combination of method, sample size, and covariate specification, averaged across 500 simulations. “M" signifies the system ran out of memory in more than half the simulations, while “F" signifies that the optimization process failed to obtain a weighting solution in more than half the simulations. |  |  $X$ correctly specified ($X$) |  $X$ transformed ($X^{*}$)
---|---|---|---
Method $n$ | 2k | 5k | 10k | 100k | 1M | 2k | 5k | 10k | 100k | 1M
$\widehat{\psi}_{\text{mod}}$ | glm | 1.72 | 1.29 | 1.20 | 1.19 | 1.18 | 1.85 | 1.37 | 1.30 | 1.18 | 1.18
sgdglm | 1.72 | 1.29 | 1.20 | 1.19 | 1.18 | 1.86 | 1.37 | 1.31 | 1.18 | 1.19
glmnet | 2.20 | 1.18 | 0.68 | 0.23 | 0.12 | 4.57 | 1.79 | 0.67 | 0.39 | 0.35
biglasso | 2.57 | 1.08 | 0.65 | 0.22 | 0.12 | 5.23 | 1.68 | 0.68 | 0.39 | 0.35
oem | 3.37 | 1.27 | 0.79 | 0.26 | 0.12 | 6.47 | 2.97 | 0.79 | 0.40 | 0.35
admm.lasso | 3.83 | 1.59 | 0.84 | 0.29 | 0.15 | 7.29 | 3.04 | 1.07 | 0.70 | 0.51
ranger | 0.71 | 0.33 | 0.22 | 0.15 | 0.09 | 0.97 | 0.41 | 0.36 | 0.27 | 0.22
kernel.ridge | 1.87 | 1.19 | 1.01 | 0.89 | M | 6.26 | 6.07 | 5.59 | 4.87 | M
bam | 2.84 | 1.31 | 1.20 | 1.19 | 1.18 | 3.37 | 1.95 | 1.25 | 1.19 | 1.18
$\widehat{\psi}_{\text{bal}}$ | quadprog | 0.31 | 0.20 | 0.15 | F | F | 0.51 | 0.34 | 0.33 | F | F
qpoases | 0.32 | 0.22 | 0.15 | F | F | 0.65 | 0.34 | 0.24 | F | F
gurobi | 0.19 | 0.13 | 0.09 | 0.04 | 0.02 | 0.31 | 0.19 | 0.15 | 0.08 | 0.05
Rmosek | 0.19 | 0.13 | 0.10 | 0.04 | 0.02 | 0.31 | 0.20 | 0.15 | 0.08 | 0.05
osqp | 0.19 | 0.03 | 0.10 | 0.04 | 0.02 | 0.31 | 0.19 | 0.16 | 0.08 | 0.05
pogs | 0.31 | 0.20 | 0.15 | M | M | 0.59 | 0.36 | 0.29 | M | M
scs | 0.25 | 0.22 | 0.19 | 0.15 | 0.10 | 0.55 | 0.38 | 0.25 | 0.18 | 0.15

Table 3 shows estimation accuracy under weak overlap. The proposed kernel
balancing estimators $\widehat{\psi}_{\text{bal}}$ substantially outperform
the modeling estimators $\widehat{\psi}_{\text{mod}}$ at all considered sample
sizes, especially in small to moderate ones, regardless of covariate
specification; however, gaps between $\widehat{\psi}_{\text{bal}}$ and
$\widehat{\psi}_{\text{mod}}$ are smaller with transformed covariates $X^{*}$.
Performance varied greatly within the class of modeling approaches; in
particular, random forests were most successful at larger sample sizes. We do
not report results for kernel-based stable balancing weights without the
Nyström approximation, since no methods were successful for samples of 100,000
or more.

Because the average CPU time did not appreciably vary across simulation
settings, we present only their averages in Table 4. For
$\widehat{\psi}_{\text{bal}}$, forming the kernel basis requires less
computation time than solving the QP, demonstrating the efficiency of the
Nyström approximation in the kernel balancing algorithm. When solving the same
QP (12), OSQP was faster than the commercial solvers MOSEK and Gurobi, by
about an order of magnitude in the case of the latter. The GLM and OSQP
methods were fastest, taking only seconds to find the optimal weights even in
the largest sample size of one million. In contrast, the lasso and non-
parametric methods required much more time at that sample size, in some cases
hundreds of times more than GLM and OSQP. Although not presented here, the
OSQP and GLM-based estimators took under a minute on average to handle a
simulated dataset of 10 million observations, whereas the slowest methods took
up to days.

Table 4: Average CPU time across simulation settings. Each cell contains the CPU time for each combination of method, sample size, and covariate specification, averaged across 500 simulations. All values are reported in seconds unless specified as being in hours (“hr"). “M" signifies the system ran out of memory in more than half the simulations, while “F" signifies that the optimization process failed to obtain a weighting solution in more than half the simulations. |  | Time (in seconds unless otherwise noted)
---|---|---
Method $n$ | 2k | 5k | 10k | 100k | 1M
$\widehat{\psi}_{\text{mod}}$ | glm | <0.1 | <0.1 | <0.1 | 0.4 | 8.6
sgdglm | <0.1 | <0.1 | <0.1 | 0.3 | 5.0
glmnet | 1.0 | 1.8 | 4.4 | 18 | 362
biglasso | 9.2 | 12 | 20 | 266 | 0.8hr
oem | 38 | 81 | 224 | 0.6hr | 5hr
admm.lasso | 1.3 | 2.7 | 6.3 | 15 | 311
ranger | 8.7 | 30 | 128 | 0.4hr | 3hr
kernel.ridge | 1.1 | 3.6 | 29 | 0.3hr | M
bam | 0.7 | 1.6 | 3.8 | 21 | 210
$\widehat{\psi}_{\text{bal}}$ | quadprog | 0.9 | 15 | 119 | F | F
qpoases | 48 | 127 | 392 | F | F
gurobi | <0.1 | <0.1 | 0.5 | 6.3 | 94
Rmosek | <0.1 | <0.1 | <0.1 | 4.2 | 37
osqp | <0.1 | <0.1 | <0.1 | 3.2 | 9.6
pogs | 3.7 | 46 | 305 | M | M
scs | <0.1 | 5.3 | 14.5 | 27 | 169

Several methods faced difficulties at larger sample sizes, particularly
$\widehat{\psi}_{\text{bal}}$ approaches. At samples $\geq$ 100,000, the open-
source solvers failed to converge, while the system ran out of memory when
attempting to execute POGS, perhaps because it cannot accommodate sparse
matrices. Among $\widehat{\psi}_{\text{mod}}$ approaches, kernel ridge
regression ran out of memory at 1 million observations.

Results for RMSE, CPU time, convergence, and memory exhaustion in the strong
covariate overlap setting were generally consistent with those in weak
overlap, but with smaller performance gaps between
$\widehat{\psi}_{\text{mod}}$ and $\widehat{\psi}_{\text{bal}}$ (see Appendix
C for details).

As noted in Section 3.2, the OSQP algorithm facilitates factorization caching
and warm starting for even better computational efficiency when multiple
iterations are required, such as for data-driven selection of $\delta$. To
demonstrate, we repeatedly solved a sequence of OSQP-based weighting problems
in the same simulation settings using $100$ different values of $\delta$. The
results showed that the OSQP algorithm achieved an additional $2.4$-fold
improvement in CPU time on average, whereas the commercial solvers did not
exhibit such efficiency gains (see Appendix D).

##  5 A national study of heart attack care by hospital profit status

Heart attacks are generally treated with either interventional cardiology
through percutaneous coronary intervention (PCI), or alternatively medical
management. In the more severe subtype known as ST-elevation MIs, PCI is
generally recommended. However, older patients are much more likely to have
non-ST-elevation MIs, and guidelines are less clear on PCI use in this subtype
and more generally among older patients with complex comorbid conditions
(Lawton et al. 2022). Thus, clinical decision-making in this setting can be
affected by both individual physician practice style and the norms of the
institutions where they work. Medicare reimburses hospitals more when patients
receive interventional cardiology, possibly incentivizing for-profit hospitals
to do so more often. We asked whether for-profit hospitals were more likely to
use interventional cardiology on heart attack patients than other hospitals,
and whether patients outcomes. differed Two previous studies reported opposing
results on relationships between profit status, PCI use, and outcomes. Sloan
et al. (2003) found for-profit hospitals had higher rates of PCI than other
hospitals yet similar mortality rates. In contrast, Shah et al. (2007) found
similar rates of PCI use among non-ST elevation MI cases. However, these
studies had fewer than 160,000 patients and used data from small numbers of
hospitals in voluntary reporting systems; for example, Shah et al. (2007)
studied only 58 for-profit hospitals.

To build on previous work by applying our proposed methods to a much larger
and more representative sample, we studied a national dataset of Medicare
administrative claims from the Centers for Medicare & Medicaid Services (CMS)
for patients admitted to short-term acute-care hospitals for initial treatment
of heart attack. We studied 1,273,636 Medicare beneficiaries hospitalized for
heart attack between October 1, 2010 and November 30, 2019. We defined
treatment as admission to a for-profit hospital, and control as admission to
other types of hospitals.

A fair comparison of for-profit and control hospitals requires careful
adjustment for differences in potentially confounding covariates. For example,
heart attack treatment and outcomes can be affected by patient demographics
and contemporary clinical guidelines. Thus, we adjusted for patient age on
admission, sex, race, and admission year. The type of heart attack also
strongly affects treatment and prognosis, as ST-elevation Myocardial
Infarctions (STEMIs) are much more likely to receive PCI but are also much
more severe than Non-ST-Elevation Myocardial Infarctions (NSTEMIs). Thus, we
defined these by checking the inpatient record for corresponding International
Classification of Diseases, 9th Revision, Clinical Modifictation (ICD-9-CM) or
ICD-10-CM principal diagnosis codes. We then included covariates for these
types of heart attack in the adjustment model. We also adjusted for emergent
admission type and history of PCI or coronary artery bypass grafting (CABG)
surgery. As comorbid risk factors such as history of diabetes, renal disease,
or stroke can affect both prognosis and likelihood of receiving PCI, we also
checked the inpatient records for evidence of a set of 43 comorbidities drawn
from risk models maintained by CMS, concentrating on conditions noted as
present on admission (POA). This screening only counted diagnoses explicitly
documented as already existing upon the patient’s admission to the hospital
rather than those noted as having occurred after the patient was admitted.

Patients may choose where to live in part on local healthcare resources, which
can affect both selection into hospitals and the outcome. Thus, we also
defined county-level measures of hospital characteristics such as academic
affiliation, nurse staffing levels, size in number of beds, available cardiac
technologies, and urban location, for a total of 17 hospital characteristics.
This partially addresses selection bias due to textcolorredresidential choice
based on potential medical need.

While the simulations included 6 adjustment covariates, our case study
includes 70 patient covariates and 17 area-level hospital characteristics for
a total of 87 covariates. Given the additional complexity arising from this
much number of covariates and the size of the dataset (over 1 million control
patients and 1.27 million in total), we proceeded with the stable kernel
balancing approach as it showed excellent computational efficiency without
diminished estimator accuracy.

To make the assumption of conditional exchangeability more plausible in this
demonstration of the method, we were careful to leverage the Medicare data as
described above to adjust for as many observable pre-treatment covariates as
possible such as demographics, heart attack type, risk factors, and hospital
market characteristics. However, as this is an observational study, it is
possible for there to be omitted unobserved confounders.

We used the proposed weighting methods to balance the 87 covariates and
functions thereof, in order to estimate the average treatment effect on the
treated; i.e., the effect of for-profit admission on the likelihood of
receiving PCI. We assessed balance using target absolute standardized mean
differences (TASMD) (Chattopadhyay et al. 2020), considering the for-profit
and control groups balanced if no TASMD exceeded 0.1. Because the same set of
weights can in principle be used to examine multiple endpoints, we used the
weights to assess PCI use and several outcomes: mortality within 30 days of
admission, length of stay, and readmission within 30 days of discharge.

###  5.1 Kernel balance

Table 5 shows selected characteristics of 206,948 patients from 775 for-profit
and 1,066,780 patients from 2,763 control hospitals. See Appendix E for a
balance table showing all 87 covariates.

Table 5: Selected Covariates and Sample Size before Weighting. The table shows selected unweighted covariate means and sample sizes among all patients admitted to for-profit hospitals or other hospitals, along with the target absolute standardized mean difference (TASMD) before weighting. Covariate (% unless noted) / Sample Size | For-Profit | Controls | TASMD Before
---|---|---|---
_Demographic Characteristics_ |  |  |
Age at admission (years, mean) | 79.0 | 79.3 | 0.03
Female | 47.8% | 47.7% | 0.00
Non-Hispanic White | 79.0% | 84.2% | 0.13
Black | 8.3% | 7.8% | 0.02
Hispanic | 8.6% | 4.4% | 0.15
_Comorbidities Present on Admission_ |  |  |
Diabetes | 37.6% | 37.1% | 0.01
Heart Failure | 43.5% | 43.8% | 0.01
Valvular Disease | 15.6% | 19.1% | 0.09
_County of Residence Hospital Characteristics_ |  |  |
Small Size ($<$ 251 beds) | 79.5% | 72.9% | 0.44
Nonteaching | 87.6% | 81.5% | 0.41
_Sample Size_ |  |  |
Patients | 206,948 | 1,066,670 |
Hospitals | 775 | 2,763 |

Because heart attack is an acute condition, the time and place of its
occurrence is not easily predictable and patients are usually brought to the
closest hospital; hence, many characteristics such as age and several
comorbidities are fairly well-balanced, although a number of imbalances
persist and require adjustment. While the occurrence of heart attack in
relation to the type of nearby hospital may be somewhat haphazard, the
distribution of patient race as well as the characteristics of the hospitals
in their counties are not. In particular, patients admitted to for-profit
hospitals were more likely to be nonwhite, and more likely to reside in
counties where hospitals tended to be smaller and nonteaching. Given the
policy importance of health equity, these community differences emphasize the
need to benchmark for-profit hospital performance.

We used the same kernel and scaling parameters as in the simulations. However,
to prioritize estimator accuracy in the case study, we used the high-accuracy
setting with a much stricter termination criterion and solution polishing
(Stellato et al. 2020, Section 4), and set $\delta=0$ rather than
$\delta=0.0005$. These made the QP more difficult to solve, but favored
estimator accuracy.

The stable kernel balancing algorithm efficiently produced excellent balance.
Figure 1 shows that despite large initial imbalances, the weighting algorithm
successfully balanced all covariates with all TASMDs after weighting well
below 0.1, and none exceeding 0.03. In our dataset of 1.27 million patients,
the kernel basis was computed in just 32.3 seconds, while the OSQP algorithm
solved for the stable kernel balancing weights in 132 seconds.

Figure 1: Each dot represents the Target Absolute Standardized Mean Difference
(TASMD) of one of 87 covariates between for-profit and control hospitals
before and after weighting.

###  5.2 Heart attack treatment and outcomes

Table 6 shows treatment and outcomes by hospital profit status after
weighting. We found that for-profit hospitals used PCI at almost the same rate
as control hospitals, which aligns with the more recent work by Shah et al.
(2007), but conflicts with the earlier study by Sloan et al. (2003).
Differences between our study and previous ones may also arise because we used
a national sample of hospitals rather than sets of hospitals in voluntary
reporting programs.

Despite similar PCI use, for-profit hospitals had significantly higher 30-day
mortality relative to control hospitals, a finding that differs from the older
study by Sloan et al. (2003). Compared to more contemporary work by Horwitz et
al. (2017), we found a larger difference in mortality; however, their analysis
pooled patients hospitalized for various reasons. Length of stay was
significantly but not meaningfully longer at for-profit hospitals. However,
readmission rates were significantly and meaningfully higher at for-profit
hospitals.

Table 6: PCI Use and Outcomes by Hospital Profit Status. The table shows outcome means among for-profit hospital patients, weighted outcome means among control hospital patients, and a 95% confidence interval for the difference. Measure (percent unless noted) | For-Profit | Control | Difference | 95% CI
---|---|---|---|---
Received PCI | 33.2 | 33.1 | 0.1 | (-0.1, 0.3)
30-day Mortality | 14.3 | 13.9 | 0.4 | (0.2, 0.6)
Length of Stay (days) | 5.76 | 5.63 | 0.13 | (0.10, 0.15)
30-day Readmission/IP Death | 19.2 | 18.2 | 0.9 | (0.7, 1.1)
30-day Readmission Only | 11.5 | 10.6 | 0.8 | (0.6, 0.9)

##  6 Summary

In this paper, we describe a weighting approach that overcomes previous
barriers to the use of the balancing approach at large scales. We show that
the rank-restricted Nyström approximation makes it computationally feasible to
find a set of general functions of the covariates in an RKHS. We describe how
the kernel basis from the Nyström approximation can be efficiently balanced
via a modern ADMM-based convex optimization method. Our extensive simulation
study shows that the proposed estimator has superior computation speed
relative to other balancing approaches and even logistic regression-based
modeling approaches. Moreover, this efficiency gain does not diminish the
proposed estimator’s accuracy, as it performed as well as or better than other
approaches in most simulation settings. Hence, our study provides new insights
into the performance of weighting estimators in large datasets.

We used the proposed approach in a large observational study of heart attack
treatment and outcomes in for-profit and control hospitals in a national
Medicare dataset. We achieved excellent balance in under 3 minutes. After
adjustment, for-profit hospitals used PCI at similar rates as others, but had
significantly higher readmission rates. Future work could examine the impact
of Medicare’s Hospital Readmissions Reduction Program on for-profit hospitals.

##  7 Software

The proposed algorithm can be implemented in R with the code provided in
https://github.com/kwangho-joshua-kim/osqp-kernel-sbw.

## References

  * Austin and Stuart (2015) Austin, P. C. and Stuart, E. A. (2015), “Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies,” Statistics in Medicine, 34, 3661–3679.
  * Ben-Michael et al. (2021) Ben-Michael, E., Feller, A., Hirshberg, D. A., and Zubizarreta, J. R. (2021), “The balancing act in causal inference,” arXiv preprint arXiv:2110.14831.
  * Boyd et al. (2011) Boyd, S., Parikh, N., and Chu, E. (2011), Distributed optimization and statistical learning via the alternating direction method of multipliers, Now Publishers, Inc.
  * Chattopadhyay et al. (2020) Chattopadhyay, A., Hase, C. H., and Zubizarreta, J. R. (2020), “Balancing vs modeling approaches to weighting in practice,” Statistics in Medicine, 39, 3227–3254.
  * Davis (2005) Davis, T. A. (2005), “Algorithm 849: A concise sparse Cholesky factorization package,” ACM Transactions on Mathematical Software (TOMS), 31, 587–591.
  * Ferreau et al. (2014) Ferreau, H. J., Kirches, C., Potschka, A., Bock, H. G., and Diehl, M. (2014), “qpOASES: A parametric active-set algorithm for quadratic programming,” Mathematical Programming Computation, 6, 327–363.
  * Fougner and Boyd (2018) Fougner, C. and Boyd, S. (2018), “Parameter selection and preconditioning for a graph form solver,” in Emerging Applications of Control and Systems Theory, Springer, pp. 41–61.
  * Friedman et al. (2010) Friedman, J., Hastie, T., and Tibshirani, R. (2010), “Regularization paths for generalized linear models via coordinate descent,” Journal of Statistical Software, 33, 1–22.
  * Gittens and Mahoney (2016) Gittens, A. and Mahoney, M. W. (2016), “Revisiting the Nyström method for improved large-scale machine learning,” The Journal of Machine Learning Research, 17, 3977–4041.
  * Goldfarb and Idnani (1983) Goldfarb, D. and Idnani, A. (1983), “A numerically stable dual method for solving strictly convex quadratic programs,” Mathematical Programming, 27, 1–33.
  * Hainmueller (2012) Hainmueller, J. (2012), “Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies,” Political Analysis, 20, 25–46.
  * Halko et al. (2011) Halko, N., Martinsson, P.-G., and Tropp, J. A. (2011), “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM review, 53, 217–288.
  * Hazlett (2020) Hazlett, C. (2020), “KERNEL BALANCING,” Statistica Sinica, 30, 1155–1189.
  * Hirshberg et al. (2021) Hirshberg, D. A., Maleki, A., and Zubizarreta, J. R. (2021), “Minimax Linear Estimation of the Retargeted Mean,” .
  * Hirshberg and Wager (2020) Hirshberg, D. A. and Wager, S. (2020), “Augmented Minimax Linear Estimation,” .
  * Horwitz et al. (2017) Horwitz, L. I., Bernheim, S. M., Ross, J. S., Herrin, J., Grady, J., Krumholz, H. M., Drye, E. E., and Lin, Z. (2017), “Hospital characteristics associated with risk-standardized readmission rates,” Medical Care, 55, 528–534.
  * Kallus (2020) Kallus, N. (2020), “Generalized Optimal Matching Methods for Causal Inference,” Journal of Machine Learning Research, 21, 1–54.
  * Kang and Schafer (2007) Kang, J. D. Y. and Schafer, J. L. (2007), “Demystifying double robustness: a comparison of alternative strategies for estimating a population mean from incomplete data (with discussion),” Statistical Science, 22, 523–539.
  * Kern et al. (2016) Kern, H. L., Stuart, E. A., Hill, J., and Green, D. P. (2016), “Assessing methods for generalizing experimental impact estimates to target populations,” Journal of Research on Educational Effectiveness, 9, 103–127.
  * Kumar et al. (2012) Kumar, S., Mohri, M., and Talwalkar, A. (2012), “Sampling methods for the Nyström method,” The Journal of Machine Learning Research, 13, 981–1006.
  * Lawton et al. (2022) Lawton, J. S., Tamis-Holland, J. E., Zwischenberger, B. A., et al. (2022), “2021 ACC/AHA/SCAI Guideline for Coronary Artery Revascularization,” Journal of the American College of Cardiology, 79, e21–e129.
  * Lehoucq et al. (1998) Lehoucq, R. B., Sorensen, D. C., and Yang, C. (1998), ARPACK users’ guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods, SIAM.
  * Li et al. (2014) Li, M., Bi, W., Kwok, J. T., and Lu, B.-L. (2014), “Large-scale Nyström kernel matrix approximation using randomized SVD,” IEEE transactions on neural networks and learning systems, 26, 152–164.
  * Little and Rubin (2019) Little, R. J. and Rubin, D. B. (2019), Single imputation methods, Wiley Online Library.
  * Micchelli et al. (2006) Micchelli, C. A., Xu, Y., and Zhang, H. (2006), “Universal Kernels.” Journal of Machine Learning Research, 7, 2651–2667.
  * Mohri and Talwalkar (2011) Mohri, M. and Talwalkar, A. (2011), “Can matrix coherence be efficiently and accurately estimated?” in Proceedings of the fourteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, pp. 534–542.
  * O’donoghue et al. (2016) O’donoghue, B., Chu, E., Parikh, N., and Boyd, S. (2016), “Conic optimization via operator splitting and homogeneous self-dual embedding,” Journal of Optimization Theory and Applications, 169, 1042–1068.
  * Pourkamali-Anaraki et al. (2018) Pourkamali-Anaraki, F., Becker, S., and Wakin, M. (2018), “Randomized clustered nystrom for large-scale kernel machines,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 32.
  * Robins and Rotnitzky (1995) Robins, J. M. and Rotnitzky, A. (1995), “Semiparametric efficiency in multivariate regression models with missing data,” Journal of the American Statistical Association, 90, 122–129.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983), “The central role of the propensity score in observational studies for causal effects,” Biometrika, 70, 41–55.
  * Rubin (1974) Rubin, D. B. (1974), “Estimating causal effects of treatments in randomized and nonrandomized studies.” Journal of Educational Psychology, 66, 688.
  * Rubin (2006) — (2006), “The design versus the analysis of observational studies for causal effects: Parallels with the design of randomized trials,” Statistics in Medicine, 26, 20–36.
  * Shah et al. (2007) Shah, B. R., Glickman, S. W., Liang, L., Gibler, W. B., Ohman, E. M., Pollack, C. V., Roe, M. T., and Peterson, E. D. (2007), “The impact of for-profit hospital status on the care and outcomes of patients with non–ST-segment elevation myocardial infarction: results from the CRUSADE Initiative,” Journal of the American College of Cardiology, 50, 1462–1468.
  * Sloan et al. (2003) Sloan, F. A., Trogdon, J. G., Curtis, L. H., and Schulman, K. A. (2003), “Does the ownership of the admitting hospital make a difference? Outcomes and process of care of Medicare beneficiaries admitted with acute myocardial infarction,” Medical Care, 41, 1193–1205.
  * Stellato et al. (2020) Stellato, B., Banjac, G., Goulart, P., Bemporad, A., and Boyd, S. (2020), “OSQP: an operator splitting solver for quadratic programs,” Mathematical Programming Computation, 12, 637–672.
  * Tran et al. (2015) Tran, D., Toulis, P., and Airoldi, E. M. (2015), “Stochastic gradient descent methods for estimation with large data sets,” arXiv preprint arXiv:1509.06459.
  * Vu et al. (2021) Vu, T., Chunikhina, E., and Raich, R. (2021), “Perturbation expansions and error bounds for the truncated singular value decomposition,” Linear Algebra and its Applications, 627, 94–139.
  * Wahba (1990) Wahba, G. (1990), Spline models for observational data, SIAM.
  * Wang et al. (2019) Wang, S., Gittens, A., and Mahoney, M. W. (2019), “Scalable kernel K-means clustering with Nyström approximation: relative-error bounds,” The Journal of Machine Learning Research, 20, 431–479.
  * Wang and Zubizarreta (2020) Wang, Y. and Zubizarreta, J. R. (2020), “Minimal dispersion approximately balancing weights: asymptotic properties and practical considerations,” Biometrika, 107, 93–105.
  * Wong and Chan (2018) Wong, R. K. and Chan, K. C. G. (2018), “Kernel-based covariate functional balancing for observational studies,” Biometrika, 105, 199–213.
  * Wood et al. (2015) Wood, S. N., Goude, Y., and Shaw, S. (2015), “Generalized additive models for large data sets,” Journal of the Royal Statistical Society: Series C, 64, 139–155.
  * Xiong et al. (2016) Xiong, S., Dai, B., Huling, J., and Qian, P. Z. (2016), “Orthogonalizing EM: A design-based least squares algorithm,” Technometrics, 58, 285–293.
  * Zeng and Breheny (2017) Zeng, Y. and Breheny, P. (2017), “The biglasso package: A memory-and computation-efficient solver for lasso model fitting with big data in R,” arXiv preprint arXiv:1701.05936.
  * Zhao (2019) Zhao, Q. (2019), “Covariate balancing propensity score by tailored loss functions,” Annals of Statistics, 47, 965 – 993.
  * Zhou et al. (2017) Zhou, L., Pan, S., Wang, J., and Vasilakos, A. V. (2017), “Machine learning on big data: Opportunities and challenges,” Neurocomputing, 237, 350–361.
  * Zubizarreta (2015) Zubizarreta, J. R. (2015), “Stable weights that balance covariates for estimation with incomplete outcome data,” Journal of the American Statistical Association, 110, 910–922.

Supplementary Material for “Scalable kernel balancing weights in a nationwide
observational study of hospital profit status and heart attack outcomes" by
Kwangho Kim, Bijan A. Niknam, and José R. Zubizarreta

## Appendix A. Proof of Theorem 1

We use three matrix norms for our result: for a matrix $A$ with the
$(i,j)$-entry $a_{ij}$, we define

Frobenius Norm:
$\displaystyle\|A\|_{F}=\sqrt{\sum_{i,j}a_{ij}^{2}}=\sqrt{\sum_{i}\sigma_{i}^{2}(A)}$
Spectral Norm:
$\displaystyle\|A\|_{2}=\underset{\|x\|_{2}=1}{\max}\|Ax\|_{2}=\sigma_{1}(A)$
Trace Norm: $\displaystyle\|A\|_{*}=\sum_{i}\sigma_{i}(A).$

###### Proof.

Recall that from Section 3.1,

$\displaystyle\widetilde{\bm{K}}=CW^{+}C^{\top},$
$\displaystyle\widetilde{\bm{K}}_{l}=CW^{-1}_{l}C^{\top},$
$\displaystyle(\widetilde{\bm{K}}_{l})_{s}=DD^{\top}\quad\text{where}\quad
D=CU_{W,l}\Lambda_{W,l}^{-1/2}\tilde{V}_{R,s}.$

Let $w_{t}=\frac{1}{n_{t}}\mathbb{1}_{n_{t}}$. Then by the Cauchy–Schwarz and
triangle inequalities, we have that

$\displaystyle\left|\left(\widehat{w}_{c}^{\top}U_{c}-w_{t}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha\right|$
$\displaystyle\leq\left|\left(\widehat{w}_{c}^{\top}U_{c}-w_{t}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha-\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha+\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha\right|$
$\displaystyle\leq\|\widehat{w}_{c}\|_{2}\left\|U_{c}\Lambda
U^{\top}-D_{c}D^{\top}\right\|_{2}\|\alpha\|_{2}+\|w_{t}\|_{2}\left\|U_{t}\Lambda
U^{\top}-D_{t}D^{\top}\right\|_{2}\|\alpha\|_{2}$
$\displaystyle\quad+\left|\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha\right|$
$\displaystyle\leq
2\left\|\bm{K}-(\widetilde{\bm{K}}_{l})_{s}\right\|_{2}\|\alpha\|_{2}+\left|\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha\right|$
$\displaystyle\leq
2\left\\{\left\|\bm{K}-(\widetilde{\bm{K}})_{s}\right\|_{2}+\left\|(\widetilde{\bm{K}})_{s}-(\widetilde{\bm{K}}_{l})_{s}\right\|_{2}\right\\}\|\alpha\|_{2}+\left|\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha\right|.$

When
$\|\widetilde{\bm{K}}-\widetilde{\bm{K}}_{l}\|_{2}\leq\frac{\sigma_{l}-\sigma_{l+1}}{2}$,
by Vu et al. (2021, Theorem 1), we have the following first-order
approximation

$\displaystyle\|(\widetilde{\bm{K}})_{s}-(\widetilde{\bm{K}}_{l})_{s}\|_{2}$
$\displaystyle\lesssim\|\widetilde{\bm{K}}-\widetilde{\bm{K}}_{l}\|_{2}+\sum_{i=1}^{l}\sum_{j=l+1}^{n}\left(\frac{\sigma_{j}^{2}}{\sigma_{i}^{2}-\sigma_{j}^{2}}+\frac{\sigma_{i}\sigma_{j}}{\sigma_{i}^{2}-\sigma_{j}^{2}}\right)\|\widetilde{\bm{K}}-\widetilde{\bm{K}}_{l}\|_{2}+\|\widetilde{\bm{K}}-\widetilde{\bm{K}}_{l}\|_{F}^{2}$
$\displaystyle\lesssim\|W^{+}-W^{-1}_{l}\|_{2}+\|W^{+}-W^{-1}_{l}\|_{F}^{2}.$

Also, by Theorem 1 of Wang et al. (2019), we have

$\displaystyle\left\|\bm{K}-(\widetilde{\bm{K}})_{s}\right\|_{2}$
$\displaystyle\leq\left\|\bm{K}-(\widetilde{\bm{K}})_{s}\right\|_{*}$
$\displaystyle\leq(1+\epsilon)\left\|\bm{K}-\bm{K}_{s}\right\|_{*}$

with high probability at least $0.9$.

For the last term in the last display, by the Cauchy–Schwarz inequality it
follows that

$\displaystyle\underset{\alpha}{\sup}\left|\left(\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right)D^{\top}\alpha\right|$
$\displaystyle\leq\left\|\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right\|_{2}\underset{\alpha}{\sup}\left\|D^{\top}\alpha\right\|_{2}$
$\displaystyle\lesssim\left\|\widehat{w}_{c}^{\top}D_{c}-w_{t}^{\top}D_{t}\right\|_{1}$
$\displaystyle\lesssim\delta.$

Hence, provided that $\underset{\alpha}{\sup}\|\alpha\|_{2}<\infty$, putting
the two pieces together we have the desired probabilistic bound

$\displaystyle\underset{\alpha}{\sup}\left|\left(\widehat{w}_{c}^{\top}U_{c}-w_{t}^{\top}U_{t}\right)\Lambda
U^{\top}\alpha\right|$
$\displaystyle=(1+\epsilon)\left\|\bm{K}-\bm{K}_{s}\right\|_{*}+O\left(\|W^{+}-W^{-1}_{l}\|_{2}+\|W^{+}-W^{-1}_{l}\|_{F}^{2}+\delta\right).$

∎

## Appendix B. More details on ADMM and OSQP

Convex QPs such as (4) are well-studied and various solution methods exist.
Some commonly used solution methods include active-set, interior-point, and
operator splitting methods. Active-set methods are the traditional algorithms
for solving QPs. These algorithms explore the feasible region containing the
solution by iteratively adding and dropping active constraints until the
solution is found. Interior-point methods, the default algorithms used by many
modern commercial solvers, solve unconstrained optimization problems for
varying barrier functions until the optimum is achieved.

Typically, active-set methods outperform interior-point methods when the
number of covariates or constraints is small, or when warm starts are used.
Otherwise, interior-point methods generally outperform active-set methods in
terms of scalability and speed. That said, in general the iteration costs of
interior-point methods grow non-linearly with the dimensionality of the
problem, rendering them impractical for very large-scale settings.
Alternatively, operator splitting methods obtain an optimal solution using
only first-order information about the cost function (i.e., the subgradients).
They reduce the optimization problem to that of finding zeros of the sum of
monotone operators. Here, we focus on a particular class of operator splitting
methods known as the alternating direction method of multipliers (ADMM) (Boyd
et al. 2011)) which has been shown to extend very well to large-scale convex
optimization problems.

ADMM. The _Alternating direction method of multipliers_ is likely the most
commonly used operator splitting method due to its simplicity and ability to
accommodate large problems. ADMM solves problems of the form

$\displaystyle{\text{minimize}}\quad g(\tilde{v})+h(v)$ (16)
$\displaystyle\text{subject to}\quad G\tilde{v}+Hv=q$

for two sets of variables $\tilde{v}$ and $v$. Following Stellato et al.
(2020), we introduce auxiliary variables $\tilde{w}$ and $\tilde{z}$ and write
Problem (12) as

$\displaystyle{\text{minimize}}\quad\tilde{w}^{\top}\mathbb{I}_{n_{c}}\tilde{w}-\frac{1}{n_{c}}\tilde{w}^{\top}\mathbb{1}_{n_{c}}+\mathcal{I}_{Qw=z}(\tilde{w},\tilde{z})+\mathcal{I}_{\mathcal{C}(\delta)}(z)$
(17) $\displaystyle\text{subject
to}\quad\mathbb{I}_{1+2n_{c}+Kd}[\tilde{w},\tilde{z}]^{\top}-\mathbb{I}_{1+2n_{c}+Kd}[w,z]^{\top}=0,$

where $\mathcal{I}_{S}$ is the indicator function for the set $S$. Problem
(17) is in the ADMM form given by (16), where
$g(\tilde{w},\tilde{z})=\tilde{w}^{\top}\mathbb{I}_{n_{c}}\tilde{w}-\frac{1}{n_{c}}\tilde{w}^{\top}\mathbb{1}_{n_{c}}+\mathcal{I}_{Qw=z}(\tilde{w},\tilde{z})$
and $h(w,z)=\mathcal{I}_{\mathcal{C}(\delta)}(z)$. Problem (17) can be solved
via the following ADMM iterations

$\displaystyle(w^{k+1},\tilde{z}^{k+1})\leftarrow\underset{Qw=z}{\operatorname*{arg\,min}}\quad{w}^{\top}\mathbb{I}_{n_{c}}{w}-\frac{1}{n_{c}}{w}^{\top}\mathbb{1}_{n_{c}}+\frac{\sigma}{2}\left\|x-w^{k}\right\|^{2}+\frac{\rho}{2}\left\|z-z^{k}+\frac{y^{k}}{\rho}\right\|^{2}$
(18) $\displaystyle
z^{k+1}\leftarrow\Pi_{\mathcal{C}(\delta)}\left(\tilde{z}^{k+1}+\frac{y^{k}}{\rho}\right)$
(19) $\displaystyle y^{k+1}\leftarrow
y^{k}+\rho\left(\tilde{z}^{k+1}-z^{k+1}\right)$ (20)

where $y$ is a dual covariate, $\sigma,\rho>0$ are step-size parameters, and
$\Pi_{\mathcal{C}(\delta)}$ denotes the Euclidean projection onto
$\mathcal{C}(\delta)$ (see Boyd et al. 2011 for details).

ADMM blends the advantages of dual decomposition and the method of multipliers
(or the augmented Lagrangian method). First, as a method of multipliers, ADMM
exhibits excellent convergence properties. In fact, Problem (16) does not
require any assumptions for convergence beyond convexity of $g$ and $h$ and
that its augmented Lagrangian has a saddle point, an advantage not shared by
other operator splitting methods. Second, as described by the steps in
equations (18) and (19), in ADMM $w^{k}$ and $z^{k}$ are updated sequentially
or alternatingly, where each step is typically far less expansive than the
original minimization problem jointly over $(w,z)$. Indeed, the inner QP (18)
is an equality-constrained problem which can be reduced to solving a linear
system, and (19) is simply a box projection. This decomposability therefore
substantially reduces the computational burden required for execution, and
allows ADMM to find solutions much faster than conventional methods.

Despite these appealing features, ADMM has some limitations. One potential
limitation is its inability to detect infeasibility. Moreover, it loses its
speed advantage when aiming to find high-accuracy solutions, and its
convergence rates also highly depend on the parameter selection and the data
setting. Despite these challenges, significant progress has recently been made
in adapting ADMM to very large-scale convex optimization problems. O’donoghue
et al. (2016) and Fougner and Boyd (2018) proposed the splitting-based conic
solver (SCS) and proximal operator graph-form solver (POGS), respectively.
Based on their own parameter selection heuristics and data preconditioning,
they showed that both SCS and POGS scale well to extremely large problems,
with SCS also able to detect infeasibility. More recently, Stellato et al.
(2020) developed the operator splitting solver for quadratic programs (OSQP),
a state-of-the-art ADMM-based solver particularly tailored to general QPs.
OSQP is the first generic QP solver that addresses all three challenges of
detecting infeasibility, finding high-accuracy solutions, and adapting to
different data settings.

input : $Q,l(\delta),u(\delta)$

parameter : $\rho>0,\sigma>0,\upalpha\in(0,2)$

1 Initialize $w^{0},y^{0},z^{0},\nu^{0}$ and set $k=0$

2 repeat

3 $(\widetilde{w}^{k+1},{\nu}^{k+1})\leftarrow$ solve
$\begin{bmatrix}-(1+\sigma)\mathbb{I}_{n_{c}}&-Q^{\top}\\\
-Q&\rho^{-1}\mathbb{I}\end{bmatrix}\begin{bmatrix}\widetilde{w}^{k+1}\\\
\nu^{k+1}\end{bmatrix}=-\begin{bmatrix}\sigma
w^{k}-1/{n_{c}}\mathbbm{1}_{n_{c}}\\\ z^{k}-\rho^{-1}y^{k}\end{bmatrix}$ ;

4 $\widetilde{z}^{k+1}\leftarrow z^{k}+\rho^{-1}(\nu^{k+1}-y^{k})$

5 $w^{k+1}\leftarrow\upalpha\widetilde{w}^{k+1}+(1-\upalpha)w^{k}$

6
$z^{k+1}\leftarrow\Pi_{\mathcal{C}(\delta)}\left(\upalpha\widetilde{z}^{k+1}+(1-\upalpha)z^{k}+\rho^{-1}y^{k}\right)$

7 $y^{k+1}\leftarrow
y^{k}+\rho\left(\upalpha\widetilde{z}^{k+1}+(1-\upalpha)z^{k}-z^{k+1}\right)$

8 $k\leftarrow k+1$

9until  _termination criterion satisfied_ ;

Algorithm 2 OSQP for SBW

OSQP for SBW. For convenience, the OSQP algorithm for SBW is presented in
Algorithm 2 separately from Figure 1. Steps 3 and 4 correspond to the reduced
Karush-Kuhn-Tucker (KKT) system for the inner QP (18), and Steps 6 and 7
correspond to (19) and (20), respectively. In Steps 5, 6, and 7, the
relaxation method is used for the _w_ -, _z_ -, and _y_ -updates with the
relaxation parameter $\upalpha$. The termination criterion is satisfied when
both the primal and dual residuals are smaller than some prespecified
tolerance, or if infeasibility is detected with some prespecified level of
accuracy. Although the method of choosing the optimal parameters
$(\rho,\sigma,\upalpha)$ may vary, Stellato et al. (2020) provides guidance
for doing so with OSQP. The parameter $\sigma$ behaves like a regularization
term and can be as small as possible, so long as it ensures that a unique
solution of the linear system in Step 3 of Algorithm 2 exists. Also, through
extensive empirical research, it has been shown that setting the relaxation
parameter $\upalpha$ in the range $[1.5,1.8]$ generally improves the
convergence rate. To set the step size $\rho$, OSQP adopts heuristics that
assign higher values for active constraints and lower values for inactive
constraints. Furthermore, to ensure fast convergence over a wide range of
problems, the algorithm uses the ratio between the primal and dual residuals
to adaptively update $\rho$ from one iteration to the next. Stellato et al.
(2020) showed that the proposed parameter selection heuristics in OSQP exhibit
solid improvements in terms of convergence and stability.

We remark that POGS and SCS are alternative ADMM-based methods that could be
applied; however, these approaches require reformulation of our QP, which
introduces many extra coefficients into the optimization procedure and
therefore may not be computationally efficient (Stellato et al. 2020).

## Appendix C. RMSE under strong overlap

Results for under strong overlap are presented in Web Table 7.

|  |  $X$ correctly specified ($X$) |  $X$ transformed ($X^{*}$)
---|---|---|---
Method $n$ | 2k | 5k | 10k | 100k | 1M | 2k | 5k | 10k | 100k | 1M
$\widehat{\psi}_{\text{mod}}$ | glm | 1.13 | 0.72 | 0.48 | 0.40 | 0.40 | 1.31 | 0.80 | 0.56 | 0.41 | 0.40
sgdglm | 1.13 | 0.72 | 0.48 | 0.40 | 0.40 | 1.31 | 0.80 | 0.56 | 0.41 | 0.40
glmnet | 1.10 | 0.87 | 0.58 | 0.17 | 0.07 | 1.30 | 1.23 | 0.51 | 0.21 | 0.13
biglasso | 1.25 | 0.82 | 0.55 | 0.16 | 0.08 | 1.30 | 1.02 | 0.50 | 0.21 | 0.12
oem | 1.23 | 0.82 | 0.58 | 0.17 | 0.08 | 1.33 | 1.49 | 0.71 | 0.22 | 0.13
admm.lasso | 1.78 | 0.87 | 0.58 | 0.22 | 0.14 | 2.23 | 1.12 | 0.56 | 0.27 | 0.16
ranger | 1.22 | 0.86 | 0.66 | 0.22 | 0.09 | 1.27 | 0.85 | 0.52 | 0.20 | 0.07
kernel.ridge | 1.18 | 0.84 | 0.67 | 0.29 | M | 1.53 | 0.91 | 0.67 | 0.32 | M
bam | 2.47 | 0.67 | 0.46 | 0.41 | 0.40 | 3.37 | 1.95 | 1.25 | 1.19 | 1.18
$\widehat{\psi}_{\text{bal}}$ | quadprog | 0.07 | 0.04 | 0.03 | F | F | 0.41 | 0.24 | 0.18 | F | F
qpoases | 0.07 | 0.05 | F | F | F | 0.41 | 0.28 | F | F | F
gurobi | 0.05 | 0.03 | 0.02 | 0.01 | 0.01 | 0.38 | 0.21 | 0.19 | 0.09 | 0.03
Rmosek | 0.05 | 0.03 | 0.02 | 0.01 | 0.01 | 0.38 | 0.21 | 0.19 | 0.09 | 0.03
osqp | 0.05 | 0.03 | 0.02 | 0.01 | 0.01 | 0.39 | 0.21 | 0.19 | 0.10 | 0.03
pogs | 0.7 | 0.04 | 0.03 | M | M | 0.40 | 0.25 | 0.25 | M | M
scs | 0.11 | 0.10 | 0.09 | 0.05 | 0.02 | 0.39 | 0.25 | 0.20 | 0.13 | 0.10
Table 7: RMSE under weak overlap

## Appendix D. Time reduction via factorization caching and warm starting

Figure 2: Average cumulative CPU time over multiple iterations of SBW calls

Since the vectors $l,u$ are the only parameters that vary in $\delta$ in
Algorithm 2, we can speed up repeated OSQP calls via factorization caching and
warm starting. Specifically, although the tuning procedure requires a number
of iterations, because the coefficient matrix in the relevant linear system is
fixed across iterations and does not depend on $\delta$, factorization caching
can be used, storing the factorization from the first iteration for use in
subsequent ones. Also, because small changes in $\delta$ do not substantially
affect weighting solutions, warm starting can be used for added efficiency. To
illustrate this computational gain, we solve a sequence of OSQP-based SBW
problems using $100$ values of $\delta$ equally spaced between $0.001$ and
$0.1$ on the same data, and compare its computing time against two popular
commercial solvers gurobi and Rmosek (with the default setting of interior-
point methods). Data are generated by the same simulation setup used in
Section 4 under weak overlap. We measure the cumulative time for each solver
at $m=1,20,40,60,80,100$ problem instances across different values of
$\delta$. We repeat the simulation $100$ times and present their average
values in Web Figure 2.

For gurobi and Rmosek, the cumulative time increases linearly with the number
of iterations as would be expected; i.e., $m$ iterations took roughly $m$
times as long as a single execution. However, $m$ iterations of the OSQP
algorithm took much less time than $m\times\text{(the execution time for a
single run)}$; in fact, it is on average $2.4$ times faster. This $2.4$-fold
time improvement demonstrates the benefits of factorization caching and warm
starting via OSQP when implementing SBW.

## Appendix E. Complete case study balance table before and after stable
kernel balancing

Covariate |  For-Profit |  Unweighted Controls |  Weighted Controls |  TASMD Before |  TASMD After
---|---|---|---|---|---
Age at Admission |  79.043 |  79.322 |  79.035 |  0.034 |  0.001
Age 66-74 |  0.369 |  0.360 |  0.367 |  0.019 |  0.005
Age 75-84 |  0.371 |  0.365 |  0.371 |  0.012 |  0.002
Age 85+ |  0.260 |  0.275 |  0.262 |  0.034 |  0.004
Admission Year |  2014.771 |  2014.668 |  2014.778 |  0.039 |  0.002
Admission 2010 |  0.029 |  0.030 |  0.029 |  0.008 |  0.000
Admission 2011 |  0.110 |  0.118 |  0.109 |  0.027 |  0.002
Admission 2012 |  0.112 |  0.117 |  0.112 |  0.016 |  0.001
Admission 2013 |  0.107 |  0.112 |  0.108 |  0.016 |  0.001
Admission 2014 |  0.109 |  0.108 |  0.108 |  0.004 |  0.004
Admission 2015 |  0.111 |  0.108 |  0.111 |  0.010 |  0.000
Admission 2016 |  0.113 |  0.110 |  0.113 |  0.009 |  0.001
Admission 2017 |  0.114 |  0.108 |  0.116 |  0.019 |  0.003
Admission 2018 |  0.099 |  0.098 |  0.100 |  0.003 |  0.003
Admission 2019 |  0.094 |  0.089 |  0.094 |  0.018 |  0.002
Sex (female) |  0.478 |  0.477 |  0.477 |  0.003 |  0.002
Non-Hispanic White |  0.790 |  0.842 |  0.794 |  0.128 |  0.008
Black |  0.083 |  0.078 |  0.082 |  0.019 |  0.004
Hispanic |  0.086 |  0.044 |  0.085 |  0.148 |  0.003
Asian/Pacific Islander |  0.025 |  0.019 |  0.025 |  0.039 |  0.001
Native American |  0.005 |  0.005 |  0.004 |  0.007 |  0.019
Other Race |  0.011 |  0.012 |  0.011 |  0.013 |  0.001
Emergent Admission |  0.837 |  0.838 |  0.836 |  0.002 |  0.003
Anterior ST-elevation MI 1 |  0.069 |  0.077 |  0.069 |  0.032 |  0.002
Other ST-elevation MI |  0.149 |  0.156 |  0.150 |  0.018 |  0.002
History of PCI |  0.150 |  0.159 |  0.152 |  0.024 |  0.004
History of CABG |  0.119 |  0.119 |  0.119 |  0.001 |  0.001
Hx Sepsis |  0.015 |  0.014 |  0.015 |  0.005 |  0.001
Hx Infection |  0.002 |  0.002 |  0.002 |  0.001 |  0.001
Hx Metastatic Cancer |  0.011 |  0.013 |  0.011 |  0.024 |  0.001
Hx Metastatic or Severe Cancer |  0.022 |  0.025 |  0.022 |  0.024 |  0.001
Hx Lung or Severe Cancer |  0.014 |  0.016 |  0.014 |  0.015 |  0.001
Hx Non-metastatic Cancer |  0.029 |  0.034 |  0.030 |  0.027 |  0.003
Hx Cancer, Other |  0.003 |  0.004 |  0.002 |  0.009 |  0.015
Hx Diabetes |  0.376 |  0.371 |  0.375 |  0.010 |  0.002
Hx Diabetes w/ Complications |  0.129 |  0.134 |  0.130 |  0.015 |  0.002
Hx Diabetes w/o Complications |  0.083 |  0.072 |  0.083 |  0.040 |  0.002
Hx Malnutrition |  0.042 |  0.036 |  0.041 |  0.031 |  0.005
Hx Morbid Obesity |  0.677 |  0.695 |  0.679 |  0.037 |  0.003
Hx Liver Disease |  0.016 |  0.017 |  0.016 |  0.006 |  0.000
Hx Dementia |  0.139 |  0.131 |  0.138 |  0.022 |  0.002
Hx Substance Use |  0.117 |  0.113 |  0.117 |  0.012 |  0.000
Hx Major Psychiatric Disorder |  0.018 |  0.017 |  0.019 |  0.010 |  0.001
Hx Cardio-respiratory Failure or Shock |  0.160 |  0.159 |  0.159 |  0.002 |  0.003
Hx Heart Failure |  0.435 |  0.438 |  0.436 |  0.006 |  0.002
Hx AMI |  0.006 |  0.007 |  0.006 |  0.007 |  0.001
Hx Unstable Angina |  0.065 |  0.062 |  0.065 |  0.011 |  0.000
Hx Valvular Disease |  0.156 |  0.191 |  0.159 |  0.095 |  0.006
Hx Hypertensive Heart Disease |  0.021 |  0.016 |  0.021 |  0.035 |  0.000
Hx Hypertension |  0.476 |  0.477 |  0.476 |  0.002 |  0.000
Hx Arrhythmia, Specified |  0.285 |  0.289 |  0.286 |  0.008 |  0.003
Hx Arrhythmia, Other |  0.129 |  0.152 |  0.130 |  0.070 |  0.005
Hx Stroke |  0.007 |  0.007 |  0.005 |  0.003 |  0.022
Hx Cerebrovascular Disease |  0.040 |  0.045 |  0.039 |  0.026 |  0.002
Hx Precerebral Occlusion |  0.032 |  0.036 |  0.032 |  0.025 |  0.000
Hx Peripheral Vascular Disease |  0.058 |  0.062 |  0.059 |  0.014 |  0.001
Hx Lung Fibrosis |  0.017 |  0.019 |  0.017 |  0.019 |  0.001
Hx Asthma |  0.028 |  0.033 |  0.029 |  0.027 |  0.002
Hx COPD |  0.215 |  0.202 |  0.213 |  0.032 |  0.004
Hx Pneumonia |  0.100 |  0.088 |  0.099 |  0.039 |  0.002
Hx Pneumothorax |  0.019 |  0.018 |  0.019 |  0.007 |  0.001
Hx Other Respiratory Disease |  0.064 |  0.070 |  0.065 |  0.023 |  0.003
Hx Eye, Other |  0.010 |  0.012 |  0.011 |  0.015 |  0.003
Hx Renal Failure |  0.364 |  0.359 |  0.364 |  0.009 |  0.000
Hx UTI |  0.089 |  0.078 |  0.088 |  0.039 |  0.004
Hx Male Genital Disorder |  0.087 |  0.092 |  0.088 |  0.018 |  0.004
Hx Skin Ulcer |  0.022 |  0.024 |  0.022 |  0.016 |  0.002
Hx Skin Disease, Other |  0.010 |  0.013 |  0.010 |  0.027 |  0.000
Hx Iatrogenic Complication |  0.035 |  0.037 |  0.035 |  0.010 |  0.001
Hx Trauma |  0.023 |  0.024 |  0.024 |  0.006 |  0.003
_County Hospital Characteristics_ |  |  |  |  |
Urban Hospital |  0.881 |  0.877 |  0.883 |  0.010 |  0.006
Large Size (501+ Beds) |  0.052 |  0.077 |  0.052 |  0.338 |  0.007
Medium Size (251-100 Beds) |  0.153 |  0.194 |  0.154 |  0.305 |  0.006
Small Size (Under 251 Beds) |  0.795 |  0.729 |  0.793 |  0.442 |  0.009
Resident-to-Bed Ratio |  0.021 |  0.034 |  0.021 |  0.336 |  0.002
Major Teaching (RTB $\geq$0.25) |  0.026 |  0.039 |  0.026 |  0.195 |  0.006
Minor Teaching (0 $\leq$RTB $\leq$0.25) |  0.098 |  0.146 |  0.098 |  0.390 |  0.001
Nonteaching (RTB=0) |  0.876 |  0.815 |  0.876 |  0.407 |  0.001
Nurse-to-Bed Ratio |  0.971 |  1.063 |  0.977 |  0.120 |  0.008
High NTB (NTB$\geq$1.0) |  0.376 |  0.421 |  0.377 |  0.214 |  0.002
Moderate NTB (0.9 $\leq$NTB $<$1.0) |  0.053 |  0.052 |  0.053 |  0.009 |  0.005
Low NTB (NTB$<$0.9) |  0.571 |  0.527 |  0.571 |  0.208 |  0.000
Cardiac Technology Index |  0.921 |  1.022 |  0.921 |  0.263 |  0.000
Cardiac Surgery Provided |  0.246 |  0.280 |  0.246 |  0.234 |  0.004
Coronary Care Unit |  0.411 |  0.480 |  0.412 |  0.330 |  0.005
Cardiac Catheterization Lab |  0.265 |  0.262 |  0.263 |  0.012 |  0.009
All Three Cardiac Services Available |  0.140 |  0.154 |  0.137 |  0.112 |  0.022