COVID-19 研究中的样本量计算

机器算法验证 规模效应 统计能力 渔民精确测试
2022-04-07 03:46:55

来自Boulware 等人的一项 Hydroxychloroquine 作为 Covid-19 暴露后预防的随机试验。在《新英格兰医学杂志》(https://www.nejm.org/doi/full/10.1056/NEJMoa2016638)中,当涉及到适当的样本量计算时,我对以下句子感到好奇:

使用具有 50% 相对效应大小以减少新的症状性感染、双边 α 为 0.05 和 90% 功效的 Fisher 精确方法,我们估计每组需要招募 621 人。

我对如何执行此计算感兴趣。我从未听说过在 Fisher 精确检验的背景下使用“效应大小”(我熟悉 Coehn 的),而且我不确定在这种情况下功率计算将如何工作(什么是适当的替代假设?)。d

请记住,我在临床试验方面的专业知识为零。我对 Casella 和 Berger 文本级别的统计数据非常满意。

教科书和期刊文章对进一步学习非常有帮助。

4个回答

我知道我迟到了几个月,但只想回复其他答案。所有答案都使用模拟和/或声称精确的 Fisher 计算计算量太大。如果你有效地编码,你可以很快得到精确的计算。下面是示例代码fisherpower()函数与power.exact.test()Exact R 包中函数的比较时间:

    > system.time(power1 <- fisherpower(0.1,0.05,621))
       user  system elapsed 
     698.23    0.93  700.23 
    > system.time(power2 <- Exact::power.exact.test(n1=621, 
                   n2=621, p1=0.1, p2=0.05, 
                   method="Fisher")$power)
       user  system elapsed 
       0.32    0.00    0.33 
    
    > power1
    [1] 0.9076656
    > power2
    [1] 0.9076656

使用函数计算只需要0.33s,而使用power.exact.test()函数计算需要700s fisherpower()请注意,该power.exact.test()函数无需模拟即可计算确切的功率,因此没有不确定性,并且比模拟更快。我还强烈建议使用 Barnard 精确检验而不是 Fisher 精确检验来比较两个比例。以下是随着组样本量增加的功效计算:

    nGroup <- 570:630
    powerFisher <- vapply(nGroup,
                          FUN = function(xn) {
                            Exact::power.exact.test(n1=xn, 
          n2=xn, p1=0.1, p2=0.05, method="Fisher")$power
                      }, numeric(1) )
powerBarnard <- vapply(nGroup,
                      FUN = function(xn) {
                        Exact::power.exact.test(n1=xn, 
      n2=xn, p1=0.1, p2=0.05, method="Z-pooled")$power
                          }, numeric(1) )
    
    plot(NA, xlim=range(nGroup), ylim = c(0.85,0.95), 
     xlab="Sample Size per Group", ylab = "Power")
    lines(nGroup, powerFisher, col='red', lwd=2)
    points(nGroup, powerFisher, pch = 21, col = 'red', bg = 
     "red", cex = 0.8)
    lines(nGroup, powerBarnard, col='blue', lwd=2)
    points(nGroup, powerBarnard, pch = 21, col = 'blue', bg = 
     "blue", cex = 0.8)
    
    abline(h=0.9, lty=2)
    abline(v=c(579, 606), col=c('blue', 'red'))
    legend(610, 0.875, c("Barnard", "Fisher"), col = c('blue', 
     'red'), lty = 1, pch=21, pt.bg=c('blue', 'red'), cex=1.2)

在此处输入图像描述

@heropup 正确的是组样本大小应该是 606(不是 621),如图所示。然而,Barnard 的检验更强大,使用“Z-pooled”检验统计量,每组只需要 579 名参与者。由于这是一个罕见的事件,因此可能需要使用 Berger 和 Boos (1994) 区间方法,该方法将样本量减少到 573 名参与者(代码未显示,需要一些时间)。重要的是,这些替代方案仍然控制 1 类错误率,并且简单地优于 Fisher 对 2x2 表的精确检验。对于分析数据集,我建议使用Exact::exact.test()@SextusEmpiricus 提供的示例数据集仅需 0.3 秒,而不是 Barnard::barnard.test()47 秒。但是,两者都产生相同的结果,而且我是 Exact R 包的维护者,所以可能会有偏见。

一个油嘴滑舌的答案是,他们可能只是将他们的数字插入了一个功率计算器。我附上了在 G*Power 3.1 中重新创建此功率分析的屏幕截图,G*Power 3.1 是一个免费提供的功率计算器。注意要匹配他们的 621 结果,我必须转到“选项”并选择“最大化 Alpha”。

该论文称,“我们预计 10% 的接触 Covid-19 的密切接触者会出现与 Covid-19 相容的疾病”以及“50% 的相对影响大小”。我将第二部分解释为他们假设治疗的效果会将患病率从 10% 降低到 5%。

这导致比例 p1 和 p2 的值分别为0.050.1

遗憾的是,我不知道 G*Power 是如何进行这个计算的,但我至少可以尝试解释一下这个想法。

我们的比例为 0.1 和 0.05。对于给定的样本大小,我们可以通过从两个二项式随机变量中抽样来随机抽样 2x2 列联表。功效计算询问“Fischer 精确检验多久会拒绝使用此过程创建的列联表的原假设?”。n

特别是,我们希望找到最小的,这样 Fischer 的检验将在至少 90% 的时间内拒绝原假设。n

一种近似的方法是模拟。对于给定的,假设有 10,000 个列联表,运行 Fischer 检验,并查看 p 值低于 0.05 的频率。继续增加直到 p 值在 90% 或更多的情况下低于 0.05...nn

G*Power 3.1 截图

他们使用了Fisher精确检验,该检验涉及无放回抽样。

但实际上这并不完全如此,它更像是二项式分布式数据。

对于这种情况,您会得到以下信息:

  • 对于零假设,它是在人们感染 covid-19 的概率相等的情况下进行抽样,无论他们是在安慰剂组还是效果组。

  • 他们计算功效的替代假设是,安慰剂组有 10% 的概率获得 covid-19,而治疗组有 5% 的概率(因此治疗将概率降低了 50%)。


精确计算功率

您可以简单地通过尝试所有可能性来计算在给定特定样本大小和概率的情况下拒绝原假设的概率,并查看哪些导致负/正 Fisher 检验。然后你对概率求和,得到你拒绝测试的情况。

P(reject)=over all i,jwhere Fisher test is rejectedP(i placebo cases and j treatment cases)

下面是一个代码示例

fisherpower <- function(p1, p2, n) {
  pf <- 0
  for (i in 1:n) {
    for (j in 1:n) {
      M <- matrix(c(i,n-i,j,n-j),2)
      if (fisher.test(M)$p.value <= 0.05) {
        pf <- pf + dbinom(i,n,p1)*dbinom(j,n,p2)
      }
    }
  }
  pf
}

这使

> fisherpower(0.1,0.05,621)
[1] 0.9076656

然而,这种方法消耗了大量的计算能力。你需要尝试 621 次 621 种可能性。上面的实现可以改进很多(您不需要计算所有 621 x 621 的情况),但它仍然很慢,所以这就是 R 中的标准实现使用模拟的原因。上述内容的快速实现在 Peter Calhoun 的 R 包Exact中,他在此处的回答中对此进行了解释。


模拟计算

您多次计算假设结果,并针对该结果确定 5% 假设检验是否会失败。

作为您获得的样本量的函数:

  • 如果原假设为真,那么您将始终获得 5% 的拒绝概率。

实际上,这并不完全正确,并且当条件不正确时,Fisher 精确检验略微保守。即使原假设为真(在我们不带放回抽样的情况下),Fisher 精确检验的拒绝率也低于 5%。在下面的示例图中,我们计算时的拒绝概率(在这种情况下,null 为真)。p1=p2=0.1

  • 如果原假设为假,并且概率不相等。然后,当样本量较大时,您会获得更大的概率来拒绝原假设。

显示功率为 n 的函数

### computing 
set.seed(1)
n <- seq(100,1000,20) 
power <- sapply(n, 
                FUN = function(xn) {
                  statmod::power.fisher.test(0.1,0.05,xn,xn, nsim = 10000)
                } )
type1 <- sapply(n, 
                FUN = function(xn) {
                  statmod::power.fisher.test(0.1,0.1,xn,xn, nsim = 10000)
                } )

### plotting of results
plot(n,power, type = "l", ylim = c(0,1),
     ylab = "reject probability")
lines(n,type1, col =2)
points(n,power, pch = 21, col = 1, bg = "white", cex = 0.7)
points(n,type1, pch = 21, col = 2, bg = "white", cex = 0.7)

# lines at 0.05 and 0.9
lines(c(0,2000),c(0.05,0.05), col = 2, lty = 2)
lines(c(0,2000),c(0.9,0.9), col = 1, lty = 2)

# legend
legend(1000,0.6,c("if p1 = p2 = 0.1",
                  "if p1 = 0.1, p2 = 0.05"), title = "reject probability",
       col = c(2,1), lty = 1, cex = 0.7, xjust = 1
      )

替代测试

还有很多其他的方式来看待它。我们还可以进行 Barnards 测试

> Barnard::barnard.test(49,58,414-49,407-58)

Barnard's Unconditional Test

           Treatment I Treatment II
Outcome I           49           58
Outcome II         365          349

Null hypothesis: Treatments have no effect on the outcomes
Score statistic = 1.02759
Nuisance parameter = 0.012 (One sided), 0.986 (Two sided)
P-value = 0.16485 (One sided), 0.320387 (Two sided)

或使用 GLM 模型

> summary(glm(cbind(c(49,58),c(414-49, 407-58)) ~ 1+c("chloroquine", "placebo"), family = binomial(link="identity")))

Call:
glm(formula = cbind(c(49, 58), c(414 - 49, 407 - 58)) ~ 1 + c("chloroquine", 
    "placebo"), family = binomial(link = "identity"))

Deviance Residuals: 
[1]  0  0

Coefficients:
                                   Estimate Std. Error
(Intercept)                         0.11836    0.01588
c("chloroquine", "placebo")placebo  0.02415    0.02350
                                   z value Pr(>|z|)    
(Intercept)                          7.455 8.98e-14 ***
c("chloroquine", "placebo")placebo   1.028    0.304    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.0568e+00  on 1  degrees of freedom
Residual deviance: 2.4780e-13  on 0  degrees of freedom
AIC: 15.355

Number of Fisher Scoring iterations: 2

> 

这些方法中的每一种都或多或少地显示相同的东西,结果 58 vs 49 不是异常(而且,效果需要达到 50% 或更多才能使我们有至少 90% 的概率检测到异常用这个测试)。

您错过了文章在引用之前引用的一条关键信息:

我们预计 10% 的接触 Covid-19 的密切接触者会出现与 Covid-19 相容的疾病。

这是在备择假设下假设的对照组发生率;50% 的相对效应量是指治疗组中 Covid-19 感染发生率的降低,即,由此得出,在备择假设下。πc=0.1πt/πc=0.5πt=0.05

但是,当我将这些(连同)输入 EAST 6 时,我没有得到每臂我得到,根据我的模拟,我相信后一个值是正确的。αβn=621n=606