Multivariate Tie-breaker Designs

There are many important settings where a costly or scarce intervention can be made for some but not all subjects. Examples include giving a scholarship to some students (Angrist et al., 2014), intervening to prevent child abuse in some families (Krantz, 2022), programs to counter juvenile delinquency (Lipsey et al., 1981) and sending some university students to remedial English classes (Aiken et al., 1998). There are also lower-stakes settings where a company might be able to offer a perk such as a free service upgrade, to some but not all of its customers. In settings such as these, there is usually a priority ordering of subjects, perhaps based on how deserving they are or on how much they or the investigator might gain from the intervention. We can represent that priority order in terms of a real-valued running variable $x_{i}$ for subject $i$.

To maximize the immediate short-term value from the limited intervention, one can assign it only to subjects with $x_{i}\geqslant t$ for some threshold $t$. The difficulty with this “greedy” solution is that such an allocation makes it difficult to estimate the causal effect of the treatment. It is possible to use a regression discontinuity design (RDD) in this case, comparing subjects with $x_{i}$ somewhat larger than $t$ to those with $x_{i}$ somewhat smaller than $t$. For details on the RDD, see the comprehensive recent survey by Cattaneo and Titiunik (2022). The RDD can give a consistent nonparametric estimate of the treatment effect at $x=t$ (Hahn et al., 2001), but not at other values of $x$. The treatment effect can be studied at other values of $x$ under an assumed regression model, but then the variance of the regression coefficients can be large due to the strong dependence between the running variable and the treatment (Gelman and Imbens, 2019).

The RDD is commonly used to analyze observational data for which the investigator has no control over the treatment cutoff. In the settings we consider, the investigators assign the treatment and can therefore employ some randomization. The motivating context makes it costly or even ethically difficult to use a randomized controlled trial (RCT) where the treatment is assigned completely at random without regard to $x_{i}$. The tie-breaker design (TBD) is a compromise between the RCT and RDD. It is a triage where subjects with large values of $x_{i}$ get the treatment, those with small values get the control level and those in between are randomized to either treatment or control.

While the tie-breaker design has been known since Campbell (1969) it has not been subjected to much analysis. Most of the theory for the TBD has been for the setting with a scalar $x$ and models that are linear in $x$ and the treatment (Owen and Varian, 2020; Li and Owen, 2022). In this paper we study the TBD for a vector-valued predictor. Our first motivation is that many use cases for TBDs will include multiple covariates. Second, although multivariate nonparametric regression models are out of the scope of this paper, we believe that TBD regressions are a useful first step in that direction. Third, Gelman et al. (2020) counsel against RDDs that do not adjust for pre-treatment variables. The multivariate regression setting supports such adjustments.

In a tie-breaker design we have a vector of covariates for subject $i$ and must choose a two-level treatment variable. We then face an atypical experimental design problem where some of the predictors are fixed and not settable by us while one of them is subject to randomization. This is known as the ‘marginally restricted design’ problem after Cook and Thibodeau (1980) who studied $D$-optimality in such a setting. We consider two such settings. In one setting, the investigator already has the covariate values. The second setting is one step removed, where we consider random covariates that an investigator might later have.

The paper is organized as follows. Section 2 outlines our notation and the regression model. Section 3 introduces our notions of efficiency and D-optimality in the multivariate regression model. Theorem 1 shows that D-optimality for the treatment effect parameters is equivalent to D-optimality for the full regression. To study efficiency for future subjects, we use a prospective D-optimality criterion, adapted from Bayesian optimal design, that maximizes expected future information. Theorem 2 then shows that the RCT is prospectively D-optimal. We also discuss an example involving Gaussian covariates and a symmetric design, which provides relevant intuition. Section 4 finds an expression for the expected short-term gain, with particular focus on this Gaussian case. When the running variable is linear in the covariates, the best linear combination for statistical efficiency is the “last” eigenvector of the covariance matrix of covariates while the best linear combination for short-term gain is the true treatment effect. Section 5 presents a design strategy based on convex optimization to choose treatment probabilities for given covariates and compares the effects of various practical constraints. In particular, we show that a monotonicity constraint in the treatment probabilities yields solutions with a few distinct treatment probability levels. This is consistent with some past results in constrained experimental design (Cook and Fedorov, 1995), but both our proof and the theorem particulars are distinct. Section 6 illustrates this procedure on a hospital data set from MIMIC-IV-ED about which emergency room patients should receive intensive care. Section 7 has a brief discussion of some additional context for our results.

