# Designing Experiments Toward Shrinkage Estimation

Evan T. R. Rosenman, Luke Miratrix

###### Abstract

We consider how increasingly available observational data can be used to
improve the design of randomized controlled trials (RCTs). We seek to design a
prospective RCT, with the intent of using an Empirical Bayes estimator to
shrink the causal estimates from our trial toward causal estimates obtained
from an observational study. We ask: how might we design the experiment to
better complement the observational study in this setting?

We propose using $\boldsymbol{\kappa}_{2}$ – which shrinks each component of
the RCT estimator toward its observational counterpart by a factor
proportional to its variance – as our estimator. First, we show that the risk
of $\boldsymbol{\kappa}_{2}$ can be computed efficiently via numerical
integration. We then propose algorithms for determining the best allocation of
units to strata (the best “design”). We consider three options: Neyman
allocation; a “naïve” design assuming no unmeasured confounding in the
observational study; and a “defensive” design accounting for the imperfect
parameter estimates we would obtain from the observational study with
unmeasured confounding.

We also incorporate results from sensitivity analysis to establish guardrails
on the designs, so that our experiment could be reasonably analyzed with and
without shrinkage. We demonstrate the superiority of these experimental
designs with a simulation study involving causal inference on a rare, binary
outcome.

###### Contents

  1. 1 Introduction
  2. 2 Set-up
    1. 2.1 Notation
    2. 2.2 Assumptions and Loss Function
    3. 2.3 Related Problems
    4. 2.4 Principles Guiding Estimator Choice and Experimental Design
  3. 3 Designing Towards Shrinkage
    1. 3.1 Choice of Estimator: $\boldsymbol{\kappa}_{2}$
    2. 3.2 Exact Risk Calculation Under Known Parameters
    3. 3.3 Design Options
      1. 3.3.1 Neyman Allocation
      2. 3.3.2 Heuristic Optimization Assuming $\boldsymbol{\xi}\boldsymbol{=0}$
      3. 3.3.3 Heuristic Optimization Assuming Worst-Case Error Under $\boldsymbol{\Gamma}$-Level Unmeasured Confounding
    4. 3.4 Imposing Guardrails on Designs
  4. 4 Simulation Study
    1. 4.1 Simulation Set-Up
    2. 4.2 Ideal Case: No Unmeasured Confounding
    3. 4.3 Practical Case: Unmeasured Confounding
  5. 5 Discussion
  6. A Alternatives Shrinkers to $\boldsymbol{\kappa}_{2}$

##  1 Introduction

Recent years have seen increased interest in methods to integrate
observational data with experimental data (Colnet et al.,, 2020). Such methods
have been used to identify average causal effects in target populations
(Bareinboim and Pearl,, 2016; Kallus et al.,, 2018), identify heterogeneous
treatment effects (Peysakhovich and Lada,, 2016), and improve precision in
causal estimation (Gagnon-Bartsch et al.,, 2021)

This surge in methodological development is motivated, at least in part, by
the increasing proliferation of observational databases. Such repositories
provide statisticians with rich new data sources from which to learn. Yet the
lurking danger of unmeasured confounding yields rightful trepidation about
incorporating these data into estimation procedures (Colnet et al.,, 2020).
Rosenman et al., (2020) developed a procedure for shrinking causal estimates
from a stratified experiment toward the analogous estimates from an
observational study. Shrinkage estimators are attractive in that they allow
researchers to use observational data in tandem with experimental data, while
protecting the integrity of the randomization of the experiment. They provide
a guaranteed reduction in expected loss, relative to using the experimental
data alone.

Separately, Rosenman and Owen, (2021) proposed a method to design more
powerful stratified experiments by utilizing information from an observational
study to inform decisions about how to allocate a sample across strata and
treatment arms. Risk reductions from this method are more modest, owing to the
fact that the observational data is used only for design and not for
inference.

Here, we combine the two approaches of design and shrinkage. We seek to design
a prospective randomized trial. We optimize the design with two objectives in
mind. Our intent is to shrink the causal estimates from our designed trial
toward those of an observational study that is already completed. First, we
want to design the experiment as a better complement to the observational
data, allowing for significant gains in estimation precision. At the same
time, we want the experiment to be usable in its own right, and seek to impose
guardrails on the design such that the experiment will yield sufficiently
precise causal estimates even if we do not elect to include the observational
study in the final estimation.

The remainder of this paper proceeds as fellow. In Section 2, we define our
problem and introduce notation and assumptions. Section 3 introduces our
choice of shrinkage estimator, $\boldsymbol{\kappa}_{2}$, and demonstrates how
to compute its risk efficiently. This section also discusses three different
heuristics under which analysts can design prospective experiments with the
intent of leveraging $\boldsymbol{\kappa}_{2}$ on the final results, while
also protecting the utility of the experiment on its own. Section 4 contains
two simulation studies, highlighting the risk improvements that can be
attained by designing toward shrinkage. Section 5 concludes.

##  2 Set-up

###  2.1 Notation

We operate in the stratified setting, with fixed strata $k=1,\dots,K$ defined
by covariates. Such strata could be defined by subject matter knowledge, or
they could be defined based on application of a modern machine learning method
for uncovering heterogeneous treatment effects (Hill,, 2011; Wager and Athey,,
2018). We begin by setting up the problem in full generality, such that we
have access to a pilot dataset that could be either observational or
experimental.

The pilot dataset comprises $n_{\pi}$ total units, indexed by $j$. With each
unit, we associate two binary potential outcomes, $Y_{j}(0)$ and
$Y_{j}(1)\in\\{0,1\\}$, representing the unit’s outcome in the presence and
absence of treatment. We also define a treatment variable $W_{j}\in\\{0,1\\}$,
representing whether unit $j$ is treated or not; and a vector
$\boldsymbol{X}_{j}\in\mathbb{R}^{p}$ representing measured covariates.

The quartets $(\boldsymbol{X}_{j},W_{j},Y_{j}(0),Y_{j}(1))$ are sampled i.i.d.
from some population distribution $F_{\pi}$. Denote as
$\mathbb{E}_{\pi},\mathrm{var}_{\pi},$ and $\mathbb{P}_{\pi}$ the expectation,
variance, and probability operators under $F_{\pi}$. Covariate membership is
defined by a variable $S_{j}$ such that
$S_{j}=k\iff\boldsymbol{X}_{j}\in\mathcal{X}_{k}$ for some set of covariate
values $\mathcal{X}_{k}$.

We design a future blocked experiment of the same treatment. In the
experiment, we will recruit a fixed total of $n_{rk}$ units for each stratum
$k$, and we will then randomize $n_{rkt}$ of those units to receive the
treatment and $n_{rkc}=n_{rk}-n_{rkt}$ to receive the control. We call the
list of tuples $\boldsymbol{d}=\\{n_{rkt},n_{rkc}\\}_{k=1}^{K}$ the “design.”
We will assume a sample size constraint such that

$\sum_{k}n_{rkt}+n_{rkc}=n_{r}\,,$

for some fixed number $n_{r}$.

Denote as $F_{r}$ the sampling distribution for the triplets
$(X_{i},Y_{i}(0),Y_{i}(1))$ among units $i$ in the experimental population. We
refer to the conditional sampling distribution for units in stratum $k$ as
$F_{r\mid S_{i}=k}$. The experiment can then be understood as drawing $n_{rk}$
units from $F_{r\mid S_{i}=k}$ and then assigning $W_{i}$ by choosing a simple
random sample of size $n_{rkt}$ from the set of $n_{rk}$ recruited units. We
define as $\mathbb{E}_{r},\mathrm{var}_{\mathrm{r}},$ and $\mathbb{P}_{r}$ the
expectation, variance, and probability operators over both the sampling and
treatment randomization in the experiment.

###  2.2 Assumptions and Loss Function

We make the following standard assumption.

######  Assumption 1 (Consistency).

For each unit $\ell$ in the pilot study or the future experiment, the observed
outcome $Y_{\ell}\in\\{0,1\\}$ is given by

$Y_{\ell}=W_{\ell}Y_{\ell}(1)+(1-W_{\ell})Y_{\ell}(0)\,,$

that is, there is only one “version” of the treatment.

A key assumption is that the causal effects across strata are shared between
the pilot and future RCT datasets, i.e.

######  Assumption 2 (Common Treatment Effects).

For all strata $k=1,\dots,K$,

$\mathbb{E}_{\pi}\left(Y(1)-Y(0)\mid
S=k\right)=\mathbb{E}_{r}\left(Y(1)-Y(0)\mid S=k\right)=\tau_{k}\,.$

Assumption 2, sometimes known as the “transportability condition” (Colnet et
al.,, 2020), imposes a congruency constraint on the distributions $F_{\pi}$
and $F_{r}$. This assumption gives us a common target of estimation, which we
define as

$\boldsymbol{\tau}=\left(\tau_{1},\dots,\tau_{k}\right)^{\mathsf{T}}\,.$

We will typically invoke a slightly stronger version of Assumption 2,
described below.

######  Assumption 3 (Common Moments).

For all strata $k=1,\dots,K$, and $w\in\\{0,1\\}$:

