机器算法验证 - Nassim Nicholas Taleb 对丹麦面具研究数据的分析（完全分离的二项式 GLM） - 吾爱随笔录

Nassim Nicholas Taleb 对丹麦面具研究数据的分析（完全分离的二项式 GLM）

机器算法验证物流二项分布列联表渔民精确测试分离

2022-02-09 11:55:17

最近，纳西姆·尼古拉斯·塔勒布 ( Nassim Nicholas Taleb ) 发表了一篇关于最近丹麦口罩研究的帖子，这是一项随机对照试验，得出的结论是，戴口罩组和不戴口罩对照组的新诊断冠状病毒感染比例没有显着差异 (42/2392 =1.8% 与 53/2470=2.1%），基于逻辑回归的 2 尾 p 值。

Taleb 指出，如果你只关注通过 qRT PCR 测试确认感染的病例，你会得到 0/2392 与 5/2470。他写道：“现在考虑更明显的错误。随机获得 0 个 PCR 与 5 个的概率是多少？

如果平均值为 5/2470，则在 2392 中实现 0 次的概率是 0.0078518，即 127 分之 1。我们可以用 p 值重新表示它，这将是 <.05，远远超过论文显示的范围为 .21-.33。这些研究人员是如何错过这一点的，我无法理解。”

对我来说，这个计算似乎没有多大意义，因为基于二项式 PMF 的计算忽略了 5/2470 的采样不确定性。对此，Taleb 回应说“我不做 P 值：https ://arxiv.org/pdf/1603.07532.pdf ”，只是后来用 Mathematica 表示法添加了这个基于蒙特卡洛双伯努利的计算，其中他获得了1 尾 p 值

ta=Table[data1=RandomVariate[BernoulliDistribution[5/2470],2470]//Total;
         data2=RandomVariate[BernoulliDistribution[5/2470],2400]//Total;
data1-data2,{10^5}];
[[Select[ta,#<-5&]//Length] / 10^5]//N
0.03483

有人可以阐明塔勒布到底在计算什么吗？如果他在这里尝试使用蒙特卡洛完成 2 样本二项式检验，我不太明白 2400 的来源以及为什么那里没有伯努利，期望值为 0/2392（这本身就会有问题，因为那么一个将有零方差）。对于 2 样本 MC 二项式检验，我宁愿期望类似（在 R 中，并使用所有计数的 +1 调整以避免一组中的 p=0 二项式期望）：

p1=rbinom(1E8, 2470, (5+1)/(2470+1))/(2470+1)
p2=rbinom(1E8, 2392, (0+1)/(2392+1))/(2392+1)
mean(p1<=p2) # 1-tailed p = 0.03401019
2*mean(p1<=p2) # 2-tailed p = 0.06802038

但似乎他改为尝试类似的东西（我更正了 2400，这可能是一个错字，并将 > 更改为 >=）：

mean((rbinom(1E7, 2470, 5/2470)-rbinom(1E7, 2392, 5/2470))>=5) 
# 1-tailed p = 0.0811906

我认为这是错误的，对吧？如果有的话，我会发现这更合乎逻辑：

mean((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392)>=(5/2470-0/2392)) 
# 1-tailed p = 0.01446185
mean(abs((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392))>=(5/2470-0/2392)) 
# 2-tailed p = 0.03479425

尽管我不确定这是否是执行此类 2 样本二项式检验的可接受方式（这似乎是Liddell 的 2x2 列联表检验的蒙特卡洛版本）。

塔勒布本人并没有提供多大帮助，指出我们只是在计算“费舍尔的双列表联合分布”，然后用他典型的直率评论说“你似乎非常非常无知，像鹦鹉一样重复公式，不明白什么概率差不多。我已经不再和你交往了。”。

我确实告诉他，因为先验口罩也可能使事情变得更糟（例如，当使用不当时）使用 2 尾测试更安全。而且由于完全分离，您无法进行常规逻辑回归（二项式 GLM）（这是作者在论文中使用的），例如在 R 中：

summary(glm(cbind(pcrpos, pcrneg) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")),pcrpos=c(0,5), pcrneg=c(2392,2470-5))))
# 2-tailed p = 1, obviously not correct

为了解决这个问题，我们可以将 1/2 作为连续性校正添加到我们的观察中（我相信，这相当于在贝叶斯二项式 GLM 中使用杰弗里先验）：

summary(glm(cbind(pcrpos+1/2, pcrneg+1/2) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")),pcrpos=c(0,5), pcrneg=c(2392,2470-5)))) 
# 2-tailed p = 0.11

然后我指出，最好做一个完全类似的逻辑回归，例如在 R 中：

library(elrm)
fit = elrm(pcrpos/n ~ treatment, ~ treatment, r=2, iter=400000, burnIn=1000, dataset=data.frame(treatment=factor(c("masks", "control")), pcrpos=c(0, 5), n=c(2392, 2470)) )
fit$p.values # 2-tailed p value = 0.06
fit$p.values.se # standard error on p value = 0.0003

并且这也将非常接近基于超几何分布的 2 尾Fisher 精确检验的结果，这也给出了 0.06 的 2 尾 p 值：

fisher.test(rbind(c(0,2392), c(5,2470-5))) # 2-tailed p value = 0.06

或 0.03 的 1 尾 p 值：

fisher.test(rbind(c(0,2392), c(5,2470-5)), alternative="less") # 1-tailed p value = 0.03

尽管 Fisher 精确检验会假设行和列边距都是固定的，但这实际上在这里并不完全正确，因为只有行边距是固定的，这将使逻辑回归/2 样本二项式更合适。

我指出的另一种选择是Firth 的逻辑回归，它会给出 0.11 的 2 尾 p 值：

library(brglm)
summary(brglm(cbind(pcrpos, pcrneg) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")), pcrpos=c(0,5), pcrneg=c(2392,2470-5))))
# 2-tailed p = 0.11

对此他回答说：“请不要给我图书馆。请提供推导。” （不要介意即使对于二项式 GLM 的最大似然解也没有封闭形式的解）。

无论如何，这里是否有人能够对整个讨论提供一些反馈，最好是从正式的统计角度，以便可能取悦塔勒布？特别是关于完全分离的问题以及如何在 2 样本二项式检验或逻辑回归中最好地处理它，以及获得精确 p 值的最佳选择是什么。

编辑：多考虑可能的选项，比较两个独立的二项式比例的精确无条件测试可能是最正确的。例如使用 Boschloo 的测试（https://en.wikipedia.org/wiki/Boschloo%27s_test）：

library(Exact)
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Boschloo", alternative="two.sided", model="Binomial") 
# Boschloo's test, 2-tailed p = 0.06223
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Boschloo", alternative="less", model="Binomial") 
# Boschloo's test, 1-tailed p = 0.03196

尽管该exact.test函数似乎有很多不同的方法，而且我不确定哪个是最好的（特别是对于计数低的情况和二项式期望 p=0 的组），因为我还没有挖掘了解所有这些方法的细节。例如method="Z-pooled"，给出更乐观的 p 值，更接近我通过类似 Liddell 的 MC 方法获得的 p 值，以测试上述常见的比例p=5/(2392+2470)：

exact.test(rbind(c(0,2392), c(5,2470-5)), method="Z-pooled", alternative="two.sided", model="Binomial") 
# 2-tailed p-value = 0.02809
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Z-pooled", alternative="less", model="Binomial") 
# 1-tailed p-value = 0.01425

同样，使用library(exact2x2)和使用method="FisherAdj"我得到了更乐观的 p 值：

uncondExact2x2(0, 2392, 5, 2470, alternative="two.sided", method="FisherAdj") 
# 2-tailed p = 0.03417
uncondExact2x2(0, 2392, 5, 2470, alternative="greater", method="FisherAdj") 
# 1-tailed p = 0.01709

欢迎就这些测试中的哪一个在这里最合适的想法进行思考。

另一方面，如果考虑到假阴性 PCR 测试（就像 Taleb 喜欢在抗体测试中考虑假阳性一样），这可能会大大改变结论......也很确定一个需要知道哪些人接受了每种类型的测试，以及他们测试的所有其他呼吸道病毒的计数是多少。

3个回答

您的双面测试隐含地将 5% 显着性水平的一半分配给“口罩有害”（），另一半分配给“口罩有益”（）。对于像塔勒布这样的贝叶斯主义者，这可能表明您没有正确考虑您的先验，因为这意味着您接受所需的证据数量与您接受所需的证据数量完全相同，即使在直觉上更有可能（至少对我来说 - 如果你在一年前问我戴口罩是否更有可能增加或减少感染呼吸道病毒的风险，我会说减少）。 $M_-$ $M_+$ $M_-$ $M_+$ $M_+$

在我看来，由于单个自变量是二元的，使用 Firth 的逻辑回归而不是 Fisher 的精确检验并没有增加太多价值。

但你的核心观点是，在试验中戴口罩减少 PCR 感染的大部分疑问来自估计中存在不确定性。这似乎无可争辩。 $p=\frac{5}{2470}$

塔勒布说：“你不明白。顺便说一句，我添加了一个双列表联合分布，就像费舍尔一样。” 也许这是他在承认你提出了有效批评的同时保住面子的方式。

该模拟试图估计精确测试的结果，其中观察到这种极端差异的机会是

\begin{aligned} Pr (all positive results in one group) & = \frac{(\binom{2470}{5}) + (\binom{2392}{5})}{(\binom{2470 + 2392}{5})} \\ = 0.03377 + 0.02876 = 0.06253. \end{aligned}

$\begin{aligned} \Pr(\text{all positive results in one group}) &= \frac{\binom{2470}{5} + \binom{2392}{5}}{\binom{2470+2392}{5}} \\&= 0.03377 + 0.02876 = 0.06253. \end{aligned}$

假设所有名受试者被独立随机分为测试组和对照组，这一计算是合理的。分子中的项计算所有五个阳性结果可能在第一组或第二组中随机结束的方式。分母计算所有受试者组中所有可能的五元素子集，在此随机化假设下，所有这些子集的可能性均等。 $2470+2392$

对于那些坚持（不正确地，恕我直言）这里适合进行单边测试的人，我已经展示了两个单独分数中的每一个的值。

我相信这个简单的结果等同于 Fisher 检验和 Boschloo 检验。

还可能值得注意的是，塔勒布的蒙特卡洛计算中的标准误差是

se = \sqrt{0.03377 (1 - 0.03377) / 10^{5}} = 0.0005712,

$\operatorname{se} = \sqrt{0.03377(1-0.03377)/10^5} = 0.0005712,$

请记住，他似乎已经接近到（这几乎没有区别）。这使得他的结果比真实值低了大约个标准误差，完全符合人们对这种模拟的期望。 $2392$ $2400$ $0.03483$ $1.9$

如果你想默认
(i) 避免使用 p 值的方法
(ii) 为这个问题产生更“定制”的测试
（Taleb 提到的目标，不一定是你同意的）

一种解决方案是通过拒绝采样（实际上是近似贝叶斯计算）来模拟和查找参数。

因此，假设有两个关键参数：（a）不戴口罩时的感染率和（b）戴口罩时的感染率变化百分比。让我们放置两个统一的先验和：您认为面具可以将风险降低 50% 或将风险增加 50% 的先验（双尾为你建议）。 $r$ $m$
$r \sim U [.0001,.005]$ $m \sim U[.5,1.5]$

我们可以运行一个模拟，在其中绘制这两个参数，使用它们对 2392 个带掩码的观测值和 2470 个不带掩码的观测值进行采样。然而，我们只接受模拟，它在不戴口罩时准确输出 5 个正数，在戴口罩时输出 0 个正数。我们一直这样做，直到我们收集到 5000 个模拟。

当我们查看后验时，乘数分布明显向左移动并远离原始中心为 1。然而，乘数高于 1 的概率仍然是 17% 左右（远高于我们 5% 的阈值）。设置我们自己）。每个人都赢了！

library(tidyverse)
library(progress)

### hard data from the paper
WITH_MASK_OBSERVATIONS<-2392
WITHOUT_MASK_OBSERVATIONS<-2470

WITH_MASK_POSITIVES<-0
WITHOUT_MASK_POSITIVES<-5

### here assume uniform priors on 
### chances of getting sick
PRIOR_INFECTION_RISK_MASKLESS<-function(){
  runif(n=1,min=.0001,max=.005)
}
### let's assume we have a uniform prior that masks can do anything
### between cutting risks by 50% and increasing it by 50%
PRIOR_CHANGE_RISK_MASK<-
  function(){
    runif(n=1,min=.5,max=1.5)
  }


### run simple ABC where we look for parameters
### where we get exactly the observations above
accepted_runs<-list()
TARGET_ACCEPTED_RUNS<-5000
pb <- progress_bar$new(total = TARGET_ACCEPTED_RUNS) ## to watch the time go by
while(length(accepted_runs)<TARGET_ACCEPTED_RUNS)
{
  ## draw a "real" infection rate from your priors
  infection_maskless<-PRIOR_INFECTION_RISK_MASKLESS()
  mask_multiplier<-PRIOR_CHANGE_RISK_MASK()
  infection_masked<- mask_multiplier * infection_maskless
  ## observe the infection rate!
  observed_infections_maskless<- sum(
rbinom(n=WITHOUT_MASK_OBSERVATIONS,
       size=1,
       prob=infection_maskless))
  observed_infections_mask<-sum(
rbinom(n=WITH_MASK_OBSERVATIONS,
       size=1,
       prob=infection_masked))
  ## if this is EXACTLY what we observe, store the drawn infection rates
  ## they will be part of our posterior
  if(observed_infections_maskless==WITHOUT_MASK_POSITIVES &&
 observed_infections_mask==WITH_MASK_POSITIVES
  ){
pb$tick()
    accepted_runs<-
      append(accepted_runs,
             list(data.frame(maskless = infection_maskless,
                             masked = infection_masked,
                             multiplier = mask_multiplier)))
  }
  
}

其它你可能感兴趣的问题

上一篇拟合多元、自然三次样条下一篇对时间序列数据使用随机森林