Given $n$ subjects with $d$ covariates each, we let $X\in\mathbb{R}^{n\times d}$ be the design matrix and $X_{i}\in\mathbb{R}^{d}$ be the variables for subject $i$. We write $\tilde{X}=\begin{bmatrix}1&X\end{bmatrix}\in\mathbb{R}^{n\times(d+1)}$ to denote the matrix with an intercept out front, and $\tilde{X}=[1\quad X^{\top}]^{\top}\in\mathbb{R}^{d+1}$ to denote its $i$’th row. For ease of notation, we zero-index $\tilde{X}$ so that $\tilde{X}_{i0}=1$ and $\tilde{X}_{ij}=X_{ij}$, for $j=1,\dots,d$.

We are interested in the effect of some treatment $Z_{i}\in\{-1,1\}$ on a future response $Y_{i}\in\mathbb{R}$ for subject $i$, where $Z_{i}=1$ corresponds to treatment and $Z_{i}=-1$ to control. The design problem is to choose probabilities $p_{i}\in[0,1]$ and then take $\mathbb{P}(Z_{i}=1)=p_{i}$. This differs from the common experimental design framework in which the covariates $X_{i}$ can also be chosen. In Section 5 we will show how to get optimal $p_{i}$ by convex optimization.

To get more general insights into the design problem, in Section 3 we also consider a random data framework. The predictors are to be sampled with $X_{i}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}P_{X}$. This allows us to relate design findings to the properties of $P_{X}$ rather than to a specific matrix $X$. After $X_{i}$ are observed, $Z_{i}$ will be set randomly and then $Y_{i}$ observed. The analyst is unable to alter $P_{X}$ at all but can choose any function $p(X)=\mathbb{P}(Z=1\!\mid\!X)$ that satisfies the imposed constraints.

We work with the following linear model:

for $\tilde{\beta},\tilde{\gamma}\in\mathbb{R}^{d+1}$, where $\varepsilon_{i}$ are IID noise terms with mean zero and variance $\sigma^{2}>0$. We use the same notational convention of writing $\tilde{\beta}=\begin{bmatrix}\beta_{0}&\beta^{\top}\end{bmatrix}^{\top}$ and $\tilde{\gamma}=\begin{bmatrix}\gamma_{0}\quad\gamma^{\top}\end{bmatrix}^{\top}$ for $\beta,\gamma\in\mathbb{R}^{d}$ to separate out the intercept term. We consider $\tilde{\gamma}$ to be the parameter of greatest interest because it captures the treatment effect of $Z$.

Equation (2) generalizes the two-line model (1) studied by Li and Owen (2023) and Owen and Varian (2020). The latter authors describe some computational methods for the model (2) but most of their theory is for model (1).

Though a strong parametric assumption, this multi-line regression is a helpful and practical model to inform treatment assignment at the design stage. Because our scenario of interest could include regions of the covariate space with $p(x)=0$, we do not have overlap and cannot rely on classical semiparametric methods for causal inference. Nonparametric generalizations of (2) such as spline models may prove more flexible in highly nonlinear settings, but we do not explore them here.

In Section 3, we also consider multivariate tie-breaker designs (TBDs), which we define as follows: the treatment $Z_{i}$ is assigned via