$\mathbb{E}_{\pi}\left(Y(w)\mid S=k\right)=\mathbb{E}_{r}\left(Y(w)\mid
S=k\right)\hskip 5.69054pt\text{ and }\hskip
5.69054pt\mathrm{var}_{\pi}\left(Y(w)\mid
S=k\right)=\mathrm{var}_{r}\left(Y(w)\mid S=k\right)\,.$

We define

$\boldsymbol{\hat{\tau}_{\pi}}=\left(\hat{\tau}_{\pi 1},\dots,\hat{\tau}_{\pi
K}\right)^{\mathsf{T}},$

as the set of stratum causal estimates arising from the pilot study. These
estimates can be obtained using a difference-in-means estimator, or a more
complex estimator. Though not yet realized, denote as
$\boldsymbol{\hat{\tau}_{r}}$ the vector of difference-in-means estimates from
the future experiment.

Our eventual causal estimator,

$\boldsymbol{\hat{\tau}}\equiv\boldsymbol{\hat{\tau}}\left(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{\pi}}\right)=(\hat{\tau}_{1},\dots,\hat{\tau}_{K})^{\mathsf{T}},$

will be a function of both $\boldsymbol{\hat{\tau}_{\pi}}$ and
$\boldsymbol{\hat{\tau}_{r}}$. Our target of estimation in this manuscript is
_not_ an average treatment effect over the experimental population. Rather, we
seek to obtain good estimates simultaneously for all of the stratum causal
effects $\boldsymbol{\tau}$. Under Assumption 2, we can define a loss function
under which we evaluate
$\boldsymbol{\hat{\tau}}\left(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{\pi}}\right)$.
We use the simple, unweighted $L_{2}$ loss:

$\mathcal{L}(\boldsymbol{\tau},\boldsymbol{\hat{\tau}})=\sum_{k}\left(\tau_{k}-\hat{\tau}_{k}\right)^{2}\,.$

We seek to optimize the risk of $\hat{\boldsymbol{\tau}}$, conditional on
$\boldsymbol{\hat{\tau}_{\pi}}$, which yields the risk expression

$\displaystyle\mathcal{R}(\boldsymbol{\tau},\boldsymbol{\hat{\tau}})$
$\displaystyle=\mathbb{E}_{r}\left(\mathcal{L}(\boldsymbol{\tau},\boldsymbol{\hat{\tau}})\right)$
$\displaystyle=\sum_{k}\mathbb{E}_{r}\left(\left(\tau_{k}-\hat{\tau}_{k}\right)^{2}\right)\,.$

###  2.3 Related Problems

Under this framing, we can relate this problem to several problems in the
causal inference and experimental design.

First, suppose that the pilot study is an RCT, such that under $F_{\pi}$ we
have $W_{i}\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_{i}(0),Y_{i}(1)$. Then, this problem becomes
one of adaptive experimental design, in which a trial is conducted in multiple
phases and information from the first phase can be used to design later
phases. This problem has a rich history in the literature, stretching back to
the foundational work of Thompson, (1933, 1935). While Thompson sampling was
originally designed for a more generic problem involving maximizing the
expected reward, it can be used to estimate average stratum treatment effects,
as discussed in Offer-Westort et al., (2021). Adaptive experimental design is
an area of active research, though much recent work has focused on methods
that define strata in the second phase, rather than taking strata as fixed
(Tabord-Meehan,, 2018; Bai,, 2019). For modern methods that incorporate a
fixed stratification scheme, see Hahn et al., (2011) and Chambaz et al.,
(2014).

Next, suppose that the pilot study is an observational study, but that our
ultimate causal estimator $\boldsymbol{\hat{\tau}}$ will simply be
$\boldsymbol{\hat{\tau}_{r}}$, the vector of difference-in-means estimators
arising from the experiment. This problem is closely related to the survey
sampling work of Neyman, (1992). Neyman computed the optimal allocation for a
stratified survey – under a budget constraint – supposing pilot estimates of
variance could be obtained for each stratum. These ideas can easily be
extended to the causal inference setting. Here, unbiased pilot stratum
variance estimates are obtainable if $W_{i}\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_{i}(0),Y_{i}(1)\mid X_{i}$ — that is,
unconfoundedness holds in the observational study (Imbens and Rubin,, 2015).
In our prior work, Rosenman and Owen, (2021), we demonstrate that efficiency
gains are possible even in the presence of unmeasured confounding, as long as
we can bound the magnitude of the unmeasured confounding using the sensitivity
model of Tan, (2006). This is because the observational data can indicate
parts of the covariate space where variation is higher or lower – and hence,
where experimenters should over- or undersample.

Lastly, suppose that the pilot study is an observational study, but that the
experimental study is already completed, and our goal is to choose an
estimator $\boldsymbol{\hat{\tau}}$ to trade off between causal estimates
derived from the two data sources. This is an example of a so-called “data
fusion” problem (Bareinboim and Pearl,, 2016) that seeks to merge the
observational and experimental data directly. Many methods rely on
unconfoundedness in the observational study, including one of our prior papers
(Rosenman et al.,, 2018) as well as Athey et al., (2019). Other studies have
sought to weaken this condition, and frequently utilize alternative
assumptions to proceed with merged estimations. Kallus et al., (2018) assumes
that the hidden confounding has a parametric structure that can be modeled
effectively. Peysakhovich and Lada, (2016) propose a method for when the
observational data are time series and the bias preserves a unit-level rank
ordering. For an excellent overview of the available methods in this area, see
Colnet et al., (2020).

###  2.4 Principles Guiding Estimator Choice and Experimental Design

In this paper, we will consider the case where the pilot study is an
observational study, and the experimental study has yet to be implemented.
Hence, we replace the subscript $\pi$ with $o$ (e.g. defining our vector of
pilot study causal estimates as $\boldsymbol{\hat{\tau}_{o}}$ rather than
$\boldsymbol{\hat{\tau}_{\pi}}$). We will choose an estimator
$\boldsymbol{\hat{\tau}}\left(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{o}}\right)$
for combining the observational and experimental data within each stratum.
Then, we will design our experiment explicitly to minimize the risk of this
estimator. In this section, we highlight several characteristics of the
estimator and experimental design procedure that we consider ideal.

First, we would like our estimator to be _robust to unmeasured confounding_ in
the observational study. The assumption of unconfoundedness – roughly, that
all variables affecting both the treatment probability and the outcome have
been measured in the observational study – is fundamentally untestable, and
rarely holds in practice (Imbens and Rubin,, 2015). In this setting, note that
our problem is somewhat asymmetric: a simple vector of difference-in-means
estimates from the experimental study will be unbiased, and will be “good
enough” in many cases. Hence, we do not want to incorporate the observational
data unless we have strong guarantees that it will reduce statistical risk.
Thus, our chosen estimator should not be highly susceptible to bias due to
unmeasured confounding in the observational study.

Second, we would like our procedure to generate experiments that are still
valid if they are analyzed alone. We term this feature _detachability_. To
motivate this idea, consider an extreme case where we fail to sample any
treated or any control units from a particular stratum, under the assumption
that the experimental units will simply complement observational units in
those strata. Suppose that we later learn of a problem with the observational
study that undermines the validity of Assumption 2 — e.g., perhaps the data
was obtained fraudulently. In this case, we would be out of luck: the
experiment has an effective sample size of $0$, so we cannot use it on its
own. Some measure of safety could be achieved by imposing design restrictions,
such as mandating that each stratum must contain a minimum proportion of
treated or control units. But we would like our procedure to naturally
regularize toward more equal experimental allocations such that detachability
is naturally achieved.

With these features in mind, we briefly discuss the possibility of using
simpler data fusion estimators. In Rosenman et al., (2018) we propose a method
that is based on simply pooling observational and experimental data within
each stratum (termed the “spiked-in” estimator) and also methods based on
convex combinations of $\boldsymbol{\hat{\tau}_{o}}$ of
$\boldsymbol{\hat{\tau}_{r}}$ (termed the “weighted average” and “dynamic
weighted average” estimators). These estimators were designed under the
assumption that unconfoundedness holds in the observational study. Hence, they
do not exhibit the desired robustness property. Hence, as discussed in the
next section, we turn our attention toward estimators based on Empirical Bayes
shrinkage.

##  3 Designing Towards Shrinkage

###  3.1 Choice of Estimator: $\boldsymbol{\kappa}_{2}$

In Green and Strawderman, (1991) and Green et al., (2005), the authors
consider how to shrink between an unbiased estimator and a biased estimator.
Their goal is to derive Empirical Bayes estimators that guarantee a risk
reduction relative to using the unbiased estimator alone. The problem turns
out to be quite similar to James-Stein estimation.

In Rosenman et al., (2020), we propose two estimators that are extensions of
Green and Strawderman’s work. The latter estimator, $\boldsymbol{\kappa}_{2}$
is constructed such that it shrinks each component of the unbiased estimator
toward its counterpart in the biased estimator by a factor inversely
proportional to its precision. We consider this property desirable in many
applied settings. Moreover, this estimator typically outperformed competitor
estimators – including those proposed by Green and Strawderman – in our
applied data analysis. Hence, we proceed in this manuscript with the intention
of using $\boldsymbol{\kappa}_{2}$ to shrink between our RCT and observational
study estimators. We discuss the use of alternative estimators in the
Appendix.

