Integer Programming for Generalized Causal Bootstrap Designs
Authors: Jennifer Rogers Brennan
Sebastien Lahaie
Adel Javanmard
Nick Doudchenko
Jean Pouget-Abadie
Venue: ICML 2025 poster (2025)
Keywords: design uncertainty
causal bootstrap
integer programming
Primary area: general_machine_learning->causality
TL;DR: We propose an integer program for computing the worst-case copula used in the causal bootstrap, which improves the estimation of design uncertainty in experimental causal inference, generalizing the settings in which this bootstrap applies.

Abstract: In experimental causal inference, we distinguish between two sources of uncertainty: design uncertainty, due to the treatment assignment mechanism, and sampling uncertainty, when the sample is drawn from a super-population. This distinction matters in settings with small fixed samples and heterogeneous treatment effects, as in geographical experiments. The standard bootstrap procedure most often used by practitioners primarily estimates sampling uncertainty, and the causal bootstrap procedure, which accounts for design uncertainty, was developed for the completely randomized design and the difference-in-means estimator, whereas non-standard designs and estimators are often used in these low-power regimes. We address this gap by proposing an integer program which computes numerically the worst-case copula used as an input to the causal bootstrap method, in a wide range of settings. Specifically, we prove the asymptotic validity of our approach for unconfounded, conditionally unconfounded, and and individualistic with bounded confoundedness assignments, as well as generalizing to any linear-in-treatment and quadratic-in-treatment estimators. We demonstrate the refined confidence intervals achieved through simulations of small geographical experiments.

