BY 4.0 license Open Access Published by De Gruyter July 15, 2021

Designing experiments informed by observational studies

Evan T. R. Rosenman and Art B. Owen

Abstract

The increasing availability of passively observed data has yielded a growing interest in “data fusion” methods, which involve merging data from observational and experimental sources to draw causal conclusions. Such methods often require a precarious tradeoff between the unknown bias in the observational dataset and the often-large variance in the experimental dataset. We propose an alternative approach, which avoids this tradeoff: rather than using observational data for inference, we use it to design a more efficient experiment. We consider the case of a stratified experiment with a binary outcome and suppose pilot estimates for the stratum potential outcome variances can be obtained from the observational study. We extend existing results to generate confidence sets for these variances, while accounting for the possibility of unmeasured confounding. Then, we pose the experimental design problem as a regret minimization problem subject to the constraints imposed by our confidence sets. We show that this problem can be converted into a concave maximization and solved using conventional methods. Finally, we demonstrate the practical utility of our methods using data from the Women’s Health Initiative.

MSC 2010: 62K05; 90C25

1 Introduction

The past half-century of causal inference research has engendered a healthy skepticism toward observational data [1]. In observational datasets, researchers do not control whether each individual receives a treatment of interest. Hence, they cannot be certain that treated individuals and untreated individuals are otherwise comparable.

This challenge can be overcome only if the covariates measured in the observational data are sufficiently rich to fully explain who receives the treatment and who does not. This is a fundamentally untestable assumption – and even if it holds, careful modeling is necessary to remove the selection effect. The applied literature includes myriad examples of treatments that showed promise in observational studies only to be overturned by later randomized trials [2]. One prominent case, the effect of hormone therapy on the health of postmenopausal women, will be discussed in this manuscript [3].

The “virtuous” counterpart to observational data is the well-designed experiment. Data from a randomized trial can yield unbiased estimates of a causal effect without the need for problematic statistical assumptions. However, experiments are not without their own significant drawbacks. Experiments are frequently expensive, and, as a consequence, often involve fewer units than observational studies. Particularly if one is interested in subgroup causal effects, this means experimental estimates can be imprecise. Moreover, experiments sometimes involve inclusion criteria that can make them dissimilar from target populations of interest. In this way, experiments are often said to have poor “external validity” [4].

In this article, we use the observational data not for inference, but rather to influence the design of an experiment. Our method seeks to retain the possibility of unbiased estimation from the experiment, while also leveraging the ready availability of observational databases to improve the experiment’s efficiency. Because the observational data are not used to estimate causal effects, we need not make onerous assumptions about the treatment assignment mechanism. However, we do need to make some assumptions to establish comparability between the observational and experimental data – assumptions that will be less likely to hold if the experiment incorporates inclusion criteria. Furthermore, our discussion will be limited to settings with binary outcomes, in which computations are tractable. We suppose the experiment has a stratified design, and seek to determine allocations of units to strata and treatment assignments.

Suppose pilot estimates of the stratum potential outcome variances are obtained from the observational study. If the outcomes are binary, we show that recent advances in sensitivity analysis from Zhao, Small, and Bhattacharya [5] can be extended to generate confidence sets for these variances, while incorporating the possibility of unmeasured confounding. Next, we pose the experimental design problem as a regret minimization problem subject to the potential outcome variances lying within their confidence sets. We use a trick from von Neumann to convert the problem into a concave maximization. The problem is not compliant with disciplined convex programming (DCP) [6], but it can be solved using projected gradient descent. This approach can yield modest efficiency gains in the experiment, especially if there is heterogeneity in treatment effects and baseline incidence rates across strata.

The remainder of the article proceeds as follows. Section 2 briefly reviews related literature, while Section 3 defines our notation, assumptions, and loss function. Section 4 gives our main results. These include the derivation of bias-aware confidence sets for the pilot variance estimates; the formulation of the design problem as a regret minimization; and the strategy to convert that problem into a computationally tractable one. We demonstrate the practical utility of our methods on data from the Women’s Health Initiative in Section 5. Section 6 discusses future work and concludes.

2 Related work

Our focus is on using observational data for experimental design, rather than for inference. We briefly review challenges in using so-called “data fusion” methods [7] that seek to merge observational and experimental data directly.

A key question is whether researchers can assume unconfoundedness – roughly, that all variables simultaneously affecting treatment probabilities and outcomes are measured – in the observational study. Under unconfoundedness, bias can be finely controlled using statistical adjustments (see e.g. ref. [1]). Hence, observational and experimental data can be merged without the risk of inducing large biases. This is the approach used in our previous work [8]; similar assumptions are made in ref. [9]. Yet unconfoundedness is a strong and fundamentally untestable assumption, and it is unrealistic to assume in many practical settings.

Some previous studies have attempted to weaken the unconfoundedness assumption, but they frequently introduce alternative assumptions in order to proceed with merged estimation. In ref. [10], the authors assume that the hidden confounding has a parametric structure, and they suggest fitting a model to correct for the hidden confounding. In ref. [11], it is assumed the bias preserves unit-level relative rank ordering (as the authors say, “bigger causal effects imply bigger bias”). The authors consider time series data with multiple observations per unit, and they argue that their assumptions are reasonable in this setting. Yet this approach does not easily extend to the case where each unit’s outcome is observed only once.

Observational studies are also frequently included in meta-analyses, which seek to synthesize evidence across multiple studies [12]. In a recent summary of methods, Mueller et al. [13] found that recommendations for the inclusion of observational studies in systematic reviews were largely unchanged from those used for experiments. They also found little consensus on best practices for combining data. Mueller and coauthors highlight a few exceptions. Thompson et al. [14] propose estimating bias reduction based on the subjective judgment of a panel of assessors, and adjusting the observational study results accordingly. Their method requires a high degree of subject matter expertise. Prevost et al. [15] suggest a hierarchical Bayes approach in which the difference between observational and experimental results is modeled explicitly. Their results are sensitive to the choice of prior.

A number of other approaches have been suggested, such as methods that make use of Bayesian networks [16] or structural causal models [17]. Broadly, this remains an area of active research, and there is no consensus best practice for merging observational and experimental causal estimates, especially when unconfoundedness is not a tenable assumption.

We instead focus on the question of experimental design, influenced by the observational data. Many recent papers have considered a closely related problem: adaptive randomization in multi-stage trials (see e.g. [18,19]). In multi-stage trials, the pilot data (or “first-stage data”) emerges not from an observational study, but instead from a randomized controlled trial (RCT). The comparative trustworthiness of these data allows for considerable flexibility in using the data to improve the design of a subsequent experiment.

In ref. [20], Tabord-Meehan considers the problem of a two-stage RCT. Unlike the setting of this study, Tabord-Meehan does not suppose that the strata are defined ahead of time. He seeks to minimize variance in estimation of the average treatment effect (ATE), rather than an L 2 loss across strata. Leveraging the reliability of the first-stage data, he proposes estimating a stratification tree using these data. Then, the choice of stratification variables, stratum delimiters for those variables, and assignment probabilities for each individual stratum in the second stage are all determined using the first-stage data. This procedure achieves a notion of asymptotic optimality among estimators utilizing stratification trees.

Bai [21] also considers randomization procedures that are informed by pilot data. He proposes a procedure in which units are first ranked according to the sum of the expectations of their treated and untreated potential outcomes (conditional on covariates), then matched into pairs with their adjacent units, with treatment randomized to exactly one member of each matched pair. Because the ranking depends on an unknown quantity, a large pilot study is required to implement this method. Bai also discusses the case in which pilot data are unavailable, in which case he proposes using the minimax framework to choose the matched-pair design that is optimal under the most adversarial data-generating process, subject to mild shape constraints on the conditional expectations of potential outcomes given covariates.

These papers share many similar goals and analytic techniques to this manuscript. Crucially, we consider the case of an L 2 loss over a fixed stratification, rather than an estimation of the ATE. Moreover, our pilot data are assumed to come from an observational study, rather than an experiment. The data are potentially informative, but significantly less reliable than a pilot RCT.

3 Problem set-up

3.1 Sources of randomness