Let $\boldsymbol{\hat{\tau}_{r}}\in\mathbb{R}^{K}$ be the unbiased (RCT)
estimator and $\boldsymbol{\hat{\tau}_{o}}\in\mathbb{R}^{K}$ be the biased
(observational study) estimator. Denote as
$\boldsymbol{\Sigma}_{r}=\mathrm{diag}(\sigma_{rk}^{2})\in\mathbb{R}^{K\times
K}$ the diagonal covariance matrix of $\boldsymbol{\hat{\tau}_{r}}$. Under
mild conditions, we can assume
$\boldsymbol{\hat{\tau}_{r}}\sim\mathcal{N}(\boldsymbol{\tau},\boldsymbol{\Sigma}_{r})$
(see e.g. Li and Ding,, 2017), where $\boldsymbol{\tau}$ is the vector of true
causal estimates. Our estimator is constructed as

$\boldsymbol{\kappa}_{2}=\boldsymbol{\hat{\tau}_{o}}+\left(\boldsymbol{I}_{K}-\frac{\text{tr}(\boldsymbol{\Sigma}_{r}^{2})\boldsymbol{\Sigma}_{r}}{(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{2}(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})}\right)\left(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}}\right)\,.$

We can compute,
$\mathbb{E}_{r}\left(\mathcal{L}(\boldsymbol{\tau},\boldsymbol{\kappa}_{2})\right)$,
the risk of $\boldsymbol{\kappa}_{2}$, using results from Strawderman, (2003).
This yields the expression

$\displaystyle\mathcal{R}(\boldsymbol{\kappa}_{2})$
$\displaystyle=\frac{1}{K}\left(\text{tr}(\boldsymbol{\Sigma}_{r})+\text{tr}(\boldsymbol{\Sigma}_{r}^{2})\mathbb{E}_{r}\left(\frac{4(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{4}}(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})}{((\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{2}}(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}}))^{2}}-\frac{\text{tr}(\boldsymbol{\Sigma}_{r}^{2})}{(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{2}}(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})}\right)\right)\,,$
(1)

where the expectation is with respect to $\boldsymbol{\hat{\tau}_{r}}$ only,
and $\boldsymbol{\hat{\tau}_{o}}$ is treated as a constant vector. As shown in
Rosenman et al., (2020), the estimator is guaranteed to dominate
$\boldsymbol{\hat{\tau}_{r}}$ under the squared-error loss
$\mathcal{L}(\boldsymbol{\tau},\boldsymbol{\kappa}_{2})=||\boldsymbol{\kappa}_{2}-\boldsymbol{\tau}||^{2}$,
as long as the condition

$4\max_{k}\sigma_{rk}^{4}<\sum_{k}\sigma_{rk}^{4}\,$ (2)

is satisfied.

We will also frequently consider the positive part analogue of
$\boldsymbol{\kappa}_{2}$,

$\boldsymbol{\kappa}_{2+}=\boldsymbol{\hat{\tau}_{o}}+\left(\boldsymbol{I}_{K}-\frac{\text{tr}(\boldsymbol{\Sigma}_{r}^{2})\boldsymbol{\Sigma}_{r}}{(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{2}(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})}\right)_{+}\left(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}}\right)\,,$

which constrains the shrinkage estimator from “over-adjusting” by not allowing
it to put a negative coefficient on $\boldsymbol{\hat{\tau}_{o}}$. The risk of
$\boldsymbol{\kappa}_{2+}$ is guaranteed to be strictly lower than that of
$\boldsymbol{\kappa}_{2}$. In simulations in Rosenman et al., (2020),
$\boldsymbol{\kappa}_{2+}$ routinely dominated $\boldsymbol{\hat{\tau}_{r}}$
even when Condition 2 was not met.

###  3.2 Exact Risk Calculation Under Known Parameters

Our goal is to optimize the experimental design over the risk given in
Expression 1. We suppose, for a moment, that the quantities
$\boldsymbol{\Sigma}_{r}$ and
$\boldsymbol{\xi}=\mathbb{E}_{r}(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})$
are known. Expression 1 can then be rewritten as

$\mathcal{R}(\boldsymbol{\kappa}_{2})=\frac{1}{K}\left(\text{tr}(\boldsymbol{\Sigma}_{r})+\text{tr}(\boldsymbol{\Sigma}_{r}^{2})\mathbb{E}_{r}\left(\frac{4\cdot\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{5}}\boldsymbol{\nu}}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{3}}\boldsymbol{\nu})^{2}}-\frac{\text{tr}(\boldsymbol{\Sigma}_{r}^{2})}{\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{3}}\boldsymbol{\nu}}\right)\right)\,,$

where
$\boldsymbol{\nu}\sim\mathcal{N}(\boldsymbol{\Sigma}_{r}^{-1/2}\boldsymbol{\xi},\boldsymbol{\Sigma}_{r})$.

The expectation is a difference of first moments of ratios of quadratic forms
in normal random variables. Exact expressions for these integrals can be found
in Bao and Kan, (2013). In particular,

$\displaystyle\mathbb{E}_{r}\left(\frac{\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{5}}\boldsymbol{\nu}}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{3}}\boldsymbol{\nu})^{2}}\right)$
$\displaystyle=\int_{0}^{\infty}\det(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r}^{3})^{-1/2}\cdot\exp\left(\frac{1}{2}\left(\boldsymbol{\xi}^{\mathsf{T}}(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r}^{3})^{-1}\boldsymbol{\xi}-\boldsymbol{\xi}^{\mathsf{T}}\boldsymbol{\xi}\right)\right)\cdot$
(3) $\displaystyle\hskip
31.29802pt\left(\text{tr}(\boldsymbol{R})+(\boldsymbol{L}\boldsymbol{\Sigma}_{r}^{-1/2}\boldsymbol{\xi})^{\mathsf{T}}\boldsymbol{R}(\boldsymbol{L}\boldsymbol{\Sigma}_{r}^{-1/2}\boldsymbol{\xi})\right)tdt$
$\displaystyle\mathbb{E}_{r}\left(\frac{1}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{3}}\boldsymbol{\nu})}\right)$
$\displaystyle=\int_{0}^{\infty}\det(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r}^{3})^{-1/2}\cdot\exp\left(\frac{1}{2}\left(\boldsymbol{\xi}^{\mathsf{T}}(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r}^{3})^{-1}\boldsymbol{\xi}-\boldsymbol{\xi}^{\mathsf{T}}\boldsymbol{\xi}\right)\right)dt$

where

$\displaystyle\boldsymbol{L}$
$\displaystyle=(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r}^{3})^{-1/2},\hskip
14.22636pt\text{ and }$ $\displaystyle\boldsymbol{R}$
$\displaystyle=\boldsymbol{L}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{5}\boldsymbol{L}.$

The integrals in Expression 8 can be computed via numerical integration,
yielding an efficient evaluation of the risk for each possible choice of the
parameter values.

###  3.3 Design Options

In this section, we consider several heuristics for designing the experiment
in the absence of perfect knowledge of $\boldsymbol{\Sigma}_{r}$ and
$\boldsymbol{\xi}$.

####  3.3.1 Neyman Allocation

Because the observational data is not randomized, we will need to conduct some
form of statistical adjustment to remove confounding bias in
$\boldsymbol{\hat{\tau}_{o}}$. In this manuscript, we will use stabilized
inverse probability of treatment weighting (SIPW) as our adjustment method.
Briefly, this involves estimating the propensity score – the probability of
treatment in the observational study – as a function of the observed
covariates, and scaling the observed outcomes by the inverse of each unit’s
estimated propensity score. For more details, see Imbens and Rubin, (2015).

Once the adjustment is made, a simple approach to the experimental design
problem is to assume that there is no residual unmeasured confounding in the
observational study. Importantly, we can make this assumption for the purposes
of _design_ only. The unconfoundedness assumption need not be strictly true to
ensure good performance of $\boldsymbol{\kappa}_{2}$ once our experiment is
completed. The shrinkage properties of $\boldsymbol{\kappa}_{2}$ ensure that
its risk will be lower than $\boldsymbol{\hat{\tau}_{r}}$ as long as Condition
2 is met, irrespective of the presence of residual confounding in the
observational study. Hence, we retain the implicit guarantee against a risk
increase when it comes to _inference_ , even if this assumption turns out to
be incorrect.

If we are willing to suppose that Assumption 3 holds, as well as
unconfoundedness, we can use the observational study to unbiasedly estimate
the mean and variance of the potential outcomes in each stratum. Denote the
estimates obtained from the observational data as

$\displaystyle\hat{\sigma}_{kt}^{2}$
$\displaystyle=\widehat{\mathrm{var}}\left(Y(1)\mid S=k\right)\hskip
14.22636pt\text{ and }\hskip
14.22636pt\hat{\sigma}_{kc}^{2}=\widehat{\mathrm{var}}\left(Y(0)\mid
S=k\right)\,.$

The simplest design heuristic is to then use a Neyman allocation (Splawa-
Neyman et al.,, 1990) without a cost constraint, e.g.

