汤普森（Thompson）、乌拉姆（Ulam）还是高斯（Gauss）？在具有二元结局的贝叶斯响应自适应试验中，针对后验概率计算方法的多标准推荐选择

《Computational Statistics & Data Analysis》：Thompson, Ulam, or Gauss? Multi-criteria recommendations for posterior probability computation methods in Bayesian response-adaptive trials with binary endpoints

　　丹尼尔·卡达吉（Daniel Kaddaj）| 斯特夫·巴斯（Stef Baas）| 埃德温·Y.N. 唐（Edwin Y.N. Tang）| 大卫·S. 罗伯逊（David S. Robertson）| 卢卡斯·平（Lukas Pin）| 索菲娅·S. 维拉尔（Sofía S. Villar）
剑桥大学纯数学与数理统计系，数学科学中心，威尔伯福斯路（Wilberforce Road），剑桥（Cambridge），CB3 0WB，英国

**摘要**
贝叶斯自适应设计通过根据累积的数据调整试验特征，使得临床试验更加灵活。其中的贝叶斯响应自适应随机化（Bayesian Response-Adaptive Randomisation，BRAR）方法会根据中期数据将患者分配到更有前景的治疗方案中。实施BRAR需要相对快速地评估后验概率。然而，现有封闭形式解的局限性意味着试验通常依赖于计算密集型近似方法，这可能会影响准确性和探索的情景范围。虽然存在更快的高斯近似方法，但其可靠性并不保证。关键问题是，所使用的近似方法往往报告不足，且相关文献缺乏选择和比较这些方法的实际指导，特别是在计算速度、推断准确性以及对患者益处的影响之间的权衡方面。

本文聚焦于具有二元终点的BRAR试验，开发了一种新颖算法，能够高效且精确地计算这些后验概率，从而对现有的近似方法进行稳健评估。利用这些精确计算结果，建立了一个全面的基准体系，用以评估不同方法的计算速度、患者效益和推断准确性。通过一系列模拟以及对一个真实的多臂试验的重新分析，发现精确计算算法在小到中等数量的试验臂数情况下通常是最快的方法。此外，研究表明常用的近似方法可能导致显著的统计功效损失和I型错误率增加，尤其是高斯近似方法。为了解决精确多臂计算中的指数级增长问题，本文提供了一个正式的定量框架。这一实用指南为实践者提供了明确的阈值，以便在各种临床试验设置中无缝切换精确计算和安全的近似方法。

**1. 引言**
自适应设计可以通过允许根据累积数据预先计划调整来提高临床试验效率，同时保持统计有效性（Burnett等人，2020年），并已广泛普及（Lee等人，2023年）。贝叶斯自适应设计通过后验优势概率（Posterior Probability of Superiority，PPS）来驱动这些调整。早期的一个例子是贝叶斯响应自适应随机化（BRAR），它根据患者的PPS顺序将他们分配到相应的治疗方案中（Thompson，1933年）。历史上，在没有严格假设的情况下，这些概率的精确计算（即使用封闭公式）是不可行的。

目前，贝叶斯设计应用范围从剂量寻找到确认性试验。例如，持续重新评估方法（Continual Reassessment Method，CRM）（O’quigley等人，1990年）使用二元终点的后验概率来进行剂量分配；而在肿瘤学中，BRAR已被提出用于在最大耐受剂量以下的患者导向剂量分配（Pin等人，2024b年）。后验概率还用于确定疗效无效的早期终止规则（Gsponer等人，2013年）。然而，由于现有封闭形式解的局限性，试验通常依赖于计算密集型模拟来近似PPS（例如Berry等人，2010年）。

尽管BRAR经常使用二元终点——在最近的一项回顾中占响应自适应试验的54.9%（Pin等人，2024a年；Wilson等人，2025年）——但精确的PPS计算常被忽视。目前已知存在适用于常见Beta先验的精确封闭公式（参数为整数值），但仅限于最多四个治疗臂数，超过两个治疗臂数时计算变得非常繁琐。因此，尽管精确方法在小样本设置中具有关键优势，但仍未得到充分利用。此外，关于选择近似方法没有标准软件或指导，计算技术也常常报告不足：例如，Bleck等人（2013年）和Richter等人（2023年）仅按需求提供代码；Guglielmetti等人（2025年）指出需要依赖独立统计学家；Alexander等人（2018年）则未提供任何信息。

本文解决了这些计算挑战。我们提出了一种新颖算法，可以在与试验规模成线性关系的时间内精确计算任意数量治疗臂的PPS。然后，我们使用这种精确方法对三种常见的近似方法进行基准测试：蒙特卡洛模拟（例如Carlson等人（2017年）、数值积分（例如Yannopoulos等人（2020年））和高斯近似（例如Marian等人（2024年）或Kamber等人（2025年））。我们证明了我们的精确算法在计算上具有显著优势——对于两个治疗臂，其速度比高斯近似快12至60倍，比蒙特卡洛方法快257至2290倍。这种效率有助于进行稳健的试验模拟，防止不准确的功效估计（例如Conor、Elm和Broglio（2013年））和不足的I型错误率控制（Baas、Jacko和Villar，2025b年）。就准确性而言，我们发现高斯近似方法会增加I型错误率并降低功效；即使进行10,000次迭代，蒙特卡洛方法也可能导致功效损失，这凸显了适当计算方法的关键性。

本文的结构如下：第2节介绍模型和精确计算算法；第3节分析计算成本；第4节展示2臂试验的模拟研究和准确性分析；第5节评估一个实际的多臂试验和一个5臂重新设计的案例；最后第6节进行总体讨论，提出严格的方法比较框架和多臂设计的实际指导。

