Using Surrogates in Covariate-adjusted Response-adaptive Randomization Experiments with Delayed Outcomes

Thanks for the careful review and critical comments. Below we give a point-by-point review to your comments.

• I'm confused by the claim that WHO stage is available immediately, which I took to mean that the experimenters receive the WHO stage for every patient every 6 months. Does that not require the patient to visit a clinic, which might also lead to missingness?

Re: Thanks for your question. WHO stage is a measure of the severity of HIV symptoms for the patients, and it can be collected from phone call check-in. However, clinical visits are required to measure HIV viral load. We will clarify this point in the next revision.

• Request: Neurips style is for citations to be by (Author, Year) instead of numerical. Please fix this in the next revision.

Thanks for pointing this out. We will correct the style in the next revision.

• The introduction states that "we allow $D\_{it} = \\infty$" (Line 139), but the algorithm only works with finite, known upper bound on $D$. It would be helpful for the authors to address the sensitivity of their results to the misspecification of $D^\*$.

Re: We apologize for the confusion. More rigorously, Our framework actually allows the delay to take a finite number of values, instead of having a finite upper bound. Therefore, we can also add in $\infty$ as long as there is a finite support for the rest of the values. In the next revision, we will further clarify this point.

• Weakness: The empirical results are missing error bars.

Re: Thank you for your suggestion. We have added confidence bands to the bias and standard deviation plots. See Figure 2 of the attached PDF file. From the bias plot Figure (A), we observe that the complete randomization design becomes unbiased when the sample size is over 800, while our proposed design is unbiased at a much smaller sample size. This suggests that our proposed design has a smaller finite sample bias. In Figure (B), our method demonstrates a significantly smaller standard deviation compared to the other two designs. For the "No surrogate" design (the green bands in both plots), we observe that its bias is rather significant in plot (A) as the green band fails to cover 0 when the sample size is smaller than 1600. We will further clarify this point in the next revision.

Thanks very much again for your engagement with our paper. We will incorporate your comments in the next revision.

Thanks for the fast response.

• You are right about the intuition. The idea is just to learn the observed delays and do a subtraction to learn the permanent missing (or censoring) probability. This gives an easier theoretical discussion and also provides a fully nonparametric framework. Putting a parametric form works too and is also a great idea to extend the framework to allow infinite delays. On the algorithm side, it does not change the general pipeline too much except for the delay estimation part. We also need to modify the theoretical argument accordingly.

• Oh we see your point. Yes, based on the confidence band, there is no significant difference in the mean of the bias. We do not claim that the mean of bias of the proposed method has superior properties. Our point is that the proposed design and estimator leads to a reduction of variance while maintaining the asymptotical unbiasedness property, and we show this both theoretically and empirically. This is the main message we want to convey, combining both figures.

Re your request, we hope to modify the wording as follows:

Figure (A) suggests that the point estimators from all three designs have a vanishing bias as the sample size grows, validating the asymptotic unbiasedness of all strategies. Nevertheless, in terms of variance, Figure (B) demonstrates that our proposed design yields smaller standard deviations and thus higher estimation efficiency for estimating the ATE.

Thanks again for these questions. All the points are well taken and will be incorporated in our next round of discussion.

Thanks for your response.

• Thank you for pointing this out. Apologize that we did not state the definition of $D^\star$ clearly in the paper. For the definition of $D^\star$, it is the upper bound for the finite values that delay can take. In the equation following Line 134, we briefly stated the values that $D$ can take. Moreover, Line 182 of the paper is not correct - it should be stated as $\mathcal{D} = \\{0,1,\dots,D^\star\\}\cup\\{\infty\\}$.

  From a high level, when the delay takes its value on finite support (no need to be finitely bounded), we can more easily estimate the delay distribution by taking empirical means and avoid unnecessary technicality in our theoretical discussion.