$\displaystyle
n_{rkt}=\frac{n_{r}\cdot\hat{\sigma}_{kt}^{2}}{\sum_{k}\hat{\sigma}_{kt}^{2}+\hat{\sigma}_{kc}^{2}}\
\hskip 14.22636pt\text{ and }\hskip 14.22636ptn_{rkc}$
$\displaystyle=\frac{n_{r}\cdot\hat{\sigma}_{kc}^{2}}{\sum_{k}\hat{\sigma}_{kt}^{2}+\hat{\sigma}_{kc}^{2}}\,.$

Such a design would be optimal if the risk of our estimator were only
$\text{tr}(\boldsymbol{\Sigma}_{r})/K$, the first term in Expression 1. Though
the design does not directly optimize over the shrinkage portion of the risk,
it serves as a reasonable starting point for the purposes of design. As we
will see in Section 4, it also typically yields good performance for
$\boldsymbol{\kappa}_{2}$ in simulations.

####  3.3.2 Heuristic Optimization Assuming $\boldsymbol{\xi}\boldsymbol{=0}$

Per the discussion in Section 3.2, we can compute the risk exactly if both
$\boldsymbol{\Sigma}_{r}$ and
$\boldsymbol{\xi}=\mathbb{E}_{r}(\boldsymbol{\hat{\tau}_{r}})-\boldsymbol{\hat{\tau}_{o}}$
are known. Under Assumption 3 and unconfoundedness, $\boldsymbol{\Sigma}_{r}$
can be estimated unbiasedly for any choice of
$\boldsymbol{d}=\\{(n_{rkt},n_{rkc})\\}_{k}$. However, $\boldsymbol{\xi}$ may
be nonzero even if unconfoundedness holds, because we consider
$\boldsymbol{\hat{\tau}_{o}}$ to be a fixed draw from the observational
distribution, rather than a random variable.

In this section, we make the additional assumption that $\boldsymbol{\xi}=0$.
Then, we have all the necessary parameter estimates to optimize
$\mathcal{R}(\boldsymbol{\kappa}_{2})$ over the choice of
$\boldsymbol{d}=\\{(n_{rkt},n_{rkc})\\}_{k}$. This problem is directly encoded
in Optimization Problem 4:

minimize $\displaystyle\hskip 14.22636pt\mathcal{R}(\boldsymbol{\kappa}_{2})$
(4) subject to $\displaystyle\hskip
14.22636pt\sigma_{rk}^{2}=\frac{\widehat{\mathrm{var}}_{r}\left(Y(1)\mid
S=k\right)}{n_{rkt}}+\frac{\hat{\sigma}_{kc}^{2}}{n_{rkc}},\hskip
5.69054ptk=1,\dots,K,$ $\displaystyle\hskip
14.22636pt0<n_{rkt},n_{rkc},,\hskip 5.69054ptk=1,\dots,K,$
$\displaystyle\hskip 14.22636ptn_{r}=\sum_{k}n_{rkt}+n_{rkc}\,,$

Unfortunately, $\mathcal{R}(\boldsymbol{\kappa}_{2})$ is not convex in $d$.
However, Optimization Problem 4 can be approximately solved using a greedy
algorithm. Define