for parameters $u<\ell$, $p\in(0,1)$ and $\eta\in\mathbb{R}^{d}$. That is, we assign treatment to subject $i$ whenever $X_{i}^{\top}\eta$ is at or above some upper cutoff $u$, do not assign treatment whenever it is at or below some lower cutoff $\ell$, and randomize at some fixed probability $p$ in the middle. For the case $u=\ell$, we take $\mathbb{P}(Z_{i}=1\!\mid\!X_{i})=\mathds{1}\{X_{i}^{\top}\eta\geqslant u\}$ which has a mild asymmetry in offering the treatment to those subjects, if any, that have $X_{i}^{\top}\eta=\ell=u$.

In equation (3) we ignore the intercept and just consider $X_{i}$ instead of $\tilde{X}_{i}$ since the intercept would merely shift everything by a constant. Here, we use $\eta\in\mathbb{R}^{d}$ instead of $\gamma$ from (2) to reflect that the vector we treat on need not be the same as the true $\gamma$, which is unknown.

In practice, $\eta$ will encode whatever constraints on randomization a practitioner may encounter. While linear constraints do not encompass all possibilities, they are simple and flexible enough to account for a great many that one could feasibly seek to impose; for example, constraints on a heart rate below some cutoff or a total SAT score that is sufficiently high can both be covered by (3). One could also set $\ell$ and $u$ to be known or estimated quantiles of a covariate or a linear combination of them if, for example, one wants deterministic treatment assignment for the top $5\%$ of candidates. More generally, we show in Section 5 how to accommodate any convex constraints on treatment assignment.

The assignment (3) greatly generalizes the one in Owen and Varian (2020) which had $d=1$, $\ell=-u$, $p=1/2$, and $P_{X}$ either $\mathcal{U}(-1,1)$ or $\mathcal{N}(0,1)$. In analogy to the one-dimensional case, we refer to the case with $u=\ell$ as an RDD. We refer to any choice of $(\ell,u,\eta)$ for which $\mathbb{P}(X^{\top}\eta\in(\ell,u))=1$ as an RCT.

We assume that the covariates $X_{i}\in\mathbb{R}^{d}$ have yet to be observed and are drawn from some distribution $P_{X}$ with a finite and invertible covariance matrix $\Sigma$. The aim is to devise a treatment assignment scheme $p(X)$ with $Z_{i}$ assigned independently via $\mathbb{P}(Z_{i}=1\text{ }|\text{ }X_{i})=p(X_{i})$ that is optimal in some sense. Random allocation of $Z$ has advantages of fairness and supports some randomization-based inference.

Section 3.1 briefly outlines efficiency and D-optimality in the fixed $Z$ setting to introduce relevant formulas. In Section 3.2, we then outline a notion of prospective $D$-optimality used when $X$ is random and we must choose a distribution of $Z$ given $X$. We model our approach on the Bayesian optimal design literature (Chaloner and Verdinelli, 1995), in which the model parameters are treated as random, which we describe in more detail in that section.

The treatment of $X$ as random permits a theoretical discussion about what is expected to happen depending on the distributional properties of the unseen $X$ data. However, the case in which the covariates are actually fixed is also covered by the procedure in Section 3.2, since it corresponds to a point mass distribution on $X$. We return to this case in Section 5, in which we show how to use convex optimization to obtain prospectively D-optimal $p(X_{i})$.

We turn now to the other arm of the tradeoff, the short-term gain. In our motivating problems, the response is defined so that larger values of $Y_{i}$ define better outcomes. Now $\mathbb{E}[Y_{i}]=\mathbb{E}[\tilde{X}_{i}^{\top}\tilde{\beta}]+\tilde{T}_{i}$ where $\tilde{T}_{i}=\mathbb{E}[Z_{i}\tilde{X}_{i}^{\top}\tilde{\gamma}]$. The first term in $\mathbb{E}[Y_{i}]$ is not affected by the treatment allocation and so any consideration of short-term gain can be expressed in terms of $\tilde{T}_{i}$. Further, $\mathbb{E}(\tilde{T}_{i})=\mathbb{E}[Z_{i}\gamma_{0}]+T_{i}$ for $T_{i}=\mathbb{E}[Z_{i}X_{i}^{\top}\gamma]$. Now $\mathbb{E}[Z_{i}\gamma_{0}]$ only depends on our design via the expected proportion of treated subjects. This will often be fixed by a budget constraint and even when it is not fixed, it does not depend on where specifically we assign the treatment. Therefore, for design purposes we may focus on $T_{i}$.