Integer Programming for Generalized Causal Bootstrap Designs
Jennifer Brennan'
Sébastien Lahaie'
Adel Javanmard
12 Nick Doudchenko 3 Jean Pouget-Abadie
Abstract
captures the uncertainty associated with inferring statisties of population given random sample of that population. Design uncertainty, by contrast, captures the uncertainty associated with the random assignment of units to treat- ment and control conditions. Design uncertainty is associ- ated with the finite population model of causal inference, in which the universe consists of exactly N units. Some units are observed under the treatment condition while oth- ers are observed under the control. The objective in the finite population model is to estimate the average difference between treated outcomes and control outcomes for these N units, typically in the absence of distributional assumptions on these outcomes. Uncertainty arises from the fact that each unit is only observed under the treated or controlled condition, so that any statistic derived from RCT data will vary slightly from one randomization to another.
In experimental causal inference, we distinguish between two sources of uncertainty: design un- certainty, due to the treatment assignment mecha- nism, and sampling uncertainty, when the sample is drawn from a popplaaooot This distinc- tion matters settings with small fixed samples and heterogeneous treatment effects, as in geo- graphical experiments. The standard bootstrap procedure most often used by practitioners primar- ily estimates sampling uncertainty, and the causal bootstrap procedure, which accounts for design uncertainty, was developed for the completely ran- domized design and the difference-in-men in-mmns es- timator, whereas non-standard designs and esti- mators often used in these low-power regimes. We address this gap by proposing integer pro- gram which computes numerically the worst-case copula used as an input to the causal bootstrap method, in a wide range of settings. Specifically, we prove the asymptotic validity of our approach for unconfounded, conditionally unconfounded. and individualistic with bounded confoundedness assignments, as well as generalizing to any linear- in-treatment and quadratic-in tuaeraeicii-trraatment treatment estima- tors. We demonstrate the refined confidence inter- vals achieved through simulations of small geo- graphical experiments.
Sampling uncertainty has received extensive attention in the statisties literature. Well-known methods to estimate sampling uncertainty include the bootstrap (Efron, 1979), in which random draws from the sample mimic the popula- tion distribution, and Wald-type confidence intervals, which rely on asymptotic normality of the sample statistic. The methodological toolkit available quantify design uncer- tainty is considerably more sparse (Imbens & Menzel, 2021; Aronow et al., 2014). For reasons explained in Appendix 2.1, methods to estimate sampling uncertainty are typically conservative for design uncertainty. The gap widens for experiments with heterogeneous effects, small sample sizes, and a large fraction of the population included in the exper- iment. Even fewer methods apply to general designs such those that balance covariates of the treated and control groups, which are widely used precisely in the settings with small experimental samples and heterogeneous treatment effects (Harshaw et al., 2024).
1. Introduction
Randomized controlled trials (RCTs) are used to estimate causal effects across the sciences, social sciences, and indus- try. An important complement to estimating causal effects providing valid uncertainty quantification, which enables statistically sound maasinnnmaking between tested policies. We distinguish between two types of uncertainty: sampling uncertainty and design uncertuinty. Sumpling uncertainty
In this paper, we propose novel application of integer programming to identify the joint potential outcome distri- bution that maximizes the variance of the chosen estimator while being consistent with the randomization design and the observed marginal distributions of potential outcomes. Solving for this joint distribution numerically is what allows us to extend prior work (Aronow et al., 2014; Imbens & Menzel, 2021) which relies on well-known results of the optimality of the isotone copula for the difference-in-meet means estimator and completely randomized assignment. In partic-
'Google Research Uniiversity of Southem California 'Meta (work completed while at Google Rescarch). Correspondence to: Jean Pouget-Abadic <jeanpa @ google.com>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Integer Programming for Generalized Causal Bootstrap Designs
ular, our method applies all unconfounded (and known and probabilistic) assignment mechanisms, as well as any individualistic (and known and probabilistic) assignment mechanisms that verifies some bounded confoundedness. We examine conditionally unconfounded assignments as a special case. Furthermore, our method applies to any linear-in-treatment estimator, i.e. estimators of the form Z,a, + bi, and quadratic-in-treant inttretmmettt estimators, i.c. es- timators of the form E; by + E; ZzZ.ja.... which include among others, difference-in-means. linear regression on the treatment variable and covariate, and any doubly robust estimator fit out-of-sample. We provide asymptotic valid- ity results and illustrate the application of our method to a simulated geographical experiment, using real data from the International Monetary Fund (IMF). Our method provides similar coverage and generally tighter confidence intervals across variety of settings when compared to a widely used, but more conservative baseline.
et al. (2020) propose a novel variance estimator for linear regression under both design and sampling uncertainty that trades off between a classical ordinary least squares variance estimator and a novel design-based based estimator based on the fraction of the population that is sampled.
Other studies expand design uncertainty estimation beyond Bernoulli or completely randomized designs. Dasgupta et al. (2015) and Lu (2016) study design uncertainty under factorial designs. Imai (2008) extends Neyman's 's upper bound to the matched-pain design, which is a special case of covariate- rovariate-balancing designs; Fogarty (2018) further extends this result to regression adjustment for matched pairs.
Neyman's 's original variance estimator is unbiased only for homogeneous ("strictly additive") treatment effects. Muker- jee al. (2018) identify less restrictive set of constraints the potential outcomes and experimental design under which the variance of certain estimators is unbiasedly es- timable. This approach relies on a "Q decomposition" of the variance, where the matrix Q must be derived for each experimental design. The approach fails for certain designs including matched pairs, where units in the same pair cannot be assigned the same treatment. Chattopadhyay & Imbens (2024) introduce two new variance decompositions for gen- eral designs and illustrate their application to completely randomized designs.
In Section 2, we define and introduce our method for the il- lustrative setting of the completely randomized assignment and the difference-in-menn -n-e-nnn--seseeeee estimator. In Section 3, we show that the same approach can be generalized to linear- and quadratic-in-treatmee show our method applies and is asymptotically valid for unconfounded assignments, conditionally unconfounded assignments. and individualistic with bounded confounded- ness assignments. We conclude in Section 5 with practical simulations on real data.
estimators. In Section 4, we
Like in (Aronow et al., 2014), our approach relies on iden- tifying the mrxiince-aiximizing coupling of treated and control observations. We then apply the causal bootstrap method of Imbens & Rubin (2015) generate confidence intervals. Our main contribution to these works is to gener- alize to new estimators and assignment mechanisms through the formulation of an integer program, which no longer re- quires the worst-case coupling to be known in closed-form. In a related optimization framework, Ji et al. (2023) show how to estimate bounds on a functional of the joint distri- bution via optimal transport. They do not explicitly solve for the optimal coupling, but their dual approach is robust to misspecification of marginal treated and control outcome probabilities in the presence of covariates. Finally, Harshaw et al. (2021) propose a minimal variance bound among all quadratic forms through a convex optimization procedure. Their method yields interesting admissibility results but only considers quadratic bounds. Our approach, inspired by the Frechet-Hoeffding bound approach used by Aronow et al. (2014) and Imbens & Menzel (2021), avoids this limitation.
1.1. Related works
The quantification of design uncertainty traces its roots to Neyman (1 (1923). who studied the difference-ing -innnennssss-sss esti- mator under complete randomization in the finite population setting. Later work identified tight upper bounds for bi- nary (Robins, 1988), interval (Stoye, 2010), and continuous (Aronow et al., 2014) outcomes. Imbens & Menzel (2021) extend the method of Aronow et al. (2014) with a resam- pling procedure they call the causal bootstrap, providing a more practical confidence interval construction.
Tighter variance estimates can be obtained when predic- tions are paired with informative covariates. Ding et al. (2019) and Wang et al. (2020) extend the work of Aronow et al. (2014) to the difference-in-nee in-means estimator in the pres- ence of covariates, while Fan & Park (2010) propose the fewer-than-n bootstrap method in the presence of covariates. Another body of literature studies design uncertainty for linear regression estimators. Freedman (2008), Schochet (2010), Samii & Aronow (2012), and Lin (2013) character- ize the contexts in which classical sampling-basee based variance estimators are valid in the design setting, while Abadie et al. (2014) and Fogarty (2018) propose regression variance es- timators specific to the mnite-population setting. Abadie
2. Our causal bootstrap procedure
2.1. Notation and background
To build intuition, we briefly review Neyman's conservative bound and follow-up work on the difference-in-men mn-means esti- mator under a completely randomized assignment. We con-
2

---

Integer Programming for Generalized Causal Bootstrap Designs
sider a finite-population potential outcomes model in which each unit i E [N] is assigned to treatment (Zz 0) con- trol (Z,- = 1), and outcome Y,.))) is observed. The average treatment effect is defined as N EvE[y] (1) - (0). Let N1 and No be the number of treated and controlled units respectively. Define the difference in means estimator ADIM - N. Eiz No
method proposed by Imbens & Menzel (2021): we impute unobserved potential outcomes using the least-favorable copula obtained above and drawing bootstrap samples via rerandomization of umits to treatment. Confidence intervals are then constructed from the quantiles of this causal boot- strap distribution. We refer the reader to a full discussion in Appendix A.9.
-i Yi
and the completely randomized design that assigns exactly N1 units to treatment and No units to control. Neyman (1923) decomposed the variance of this estimator as
We propose to solve the least-favorable copula problem numerically, by formulating the optimization problem in Eq. (2) in solver-friendly way. To build intuition, we first consider the effference-in-mean estimator. We will explore generalizations to a wider class of estimators in Section 3.1 and 3.2. The variance of the difference-in- in-mnasss estima- tor has a quadratic closed-form for general assignments: Var z [+] {Y) (0)) + Y Y. where Y is the vector of coordinates LNN'YY(() + +No Y:()))=1..... and Σzz is the covariance matrix of the random vector (Z.)i=....,,, with coordinates (Ezz)aj = Cov(Z,, Zz). While it may be tempting to use this formulation directly, Σzz is positive semi-definite, making the objective straight- forward to minimize, but more difficult to maximize under constraints, which is our goal (2).
St N1
S
Varz [fDIM]
(1)
No
N
where S is the sample variances of the treated (= = 1) and control (z = 0) units, while S? is the sample variance of the unit-level treatment effects ((1)-Y((0) See Appendix A.8 for the full formulas. S depends on the unknown joint distribution between treated and control outcomes since it contains terms ())) making it impossible to estimate consistently from experimental data in which only one of (1) or (0) is ever observed. While Neyman used the fact that variances
nonnegative to conservatively bound S 2 O, thus proposing an upper-bound of Varz [fDIM]. this bound can be improved by identifying joint distribution that maximizes S? while matching the observed marginal distributions of (0) and (1). Aronow et al. (2014) iden- tify this aniaace-maximizing distribution to be the isotone (assortative) copula, which pairs the smallest treated out- come with the smallest control outcome, then the second- smallest, and so on. They derive analytical variance bounds using that joint distribution, while Imbens & Menzel (2021) propose bootstrap-style approach, imputing the counter- factual outcomes using the isotone copula and resampling those imputed outcomes to generate bootstrap samples.
Instead, propose the following formulation: consider the indicator variables x(ee) with a E {00 1} indexing the treat- ment assignment, i N] indexing the units in the popula- tion, and yv the k outcome of y:: x-(a := I(Y((a) = yye) sk
In particular, this notation implies that y is discrete. For continuous outcomes, y can be made discrete by replac- ing it with the support of observed outcomes. We provide asymptotic validity results of this approach in Section 4.2, though nothing prevents practitioners from expanding this discrete support of sitcoomes-byyone computational scala- bility, which we cover in Section 5-to cover other unseen outcomes. We can rewrite the expression of the variance of the Lufference-in-mean estimator under general assign- ments as a quadratic form:
2.2. Our approach
In this work, we propose to identify the variance- maximizing, joint distribution numerically using optimiza- tion instead of analytical derivations, allowing us extend the least favorable copula approach to a broader class of estimators and assignments. Our goal is to find the joint distribution of potential outcomes Fo, which maximizes the variance of some estimator f under some assignment mechanism Z and with some constraint C on the shape of this joint distribution.
Var X Y ExT'zzY Yx x XTQX,
(3)
where X E (0, and Y E Nx(N-K2 K) are a vector and matrix representation of the indicator variables and po- tential outcomes respectively; ef. Appendix A.3 for details.
2.3. Adding constraints on the potential outcomes
An important part of our proposed algorithm is imposing realistic constraints on the missing potential outcomes that do not jeopardize the asymptotic validity of our bound while reducing the value of our objective (2), hence making our bound tighter. We constrain each unit i's potential outcome (a() to only match one possible value yk. leading to the
max YO
Varz|f]
(2)
S.L.
Fo,1 CC
There are multiple ways of constructing confidence inter- vals from the solution to this optimization objective. We propose to follow a similar procedure to the causal bootstrap
constraint Σ₁ x(a) 1. Additionally, because we observe
3

---

Integer Programming for Generalized Causal Bootstrap Designs
N of the 2N potential outcomes, we set x(2)) and only if y,oha yk-
for which powerful solvers available, using standard transformations; details are provided in Appendix A.3.
A clever observation by Aronow al. (2014) is that the observed potential outcomes allow us to impose further constraints on the missing potential outcomes. While they focus entirely on the completely randomized assignment, which we use now to build intuition, we will generalize our results to other assignments in Section 4.1. I. If all units receive treatment with the same probability (different from 1 or 0), we expect the treated observed Fiehs and missing Fymis marginals to match and the controlled observed Fole and missing Fo marginals to match. Formally, for any outcome y and treatment assignment value a E (0, 1 con- sider
max X"QX for X {0, 1] 1
(4)
Vi
X (2)) x(a
1 iff Y yk
(a)
X ((a) = 0 iff yn supp(Fa)
Ya, i,
(b)
Va, i,
Σ x
(c)
k=1
Va, k,
N
(a)
ik
Ni
Zi
1
Ze No
=
(d)
1=1
Constraint (4.b) is optional and allows us to account for prac- tical assumptions we may have on the support of Fa. For example, we may be able to assume that Fa C Fahs, further reducing our estimate. Constraint (4.d) is not always feasi- ble, when taken in conjunction with the other constraints on (x(e). To ensure feasibility, we can relax these constraints for some chosen slackness parameter e > 0,
N Σ
Fobs (v)
- Yi(a) =))
= Na i=1
Fmmis
Σ 1(2 11 a, (a) y)
1=1
For a completely randomized assigmment, each observed and missing marginal is equal in expectation to the true un- derlying marginal Fa- In other words, have the following equality for any potential outcome y::
Va E {0, 1), Vk,
N
Z N1
( ik
1 1- Zi ((()) No
S..
The regime of feasibility is captured in the following Lemma, a proof of which can be found in Appendix A.2.
Eg FF (w)
Σ I(Yi(a)= = y) =)=Ez|FF" E2(FFo (y)l
Lemma 2.1. cor O, there may be solution to the constrained optimization below (4), unless No Ni. It admits a feasible solution for e 2 1/11mmmmm mio
F (a)
In our chosen parameterization, this equality on the marginal distributions can be written as a linear constraint on the variables (x (x(a)) sk
As decreases, the worst case bound on the variance de- creases, improving the performance of our estimator. We include further discussions on the choice of e in Section 4.2.
Va E 10, 1J,
Z
x
((a) 2k:
No
3. Generalizing to new estimators
3.1. Linear-in-treatment estimators
2.4. Our proposed integer linear program
So far, our results hold only for the difference-in means estimator, which has been the primary object of study of previous work. In f'ac,, they can be easily extended to slightly more general class of estimators of the form: f :=
-neamss
We now rewrite our objective in Eq. (2) as integer pro- gram. For the completely randomized assignment and the difference-in-meant estimator, we know from prior work that its solution is the isotone copula. The value of the inte- ger program formulation will become clear as we generalize beyond the completely randomized assignment, and be- yond the lffference-in-mest closed-form solutions are not always known.
Σ Z,a, + b, where a, and b; are constants with respect to the assignment mechanism and linear in the potential outcomes. We refer to such estimators as being linear-in- treatment (the second linearity assumption is then implied, for simplicity). Whether an estimator linear-in-treatmen depends on the assignment mechanism. Strictly speaking, the difference-in-meann estimator is linear-in-treateen for completely randomized assignment but not for a Bernoulli assignment, since the normalization constants N and No are not constants with respect to Z. The Horvitz- Thompson estimator (Horvitz Thompson, 1952) is linear-
setting for which
When indicated, the equations below are defined for all E (0. 1). for all i E [1, N). and for all yk V. As written, the optimization has binary variables and its constraints are all linear, but the objective is quadratic. Because the product of two binary variables also a binary variable, the objective can be rewritten as a linear function of binary variables,

---

Integer Programming for Generalized Causal Bootstrap Designs
in-treatment for either assignment mechanism:
pre-freatment-period value of the outcome as a plug-in estimator (Ahadie, 2005; Deng al., 2013). If we expect growth of this variable, we can also account for its rate of growth and rescale the plug-in estimator with its expected value; however, fitting this growth rate using data from the experimental period depends on the treatment assign- ment and should be avoided. We include an example of this approach in Section 5. Another option is to rely on sample-spliting, where split units into C' cohorts, cho- sen independently of the treatment assignment Z. In that case, we can estimate a (pair of) function(s) per cohort using data from all other cohorts.
(1)1)
+ Yi(0)
- (0:
f Σ Zi
N P N (1 (1 )P)
+ N- (1 P.)
where P. = [ = 1) is the probability that unit i receives treatment. For any assignment mechanism and linear-in- treatment estimator, Varz [+] = Dig a:ag Cov Zi, Zn] = Q where Q' = TUTTZZUU, with Y and X defined similarly to Section 2.4, and U such that = UYX, with a the vector of coordinates bel)e The existence of U is given by the linearity of a, in the potential outcomes. The constraints in our integer program in Eq. 4 do not incor- porate any estimator information, such that we can easily generalize our results to any -ntraatteettttmmettt estimator by replacing Q with Q'
3.3. Quadratic-in-treatmee estimators
Repeating a similar approach to Section 3.1 for estimators of the form f = >; 6y + Σ; Z,Z, we obtain the following quadratic form in the coefficients aıj when computing their variance: Varz[f] - Enj.e. (ijuarCov [Z,Z,, Zi We can then repeat the approach of Section 3.1 as long as the coefficients dg continue to be constant with respect to the treatment assignment and linear in the potential out- comes. In particular, we show in Appendix A.6 that linear regression fitting observed outcome Y to the treatment variable and a scalar covariate X can be written in this form, with coefficients dij that are linear in the potential outcomes. One advantage of this approach is that we are no longer required to do sample-spliting or use solely pre- period information to incorporate covariate information and reduce the size of our confidence intervals.
3.2. Fitting regressions out-of-sample
In most practical applications, we also have access to covari- ates, which improve the precision of effect estimates when they are predictive of the outcomes. In this section and the next, we consider ways to extend our integer program to estimators that include covariate information.
One straightforward way to extend our previous analysis is to consider aatcome-to-covariate functions that are fit out-of-sample. Consider the doubly-robust estimator as illustrative example, where we replace the propensity scores by the true probabilities of each unit being treated, since we are motivated by randomized experiments where the probabilities are known. Let fi and fo be the predicted outcomes using functions that are fit out-of-sampil samls and do not depend Z, and E Y; fi (W,) and Eo- == Y, (W,) the corresponding residuals, such that we write F: = N I + Zá6 P. fo..))
4. Generalizing to new mechanisms
We now show that our integer programming formulation in Eq. 4 can be generalized to a wider class of probabilistic as- signment mechanisms for any of the estimators in Section 3.
(1-2,66) I-P. T-P
4.1. Known and Probabilitistic assignments
The variance of this estimator is equal to the variance of the Morvitz-Thomppe
estimator applied to the residuals Y < Y 2.f.(W)) - (1 (W.)).W.. When taking the variance of with respect to and replacing Y! +
Of the constraints in Eq. requiring generalization, only (4.d) is restricted assignments having P. = p Vi. In the more general setting of all known and probabilistic as- signment, we no longer expect Fini and Fobs to match exactly-or even approximately. If the treatment proba- bilities are known however, we can construct an unbiased estimator of the true marginal F, from both observed and missing potential outcomes. For a E {0, 1} and yk E V, consider the marginal estimators:
Y, Z.f,(W.)- - (1 (W.)). we have:
N Σ
Z.Y' P.
1
N
(1-
(1-
(1
Varz(f] = Varz
N
N
1-P.
We recognize the Horvitz- Thommpsn estimator applied to the residuals Y/ We can apply our method to Y + Y/. In particular, the constraints on the marginal distributions now constrain the residual marginal distributions to match.
1 71 - 7,) 1-a (a) and
Fobe robb((x))) == N PO P.) 1-a
Σ
There are many practical settings in which fitted-out-of-- sample functions are available. For example, such functions may be fitted on the same cohort before the start of the ex- periment. A typical and straightforward example to use
1 Z 2,)0 (0)
fomis
famis(ya) (yk) == N
p "(1 P.)a
5

---

Integer Programming for Generalized Causal Bootstrap Designs
Through iterated expectations, it holds that El E (yne)] = El R[ [y+)]. We can thus reformulate constraints (d) as follows, for all a, k. b E {0, x |y| x 10,
across E [n], and suppose that
PiZ = a (0) a)
cd.
N (-1)
b
( sk
Z
1 1- Z
1 - P. SeN.
for some 8 0, such that 5 controls the effect of confound- edness, with d 0 corresponding unconfoundedness. If the slackness satisfies c 2 6, then
P;
The resulting integer program applies to any known, and probabilistic assignment mechanism.
P(V* > Varz|f)> 2 1 3,
where 3 ;= 11 - 8e 2 -8 -NP2/// 5D 240+c
(c-8)P(1+8)
6-6 and E1 = 240+c*
4.2. Asymptotic Validity
Just because can formulate (and solve) an integer pro- gram for any known and probabilistic assignment mech- anism, does not mean our bound is asymptotically valid. Our procedure is guaranteed to upper-bound the true vari- ance only if the true joint distribution is also a feasible solution of our integer program. This could be violated if (a) the discrete support we choose, equal to the support of observed outcomes, does not also encompass the support of the missing potential outcomes, or (b) the realized dis- tributions Fimss and Fobs do not match exactly. Recall that we proposed relaxing constraint (4.d) to hold with slack- ness to ensure feasibility of our integer program. In fact, by bounding the probability that the observed and missing marginal distributions are uniformly within an e distance of each other, we can provide rate for the probability that our procedure upper-bounds the true variance as N grows:
Details can be found in Appendix A.5. We leave the gener- alization to general confounded assignments to future work.
4.3. Conditionally unconfounded assignments
It is possible for the imbalance in Fimss and Fobe to be fully captured by an observed covariate W,, as is the case for conditionally unconfounded assignments, which verify Z, 1 (1)), ((()) In that case, we have at least two op- tions for formulating the integer program. The first option, and the most straightforward one, to re-use the integer pro- gram from Section 4.1 with P. P(Z, WW; = W.). While our previous results follow through trivially, this approach does not leverage all the covariate information available.
A second approach that leverages the covariate information is to impose equality of the conditional marginal distribu- tions, reweighted by the suitable inverse probabilities. In other words, for all a E {0, 1).1,.. yk y and values of covariate W, consider the conditional marginal estimators:
Theorem 4.1. For any probabilistic unconfounded assign- ment mechanism, P (V* 2 Var (f)z(())> 2 1 1-3, where V* is the output of our integer program, / = 8 exp (
(1- 2,)²=
x ((a)
N= ENJE(NI Cov(Z,, Zj). is the slackness parameter in (4.d). P == min(P. P), and P N E.
frabs (yu (u]ww)= == Σ
::W: Pı(w)a(1- - P. P.(w)) [w]
2 (1- 2.
X(
As expected, as the sample size N and the slackness e that we choose grows, the probability of upper-boundie the true variance increases. We can therefore justify our approach theoretically for any unconfounded (and known and prob- abilitistic) assignments. We will examine the special case of conditionally unconfounded assignments in Section 4.3. We can further refine the result above if the assignment mechanism is also individualistic:
Ermie (yk|w)= =
i: = P.(w) 1-a Pi(ww)) KwW
Here, P = 1|W, = w) is the probability that unit receives treatment Z, = a given covariate W; = w. and [wH denotes the number of units with covariate W;- = w. As a result, the following equalities hold for all (a, Mk, w) triplets: Ez Foobe (vxlw)) = Ez famis (yx (u))wr))
Corollary 4.2. For any probabilistic, unconfounded, and individualistic assignment mechanism. Theorem 4.1 holds with 3 = 8 8exp f (N)) + 8 8xxp (
These equalities can be expressed in solver-friendly way. similar to the formulation of the (4.d) constraints in Sec- tion 1. Doing so allows us to incorporate covariate infor- mation into our variance estimator for tighter confidence intervals, since it adds additional constraints the optimiza- tion problem. On the other hand, the empirical conditional marginal distribution F. (2)(ww) converges to its expectation more slowly than the unconditional marginal, which will negatively affect the performance of our variance bound in small samples. This is especially true when W is high- dimensional or continuous.
What of confounded assignments? In fact, we can also prove our approach is asymptotically valid if the assignment is individualistic with bounded confoundedness, as defined by (Tan, 2006; De Bartolomeis et al., 2024). Adopting similar notation to Theorem 4.1,
Theorem 4.3. Suppose that (Y:(0), Y((1), Z,) are i.dd
6

---

Integer Programming for Generalized Causal Bootstrap Designs
5. Simulations on real data
tor distribution. We also compare to the causal bootstrap method suggested by (Aronow et al., 2014; Imbens & Men- zel, 2021), which resamples from the isotone copula. This method is only theoretically justified in the complete ran- domization setting for the difference-in-means estimator. All reported confidence intervals have 95% nominal cov- erage. We also report the "true confidence intervals" as defined by the 2.5% and 97.5% quantiles of the distribution of the estimate when resampling the assignment according to each design 500 times.
5.1. Empirical Results
Our approach is well-motivated for geographical experi- ments, which assign treatment and control to, e.g. countries or Designated Marketing Areas. This is due to the fact that geographical regions are fixed in number, meaning design uncertainty likely dominates sampling uncertainty. Addi- tionally, there are often few available regions (on the order of hundreds) with significant heterogeneity in potential out- comes. As a result, (a) sample uncertainty quantification methods are poor approximations of the design uncertainty of our estimator, and (b) practitioners rely on non-standard designs estimators to achieve statistical power. This is precisely the setting in which existing design uncertainty quantification methods fall short.
A Analysis more complete of results can be found in Appendix A. including implementation details and reference to the open- source solver. For any experiment with complete random- ization, our optimal causal bootstrap mechanism and the isotone copula bootstrap yield the same results down to each simulation as expected: our optimal causal bootstrap recovers the well-known suriance-maximiziin joint distribu- tion exactly. Across all three inference methods and for the ground truth confidence intervals (CIs), the CIs are greatly reduced by incorporating pre-period information. These CIs are further reduced by the matched -piiss design, with more discussion in Appendix A.1
Select results are summarized in Tables 1-2. A
Dataset To work with realistic dataset, we look to the International Monetary Fund's publicly-avaiable Gross Do- mestic Product (GDP) country-level report for the years 2017-2019, restricted to the top 50 countries by GDP (IMF, 2024). While a coordinated worldwide randomized experi- ment affecting GDP is implausible, we can make the more reasonable assumption that some global company's 's metric- of-interest in each country is proportional to that country's s GDP, such that the simulation results would approximately hold if such a company were to run a country-level ex- periment. Our outcome of interest is GDP in 2019, with 2017-2018 GDP acting as pre-period covariates.
All three bootstrap methods achieve lower-than-nomee coverage for the completely randomized design, due the heavy-tailed nature of GDP data, which causes F. to con- verge slowly. This problem is particularly acute for the completely randomized design, which has high likelihood of assigning the largest GDP countries to the same treatment group. By contrast, the matched-pail design balances the treated and controlled outcomes, improving the convergence of Fa and thus the coverage of the sampling bootstrap and our method. We provide more discussion of the observed undercoverage in Section A.
Settings To showcase the validity and usefulness of our method, we compare two assignment strategies: the stan- dard complete randomization design ("Compl. R.") treat- ing exactly half of the countries, and a matched-pairs de- sign ("Matched Pairs"). The matched-pairs design is con- structed by matching the country with the largest GDP in 2018 with the country with the second largest GDP for that year, the third largest with the fourth largest, and so on. We consider two estimators: first, the standard "no- covariates" difference-in- means estimator, and second, the doubly-robus estimator from Section 3.2, using pre-period covariates from 2018 GDP adjusted for growth using the historical change between 2017 and 2018, as described in Section 3.2.
The conservative variance estimators fare well at achieving dlose-to-nominll coverage. For the doubly robust estimator under matched-paire -aiidd--aisshhee the setup which achieves the best power here- the conservative estimator has the narrowest CI widths for small effect sizes (less than 2.5%), but for large effect size (over 10%), its CI widths become larger than those of the sampling bootstrap. Our optimal causal bootstrap maintains a stable of improvement of at least 11% in CI widths relative to the sampling bootstrap, across effects sizes. We emphasize that analytical variance bounds are not typically available for more complex designs, whereas our procedure extends as long as the treatment covariance matrix is available.
Basclines Finally, we compare our method to three base- lines. Difference-in-ment (and its doubly-robust variant) admits conservative, analytical bounds for its sampling vari- ance under both complete randomization and matched-pain see Appendix A.7. We form symmetric, Neyman-style con- fidence intervals based on these bounds. We next consider the standard sampling bootstrap, which samples with re- placement N units 1, 000 times, and reports the empirical confidence intervals as quantiles of the resulting estima-
For the matched-pairs design, the isotone copula bootstrap is no longer known to be conservative estimator. In prac- tice, we find its resulting Cls are unrealistically narrow: in Table 1, they never intersect with zero, leading to a coverage of 0%. This is not surprising in a well-matched A/A setting
7

---

Integer Programming for Generalized Causal Bootstrap Designs
Covariates
True
Inference Method
Coverage
Median CI
Compl.
Matched Pairs
Compl. R.
Matched Pairs
Compl.
Matched Pairs
Sampling Boot.
90%
100%
3967
4110
Conservative Var.
98%
100%
3991
1152
Nome
3574
762
Isotone Copula
87%
0%
3348
44
Opt. Causal Boot.
87%
100%
3348
2233
Sampling Boot.
91%
100%
141
142
Conservative Var.
94%
95%
143
49
GDP
142
46
Isotone Copula
90%
91%
133
41
Opt. Causal Boot.
90%
100%
133
102
Table 1: "True CT" and "Median CT" correspond to 95% CI widths in an A/A test under two designs and two estimators.
Covariates
Effect
True CI
Inference Method
Median CI
Power
Compl. R.
Matched Pairs
Compl. R.
Matched Pairs
Compl. R.
Matched Pairs
Sampling Boot.
4015
4160
0.10
0.00
Conservative Var.
4040
1215
0,02
0.00
2.5%
3618
772
Isotone Copula
3395
0.13
Opt. Causal Boot.
3395
2256
0.13
0.00
None
Sampling Boot.
4181
4306
0.10
0.00
Conservative Var.
4149
1464
0.02
0.00
10%
3752
800
Isotone Copula
3540
0.13
Opt. Causal Boot.
3540
2341
0.13
0.00
Sampling Boot.
111
109
0.32
0.20
Conservative Var.
117
84
0.24
0.46
2.5%
101
45
Isotone Copula
99
0.39
Opt. Causal Boot.
99
97
0.39
0.25
2018 GDP
Sampling Boot.
185
1.00
1.00
Conservative Var.
189
310
1.00
0.93
10%
90
64
Isotone Copula
154
1.00
Opt. Causal Boot.
154
203
1.00
1.00
Table 2: "True CI" and "Median CT" correspond to 95% CI widths for multiplicative effects under two designs and two estimators.
where the difference of outcomes across pairs larger than the difference within pairs. In that case, the isotone copula imputes the counterfactual outcome of each unit with the observed outcome of its paired unit. No matter how many times one samples from this imputed copula, the estimated effect within each pair will be the same, leading to an es- timated variance of zero. However, the point estimate of our estimator can be non zerze even for well-matee well-matched A/A test, if the paired units not exactly identical, leading zero coverage. As result, we do not report the median CIs of the isotone copula bootstrap in Table for the matched- pairs design. On the other hand, both the sampling bootstrap and our optimal causal bootstrap have satisfactory coverage both complete randomization and matched pair designs. Across experiments, our method reduces CI widths by at least 10% compared the naive bootstrap, with no loss in power coverage, with an exceptional 50% reduction in CI widths for the difference-in-meant estimator with injected multiplicative effects.
5.2. Scalability Analysis
As mentioned, we focus on the top 50 countries by GDP for our simulation analysis. For dataset of this size, our integer program took 671 seconds average for the matched pairs design, and 1057 seconds on average for complete random- ization. Matched pairs is faster because of the sparsity of its covariance matrix. However, complete randomization also admits a linear programming formulation for the worst-case copula, which took only seconds on average to solve. We used this faster implementation in our experiments for the sake of expedience.
To assess the scalability of our formulation, we solved for the worst-case copula under null effect, varying the num- her of units. Figures 1-2 presents the results of our scaling analysis, where each data point in the plot represents an av- erage over 10 simulation runs, along with 95% confidence intervals based standard errors.
8

---

Integer Programming for Generalized Causal Bootstrap Designs
180
Averppn
123
Figure Runtime scaling for matched pairs, integer pro- gram formulation.
Figure 2: Runtime scaling for complete randomization, linear program formulation.
6. Conclusion
2005.
We propose integer program approach design- uncertainty quantification for linear- and quadratic-in- treatment estimators for known and probabilistic assign- ments, that are unconfounded, conditionally unconfounded, or individualistic with bounded confoundedness. Future work could look at new classes of estimators or provide asymptotic validity results for general confounded mecha- nisms. There is also the possibility of unknown probabilis- tic conditionally unconfounded assignmentsm, in which case the probabilities P. would be replaced with estimated propensity scores. The validity of our approach would then depend on the quality of these estimates. The limitations of our work are around computational scaling, which is lim- ited to tutbut also best motivated for-small datasets with hundreds of units.
Abadie, A., Imbens, G. W., and Zheng, F. Inference for misspecified models with fixed regressors. Journal of the American Statistical Association, 109(508):1601-1114 2014.
Abadie, A., Athey, S., Imbens, G. W., and Wooldridge, J. M. Sampling-based versus design- based uncertainty in regression analysis. Econometrica, 88(1) (1:1)))2229999) 2020.
Aronow, P. M., Green, D. P., and Lee, D. K. Sharp bounds on the variance in randomized experiments. The Annals Statistics, 42(3):850-877 2014.
Chattopadhyay, A. and Imbens, G. W. Neymanian in- ference in randomized experiments. arXiv preprint arXiv:2409. 12498, 2024.
Impact Statement
Dasgupta, T., Pillai, N. S., and Rubin, D. B. Causal in- ference from 2k factorial designs by using potential out- comes. Journal of the Royal Statistical Society Series B: Statistical Methodology, 77(4):727-733, 2015.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences our work, none of which we feel must specifically highlighted here.
De Bartolomeis, P., Martinez, J. A., Donhauser, K., and Yang, F. Hidden yet quantifiable: lower bound for confounding strength using randomized trials. Interna- tional Conference on Artificial Intelligence and Statistics. pp. 1045-105 PMLR, 2024.
Acknowledgments
Adel Javanmard is supported in part by the NSF Award DMS-231 1024, the Sloan fellowship in Mathematics, an Adobe Faculty Research Award, an Amazon Faculty Re- search Award, and an iORB grant form USC Marshall School Business. The authors are grateful anonymous reviewers for their feedback improving this paper.
Deng, A., Xu, Y., Kohavi, R., and Walker, T. Improving the sensitivity of online controlled experiments by utilizing pre-experiment data. In Proceedings of the sixth ACM international conference on Web search and data mining, pp. 123-132, 2013.
References
Ding, P., Feller, A., and Miratrix, L. Decomposing treat- ment effect variation. Journal of the American Statistical Association, 114(525):004-317 2019.
Abadie, A. Semiparametric difference-in-ddifeeennss es- timators. The review of economic studies, 72(1):1-19,

---

Integer Programming for Generalized Causal Bootstrap Designs
Efron, B. Bootstrap methods: Another look the jackknife. The Annals of Statistics, 7(1):1-26, 1979.
Mukerjee, R., Dasgupta, T., and Rubin, D. B. Using stan- dard tools from finite population sampling improve causal inference for complex experiments. Journal of the American Statistical Association, 113(522):868-888 881, 2018.
Fan, Y. and Park, S. Sharp bounds on the distribution of treatment effects and their statistical inference. Econo- metric Theory, 26(3):931-9) 191111911112 2010.
Neyman, S. On the application of probability theory agricultural experiments. Essay on principles. Section 9, Statistical Science, pp. 465 472, 1923. 1990 reprint of the original 1923 paper.
Fogarty, C. B. Regression-assisted inference for the average treatment effect in paired experiments. Biometrika, 105 (4):994-1000, 2018.
Freedman, D. A. On regression adjustments to experimental data. Advances in Applied Mathematics, 40(2):180-192 2008.
Robins, J. M. Confidence intervals for causal parameters. Statistics in medicine, 7(7):773-7855 1988.
Samii, C. and Aronow, P. M. On equivalencies between design-based and regression-basee variance estimators for randomized experiments. Statistics & Probability Letters, 82(2):365-315 3700, 2012.
Hall, P. The bootstrap and Edgeworth expansion. Springer Science & Business Media, 2013.
Harshaw, C., Middleton, A., and Sävje, F. Optimized variance estimation under interference and complex perimental designs. arXiv preprint arXiv:21/2.01709, 2021.
Schochet, Z.
Is regression adjustment supported by neyman model for causal inferenee? Journal of Statistical Planning and Inference, 140(1):246-259, 2010.
Harshaw, C., Sävje, F., Spielman, A., and Zhang, Balancing covariates in randomized experiments with the gram-schmidt walk design. Journal of the American 2024.
Stoye, J. Partial identification spread parameters. Quanti- tative Economics, 1(2):323-357, 2010.
Statistical Association, pp. 1-24.
Tan, Z. A distributional approach for causal inference using propensity scores. Journal of the American Statistical Association, 101(476):161--163) 2006.
Horvitz, D. G. and Thompson, J. A generalization of sam- pling without replacement from a finite universe. Journal of the American statistical Association, 476260):663-685, 1952.
Wang. R., Wang, Q., Miao, W., and Zhou, X. Sharp bounds for variance of treatment effect estimators in the finite population in the presence of covariates. arXiv preprint arXiv:2011.09829, 2020.
Imai, Variance identification and efficiency analysis in randomized experiments under the matched- pair design. Statistics in medicine, 27(24)4485744873 2008.
Imbens, G. and Menzel, K. A causal bootstrap. Annals of Statistics, 4963):1400-1445 June 2021.
Imbens, G. W. and Rubin, D. B. Causal inference in statis- tics, social, and biomedical sciences. Cambridge Univer- sity Press, 2015.
IMF. World economic outlook database, 2024.
https: http::://wwwwnnm/ wwtp/////twiiaatttieeooooo http://oooaataaseee22222 Aprril
Ji, W., Lei, L., and Spector, A. Model- Model-agnostic covariate- assisted inference on partially identified causal effects. arXiv preprint arXiv:2310.08115 2023.
Lin, W. Agnostic notes on regression adjustments to ex- perimental data: Reexamining freedman's critique. THE ANNALS of APPLIED STATISTICS, pp. 295-318, 2013.
Lu, J. On randomization based and regression-basee infer- ences for 2k factorial designs. Statistics & Probability Letters, 1 112:72-78, 2016.
10