$\boldsymbol{d}_{j}=\\{(n_{rkt}^{(j)},n_{rkc}^{(j)}\\}_{k}$

as the allocation of RCT units to strata and treatment level at iteration $j$
of the algorithm. Next, define

$D_{j}=\\{\boldsymbol{d}\mid\text{ $\boldsymbol{d}$ swaps exactly one unit
across strata or treatment level from }\boldsymbol{d}_{j}\\}\,.$

Because there are $K$ strata and two treatment levels, the “swap set” $D_{j}$
will typically contain $2K\times(2K-1)$ possible allocations. This cardinality
can be reduced by imposing additional constraints, as will be discussed in
Section 3.4.

Define $R(\boldsymbol{d},\boldsymbol{V},\boldsymbol{\xi})$ as the value of
$\mathcal{R}(\boldsymbol{\kappa}_{2})$ evaluated under the design
$\boldsymbol{d}$ with estimated stratum potential outcome variances
$\boldsymbol{V}$ and error vectors $\boldsymbol{\xi}$. We will evaluate
$R(\boldsymbol{d},\boldsymbol{V},\boldsymbol{\xi})$ under estimated variances
$\boldsymbol{V}^{\star}=\left(\hat{\sigma}_{kt}^{2},\hat{\sigma}_{kc}^{2}\right)_{k=1}^{K}$
and $\boldsymbol{\xi}=0$.

Now, we can approximately solve Optimization Problem 4 using Algorithm 5:

Start with design $\boldsymbol{d}_{0}=\\{(n_{rkt}^{(0)},n_{rkc}^{(0)}\\}_{k}$.
(5) For iteration $j=1,2,\dots$: For each design $\boldsymbol{d}$ in
$D_{j-1}$: Compute $R(\boldsymbol{d},V^{\star},0)$. Set
$\boldsymbol{d}_{j}=\operatorname*{arg\,min}_{\boldsymbol{d}\in
D_{j-1}}R(\boldsymbol{d})$ If
$R(\boldsymbol{d}_{j},V^{\star},0)>=R(\boldsymbol{d}_{j-1},V^{\star},0)$
Return $\boldsymbol{d}_{j-1}.$

In words, Algorithm 5 will continue to swap units between strata and treatment
levels until no swap will further reduce the estimated risk of the shrinkage
estimator. The algorithm naturally enforces the sample size constraint and
ensures that the returned values will be integers.

Because the objective is non-convex, it is plausible that Algorithm 5 could
get stuck at a local optimum. Practically, we recommend running the algorithm
at a few different starting points (e.g. an equally allocated design, the
Neyman allocation, and several randomly chosen designs), and choose the design
which achieves the minimum final value of the risk. While this approach is not
guaranteed to find an optimum, it will nonetheless find a point with a reduced
value of the objective function. If the relevant assumptions hold, this point
will also be guaranteed to improve efficiency relative to a Neyman allocation.

####  3.3.3 Heuristic Optimization Assuming Worst-Case Error Under
$\boldsymbol{\Gamma}$-Level Unmeasured Confounding

The assumption that $\boldsymbol{\xi}=0$ is fundamentally optimistic: it is
unlikely to hold even in the absence of unmeasured confounding. If there are
unmeasured confounders, it may be far from the truth. We can take a more
defensive approach by imposing a sensitivity model on the observational study,
and optimizing under the worst-case choice of $\boldsymbol{\xi}$.

As in Rosenman and Owen, (2021), we constrain the magnitude of the unmeasured
confounding by imposing the marginal sensitivity model of Tan, (2006). Under
this model, a key odds ratio – between the treatment probability conditional
on the potential outcomes and covariates and the treatment probability
conditional on covariates only – is bounded between $1/\Gamma$ and $\Gamma$,
for a user-chosen parameter $\Gamma\geq 1$. The Tan model can be seen as
extending the popular Rosenbaum sensitivity model (Rosenbaum,, 1987) to the
setting of inverse probability weighting. In practice, $\Gamma$ can be chosen
by computing the treatments odds ratio with and without conditioning on each
of the _measured_ covariates, and choosing the maximum.

Under Assumption 3, we have

$\displaystyle\mu_{kt}$ $\displaystyle\equiv\mathbb{E}_{o}(Y(1)\mid
S=k)=\mathbb{E}_{R}(Y(1)\mid S=k),\hskip 14.22636pt\text{ and }$
$\displaystyle\mu_{kc}$ $\displaystyle\equiv\mathbb{E}_{o}(Y(0)\mid
S=k)=\mathbb{E}_{R}(Y(0)\mid S=k),$

for $k=1,\dots,K$. For any choice of Type I error bound $\alpha\in(0,1)$, we
can use the method of Zhao et al., (2019) to obtain intervals
$\left(\ell_{kt}^{(\Gamma,\alpha)},u_{kt}^{(\Gamma,\alpha)}\right)$ and
$\left(\ell_{kc}^{(\Gamma,\alpha)},u_{kc}^{(\Gamma,\alpha)}\right)$ such that
parameters $\mu_{kt}$ and $\mu_{kc}$ reside within the intervals with at least
$1-\alpha$ probability as long as the true confounding structure lies within
the sensitivity model parameterized by $\Gamma$. The method relies on convex
optimization and the bootstrap in order to generate valid confidence sets.

Because our outcome is binary, we can then use the confidence sets on
$\mu_{kt}$ and $\mu_{kc}$ to obtain valid confidence sets for the stratum
potential outcome variances $\sigma_{kt}^{2}$ and $\sigma_{kc}^{2}$ via the
relations

$\sigma_{kt}^{2}=\mu_{kt}\cdot(1-\mu_{kt})\hskip 14.22636pt\text{ and }\hskip
14.22636pt\sigma_{kc}^{2}=\mu_{kc}\cdot(1-\mu_{kc})\,.$

For a full justification, see Rosenman and Owen, (2021).

Putting these ideas together: under our calibrated choice of $\Gamma$ and a
reasonable choice of $\alpha$, we set

$\displaystyle\xi_{k}$
$\displaystyle=\max\left(\left|u_{kt}^{(\Gamma,\alpha)}-\ell_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|,\left|\ell_{kt}^{(\Gamma,\alpha)}-u_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|\right),$

the worst-case value of the error under our sensitivity model. We collect
these quantities into a vector $\boldsymbol{\xi}^{\prime}$. Next, we collect
the corresponding values of the variances, e.g.

$\hat{\sigma}_{kt}^{2}=\left\\{\begin{array}[]{cc}u_{kt}^{(\Gamma,\alpha)}\cdot(1-u_{kt}^{(\Gamma,\alpha)})&\text{if
}\left|u_{kt}^{(\Gamma,\alpha)}-\ell_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|>\left|\ell_{kt}^{(\Gamma,\alpha)}-u_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|\\\
\ell_{kt}^{(\Gamma,\alpha)}\cdot(1-\ell_{kt}^{(\Gamma,\alpha)})&\text{otherwise}\end{array}\right.$

and

$\hat{\sigma}_{kc}^{2}=\left\\{\begin{array}[]{cc}\ell_{kc}^{(\Gamma,\alpha)}\cdot(1-\ell_{kc}^{(\Gamma,\alpha)})&\text{if
}\left|u_{kt}^{(\Gamma,\alpha)}-\ell_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|>\left|\ell_{kt}^{(\Gamma,\alpha)}-u_{kc}^{(\Gamma,\alpha)}-\hat{\tau}_{ok}\right|\\\
u_{kc}^{(\Gamma,\alpha)}\cdot(1-u_{kc}^{(\Gamma,\alpha)})&\text{otherwise}\end{array}\right.$

into a matrix $\boldsymbol{V}^{\prime}$.

Finally, we can evaluate our function
$R(\boldsymbol{d},\boldsymbol{V}^{\prime},\boldsymbol{\xi}^{\prime})$ to
obtain the risk of $\boldsymbol{\kappa}_{2}$ for any experimental design
$\boldsymbol{d}$ under these parameters. The procedure is henceforth exactly
analogous to the one used in the prior section: we run Algorithm 5,
substituting
$R(\boldsymbol{d},\boldsymbol{V}^{\prime},\boldsymbol{\xi}^{\prime})$ for
$R(\boldsymbol{d},\boldsymbol{V},0)$, and obtain the design that yields the
lowest value of the risk.

###  3.4 Imposing Guardrails on Designs

Algorithm 5 can easily incorporate further restrictions on the set of possible
designs. The simplest method is to set the risk computation for any invalid
design equal to infinity. This will force the algorithm not to choose that
design on the next iteration.

Practically speaking, we suggest incorporating three constraints on the set of
possible designs. First, we suggest imposing a _minimum sample size
constraint_ such that the allocations to any stratum and treatment group
cannot be lower than some value $\text{SS}_{\text{min}}$. This constraint
serves two purposes. First, the risk expression in Equation 1 is dependent on
the normality of $\boldsymbol{\hat{\tau}_{r}}$. In simulations, we find that
the risk expression is quite robust to modest deviations from normality, but
sample sizes could plausibly become so small that a Central Limit Theorem need
not hold even approximately. Second, this constraint naturally helps improve
detachability, since it prevents the variance of any entry of
$\boldsymbol{\hat{\tau}_{r}}$ from growing too large.

Second, we suggest imposing an explicit _detachability_ constraint. Inherent
in this constraint is the idea that an analyst wants to ensure that the
experiment can be analyzed on its own in case stakeholders ultimately decide
not to use the observational data at all. One way to impose this constraint is
to assume the observational study point estimates of the strata potential
outcome variances are correct, and to choose some baseline design
$\boldsymbol{\tilde{d}}=\\{\tilde{n}_{rkt},\tilde{n}_{rkc}\\}_{k}$, e.g. equal
allocation or Neyman allocation. We then invalidate any design
$\boldsymbol{d^{\prime}}=\\{n_{rkt}^{\prime},n_{rkc}^{\prime}\\}_{k}$ for
which

$\sum_{k}\frac{\hat{\sigma}_{kt}^{2}}{n_{rkt}^{\prime}}+\frac{\hat{\sigma}_{kc}^{2}}{n_{rkc}^{\prime}}\geq\delta_{d}\sum_{k}\frac{\hat{\sigma}_{kt}^{2}}{\tilde{n}_{rkt}}+\frac{\hat{\sigma}_{kc}^{2}}{\tilde{n}_{rkc}}\,,$

where $\delta_{d}\geq 1$ is some user-chosen tolerance parameter. In words,
this approach invalidates any design $\boldsymbol{d^{\prime}}$ for which the
estimated risk of $\boldsymbol{\hat{\tau}_{r}}$ under design
$\boldsymbol{d^{\prime}}$ is larger than the estimated risk under the default
design by a factor greater than $\delta_{d}$.

A more robust version of the detachability constraint can be imposed if the
analyst assumes the Tan sensitivity model, as discussed in Section 3.3.3.
Under a given choice of $\Gamma$ and $\alpha$, we can obtain bounds on the
potential outcome means in each stratum. We reject any design such that

$\max_{\mu_{kc}\in\left(\ell_{kc}^{(\Gamma,\alpha)}u_{kc}^{(\Gamma,\alpha)}\right),\atop\mu_{kt}\in\left(\ell_{kt}^{(\Gamma,\alpha)},u_{kt}^{(\Gamma,\alpha)}\right)}\frac{\sum_{k}\frac{\mu_{kt}(1-\mu_{kt})}{n_{rkt}^{\prime}}+\frac{\mu_{kt}(1-\mu_{kt})}{\tilde{n}_{rkc}}}{\sum_{k}\frac{\mu_{kt}(1-\mu_{kt})}{n_{rkt}^{\prime}}+\frac{\mu_{kc}(1-\mu_{kc})}{\tilde{n}_{rkc}}}\geq\delta_{d}\,.$
(6)

In words, this means that we are rejecting any design
$\boldsymbol{d^{\prime}}$ such that the risk of $\boldsymbol{\hat{\tau}_{r}}$
under design $d^{\prime}$ is larger than the estimated risk under the default
design by a factor greater than $\delta_{d}$, for _any_ configuration of the
potential outcome means consistent with our sensitivity bounds. Practically
the left-hand-side of Inequality 6 can be reduced to a quadratic fractional
programming problem and solved via Dinkelbach’s method (Dinkelbach,, 1967).

Third, we suggest imposing a _risk reduction_ constraint. This means that we
enforce Condition 2 explicitly by not allowing the optimization algorithm to
choose allocations under which it fails to hold. Like the detachability
constraint, we can impose either a point estimate version or a robust version.
The point estimate version invalidates any design design
$\boldsymbol{d^{\prime}}=\\{n_{rkt}^{\prime},n_{rkc}^{\prime}\\}_{k}$ for
which

$\displaystyle 4\max_{k}$
$\displaystyle\left(\frac{\hat{\sigma}_{kt}^{2}}{n_{rkt}^{\prime}}+\frac{\hat{\sigma}_{kc}^{2}}{n_{rkc}^{\prime}}\right)^{2}>\sum_{k}\left(\frac{\hat{\sigma}_{kt}^{2}}{n_{rkt}^{\prime}}+\frac{\hat{\sigma}_{kc}^{2}}{n_{rkc}^{\prime}}\right)^{2}.$

The robust version incorporates the parameter bounds from the Tan sensitivity
model. We reject any design such that

$\displaystyle 4\max_{k}$
$\displaystyle\min_{\mu_{kc}\in\left(\ell_{kc}^{(\Gamma,\alpha)}u_{kc}^{(\Gamma,\alpha)}\right),\atop\mu_{kt}\in\left(\ell_{kt}^{(\Gamma,\alpha)},u_{kt}^{(\Gamma,\alpha)}\right)}\left(\frac{\mu_{kt}(1-\mu_{kt})}{n_{rkt}^{\prime}}+\frac{\mu_{kc}(1-\mu_{kc})}{n_{rkc}^{\prime}}\right)^{2}>$
$\displaystyle\max_{\mu_{kc}\in\left(\ell_{kc}^{(\Gamma,\alpha)}u_{kc}^{(\Gamma,\alpha)}\right),\atop\mu_{kt}\in\left(\ell_{kt}^{(\Gamma,\alpha)},u_{kt}^{(\Gamma,\alpha)}\right)}\sum_{k}\left(\frac{\mu_{kt}(1-\mu_{kt})}{n_{rkt}^{\prime}}+\frac{\mu_{kc}(1-\mu_{kc})}{n_{rkc}^{\prime}}\right)^{2}.$

##  4 Simulation Study

###  4.1 Simulation Set-Up

We pattern the simulations on those in Rosenman et al., (2018), considering a
situation with a relatively rare, binary outcome and a modest effect size.

We first generate an observational super-population of $1\times 10^{6}$ units,
and an experimental super-population of the same size. We suppose that the
completed observational study comprises $20,000$ units, sampled a single time
from the corresponding super-population. The prospective RCT comprises $1,000$
units, which will be sampled repeatedly from the experimental super-
population.

We define $j\in\mathcal{O}$ as the indexing variable for the observational
super-population and $i\in\mathcal{E}$ as the indexing variable for the
experimental super-population. We use $\ell$ as an index over both
populations. For each unit $\ell\in\mathcal{O}\cup\mathcal{E}$, we suppose
there is a covariate vector $\boldsymbol{X}_{\ell}\in\mathbb{R}^{5}$ where
$\boldsymbol{X}_{\ell}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathcal{N}\left(0,\Sigma\right)$
for $\Sigma$ such that each covariate has unit variance and roughly a quarter
of the covariances are $+0.1$, roughly a quarter are $-0.1$, and the remainder
are $0$. Such a covariance structure is roughly consistent with the applied
data analysis from the Women’s Health Initiative, as used in Rosenman et al.,
(2018).

The untreated potential outcomes $Y_{\ell}(0)$ are sampled as independent
Bernoulli random variables with

$\text{Pr}(Y_{\ell}(0)=1\mid{x}_{\ell})=\frac{1}{1+e^{-\alpha-\beta^{\mathsf{T}}{x}_{\ell}+\varepsilon_{\ell}}},\quad\text{for
$\beta=(1,1,1,1,1)^{\mathsf{T}}$}$

where $\varepsilon_{\ell}$ are generated as IID $\mathcal{N}(0,1)$ random
variables and $\alpha$ is chosen such that the average incidence rate is 10%.
The treatment variables in the observational study are independent Bernoulli
random variables with

$\text{Pr}(W_{j}=1\mid{x}_{j})=\frac{1}{1+e^{-\gamma^{\mathsf{T}}{x}_{j}}},\quad\text{for
$\gamma=(\sqrt{2},\sqrt{2},\sqrt{2},0,0)^{\mathsf{T}}$.}$

The nonzero inner product between $\beta$ and $\gamma$ induces strong
selection bias in the observational study, with treated units likelier to have
larger untreated potential outcomes.

For both datasets, we suppose the data is split into twelve strata based on
the first and second covariates. The strata boundaries are defined by the 25th
and 50th quantiles of the first covariate and by the quartiles of the second
covariate. We seek to obtain estimates of the average causal effect in each of
the resultant strata.

In both datasets, we assign individual treatment effects to match three
different structures. We number the strata from $k=1,\dots,12$ such that
stratum $1$ corresponds to the lowest stratum on both covariates, stratum 2
corresponds to the lowest stratum on the first covariate and the second
quartile of the second covariate, stratum $3$ corresponds to the lowest
stratum on the first covariate and the third quartile of the second covariate,
etc. The values of $\tau_{k}$, the stratum average treatment effects, in the
constant, linear, and quadratic strata-level treatment effect models are

$\displaystyle\tau_{k}=T,\quad\tau_{k}=-T\times\frac{k}{K},\quad\text{and}\quad\tau_{k}=T\times\left(\frac{k}{K}\right)^{2}$
(7)

respectively. In each case we choose the scale $T>0$ so that Cohen’s $d$
(Cohen,, 1988) in the observational study precisely equals $0.5$, which Cohen
calls a medium effect size.

Denote as $n_{k}$ the total number of super-population units that fall into a
given stratum. Given the $\tau_{k}$, the treated potential outcomes are
assigned by randomly selecting $\tau_{k}\times n_{k}$ units for which
$Y_{\ell}(0)=0$, and setting $Y_{\ell}(1)=1$ for those units. For all of the
remaining units, $Y_{\ell}(0)=Y_{\ell}(1)$. Because we use the same process
across both super-populations, we enforce Assumption 3, such that the stratum-
specific causal effects (as well as potential outcome means and variances) are
shared under the observational and RCT data distributions.

The observational data is sampled a single time. Next, leveraging the
observational sample, we compute the allocations of units to strata in the RCT
(the RCT “design”) using each of the methods discussed in the prior section.
For each possible design, we sample the RCT units from the super-population
$5,000$ times. For each iteration, for each stratum $k$, we assume treatment
is assigned via a simple random sample of $n_{rkt}$ units out of the
$n_{rkt}+n_{rkc}$ units recruited for the stratum. Once the units are drawn
and treatments are assigned, we compute the estimators
$\boldsymbol{\hat{\tau}_{r}}$, $\boldsymbol{\kappa}_{2}$, and
$\boldsymbol{\kappa}_{2+}$. We compute the $L_{2}$ distance between each
estimate and the true treatment effects $\tau$, and take the average over all
$5,000$ simulations.

###  4.2 Ideal Case: No Unmeasured Confounding

We begin with the simplest case: we suppose all of the covariates are measured
in the observational study, so there is no residual unmeasured confounding.
This is an idealized case, in which all of the selection bias in the
observational study can be removed with a statistical adjustment. We fit a
propensity score to the observational study data and use stabilized inverse
probability of treatment weighting (SIPW) to compute the observational causal
estimates.

In Table 1, we show the average $L_{2}$ errors over $5,000$ simulations. We
consider the equal allocation and Neyman allocation designs, as well as the
“naive” design assuming $\boldsymbol{\xi}=0$. We also consider the defensive
approach discussed in Section 3.3.3, and compute the optimal design under
errors computed with the Tan sensitivity model parameter set to 1.0, 1.1, 1.2,
and 1.5. Lastly, we consider an “oracle” design, in which we run Algorithm 5
but provide it with the true values of the potential outcome variances and the
error in the observational study. Results are given for the three different
treatment effect models: constant (c), linear ($\ell$), and quadratic (q).
Risk estimates are expressed as percentages of the risk of
$\boldsymbol{\hat{\tau}_{r}}$ when using an equal allocation for the given
treatment effect model.

|  |  |  |  |  Max Bias, $\boldsymbol{\Gamma}$ Value |
---|---|---|---|---|---|---
Est. | Trt. | Equal | Neyman | Naïve | 1.0 | 1.1 | 1.2 | 1.5 | Oracle
$\boldsymbol{\hat{\tau}_{r}}$ | c | 100% | 87% | 91% | 100% | 96% | 94% | 94% | 96%
$\boldsymbol{\kappa}_{2}$ | 82% | 48% | 44% | 52% | 48% | 47% | 50% | 42%
$\boldsymbol{\kappa}_{2+}$ | 38% | 28% | 26% | 26% | 26% | 26% | 28% | 23%
$\boldsymbol{\hat{\tau}_{r}}$ | $\ell$ | 100% | 89% | 92% | 95% | 94% | 95% | 97% | 104%
$\boldsymbol{\kappa}_{2}$ | 93% | 66% | 58% | 58% | 57% | 60% | 64% | 50%
$\boldsymbol{\kappa}_{2+}$ | 59% | 51% | 45% | 43% | 45% | 47% | 49% | 33%
$\boldsymbol{\hat{\tau}_{r}}$ | q | 100% | 86% | 91% | 95% | 98% | 94% | 92% | 91%
$\boldsymbol{\kappa}_{2+}$ | 81% | 47% | 45% | 52% | 52% | 50% | 48% | 41%
$\boldsymbol{\kappa}_{2+}$ | 37% | 29% | 27% | 28% | 28% | 30% | 29% | 25%
Table 1:  Risk over $5,000$ iterations of
$\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\kappa}_{2}$, and
$\boldsymbol{\kappa}_{2+}$ under various experimental designs, in the case of
no unmeasured confounding in the observational study. Risks are expressed as a
percentage of the risk of $\boldsymbol{\hat{\tau}_{r}}$ using an equally
allocated experiment, for each of the three treatment effect models. The
minimum non-oracle risk in each row is denoted with an underline.

Across the three treatment effect models, we see that results are relatively
consistent in the ordering of the estimators. The Neyman allocation always
performs best when using $\boldsymbol{\hat{\tau}_{r}}$ alone, realizing an
8$-$14% error reduction relative to the equal allocation. The other
allocations typically achieve more modest error reductions when using
$\boldsymbol{\hat{\tau}_{r}}$.

If we stick with the Neyman allocation but switch to using either the
shrinkage estimator $\boldsymbol{\kappa}_{2}$ or its positive part analogue
$\boldsymbol{\kappa}_{2+}$, we can realize massive additional error reductions
on the order of 45 to 60%, depending on the treatment effect model. However,
the naïve allocation is typically even better, allowing us to realize
improvements of a few additional percentage points in error. The naïve
allocation yields the best performance with $\boldsymbol{\kappa}_{2}$ and
$\boldsymbol{\kappa}_{2+}$ in the linear and quadratic treatment effect
models, with the robust allocations typically doing nearly as well. In the
linear treatment effect model, we actually do best using a $\Gamma=1.1$ robust
design for $\boldsymbol{\kappa}_{2}$ and a $\Gamma=1.0$ robust design for
$\boldsymbol{\kappa}_{2+}$, with the naïve allocation performing nearly as
well.

###  4.3 Practical Case: Unmeasured Confounding

The assumption of unconfoundedness is not defensible in many practical
settings, in which a relatively sparse set of covariates are measured in the
observational study. Moreover, the assumption is not testable. Hence, we run a
second set of simulations in which we induce confounding by assuming that the
third entry in ${x}_{i}$ is not measured. This third covariate affects both
treatment probabilities and outcomes, so bias in the observational study
cannot be fully corrected with a propensity score adjustment.

Results are given in Table 2. Results differ slightly from the simulations
assuming no unmeasured confounding. When using $\boldsymbol{\kappa}_{2}$ and
$\boldsymbol{\kappa}_{2+}$, the best non-oracle designs are the robust designs
assuming $\Gamma$-level unmeasured confounding for
$\Gamma\in\\{1.0,1.1,1.2\\}$. The design under $\Gamma=1$ does best in the
constant and linear case, while the design under $\Gamma=1.2$ does best in the
quadratic case. The naïve allocation still performs quite well, with risk
typically only a few percentage points worse than the best performer.

|  |  |  |  |  Max Bias, $\boldsymbol{\Gamma}$ Value |
---|---|---|---|---|---|---
Est. | Trt. | Equal | Neyman | Naïve | 1.0 | 1.1 | 1.2 | 1.5 | Oracle
$\boldsymbol{\hat{\tau}_{r}}$ | c | 100% | 90% | 90% | 90% | 92% | 93% | 95% | 102%
$\boldsymbol{\kappa}_{2}$ | 102% | 81% | 74% | 72% | 72% | 72% | 77% | 69%
$\boldsymbol{\kappa}_{2}$ | 96% | 80% | 74% | 71% | 72% | 72% | 76% | 67%
$\boldsymbol{\hat{\tau}_{r}}$ | $\ell$ | 100% | 93% | 93% | 94% | 95% | 96% | 96% | 104%
$\boldsymbol{\kappa}_{2}$ | 102% | 85% | 77% | 75% | 76% | 77% | 79% | 73%
$\boldsymbol{\kappa}_{2+}$ | 98% | 84% | 77% | 75% | 76% | 76% | 79% | 71%
$\boldsymbol{\hat{\tau}_{r}}$ | q | 100% | 89% | 90% | 93% | 92% | 91% | 96% | 96%
$\boldsymbol{\kappa}_{2+}$ | 101% | 74% | 69% | 68% | 68% | 67% | 73% | 66%
$\boldsymbol{\kappa}_{2+}$ | 88% | 72% | 67% | 66% | 66% | 65% | 71% | 63%
Table 2:  Risk over $5,000$ iterations of
$\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\kappa}_{2}$, and
$\boldsymbol{\kappa}_{2+}$ under various experimental designs, in the case of
unmeasured confounding in the observational study via failure to measure the
third covariate. Risks are expressed as a percentage of the risk of
$\boldsymbol{\hat{\tau}_{r}}$ using an equally allocated experiment, for each
of the three treatment effect models. The minimum non-oracle risk in each row
is denoted with an underline.

These results point to a few practical guidelines for designing toward
shrinkage. First, we find that the naïve allocation is quite robust, even when
$\boldsymbol{\xi}$ is far from zero. This is evident from the fact that the
naïve allocation always yields significant performance gains over the equal
and Neyman allocations when using a shrinkage estimator, even when unmeasured
confounding is present.

Second, in the presence of unmeasured confounding, one may achieve modest
further improvements from enforcing a robust design at a relatively low value
of $\Gamma$. In this example, we have _not_ constructed the confounding such
that it matches the form of our sensitivity model. Exclusion of the third
covariate induces enormous discrepancies in the treatment odds between the
true and estimated propensity scores in a small proportion of individuals: for
about 0.2% of the super-population units, the difference exceeds a
multiplicative factor of 100. For most of the population, however, the true
and estimated treatment odds differ by a much smaller factor. The Tan model
imposes a worst-case bound on the deviation between the true and estimated
odds of treatment, so no value of $\Gamma$ between 1.0 and 2.0 is large enough
to account for our most extreme deviations.

Nonetheless, the magnitude of the worst-case error under the Tan model
correlates reasonably well with the true values of $\boldsymbol{\xi}$ when
choosing $\Gamma=1.0,1.1,$ or $1.2$, offering one explanation for the strong
performance of the allocations designed under these schemes. In a more general
sense, the robust allocations under $\Gamma$ serve as a form of
regularization, bringing the design closer to an equal allocation to hedge
against the possibility that $\boldsymbol{\hat{\tau}_{o}}$ is far off from
$\boldsymbol{\tau}$. In doing so, the robust designs may improve empirical
performance even when the Tan model does not accurately characterize the form
of unmeasured confounding.

##  5 Discussion

We have considered the problem of designing a stratified experiment when we
plan to use an Empirical Bayes estimator, $\boldsymbol{\kappa}_{2}$, to shrink
the stratum causal estimates toward estimates previously obtained from an
observational study. As we have shown, the risk of $\boldsymbol{\kappa}_{2}$
can be computed explicitly if the stratum potential outcome variances and the
stratum-specific errors of the observational study estimates,
$\boldsymbol{\xi}$, are known. In the absence of such information, we have
proposed three heuristics – Neyman allocation, “naïve” allocation assuming
$\boldsymbol{\xi}=0$, and robust allocation under the worst-case value of
$\boldsymbol{\xi}$ given a sensitivity model – for designing the experiment.
We have also emphasized “detachability,” the ability to do good causal
inference using the RCT on its own, and suggested imposing constraints on the
experimental designs to ensure detachability is preserved.

We simulated a realistic scenario in which we have access to a large
observational database and are interested in a relatively rare outcome. We
considered a stratification yielding twelve strata, in which the outcome
frequency varies dramatically across the strata. In this setting, there were
significant risk improvements to be realized by using a shrinkage estimator.
Additional gains were possible when switching from an equally allocated RCT to
one designed explicitly for use with a shrinkage estimator. In simulations
with and without unmeasured confounding in the observational study, we saw
that the naïve allocation typically performed well, exhibiting surprising
robustness to the case when $\boldsymbol{\xi}$ was far from zero. Our robust
allocations outperformed in the case when unmeasured confounding was present
in the data.

There are many plausible extensions to this line of work. We have assumed a
rigid model for treatment effect heterogeneity in this paper: treatment
effects vary according to a known stratification on observed covariates. More
modern work (Lee et al.,, 2021; Nie and Wager,, 2021) focuses on estimating
heterogeneous treatment effects empirically, allowing for greater flexibility
in how the effects differ across units. Incorporating such ideas into this
work, we might imagine using some of the observational data in a first stage
to estimate a stratification scheme, and using the remaining data for
shrinkage in a second stage. Alternatively, we could consider modeling the
causal effect explicitly as a function of the covariates, e.g. defining
$\tau(X_{i})$ rather than a vector of true treatment effects
$\boldsymbol{\tau}$. We could then define a flexible shrinker to trade off
between estimates not within strata, but within nearby values of the
covariates themselves.

## References

  * Athey et al., (2019) Athey, S., Chetty, R., Imbens, G. W., and Kang, H. (2019).  The surrogate index: Combining short-term proxies to estimate long-term treatment effects more rapidly and precisely.  Technical report, National Bureau of Economic Research.
  * Bai, (2019) Bai, Y. (2019).  Optimality of matched-pair designs in randomized controlled trials.  Available at SSRN 3483834.
  * Bao and Kan, (2013) Bao, Y. and Kan, R. (2013).  On the moments of ratios of quadratic forms in normal random variables.  Journal of Multivariate Analysis, 117:229–245.
  * Bareinboim and Pearl, (2016) Bareinboim, E. and Pearl, J. (2016).  Causal inference and the data-fusion problem.  Proceedings of the National Academy of Sciences, 113(27):7345–7352.
  * Chambaz et al., (2014) Chambaz, A., van der Laan, M. J., and Zheng, W. (2014).  Targeted covariate-adjusted response-adaptive lasso-based randomized controlled trials.  Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, pages 345–368.
  * Cohen, (1988) Cohen, J. (1988).  Statistical power analysis for the behavioral sciences.  Lawrence Erlbaum Associates, Hillsdale, NJ, 2nd edition.
  * Colnet et al., (2020) Colnet, B., Mayer, I., Chen, G., Dieng, A., Li, R., Varoquaux, G., Vert, J.-P., Josse, J., and Yang, S. (2020).  Causal inference methods for combining randomized trials and observational studies: a review.  arXiv preprint arXiv:2011.08047.
  * Dinkelbach, (1967) Dinkelbach, W. (1967).  On nonlinear fractional programming.  Management Science, 13(7):492–498.
  * Gagnon-Bartsch et al., (2021) Gagnon-Bartsch, J. A., Sales, A. C., Wu, E., Botelho, A. F., Erickson, J. A., Miratrix, L. W., and Heffernan, N. T. (2021).  Precise unbiased estimation in randomized experiments using auxiliary observational data.  arXiv preprint arXiv:2105.03529.
  * Green and Strawderman, (1991) Green, E. J. and Strawderman, W. E. (1991).  A James-Stein type estimator for combining unbiased and possibly biased estimators.  Journal of the American Statistical Association, 86(416):1001–1006.
  * Green et al., (2005) Green, E. J., Strawderman, W. E., Amateis, R. L., and Reams, G. A. (2005).  Improved estimation for multiple means with heterogeneous variances.  Forest Science, 51(1):1–6.
  * Hahn et al., (2011) Hahn, J., Hirano, K., and Karlan, D. (2011).  Adaptive experimental design using the propensity score.  Journal of Business & Economic Statistics, 29(1):96–108.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction.  Cambridge University Press, New York, NY, USA.
  * Kallus et al., (2018) Kallus, N., Puli, A. M., and Shalit, U. (2018).  Removing hidden confounding by experimental grounding.  In Advances in Neural Information Processing Systems, pages 10888–10897.
  * Lee et al., (2021) Lee, K., Small, D. S., and Dominici, F. (2021).  Discovering heterogeneous exposure effects using randomization inference in air pollution studies.  Journal of the American Statistical Association, 116(534):569–580.
  * Li and Ding, (2017) Li, X. and Ding, P. (2017).  General forms of finite population central limit theorems with applications to causal inference.  Journal of the American Statistical Association, 112(520):1759–1769.
  * Neyman, (1992) Neyman, J. (1992).  On the two different aspects of the representative method: the method of stratified sampling and the method of purposive selection.  In Breakthroughs in Statistics, pages 123–150. Springer.
  * Nie and Wager, (2021) Nie, X. and Wager, S. (2021).  Quasi-oracle estimation of heterogeneous treatment effects.  Biometrika, 108(2):299–319.
  * Offer-Westort et al., (2021) Offer-Westort, M., Coppock, A., and Green, D. P. (2021).  Adaptive experimental design: Prospects and applications in political science.  American Journal of Political Science.
  * Peysakhovich and Lada, (2016) Peysakhovich, A. and Lada, A. (2016).  Combining observational and experimental data to find heterogeneous treatment effects.  arXiv preprint arXiv:1611.02385.
  * Rosenbaum, (1987) Rosenbaum, P. R. (1987).  Sensitivity analysis for certain permutation inferences in matched observational studies.  Biometrika, 74(1):13–26.
  * Rosenman et al., (2020) Rosenman, E., Basse, G., Owen, A., and Baiocchi, M. (2020).  Combining observational and experimental datasets using shrinkage estimators.  arXiv preprint arXiv:2002.06708.
  * Rosenman et al., (2018) Rosenman, E., Owen, A. B., Baiocchi, M., and Banack, H. (2018).  Propensity score methods for merging observational and experimental datasets.  arXiv preprint arXiv:1804.07863.
  * Rosenman and Owen, (2021) Rosenman, E. T. R. and Owen, A. B. (2021).  Designing experiments informed by observational studies:.  Journal of Causal Inference, 9(1):147–171.
  * Splawa-Neyman et al., (1990) Splawa-Neyman, J., Dabrowska, D. M., and Speed, T. P. (1990).  On the application of probability theory to agricultural experiments. essay on principles. section 9.  Statistical Science, pages 465–472.
  * Strawderman, (2003) Strawderman, W. E. (2003).  On minimax estimation of a normal mean vector for general quadratic loss.  In Mathematical Statistics and Applications: Festschrift for Constance Van Eeden, pages 3–14. Institute of Mathematical Statistics.
  * Tabord-Meehan, (2018) Tabord-Meehan, M. (2018).  Stratification trees for adaptive randomization in randomized controlled trials.  arXiv preprint arXiv:1806.05127.
  * Tan, (2006) Tan, Z. (2006).  A distributional approach for causal inference using propensity scores.  Journal of the American Statistical Association, 101(476):1619–1637.
  * Thompson, (1933) Thompson, W. R. (1933).  On the likelihood that one unknown probability exceeds another in view of the evidence of two samples.  Biometrika, 25(3/4):285–294.
  * Thompson, (1935) Thompson, W. R. (1935).  On the theory of apportionment.  American Journal of Mathematics, 57(2):450–456.
  * Wager and Athey, (2018) Wager, S. and Athey, S. (2018).  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, 113(523):1228–1242.
  * Zhao et al., (2019) Zhao, Q., Small, D. S., and Bhattacharya, B. B. (2019).  Sensitivity analysis for inverse probability weighting estimators via the percentile bootstrap.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(4):735–761.

##  Appendix A Alternatives Shrinkers to $\boldsymbol{\kappa}_{2}$

We have supposed that, after the experiment’s conclusion, we will use
$\boldsymbol{\kappa}_{2}$ or $\boldsymbol{\kappa}_{2+}$ to shrink causal
estimates from $\boldsymbol{\hat{\tau}_{r}}$ toward
$\boldsymbol{\hat{\tau}_{o}}$. However, our results are not dependent on the
form of the shrinkage estimator. Our goal is to apply Theorem 3.1 in
Strawderman, (2003) to compute the exact risk of the estimator. Hence, we can
use any estimator $\boldsymbol{\theta}$ of the form

$\boldsymbol{\theta}=\boldsymbol{\hat{\tau}_{r}}+\boldsymbol{\Sigma}_{r}g_{\theta}(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{o}})$

where $g_{\theta}(x,y)$ is weakly differential in $x$ and
$\mathbb{E}_{r}\left(||g_{\theta}||^{2}\right)<\infty$.

A plausible alternative estimator proposed in Rosenman et al., (2020) is

$\boldsymbol{\kappa}_{1}=\boldsymbol{\hat{\tau}_{r}}+\left(\frac{\text{tr}(\boldsymbol{\Sigma}_{r})}{(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})}\right)\left(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}}\right)\,,$

