r - 贝叶斯自适应随机化：如何计算参加治疗组的平均人数

Question

这是我正在处理的问题。这是贝叶斯自适应随机化设计。符合条件的患者已被随机分配至标准治疗或新的靶向药物。试验的主要终点是 4 周的临床反应状态。设计一个 II 期试验，假设两个手臂：

ARM 1：患者接受标准化疗，已知有效率 p1 = 0.2（成功的后验概率 = thetaA）

ARM 2：给予患者新的新辅助化疗，这应该使反应率翻倍 p2 = 0.4（成功的后验概率 = thetaB）

前 40 名患者同样随机分配到新组 (20 + 20)，然后根据成功的后验概率自适应随机分配。

为了让可能性支配后验分布，对 p1 和 p2 采用无信息分布（最简单的情况是在同一患者身上测试两种治疗方法，因此这意味着要考虑成功/失败的先验分布概率等于 0.5 + 0.5 = 1)。我被要求进行 100 次模拟，假设参数 c 为最优，计算 Nmax = 100（最大患者数）的随机概率，并最终认为：

the mean number of patients enrolled in ARM 1 and ARM 2

我还假设 0.01 和 0.98 分别作为无效和有效性提前停止的阈值，如果所有 100 名患者都经过评估，则宣布有效的阈值为 0.9。

因此，我为模拟提出的代码如下（设置参数，如您在#SET PARAMETERS 处所读到的）

AR = function(priorA, priorB, c=1){
  posteriorA = rbeta(10000, priorA[1], priorA[2])
  posteriorB = rbeta(10000, priorB[1], priorB[2])
  pab = posteriorA-posteriorB
  
  pA.less.B = length(pab[pab<0])/length(pab)
  pA.great.B = length(pab[pab>0])/length(pab)
  pB.great.A = length(pab[pab<0])/length(pab)
  
  rb = (pA.less.B^c)/(pA.less.B^c + pA.great.B^c)
  ra = 1- rb
  ar=list(ra=ra, rb=rb, pA.less.B=pA.less.B, 
          pA.great.B = pA.great.B, pB.great.A = pB.great.A)
  ar
}

update.prior = function(priorA, patientA, priorB, patientB){
  priorA=priorA+patientA
  priorB=priorB+patientB
  prior=list(priorA=priorA, priorB=priorB)
  prior
  
}

#SET PARAMETERS
set.seed(123)
thetaA = 0.2
thetaB = 0.4
priorA = c(0.5,0.5) 
priorB = c(0.5,0.5) 
inits = 20
c = 'optim'
futility = 0.01
efficacy = 0.98
nsims = 100
Nmax = 100

Na=NULL
Nb=NULL
result = NULL
Rb = NULL


for(j in 1:nsims){
  Nmax=Nmax
  # yA = rbinom(n=Nmax, size=1, prob=0.25)
  # yB = rbinom(n=Nmax, size=1, prob=0.25)
  
  priorA = c(0.5,0.5) 
  priorB = c(0.5,0.5) 
  treatment=NULL
  posterior=NULL
  success = NULL
  
  for (k in 1:inits){
    patientA = rbinom(1, size=1, prob=thetaA)
    patientB = rbinom(1, size=1, prob=thetaB)
    prior=update.prior(priorA, c(patientA,1-patientA), priorB, c(patientB,1-patientB))
    priorA=prior$priorA
    priorB=prior$priorB
    success=c(success,patientA,patientB)
  }
  
  A=inits
  B=inits
  treatment[1:(A+B)]=rep(c('A','B'), inits)
  
  if(c=='optim'){
    c=(A+B)/(2*Nmax)}else{c=c}
  ar= AR(priorA, priorB, c=c)
  
  #check for futility
  if (ar$pA.less.B<=futility) {print('stop for futility') 
    
    result[j] = 'stop for futility'
    Rb[j] = 'stop for futility'}
  
  #check for efficacy
  if (ar$pA.less.B>=efficacy) {print('stop for efficacy')
    result[j] = 'stop for efficacy'
    Rb[j] = 'stop for efficacy'}
  
  ra=ar$ra
  rb=ar$rb
  result[j] = ar$pA.less.B
  Rb[j] = ar$rb
  
  i=A+B
  while ((A+B)<Nmax & ar$pA.less.B<efficacy & ar$pA.less.B>futility ){
    i=i+1
    randomize = sample(c('A','B'), 1, prob = c(ra, rb))
    if (randomize=='A'){
      A=A+1
      patientA = rbinom(1, size=1, prob=thetaA)
      success[i] = patientA
      patientA = c(patientA, 1- patientA)
      patientB = c(0,0)
      treatment[i]='A'
      # success[i]=yA[A]
    }
    
    if (randomize=='B'){
      B=B+1
      patientB = rbinom(1, size=1, prob=thetaB)
      success[i] = patientB
      patientB = c(patientB, 1-patientB)
      patientA = c(0,0)
      treatment[i]='B'
    }
    
    if(c=='optim'){
      c=(A+B)/(2*Nmax)}else{c=c}
    prior=update.prior(priorA, patientA, priorB, patientB)
    priorA=prior$priorA
    priorB=prior$priorB
    ar=AR(priorA, priorB, c=c)
    #print(ar)
    ra=ar$ra
    rb=ar$rb
    Rb[j]=ar$rb
    posterior[i]=ar$pA.less.B
    result[j] = posterior[i]
    
  }
  Na[j]=A
  Nb[j]=B
}

我最终计算的手段是

> mean(Na[1:j])
[1] 30.4

> mean(Nb[1:j])
[1] 34.89

这似乎是错误的。您认为我如何解决此错误。如果我没有在代码中输入的所有 100 名患者都经过评估（0.90），您是否认为这可能取决于宣布疗效的阈值？任何人都可能知道如何建立一个正确的命令，如果那是差距？

谢谢

标签： rmodelmeanbayesiantrial