• Apologize for the ambiguity in our wording. By "unbiased" we want to make the claim of asymptotic unbiasedness of the three approaches. We fully agree that all three methods are asymptotically unbiased or asymptotically consistent for the truth, or in your words, "the mean of the bias is identical". Our claim wants to highlight a bit on the finite sample performance of the estimators in this example; or, how fast can these bias converge to zero. The confidence band suggests the "No surrogate" approach converges slightly slower, which makes sense because the surrogate information is not used in the implementation to mitigate the missingness in primary outcomes. The proposed method and complete randomization behave quite similar. Therefore, while asymptotically all three methods are unbiased, incorporating surrogates can slightly improve the finite sample convergence rate of the final estimator.

Hope the above clarification resolves your concerns. We will also incorporate the updates in the next revision.

Thanks for your careful review and critical comments. Below we add a point-by-point response to the weaknesses you listed.

• The presentation is too messy and lacks explanations in multiple places. It is unclear why EIF and variance bound are presented at the beginning of section 3.

Re: We apologize for the confusion. The semiparametric efficiency bound is the building block for our adaptive design proposal, so we use the majority of Section 3 to present it. We agree that we should try to include more intuition on why we introduce the EIF and variance bound before diving into the technical details. In the next revision, we plan to add the following explanation in Section 3:

The idea of adaptation is to adaptively collect data and adjust the allocation probability across multiple experiment stages to minimize the variance of a certain estimator for ATE. One question is: what estimator should be used for treatment effect evaluation at the end of the study? For ATE estimation, one popular choice is the semiparametric efficient estimator, whose variance matches the efficiency lower bound across all semiparametric estimators. Therefore, it is a natural idea to set the semiparametric efficiency bound as the design objective and to optimize this objective with adaptive treatment allocation. As the building block, we first present the EIF and semiparametric efficiency bound in Theorem 1 to serve as the building block for our design strategy.

• Section 4 is also presented in a poor way. The authors go through the design steps with dense notation. The authors should provide more explanation of the steps taken in the design.

Re: We apologize for the confusion. In the next revision, we will try to improve the presentation of Section 4. Here is a plan for modifying Section 4:

1. We want to add more discussion on the intuition. From a high-level perspective, our proposed design consists of four steps. The first step is from Stage $1$ to $D^\star+1$, a parameter exploration step that learns the nuisance parameters such as the delay distribution, outcome function, etc. The second step is the policy optimization step after Stage $D^\star + 1$, where the experimenter uses the estimated parameter information to compute the optimal treatment allocation strategy. The third step is the policy calibration step from Stage $D^\star + 1$ to $T$, where we calibrate the treatment probability to match the optimal treatment allocation strategy computed from the policy optimization step. The final step is constructing point and variance estimators based on the collected data.

2. We can also simplify the notation by introducing some high-level explanations for the quantity estimation part. The nuisance parameters such as outcome regressions $\hat{\tau}^{(t)}(a,x,s)$ and $\hat{\tau}^{(t)}(a,x)$ as well delay distribution $\hat{\rho}_a^{(t)}(0 | x)$ can be obtained by taking a stratified average over each combination of the treatment and covariate levels. Also, we can deliver the expression of the point estimator $\hat{\tau}$ and the variance estimator $\hat{\mathbb{V}}$ to the Appendix.

Thanks very much again for your engagement with our paper. We will incorporate your comments in the next revision.

Thanks for your careful review and critical comments. Below we add a point-by-point response.

• (W1): We agree that our efficiency bound may share some similarities with [1], but also hope to emphasize that deriving a neat decomposition of the tangent space is not straightforward when data come in multi-stages with delays (Step 2 in the proof of Theorem 1). This step is crucial for score function projection and finding the efficient influence function and requires more careful treatment in our setting than [1].

On the other hand, as you have correctly pointed out, this result is the building block for the proposed adaptive design strategy and motivates point and variance estimators for ATE. Grounded by Theorem 1, Theorem 2 and 3 further provide guarantees for the optimality of design and validity of inference. We will further clarify these points in the next revision of our manuscript.

