Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,838 characters · 16 sections · 59 citation commands
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 Angrist2014leveling, intervening to prevent child abuse in some families kran:2022, programs to counter juvenile delinquency lips:cord:berg:1981 and sending some university students to remedial English classes aike:west:schw:carr:hsiu: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 cattaneo2022regression. The RDD can give a consistent nonparametric estimate of the treatment effect at $x=t$ hahntodd, 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 gelman2019high.
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 camp: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 owenvarian, li:owen:2022:tr. 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, gelm:hill:veht: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 cookthibodeau 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 (ref) outlines our notation and the regression model. Section (ref) introduces our notions of efficiency and D-optimality in the multivariate regression model. Theorem (ref) 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 (ref) 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 (ref) 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 (ref) 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 cookfedorov, but both our proof and the theorem particulars are distinct. Section (ref) illustrates this procedure on a hospital data set from MIMIC-IV-ED about which emergency room patients should receive intensive care. Section (ref) has a brief discussion of some additional context for our results.
The earliest tie-breaker design of which we are aware camp:1969 had a discrete running variable, and literally broke a tie by randomizing the treatment assignments for those subjects with $x=t$. For an imperfectly measured running variable, we might consider $|x-t|<\Delta$ to not be meaningfully large for some $\Delta>0$, so randomizing in that window is like breaking ties. This is the approach from boru:1975.
When we randomize for $x$ with $|x-t|< \Delta$, then setting $\Delta=0$ provides an RDD while $\Delta\to\infty$ provides an RCT. The TBD transitions smoothly between these extremes. Many of the early comparisons among these methods use the two-line regression model
Here, $Y_i$ is the response for subject $i$, $Z_i \in \{-1, 1\}$ is a binary treatment variable, and $\varepsilon_i$ is an IID error term. goldberger finds for $X_i\sim\mathcal{N}(0,1)$ that the RDD has about 2.75 times the variance of an RCT. jacob2012practical consider polynomials in $X$ up to third degree, with or without interactions with $Z$ for Gaussian and for uniform random variables. Neither paper considers tie-breakers. cappelleri1994power study the two-line model above (absent the $X$-$Z$ interaction) and compare the RDD, RCT and three intermediate TBD designs. They study the sample size required to reject $\gamma=0$ with sufficient power for those designs and three different sizes of $\gamma$. Larger $\Delta$ brings greater power.
owenvarian find an expression for the variance of the parameters in the two-line model for $x$ with a symmetric uniform distribution about a threshold $t=0$, and also for a Gaussian case. The randomization window is $|x|\leqslant\Delta$, with a larger $\Delta$ yielding smaller variance but a lesser short-term gain. Their designs have only three levels (0, 50 and 100 percent) for the treatment probability given $x$. They show that there is no advantage to using other levels or a sliding scale when $x$ has a symmetric distribution and $50$% of the subjects will get the treatment.
li:owen:2023 study the two-line regression for a general distribution of $x$ and an arbitrary proportion of treated subjects. They constrain the global treatment probability and a measure of short-term gain. They find that there exists a D-optimal design under these constraints with a treatment probability that is constant within at most four intervals of $x$ values. Moreover, with the addition of a monotonicity constraint, there exists an optimal solution with just two levels corresponding to $x < t'$ and $x > t'$ for some $t'$.
klugerowen consider tie-breaker design for nonparametric regression for real-valued treatments $x$. A tie-breaker design can support causal inference about the treatment effect at points within the randomization window, not just at the threshold $x=t$. At the threshold, the widely use kernel regression methods for RDD become much more efficient if one has sampled from a tie-breaker, because they can use interpolation to $x=t$ instead of the extrapolation from the left and right sides that an RDD requires.
The tie-breaker setup is a special case of marginally constrained experimentation of cookthibodeau. Some more general constraints are considered in lopez2004optimal. Some of our findings are special cases of more general results in nachtsheim1989design. heavlin1998columnwise consider a sequential design setting in which the columns of a design matrix correspond to steps in semi-conductor fabrication, with each step constrained by the prior ones. The marginally constrained literature generally considers choosing $n_i$ values of the settable variables for level $i$ of the fixed variables. In our setting, we cannot assume that $n_i>1$.
The TBD setting often has a monotonicity constraint on treatment probabilities that we have not seen in the marginally constrained design literature. Under that constraint a more “deserving” subject should never have a lower treatment probability than a less deserving subject has. This and other constraints in the TBD will often lead to optimal designs with a small number of different treatment levels.
We anticipate that the tie-breaker design could be used in sequential experimentation. In our motivating applications, the response is measured long enough after the treatment (e.g., six years in some educational settings) that bandit methods slivkins2019introduction are not appropriate. There is related work on adaptive experimental design. See for instance metelkina. The problem we focus on is designing the first experiment that one might use, and a full sequential analysis is not in the scope of this paper.
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} =
\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 (ref) we will show how to get optimal $p_i$ by convex optimization.
To get more general insights into the design problem, in Section (ref) we also consider a random data framework. The predictors are to be sampled with $X_i\stackrel{\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=
^{\top}$ and $\tilde{\gamma} =
^{\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 (ref) generalizes the two-line model (ref) studied by li:owen:2023 and owenvarian. The latter authors describe some computational methods for the model (ref) but most of their theory is for model (ref).
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 (ref) such as spline models may prove more flexible in highly nonlinear settings, but we do not explore them here.
In Section (ref), 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 (ref) 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 (ref) 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 (ref). 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 (ref) how to accommodate any convex constraints on treatment assignment.
The assignment (ref) greatly generalizes the one in owenvarian 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 (ref) briefly outlines efficiency and D-optimality in the fixed $Z$ setting to introduce relevant formulas. In Section (ref), 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 chaloner1995bayesian, 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 (ref), since it corresponds to a point mass distribution on $X$. We return to this case in Section (ref), in which we show how to use convex optimization to obtain prospectively D-optimal $p(X_i)$.
In this brief section, we outline our notions of efficiency and D-optimality in the setting in which $X$ is fixed and $Z$ are assigned non-randomly. While not the focus of this paper, it will allow us to motivate various definitions and derive a property of D-optimality in the multivariate model that will be useful going forward.
We begin by conditioning on $X$ and $Z$. Let $D \in \mathbb{R}^{n \times n}$ be the diagonal matrix whose diagonal entries are $D_{ii} = Z_i$. We can write the linear model (ref) in matrix form as $Y = U \delta + \varepsilon$, where $U =
$ and $\delta =
^{\top}$.
In the general model (ref), conditionally on $X$ and $Z$ we have
where we assume here that $U$ is full rank. Because $\sigma^2$ is merely a multiplicative factor independent of all relevant parameters, it is no loss of generality to take $\sigma^2 = 1$ going forward for simplicity. The treatment effect vector $\tilde{\gamma}$ is our parameter of primary interest, so we want to minimize a measure of the magnitude of $((U^{\top}U)^{-1})_{22}$. When $X$ is fixed, a standard choice would be to assign $\{Z_1, \ldots, Z_n\}$ to minimize the $D$-optimality criterion of $$\det\bigl(((U^{\top}U)^{-1})_{22}\bigr)=\prod_{j=1}^{d + 1}\lambda_j(\mathrm{Var}(\hat\gamma \!\mid\! Z)),$$ where $\lambda_j(\cdot)$ is the $j$th eigenvalue of its argument. D-optimality is the most studied design choice among many alternatives atkinson2007optimum. It has the convenient property of being invariant under reparametrizations. This criterion is actually a particular case of D$_S$-optimality, in which only the parameters corresponding to a subset $S$ of the indices are of interest. See Section 10.3 of atkinson2007optimum. Much of our theory and discussion will generalize to a broader class of criteria that we discuss in Section (ref).
Before proceeding to the setting in which $X$ and $Z$ are random, we note a helpful result. Under the model (ref), there is a convenient property of D-optimality in this setting, which we state as the following simple theorem.
The equivalence (ref) is well-known in the D$_S$-optimality literature (see, e.g., Equation 10.7 of atkinson2007optimum). In our setting, only the right half of the columns of $U$ can be changed. Then the numerator in Equation (ref) is conveniently fixed and so optimizing the denominator alone optimizes the ratio.
The simple structure of the model (ref) has made $D$-optimal estimation of $\tilde\gamma$ equivalent to D-optimality for $\delta$ and for $\tilde\beta$. By the same token, a design that is A-optimal for $\tilde\gamma$, minimizing $\mathrm{tr}(\mathrm{Var}(\hat{\tilde\gamma}))$ is also A-optimal for $\delta$ and for $\tilde\beta$. Lemma 1 of nachtsheim1989design includes our setting and shows that $D$-optimality for $\delta$ is equivalent to $D$-optimality for $\tilde\gamma$. It does not apply to $\tilde\beta$ for our problem nor does that result consider other criteria such as A-optimality.
The theory of marginally restricted D-optimality in cookthibodeau describes some settings where the D-optimal design for all variables is a tensor product of the given empirical design for the fixed variables and a randomized design for the settable variables. Such a design is simply an RCT on the settable variables. By their Lemma 1, this holds when the regression model is a Kronecker product of functions of fixed variables times functions of settable variables. In their Lemma 3, this holds when the regression model has an intercept plus a sum of functions of settable variables and a sum of functions of fixed variables. Neither of those apply to model (ref) but Lemma 2 of nachtsheim1989design does. The TBD designs we consider usually have constraints on the short-term gain or monotonicity constraints, and those generally make RCTs non-optimal.
In this section, we modify the approach in Section (ref) to account for the randomness in $X$ and $Z$. To do so, we adopt a prospective D-optimality criterion to apply to the setting where both $X_i$ and $Z_i$ have yet to be observed. We derive our approach from ideas in Bayesian optimal design, which we briefly summarize below based on Chapter 18.2 of atkinson2007optimum, omitting details that will not play a role in our setting.
Bayesian optimal design often arises when the variance of the parameter estimate depends on the unknown true value of the parameter $\theta$, as is the usual case for models where the expected response is a nonlinear function of $\theta$. In this case, the information matrix $M=M(\theta)$ is usually constructed by linearizing the model form and then taking the expected outer product of the gradient with itself under a design measure on the predictors. This $M(\theta)$ is random because $\theta$ has a prior distribution. In our case, $M = U^{\top}U$ since this is the information matrix for the multivariate regression.
The favored approach in Bayesian optimal design is to choose the design in order to minimize $\mathbb{E}( \log \det (M^{-1}))$ where the expectation is over the prior distribution on the parameters chaloner1995bayesian, dette1996. That can be quite expensive to do. Many of the examples in the literature optimize the design over a grid which is reasonable when the dimensions of $X$ and $\theta$ are both small, but those methods do not scale well to larger problems.
Table 18.1 of atkinson2007optimum lists four additional criteria along with the above one. They are $\log\mathbb{E}(\det(M^{-1}))$, $\log\det(\mathbb{E}(M^{-1}))$, $\log(\mathbb{E}(\det(M))^{-1})$, and $\log(\det\mathbb{E}(M))^{-1}$. The objective becomes more tractable each time a nonlinear operation is taken out of the expectation. When the logarithm is the final step, the criterion is equivalent to not taking that logarithm.
We choose the last of those four quantities (choice V in their Table 18.1) for our definition of prospective $D$-optimality.
We could analogously define prospective D-optimality for $\tilde{\beta}$ or $\tilde{\gamma}$ as minimizing $\det((\mathbb{E}[U^{\top}U]^{-1})_{11})$ or $\det((\mathbb{E}[U^{\top}U]^{-1})_{22})$, respectively. By Theorem (ref), prospective D-optimality in the sense of Definition (ref) is equivalent to these conditions in our model, so the three notions all align.
Our choice is known as EW D-optimality in Bayesian design for generalized linear models (GLMs) YMM2016, Bu2020. The E is for expectation and the $W$ refers to a weight matrix arising in GLMs. It is valued for its significantly reduced computational cost. While a computationally efficient design is not necessarily the best one, YMM2016 and YTM2017 find in a range of simulation studies that EW D-optimal designs tend to have strong performance under the more standard Bayesian D-optimality criterion as well.
Since our prospective D-optimality criterion uses the expected value of $U^{\top}U$, it depends only on the $p_i$ terms and not on the joint distribution of the $Z_i$. It is therefore important that the $Z_i$ be independent in order to distinguish, say, an RCT with $\mathbb{P}(Z_i=1)=1/2$ from an allocation that takes all $Z_i=Z_1$ where $Z_1=1$ with probability $1/2$. In addition, morr:etal:2024 show in a similar problem that, when the $Z_i$ are sampled independently, the design that minimizes $\det(\mathbb{E}[M])$ also minimizes $\mathbb{E}[\det(M)]$ asymptotically as $n \to \infty$. klugerowen consider some stratified sampling methods that incorporate negative correlations among the $Z_i$ which makes $U^{\top}U$ come even closer to its expectation than it does under independent sampling.
Under sampling with $X_i\sim P_X$, define
where the bullet subscript denotes an arbitrary subject with $X_\bullet\sim P_X$ and $\mathbb{P}(Z_{\bullet} = 1 \text{ } | \text{ } X_{\bullet}) = p(X_{\bullet})$. In addition,
Now let $N$ be the matrix with
Under our sampling assumptions
The right-hand side of (ref) represents the expected information per observation in our tie-breaker design. Using these formulas, we obtain the following desirable result for prospective D-optimality.
The proofs of Theorem (ref) and of all subsequent theorems are presented in Appendix (ref). Theorem (ref) does not require $X_i$ to be independent though that would be the usual model.
The theorem establishes that the RCT is prospectively D-optimal among any randomization scheme $\mathbb{P}(Z=1\!\mid\! X_\bullet) =p(X_\bullet)\in[0,1]$. It is not necessarily the unique optimum in this larger class. For instance, if $$\mathbb{E}[ X_{\bullet j}X_{\bullet k}(2p(X_\bullet)-1)]=0$$ for all $j$ and $k$, then the function $p(\cdot)$ would provide the same efficiency as an RCT since it would make the matrix $N$ in the above proof vanish.
Though we frame Theorem (ref) via prospective D-optimality, the result also holds for a broader class of criteria that we call prospectively monotone. The basic idea of such criteria is that they only depend on the bottom-right submatrix of $\mathbb{E}[(U^{\top}U)^{-1}]$ and that they encourage this submatrix to be small in the standard ordering on positive semi-definite matrices (i.e., $\Sigma_1 \preceq \Sigma_2$ if and only if $\Sigma_2 - \Sigma_1$ is PSD). The precise definition is as follows.
Because these are the only properties of prospective D-optimality we use in Theorem (ref), the same result holds immediately for any prospectively monotone criterion. Examples include prospective A-optimality, which minimizes the quantity $\text{Tr}(\mathbb{E}[\text{Var}(\hat{\tilde{\gamma}} \!\mid\! Z)^{-1}])$, or prospective C-optimality, which minimizes $c^{\top}\mathbb{E}[\text{Var}(\hat{\tilde{\gamma}} \!\mid\! Z)^{-1}]c$ for some preset vector $c$ (the latter of use if a particular linear combination of elements of $\hat\gamma$ is of primary interest).
We can gain particular insights, both theoretical and practical, by considering a special case satisfying two conditions. First, $P_X$ has a symmetric density, i.e., $f_X(\vec{x}) = f_X(-\vec{x})$ for $\vec{x} \in \mathbb{R}^d$. This includes the special case of Gaussian covariates, which we consider in more detail as well. Secondly, we will further assume that $p = 1/2$ and the randomization window is symmetric about zero with width $\Delta\geqslant0$, which we call a symmetric design. That is, we restrict (ref) to simply
For $j=k=0$ (i.e., both terms are intercepts), equation (ref) reduces to $$N_{00} = \mathbb{E}[(\mathds{1}\{X_{\bullet}^{\top} \eta \geqslant \Delta\} - \mathds{1}\{X_{\bullet}^{\top} \eta \leqslant -\Delta\})] = 0$$ since we are integrating an odd function with respect to a symmetric density. Likewise, when both $j,k\geqslant1$ we have $N_{jk} = 0$. The only cases that remain are the first row and first column of $N$, besides the top-left entry. Thus, we can write
where $\alpha \in \mathbb{R}^d$ with
We note that $\alpha = \alpha(\Delta, \eta)$ depends on the width $\Delta$ and the treatment assignment vector $\eta$, but we suppress that dependence for notational ease. From (ref), we can compute explicitly that $$ N \tilde{\Sigma}^{-1} N =
,$$ so our criterion becomes
In the last line we use the formula $\det(A + cd^{\top}) = \det(A) (1 + d^{\top} A^{-1} c)$ for the determinant of a rank-one update of an invertible matrix and we also note that $\det(\tilde \Sigma)=\det(\Sigma)$. Let $W = \Sigma^{1/2}$ so that $\mathrm{Var}(W^{-1}x) = I$. The efficiency therefore only depends on $\alpha$ through $\alpha^{\top}\Sigma^{-1}\alpha = \Vert W^{-1}\alpha\Vert^2$.
We could also ask whether we can do better by changing our randomization scheme to allow
for some other $p \neq 1/2$. While this may be a reasonable choice in practice when treatment cannot be assigned equally, it cannot provide any efficiency benefit in the symmetric case, as shown in Theorem (ref) below. Just as an RCT is most efficient globally, if one is using the three level rule (ref) then the best choice for the middle level is $1/2$ and that choice is unique under a reasonable assumption.
An informative example is the case in which $P_X = \mathcal{N}(0, \Sigma)$ for some covariance matrix $\Sigma$. In this case, we can compute the efficiency explicitly as a function of $\Delta$.
From Theorem (ref) we find that the efficiency ratio between $\Delta = \infty$ (the RCT) and $\Delta = 0$ (the RDD) is $(1 - {2}/{\pi})^{-2} \approx 7.57$. The result in goldberger gives a ratio of $(1-2/\pi)^{-1}$ for the variance of the slope in the case $d = 1$. Our result is the same, though we pick up an extra factor because our determinant criterion incorporates both the intercept and the slope. Their result was for $d=1$; here we get the same efficiency ratio for all $d\geqslant1$.
In this multivariate setting we see that for any fixed $\Delta>0$, the most efficient design minimizes $\eta^{\top} \Sigma \eta$, so it is the eigenvector corresponding to the smallest eigenvalue of $\Sigma$. This represents the least “distribution-aware” choice, i.e., the last principal component vector, which aligns with our intuition that we gain more information by randomizing as much as possible.
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_iX_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 (ref) with treatment assignment from (ref),
If $\eta = \gamma$, so that we assign treatment using the true treatment effect vector, then equation (ref) 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 (ref), 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 (ref) for $\alpha$ in the Gaussian case, this is simply $$\mathbb{E}[T_{\bullet}] = \sqrt{\frac{2}{\pi}} \frac{\gamma^{\top} \Sigma \eta}{\sqrt{\eta^{\top} \Sigma \eta}} \text{ } e^{\frac{-u^2}{2 \eta^{\top}\Sigma\eta}}.$$
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\geqslant0$, and then $$\mathbb{E}[T_{\bullet}] = \sqrt{\frac{2}{\pi}} \frac{\gamma^{\top} \Sigma \eta}{\sqrt{\eta^{\top} \Sigma \eta}} \text{ } e^{{-u_0^2}/{2}}.$$ Let $\gamma_w=\Sigma^{1/2}\gamma$ and $\eta_w=\Sigma^{1/2}\eta$, using the same matrix square root in both cases. Then $$ \frac{\gamma^{\top}\Sigma\eta}{\sqrt{\eta^{\top}\Sigma\eta}} =\frac{\gamma_w^{\top}\eta_w}{\Vert\eta_w\Vert} $$ 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 (ref) and (ref) 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 (ref) is $$U =
, \quadfor\quad u_i(1) =u_{i+}\equiv
\quadand\quad u_i(-1) =u_{i-}\equiv
. $$ 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 boyd2004convex 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 2 p - 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 (ref). The constrained optimization of (ref) can be cast as both a semi-definite program boyd2004convex and a mixed integer second-order cone program sagnolharman.
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.
We outline a few reasonable (and convex) constraints that one could impose on this problem:
Budget constraint: In practice we have a fixed budget for the treatments. For instance the number of scholarships or customer perks to give out may be fixed for economic reasons. We can impose this constraint in expectation by setting $(1/n) \sum_{i = 1}^{n} p_{i} = \mu$, where $\mu$ is some fixed average treatment rate.
Monotonicity constraint: It may be reasonable to require that $p_{i}$ is nondecreasing in some running variable $r_i = f(\tilde{x}_i)$. For example, a university may require that an applicant's probability of receiving a scholarship can only stay constant or increase with their score, which is some combination of applicant variables. We can encode this as a convex constraint by first permuting the data matrix so that $r_{(1)} \leqslant r_{(2)} \leqslant \cdots \leqslant r_{(n)}$ and then forcing $p_{(1)} \leqslant p_{(2)} \leqslant \cdots \leqslant p_{(n)}$. Note that the formulation (ref) satisfies this monotonicity constraint, in which case $r_i = \tilde{x}_i^{\top}\eta$.
Gain constraint: One may also want to impose that the expected gain is at least some fraction of its highest possible value, i.e.
The left-hand side of (ref) is the expected gain for this choice of $p_i$, whereas the right-hand side is the highest possible gain, which corresponds to the RDD $\mathbb{P}(Z_i = 1) = \mathds{1}\{\tilde{x}_i^{\top}\eta \geqslant 0\}$. Because $\gamma$ is typically not known exactly, (ref) computes the anticipated gain under the sampling direction $\eta$ we use. If an analyst has access to a better estimate of gain, such as a set of estimated treatment effects $\{\hat{\tau}(\tilde{x}_1),\ldots, \hat{\tau}(\tilde{x}_n)\}$ fit from prior data, they could replace (ref) by $\sum_{i = 1}^{n} p_i \hat{\tau}(\tilde{x}_i) \geqslant \rho \sum_{i = 1}^{n} |\hat{\tau}(\tilde{x}_i)|$, which is again a linear constraint on $q$.
Covariate balance constraint: The expected total value of the $j$th variable for treated subjects is $\sum_i p_i \tilde{x}_{ij}$. The constraint $$\sum_ip_i(\tilde{x}_{ij}-\Delta_j) \leqslant 0$$ keeps the expected sample average of the $j$th variable among treated subjects to at most $\Delta_j$ regardless of the actual number of treated subjects. We note that sufficiently stringent covariate balance constraints may not be compatible with gain or monotonicity constraints.
Though we delay discussion of a particular example to Section (ref), we prove here an empirically-observed phenomenon regarding the monotonicity constraint. In particular, optimal solutions when a monotonicity constraint is imposed display only a few distinct levels in the assigned treatment probability. This is consistent with the one-dimensional theory of li:owen:2023, though interestingly that paper observed the same behavior even with only a budget constraint. It also used an approach that does not obviously extend to $d>1$.
These upper bounds closely resemble those in cookfedorov, who study constraints on a design of $\tilde{x}_1,\ldots,\tilde{x}_n$ themselves rather than on a treatment variable. However, theirs is an existence result for some solution with this level of sparsity, whereas our result holds for any solution with probability one. Moreover, their result derives from Caratheodory's theorem, whereas ours comes from a close analysis of the KKT conditions.
The upper bounds in Theorem (ref) are not tight. In our experience, there have typically been no more than five or six distinct levels even for $d$ as high as ten.
In this section we detail a simulation based on a real data set of emergency department (ED) patients. The MIMIC-IV-ED database mimic provided via PhysioNet physionet 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 phua2020lessmore. It is also of interest chang2017icu 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 kumar2018healthcare, vranas2018lowrisk. 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 cvxr and the MOSEK solver mosek.
To carry out a full analysis of the sort described in this paper, we need a vector $\eta$, as in (ref). In practice, one could assume a model of the form (ref) 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 (ref).
Figure (ref) presents the efficiency/gain tradeoff as we vary the size of the randomization window $\Delta$ in (ref). For ease of visualization, we plot the logs of both quantities. As expected, we get a clear monotone increase in efficiency and decrease in gain as we increase $\Delta$, moving from an RDD to an RCT. It should be noted that our efficiency criterion, because it only uses information in the $X$ values, is robust to a poor choice of $\eta$, whereas our gain definition is constrained by the assumption that $\eta$ is a reasonably accurate stand-in for the true treatment effect $\gamma$.
In practice, it is hard to interpret what a “good” value of efficiency is because of our D-optimality criterion. Hence, as in owenvarian, a pragmatic approach is to first stipulate that the gain is at least some fraction of its highest possible value, and then pick the largest $\Delta$ for this choice to maximize efficiency. A more qualitative choice based on results like Figure (ref), such as picking the right endpoint of a sharp efficiency jump or the left endpoint of a sharp gain decline, would also be sensible.
As we see in Figure (ref), the treatment constraint causes most $p_i$ to be at or near zero or one. Adding the gain constraint pushes most of the treatment probabilities to zero for low values of the running variable and one for high values. This scenario most closely resembles the RDD, with some deviations to boost efficiency. Indeed, the optimal solution would necessarily tend towards the RDD solution as the gain constraint increased. Finally, the monotonicity constraint further pushes the higher values of $p$ to the positive values of the running variable and vice-versa, since we lose the opportunity to counterbalance some high and low probabilities at the extreme with their opposites. The right two panels display a few discrete levels in the treatment probability, consistent with Theorem (ref).
In this paper, we add to a growing body of work demonstrating the benefits of tie-breaker designs. Though RCTs are often infeasible, opportunities for small doses of randomization may present themselves in a wide variety of real-world settings, in which case treatment effects can be learned more efficiently. This phenomenon is analogous to similar causal inference findings about merging observational and experimental data rct+odb,rosenman2023combining,colnet2024causal.
The convex optimization framework in Section (ref) is more general and conveniently only relies on knowing sample data rather than population parameters. It is also simple to implement and allows one to incorporate natural economic and ethical constraints with ease.
Multivariate tie-breaker designs are a natural option in situations in which there is no clear univariate running variable. For example, subjects may possess a vector of covariates, many of which could be associated with heterogeneous treatment effects in some unknown way of interest. Of course, two-line models and their multivariate analogs are not nearly as complicated as many of the models found in practice. Our view is to use them as a working model by which to decide on treatment allocations, in which case a more flexible model could be used upon full data acquisition as appropriate.
We thank John Cherian, Anav Sood, Harrison Li and Dan Kluger for helpful discussions. We also thank Balasubramanian Narasimhan for helpful input on the convex optimization problem, Michael Baiocchi and Minh Nguyen of the Stanford School of Medicine for discussions about triage to hospital intensive care units, and an anonymous reviewer for helpful comments. This work was supported by the NSF under grants IIS-1837931 and DMS-2152780. T.\ M.\ is supported by a B.\ C.\ and E.\ J.\ Eaves Stanford Graduate Fellowship.