We assume we have access to an observational study with units i in indexing set O such that O = n o . We associate with each unit i O a pair of unseen potential outcomes ( Y i ( 0 ) , Y i ( 1 ) ) ; an observed covariate vector X i , where X i R p ; and a propensity score p i ( 0 , 1 ) denoting that probability of receiving treatment. We also associate with each i a treatment indicator Z i and an observed outcome defined by Y i = Z i Y i ( 1 ) + ( 1 Z i ) Y i ( 0 ) .

There are multiple perspectives on randomness in causal inference. In the setting of ref. [22] – as in much of the early potential outcome literature – all quantities are treated as fixed except the treatment assignment Z i . More modern approaches sometimes treat the potential outcomes Y i ( 0 ) and Y i ( 1 ) and covariates X i as random variables (see e.g. ref. [23]). Similarly, some authors treat all of the data elements (including the treatment assignment Z i ) as random draws from a super-population (see e.g. ref. [1]). Per the discussion in ref. [24], these subtleties often have little effect on the choice of estimators, but they do affect the population to which results can be generalized.

In our setting, we assume that the experimental data have not yet been collected, so it does not make sense to talk about fixed potential outcomes. More naturally, we treat the potential outcomes and covariates as random for both the observational and experimental datasets. Thus, for units i O , we view ( Y i ( 0 ) , Y i ( 1 ) , X i ) as drawn from a joint distribution F O . Similarly, the experimental data will be denoted ( Y i ( 0 ) , Y i ( 1 ) , X i ) for i , sampled from a joint distribution F R . Because we are treating the potential outcomes as random variables, we can reason about their means and variances under the distribution F R .

3.2 Stratification and assumptions

We suppose we have a fixed stratification scheme based on the covariates X i . This can be derived from substantive knowledge or from applying a modern machine learning algorithm on the observational study to uncover treatment effect heterogeneity (see e.g. ref. [25,26]). The stratification is such that there are k = 1 , , K strata and each has an associated weight w 1 , , w K , where w k > 0 for all k and k w k = 1 . The w k define the relative importance of the strata and thus ordinarily reflect their prevalence in a population of interest.

Using the stratification on the observational study, we define indexing subsets O k (with cardinalities n o k ) to identify units in each stratum. For each stratum, define k as the set of covariate values defining the stratum, such that X i k i O k .

Suppose we have a budget constraint such that we can recruit only n r total units for the experiment, which we will also refer to as an RCT. One goal of our procedure is to decide the number of units n r k recruited for each stratum, subject to the constraint k n r k = n r . Once the experimental units are recruited, we will identically define indexing subsets k such that X i k i k . Within each stratum k , a second goal of our procedure will be to decide the count of units we will assign to the treatment vs. control conditions, such that the associated counts n r k t and n r k c sum to n r k . Hence, our variables of interest will be { ( n r k t , n r k c ) } 1 K .

We will make the following assumption about allocation to treatment.

Assumption 1

(Allocations to treatment) For each observational unit i O , treatment is allocated via an independent Bernoulli trial with success probability p i ( 0 , 1 ) . For the experimental units, treatment is allocated stratum-wise by drawing a simple random sample of size n r k t treated units from the n r k total units within stratum k .

Under Assumption 1, the experiment is a stratified randomized experiment [1], and the number of treated units in each stratum is fixed ahead of time.

Define E R , Var R , E O , and Var O as expectations and variances under the distributions F R and F O , respectively. We will need two further assumptions for our derivations.

Assumption 2

(Common potential outcome means) Conditional on the stratum, the potential outcome averages for the two populations are equal. In other words,

E R ( Y i ( 0 ) X i k ) = E O ( Y i ( 0 ) X i k ) and E R ( Y i ( 1 ) X i k ) = E O ( Y i ( 1 ) X i k )

for all k 1 , , K . We denote these shared quantities as μ k ( 0 ) and μ k ( 1 ) , respectively.

Assumption 3

(Common potential outcome variances) Conditional on the stratum, the potential outcome variances for the two populations are equal. In other words,

Var R ( Y i ( 0 ) X i k ) = Var O ( Y i ( 0 ) X i k ) and Var R ( Y i ( 1 ) X i k ) = Var O ( Y i ( 1 ) X i k )

for all k 1 , , K . We denote these shared quantities as σ k 2 ( 0 ) and σ k 2 ( 1 ) , respectively.

Assumptions 2 and 3 establish commonality between the observational and experimental datasets. Assumption 3 is needed explicitly to relate the optimal experimental design to quantities estimated from the observational study. These assumptions are not testable, though they need not hold exactly for the proposed methods to generate improved experimental designs. Researchers must apply subject matter knowledge to assess their approximate viability. For example, in cases in which the RCT units are sampled from the same underlying population as the observational units, these assumptions are likelier to hold. However, if the experiment incorporates onerous inclusion criteria such that the covariate distributions within stratum differ significantly between experimental and observational datasets, Assumptions 2 and 3 may be less plausible.

3.3 Loss and problem statement

Given Assumption 2, we can define a mean effect,

τ k = E R ( Y i ( 1 ) Y i ( 0 ) X i k ) = E O ( Y i ( 1 ) Y i ( 0 ) X i k ) = μ k ( 1 ) μ k ( 0 )

for each k 1 , , K . We can collect these values into a vector τ .

Denote the associated causal estimates derived from the RCT as τ ˆ r k for k = 1 , , K . We can collect these estimates into a vector τ ˆ r . We use a weighted L 2 loss when estimating the causal effects across strata,

( τ , τ ˆ r ) = k w k ( τ ˆ r k τ k ) 2 .

Our goal will be to minimize the risk, defined as an expectation of the loss over both the treatment assignments and the potential outcomes. For simplicity, we suppress the subscript and write

(1) ( τ , τ ˆ r ) = E R k w k ( τ ˆ r k τ k ) 2 = k w k σ k 2 ( 1 ) n r k t + σ k 2 ( 0 ) n r k c .

4 Converting to an optimization problem

4.1 Decision framework

Were ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) k = 1 K known exactly, it would be straightforward to compute optimal allocations in the RCT. The optimal choice from minimizing (1) is simply:

(2) n r k t = n r w k σ k ( 1 ) j w j ( σ j ( 1 ) + σ j ( 0 ) ) , n r k c = n r w k σ k ( 0 ) j w j ( σ j ( 1 ) + σ j ( 0 ) ) ,

which yields a risk of

1 n r k w k ( σ k ( 1 ) + σ k ( 0 ) ) 2 .

Note that the expressions in (2) are closely related to the well-known Neyman allocation formulas for stratified sampling [27]. In our setting, we are allowing for arbitrary stratum weights, but we are imposing a sample size constraint rather than a cost constraint, as is frequently used in the Neyman allocations. We will continue using a sample size constraint for the remainder of the article. It is straightforward to extend this work to the setting in which the treated and control arms have different costs, and the constraint is imposed in terms of cost rather than sample size. These formulas are computed explicitly in Appendix D.

Assumption 3 guarantees shared potential outcome variances across the observational and RCT datasets. One approach would be to obtain pilot estimates of σ k 2 ( 1 ) and σ k 2 ( 0 ) from the observational study and then plug them into the expressions in (2) to determine the allocation of units in the RCT. We refer to this approach as the “naïve allocation.” However, any estimate of the variances derived from the observational study should be treated with caution. Our assumptions do not preclude the possibility of unmeasured confounding, which can introduce substantial bias into the pilot estimation step. Hence, we would be better served by a framework that explicitly accounts for uncertainty in the pilot estimates.

A number of heuristic approaches are appealing. The experimenter might, for example, take a weighted average between the naïve allocation and a design that allocates units equally across strata and treatment arms. Such an approach would rely on a subjective weighting to account for the possibility of unmeasured confounding, but would be difficult to calibrate in practice. Alternatively, the experimenter might seek to develop confidence regions for the pilot estimates of σ k 2 ( 1 ) and σ k 2 ( 0 ) and solve for the best possible allocation consistent with these regions. But such an approach would be fundamentally optimistic and would ignore the possibility that σ k 2 ( 1 ) and σ k 2 ( 0 ) could take on more adversarial values.

