Integer Programming for Generalized Causal Bootstrap Designs

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.

We thank the reviewer for their thoughtful engagement with our work.

Responses to Questions

1. Thank you for this question; please see the discussion of scalability in our response to Reviewer bwtT.

2. 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.

3. 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.

4. 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.

5. 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.

6. 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.

7. 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.

8. 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.

9. 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.

10. 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.

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.

We thank the reviewer for their careful and positive review of our paper. Below, we address the two questions raised.

Q1: 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].

Bounded 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).

We would be happy to expand on these assumptions and include them in the paper if the reviewer would find it helpful.

[1] Sävje, Fredrik. "Randomization does not imply unconfoundedness." arXiv preprint arXiv:2107.14197 (2021).

Q2: 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.

The 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.

To scale the IP beyond hundreds of units, we see potential to apply relaxation or approximation algorithms, as also suggested by Reviewer XXTL.

For 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.

For 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.