**2. 方法**
**2.1. 介绍性符号**
使用Robertson等人（2023年）和Baas等人（2025b年）的符号，我们考虑最多有n名患者的临床试验，患者依次被分配到k个治疗臂中的一个（编号为0到k?1）。患者i被分配到臂j的概率用指示符Ai,j表示，他们的二元反应用Yi～Bernoulli(pai)表示，其中ai是患者i被分配到的治疗。直到患者i的所有试验历史记录用Hi表示：=(A1,Y1,…,Ai,Yi)，Ni,j和Si,j分别表示臂j的总患者数和阳性反应数。

一种简单的贝叶斯响应自适应随机化（S-BRAR）程序是一种随机分配程序（Thompson，1933年），它根据当前数据将患者按给定治疗是最优选项的概率比例分配。给定响应概率的先验p：Q:=Q0×…×Qk?1，以及给定试验历史Hi的p的后验分布Qi，分配概率为：
(1) πS-BRAR,jQ(Hi) = Qi * (pj > maxj′ ≠ pj′)，
其中i*是决定随机化概率更新频率的批次大小b的最大整数倍，且i < i。实际上，Qi可以很好地分解，仅依赖于Hi、Ni,j和Si,j，从而简化了计算（见引理A.1）。注意，这些随机化概率的更新频率……

我们将检验无唯一最优治疗的零假设，对抗个别更优治疗的备择假设，即
H0: |argmax{p0,…,pk?1}| > 1；
HAj: pj > maxj′ ≠ pj′。

如Tang等人（2025年）在BRAR试验中常见的做法，我们使用后验概率作为早期终止和最终分析的检验统计量。我们定义一个拒绝规则R=?i=1≤n?j=0≤k?1Rji，其中Rji?{0,1}(k+1)，i是在患者i之后试验停止并支持假设HAj的试验历史集合。这些区域基于检验统计量Tj(Hi):=Qi(pj > maxj′ ≠ pj′)定义。对于一组预先计划的分析点P?{1,?,n}，如果maxi∈PTj(Hi) > c，则试验提前终止并支持假设HAj，其中P?{1,?,n}，c是选定的阈值概率（详见附录A.5.1）。此类规则在实践中已被广泛应用（Bleck等人，2013年；Carlson等人，2017年）。由于患者i之后的检验统计量正好与患者i+1的随机化概率相匹配，我们可以简单地通过为一个患者扩展试验来计算所有必要的检验统计量。

假设Beta先验是独立的，后验分布仍然是仅依赖于先验和观察到的反应的独立Beta分布（引理A.2）。为了简化符号，假设所有j的aj、bj都是整数值，我们希望计算的后验概率可以写成P(x;{i}):=P(Xi > maxj ≠ iXj)，其中Xj～Beta(x2j,x2j+1)对所有j都是独立的，且x∈N2k。这种符号的合理性在于，我们希望表示当S?{0,1,?,k?1}且xS′是从S中移除了x2s和x2s+1的向量时，通过合并S中的臂获得的复合臂更优的概率：
P(x;S):=P((∑s∈Sx2s, ∑s∈Sx2s+1,xS′);{0})。

**2.2. 计算随机化概率**
计算随机化概率通常不直接，往往需要近似。下面我们分别评估精确方法和近似方法。

**2.2.1. 精确评估**
对于k=2的情况，Liebermeister（1877年）和Thompson（1933年）推导出了P(x; {j})的精确公式，Miller（2014年）最近提出了一个计算上更易处理的版本：
(2) P(x;{1})=∑i=0x2?1B(x0+i,x1+x3)(x3+i)B(i,x3)B(x0,x1)，
其中B是Beta函数。利用对称性，该求和可以简化为仅包含min(x0, x1, x2, x3)项的形式。

对于任何k≥2的情况，我们提出了一种新颖的递归方法（算法1），用于高效精确计算这些概率：
**定理2.1**
设P(x;{0,?,k?1}):=1。给定一个N2k中的单坐标增量序列{x0,?,xn}，算法1在Θ(k2k(n+∑i=02k?1xi0))次操作内计算{P(xn′;{j}):j∈{0,1,?,k?1},n′∈{0,1,?,n}}（包括评估Beta函数）。这些正是运行一个n患者S-BRAR试验所需的概率，前提是先验为x0，状态路径等于x0,x1,…,xn。

**下载：下载高分辨率图像（133KB）**
**下载：下载全尺寸图像**

算法1用于精确计算具有整数Beta先验的BRAR试验中的所有后验概率。证明见附录A.2。该计算自然适用于批量随机化和初始均衡分配阶段。

通过一个简单的例子可以理解算法的工作原理。考虑一个具有均匀Beta(1, 1)先验的双臂试验。初始状态可以表示为向量x0=(α0,β0,α1,β1)=(1,1,1,1)。此时，任一臂更优的概率正好是0.5。假设前三个观察到的患者结果是：患者1被分配到臂0并成功，患者2被分配到臂1并失败，患者3被分配到臂0并失败。那么，后续状态只需将相关坐标递增1：
(1,1,1,1)→P(x;{1})=1/2→(2,1,1,1)→P(x;{1})=1/3→(2,1,1,2)→P(x;{1})=1/6→(2,2,1,2)→P(x;{1})=3/10。

在每一步中，算法1使用先前状态存储的后验概率进行更新，而不是从头开始重新计算。例如，在患者1在臂0上成功后，我们可以使用先验状态计算臂1更优的更新概率：
P((2,1,1,1);{1})?P((1,1,1,1);{1})=?b0,1((1,1,1,1))x0P((1,1,1);{0,1})，
因为k=2，P(x;{0,1})=1。项b0,1来自Beta函数，对于初始状态等于1/6，且x0=1。因此，新的概率只需基本算术：P((2,1,1);{1})=0.5?(1/6)。因此，臂0的更优概率为2/3。重复这种一坐标更新可以高效地得到整个试验的概率。

