---
author:
- |
    Steffen Ventz$^{1,2,\star}$[^1], Matteo Cellamare$^{1,2,\star}$[^2], Sergio Bacallado$^{3}$[^3], and Lorenzo Trippa$^{1,2}$[^4]\
    $^{1}$[Dana-Farber Cancer Institute, US]{}, $^{2}$[Harvard T.H. Chan School of Public Health, US]{},\
    $^{3}$[University of Cambridge, UK]{}\
    $\star$[co-first authorship]{}
bibliography:
- 'Bibliography.bib'
title: 'Supplementary Material - Bayesian Uncertainty Directed Trial Designs'
---

Additional Simulation Results {#Additional:Sim}
=============================

[ Consider two studies, A and B, [ both]{} with five treatment arms and normally distributed outcomes with variance $\sigma^2=1$ and unknown means $\mu_a$, $a=1, \cdots, 5$. Study A seeks to identify $a^\star = \arg \max_{a} \mu_a$, and attempts to minimize the discrepancy between $\mu_{a^\star}$ and $\mu_{\widehat{a}^\star}$, where $\displaystyle \widehat{a}^\star=  \arg \min_{a=1, \cdots, 5} E[ (  \mu_{a} - \mu_{a^\star} ) ^2 | \Sigma_t ]$. An uncertainty metric consistent with this goal is $E[(\mu_{\hat{a}^\star} - \mu_{a^\star} )^2| \Sigma_t ].$ In contrast, the goal of study B is the estimation of the means $\mu_a$, $a=1,\cdots, 5$. The goals are different and an appropriate design for one study might be inadequate for the other. For study B a balanced design is consistent with the uncertainty measure $\sum_a \mbox{Var}(\mu_a|  \Sigma_t   )$. But the balanced design might be inadequate for study A. Consider for example a scenario with means $( 1.2,1, 0.2, 0, -0.2)$. For study A, in this case, a design that progressively eliminates candidate arms performs considerably better. With $80$ patients, a sequential design which after $10, 20, \dots, 70$ enrollments drops $a$ if $\displaystyle p({\mu}_a < \max_{1\leq \ell \leq 5} {\mu}_\ell | \Sigma_t ) > 0.95$, where [*a priori*]{} $\mu_a \overset{iid}{\sim}N(0, 1)$, selects $a^\star$ with a probability 5% higher than the balanced design. ]{}

Simulation Results for Normal Distributed Outcomes {#Sec:Normal:Example .unnumbered}
--------------------------------------------------

We present additional simulation results for a multi-arm trial with a BUD design as described in Example 1 of Section 3.1 in the main manuscript. Outcomes for patients assigned to arms $a=0, \cdots, K$ are $N(\theta_a, \sigma^2_a)$ distributed with known variances $ \sigma^2_a$ and unknown means $\theta_a \sim N(0, \sigma^2_0)$. The BUD information measure coincides with the sum of posterior variances $u(\Sigma_t)= -\sum_{k=0}^K Var(\theta_a|\Sigma_t)$. Table \[Tab:Normal\] summarizes selected operating characteristics for a trial with $K=3$ experimental arms and $T=150$ patients in four scenarios with $h(t)=3$. We consider two scenarios with identical variances $\sigma^2_a=1$ across arms and two scenarios with different variances $\sigma^2_a=2,2,1.5, 0.5$ across arms. For each scenario Table \[Tab:Normal\] shows the asymptotic approximation of the patient allocation (expression \[limit\] in Section \[Sec:Asymptotic\] multiplied by $T=150$). Table \[Tab:Normal\] indicates that the average allocation per arm agrees closely with the asymptotic approximation.

Simulation Results for Section \[Sec:Multi:Arm1\] {#simulation-results-for-section-secmultiarm1 .unnumbered}
-------------------------------------------------

BUD designs for controlled multi-arm trials can be defined for studies with different aims. For instance one may consider the minimization of uncertainty on the treatment effect of (i) the most effective experimental arm, (ii) all experimental arms with positive treatment effects or (iii) all experimental arms. We compared BUD designs with information functions $u_i, i=1,2,3$ ($BUD_i$) defined by distinct study aims, $$\begin{aligned}
u_1(\Sigma_t) &= - \mbox{Var}( \max_a \gamma_{a}    \mid  \Sigma_t ) \label{s1},\\
u_2(\Sigma_t) &= \sum_{a=1}^K \Big[ \mbox{Var}(\gamma_a I(\gamma_a >0) ) - \mbox{Var}( \gamma_a I(\gamma_a>0)   \mid  \Sigma_t ) \Big], \label{s2}\\
u_3(\Sigma_t) &= \sum_{a=1}^K \Big[ \mbox{Var}(\gamma_a)  - \mbox{Var}( \gamma_a   \mid  \Sigma_t ) \Big], \label{s3}\end{aligned}$$ for a trial with two experimental arms $a=1,2$. In the four scenarios that we consider, the probabilities of response $(\theta_0, \theta_1, \theta_2)$ to arms $a=0,1,2$ equal $(0.6, 0.5, 0.8)$, $(0.6, 0.5, 0.7), (0.6, 0.4, 0.9)$ and $(0.5, 0.6, 0.7).$ The overall sample size $T=3 \times 97$ was selected based on targeted arm-specific type I and II error bound of $0.05$ and $0.2$ for $(\theta_0,\theta_a)=(0.6, 0.8)$.

Table \[TAB:Inferior:utility\] summarizes selected operating characteristics of the three BUD designs. Compared to $BUD_3$, a BUD design defined by $u_2$ focuses on positive treatment effects and reduces the average number of assignments to inferior agents (Scenarios 1 to 3 of Table \[TAB:Inferior:utility\]) by 58 to 77 patients (34, 37, 20 enrollments on average for $BUD_2$ compared to 97 under BR). A BUD design defined by $u_1$ focuses on the minimization of uncertainty on the treatment effect of the most effective arm, and, compared to $BUD_2$, reduces the number of patient assignments to experimental arm one on average by 1 to 9 patients in scenarios 1-3. But when multiple effective experimental arms exists (scenario 4) $BUD_1$ has $10\%$ lower power for arm 1 than $BUD_2$.

Additional Figures and Tables
=============================

Supplementary Material for Section \[Sec:Multi:Arm1\]
-----------------------------------------------------

![ Multi-arm trial, as discussed in Section \[Sec:Multi:Arm1\] of the manuscript. Regret function $E_{d^\star }[u(\Sigma_T)] - E_{d}[u(\Sigma_T)]$ between the optimal sequential Bayesian design $d^\star$, obtained through backwards induction (BI), and either a Bayesian uncertainty directed design $d=d_{BUD}$ (green curve), a balanced randomized design $d=d_{BR}$ (blue curve), Bayesian adaptive randomization $d=d_{BAR}$ (red curve), or the randomized play-the winner design $d=d_{RPW}$ (purple curve). Results are based on the exact utility $E_{d^\star }[u(\Sigma_T)]$ of $d^\star$ under a uniform prior for the outcome probabilities $\theta = (\theta_1, \dots, \theta_4)$. For $d=d_{BUD}, d_{BAR}, d_{BR}, d_{RPW}$ we computed the expected utility $\mathbb E_{d}[u(\Sigma_T)]$ based on 100,000 simulated trials with outcome probabilities $\theta = (\theta_1, \dots, \theta_4)$ generated under the prior. Exact computations based on BI become computationally unfeasible for larger sample sizes in the order of hundreds of patients. []{data-label="Bandit"}](figures/FigSec3p1Regret){width="15cm" height="8cm"}

Supplementary Material for Section \[Section:FindBestArm\]
----------------------------------------------------------

Supplementary Material for Section \[Sec:SubgroupTrial\]
--------------------------------------------------------

Proofs
======

We show that for any $\varepsilon>0$, the process $ \widehat p_{0,t}$ can only visit the interval $[\rho_0+\varepsilon,1]$ finitely often. A symmetric argument can be applied to intervals lying below $\rho_0$. Define the martingale [$$\begin{aligned}
M_{t} =  \widehat p_{0,t} - \sum_{i=1}^t {\mathbb{E}}[\widehat p_{0,i}-\widehat p_{0,i-1}\mid \Sigma_{i-1}] .\end{aligned}$$]{} As $M_{t+1}-M_{t}\leq \mathcal O(1/t)$, Doob’s $L^2$ martingale convergence theorem implies that $M_{t}\to M_{\infty}<\infty$ a.s. Furthermore, a simple derivation yields [$$\begin{aligned}
{\mathbb{E}}[\widehat p_{0,t+1}- \widehat p_{0,t}\mid \Sigma_t]= 
-\frac{ \widehat p_{0,t}}{t+1} + \frac{1}{t+1}\frac{\Delta_{t}^h(0)}{\Delta_{t}^h(0)+\Delta_{t}^h(1)} = \frac{F_t}{t+1}.\end{aligned}$$]{} In order to prove that $ \widehat p_{0,t}=M_{t}+ \sum_{i=1}^t {\mathbb{E}}[\widehat p_{0,i}-\widehat p_{0,i-1}\mid \Sigma_{i-1}] $ only visits $[\rho_0+\varepsilon,1]$ finitely often, we will use the convergence of the first term, and the fact that the summands in the second term define a drift toward $\rho_0$ for large values of $i$.

On a set of probability 1, there is a number $N$ such that for any $t>N$, $| M_{t}-M_{\infty}|<\varepsilon/4$, and $F_t<-c<0$ when $ \widehat p_{0,t}\in[\rho_0+\varepsilon/2,1]$. Therefore, the trajectory $(\widehat p_{0,N+k})_{k\geq 1}$ cannot increase more than $\varepsilon/4$ while staying in $[\rho_0+\varepsilon/2,1]$. As $\sum_{t=1}^\infty(t+1)^{-1}=\infty$, every time the process enters this interval, it will eventually exit below $\rho+\varepsilon/2$. Finally, observe that $|F_t|\leq 2$ and thus $|F_t/(t+1)|\to 0$, so for $t$ large enough, the process cannot reenter the interval $[\rho+\varepsilon/2,1]$ above $\rho+3\varepsilon/4$, so it can only enter $[\rho+\varepsilon,1]$ finitely often.

Due to the conjugacy between the Bernoulli likelihood and the Beta prior, the variables $\theta_0,\dots,\theta_K$ are independent and Beta a posteriori. This makes it possible to write down the information gain in closed form. Letting $\hat\theta_{a,t}$ be the proportion of successes in arm $a$ by time $t$, [$$\begin{aligned}
\Delta_{t}(a) =\;& -\frac{(\alpha+t  \widehat p_{a,t}\hat\theta_{a,t}+1)(\beta+t  \widehat p_{a,t}(1-\hat\theta_{a,t}))}{(\alpha+\beta+t  \widehat p_{a,t}+1)^2(\alpha+\beta+t  \widehat p_{a,t}+2)} 
\left( \frac{t \widehat p_{a,t}\hat\theta_a+\alpha}{t \widehat p_{a,t}+\alpha+\beta} \right)\\
& - \frac{(\alpha+t  \widehat p_{a,t}\hat\theta_{a,t})(\beta+t  \widehat p_{a,t}(1-\hat\theta_{a,t})+1)}{(\alpha+\beta+t  \widehat p_{a,t}+1)^2(\alpha+\beta+t  \widehat p_{a,t}+2)} 
\left( \frac{t \widehat p_{a,t}(1-\hat\theta_a)+\beta}{t \widehat p_{a,t}+\alpha+\beta} \right)\\
& +\frac{(\alpha+t  \widehat p_{a,t}\hat\theta_{a,t})(\beta+t  \widehat p_{a,t}(1-\hat\theta_{a,t}))}{(\alpha+\beta+t  \widehat p_{a,t})^2(\alpha+\beta+t  \widehat p_{a,t}+1)} \\
=\;& \frac{1}{t^2}\left(\frac{\hat\theta_{a,t}(1-\hat\theta_{a,t})}{ \widehat p_{a,t}^2}\right) + \mathcal O((t\widehat p_{a,t})^{-3}). \\
=\;& \frac{1}{t^2}\left(\frac{\theta_{a,t}(1-\theta_{a,t})}{ \widehat p_{a,t}^2}\right) +\mathcal O(|\hat\theta_{a,t}-\theta_a|(t\widehat p_{a,t})^{-2})+ \mathcal O((t\widehat p_{a,t})^{-3}).\end{aligned}$$]{} where the last two equalities hold due to Taylor’s theorem. Now, we can check the condition of Lemma \[asymptotics::lemma\] by the same argument used in the normal outcomes case. Writing $\sigma_a=\theta_a(1-\theta_a)$, the Taylor series approximation implies [$$\begin{aligned}
F_t = - \widehat p_{0,t}+\frac{ \left(\widehat p_{0,t}^{-2}\sigma_0^{2}
+\mathcal O(|\hat\theta_{0,t}-\theta_0|\widehat p_{0,t}^{-2})+ \mathcal O(t^{-1}\widehat p_{0,t}^{-3})\right)^h}
{\left(\widehat p_{0,t}^{-2}\sigma_0^{2}
+\mathcal O(|\hat\theta_{0,t}-\theta_0|\widehat p_{0,t}^{-2})+ \mathcal O(t^{-1}\widehat p_{0,t}^{-3})\right)^h
+ \left(\widehat p_{1,t}^{-2}\sigma_1^{2}
+\mathcal O(|\hat\theta_{1,t}-\theta_1|\widehat p_{1,t}^{-2})+ \mathcal O(t^{-1}\widehat p_{1,t}^{-3})\right)^h }.\end{aligned}$$]{} We show that the residual terms in $\mathcal O(\cdot)$ are negligible in the limit, which is trivial when $\widehat p_{0,t}$ and $\widehat p_{1,t}$ are of constant order. For any $\delta\in(0,1/4)$, on the subsequence where $\widehat p_{a,t}> \max\{(t\widehat p_{a,t})^{-1/4+\delta} \allowbreak (\log\log (t\widehat p_{a,t}))^{1/2}, \allowbreak t^{-1/3+\delta}\}$ for $a=0,1$, by the law of the iterated logarithm, the residual terms in the above expression go to 0 a.s. as $t\to\infty$. Otherwise, one of $\widehat p_{0,t}$ or $\widehat p_{1,t}$ goes to $0$ as $t\to\infty$. We conclude that, on a set of probability 1, $F_t$ converges to [$$\begin{aligned}
\tilde F_t = - \widehat p_{0,t}+\frac{ \widehat p_{0,t}^{-2h}\sigma_0^{2h}}
{\widehat p_{0,t}^{-2h}\sigma_0^{2h}
+ (1-\widehat p_{0,t})^{-2h}\sigma_1^{2h}.}\end{aligned}$$]{} The argument used in the normal model allows us to apply Lemma \[asymptotics::lemma\] and Remark \[multiarm::remark\], which imply the a.s. convergence $\widehat p_{a,t}\to\rho_a\propto \sigma_a^{2h/(2h+1)}$ for a BUD design with $K$ arms.

 Computational Details  {#Sec:Computational:Details}
=======================

The computing time to simulate a BUD trial increases linearly with respect to the sample size $T$ and with respect to the number of possible actions. In our examples the computing time for BUD simulations, as expected, is very similar to BAR, and significantly faster than computing the optimal Bayesian design with dynamic programming.

[**Computational details for Section \[Exa1\]:** ]{} In the first multi-arm trial designs of Section \[Exa1\], the utility $
u(\Sigma_t) = \sum_{a=1}^K ( v_a - {\textnormal{Var}}(\gamma_a | \Sigma _t)  )
$, with independent and conjugate priors Beta$(1, 1 )$ for the response rates of each arm, has a closed form. The expected utility conditional on $A_{t+1}=a$ is computed with the formula $E[ u(\Sigma_{t+1}) \mid \Sigma_t, A_{t+1}=a ] = 
u_{a,1}\widehat \theta _a(t) + 
u_{a,0}(1- \widehat \theta _a(t) ),$ where $\widehat \theta _a (t) = E[\theta_a \mid \Sigma _t]$, and $u_{a,0}$ and $u_{a,1}$ are the utilities at time $t+1$ with a negative and positive outcome, respectively. In the second example we computed the negative entropy $u(\Sigma_t) = E[ \log p_{\theta_{a^\star}} (\theta_{a^\star} \mid \Sigma_t)\mid \Sigma_t ]$ for the posterior of the most effective arm $a^\star$ using numerical integration.

[**Computational details for Section \[Sec:SubgroupTrial\]:** ]{} The prior specified on the response rates $\theta$ makes it possible to compute the posterior probability of treatment effects $E_a=(E_{0,a},E_{1,a})$ in arm $a$ in closed form, which yields an exact expression for the information gain $\Delta_a(t)$. For simplicity, consider the case with one control and one experimental arm, $a=0,1$. Here, [$$\begin{aligned}
p(E_1\mid Y_{1:t}) = \int p(Y_{1:t}\mid \theta) p(\theta\mid E_1) p(E_1) d\theta\end{aligned}$$]{} where $\theta= (\theta_{0,0},\theta_{1,0},\theta_{0,1},\theta_{1,1})$. Inside the integral, $p(E_1)$ is a simple function of $\pi$ and $\lambda$, while $p(\theta \mid E_1)$ is a product Beta density with integer parameters restricted to a set $R\subset [0,1]^{4}$. The normalizing constant of this truncated density can be computed exactly using a well-known combinatorial formula, derived by representing Beta random variables as order statistics of uniform random variables. By conjugacy, the product $p(Y_{1:t}\mid \theta)p(\theta\mid E_1)$ is proportional to a product of beta densities, and its integral over $R$ is obtained using the combinatorial formula mentioned.

For each treatment $a=1, \dots, K$ we tested treatment effects in the biomarker-targeted population (all patients with biomarker characteristic $X_{t,b_a}=1$) and in the biomarker-negative subgroup (patients with $X_{t,b_a}=0$) using a parametric bootstrap test (see [@ventz2017combining] for details).

[**Computational details for Section \[Sec:CoPrimary\]:** ]{} In the BUD design for multi-arm trials with two co-primary endpoints, we estimated the posterior probability of treatment effects by Monte Carlo integration, drawing independent samples from the Dirichlet posterior of the response rates $\theta$. The expected utility is then computed by $E[ u(\Sigma_{t+1}) \mid \Sigma_t, A_{t+1}=a ] =  \sum_{y \in \{0,1\}^2} u_{y, a}(t+1) \times \widehat \theta _{y, a}(t),$ where $\widehat \theta _{y, a}(t) = E[ \theta _{y, a}\mid  \Sigma_t ]$ and $u_{y, a}(t+1)$ indicates the expected utility $u(\Sigma_{t+1})$ given $\Sigma_t$ and $(A_{t+1}, Y_{t+1})=(a,y)$.

[^1]: steffen@jimmy.harvard.edu

[^2]: matteoc@jimmy.harvard.edu

[^3]: sb2116@cam.ac.uk

[^4]: ltrippa@jimmy.harvard.edu