• (W2): In HIV trials, delays can occur for many reasons. For instance, viral load is typically measured every six months due to the protocol, which are independent of outcomes. Also, factors like economic status or distance to clinics are confounders for delay and outcome. Once adjusted, these delays become non-informative. Meanwhile, delays due to severity of HIV, as you mentioned, are informative for the outcome. In reality, the interplay among these factors is complex. In our simulation, we simplify the setup by considering a hypothetical setting where delays are due to protocol or structural factors.

  In some cases, the delay is conditionally independent of the outcome. For example, if we are testing the effect of a vitamin supplement on blood vitamin levels and ask patients for regular check-ins, the delay is less likely to depend on the outcome since it does not cause severe symptoms like HIV. In general, there is a complicated delay-outcome interaction. We leave the extension as future work.

• (W3): Apologize for the confusion. The optimal allocation means the sequence of propensity scores that gives the minimal asymptotic variance matching the efficiency bound. So optimal allocation implies optimal variance. Nevertheless, the optimal treatment allocation that achieves the efficiency bound might not be unique, and one can add additional constraints or penalization (say $\ell\_2$, entropy, and so on). We will clarify this point in the next revision.

• (W4): Thanks for the suggestion. In the next revision, we will incorporate your suggestions to improve the presentation. Here is a modification plan:

  1. More discussion on the intuition. Our proposal has four steps. Step 1 is from Stage $1$ to $D^\star+1$, an exploration step that learns parameters like the delay distribution, outcome model, etc. Step 2 is the policy optimization step after Stage $D^\star + 1$, to compute the optimal allocation. Step 3 is policy calibration from Stage $D^\star + 1$ to $T$, where we calibrate the allocation to match the optimal strategy. Step 4 is constructing point and variance estimators based on the data.

  2. Simplify the notation. The elements such as outcome models $\hat{\tau}^{(t)}(a,x,s)$ and $\hat{\tau}^{(t)}(a,x)$, and delay distribution $\hat{\rho}_a^{(t)}(0 | x)$ can be obtained by taking a stratified average over each combination of the $(A,X,S)$ levels. Also, we will defer the form of $\hat{\tau}$ and $\hat{\mathbb{V}}$ to the Appendix.

  Moreover, your intuition behind variance reduction is correct. In our setting, the missing probability is different across stages and combinations of $(X, A)$ levels. We need enough data in the early stage to learn the mean and variance across arms, especially for those that are more likely to delay, so that we can optimize variance in later stages. As a toy example, consider a two-stage example, with the same $n_t$ for each stage. We omit surrogates here for simplicity. The bound in Theorem 1 becomes:

$

& \mathbb{V}\_{SPEB}
=\mathbb{E}[\big(\tau(1,X) -\tau(0,X) - \tau\big)^2]\\\\
& + 2\mathbb{E}\Big[\frac{\sigma\_1^2(X)}{ e\_1(X) \rho_1(1|X)}  + \frac{\sigma\_0^2(X)}{(1-e\_1(X)) \rho_0(1|X)}\Big]\\\\
& + 2\mathbb{E}\Big[\frac{\sigma\_1^2(X)}{ e\_2(X) \rho_1(0|X)}  + \frac{\sigma\_0^2(X)}{(1-e\_2(X)) \rho_0(0|X)}\Big]. 

$

The first stage is exploration for estimating the nuisance. We set $e_1(X) = 1/2$ to guarantee each $(x,a)$ level has enough sample. For the second stage, we minimize $\mathbb{V}_{SPEB}$ over $e_2(X)$, which gives

$

   e_2^\star(x) = \frac{\sigma_1(x)/\sqrt{\rho_1(0\mid x)}}{\sigma_1(x)/\sqrt{\rho_1(0\mid x)} 
   + 
   \sigma_0(x)/\sqrt{\rho_0(0\mid x)}}

$

Therefore, the allocation depends on $\sigma_a(x)/\sqrt{\rho_a(0\mid x)}$. If some $(x,a)$ levels lead to larger variance and fewer delays, the design will sample more from those levels to facilitate variance reduction.