and its positive part analogue, $\boldsymbol{\kappa}_{1+}$. Two other standard
alternatives, introduced in Green et al., (2005), are

$\boldsymbol{\delta}_{1}=\boldsymbol{\hat{\tau}_{r}}+\left(\frac{K-2}{(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{r}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{-1}\left(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}}\right)}\right)\left(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}}\right)$

and

$\boldsymbol{\delta}_{2}=\boldsymbol{\hat{\tau}_{r}}+\left(\frac{(K-2)\boldsymbol{\Sigma}_{r}^{-1}}{(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{-2}(\boldsymbol{\hat{\tau}_{r}}-\boldsymbol{\hat{\tau}_{o}})}\right)(\boldsymbol{\hat{\tau}_{o}}-\boldsymbol{\hat{\tau}_{r}})\,.$

Each also has a positive part version that is straightforward to define.

Once a suitable estimator $\boldsymbol{\theta}$ is chosen, we can redefine our
procedures to work with that estimator. The first step is to apply Theorem 3.1
from Strawderman, (2003) to compute its risk conditional on
$\boldsymbol{\hat{\tau}_{o}}$,

$\mathcal{R}(\theta)=\frac{1}{K}\left(\text{tr}(\boldsymbol{\Sigma}_{r})+\mathbb{E}_{r}\left(\sum_{k}\sigma_{rk}^{4}\left(g_{\theta,k}^{2}(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{o}})+2\frac{\partial
g_{\theta,k}(\boldsymbol{\hat{\tau}_{r}},\boldsymbol{\hat{\tau}_{o}})}{\partial\tau_{rk}}\right)\right)\right)\,.$