We argue that the problem is somewhat asymmetric. Were the experimenter to ignore the observational data and use a sensible default allocation – e.g., equal allocation – they might lose some efficiency, but they would likely obtain a fairly good estimate of τ . Hence, we argue that one should incorporate the observational data somewhat cautiously and seek a strong guarantee that doing so will not make the estimate worse. Decision theory provides an attractive framework in the form of regret minimization [28,29]. In this framework, a decision-maker chooses between multiple prospects and cares about not only the received payoff but also the foregone choice. If the foregone choice would have yielded higher payoff than the chosen one, the decision-maker experiences regret [30]. Decisions are made to minimize the maximum possible regret.

In our case, the decision is on how to allocate units in our RCT. One choice is an allocation informed by the observational study. The other is a “default” allocation against which we seek to compare. Denote the default values as n ˜ r k t and n ˜ r k c , where a common choice would be equal allocation, n ˜ r k t = n ˜ r k c = n r / 2 K for all k ; or weighted allocation n ˜ r k t = n ˜ r k c = w k n r / 2 for all k . Regret is defined as the difference between the risk of our chosen allocation and the default allocation,

Regret ( { n r k t , n r k c } k = 1 K ) = k w k σ k 2 ( 1 ) 1 n r k t 1 n r k t ˜ + σ k 2 ( 0 ) 1 n r k c 1 n r k c ˜ .

Choosing this as our objective, we can now begin to formulate an optimization problem.

Suppose we can capture our uncertainty about ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) via a convex constraint, indexed by a user-defined parameter Γ ,

( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K ,

where A k ( Γ ) R 2 . We could then obtain the regret-minimizing unit allocations as the solution to

(3) min n r k t , n r k c max σ k 2 ( 1 ) , σ k 2 ( 0 ) k w k σ k 2 ( 1 ) 1 n r k t 1 n r k t ˜ + σ k 2 ( 0 ) 1 n r k c 1 n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K k n r k t + n r k c = n r .

Crucially, observe that the objective in Problem (3) can be set to zero by choosing n r k t = n ˜ r k t and n r k c = n ˜ r k c for k = 1 , , K , and this allocation must satisfy the sample size constraint by definition. Hence, the problem will only return an allocation other than the default in the case that such an allocation outperforms the default under all constraint-satisfying possible values of the variances σ k 2 ( 1 ) , σ k 2 ( 0 ) , k = 1 , , K . This captures our intuition about the asymmetry of the problem.

Defining and solving Optimization Problem (3) will be the goal of the remainder of this article.

4.2 Tractable case: binary outcomes

To construct our confidence regions A k , k = 1 , , K , we will extend recent sensitivity analysis results from Zhao et al. [5].

The authors consider the case of causal estimation via inverse probability of treatment weighting (IPW). They focus on observational studies and consider the case where unmeasured confounding is present. To quantify this confounding, they rely on the marginal sensitivity model of Tan [31]. In this model, the degree of confounding is summarized by a single researcher-chosen value, Γ 1 , which bounds the odds ratio of the treatment probability conditional on the potential outcomes and covariates and the treatment probability conditional only on covariates. The Tan model extends the widely used Rosenbaum sensitivity model [32] to the setting of IPW.

Zhao and co-authors focus on developing valid confidence intervals for the ATE even when Γ -level confounding may be present. They offer two key insights. First, they demonstrate that for any choice of Γ , one can efficiently compute upper and lower bounds on the true potential outcome means via linear fractional programming. These bounds, referred to as the “partially identified region,” quantify the possible bias in the point estimate of the ATE. Second, the authors show that the bootstrap is valid in this setting. Hence, they propose drawing repeated bootstrap replicates; computing extrema within each replicate using their linear fractional programming approach; and then taking the relevant α -level quantiles of these extrema. This procedure yields a valid α -level confidence region for the ATE.

Practically speaking, the choice of Γ is crucial in establishing the appropriate width of the confidence intervals. A common approach is to calibrate the choice of Γ against the disparities in treatment probability caused by omitting any of the observed variables [33,34]. The central logic to this approach is that unobserved covariates are unlikely to have affected the treatment probability more than any of the relevant measured covariates that are available in the dataset. A broader treatment on how to choose sensitivity parameters can be found in the study by Hsu and Small [35].

We adapt this approach to our setting in the case of binary outcomes. Note that if Y i { 0 , 1 } , then potential outcome variances can be expressed directly as a function of potential outcome means, via

σ k 2 ( 1 ) = μ k ( 1 ) ( 1 μ k ( 1 ) ) and σ k 2 ( 0 ) = μ k ( 0 ) ( 1 μ k ( 0 ) ) .

In this setting, note also that Assumption 2 implies Assumption 3.

As the work of Zhao et al. provides the necessary machinery to bound mean estimates, we can exploit this relationship between the means and variances to bound variance estimates. In particular, we can show that the bootstrap is also valid if our estimand is μ k ( e ) ( 1 μ k ( e ) ) , rather than μ k ( e ) , for e { 0 , 1 } and k = 1 , , K . Computing the extrema is also straightforward. Note that the function f ( x ) = x ( 1 x ) is monotonically increasing in x if 0 < x < 0.5 and monotonically decreasing in x if 0.5 < x < 1 . Hence, if we use the method used in ref. [5] to solve for a partially identified region for μ k ( 1 ) and μ k ( 0 ) , we can equivalently compute such intervals for σ k 2 ( 1 ) and σ k 2 ( 0 ) .

Denote as μ ˆ k U ( e ) the upper bound and μ ˆ k L ( e ) the lower bound computed for a mean for e { 0 , 1 } . Denote ( σ ˆ k 2 ( e ) ) U and ( σ ˆ k 2 ( e ) ) L as the analogous quantities for variance. We apply the following logic:

  • If μ ˆ k U ( e ) 0.5 , set

    (4) ( σ ˆ k 2 ( e ) ) L = μ ˆ k L ( e ) ( 1 μ ˆ k L ( e ) ) and ( σ ˆ k 2 ( e ) ) U = μ ˆ k U ( e ) ( 1 μ ˆ k U ( e ) ) .

  • If μ ˆ k L ( e ) 0.5 , set

    (5) ( σ ˆ k 2 ( e ) ) L = μ ˆ k U ( e ) ( 1 μ ˆ k U ( e ) ) and ( σ ˆ k 2 ( e ) ) U = μ ˆ k L ( e ) ( 1 μ ˆ k L ( e ) ) .

  • If μ ˆ k L ( e ) < 0.5 < μ ˆ k U ( e ) , set

    (6) ( σ ˆ k 2 ( e ) ) L = min ( μ ˆ k L ( e ) ( 1 μ ˆ k L ( e ) ) , μ ˆ k U ( e ) ( 1 μ ˆ k U ( e ) ) ) and ( σ ˆ k 2 ( e ) ) U = 0.25 .

Hence, we propose the following procedure for deriving valid confidence regions for ( σ k 2 ( 0 ) , σ k 2 ( 1 ) ) for each choice of k :

  1. Draw B bootstrap replicates from the units i O k .

  2. For each replicate:

    1. Compute μ ˆ k U ( e ) , μ ˆ k L ( e ) for e { 0 , 1 } using Zhao and co-authors’ linear fractional programming approach.

    2. Determine ( σ ˆ k 2 ( e ) ) U and ( σ ˆ k 2 ( e ) ) L for e { 0 , 1 } using the approach described in (4), (5), and (6).

  3. Each replicate can now be represented as a rectangle in [ 0 , 1 ] × [ 0 , 1 ] , where one axis represents the value of ( σ ˆ k 2 ( 1 ) ) , and the other the value of ( σ ˆ k 2 ( 0 ) ) , and the vertices correspond to the extrema. Any set such that a 1 α proportion of the rectangles have all four corners included in the set will asymptotically form a valid α -level confidence interval.

A full proof of the validity of this method can be found in Appendix B.

Note that the final step does not specify the shape of the confidence set (it need not even be convex). For simplicity, we compute the minimum volume ellipsoid containing all vertices, then shrink the ellipsoid toward its center until only B ( 1 α ) of the rectangles have all four of their vertices included. For details on constructing the ellipsoids (sometimes known as Löwner–John ellipsoids), see ref. [36]. Observe that this is by no means the smallest valid confidence set, but it is convex and easy to work with numerically. In Appendix E, we briefly discuss the use of rectangular confidence regions, finding that results are substantively similar.

In Figure 1, we demonstrate this procedure on simulated data using Γ = 1.2 . We suppose there are four strata, each containing 1,000 observational units. The strata differ in their treatment probabilities with 263, 421, 564, and 739 treated units in each stratum, respectively. The large black dot at the center of each cluster represents the point estimate ( σ ˆ k 2 ( 0 ) , σ ˆ k 2 ( 1 ) ) . In purple, we plot the rectangles corresponding to the extrema computed in each of 200 bootstrap replicates drawn from the data. The dashed ellipsoids represent 90% confidence sets. In the cases of strata 2 and 3, the ellipsoids extend beyond the upper bound of 0.25 in at least one direction, so we intersect the ellipsoids with the hard boundary at 0.25. The resulting final confidence sets, A 1 , A 2 , A 3 , and A 4 , are all convex.

Figure 1 
                  Simulated example of confidence regions in four strata under 
                        
                           
                           
                              Γ
                              =
                              1.2
                           
                           \Gamma =1.2
                        
                     .

Figure 1

Simulated example of confidence regions in four strata under Γ = 1.2 .

The objective is convex in n r k t , n r k c and affine (and thus concave) in σ k 2 ( 1 ) , σ k 2 ( 0 ) . Now, having obtained convex constraints, we can invoke von Neumann’s minimax theorem [37] to switch the order of the minimization and maximization. Hence, the solution to Problem (3) is equivalent to the solution of

max σ k 2 ( 1 ) , σ k 2 ( 0 ) min n r k t , n r k c k w k σ k 2 ( 1 ) 1 n r k t 1 n r k t ˜ + σ k 2 ( 0 ) 1 n r k c 1 n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K k n r k t + n r k c = n r .

But the inner problem has an explicit solution, given by the expressions in (2). Plugging in these expressions, we arrive at the simplified problem:

(7) max σ k 2 ( 1 ) , σ k 2 ( 0 ) 1 n r k w k ( σ k ( 1 ) + σ k ( 0 ) ) 2 k w k σ k 2 ( 1 ) n r k t ˜ + σ k 2 ( 0 ) n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K .

Problem (7) is concave. See Appendix C for a detailed proof. The solution is non-trivial, owing to the fact that the problem is not DCP-compliant. Nonetheless, a simple projected gradient descent algorithm is guaranteed to converge under very mild conditions given the curvature [38]. Similarly, under mild conditions, the convergence rate can be shown to be linear (see e.g. ref. [39]), meaning that distance to the optimum declines at a rate of O ( 1 / m ) , where m is the number of steps taken by the algorithm. Hence, we can efficiently solve this problem.

5 Application to the data from the Women’s Health Initiative

5.1 Setting

To evaluate our methods in practice, we make use of data from the Women’s Health Initiative (WHI), a 1991 study of the effects of hormone therapy on postmenopausal women. The study included both an RCT and an observational study. A total of 16,608 women were included in the trial, with half randomly selected to take 625 mg of estrogen and 2.5 mg of progestin, and the remainder receiving a placebo. A corresponding 53,054 women in the observational component of the WHI were deemed clinically comparable to women in the trial. About a third of these women were using estrogen plus progestin, while the remaining women in the observational study were not using hormone therapy [40].

We investigate the effect of the treatment on incidence of coronary heart disease. We split the data into two non-overlapping subsets, which we term the “gold” and “silver” datasets. We estimate the probability of treatment for observational units via fitted propensity scores. The data split is the same as the one used in ref. [8]. Details on the construction of these data elements can be found in Section A.2, while further details about the WHI can be found in Section A.1.

To choose our subgroups for stratification, we utilize the clinical expertise of researchers in the study’s writing group. The trial protocol highlights age as an important subgroup variable to consider [41], while subsequent work considered a patient’s history of cardiovascular disease [42]. To evaluate the impact of a clinically irrelevant variable, we also consider langley scatter, a measure of solar irradiance at each woman’s enrollment center, which is not plausibly related to baseline incidence or treatment effect. Langley scatter exhibits no association with the outcome in the observational control population: a Pearson’s Chi-squared test yields a p-value of 0.89. The analogous tests for age and history of cardiovascular disease have p-values below 1 0 5 .

The age variable has three levels, corresponding to whether a woman was in her 50s, 60s, or 70s. The cardiovascular disease history variable is binary. The langley scatter variable has five levels, corresponding to strata between 300 and 500 langleys of irradiance. We provide brief summaries of these variables in Tables A.7, A.8, and A.9 in Section A.3.

The RCT gold dataset is used to estimate “gold standard” stratum causal effects. We suppose that the observational study is being used to assist the design of an experiment of size n r = 1,000 units. In all cases, the default allocation is an equal allocation across strata and treatment statuses.

We face the additional challenge of choosing the appropriate value of Γ . The WHI provides a very rich set of covariates, and our propensity model incorporates more than 50 variables spanning the demographic and clinical domains (see details in Section A.2). Hence, we will run our algorithm at values of Γ = 1.0 (reflecting no residual confounding) as well as Γ = 1.1 , 1.5 , and 2.0 (reflecting a modest amount).

5.2 Detailed example: Γ = 1.5 , fine stratification

We show one example in detail, in which we choose Γ = 1.5 and stratify on all three subgroup variables: age, history of cardiovascular disease, and langley scatter. The cross-product of these variables yields 30 strata, which we suppose are weighted equally. We number these groups from 1 through 30.

In the top panel of Figure 2, we show a naïve RCT allocation based purely on the pilot estimates of the stratum potential outcome variances from the observational study. In the bottom panel, we show the regret minimizing allocations. Visually, it is clear that we have heavily shrunk the allocations toward an equally allocated RCT, but there remain some strata where we recommend over- or under-sampling. Note, too, that the shrinkage is not purely reflective of the magnitude of the pilot estimate, since the number of observational units from each stratum and treatment status also influences the width of our confidence region for each of the pilot estimates.

Figure 2 
                  Allocation of units to strata under naïve allocation and regret-minimizing allocation.

Figure 2

Allocation of units to strata under naïve allocation and regret-minimizing allocation.

To investigate the utility of our regret-minimizing allocations, we sample pseudo-experiments of 1,000 units from the RCT silver dataset 1,000 times with replacement. We do so under three designs: equal allocation by strata; naïve allocation based on the pilot estimates; and the regret-minimizing allocations under Γ = 1.5. We compute the average L 2 loss when compared against the gold standard estimates derived from the RCT gold dataset. Results are shown in Figure 3. Our method yields a modest reduction in average loss (3.6%) relative to the naïve design. It also outperforms the equal design, though by a slimmer margin (1.6%).

Figure 3 
                  Average loss over 1,000 resamples of 1,000-unit experiments under equal allocation, naïve allocation, and regret-minimizing allocation designs.

Figure 3

Average loss over 1,000 resamples of 1,000-unit experiments under equal allocation, naïve allocation, and regret-minimizing allocation designs.

5.3 Performance over multiple conditions

We now simulate with all possible combinations of the stratification variables. For each choice of a stratification, we sample 1,000 units from the RCT silver dataset with replacement, under equal allocation, naïve allocation, and regret-minimizing allocation with Γ = 1.0 , 1.1 , 1.5 , and 2.0. We then compute the L 2 loss versus the “gold standard” estimates derived from the RCT gold dataset.

In Table 1, we summarize the loss of the regret-minimizing allocations relative to equal allocation. We see immediately that the entries are all non-positive. This makes some intuitive sense: the objective in Problem (3) can always be set to 0 by choosing n r k t = n ˜ r k t and n r k c = n ˜ r k c for all k ; hence, the algorithm is designed to guarantee that we cannot do worse than allocating equally. By the same token, many of the gains we see are modest, owing to the conservatism of the regret-minimizing approach. Notably, we seem to achieve the greatest gains when we are stratifying only on clinically relevant variables and using a relatively low value of Γ . We achieve a 5–6% risk reduction at low values of Γ in the fourth row of the table, in which we stratify on the clinically relevant age and cardiovascular disease variables. On the other hand, the algorithm quickly defaults to recommending equal allocation when variables are not clinically relevant. In the third row, in which we stratify only on the irrelevant langley scatter variable, the starred entries correspond to cases in which the regret-minimizing allocation is equal allocation.

Table 1

L 2 loss comparisons for regret-minimizing allocations relative to equal allocation

Subgroup Var(s) Equal alloc loss Loss relative to equal allocation
Γ = 1 (%) Γ = 1.1 (%) Γ = 1.5 (%) Γ = 2 (%)
Age 0.000517 2.0 1.9 2.0 0.0
CVD 0.000498 2.3 2.0 1.5 0.0
Langley 0.000841 0.0 0 . 0 0 . 0 0 . 0
Age, CVD 0.001541 5.5 5.6 3.8 2.3
Age, langley 0.003417 1.6 1.6 0.7 0.1
CVD, langley 0.002495 1.7 1.2 0.8 0.2
Age, CVD, langley 0.008395 1.9 2.1 1.6 0.7

For starred entries, the regret-minimizing allocation defaults to equal allocation.

In Table 2, we summarize the loss relative to naïve allocation. our method can underperform a naïve allocation derived from the observational study pilot variance estimates. This can be seen most clearly in the first row of the table, in which we stratify only on the age variable. However, there are two clear trends in the results. First, when we stratify on a variable that turns out not to be clinically relevant, like langley scatter, the naïve allocation is essentially recommending an allocation based on noise from the data; as a result, our regret-minimizing allocations uniformly outperform naïve allocations. Second, the regret-minimizing allocations tend to outperform the naïve allocations as the number of strata grows. We significantly outperform naïve allocation in the final row, which corresponds to stratification on all three variables and a total of 30 strata.

Table 2

L 2 loss comparisons for regret-minimizing allocations relative to naïve allocation

Subgroup Var(s) Naïve alloc loss Loss relative to naïve allocation
Γ = 1 (%) Γ = 1.1 (%) Γ = 1.5 (%) Γ = 2 (%)
Age 0.000501 1.2 1.2 1.1 3.2
CVD 0.000488 0.3 0.0 0.6 2.1
Langley 0.000852 1.1 1.3 1.3 1.3
Age, CVD 0.001484 1.8 1.9 0.1 1.5
Age, langley 0.003393 0.9 0.9 0.0 0.6
CVD, langley 0.002481 1.1 0.7 0.3 0.3
Age, CVD, langley 0.008574 3.9 4.1 3.6 2.8

Recall that as Γ rises, the feasible set of Optimization Problem (3) grows larger. Hence, we expect the allocation to be closer to the naïve allocation for smaller values of Γ , but to be regularized more toward the default allocation for larger values of Γ . For large Γ , we would thus expect the loss to converge to the equal allocation loss. This is precisely what we see in Table 1: for each possible stratification, the performance is closest to that of the default allocation when Γ = 2 . However, in Table 2, we do not see the inverse pattern – that is, performance is not uniformly closest to that of the naïve allocation when Γ = 1 . This is because the confidence set does not collapse to a single point at Γ = 1 ; rather, it incorporates the possibility of variance but not bias in the pilot estimation. More broadly, we do not expect a monotone relationship between Γ and the average loss. In many cases, the pilot estimates will be somewhat informative, but incorporate some bias. Hence, we may see the lowest average loss at intermediate values of Γ , which encourage the algorithm to extract some relevant information from the pilot data without relying too heavily on these estimates.

While these simulation results show modest performance gains, they are encouraging. A wise analyst would be cautious about designing an RCT exclusively using observational study pilot estimates of potential outcome variances. Because such pilot estimates can have both bias and variance, relying too heavily upon them might waste resources. Our framework allows data from the observational study to be incorporated into the RCT design while guarding against the possibility of underperforming a default allocation.

6 Extensions

We briefly discuss potential extensions of this work.

One natural consideration is the case of multiple treatment levels, rather than the binary setting of treatment versus control. The machinery discussed in this manuscript naturally extends to the multilevel case. If we suppose there are L treatment levels, then we instead optimize over sample sizes n r k and stratum potential outcome variances σ k 2 ( ) for { 1 , , L } . The optimization problem becomes:

(8) min { n r k } , k max { σ k 2 ( ) } , k k =1 K w k =1 L σ k 2 ( ) 1 n r k 1 n ˜ r k subject to ( σ k 2 (1), , σ k 2 ( L ) ) A k ( Γ ) , k = 1 , , K k =1 K =1 L n r k = n r .

The curvature of Problem (8) is unchanged from that of Problem (3), so we can use the same von Neumann trick to obtain a readily solvable concave maximization problem. The only remaining complexity is the construction of the confidence sets A k ( Γ ) . The procedure described in Section 4.2 can be easily generalized to the multilevel case, with the bounds derived from each bootstrap replicate now represented as an L -dimensional box rather than a rectangle. The proof in Appendix B does not depend on the problem’s dimensionality, so we can again obtain asymptotic α -level validity for any confidence set drawn to include a 1 α proportion of the boxes. The method of drawing Löwner–John ellipsoids also generalizes to dimensions greater than two, so we can use this exact procedure to obtain our confidence sets.

Another obvious extension is to the more general case of Y i R . In keeping with the theme of IPW estimation, we consider estimators of the form

(9) σ ˆ k 2 ( 1 ) = i O k Y i 2 Z i p i i O k Z i p i i O k Y i Z i p i i O k Z i p i 2 σ ˆ k 2 ( 0 ) = i O k Y i 2 1 Z i 1 p i i O k 1 Z i 1 p i i O k Y i 1 Z i 1 p i i O k 1 Z i 1 p i 2 ,

where p i are the true treatment probabilities. Such estimators are asymptotically unbiased.

We suppose we estimate p i with fitted propensity scores, π ˆ i , defined as

π ˆ i = 1 1 + e g ˆ ( X i ) .

We typically use logistic regression to estimate the propensity scores, such that g ˆ ( X i ) = β ˆ T X i .

We account for the possibility of Γ -level unmeasured confounding by allowing the true probability p i to satisfy

p i 1 1 + z i e g ˆ ( X i ) 1 Γ z i Γ .

We redefine the problem in terms of the v i = p i 1 , an affine function of the z i . We define two vectors v t = ( v i ) i : Z i = 1 and v c = ( v i ) i : Z i = 0 , and analogously define vectors Y t = ( Y i ) Z i = 1 and Y c = ( Y i ) Z i = 0 . Now, we can express the equations in (9) as quadratic fractional programs, e.g.,

σ ˆ k 2 ( 1 ) = v t T Θ t v t v t T 1 1 T v t , σ ˆ k 2 ( 0 ) = v c T Θ c v c v c T 1 1 T v c ,

where

Θ t = Y t 2 1 T Y t Y t T and Y c 2 1 T Y c Y c T .

We have few guarantees on the curvature of the problem: the numerators will be neither convex nor concave in the v e terms, e { 0 , 1 } , as long as the vectors 1 , Y t , and Y t 2 are linearly independent. The denominators will be convex in the v e terms. This poses a major challenge. Quadratic fractional programming problems can be solved efficiently in some special cases, but are, in general, NP-hard [43].

One avenue is to apply Dinkelbach’s method to transform the quadratic fractional problem to a series of quadratic programming problems [44]. This would not immediately yield a solution because of the indefinite numerator, but it would potentially allow one to make use of considerable recent work on solution methods in quadratic programming (see e.g. ref. [45]). This path represents a possible future extension of this work.

Acknowledgments

We thank Mike Baiocchi, Guillaume Basse, and Luke Miratrix for their useful comments and discussion.

  1. Funding information: Evan Rosenman was supported by Google, and by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program. This work was also supported by the NSF under grants DMS-1521145, DMS-1407397, and IIS-1837931.

  2. Conflict of interest: The authors have no conflicts of interest to declare.

Appendix

A.1 Further details about the Women’s Health Initiative

We evaluate our estimators on data from the Women’s Health Initiative to estimate the effect of hormone therapy on coronary heart disease (CHD). The Women’s Health Initiative is a study of postmenopausal women in the United States, consisting of RCT and observational study components with 161,808 total women enrolled [40]. Eligibility and recruitment data for the WHI can be found in the results of previous studies [3,46]. Participants were women between 50 and 79 years old at baseline, who had a predicted survival of at least 3 years and were unlikely to leave their current geographic area for 3 years.

Women with a uterus who met various safety, adherence, and retention criteria were eligible for a combined hormone therapy trial. A total of 16,608 women were included in the trial, with 8,506 women randomized to take 625 mg of estrogen and 2.5 mg of progestin, and the remainder receiving a placebo. A corresponding 53,054 women in the observational component of the Women’s Health Initiative had an intact uterus and were not using unopposed estrogen at baseline, thus rendering them clinically comparable [40]. About a third of these women were using estrogen plus progestin, while the remaining women in the observational study were not using hormone therapy [40].

Participants received semiannual contacts and annual in-clinic visits for the collection of information about outcomes. Disease events, including CHD, were first self-reported and later adjudicated by physicians. We focus on outcomes during the initial phase of the study, which extended for an average of 8.16 years of follow-up in the RCT and 7.96 years in the observational study.

The overall rate of coronary heart disease in the trial was 3.7% in the treated group (314 cases among 8,472 women reporting) versus 3.3% (269 cases among 8,065 women reporting) for women not randomized to estrogen and progestin. In the observational study, the corresponding rates were 1.6% among treated women (706 out of 17,457 women reporting) and 3.1% among control women (1,108 out of 35,408 women reporting). Our methodology compares means and not survival curves. In the initial follow-up period, death rates were relatively low in both the observational study (6.4%) and the randomized trial (5.7%). Hence, we do not correct for the possibility of these deaths censoring coronary heart disease events.

A.2 Propensity score construction, covariate balance, and gold standard effects

The Women’s Health Initiative researchers collected a rich set of covariates about the participants in the study. For the purposes of computational speed, we narrow to a set of 684 variables, spanning demographics, medical history, diet, physical measurements, and psychosocial data collected at baseline.

The most meaningful measure of covariate imbalance can be found by looking at clinically relevant factors. Prentice et al. [40] found that hormone therapy users in the observational study were more likely to be Caucasian or Asian/Pacific Islander, less likely to be overweight, and more likely to have a college degree. These imbalances strongly suggest that applying a naïve differencing estimate to the observational data will yield an unfairly rosy view of the effect of hormone therapy on CHD.

To generate our estimators for this dataset, we need a propensity model e ( x ) to map the observed covariates to an estimated probability of receiving the treatment in the observational study. We used a logistic regression to generate an expressive model while limiting overfit. A heuristic procedure was used for careful construction of the propensity scores. Full details can be found in ref. [8].

Matching on the propensity score should reduce imbalances on clinically relevant covariates. We provide a summary of these diagnostics, with additional details again available in ref. [8]. We use standardized differences to measure imbalance, as advocated by Rosenbaum [47].

We compute the standardized differences between treated and control on risk factors listed in ref. [40], before and after adjusting for the propensity score by grouping the units into ten equal-width propensity score strata. With the exception of the physical functioning score, all evaluated covariates were included in the propensity model. Imbalance measures for the continuous covariates can be found in Table A1. The stratification procedure reduces all standardized differences to less than 0.05 in absolute value, representing very good matches between the populations.

Table A1

Standardized differences (SD) between treated and control populations in the observational dataset, before and after stratification on the propensity score, for clinical risk factors for coronary heart disease

Before stratifying After stratifying
Variable Test Ctrl SD Test Ctrl SD
Age 60.78 64.72 0.56 63.06 63.33 0.04
BMI 25.55 27.11 0.25 26.71 26.62 0.00
Physical functioning 85.23 79.58 0.26 81.15 81.23 0.03
Age at menopause 50.49 50.19 0.06 50.35 50.33 0.02

For categorical variables, the stratification procedure reweights individual women, such that the effective proportion of women in each category changes after stratifying on the propensity score. Standardized differences can also be computed for categorical variables, using the procedure described in Graziano et al. [48]. We achieve similar balance on two significant categorical variables – ethnicity and smoking status – in Tables A2 and A3.

Table A2

Standardized differences (SD) between treated and control populations in the observational database, before and after stratification on the propensity score, for ethnicity category

White (%) Black (%) Latino (%) AAPI (%) Native American (%) Missing/others (%) SD
Before stratifying Treated 89.0 2.7 2.9 4.0 0.2 1.1 0.26
Control 83.1 8.1 3.9 2.8 0.4 1.5
After stratifying Treated 83.4 6.9 4.3 3.6 0.5 1.4 0.05
Control 84.8 6.4 3.6 3.4 0.4 1.4

Table A3

Standardized differences (SD) between treated and control populations in the observational database, before and after stratification on the propensity score, for smoking category

Never smoked (%) Past smoker (%) Current smoker (%) SD
Before stratifying Treated 48.7 46.2 5.1 0.11
Control 52.3 41.1 6.6
After stratifying Treated 50.9 42.5 6.6 0.01
Control 51.0 42.7 6.3

Finally, we consider estimation of the “gold standard” causal effect. We randomly partition the randomized trial data into two subsets of equal size, such that each contains the same number of treated and control women. We select one of these subsets and refer to it as our “gold” dataset, to be used for estimating the true causal effect. The remaining subset is referred to as the “silver” dataset and is used for evaluating our estimators.

Because of the randomization, we find that treated and control are already well balanced on the coronary heart disease risk factors in the gold dataset, as summarized in Tables A4, A5, and A6.

Table A4

Standardized differences (SD) between treated and control populations in RCT gold dataset, for clinical risk factors for coronary heart disease

Variable Treated Control SD
Age 63.24 63.41 0.02
BMI 28.33 28.38 0.01
Physical functioning 80.97 81.11 0.01
Age at menopause 44.97 46.33 0.09

Table A5

Standardized differences (SD) between treated and control populations in RCT gold dataset, for the measured ethnicity variable

White (%) Black (%) Latino (%) AAPI (%) Native American (%) Missing/others (%) SD
Treated 84.1 6.5 5.5 2.1 0.26 1.6 0.05
Control 84.6 6.8 5.1 1.9 0.40 1.2

Table A6

Standardized differences (SD) between treated and control populations in the RCT gold dataset, for the measured smoking variable

Never smoked (%) Past smoker (%) Current smoker (%) SD
Treated 50.1 38.7 11.2 0.03
Control 50.6 39.1 10.2

A.3 Stratification variable distributions

In Tables A7, A8, and A9, we provide the distributions for variables with which we stratify in the main text.

Table A7

Distribution of age variable values in the observational study, RCT, and RCT silver datasets

Age Observational study RCT RCT silver dataset
50–59 17,447 (33.0%) 5,491 (33.2%) 2,806 (33.9%)
60–69 23,030 (43.6%) 7,473 (45.2%) 3,689 (44.6%)
70–79 12,388 (23.4%) 3,573 (21.2%) 1,774 (21.5%)

Table A8

Distribution of history of cardiovascular disease in the observational study, RCT, and RCT silver datasets

History of CVD Observational study RCT RCT silver dataset
Yes 8,709 (16.5%) 1,828 (11.1%) 900 (10.9%)
No 44,156 (83.5%) 14,709 (88.9%) 7,369 (89.1%)

Table A9

Distribution of langley scatter (g-cal/ cm 2 ) variable across categories in the observational study, RCT, and RCT silver datasets

Langley scatter Observational study RCT RCT silver dataset
300–325 15,599 (29.5%) 4,854 (29.4%) 2,411 (29.2%)
350 12,521 (23.7%) 3,917 (23.7%) 1,935 (23.4%)
375–380 5,841 (11.0%) 1,858 (11.2%) 934 (11.3%)
400–430 8,216 (15.5%) 2,585 (15.6%) 1,310 (15.8%)
475–500 10,688 (20.2%) 3,323 (20.1%) 1,679 (20.3%)

B Proof of validity of confidence regions

We hew closely to the proofs provided in ref. [5]. The primary proofs consider the missing data problem, which is equivalent to estimating the mean of either of the potential outcomes in our setting. We begin by providing a summary of key proofs. All references to Remarks, Lemmas, etc. in Section B.1 refer to the text of ref. [5].

B.1 Review of proofs in Zhao et al. [5]

The authors define e ( x , y ) as the probability of treatment given covariates X = x X and outcome Y = y R and compare it against the marginal treatment probability e ( x ) . They use A rather than W to denote a treatment indicator, so in keeping with their notation:

e ( x , y ) = P 0 ( A = 1 X = x , Y = y ) and e ( x ) = P 0 ( A = 1 X = x ) .

Then, for any choice of Λ > 1 , they define a collection of sensitivity models

( Λ ) = 0 e ( x , y ) 1 1 Λ OR ( e ( x , y ) , e ( x ) ) Λ , for all x , y ,

where OR ( p 1 , p 2 ) = [ p 1 / ( 1 p 1 ) ] / [ p 2 / ( 1 p 2 ) ] is the odds ratio. This model was originally introduced by Tan [31]. Per Proposition 7.1, it is related to the widely used Rosenbaum sensitivity model. In keeping with that model, we use Γ rather than Λ to denote our sensitivity parameter in the text, but retain the notation Λ throughout this proof.

Via Remark 3.2, Zhao and co-authors reparameterize the problem such that each model in ( Λ ) corresponds to a choice of h ( x , y ) , the logit-scale difference of the observed probability e ( x ) , and the complete data selection probability e ( x , y ) . So, we can alternatively write:

( Λ ) = { e ( h ) ( x , y ) h ( λ ) } ,

where λ = log ( Λ ) and ( λ ) = { h : X × R h λ } . In words: every choice of h ( λ ) defines, at each possible value of X and Y , a discrepancy between e ( x ) and e ( x , y ) . The choice of ( λ ) bounds the maximum of those discrepancies. So, as Λ grows, we are allowing for greater and greater discrepancies in these probabilities.

For each choice of h , they define a “shifted estimand,”

μ ( h ) = E A e ( h ) ( x , y ) 1 E A Y e ( h ) ( x , y ) ,

where A is the treatment indicator and the expectation is over the joint distribution of X , Y , A . The corresponding “shifted estimator” is given by

μ ˆ ( h ) = 1 n i = 1 n A i e ( h ) ˆ ( X i , Y i ) 1 1 n i = 1 n A i Y i e ( h ) ˆ ( X i , Y i ) .

The sum is over a sample of points ( X i , Y i , A i ) drawn i.i.d. from their joint distribution. The quantity in the denominators, e ˆ ( h ) ( X i , Y i ) , is obtained by estimating P ( A = 1 X = x ) and then shifting the estimate by h ( x i , y i ) for all units i such that X i = x i and Y i = y i .

Now, the proof of the validity of their approach proceeds in several stages.

  1. First, they consider the case where data-dependent intervals [ L ( h ) , U ( h ) ] are asymptotically guaranteed to contain μ ( h ) with 1 α probability. They argue that taking L = inf h ( λ ) L ( h ) and U = sup h ℋ(λ) U ( h ) yields an interval [ L , H ] with asymptotic 1 α coverage for every value of μ ( h ) for which h ( λ ) (Proposition 4.1).

  2. For each choice of h ( λ ) , they establish that the bootstrap is valid (Theorem 4.2).

    1. First, they use the general theory of Z -estimators to show that μ ˆ ( h ) and its bootstrap analogue, μ ˆ ˆ ( h ) , are asymptotically normal with the same mean and variance (Theorem C.1 and Corollary C.2).

    2. Then, they conclude that by defining L B ( h ) as the α / 2 bootstrap quantile, they have

      P ( μ ( h ) < L B ( h ) ) α 2 ,

      where the expectation is taken under the joint distribution of X , Y , and A . Analogous results holds for U B ( h ) , the 1 α / 2 bootstrap quantile (Section C.3).

  3. They argue that the quantile and infimum/supremum functions can be interchanged, such that

    Q α / 2 ( inf h ( λ ) μ ˆ ˆ ( h ) ) inf h ( λ ) L ( h )

    and

    Q 1 α / 2 ( sup h ( λ ) μ ˆ ˆ ( h ) ) sup h ( λ ) U ( h )

    via Lemma 4.3.

B.2 Extension to design case

Our challenge is to extend this argument to the case where our estimand of interest is not a single μ but rather the pair ( σ k 2 ( 0 ) , σ k 2 ( 1 ) ) = ( μ k ( 0 ) ( 1 μ k ( 0 ) ) , μ k ( 1 ) ( 1 μ k ( 1 ) ) . Crucially, we will now have two h functions h 0 and h 1 , corresponding to each of the potential outcomes, but they both lie within ( λ ) . The definition of the shifted estimand under h given above generalizes to the case of two shifted estimands in a straightforward way. We extend Proposition 1 from ref. [5] in the following argument.

Proposition 1

Suppose there exists a data-dependent region β k ( h 0 , h 1 ) R 2 such that

lim inf n P ( ( σ k ( h 0 ) ( 0 ) 2 , σ k ( h 1 ) ( 1 ) 2 ) β k ( h 0 , h 1 ) ) 1 α

holds for every ( h 0 , h 1 ) ( λ ) × ( λ ) , where σ k ( h 0 ) ( e ) 2 = μ k ( h 1 ) ( e ) ( 1 μ k ( h ) ( e ) ) for e { 0 , 1 } , and n is the sample size. Under these conditions, the set

β k = h 0 , h 1 ( λ ) β k ( h 0 , h 1 )

is an asymptotic confidence set of ( σ k 2 ( 0 ) , σ k 2 ( 1 ) ) with at least 1 α coverage if h 0 , h 1 ( λ ) .

Proof

This follows from the fact that, by assumption, the true data-generating distribution satisfies h 0 , h 1 ( λ ) .□

Next, we must show that the bootstrap is valid in our setting. We adopt the same model and regularity conditions of Theorem 4.2 in ref. [5]. In their proof of Corollary 5.1, the authors show that the pairs ( μ ˆ k ( h 0 ) ( 0 ) , μ ˆ k ( h 1 ) ( 1 ) ) and ( μ ˆ ˆ k ( h 0 ) ( 0 ) , μ ˆ ˆ k ( h 1 ) ( 1 ) ) are both jointly asymptotically normal, with the same limiting distribution. We define the function

f ( x , y ) = ( x ( 1 x ) , y ( 1 y ) ) .

We can see that applying f ( ) to the tuple of potential outcome means will yield the potential outcome variances, and the same logic holds for applying f ( ) to any estimator of the potential outcome means. Moreover, because f ( ) is continuously differentiable, we can use the Delta method to observe immediately that ( σ ˆ k ( h 0 ) ( 0 ) 2 , σ ˆ k ( h 1 ) ( 1 ) 2 ) and ( σ ˆ ˆ k ( h 0 ) ( 0 ) 2 , σ ˆ ˆ k ( h 1 ) ( 1 ) 2 ) have the same asymptotic distribution, and thus the bootstrap is valid [49].

Finally, we generalize Lemma 4.3 to our setting. For each possible bootstrap replicate b { 1 , , N } where N = n n , define the quartet of points

R ˆ ˆ k , b = ( inf h 0 ( λ ) σ ˆ ˆ k , b ( h 0 ) (0) 2 , inf h 1 ( λ ) σ ˆ ˆ k , b ( h 1 ) (1) 2 ) , ( inf h 0 ( λ ) σ ˆ ˆ k , b ( h 0 ) (0) 2 , sup h 1 ( λ ) σ ˆ ˆ k , b ( h 1 ) (1) 2 ) , ( sup h 0 ( λ ) σ ˆ ˆ k , b ( h 0 ) (0) 2 , inf h 1 ( λ ) σ ˆ ˆ k , b ( h 1 ) (1) 2 ) , ( sup h 0 ( λ ) σ ˆ ˆ k , b ( h 0 ) (0) 2 , sup h 1 ( λ ) σ ˆ ˆ k , b ( h 1 ) (1) 2 ) .

In words, R ˆ ˆ k , b contains the vertices of a rectangle in R 2 which defines the extrema of the potential outcome variances consistent with h 0 , h 1 ( λ ) .

Denote as Conv ( ) the standard convex hull operator. Define a related operator,

Conv ( S , ) = Conv ( b S b ) ,

which takes in a set S of cardinality N S as well as a set { 1 , , N S } . The function returns the convex hull of the points contained in the entries in S indexed by .

We choose a set α { 1 , 2 , , N = n n } such that α = ( 1 α ) N , and we define the set

A k = Conv ( { R ˆ ˆ k , b } , α ) .

Lemma 1

The set A k is an asymptotically valid confidence set.

Proof

For 1 b N , where N = n n is the total number of possible bootstrap samples, we have that for every h 0 , h 1 ( λ ) ,

( σ ˆ ˆ k , b ( h 0 ) ( 0 ) 2 , σ ˆ ˆ k , b ( h 1 ) ( 1 ) 2 ) Conv ( R ˆ ˆ k , b ) , for all 1 b N ,

Since this holds entrywise, it follows that any set containing a fixed proportion of the sets on the RHS must contain at least that proportion of points on the LHS, and hence

Conv ( { ( σ ˆ ˆ k , b ( h 0 ) ( 0 ) 2 , σ ˆ ˆ k , b ( h 1 ) ( 1 ) 2 ) } , α ) Conv ( { Conv ( R ˆ ˆ k , b ) } , α ) .

Since this holds for every h 0 , h 1 ( λ ) , we can take the union on the LHS to observe

h 0 , h 1 ( λ ) Conv ( { ( σ ˆ ˆ k , b ( h 0 ) ( 0 ) 2 , σ ˆ ˆ k , b ( h 1 ) ( 1 ) 2 ) } , α ) Conv ( { Conv ( R ˆ ˆ k , b ) } , α ) .

Observe that the RHS is simply A k , since any convex set containing the vertices of a rectangle will contain the convex hull of those vertices as well.

On the LHS, we can make use of our bootstrap validity result to observe

liminf n P ( ( σ k ( h 0 ) ( 0 ) 2 , σ k ( h 1 ) ( 1 ) 2 ) Conv ( { ( σ ˆ ˆ k , b ( h 0 ) ( 0 ) 2 , σ ˆ ˆ k , b ( h 1 ) ( 1 ) 2 ) } , α ) ) 1 α .

It follows from Proposition 1 that the LHS is a valid 1 α level confidence region. Hence, the right-hand side must be as well.

To conclude, we observe that our ellipsoid method must necessarily comprise a superset of a convex hull for some choice of α . Hence, our method will indeed generate valid confidence regions for the potential outcome variances.□

C Proof of concavity of minimax problem

We begin with the unweighted case and demonstrate concavity by direct computation of the Hessian. Define

f ( { σ k 2 (0), σ k 2 (1)} k ) = 1 n r k σ k (1) + σ k (0) 2 k σ k 2 (1) n ˜ r k t + σ k 2 (0) n ˜ r k c .

The Hessian is given by

2 f = 1 2 n ( H + v v T ) ,

where

H = diag j σ j (0) + σ j (1) σ k 3 ( e ) k , e and v = 1 σ 1 (0) , 1 σ 1 (1) , , 1 σ K (0) , 1 σ K (1) T .

We want to consider the eigenvalues of H + v v T . First, observe that at most one eigenvalue can be nonnegative. This follows from the famed Weyl Inequalities [50]. H has all strictly negative eigenvalues, while v v T , being an outer product, has one positive eigenvalue, v T v , with all other eigenvalues 0. Denoting as λ i ( G ) the ith largest eigenvalue of matrix G , the Weyl Inequalities tell us that

λ 2 ( H + v v T ) λ 1 ( H ) + λ 2 ( v v T ) = λ 1 ( H ) < 0 .

Hence, only one non-negative eigenvalue is possible.

Next, we can use the matrix determinant lemma to observe that

det ( H + v v T ) = ( 1 + v T H 1 v ) det ( H )

and direct computation tells us that

v T H 1 v = 1 .

Hence, the determinant is 0, meaning at least one of our eigenvalues must be zero. Combined with our prior result, this means our maximum eigenvalue must be zero and we conclude the Hessian is negative semidefinite. Thus, f is indeed concave.

Finally, note that the extension to the weighted case is straightforward. We can simply define new variables σ ˜ k ( e ) = w k σ k ( e ) for e { 0 , 1 } , and then repeat the proof above using the σ ˜ k ( e ) variables. Since σ k ( e ) is simply an affine transformation of σ ˜ k ( e ) , concavity in the former follows from concavity in the latter.

D Extension to cost-constrained case

The existing results can easily be extended to the case in which costs vary by treatment status, and we have a constraint in terms of cost rather than sample size. In this setting, we associate with each stratum k and treatment status { t , c } a constant c k e , e { t , c } , which represents differential cost per unit. With a budget constraint C , the regret minimization problem now becomes

(10) min n r k t , n r k c max σ k 2 ( 1 ) , σ k 2 ( 0 ) k w k σ k 2 ( 1 ) 1 n r k t 1 n r k t ˜ + σ k 2 ( 0 ) 1 n r k c 1 n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K k c k t n r k t + c k c n r k c = C .

Crucially, our constraints are still affine, so we can again invoke von Neumann’s minimax theorem to switch the order of the minimization and maximization, yielding the problem

max σ k 2 ( 1 ) , σ k 2 ( 0 ) min n r k t , n r k c k w k σ k 2 ( 1 ) 1 n r k t 1 n r k t ˜ + σ k 2 ( 0 ) 1 n r k c 1 n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K k c k t n r k t + c k c n r k c = C .

But the inner problem has an explicit solution, given by

n r k t = C w k σ k (1)/ c k t k w k ( c k t σ k (1) + c k c σ k (0) ) , n r k c = C w k σ k (0)/ c k c k w k ( c k t σ k (1) + c k c σ k (0) ) .

Plugging this in yields the simplified problem

(11) max σ k 2 ( 1 ) , σ k 2 ( 0 ) 1 C k w k ( c k t σ k ( 1 ) + c k c σ k ( 0 ) ) 2 k w k σ k 2 ( 1 ) n r k t ˜ + σ k 2 ( 0 ) n r k c ˜ subject to ( σ k 2 ( 1 ) , σ k 2 ( 0 ) ) A k ( Γ ) , k = 1 , , K .

Using the same logic as the final paragraph in Appendix C, we immediately see that Problem (11) is concave. Hence, we can use the same projected gradient descent approach to efficiently solve this problem.

E Results with rectangular confidence sets

We consider the case of using rectangular confidence sets, rather than ellipsoids. The implementation is straightforward: per the discussion in Section 4.2, we simply draw B bootstrap replicates and compute the associated extremal rectangles; then draw the minimum rectangular set that contains all of them; and then shrink the rectangle proportionally toward its center until a ( 1 α ) proportion of the rectangles have all four corners included in the set. Per the results in Appendix B, this procedure yields asymptotically valid α -level confidence sets.

We use this procedure on the WHI data and report the results in Tables A10 and A11. We see that results are substantively similar to those in Tables 1 and 2. In particular, all the entries in Table A10 are non-positive. The rectangular confidence sets actually seem to yield slightly larger performance improvements in the case of stratification on a single variable (yielding fewer total strata). However, the performance is slightly poorer when stratifying on multiple variables.

Table A10

L 2 loss comparisons for regret-minimizing allocations versus equal allocation, rectangular confidence regions

Subgroup Var(s) Equal alloc loss Loss relative to equal allocation
Γ = 1 (%) Γ = 1.1 (%) Γ = 1.5 (%) Γ = 2 (%)
Age 0.000691 3.8 3.2 3.5 2.3
CVD 0.000649 2.4 2.5 1.3 1.4