• (W5): We apologize for the confusion. More rigorously, our framework allows the delay to take a finite number of values, instead of having a finite upper bound. Therefore, we can also add in $\infty$ as long as there is a finite support for the rest of the values. In the next revision, we will further clarify this point.

• (W6): One role of surrogates here is, as you have pointed out, to facilitate the derivation of the point estimator variance bound, which suggests that incorporating surrogates reduces variance (Remark 1). Another role is that it helps with policy optimization in the early stages when many primary outcomes are missing. It is interesting to think about the interplay between delay and surrogates. As a special setup, Kallus & Mao (2020) discussed missingness that depends on surrogates. We believe such an extension also works in our case and will leave this as a future work.

• Minor nits: will correct in the next revision.

Thanks again for the insightful comments! We will incorporate these points in the next revision of our manuscript.

Thank you for the thoughtful response, a quick point-by-point of my reactions. Overall I'm positive on the paper, and I'm going to raise my score to a 6.

The one point I'm still fuzzy on is (W3):

Do you need any other conditions to say something about the asymptotic variance being optimal? To clarify my question with an example: If I want to estimate $E[g(X)]$ for some deterministic function g, then obviously the empirical mean $E_n[g(X)]$ has optimal variance, but it's not usually the case that the plug-in $E_n[\hat{g}(X)]$ has optimal variance, even if $\hat{g}$ is consistent/converges in probability to $g$. Similarly, you're saying that the allocation converges in probability to the optimal one, and that if you had the optimal allocation, it would minimize variance, but I struggle to see how that's different from my example, where I could say "$\hat{g}$ converges to $g$, and using $g$ gives optimal variance, therefore my procedure is also optimal" but obviously that's incorrect. That's why I was wondering if there are other conditions, e.g., in my example above you might require that the stronger condition that $\hat{g}$ converges to $g$ at o_p(n^{-1/2)) rates so it doesn't contribute to the asymptotic variance, etc.

For the remaining points, thank you for the clarifications, my quick takes:

• (W1) That seems reasonable enough as a distinction, and I would clarify that point in the paper, as you state.
• (W2) Fair enough, I would suggest some discussion in the paper as space permits, to at least highlight reasons why the assumption might be considered a strong one.
• (W4) I like the example! Would be nice to include somewhere for intuition, if space permits.
• (W5-6) Makes sense, understood, good clarifications all around.

Got it. We agree with you about the variance claim - if we have a general $\hat{g}$ there is an additional layer of variance due to the estimation of $g$. Nevertheless, there is some specialty in the IPW and AIPW cases. Here the propensity score and the outcome models are nuisance parameters. For the outcome model, we can achieve $O(N^{-1/2})$ convergence because of the bounded moment condition on the potential outcomes and discrete covariates. Also for the delay distribution $\rho(t\mid X)$, we can achieve $O(N^{-1/2})$ because the delay indicators are bounded, and the covariate and delay have finite support. Then for the propensity score modeling, we only require $o_p(1)$ convergence. The proof is similar to the techniques adopted in double machine learning literature. To see this, consider the following math in the classical AIPW case:

$

&\frac{1}{n}\sum_{i} \frac{(Y_i - \hat{\mu}(X_i)) A_i}{\hat{e}(X_i)} + \frac{1}{n}\sum_{i} \hat{\mu}(X_i)-\mu\\ =&\frac{1}{n}\sum_{i} \frac{{e}(X_i)}{\hat{e}(X_i)} \frac{(Y_i - {\mu}(X_i)) A_i}{{e}(X_i)} + \frac{1}{n}\sum_{i}\mu(X_i) - \mu \tag{1}\\ +&\frac{1}{n}\sum_{i} \{\hat{\mu}(X_i) - \mu(X_i) \}\{ 1 - \frac{A_i}{\hat{e}(X_i)}\}. \tag{2} \\

$

For (1), we show this part gives the asymptotic minimal variance. We have

$