Next, we can deploy the method from Section 3.2 to obtain an exact integral
expression for the the risk of the estimator. Bao and Kan, (2013) contains a
detailed explanation of how to construct each component of the integral. For
example, if we use $\boldsymbol{\kappa}_{1+}$, we obtain

$\mathcal{R}(\boldsymbol{\kappa}_{1+})=\frac{1}{K}\left(\text{tr}(\boldsymbol{\Sigma}_{r})+\text{tr}(\boldsymbol{\Sigma}_{r})\mathbb{E}\left(\frac{4\cdot\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{2}}\boldsymbol{\nu}}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}\boldsymbol{\nu})^{2}}-\frac{\text{tr}(\boldsymbol{\Sigma}_{r})}{\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}\boldsymbol{\nu}}\right)\right)\,,$

and the relevant integrals can be computed as

$\displaystyle\mathbb{E}\left(\frac{\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{{}^{2}}\boldsymbol{\nu}}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}\boldsymbol{\nu})^{2}}\right)$
$\displaystyle=\int_{0}^{\infty}\det(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma})^{-1/2}\cdot\exp\left(\frac{1}{2}\left(\boldsymbol{\xi}^{\mathsf{T}}(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r})^{-1}\boldsymbol{\xi}-\boldsymbol{\xi}^{\mathsf{T}}\boldsymbol{\xi}\right)\right)\cdot$
(8) $\displaystyle\hskip
31.29802pt\left(\text{tr}(\boldsymbol{R})+(\boldsymbol{L}\boldsymbol{\Sigma}_{r}^{-1/2}\boldsymbol{\xi})^{\mathsf{T}}\boldsymbol{R}(\boldsymbol{L}\boldsymbol{\Sigma}_{r}^{-1/2}\boldsymbol{\xi})\right)tdt$
$\displaystyle\mathbb{E}\left(\frac{1}{(\boldsymbol{\nu}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}\boldsymbol{\nu})}\right)$
$\displaystyle=\int_{0}^{\infty}\det(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r})^{-1/2}\cdot\exp\left(\frac{1}{2}\left(\boldsymbol{\xi}^{\mathsf{T}}(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r})^{-1}\boldsymbol{\xi}-\boldsymbol{\xi}^{\mathsf{T}}\boldsymbol{\xi}\right)\right)dt$

where

$\displaystyle\boldsymbol{L}$
$\displaystyle=(\boldsymbol{I}_{K}+2t\boldsymbol{\Sigma}_{r})^{-1/2},\hskip
14.22636pt\text{ and }$ $\displaystyle\boldsymbol{R}$
$\displaystyle=\boldsymbol{L}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{2}\boldsymbol{L}.$

This process can be repeated for any estimator $\boldsymbol{\theta}$. Once
integral expressions for the risk have been obtained, the design heuristics
described in Section 3.3 can be deployed, using
$\mathcal{R}(\boldsymbol{\theta})$ as the objective function rather than
$\mathcal{R}(\boldsymbol{\kappa}_{2})$.