Under the model (2) with treatment assignment from (11),

If $\eta=\gamma$, so that we assign treatment using the true treatment effect vector, then equation (12) shows that the best expected gain comes from taking $u=\ell=0$, which is an RDD centered at the mean. Moreover, the expected gain decreases monotonically as we increase $u$ beyond $0$ or decrease $\ell$ below $0$. This matches our intuition that we must sacrifice some short-term gain to improve on statistical efficiency. Ordinarily $\eta\neq\gamma$, and a poor choice of $\eta$ could break this monotonicity.

In the Gaussian case considered in Section 3.3, we can likewise derive an explicit formula for the expected gain as a function of $u$, $\eta$, and $\gamma$. Letting $T_{\bullet}=Z_{\bullet}X_{\bullet}^{\top}\gamma$, we have

Using the formula (23) for $\alpha$ in the Gaussian case, this is simply

We expect intuitively that $\eta=\gamma$ will maximize $\mathbb{E}[T_{\bullet}]$. To verify this we start by choosing $u=u(\eta)$ in a way that keeps the proportion of data in the three treatment regions constant. We do so by taking $u=u(\eta)=u_{0}\sqrt{\eta^{\top}\Sigma\eta}$ for some $u_{0}\geqslant 0$, and then

Let $\gamma_{w}=\Sigma^{1/2}\gamma$ and $\eta_{w}=\Sigma^{1/2}\eta$, using the same matrix square root in both cases. Then

is maximized by taking $\eta_{w}=\gamma_{w}$ or equivalently $\eta=\gamma$. Any scaling of $\eta=\gamma$ leaves this criterion invariant.

Working under the normalization $\eta^{\top}\Sigma\eta=1$, we can summarize our results in the Gaussian case as

With our normalization, $u^{2}=u_{0}^{2}\eta^{\top}\Sigma\eta=u_{0}^{2}$. Equations (13) and (14) quantify the tradeoff between efficiency and short-term gain, that come from choosing $u_{0}$. Greater randomization through larger $u_{0}$ increases efficiency, and, assuming that the sign of $\eta$ is properly chosen, decreases the short-term gain.

In this section, we return to the setting where $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ are fixed values but $Z_{i}$ are not yet assigned. We use lowercase $x_{i}$ to emphasize that they are non-random. We assume that any subset of the $x_{i}$ of size $d$ has full rank, as would be the case almost surely when they are drawn IID from a distribution whose covariance matrix is of full rank. The design problem is to choose $p_{i}=\mathbb{P}(Z_{i}=1)$.

For given $x_{i}$, the design matrix in (2) is

Introducing $p_{i+}=p_{i}$ and $p_{i-}=1-p_{i}$ we get

Our design criterion is to choose $p_{i\pm}$ to minimize

This problem is only well-defined for $n\geqslant d$, since otherwise the matrix does not have full rank for any choice of $p\in\mathbb{R}^{n}$ and the determinant is always zero. This criterion is convex in $\{p_{is}\!\mid\!1\leqslant i\leqslant n,s\in\{+,-\}\}$ by a direct match with Chapter 7.5.2 of Boyd and Vandenberghe (2004) over the convex domain with $0\leqslant p_{i\pm}\leqslant 1$ and $p_{i+}+p_{i-}=1$ for all $i$.

It will be simpler for us to optimize over $q\equiv 2p-1\in[-1,1]^{n}$, in which case

Absent any other constraints, we have seen that the RCT ($q_{i}=0$, $p_{i}={1}/{2}$ for all $i\leqslant n$) always minimizes (16). The constrained optimization of (16) can be cast as both a semi-definite program (Boyd and Vandenberghe, 2004) and a mixed integer second-order cone program (Sagnol and Harman, 2015).