&\frac{1}{n}\sum_{i} \frac{{e}(X_i)}{\hat{e}(X_i)} \frac{(Y_i - {\mu}(X_i)) A_i}{{e}(X_i)} + \frac{1}{n}\sum_{i}\mu(X_i) - \mu \\ =& \frac{1}{n}\sum_{i} (1 + o_p(1)) \frac{(Y_i - {\mu}(X_i)) A_i}{{e}(X_i)} + \frac{1}{n}\sum_{i}\mu(X_i) - \mu \\ \asymp & \frac{1}{n}\sum_{i} \frac{(Y_i - {\mu}(X_i)) A_i}{{e}(X_i)} + \frac{1}{n}\sum_{i}\mu(X_i) - \mu + o_p(n^{-1/2}).

$

For (2), we show it's of small order $o(n^{-1/2})$.

$

&\frac{1}{n}\sum_{i} \{\hat{\mu}(X_i) - \mu(X_i) \}\{ 1 - \frac{A_i}{\hat{e}(X_i)}\}\\ = & \sum_{x\in\mathcal{X}}\frac{\hat{\mu}(x) - \mu(x) }{\hat{e}(x)}\frac{1}{n}\sum_{i} \{A_i - \hat{e}(x)\}1\{X_i = x\}\\ = & \sum_{x\in\mathcal{X}}\frac{\hat{\mu}(x) - \mu(x) }{\hat{e}(x)} \cdot \frac{1}{n}\sum_{i} \{A_i - {e}(x)\}1\{X_i = x\} +\sum_{x\in\mathcal{X}}\frac{(\hat{\mu}(x) - \mu(x))(e(x) - \hat{e}(x)) }{\hat{e}(x)}\frac{1}{n}\sum_{i}1\{X_i = x\}\\ = & \sum_{x\in\mathcal{X}} \frac{O_p(n^{-1/2})}{e(x) + o_p(1)} \cdot O_p(n^{-1/2}) +\sum_{x\in\mathcal{X}} \frac{O_p(n^{-1/2}) \cdot o_p(1)}{e(x) + o_p(1)} \cdot \{p(x) + O_p(n^{-1/2})\}\\ = & o_p(n^{-1/2}).

$

Therefore, we can see that (1) plays the dominant role in the asymptotic regime and gives minimal variance. Again, we highlight that the specialty is due to we can estimate the outcome model in square root order so we do not need this order for the propensity score models.

Hope this clarifies your concern regarding the variance part.

Thanks so much for your careful review. We add a point-by-point response to your questions below.

• It would be useful to have a finite sample analysis.

Re: Thanks for the suggestion. We pursue asymptotic analysis to establish the distributional convergence and construct confidence intervals based on normal distribution approximations. In settings like A/B testing, many users are involved and asymptotic analysis suffices to help with inference. Nevertheless, we agree that in other settings such as clinical trials, sample size becomes a more important constraint which calls for a more delicate finite sample quantification. To translate the asymptotic results into a finite sample analysis, we can utilize the martingale concentration inequalities to quantify the tail of the estimates and also use Berry-Esseen bounds to obtain rates of distribution convergence. To avoid an overly lengthy paper, we hope to save these for future work.

• It would have also been nice to have seen a slightly larger set of empirical results.

Re: Thanks for the suggestion. In the rebuttal, we have added a new experiment to compare different tuning strategies for our design optimization program (Figure 1 in the attached pdf). Also, we have refined the presentation of the previous experiments by adding confidence bands to the curves (Figure 2 in the attached pdf).

• It wasn’t clear to me what is being assumed about the size of each wave of participants.