为了计算孤立的后验概率P(x; {j})，我们可以将问题视为找到具有均匀先验的S-BRAR试验的最终状态。通过将目标参数视为从全1向量开始的试验路径的终点，应用定理2.1得到：
**推论2.1**
P(x; {j})可以在Θ(k2k(∑i=02k?1xi?2k))次操作内精确计算，对于任何k≥2且x∈N2k。

值得注意的是，对于k=2，这导致运算复杂度为Θ(x0+x1+x2+x3)。相比之下，方程(2)的改编实现的复杂度为Θ(min(x0, x1, x2, x3)。这些结果的递归结构基于以下引理，该引理扩展了Cook（2003年）的方法：
**引理2.1**
如果0≤j,j′≤k?1, k≥2, s∈{0, 1}, x∈N2k，则
(3) P(x+ej,s;{j′})?P(x;{j′})={
(?1)s∑j″≠jbj,j″(x)x2j+sP(x;{j,j″})
如果j=j′，
(?1)s+1bj,j′(x)x2j+sP(x;{j,j′})
如果j≠j′，
其中el,j,s:=1(l=2j+s)表示递增1，bj,j′(x):=B(x2j+x2j′,x2j+1+x2j′+1)B(x2j,x2j+1)是为了符号方便引入的量。

引理2.1的证明见附录A.3。这个看似技术性的结果可以通过类比b=1的S-BRAR试验来直观理解。假设一个患者根据已知的当前概率{P(x;{j′}):j′∈{0,?,k?1}}被随机分配到治疗j，并表现出反应s。然后根据更新后的概率{P(x+ej,s;{j′}):j′∈{0,?,k?1}}对下一个患者进行随机分配。引理2.1表明，计算这些更新后的概率相当于计算集合{P(x; {j, j′}): j′≠j}，这实际上代表了k?1臂试验的后验概率。虽然可以通过递归减少到双臂情况（可以显式评估），但算法1通过存储所有非空S?{0,?,k?1}的P(x; S)采用了更高效的方法。因为(3)自然适用于任何非空子集S（而不仅仅是单元素集合S={j′}），我们可以使用这些先前存储的值更新整个概率集{P(x+ej,s;S)}，从而大幅减少计算负担，仅需Θ(k2k)次操作即可。

虽然我们专注于整数先验，但算法1也适用于连续Beta参数。仅需要用对第一个随机化概率的准确前期评估来替换第1-10行，这对于小的先验值来说是容易实现的（Cook，2003）。然而，研究人员应注意，非整数先验可能会强烈影响抽样过程，因此在临床环境中可能不如均匀先验受欢迎。2.2.2. 近似方法另一种方法是使用以下方法之一来近似P(x; {j})：•数值积分（NI），例如在ARREST试验中使用的（Yannopoulos等人，2020），它利用了由于各组独立性而产生的后验概率的一维积分表示。设αj:=aj+Si,j和βj:=bj+Ni,j?Si,j，并分别用fα,β和Fα,β表示Beta(α, β)的密度和分布函数，我们有(4)P(x;{j})=∫01fαj,βj(u)∏j′≠jFαj′,βj′(u)du。我们使用m节点高斯-勒让德规则在[0, 1]上近似这个积分，通过log?B(αj, βj)来评估fαj,βj(u)，并通过正则化不完全贝塔函数来评估每个Fαj′,βj′(u)，并将近似量表示为PNI(x; {j})。m的选择是从m=20开始，每次增加20，直到评估积分的变化小于10?6。•重复抽样（RS）最初由Ulam开发（Metropolis和Ulam，1949），并应用于如ADORE（Carlson等人，2017）这样的试验，这是一种直观的蒙特卡洛方法。它涉及从联合后验P:=?j′=0k?1Beta(x2j′,x2j′+1)中抽取K个独立样本X1,…,XK，并计算其中组j是最优组的样本比例，即(5)P^RS(x;{j}):=1K∑K′=1K1(Xj,K′>maxj′≠jXj′,K′)。这种方法需要Θ(Kk)次操作。•高斯近似（GA）由Cook（2012）提出，并用于Marian等人（2024）的贝叶斯重新设计中，它用均值和方差匹配的正态分布来近似Beta分布。对于一般情况k?≥?2，这产生：(6)PGA(x;{j})=Fj(0,…,0)，其中Fj是多变量（k?1维）正态分布N(μ(j),Σ(j))的累积分布函数，μ(j)=μ?j?mj1k?1，Σ(j)=D?j+sj21k?11k?1?，μ?j是移除了mj的向量[m0,?,mk?1]，1k?1是长度为k?1的全1向量，D?j是对角线为[s02,?sk?12]的对角矩阵，且移除了sj2，对于所有j′，都有mj′=x2j′x2j′+x2j′+1,sj′2=x2j′x2j′+1(x2j′+x2j′+1)2(x2j′+x2j′+1+1)。Altham（Altham，1969）提出了另一种2组GA方法，但它不容易推广到k?>?2的情况，据我们所知，在实践中很少使用。因此，我们的一般GA的计算成本为Θ(CMVN(k?1))，其中CMVN(k)是k维正态CDF评估的成本。这个成本会根据分配的患者及其结果产生的协方差结构而有很大差异。2.3. 后验概率的双重作用：推断和随机化首先，我们评估单独近似单一概率P(x; {j})的准确性和计算复杂性之间的基本权衡。在此基础上，我们在试验中后验概率发挥的两种主要作用——推断和随机化方面对这两种近似方法进行基准测试。当用作检验统计量时，我们量化近似误差如何影响关键推断属性，如I型错误率和功效。相反，当这些概率用于指导整个试验中的患者随机化时，我们研究近似误差如何随时间累积并影响患者分配以及其他关键操作特性。由于概率在测试和随机化中的使用方式不同，每种任务可能需要不同的计算方法。为了实施贝叶斯拒绝规则进行假设检验，必须明确评估后验概率。然而，对于标准随机化，明确定义πS-BRARQ(Hi)的值并不是严格必要的，可以直接从后验Qi*中抽样，并将患者i+1分配给产生最大样本的组（Russo等人，2018）。尽管如此，许多正则化的S-BRAR版本仍然需要显式的概率。调整程序用于缩放随机化概率，从而减轻极端或高度可变分配的影响。我们通过应用调整函数t:[0,1]→R+k来定义一个调整后的BRAR（T-BRAR）程序：(7)πT-BRAR,jQ(Hi):=tj(πS-BRAR,jQ(Hi))∑j′=0k?1tj′(πS-BRAR,j′Q(Hi))。Bleck等人（2013）使用的一种标准调整方法不仅用于优先考虑预期治疗效果较高的组，还优先考虑方差较大的组。考虑方差对于防止过早锁定和极端分配不平衡至关重要。对于参数m∈N，这种方差缩放方法定义为：(8)tj,VS(m)(πS-BRAR,jQ(Hi)):=(πS-BRAR,jQ(Hi)VarQi*+1(pj)Ni,j+1)1m。2.4. 假设检验和误差控制在设计具有拒绝规则R(P,c)的试验时，阈值c通常被选择为将I型错误率控制在名义显著性水平α以下。假设零假设可以简化为全局零假设p0=…=pk?1，一个无条件的精确检验UX(P,α)通过将阈值设置为c(P,α)来实现这一点，c(P,α)被定义为保持最大I型错误率低于α的最小c。确定这个精确阈值在计算上非常耗费资源，这是高效概率计算的另一个关键应用。这通常首先通过定义校准后的后验概率检验PP(P,p,α)来进行。该检验找到在零假设下控制特定响应概率p的错误率的最小c(P,p,α)。因为c(P,α)=supp∈[0,1]c(P,p,α)，所以可以通过最大化这些逐点阈值来推导出精确的无条件阈值（参见Baas等人（2025b）中的算法2）。计算c(P,p,α)需要评估所有可能试验状态x下的所有组j的后验概率P(x; {j})，直到最大样本量n（算法1，Baas等人（2025b））。对于均匀先验，这需要计算Θ(n2k?1)个不同的概率。因此，随着n或k的增长，这些评估的计算效率变得至关重要。最终，准确的概率计算不仅支撑患者分配和检验统计量的计算，还支撑精确检验本身的推导。如果一致使用如GA这样的确定性近似来计算序列检验统计量和阈值c(P,p,α)，I型错误率将保持在指定的点p。然而，近似仍然可能无意中改变更广泛的I型错误分布和试验的统计功效。2.5. 操作特性我们考虑了一系列操作特性（OCs）。设Tˉ为试验停止的患者编号，根据Baas等人（2025a）的定义：•I型错误率是在零假设为真时我们拒绝零假设的概率，即Pπ(maxj∈{1,?,k}Tj(HTˉ)>c)。•功效是我们正确拒绝零假设并支持HAj的概率，即Pπ(Tj(HTˉ)>c)，当j是更优的治疗组时。•当j是更优的治疗组时，分配到更优组的预期患者数量（EPASA）由Eπ[Nj,Tˉ+(n?Tˉ)1(Tj(HTˉ)>c)给出。这考虑了在可能提前停止之前和之后的患者分配。这确保了为了支持最佳治疗而提前停止不会受到惩罚。•当j是更优的治疗组时，分配到更优组的患者方差（VPASA）由Vπ(Nj,Tˉ+(n?Tˉ)1(Tj(HTˉ)>c)给出。这旨在捕捉患者收益的变异性；否则，高的EPASA可能会掩盖将许多患者分配到较差组的风险。因为BRAR试验已被证明会产生有偏的治疗效果估计（Bowden和Trippa，2017），它们更适合区分感兴趣的检验假设，而不是提供精确的效果估计。注意，所有这些OC都可以表示为HTˉ的函数的期望值（VPASA通过标准方差公式得出）。如果我们进一步指定拒绝规则的形式为R(P,c)，这些OC可以使用Baas等人（2025b）中概述的框架来精确计算。如果使用如GA这样的确定性近似来随机化患者或近似检验统计量，也可以进行这种精确计算。然而，如果近似是随机的（如RS），或者患者和治疗组的状态空间太大，则精确计算变得不可行。在这种情况下，我们可以通过模拟来估计OCs。具体来说，我们运行大量K次试验复制，并取平均结果来估计目标OC。虽然我们的基线分析主要使用K=104个样本来平衡计算可行性和I型错误膨胀，但实际上试验有时会使用更大的样本量，例如Carlson等人（2017）使用的106个样本。为了全面描述这种准确性-速度权衡，我们还系统地研究了K的变化（从104到106）如何影响计算时间和结果的推断操作特性。3. 计算成本本节分析了使用精确和近似方法评估单个优效后验概率P(x; {j})（对于某些j?≥?0和x∈N2k）的计算成本，然后检查这在BRAR试验过程中的变化。3.1. 评估单个优效后验概率在第2节中，我们描述了不同方法进行PPS评估所需的预期理论操作规模。然而，对于GA，这涉及一个黑箱函数CMVN(k)，我们预计在组分配不均匀时它会更小。因此，我们在两种设置下对计算成本进行基准测试：一种是“平衡”设置，其中各组的成功、失败和分配是相等的；另一种是“典型”设置，在S-BRAR试验中评估第200个患者的试验状态的中位数，b=1（假设一个组的响应概率为0.6，其他组为0.5，并且分配概率每20个患者更新一次）。同样，对于k=2的精确计算，方程7在经过对称调整后比算法1计算更快，但在“平衡”设置中也比“典型”设置成本更高，因此需要这种区分。因此，图1(a)绘制了k=2时10次单PPS计算的平均时间（在600次运行中的中位数），并以对数尺度表示。它还包括最小二乘预测：精确方法使用n的线性模型（可能具有非零截距以捕捉小n时的固定效应），而近似方法使用常数模型。固定n=200，图1(b)绘制了k变化时的平均时间（在60次运行中的中位数），以及第2节中概述的缩放下的最小二乘预测（截距设置为0，除了RS，需要截距来捕捉实证行为）。因为GA根据k表现出特定的校准，其最小二乘拟合仅针对k?≥?8进行计算，并进行了外推。预测和数据之间的强烈吻合证实了理论缩放律在实践中的有效性。下载：下载高分辨率图像（332KB）下载：下载全尺寸图像图1. 使用不同方法的10次单PPS计算的平均时间（在600次运行中的中位数）以及它们的预测，其中(a) k=2且n变化；(b) n=200且k变化。对于k=2和变化的n，NI的经验计算成本在找到适当的m之后进行评估，与算法1相似，尽管对于较大的n速度稍快，但对于较小的n速度较慢，因此在这种情况下它是最慢的近似方法，同时也没有下一小节讨论的精确计算的顺序效率优势。然而，对于n=200和变化的k，选择的m并不随k增加而增加，因此它的成本与RS在K=104时相似；RS在k?≤?6时稍快，而对于更大的k，NI更快。然而，我们在本文的其余部分将重点放在RS上，有两个原因。首先也是最重要的是，NI并不直接承认n的理论缩放律，不像所有其他方法那样，而且像GA和k那样对每个n进行基准测试是不可行的，因为可能的n范围要大得多，这限制了在没有运行完整试验模拟的情况下计算预期计算成本的可能性。其次，这里评估的特定NI方法只是许多潜在的求积选项中的一种，没有任何一种在临床试验实施中成为标准。对于k=2和n?≤?400，通过方程（2）进行的精确计算是“典型”场景中最快的方法，而GA对于较大的n更快。当n固定而k变化时，算法1和GA都表现出接近指数的缩放；值得注意的是，在“平衡”场景中，GA虽然较慢但更稳定。至关重要的是，虽然使用算法1进行精确计算通常对于单个PPS效率低下（即使对于较小的k也是如此），但这种静态分析并没有反映出BRAR试验的顺序性质，其中多个PPS是依次评估的。正如我们在下一节中所示，正是在顺序更新的动态背景下，算法1展现了强大的比较优势。3.2.### 评估整个试验中优势的后验概率

在典型的BRAR试验中，初始的燃烧期通过均匀随机化将B名患者分配到每个组中，然后每b名患者更新一次概率。假设调整程序的计算成本可以忽略不计，那么S-BRAR和T-BRAR的预测总成本分别为：
- \(T_{\text{Exact}} = c_{A1} \cdot n \cdot k \cdot 2^k\)
- \(T_{\text{TRS}} = c_{RS} \cdot K \cdot k \cdot U(n,k,B,b)\)
- \(T_{\text{GA}} = k \cdot c_{GA} \cdot C_{MVN}(k-1) \cdot U(n,k,B,b)\)

其中\(c_{A1}\)、\(c_{RS}\)、\(c_{GA}\)是预先校准的常数，\(U(n,k,B,b):=\lceil n - kBb \rceil\)表示燃烧期后进行的中间分析次数。关键的是，算法1能够一次性处理所有患者的数据序列。它维护所有\(2^k\)个子集概率的运行状态，自然地能够在试验的每一步产生所需的后验概率（忽略\(n\)成本中的截距项）。因此，其总计算成本与中间分析的频率无关。相比之下，RS和GA必须在每次中间分析时从头开始重新计算后验概率，使得它们的总成本与\(U\)成正比。注意，在每次中间分析时都需要所有的\(k\)个后验概率。对于RS，从\(k\)个Beta分布中抽取一批样本可以同时得到所有\(k\)个值（通过计算不同的argmax结果），因此其每次中间分析的成本保持在\(\Theta(Kk)\)。然而，对于GA，每个组都需要单独进行\(k-1\)维的MVN累积分布函数(CDF)评估，这导致了(11)式中的额外因子\(k\)。为了说明这种精确计算的递归优势，当\(k=2\)时，它的速度可以是RS的257倍（随着\(K\)的增加这个优势线性增长），并且是GA的60倍。

这些成本模型定义了两个计算阈值。由于RS和GA的每次中间分析成本都有一个共同的因子\(k\)，它们的切换仅取决于\(K\)：
- \(k_1(K) = \min\{k \geq 2 : c_{RS}K \leq c_{GA} \cdot C_{MVN}(k-1)\}
- 当\(k < k_1\)时，GA是更快的近似方法；当\(k \geq k_1\)时，RS是更快的方法。

由于其精确性，除非算法1明显较慢，否则优先选择它。通过引入一个容忍度\(r_{\text{mix}} > 1\)，当\(T_{\text{Exact}} \leq r_{\text{mix}} \cdot \min(T_{\text{TRS}}, T_{\text{GA}}\)时，推荐使用精确计算。对于\(k \geq k_1\)，忽略进一步的目标精度考虑，只要满足\(c_{A1} \cdot n \cdot 2^k \leq r_{\text{mix}} \cdot c_{RS} \cdot K \cdot U(n,k,B,b)\)，就应该优先选择精确计算。

### 3.3 实现、稳定性和内存扩展

本文中的所有基准测试都是在14英寸MacBook Pro（M3，18 GB RAM，macOS Sonoma 14.1）上使用Python 3.11在VSCode中执行的。实现依赖于NumPy进行向量化数组操作，以及SciPy进行特殊函数（如\(scipy.special.betaln\)用于对数Beta值评估，\(scipy.special.betainc\)用于规则化不完整Beta函数）和多变量正态CDF评估（\(scipy.stats.multivariate_normal\)）。没有使用编译扩展或GPU加速。为了确保可重复的单线程基准测试，所有的BLAS和OpenMP线程数都限制为一个。

虽然绝对的墙钟时间和特定的切换点（例如，当\(k=2\)时精确计算优于GA，对于\(n \leq 400\)）是硬件特定的，但理论上的操作复杂性——精确计算的\(\Theta(nk^2k)\)、RS的每次评估\(\Theta(Kk)\)和GA的每次评估\(\Theta(C_{MVN}(k-1)\)是完全与硬件无关的。因为所有方法都使用相同的编译数值库，硬件升级应该会按比例加速它们，保持常数比率\(c_{A1}/c_{RS}\)和\(c_{A1}/c_{GA}\)。重要的是，我们在k上的指数与多项式差异推荐仍然对常数因子执行速度的变化具有鲁棒性。

算法1需要存储2^k?1个浮点值，用于子集概率\(\{P(x; S): \emptyset \neq S \subseteq \{0,\ldots, k-1\}\)，以及相同大小的辅助数组，用于累积参数之和和缓存的log-Beta值，总内存占用为\(\Theta(2^k)\)。实际上，每个双精度浮点数占用8字节，这意味着当\(k=10\)时大约为32 KB，\(k=15\)时为1 MB，\(k=20\)时为32 MB。相比之下，RS需要\(\Theta(Kk)\)的内存来存储样本矩阵（例如，当\(K=10^5\)和\(k=15\)时大约为15 MB），而GA只需要\(\Theta(k^2)\)的内存来存储协方差矩阵。因此，精确方法的内存限制不如其计算成本严格：当\(k=25\)时，大约1 GB的内存需求仍然是可管理的，而\(\Theta(nk^2k)\)计算时间对于合理的\(n\)来说可能是不可行的。

### 4. 准确性：研究双臂S-BRAR

在评估计算复杂性后，我们评估准确性以确定近似方法何时会失败，并突出未验证近似方法的风险。为了特别清楚地展示高度顺序（\(b=1\)）、零燃烧期（\(B=0\)设计的近似方法的陷阱，我们分析了一个双臂S-BRAR试验，假设先验是均匀的，并且不提前停止。

#### 4.1 评估近似单个后验概率的效果

首先，我们评估近似单个后验概率的准确性。用\(P(a, b, c, d)\)表示\(P((a+1,b+1,c+1,d+1); \{1\})\)。作为独立同分布(Bernoulli)随机变量的均值，有\(10^4P^RS \sim Bin(10^4,P)\)。使用Blyth (1980)的公式，我们发现：
\[E[|P^RS - P|] = 2(10^4 - 1)^{l-1} P^l(1 - P)^{10^4 - l + 1}\]
其中\(l\)满足\(l-1 \leq P \times 10^4 < l\)。这个函数在\(P=0.5\)时达到最大值3.99×10^?3。根据经验，对于任何固定的\(n \leq 200\)的分配，最大均值绝对误差\(E[|P^RS(a,b,c,d) - P(a,b,c,d)|\)大约为3.99×10^?3。这提供了一个普遍的误差界限，适用于任何\(j\)、\(k\)和\(x\)，使我们能够根据需要选择\(K\)：
\[E[|P^RS(x; \{j\}) - P(x; \{j\})|] \leq (K-1)^{\lceil K/2 \rceil}^2 - K\]

相比之下，GA产生的误差更大，提供的准确性保证非常弱。对于\(n \geq 113\)，GA的绝对误差可以超过0.1。反直觉的是，如果成功/失败计数保持较小，较大的\(n\)可能会导致更大的误差。此外，RS的误差在极端分配时达到最大（如果一个组没有患者，则超过0.003），这对于像S-BRAR这样的RA程序来说是令人担忧的。然而，即使在平均情况下，GA的误差也比使用\(K=10^4\)的RS大一个数量级（见图B.3）。

#### 4.2 评估在整个试验中近似后验概率的效果

我们区分了近似随机化概率和测试统计量的影响。尽管试验通常会近似两者，但可以混合使用这些方法来优化速度-准确性权衡。对于较小的\(k\)（\(k=2\)，GA比RS快得多；因此，将GA用于随机化，RS用于测试可以优化速度-准确性权衡。然而，随着设计中中间分析的频率增加，这种优势会减弱。

首先，考虑试验内患者收益的操作特性(OCs)。仅使用GA进行随机化对预期患者收益(EPASA)的影响有限，变化范围从-4.17%到+0.76%（绝对变化从-8.00到1.49）。这种影响在大部分参数空间中是正的，在响应概率不同（例如，\(p_0, p_1 \approx 0.03, 0.4\)）时达到峰值，只有在两者都接近1时才达到最小值。然而，GA引入了明显更高的可变性。方差(VPASA)在几乎整个参数空间中增加——最多增加55.0%——只有在极端概率接近1时才会减少（最多减少26.1%）。相应地，分配给优效治疗的比率的标准差在-0.0407到0.0107之间变化。总体而言，用GA代替精确计算进行随机化会导致平均患者收益略有增加，但代价是方差增大。相反，RS用于随机化不会改变OCs。因为没有调整程序，从每个Beta分布中抽样并选择样本量最大的组是完全精确的。同样，仅近似测试统计量概率也不会影响试验内患者收益的OCs。

关于I型错误和功效，用RS近似随机化没有影响，无论\(K\)是多少。然而，对于GA，影响是显著的（见图B.4）。假设\(n=200\)，\(P=\{n}\)，并且在本节的其余部分\(\alpha=0.05\)，GA根据测试和响应概率\(p\)导致I型错误膨胀或缩减。最大I型错误膨胀为PP(0.6)的0.0127，UX的0.00369。值得注意的是，在GA下，UX的I型错误率严格保持在0.05以下，这意味着GA没有损害其错误控制。然而，GA在大部分参数空间中都会导致功效损失。这种损失可能很严重，对于PP(0.6)在某些\(p_0, p_1\)组合下超过0.42。总体而言，虽然GA对PP(0.6)和UX的I型错误膨胀最小，但它以显著降低功效为代价，仅提供微小的功效增益。

如果只近似测试统计量\(T_1=1 - T_0\)，对PP(0.6)和UX的I型错误和功效的影响比近似随机化要严重得多（见图2）。如图2(a)所示，使用GA显著改变了I型错误率的定性行为。对于这两种测试，大的\(p\)值会导致严重的误差减少，而适度的\(p\)值会导致较大的误差膨胀（PP(0.6)达到0.048，UX达到0.031）。相反，使用RS（\(K=10^4\)进行测试统计量只会导致与精确I型错误率的最小偏差（见图2(b)），尽管这种差异在\(p\)接近1时最大。

### 下载：
- 下载高分辨率图像（634KB）
- 下载全尺寸图像

**图2.** 对于\(n=200\)和不同的\(p_0\)和\(p_1\)值：(a) 如果使用精确概率和GA进行测试，PP(0.6)和UX的I型错误率；(b) 如果使用RS或精确概率进行测试，I型错误率的估计差异；(c) 对于PP(0.6)的测试；(d) 对于UX的测试；(e) 对于PP(0.6)的测试；(f) 对于UX的测试，使用GA进行测试时的功效变化。

使用GA进行推断会显著改变功效（见图2(c)–(d)）。对于这两种测试，功效通常在参数空间中增加，但当任一响应概率较大时急剧下降。在PP(0.6)下，功效最多增加0.65（接近\(p_0, p_1 \approx 0, 0.2\)），而在\(p_0, p_1 \approx 1, 0.9\)时减少0.75。UX测试的波动更大：最大增加0.81（接近0, 0.4），最大减少0.84（接近1, 0.8）。有趣的是，当使用GA进行随机化时提高功效的响应概率，在使用GA进行测试时往往会导致功效降低。

相反，使用RS进行测试会导致较小的但非微不足道的功效变化（见图2(e)–(f)）。一般来说，较小的响应概率会带来功效增益，而较大的响应概率会导致功效损失。PP(0.6)的最大增益约为0.02，最大损失为0.008，而UX则呈现相反的趋势。即使考虑到估计值的95%置信区间（±0.0043），RS对功效的影响也明显小于GA。为了量化蒙特卡洛预算\(K\)对RS测试的影响，我们在相同的100×100网格上计算了最大绝对扭曲。对于PP(0.6)，最大扭曲分别为K=103、104、105和106时的0.0372、0.0162、0.0122和0.00664。UX的值类似：0.0660、0.0190、0.0123和0.00645。因此，增加K可以严格减少最坏情况下的扭曲，尽管没有简单的通用缩放规律。最终，如果保持原始试验设计的功效至关重要，RS在计算测试统计量方面明显优于GA，其可管理的计算成本是合理的。因为微小的测试统计量错误可以直接改变最终的试验决策，即使RS也会引入一些在近似随机化时看不到的偏差。虽然GA的速度通常伴随着严重的推断扭曲，但UX测试是一个显著的例外：GA保持了I型错误控制，并提高了小响应概率的功效，因此在\(p_0\)和\(p_1\)不太可能较大的情况下是可取的。因此，近似方法的适用性在很大程度上取决于特定的贝叶斯测试统计量和试验的潜在目标，以及我们预期在响应概率空间中的位置。

### 5. 案例研究：ESET试验

实际试验很少像之前讨论的设置那样简单。尽管如此，这些更简单的模型为研究人员如何系统地评估和选择控制患者分配和统计测试的概率的计算方法提供了一个有价值的框架。例如，在实践中，到达终点观察的时间可能有很大差异，新患者的累积速度通常比这些快，因此在分配后续患者之前观察所有响应（\(b=1\)通常是不现实的。在本节中，我们还考虑了更大的\(b\)范围，但请注意，新块中第一个患者的分配与确定此分配所需的中间分析之间的延迟是BRAR方法的一般性后勤限制。

基于前一节的逻辑，我们推导出针对现实场景的实际建议，现在我们分析了一个具有多个组、块随机化、初始等分配燃烧期、提前停止和调整的试验。已建立的癫痫治疗(ESET)试验(Bleck等人，2013)比较了三种活性治疗：磷苯妥英、左乙拉西坦和丙戊酸（分别标记为0、1和2）。假设有720名可评估的患者，原始设计将前300名患者平均分配到所有组中（每个组100名的燃烧期）。剩余的患者使用T-BRAR进行分配，方差缩放（\(m=2\)，每100名患者更新一次随机化概率）。

该试验旨在确定唯一最佳和/或最差的治疗。原假设和备择假设为：
H0: p0 = p1 = p2;
Hj?: pj > pj′ 对于所有 j′ ≠ j；
Hj?: pj < pj′ 对于所有 j′ ≠ j；
Hj?, j″: pj > pj″ > pj 对于所有 j″ ≠ j, j′。

我们定义治疗 j 最差的检验统计量为 Tj′(Hi) = ∑Qi(pj < minj′ ≠ pj′)。这可以与计算 Tj 的方法相同，因为 Qi(pj < minj′ ≠ pj′) = Qi(1 ? pj) > maxj′ ≠ j)。

在试验结束时，如果对于剩余的不同处理组 j (最佳) 和 j′ (最差)，有 Hj?, j″ 且 Tj(Hn) > 0.975 以及 Tj′′(Hn) > 0.975，则拒绝 H0。否则，如果 Tj(Hn) > 0.975，则拒绝 Hj；如果 Tj′(Hn) > 0.975，则拒绝 Hj。

在每 100 个患者完成一个治疗组后进行中期分析（从 N=400 开始）。如果 Tj(Hi) > 0.975，则提前终止试验。此外，如果 Qi(pj < 0.25) ≥ 0.95，则放弃该治疗组；如果所有治疗组都被放弃，则完全终止试验。（注：为了严格分离 PPS 计算的计算成本，这里省略了原始设计中的额外预测性无效规则。）

我们评估了每个治疗组的 burn-in 长度 (B) 和中期分析频率 (b) 的影响。原始配置是 B = b = 100。我们还考虑了一个扩展的 k=5 治疗组设计。因为在随机化每个新治疗组之前已经计算了检验统计量，所以随机化概率本质上是重新使用的检验统计量。因此，仅为了随机化而使用更快但不更准确的方法并没有计算上的优势。

表 1 详细说明了在假设没有提前终止的情况下模拟 100 次试验所需的最大时间。主要发现包括：
- 减少 B 或 b 对精确计算时间的影响很小，因此精确计算显示出很高的稳定性。
- 随着 B 和 b 的减少，遗传算法 (GA) 和 Robbins-Sampson (RS) 的速度大幅降低。对于高度顺序化的试验 (b=1, k=3)，GA 和 RS 的时间长达 b=100 时的 61 倍。因此，在顺序化设计中，近似方法会受到影响。
- GA 显示出重大的计算瓶颈：对于 k=5 和完全顺序化的设置 (b=1)，准确的操作特征 (OC) 计算变得几乎不可能。在 B=100 时，通过 GA 评估第一类错误率大约需要 76 小时。因此，排除了 B=0 或 b=1 且 k=5 的配置，不进行进一步的 OC 分析。

表 1. 计算 1000 次 ESET 试验所需的最大时间（以秒为单位），考虑了不同的后验概率计算方法、burn-in 周期长度 B、中期分析频率 b 和治疗组数量 k。
| 空单元 | k=3 | k=5 |
| ---- | ---- | ---- |
| b=100 | 33.44 | 33.42 |
| 35.50 | 45.33 | 144.21 |
| 45.11 | 44.01 | 49.55 |
| 34.33 | 35.41 | 38.84 |
| 49.39 | 144.01 | 41.31 |
| 44.51 | 55.50 | 33.81 |
| 32.90 | 35.79 | 49.04 |
| 139.91 | 41.51 | 38.21 |
| 49.8 | | |

表 2 报告了基于 105 次模拟的第一类错误率和功效。使用 McNemar 的配对检验基于共同随机数来评估拒绝率之间的差异。

• 第一类错误率 (p0 = p1 = p2 = 0.5)：与精确计算相比，GA 在每个评估的配置中都表现出统计上显著的第一类错误率膨胀（p < 10^-6）。虽然在 B=100 时膨胀幅度较小（0.13–0.34 个百分点），但随着 B 和 b 的减小，这种膨胀会急剧增加。在 B=0, b=1 (k=3) 时，GA 的错误率为 18.9%，而精确计算的错误率为 13.7%——相对增加了 38%。相反，RS (K=104) 在 18 个配置中成功控制了第一类错误率，与精确计算相比没有统计上显著的膨胀。
• 效效：在最佳治疗组的情景下，所有方法的功效几乎相同；GA 和 RS 的估计值都在精确计算的模拟不确定性范围内。在最差治疗组的情景下，GA 在极端配置下显示出轻微的、非系统性的功效波动（最高达 1.6 个百分点），而 RS 则精确地反映了精确计算的结果。

表 2. 基于 105 次模拟的 ESET 试验的第一类错误率和功效，考虑了不同的随机化和检验计算方法、burn-in 周期长度 B 和中期分析频率 b。还报告了 95% Clopper–Pearson 精确置信区间的半半径。

表 3 提供了使用不同后验概率计算方法、burn-in 周期长度 B、中期分析频率 b 和治疗组数量 k 进行 1000 次 ESET 试验所需的最大时间（以秒为单位）。
| 空单元 | k=3 | k=5 |
| ---- | ---- |
| b=100 | 33.44 | 33.42 |
| 35.50 | 45.33 | 144.21 |
| 45.11 | 44.01 | 49.55 |
| 34.33 | 35.41 | 38.84 |
| 49.39 | 144.01 | 41.31 |
| 44.51 | 55.50 | 33.81 |
| 32.90 | 35.79 | 49.04 |
| 139.91 | 41.51 | 38.21 |
| 49.8 | | |

表 2 报告了基于 105 次模拟的第一类错误率和功效。差异使用 McNemar 的配对检验基于共同随机数进行评估。