This setting is close to the usual design measure relaxation. Instead of choosing $n_{i}=1$ point $(\tilde{x}_{i},Z_{i})$ for observation $i$ we make a random choice between $(\tilde{x}_{i},1)$ and $(\tilde{x}_{i},-1)$ for that point. The difference here is that we have the union of $n$ such tiny design problems.

In this section we detail a simulation based on a real data set of emergency department (ED) patients. The MIMIC-IV-ED database (Johnson et al., 2021) provided via PhysioNet (Goldberger et al., 2000) includes data on ED admissions at the Beth Israel Deaconess Medical Center between 2011 and 2019.

Emergency departments face heavy resource constraints, particularly in the limited human attention and beds available. It is thus important to ensure patients are triaged appropriately so that the patients in most urgent need of care are assigned to intensive care units (ICUs). In practice, this is often done via a scoring method such as the Emergency Severity Index (ESI), in which patients receive a score in $\{1,2,3,4,5\}$, with $1$ indicating the highest severity and $5$ indicating the lowest severity. MIMIC-IV-ED contains these values as acuity scores, along with a vector of vital signs and other relevant information about each patient.

Such a setting provides a very natural potential use case for tie-breaker designs. Patients arrive with an assortment of covariates, and hospitals acting under resource constraints must decide whether to put them in an ICU. A hospital or researcher may be interested in the treatment effect of an ICU bed; for example, a practical implication of such a question is whether to expand the ICU or allocate resources elsewhere (Phua et al., 2020). It is also of interest (Chang et al., 2017) to understand which types of patient benefit more or less from ICU beds to improve this costly resource allocation. Obviously, it is both unethical and counterproductive to assign some ICU beds by an RCT. Patients with high acuity scores must be sent to the ICU, and those with low acuity scores may be exposed unnecessarily to bad outcomes such as increased risk of acquiring a hospital-based infection (Kumar et al., 2018; Vranas et al., 2018). However, it may be possible to randomize “in the middle”, e.g., by randomizing for patients with an intermediate acuity scores such as $3$, or with scores in a range where there is uncertainty as to whether the ICU would be beneficial for that patient. Because such patients are believed to have similar severities, this would limit ethical concerns and allow for greater information gain.

The triage data set contains several vital signs for patients. Of these, we use all quantitative ones, which are: temperature, heart rate (HR), respiration rate (RR), oxygen saturation ($O_{2}$ Sat.), and systolic and diastolic blood pressure (SBP and DBP). There is also an acuity score for each patient, as described above. The data set contains 448,972 entries, but to arrive at a more realistic sample for a prospective analysis, we randomly select $200$ subjects among those with no missing or blatantly inaccurate entries. Our optimization was done using the CVXR package (Fu et al., 2020) and the MOSEK solver (ApS, 2022).

To carry out a full analysis of the sort described in this paper, we need a vector $\eta$, as in (3). In practice, one could assume a model of the form (2) and take $\eta=\hat{\gamma}$ for some estimate of $\gamma$ formed via prior data. Since we do not have any $Y$ values (which in practice could be some measure of survival, length of stay, or subsequent readmission), we will construct $\eta$ via the acuity scores, using the reasonable assumption that treatment benefit increases with more severe acuity scores.

We collapse acuity scores of $\{1,2\}$ into a group ($Y=1$) and acuity scores of $\{3,4,5\}$ into another ($Y=0$) and perform a logistic regression using these binary groups. The covariates used are the vital signs and their squares, the latter to allow for non-monotonic effects, e.g., the acuity score might be lower for both abnormally low and abnormally high heart rates. All covariates were scaled to mean zero and variance one. For pure quadratic terms the squares of the scaled covariates were themselves scaled to have mean zero and variance one. We also considered an ordered categorical regression model but preferred the logistic regression for ease of interpretability. Our estimated $\hat{\eta_{j}}$ are in Table 1.