Re: In our setting, we assume that in each stage, the portions of participants to be recruited in each stage (the $r_t$'s) are decided a priori. This aligns with the practice in many real-world applications. For example, in clinical trials, researchers need to decide the number of patients to enroll as part of the protocol before the trial starts. Theoretically, this assumption makes it easier to discuss asymptotic regimes when $N$ goes to $\infty$.

• Also these participants are arriving i.i.d. over time?

Re: In our setting, we assumed that the potential outcomes of the participants across stages are i.i.d., but the truly observed outcomes are not. The i.i.d. assumption is only imposed on the potential outcomes, which can be interpreted as the patients are drawn from the same target population. Nevertheless, due to the adaptive recruitment process, the observed outcomes are decided by both the distribution of the potential outcomes and the history-dependent data collection policy and thus are not i.i.d. anymore.

• Does the proposed procedure give any guidance to the setting where the experimenter can choose the number of participants to enroll at each stage?

Re: Thanks for this insightful point. Our results can provide some guidance in the following two aspects:

1. For Stage $1$ to $D^\star$, Condition (4) in Theorem 3 provides a sufficient condition to quantify the portion of people to recruit for the first $D^\star$ stages to achieve the optimal treatment allocation strategy, which depends on the variance and the delay mechanism. Therefore, a promising strategy is to use the first stage to estimate the size of the quantity on the right-hand side of Condition (4), then adjust the portions in the following stages to satisfy the constraint.

2. For Stage $D^\star + 1$ to $T$, we apply the calibrated treatment allocation probabilities $\hat{e}\_l$. To make sure $\hat{e}\_l$ is bounded away from $0$ and $1$, we could set up a sequence of $n_s$ that satisfy the following inequality:

$$\delta \le \frac{\sum\_{s=1}^{T-D^\star} n_s \cdot \tilde{e}\_{s} - \sum_{s=1}^{D^*+1} n\_s\cdot \hat{e}\_{s}}{\sum\_{s = D^* + 2}^{T-D^*}n\_{s}} \le 1 - \delta.$$

This will ensure a sufficient portion of participants in both the treatment and control groups for Stage $D^\star + 1$ to $T$.

Thanks again for the insightful comments! We will incorporate these points in the next revision of our manuscript.

Thank you very much for your careful review and critical questions. Below, we add a point-by-point response to your questions:

• In the numerics, the "No Surrogate" approach seems to perform worse than complete randomization. Is there any insight on why this adaptive approach performs worse than the complete randomization approach? Does the "No Surrogate" approach also account for delays or does it utilize an existing adaptive randomization design approach?

Re: Sorry for the confusion. The complete randomization strategy incorporates surrogates and adjusts for delays in the primary outcome. The "No Surrogate" strategy accounts for delays but does not incorporate any surrogates. In general, these two strategies are not directly comparable, and the point we want to emphasize here is that both strategies are inferior to the proposed CARA design. Between complete randomization and the "No surrogate" strategy, which performs better depends on the surrogates' quality. If the surrogates are more predictive of the outcome, then incorporating the surrogates can greatly reduce the bias and variance, leading to an estimator that outperforms the "No Surrogate" strategy. According to the summary statistics in Table 2 in the Appendix, different levels of the WHO HIV stage lead to a significant variation in the conditional mean and variance of primary outcomes, which suggests the importance of incorporating surrogates and explains the superiority of this strategy in this case.

• The proposed method suggests applying regularization to the oracle optimization program to deal with potential multiple optima. Does the procedure tune the regularization parameter so it produces calibrated treatment allocation probabilities with better estimation error? Does this align with the cross-validation set-up briefly mentioned in the paper? If so, what is the exact cross-validation setup and if not how is the the tuning parameter chosen?

Re: Thank you for your question. We provide a comparison of various tuning parameter selections in Figure 1 of the attached pdf file. Overall, our design is not sensitive to tuning parameter choices, especially when the sample size is large. Based on Figure 1, as heuristic guidance, we recommend choosing a tuning parameter between 0 and 1 because a large tuning parameter may push the treatment allocation solutions to the boundary. As demonstrated in Figure 1, when the tuning parameter is chosen to be 10, the design is less efficient than the performance under a smaller tuning. We agree that having a cross-validation procedure for tuning parameter selection would be more interesting and systematic. We hope to explore this piece of methodological development in future work.

Thanks again for your engagement with our paper. We will incorporate your comments in the next revision.