---

Integer Programming for Generalized Causal Bootstrap Designs
A. Appendix
A.1. Complete Simulation Results
Integer programs to compute the optimal causal bootstrap imputation were solved using the CP-SAT solver (https: /eveloppess google com/ com/optinnn zat wttptattt///a/aaaaaatattattttttt cp/cp cp_ sop_soover)) Experiments were parallelized on a cloud CPU cluster using 100 workers. Runtime for the final experiments was approximately hours. Smaller runs of approximately 10 minutes used debug the simulation code. Tables 4 and 3 report the complete simulation results CI width and power for multiplicative and additive effects, respectively. Table in the main body contains complete results for A/A tests.
Covariates
Effect Size
Inference Method
Median CT
Power
Compl. R.
Matched Pairs
Compl. R.
Matched Pairs
None
82
Sampling Boot.
3967
4110
0.10
0.00
None
82
Conservative Var.
3991
1152
0,02
0.00
None
82
Isotone Copula
3348
0.13
None
82
Opt. Causal Boot.
3348
2230
0.13
0.00
None
164
Sampling Boot.
3967
4110
0.10
0.00
None
164
Conservative Var.
3991
1152
0.03
0.00
None
164
Isotone Copula
3348
0.14
None
164
Opt. Causal Boot.
3348
2230
0.14
0.00
2018 GDP
Sampling Boot.
141
142
0.61
0.78
2018 GDP
82
Conservative Var.
143
49
0.61
1.00
2018 GDP
82
Isotone Copula
133
0.68
2018 GDP
82
Opt. Causal Boot.
133
102
0.68
1.00
2018 GDP
164
Sampling Boot.
141
142
1.00
1.00
2018 GDP
Conservative Var.
143
49
1.00
1.00
2018 GDP
164
Isotone Copula
133
1.00
2018 GDP
164
Opt. Causal Boot.
133
101
1.00
1.00
Table 3: Complete results for 95% CI widths and power for injected additive effects.
As footnote, we observe that the "true CT" widths increase with effect size for both designs when there are no covariates, as expected. However, they sometimes decrease with effect size with covariates. This is due to the fact that the covariates model some growth in outcomes, and these models are not fit in-sample. Therefore, it possible for variance reduction to be largest for nonzero effect sizes.

---

Integer Programming for Generalized Causal Bootstrap Designs
Covariates
Effect Size
Inference Method
Median
Power
Compl. R.
Matched Pairs
Compl. R.
Matched Pairs
None
2.5%
Sampling Boot.
4015
4160
0.10
0.00
None
2.5%
Conservative Var.
4040
1215
0.02
0.00
None
2.5%
Isotone Copula
3395
0.13
None
2.5%
Opt. Causal Boot.
3395
2256
0.13
0.00
None
5.0%
Sampling Boot.
4068
4206
0.10
0.00
None
5.0%
Conservative Var.
4079
1298
0.01
0.00
None
5.0%
Isotone Copula
3447
0.13
None
5.0%
Opt. Causal Boot.
3447
2284
0.13
0.00
None
7.5%
Sampling Boot.
4121
4260
0.10
0.00
None
7.5%
Conservative Var.
4114
1381
0.02
0.00
None
7.5%
Isotone Copula
3495
0.13
None
7.5%
Opt. Causal Boot.
3495
2314
0.13
None
10.0%
Sampling Boot.
4181
4306
0.10
0.00
None
10.0%
Conservative Var.
4149
1464
0.02
0.00
None
10.0%
Isotone Copula
3540
0.13
None
10.0%
Opt. Causal Boot.
3540
2341
0.13
2018 GDP
2.5%
Sampling Boot.
111
109
0.32
0.20
2018 GDP
2.5%
Conservative Var.
117
84
0.24
0.46
2018 GDP
2.5%
Isotone Copula
99
0.39
2018 GDP
2.5%
Opt. Causal Boot.
99
97
0.39
0.25
2018 GDP
5.0%
Sampling Boot.
122
0.93
0.98
2018 GDP
5.0%
Conservative Var.
124
148
0.89
0.75
2018 GDP
5.0%
Isotone Copula
102
0.97
2018 GDP
5.0%
Opt. Causal Boot.
102
126
0.97
0.92
2018 GDP
7.5%
Sampling Boot.
154
169
1.00
1.00
2018 GDP
7.5%
Conservative Var.
152
2018 GDP
7.5%
Isotone Copula
126
1.00
2018 GDP
7.5%
Opt. Causal Boot.
126
162
1.00
1.00
2018 GDP
10.0%
Sampling Boot.
185
235
1,00
1.00
2018 GDP
10.0%
Conservative Var.
189
310
1.00
0.93
2018 GDP
10.0%
Isotone Copula
154
1.00
2018 GDP
10.0%
Opt. Causal Boot.
154
203
1.00
1.00
Table Complete results for 95% CI widths and power for injected multiplicative effects.
12

---

Integer Programming for Generalized Causal Bootstrap Designs
A.2. Proof of Lemma 2.1 in Section 2.3
To show the first part of the Lemma 2.1, we observe that for any outcome yk setting implies Fabs (yx) pimis (yx). The observed distribution (yx) is already determined and not in our control. Further, Fy (vx) Fy (yx)
(1 - x which an integer value. Since Fi (yk) is already given, it may be that there exists an outcome y/: such that No P1 sk
(ye) not an integer, in which case the (d) constraints above are not feasible when O.
We now tackle the second statement of Lemma 2.1. The constraints corresponding corresponding to 0. Here, we focus
1 are decoupled from those and the case of a 0 can be treated similarly. Define the shorthand
(1) 2 N1
Pk 1ipbe (yk)
Σ
for every yb E supp(F("e). Note that the values of pa are given by constraint (a) and Σ. = 1. Let m = supp 170b0 /uppp and C be the controlled units (C)= =N)). For each [m]) we let Ag C C be an arbitrary set with |Axl [-Nop&] such that the sets A are disjoint. Note that
m NoPa] I No
Ael
A& = L
Σ Pk: = No- i=1
k=
k=1
Therefore, the number of controlled units not covered by any of A is at most
= 10=10u U Aul = No Σ Nopx = Σ
(Nops [Nopw]]
k=1
We next increase the size of |Ail..... |A, by adding one element of C²\ U/=11 A& to each of them (Note that this is possible since showed above m). This way, have
[& [Nopk] + 1, for k: ,,...,r, A& [Nopk], for k = I...., m,
the sets Ag are disjoint and their union covers C.
We are now ready to construct a feasible solution as follows. For i E C' (i.e., a controlled unit), we set
J1 i E A
x=( tk
otherwise.
In addition, for i E C (i.e., treated unit), set
=1 iff Y(1) = yk
(5)
Constraints (a) and (c) are clearly satisfied as x(1) are binary variables that are zero outside supp(Fofi)). Condition (a) is also satisfied by (5) and noting that for E C we have yooh = Y((1). To verify constraint (c), note that if i C, only one of ( is one and the rest are zero because the sets A are disjoint. The same holds for i E C clearly by the construction in (5).
We next verify constraint (d). We have
N
N
2,
x(1) ik
|
N
ik
Nu
=
Pk
No
=
Pk
No
1=1
1=1
IEC
For 1,...,7,, T, we have
|Ag]
[oop
+1
1
(Nopk
[Nop
Pk
No
Pk
No
No
<
No
13

---

Integer Programming for Generalized Causal Bootstrap Designs
Likewise, for k = 1, .......... m have
[& No
Nopk
Nops
NNpk
Ph
Pk
No
No
No
Combining the previous two inequalities, we obtain that constraint (d) satisfied for c c 1/No for a =
Following the same procedure for a = O, we get that there exists feasible solution to optimization (4) for 2
1/ //n
When No Ni, note that in this case we can form one-to-one matching of units in the control set to units in the treatment
set. Take any such matching, and for E let k; be the outcome of the unit that is matched to. Set X;; (1--2) 1 and
x(1-2)) = for all & + kū, Set x if y yı and x!2.) = otherwise. By construction, this solution satisfies constraints (4).a and also constraints (4).d for e 0.
A.3. Constructing the Integer Linear Program in Section 2.4
Let y(a) E RDD be the vector of potential outcomes V. Consider the following matrix and vector variables, with ( an N (N- K) matrix and (()) an (N K) vector:
(a)T
0
0
(a)
0
(
0
0
0
0
y(a)TT
( (x(aa
(a) ik
x((a
X
T
i(k+1)
NK
We further consider the following combined matrices:
Y: (:=(<t) No 1 No x
X X(1)
X
While the notation is clumsy, it allows us rewrite Varz[f] [f] in solver-friendly way, with X vector and objective the quadratic function. Recall that the variance of the difference- innecee-i-maan estimator has quadratic closed-form for general assignments: Varz[f] (Y,(1), Y,()) + Y ΣzzY Y, where Y is the vector of coordinates (N,1 + No Na'Y, (Y(O))== and Σz is the covariance matrix the random vector {Z;)i=1... with coordinates (Σzz)ig = Cov [i, Zj]. Rewriting Varz(f) in terms of X and Y, we obtain:
Varz X X -
Because Σzz is positive definite for any probabilistic assignment, it can be decomposed into the product of a matrix and its transpose RTR by the Cholesky decomposition. Thus, Q can be written as the positive semi-definite matrix Q Q= YYRTRRR..
Furthermore, in the objective (4) we have products binary variables of the form X x(a)
We replace this product with ikkje E {0, 1}, and introduce the following linear constraints to ensure the variable corresponds to
X je
a new binary variable x-(..) product of the original variables:
((ab) ik.j s (() ik
(a.b)
ik.je < (()
x
(a.b) 1k.je >
x ik
+ (b) -
is easy to check that x(aab) = if and only both X = 1 and x je consists only of binary variables, and the objective and constraints all linear.
1. After this transformation. the integer program
14

---

Integer Programming for Generalized Causal Bootstrap Designs
A.4. Derivations of Section 3.1
For any assignment mechanism where No. Mi >
the difference-in-meaan estimator can be written
(1) Y((0)
Y0(
Σ Z
Σ Z,a; + b
N +
No
No
where N1 is the number of treated units and No is the number of controlled units. Estimators of this form have expectation over all assignments:
E(f) Σ, ElZi]a; + b; and we can write their variance
N
N
Varz -arz P(Z= PP(2-= =2
Σ
n,, + b, -E(f) Σ P(Z =))
Σ a [[((s-E(Z.))
Substituting for Dj = 2 El Zz], we have Varz[f] = Varp(f) = Px = d) ( a,d,)² = Ea =
d) (Ei, dıdja;a)) = Ei, 0,a, E, such that the objective function we are trying to maximize with respect to (1) and Yi(0)
Varz(f]= Σ a,a Cov Cov[?,,24 (4,
(6)
A.5. Proof of Theorems in Section 4.2
In optimization (4) the slackness parameter was chosen allow for feasible solution. Other than feasibility, note that while Ex([""" section, we derive a finite sample analysis of optimization (4). Specifically, the next theorem provides lower bound on the probability that for given sample, the two distributions Fmis and Fobe are uniformly close to each other, namely within a distance of at most e. This also gives a rate for the probability of the optimal objective value of (4) upper bounding the quantity of interest Varz(f) as the sample size increases.
Bz (v)], for given sample the distributions Pmis and Pobs may not match exactly. In this
For the readers' convenience recall the statement of Theorem 4.1. Theorem A.1. Suppose that (i(a) ~ F. for E {0,1} and [n]. Also assume that {(Y,(0), Yi(1)) HE(N] and the assignment variables [ZiJiEIN] are independent. Let V be the optimal objective value of the problem (4), We have
P(V* 2 Varz(f))2 21-B
(7)
where
8 : 8exp
2 -NNP) +
32
Σ
Cov Zz).
4
N2
[JE[N
p= min (P,1 P),
P
1
Σ
1).
N
iE[N]
Proof. We start by finding a lower bound for the probability that the following statement holds:
N
Z
Va E (0, 10,1), Ak,
Σ
x())
(16 2
N
ik
Since is independent of Y(a), , the set of (r : [N], 2, = 1} is a set of random draws from F We condition on (DKW) inequality', we have that for every > 0,
and let N₁ = EN Ze By using the Dvoretzky Evoretzky-Kiefer-
P
sup Σ
2x (a) N1 F(()
( (
2c 222
'Note that DWK inequality provides a bound on the worst case distance of empirically determined "distribution function" Ō from its associated population "distribution function" G. For discrete variable, can write the pmf value as Fly) = G(y+ 8) G(y 5), for a small enough and apply DWK bound on G(y + 6) and separately to get similar bound F(y).
15

---

Integer Programming for Generalized Causal Bootstrap Designs
By taking expectation with respect to Z we arrive at
P
sup
N
N1
Z
x
ik
F((U))
>2 2
6
- 2E|e 121
(8)
Olk
The correlation among Z will affect the rate of convergence. By using Chebyshev's incquality, we have
P
Σ VE[M]
2
Σ iE[N]
P
P
No)
Σ
Cov (z,77)
(9)
N282
[JE[N]
with PP = 1). Therefore.
P(N
Σ P. N6)
S
1
Σ
Cov(ZiZg).
(10)
N282
iE[N]
(JE[N]
We set 6 1 2N EiE(N] P., we have
P (
S
Σ
P. EP)
S
4
Σ
Cz,
2 (E[N)
C P.Y
[JE[
Let & be the following probabilistic event:
8
8=
f Ni
1 Σ
P...
KE(N)
Using the above bound bound on IP(ES), we get
2
Ele 1
S exp
(
Σ
P.
P(S) +P(E)
4 iE[N]
S exp
Σ
P.
+
Σ
2))
(11)
4
E
P.Y
iE[M]
1.J[[M]
Combining (11) and (8) we obtain the following result:
IP
sup gk
Z
sk
(a)
> 5) S 2exp ( (- 4
2
Σ
P.)
+
C, P.J
8
Σ Cov 2i)
[E[N]
AJENN]
= 2exp
1--N- NP)
+
8
Σ 2g)
(12)
N2P2
JJE[N]
By applying a similar argument (using
Z; instead of Z₂) we also get
N
IP(
sup
Zi F2(2a)
x(a)
>
-) 2
< 2exp
(-ENI
N(1 (o(1-)))=
+
8
Σ iJE(N]
Cov 2z). (13)
No
4
N2 pLZ
In addition, by triangle inequality we have
Zx x(aa==
1 - 7. x(aa)
7x ((a)
Za
(v))
Σ
x
No
N1
No
Therefore, by union bounding over the two events in (12) and (13) and over a E {0, 1} we get that
N
1- Z
Ak, a,
se,
N
No
16

---

Integer Programming for Generalized Causal Bootstrap Designs
with probability at least where
8 = 8exp
2 NP) (-
- 4
+
32
Σ
Cov(Zi, 2z),
N2 p2
(JEIM
with P = min(P,1- P)-
Hence, is feasible solution to (8) with the same probability. addition, recalling (3) we have
=" -X-) 2 X'QX= Varz(f),
which completes the proof.
Remark A.2. Under the setting of Theorem A.1, if the assignment variables Z; are independent, then we can sharpen the probability bound using Hoeffding inequality in (9). In this case, the failure probability /3 will be given by
2² 1-=NN NP) 8exp N-2NPt) (-
8 = 8exp
Confounded independent treatment assignments. We next prove a similar result for the case that the treatment variable can be correlated with the outcomes through a confounding factor.
For the readers' convenience we recall the statement of Theorem 4.3.
Theorem A.3. Suppose that (Y)(0), Yi(1), Z4) are i..d across E [n] and (;(a) ~ F for a E {0, In addition, suppose that
PPZ ((Y())))(1))) a)
88,
for some 5 > O. Note that â controls the effect of confoundedness, with ō = 0 corresponding to unconfoundedness. Let V be the optimal objective value of the optimization problem (4), and suppose that the slackness parameter e, in constraint (d) satisfies e 2 5. We have
P(V 2 2arz(f)) 11
(14)
where
8 : 1 8c -2Nc
86 -NPe7/3
(c (c-8))((1+
c-8
p min( (1,1 P), = 5+5+c
61 2+
Proof. Recall that ik = I(Y = yl.). We first establish the following lemma.
Lemma A.4. Under the assumption of Theorem A.3 we have
2, x (8) P(Z, = a)
58.
We define the probabilistic event &, follows, for fixed arbitrary en 0:
Z P 1(z==1ix()) F (yk) S sco
E,
sup Uh
N
Note that by assumption X are i.i.d and we denote their distribution by dropping the index as x
By an application of Dyorezky-Kiefee-Woloooe (DKW) inequality, similar to the proof of Theorem A. have
P(E)) & 22N
17

---

Integer Programming for Generalized Causal Bootstrap Designs
We also define the probabilistic event & as follows:
N
-1) 1) Se P(Z P(2 -112-111 1
غ =
Σ
P(Z - 2
2:
(15)
N
1=1
Recall the notation N1 :- 2i 1 Zi and let P = P(Z 1). Then, by an application of the multiplicative Chernoff bound we have P(E) I 2e -NP//33
Therefore, on the event & n E, we have
Z
Ei=1
7,X((a)
(0)
ik
N1
ik
Σ i=1 Zz
N
24 X(a)
P(Z = 1|X ( Pa (1x) PP(Z 1|x((0)) Fa (Uk)
i=
ik
P(Z= - PP = 1) $/2 E/ 2z
60 + P(Z 1
P(1 (1)
$ 60 1 5
Fz (yn+),
((1-g)
11
where in the last step we used the result of Lemma A.4, by which P(Z = 1|xx < 1 Likewise, the event & E,E, we have
7
PíZ 1|x; (() ) (fs)) -0
X
(a)
N
ik
P(1 +61
2
1 1 8
(v:)) p(1+G))
60
+ 61
Therefore, on the event & n Nε we have
Zi
x(a)
(2k)
sup 9k
N1
ik
1=1
< max
6D
+
1+6 -1)
FT (mn):
E0 P(1 +c) + (1-1 1-0)
F (2x)
P(1 61)
1 --er E1
60
(1 0
P(1 e1)
1-61
since F. (3k) S 1. In addition, by union bounding.
P(Б n Ł) 2 - P(89) - -P(89) 2 1 2e -2Ne 2e /N
By following a similar argument (replacing the role of Z, by 1 - Zz), obtain that
N
sup 5/k
Σ
x(a)
F. (2k)
1 (1 - P) - -1)
60
+
61 5
Nu
1-
i=1
with probability at least - 2c Therefore, by triangle inequality and union bounding over a E (0, 1} we get that
2N
2c -N(1
Yk, O,
Σ
(a)
Z
1 1
Zz
5 2
60
+
6 + 8
(16)
N₂
No
P(1 - (1)
1-6 61
18

---

Integer Programming for Generalized Causal Bootstrap Designs
with probability at least where
B = 1 Se -2Nc 88 /NP
and P min(P, P).
We next ED, ei such that the right-hand side of (16) becomes e, Namely, we set
(e (e-6)
e
60 2+0+6
61 2+6+6
This way, X is a feasible solution to (8) with probability at least - 3, and so
V* = (x*) (x*) 2 X'QX= Varz(f].
Proof. (Proof of Lemma A.4) By definition, P(Z al x( = 1) P(Z. aa(Y(aa) = =yk). Therefore, from the
assumption we have
P = a|x
ik 1) 1 58.
P(Z: a)
For the other case we write
P(Z, = 1|x (1) = =0) IP(2: 1, (())
= 0)
Σesr P(Z: = 1|Y7(0) = ())P(( = ve) P(X (0) 0)
PP (a) =0)
Expanding the denominator PP(((1) (a) 0) Eijr ((V;(a)) yt). get
P(Z: = 1|x (a) =0)
Eipr PIZ: (| w.))P(Y((a)) = We)
P(Z a.
Eesr P(Z = (a)) =YE)
Define the shorthand dę P(Z, 1/Y y (Y) = =ye) and be P(Z, aY)P(Y.(a) ye). Then by assumption by - 1| S 8, or equivalently, (1 5)b₂ S $ (1 11+5) 8)b.. This implies that
1-85 PP 1|x ik =0) Eeak az P(Z, = a) Σejr be
(a)
51+5.
This concludes the proof lemma.
A.6. Derivations of Section 3.3
We first prove that the variance of estimators of the form i = E;b: + E3 ZiZzaij can be written as a sum of 4 products of asj and age terms:
Varz
Σ P(Z 2)
Σn: bu Σ ZiZzaij Σb. b: Σ B B[Z,Zzlaj
Σ P(Z =)
Σ (Z:Zz - B(2:2z)) aij
i,j
Σ k.
aiyaki Cov [Z:Zy; Zi
We now prove that the regression estimator on Z and covariate X (later rewritten to W to avoid confusion with the indicator variables) can be parameterized such a way.
19

---

Integer Programming for Generalized Causal Bootstrap Designs
Suppose that I' is linear regression with covariate X. The exact model is Y; = + B2X,, where D, Z, - x |1 gives the average treatment effect. Without loss of generality, assume x E; X; 0 (any nonzero can be absorbed into the intercept). We can analytically compute the hat matrix, which is simplified by the fact that E, D; and >. X; 0.
N O
O E, X,D, D,D,
X X =
No
0
Σ
(x())
We compute the inverse of this matrix using Gaussian elimination. Because we only need the second row of the inverse, it suffices to only partially solve the Gaussian elimination. We only compute the first two rows of the inverse:
0
(XTx)-1=
N
N
E,
=
NYNo ?
NINo
(x+)
The average treatment effect estimate 81 is given by
E Yabe yob Dj Yoba Xs
-= xxxxxxx-
N
N
E. X.D.
No
No
[x(3?)
N
Σ Y.oho Di)
N
E, XiD: E,(x?)
Σ y-obo Xi
No
N1
N
E3 X
Σ
Dj
yobbe -X,
No
Σ, x3
Next, we perform two algebraic manipulations write our estimator in terms of the potential outcomes (instead of Y,ob)).
Yobs = Z1)((1 (1 Z,)Y.(0) =: ZY
N1
D, D,Yoha
2i
(Z.Y(11)+ + (2))((0)
N
N1
N1 Yi(0) + ZY
N1
ZY((1)
ZY((1)
N
N N YY)
N
(Xx(1)
No
Ni
Zi (1(1))
Y,(0)
N
N
N
No No
Yi(1)
Y((0)
N
No
N V2(0)
Z:
Y2(0) =: 7,
V.
N
N1
No
N
N
N
where we introduce the notation Y, = Y,(1) (, and 8, Y1())
Y.0))
+
No
20

---

Integer Programming for Generalized Causal Bootstrap Designs
N No
Σ
DiY, x
N
Zi N Y,
Ni N xtx (2,XX2228,,
No
N No
E(
N N XY((
N +
Σ
Z, N Y,
x Σ (Z: +
1
NINo
No Σx((0)
N N.N.
X ΣΣ
X - X.X, (Z:Zzi + - * *r((1))
1 Σ Y,(0)
No x ΣΣ X,X)Y,(0)
1
No
N
x N Y, - -XX,, (Z.2,i + 2Y, N4 2113)
NoNo
X
follows that
f Σ
No (v((0)x,
ΣΣ
X X,X, -
N.
X
Because the double sum contains all combinations of and j. we can reindex i and j in the term N: Z,gg
by NINo
ΣΣ X
x
N N - X,X, 2,258 + ZY, No x..)
Σ b No X ΣΣ
Z2 X NoN XxX,Y)(0) m (9,)-22,,x. -22,xxxx,))
)
We see we can write f = =. + ΣiE, aij Z.j, where
1. (Y:(0) + + 1T
No
1 f N
No X X, X,X,Y,(0)) + No x fi=
aij
No X X.X,V,
N
otherwise
Substituting in our definitions of Y; and ýi we have
6,
1
1T (Y:(0)+ + Xi X
No
1 Y NoW Y2(0) - N₁ No xix X,Xx,,(0)) + 1 x (Y1(1) - Yi(0)) ifi=j otherwise
1
N
No
a.j
NoNo X XxX, (Y,(1) - (.))
A.7. Conservative Variance Bounds for Section 5
Here we make explicit the analytical, conservative variance estimators used as baselines for the complete randomization and matched pairs designs in Section Recall that No, Nj are the number of units the control and treatment sets, respectively,
21

---

Integer Programming for Generalized Causal Bootstrap Designs
that N = No + N1 is the total number of units, and that Z, {0. 1} is the treatment indicator for unit For {0, 1),
define the sample variance estimator
2
Yi(2)
N. Yi(2)
Σ
N
An estimator for the sampling variance of diff-in-meant under complete randomization is as follows; see Theorem 6.3 in (Imhens & Rubin, 2015).
No
N1
Comparing to (I), see that this estimator is conservative upper bound for the true sampling variance S 0.
For the matched pairs design, let F(j) be the difference between the outcomes for the treatment and control units in pair j. for j E [2] A conservative estimator for the sampling variance of diff-in-means under the matched pairs design is as follows; see Theorem 10.1 in (Imbens& & Rubin, 2015).
N/2
N/2 +(J)
Σ +(()) j=1
N(N - 2) J=1
N/2
A.8. Variance of the Difference-in-menne means Estimator under a Completely Randomized Assignment The variance of the difference-in-menn iennnss estimator given by
S?
S
S?
VarzlfDIMI = N1
(17)
No
N
(i(()) - N Ew v1(x)) 2
where 52 1
N-1 Emm
is the sample variances of the treated (2 = 1) and control (a = 0) units,
while S7 ;= N-1 Ei (Y/(1) -((0)--)0)-))) 1. EE (Y,(1) - -Y(O)) is the sample variance of the unit-level treatment
N
effects (1) Y:))).
A.9. Constructing confidence intervals
We propose two different ways of constructing confidence intervals from the optimal value V the quadratic program above. The first approach is to use a normal asymptotic approximation. In this case, the (1 - a) -c-ofdfidence interval is computed + (VV~/N) 21-0/2 with 21-0/2 being the (1 a/2))-qunntin of the standard normal distribution. This approach is justified in some settings. In particular, Abadie, Athey, Imbens, and Wooldridge (2020) study the asymptotic validity of this approach for linear regression. In other settings, it can be viewed as an approximation, the properties of which can be investigated by running simulations or A/A-tests. Another approach is to employ a bootstrap procedure similar to the one proposed in Imbens Menzel (2021), in which unobserved potential outcomes are imputed via the least favorable copula and bootstrap samples are drawn via rerandomization of units to treatment. Confidence intervals are constructed based on the quantiles of the bootstrap distribution f, as illustrated Algorithm 1. Many variants of the bootstrap can also be applied to a causal bootstrap; for example Imbens and Menzel (2021) bootstrap the t-statistic which is asymptotically pivotal and is known to lead to asymptotic refinement of the confidence interval in standard settings (Hall, 2013).
In Algorithm 1, (F.) is the support of the marginal distributions of outcomes Fa. If supp(Fa) is not known (e.g. continuous outcomes), we can use the support of the observed marginal distributions supp Faahs) in Algorithm 1. We revisit the validity of this assumption in 4.2, and the specific parameterization of the solution X in section 2.4.
A. Understanding the undercoverage of bootstrap-based methods
Copula-based methods are liable to suffer from undercoverage in small samples. This is true the isotone copula method proposed by Aronow et al. (2014) and used by Imbens & Menzel (2021), as well as our own method. These concerns diminish in larger samples. Consider the following scenario: The true effect Regardless of the realized randomization,
22

---

Integer Programming for Generalized Causal Bootstrap Designs
Algorithm I Our Causal Bootstrap Procedure
Require: for a E { y(a), vector representation of supp(Fa); X()), solution to Eq. 2; design P(Z). Initialize: 5 0
for E {0.1)
(a) t (x()))
end for
for some large number of repetitions do
Sample Z ~ (Z) and S + u (f(Z,
end for
Return [ya/2, aa/2 (1x is the th quantile of 5).
Unit
Pair
Yo
Y1
A
11
11
B
10
10
0
any copula that matches the marginals will lead non-zero ATE deterministiiaalyy For example, in each of the observed
potential outcome tables below,
Unit
Pair
Yo
Y1
Unit
Pair
Yo
Yi
A
11
A
11
B
10
X
B
X
10
the copula will always be imputed either:
Unit
Pair
Yo
Y1
Unit
Pair
Yo
Y1
A
10
11
A
11
10
B
10
11
B
11 10
In either case, the samples of the estimated ATE will be 0.5 -0.5 deterministically, leading coverage zero. These issues disappear in larger samples. We refer the reader asymptotic validity results for theoretical guarantees in this setting.
23


Reviews:
[{"review_id": "B2489CEhFn", "reviewer": "Reviewer_bwtT", "summary": "The paper presents a method that employs integer programming to maximize the variance of proposed estimators in randomized experimental design, addressing the issue of design uncertainty. It extends linear-in-treatment and quadratic-in-treatment estimators and generalizes assignment mechanisms using integer programming. The approach is built on a strong theoretical foundation.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and positive review of our paper. Below, we address the two questions raised. \n\nQ1:  Unconfoundedness implies that the treatment assignment is independent of any potential outcome, usually achieved through randomization. We distinguish between conditional and unconditional unconfoundedness, i.e. whether or not this independence requires conditioning on an observed covariate, though the literature does not always make this distinction. Strictly speaking, a stratified assignment is not unconditionally unconfounded and is instead conditionally unconfounded if the stratifying variable is correlated with the potential outcomes. In observational studies, we usually cannot eliminate confounders, and it is more common then to assume conditional unconfoundedness. An excellent discussion on the relationship between randomization and unconfoundedness is [1].\n\nBounded confoundedness is not a standard defined assumption like unconfoundedness is. It typically arises in the context of sensitivity analysis and assumes that the magnitude of any confoundedness is limited within some range, i.e. one where treatment probabilities do not change drastically with remaining confounders. An individualistic assignment is one where the treatment assignment of one unit is independent from the treatment assignment of other units. For example, a Bernoulli randomized design is individualistic, but a completely randomized design, strictly speaking design is not (it is unconfounded however, so it is covered by Theorem 4.1). \n\nWe would be happy to expand on these assumptions and include them in the paper if the reviewer would find it helpful.\n\n[1] Sävje, Fredrik. \"Randomization does not imply unconfoundedness.\" arXiv preprint arXiv:2107.14197 (2021).\n\n\nQ2: Regarding IP scalability, we would first like to re-emphasize that our method is primarily motivated for small-sample experiments where design uncertainty dominates. Also, we note that given the results of a particular experiment, the IP only needs to be solved once. (In our simulations, the IP had to be solved multiple times for coverage and power analysis.) Nonetheless, we agree that for a general-purpose tool, scalability is important.\n\nThe scalability of solving the IP depends on the sparsity of the treatment covariance matrix. This is why we see sub-quadratic scaling for Matched-Pairs in Table 1 (Appendix A.8). For the case of complete randomization where the off-diagonal entries to the covariance matrix are all equal, an LP formulation is possible as mentioned in Section 5. \n\nTo scale the IP beyond hundreds of units, we see potential to apply relaxation or approximation algorithms, as also suggested by Reviewer XXTL.\n\nFor a relaxation, we could relax the IP to an LP. The objective would then provide an upper-bound (i.e., an overestimate) to the worst-case variance, which can be used for Neyman-style confidence intervals (e.g., 1.96 * sqrt(variance upper bound)). However, the LP solution could not be used for the causal bootstrap because it might not be feasible for the IP.\n\nFor approximation, on the other hand, we could apply some technique like LP rounding to obtain a feasible solution, but this would now underestimate the variance. However, if we can guarantee that we have a c-approximation (so that the solution objective is 1/c of the optimal, for c >= 1), then one could scale imputed unit outcomes by sqrt(c) to ensure the variance of the causal bootstrap distribution upper bounds the true variance. Developing LP rounding techniques for our IP is still an open question.", "claims_and_evidence": "Yes, the paper provides theoretical guarantees and empirical evidence to support the claims.", "methods_and_evaluation_criteria": "Yes, the authors validated the proposed approach via simulations on real data and compared it with three baselines.", "theoretical_claims": "Yes", "experimental_designs_or_analyses": "Yes", "relation_to_broader_scientific_literature": "Building on previous work, this paper extends to linear-in-treatment and quadratic-in-treatment estimators and introduces new assignment mechanisms using integer programming.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "**Strengths:**\n1. The paper tackles the important challenge of design uncertainty in experimental design and extends previous work by incorporating linear-in-treatment and quadratic-in-treatment estimators, along with new assignment mechanisms using integer programming.\n2. The proposed approach is supported by asymptotic guarantees under different assignment mechanisms.\n**Weaknesses:**\n1. Computational complexity and scalability of the integer programming formulation are not thoroughly discussed, which may pose challenges for large-scale experimental designs.", "other_comments_or_suggestions": "No", "questions_for_authors": "1. The authors establish the asymptotic validity of their approach under unconfounded, conditionally unconfounded, and individualistic assignments with bounded confoundedness. Could the authors elaborate on the assumptions underlying each of these assignment mechanisms?\n2. Solving Integer Programming problems is computationally expensive as the problem size increases. Could the authors provide a more detailed discussion on the computational complexity and scalability of their approach?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1C2cI4o7yG", "reviewer": "Reviewer_XXTL", "summary": "This paper proposes a novel method for quantifying **design uncertainty** in causal inference settings, particularly when experiments involve small samples, heterogeneous treatment effects, or non-standard assignment mechanisms. The standard bootstrap only captures sampling uncertainty, while existing causal bootstrap methods are limited to completely randomized designs and simple estimators like difference-in-means. The authors generalize causal bootstrap by formulating the worst-case copula problem as an **integer programming** (IP) task, applicable to a broad class of estimators (linear- and quadratic-in-treatment) and assignment mechanisms (including unconfounded, conditionally unconfounded, and bounded confoundedness cases). They prove asymptotic validity of the approach, and showcase improved variance estimates and tighter confidence intervals on simulated geographical experiments using IMF data. The approach enables exact variance maximization without closed-form copulas, providing flexibility and rigor for real-world experimental designs.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful engagement with our work.\n\n## Responses to Questions\n\n1. Thank you for this question; please see the discussion of scalability in our response to Reviewer bwtT.\n\n2. We currently discretize the continuous outcomes using the grid of observed outcomes, so that the discretization is lossless. A bigger issue related to discretization is that we require the marginal distributions of the imputed joint distribution to match the observed marginals. While the observed marginal will converge to its expectation asymptotically, in finite samples there can be a significant gap between the empirical sample and the marginal distribution (characterized by the DKW inequality). This finite-sample gap can be mitigated through the epsilon slack parameter; see Section A.5 for results.\n\n3. This is a nice suggestion to improve efficiency of the IP. We have not considered it, but it would be worthwhile to investigate in future work.\n\n4. If SUTVA holds, our approach covers clustered designs, since they are usually implemented to be probabilistic and unconfounded, and the estimators they are usually paired with are linear- or quadratic-in-treatment. It is not immediately clear whether our approach yields good coverage when spillover effects are present. We agree with the reviewer that proposing a bootstrap procedure in this setting is an interesting area for future work.\n\n5. This is a great question, and we would be happy to provide more practical guidance. In the regime of feasibility, the smaller epsilon is, the tighter the upper-bound for the variance, but the stronger the coverage guarantees as provided by Theorems 4.1 and 4.3 and Corollary 4.2. Some practitioners may choose to contrast the upper-bound provided by the Integer Program against these theoretical guarantees to select epsilon. In our simulations, we solved our IP with epsilon = 0 (feasible since N_1 = N_0, cf. Lemma 2.1) and found the coverage to be acceptable. The adaptive or selection of epsilon is an interesting area for future work.\n\n6. We did not fully understand this question. For the matched-pairs design, the optimal copula is the one computed by our IP. If the reviewer was asking whether there are other (non-IP) ways to construct a copula for this design in the literature, then we are not aware of any.\n\n7. The causal bootstrap proceeds by first imputing a single (deterministic) outcome to each unit for the treatment it did not receive, which is why we need binary assignment of outcome to units (hence an IP). A convex relaxation would overestimate the variance of the worst-case copula, which could negate any power gains from the causal bootstrap. The same applies to dual OT (which would have the same objective value as the convex relaxation by strong duality). In the case of complete randomization, we found that an LP reformulation of the imputation problem was possible, but this was not the case for general treatment covariance matrices.\n\n8. Our method is also a copula-based bound; we are optimizing the least favorable copula numerically instead of deriving it analytically. Our method differs from previous work because we allow a broader class of treatment assignment mechanisms. For some concrete comparisons: (1) Aronow et al (2014) derive the assortative copula as variance-maximizing for the difference-in-means estimator under complete randomization. Our technique recovers their copula in this specific setting. (2) Fan and Park (2010) estimate quantiles of the treatment effect distribution under effect heterogeneity, assuming iid treatment assignments. They show that the assortative copula is not sharp for measuring quantiles of the distribution of the treatment effects. By contrast, our method estimates the ATE and allows assignment mechanisms other than iid.\n\n9. This is an interesting idea. The choice of estimands and estimators grows with the number of treatment options, but a priori, such an extension should be possible in most cases, and would be worthwhile to explore in future work. Beyond the asymptotic validity results which would need careful consideration, one would need to make sure that the estimator variance and potential outcome constraints under multi-valued treatment remain optimizable.\n\n10. We are deeply committed to open-sourcing valuable code and research, and hope to do so here. The feedback from practitioners has been positive: many practitioners find bootstrap methods intuitive, and shy away from complex variance formulas, which was the impetus for this work. Anecdotally, we find many users of A/B testing platforms shy away from implementing sophisticated randomized designs because they fail to see strong variance improvements when paired with the standard (but often incorrect) bootstrap confidence interval construction methods. We hope that making the construction of correct confidence intervals easier for practitioners to implement will encourage them to adopt more sophisticated designs.", "claims_and_evidence": "The main claims are:\n- Existing causal bootstrap methods are limited in scope; integer programming enables generalization to arbitrary estimators and assignment designs.\n- The proposed method provides valid, tighter confidence intervals across a range of realistic designs and estimators.\n- The method remains asymptotically valid even under bounded confoundedness and conditional unconfoundedness.\n\nThese claims are strongly supported:\n- The integer program is precisely formulated with all necessary constraints and relaxation mechanisms.\n- The theoretical results (Lemmas and Theorems 4.1–4.3) provide probabilistic bounds and convergence rates.\n- Simulations across additive and multiplicative effects, covariate scenarios, and design types (CR, matched pairs) validate the empirical advantages.", "methods_and_evaluation_criteria": "The core methodology—variance maximization via an integer program under potential outcome and assignment constraints—is well-motivated and mathematically sound. The paper systematically builds from the classic Neyman decomposition to modern extensions, supporting:\n- Arbitrary potential outcome distributions via discretization\n- Arbitrary linear- and quadratic-in-treatment estimators\n- Known and probabilistic assignment mechanisms (CR, Bernoulli, matched-pairs)\n\nEvaluation is appropriate:\n- Experiments simulate realistic geographical designs\n- Estimators include difference-in-means and doubly robust estimators\n- Baselines include sampling bootstrap, conservative variance, and isotone copula", "theoretical_claims": "The paper contains multiple rigorous theoretical results:\n- Lemma 2.1 guarantees feasibility of the IP under relaxed marginal constraints\n- Theorem 4.1–4.3 prove asymptotic validity under unconfounded, conditionally unconfounded, and bounded confoundedness regimes\n- Proofs rely on marginal balancing, inverse-probability weighting, and probabilistic bounds (Hoeffding-type inequalities)\n\nThe clarity of the mapping from copula constraints to linear conditions on binary indicator variables is a notable strength.", "experimental_designs_or_analyses": "The experiments are well designed:\n- Real GDP data is used to simulate geographical treatment effects\n- Both CR and matched-pairs designs are compared\n- Outcome models include both additive and multiplicative effects\n- Covariate information is incorporated through doubly robust estimation and growth modeling\n\nCI width, power, and coverage are all reported. Tables are clear and informative. Solver runtimes and scalability are also discussed.", "relation_to_broader_scientific_literature": "This paper fits well into the literature on causal inference and design uncertainty:\n- Extends Aronow et al. (2014), Imbens & Menzel (2021) on causal bootstrap\n- Goes beyond analytical copula assumptions by solving for optimal coupling numerically\n- Relates to optimal transport and Frechet bounds literature\n- Empirically complements recent work on balancing designs (Harshaw et al., 2024)\n\nPositioning is clear and the novelty of integer-program-based causal bootstrap is well-motivated.", "essential_references_not_discussed": "Most relevant references are covered, including:\n- Neyman, Aronow, Imbens, Robins, Harshaw, Ji et al., and causal bootstrap literature\n\nOne possible addition: more discussion of connections to **causal bounds via OT or copula methods**, such as:\n- Ji et al. (2023, AISTATS) used dual OT but not explicit worst-case coupling\n- Literature on robust causal bounds under partial identification (e.g., Manski-style or Fan & Park (2010))\n\nHowever, this is a minor omission.", "other_strengths_and_weaknesses": "**Strengths:**\n- Elegant use of integer programming to generalize causal inference tools\n- Applicable to a wide range of estimators and designs\n- Excellent theoretical backing and empirical validation\n- Readable and well-structured exposition\n\n**Weaknesses:**\n- Integer programming scalability limits application to large-scale datasets (>1000 units)\n- Requires discretization of outcome support, which may lead to approximation error\n- Some approximations (e.g., empirical matching of marginals) may be fragile in heavy-tailed or sparse data settings", "other_comments_or_suggestions": "NA", "questions_for_authors": "1. **Scalability to Large-Scale Experiments:**  \n   Your method is elegant but computationally heavy. Do you foresee any way to scale to thousands or millions of units (e.g., through greedy relaxation, LP relaxations, stochastic approximations)?\n\n2. **Robustness to Outcome Discretization:**  \n   Discretizing continuous outcomes is a strong assumption. How sensitive are your variance bounds and coverage to discretization granularity?\n\n3. **Adaptive Grid Design:**  \n   Have you considered adaptive or non-uniform binning (e.g., quantile bins) to improve efficiency or accuracy of the integer program formulation?\n\n4. **Extension to Clustered Designs:**  \n   Can your approach be extended to clustered randomizations or interference settings, where design uncertainty is entangled with spillover effects?\n\n5. **Automated Constraint Tuning (ε):**  \n   Is there a principled way to select or adaptively shrink the marginal balance slackness ε, beyond theoretical feasibility guarantees?\n\n6. **Limitations of Matched-Pair Analysis:**  \n   Your method identifies matched-pair isotone copula bootstrap as degenerate. Are there better imputations in that regime (e.g., hierarchical Bayesian)?\n\n7. **Alternative Optimization Frameworks:**  \n   Why integer programming instead of convex relaxations or dual OT (e.g., Kantorovich dual)? Is IP strictly necessary to preserve optimality?\n\n8. **Comparison to Copula-Based Bounds (Fan & Park, 2010):**  \n   How does your method relate in tightness or assumptions to classic copula-based bounds, or nonparametric bound estimators in econometrics?\n\n9. **Multi-valued Treatments and Interaction Effects:**  \n   Could the integer program be extended to handle multi-valued treatments, interactions, or factorial designs beyond binary assignments?\n\n10. **Real-World Adoption and Open-Source Tools:**  \n   Are there plans to release user-friendly packages or APIs for this method? What feedback have you received from practitioners in A/B testing platforms or policy experiments?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZrQGpt0fTZ", "reviewer": "Reviewer_bxyv", "summary": "This paper proposed a new integer program to jointly address two sources of uncertainty in causal inference, the design uncertainty due to the treatment assignment mechanism, and sampling uncertainty. Traditional methods tend to address one of the two uncertainties, but do not handle them at the same time. Motivated by this gap, the paper proposed an integer program formulation which computes numerically the worst-case copula used as an input to the causal bootstrap method. Further, the paper proved the asymptotic validity of this method for unconfounded, conditionally unconfounded, and individualistic with bounded confoundedness assignments. Numerical experiments support the effectiveness of the proposed methodology.\n\n### update after rebuttal:\nI have read through the authors' response and updated my scores accordingly.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful consideration of our work and overall positive review. The main concern seems to be with the number of datasets and baselines. We were unfortunately limited by space considerations, and relegated additional experimental results to Appendix A.1. We are excited about evaluating this method in a wide variety of contexts against many possible baselines, though we believe a truly comprehensive evaluation is best left to future work given the current space constraints.", "claims_and_evidence": "Overall the paper stated its claims and evidence clearly. \n- Theoretically, the proposed integer program aims to identify the joint potential outcome distribution that maximizes the variance of the chosen estimator, while being consistent with the randomization design and the observed marginal distributions of potential outcomes. The proposed optimization's objective and constraints were described in detail in Section 2 for the basic difference-in-means estimator.\n- On the statistical validity front, the paper provided a rate for the probability that the proposed method upper-bounds the true variance in large sample limit, by bounding the probability that the observed and missing marginal distributions are uniformly within epsilon distance. \n- In the special case of conditionally unconfounded assignments, the paper shows that one can further imposes equality of the conditional marginal distributions to fully utilize the covariate information.\nThe main claims were supported by convincing evidence.", "methods_and_evaluation_criteria": "The proposed method matters the most in settings with small fixed samples and heterogeneous treatment effects, such as in geographical experiments. Thus the paper used the GDP data report by IMF as a real-world geographical dataset. The proposed procedure were compared with three baselines including standard bootstrap and causal bootstrap proposed in prior works. The dataset and evaluation criteria are reasonable.", "theoretical_claims": "Checked through section 2's theoretical Claims, including proofs in appendix A.2 - for which I did not find issue with their correctness.", "experimental_designs_or_analyses": "I did not find particular issues with the experimental design.", "relation_to_broader_scientific_literature": "The proposed integer programming is novel and interesting as in the broader literature, which brings a gap between the design uncertainty and sampling uncertainty.", "essential_references_not_discussed": "Essential references were discussed in this work.", "other_strengths_and_weaknesses": "The experiment section is relatively weak with one single dataset and limited comparisons. It will strengthen the work if a larger dataset and more baselines can be used for stress testing and comparisons.", "other_comments_or_suggestions": "typos: \"and and\" in abstract", "questions_for_authors": "See comments above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xdiAabt6Jx", "reviewer": "Reviewer_5o6Q", "summary": "The authors proposed to numerically estimate the joint distribution with the highest variance thus leading to the ATE estimate with the lowest variance for an RCT. \n\nTheir optimization problem is modeled by all the posible choices of assignment rules but instead of optimizing over all possible assignments they optimize over all possible potential outcomes for the sample observed. This modification leads to the constraints where there is only one posible potential outcome per observation and its consistency with the ascertainment variable Z. Additionally for a random assignment the marginals of the observed vs missed outcomes must match which add an extra constraint which the authors relax a bit in order to get a more stable program. However this last constraint is relaxed to allow for general assignment probability.\n\nThe authors then proceed to prove asymptotic validity of their methods produce a true bound bye ensuring that the true distribution is always a feasible solution. They derive a high probability bound on the variance that depends on the slack factor and the sample size. the result is valid even in the case of individual confounded assignment.\n\nFinally simulators rustles are presented.\n\nWeakness\n- Is a computationally expensive approach.\n- The guarantees are for the variance, but there is no guarantees on the ATE. \n\nStrengths\n- The authors produce finite sample guarantees for their method.\n- The authors extend their method to a wide variety of estimators.\n\nOverall a very strong work", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their question about theoretical guarantees on the ATE. Our results all bound the variance of the ATE, so we believe the reviewer is asking about the bias of the ATE. Our results hold for a general class of quadratic-in-treatment estimators, which each have their own bias characteristics. This class includes estimators that are unbiased under SUTVA, such as Horvitz-Thompson. The class also includes a wide range of biased estimators, such as \"any unbiased estimator plus a constant offset.\" Because of the generality of our estimator class, we do not aim to provide general characterizations of the bias of the ATE.", "claims_and_evidence": "The theory is sound and the experiments back up the results obtained.", "methods_and_evaluation_criteria": "Refer to the summary", "theoretical_claims": "The math is sound", "experimental_designs_or_analyses": "Refer to the summary", "relation_to_broader_scientific_literature": "I am not well versed in the literature closely related to this paper.", "essential_references_not_discussed": "I am not well versed in the literature closely related to this paper.", "other_strengths_and_weaknesses": "refer to the Summary.", "other_comments_or_suggestions": "none", "questions_for_authors": "- Besides the variance are there any asymptotic guarantees for the ATE? How the computed ATE from the sample obtained as the minimizer relates with the true ATE?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduces a new method for quantifying design uncertainty in causal inference, particularly useful for small sample sizes or non-standard assignment mechanisms. Traditional bootstrap methods primarily address sampling uncertainty, while existing causal bootstrap techniques are often limited to completely randomized designs and simple estimators, such as the difference in means. The authors expand on causal bootstrap by employing integer programming to manage a wider range of estimators and assignment scenarios, including varied levels of confounding.

They validate their approach asymptotically and demonstrate its effectiveness through improved variance estimates and tighter confidence intervals in simulated geographical experiments using IMF data. This method allows for exact variance maximization without relying on closed-form copulas.

Overall, the proposed method contributes significantly to the literature on causal inference and design uncertainty. Given the importance of evaluating causal effects from confounded data, this work could have a substantial impact across various disciplines. However, since integer programming is generally computationally demanding, there are challenges in scaling the proposed method to larger datasets. Exploring efficient approximations for integer programs could be a promising future direction.
Decision: Accept (poster)