【1.3.1】疫苗的保护率

primer 阅读量：次

我们知道正态分布，二项分布等常见分布，但是对于Beta分布，可能很多人并不太熟悉。在估计疫苗保护率之前，我们先看一下Beta分布的特性。

一、Beta分布

在二项分布中，某个事件是以一个固定概率出现的，比如掷硬币，那么正面朝上的概率我们坚信是0.5，但是硬币可能并非完全对称重。比如，有可能正面朝上的概率是0.499，或者0.501。那么此时在二项分布中的事件发生概率就很难用一个固定的值来表示了，通常做法是用概率分布来替代单个的点估计。常用的概率分布就是Beta分布。也就是说Beta分布是描述事件发生概率的分布。更广泛的说，Beta分布适用于百分比数据，它的取值分为为0-1之间，累积概率密度和为1.

其概率密度函数：

$ f(x;\alpha,\beta) = \frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1} $

其中B(α，β)是beta函数，主要为了使全部累积概率为1，我们这儿暂且忽略。只考虑后面含有x的项。

Beta分布有两个正数参数：α>0和β>0，这两个参数共同决定了Beta分布的形状。

比如：

Beta分布的期望均值为：

$ E[X] = \frac{\alpha}{\alpha+\beta} \\ $

当的时候，整个概率密度分布是对称的，期望均值为0.5;

当的时候, 概率密度分布偏左（靠近0）；

当的时候, 概率密度分布偏有（靠近1）；

Beta分布在Bayes推断中可以作为一个先验概率分布。比如，我们曾经掷硬币100次，49次正面朝上，51次反面朝上，那么我们可以认为掷硬币正面朝上的概率分布为 Beta(49,51)

curve(dbeta(x,49,51))

随着掷硬币次数的增加，概率分布变得更加瘦长，更加集中。

假如上面的Beta(49,51)是我们的一个先验概率分布。现在我们再次某个硬币，一共掷了1000次，其中正面朝上只出现了100次，其余900次均为反面朝上。那么此时的后验概率分布Beta则变成了 Beta(49+100,51+900)

curve(dbeta(x,49+100,51+900))

可以看到此时掷硬币正面朝上的概率明显偏低了，平均期望概率只有：

$ \frac{\alpha}{\alpha + \beta} = \frac{149}{149+951} = 0.135 $

所以，掷同一个硬币，下次出现正面朝上的概率只有13.5%，这儿多出的3.5%是因为我们先验概率的原因。随着我们实验次数的增加，先验概率对结果的影响会越来越弱。

二、疫苗保护率

2020年12月，BNT在《新英格兰医学杂志》上发布了其疫苗实验数据，如下：

有两种观测点，我们以第一个观测为例，实验组有17411人参与，其中8人感染新冠；对照组有17511个参与，其中162人感染新冠。以此为数据，实验推断疫苗的保护率为95.0%，置信区间为[90.3%, 97.6%]，疫苗保护率超过30%的概率超过99.99%。

那么这组数据是怎么计算的来的？

三、疫苗保护率

$ Vaccine\ Efficacy = 100 \times (1-IRR) $

其中，

$ IRR = \frac{疫苗组发病率}{对照组发病率} $

按照文章中的数据，我们可以计算

$ IRR = \frac{8/17411}{162/17511} = 0.04966 $

所以疫苗保护率

$ Efficacy = 100\times(1-IRR) = 95\% $

四、 I型错误

上述95%是一个点估计值，实际疫苗的保护率应该是一个分布。假设疫苗上市要满足保护率超过30%，那么我们还要知道在保护率的分布中有多大面积是超过30%的阈值的。按照文章中给的结果 ，

$ P(VaccineEfficacy ＞ 30%) = 0.9999 $

也就是说保护率分布中只有0.01%的面积在30%保护率以下，也即I型错误的概率只有0.01%

五、先验概率

上图中，我们直接给出了有效率的后验概率分布，那么这个分布是怎么得来的？

根据贝叶斯理论，疫苗后验概率分布

$ (Posterior \ beliefs) \propto (Prior \ beliefs) \times (Likelihood \ of \ observed \ data ) $

根据文中开始说到的Beta分布，我们可以给一个保护率的先验概率分布，其服从Beta分布：

$ Prior \ beliefs = constant \times \theta^{ \alpha-1}(1- \theta)^{ \beta-1 } $

似然函数则为二项分布函数，也即观测了xx人数，有xx人发病。

$ Likelihood \ of \ observed \ data = \binom{n}{k} \theta^{k}(1-\theta)^{n-k } $

将上述两式相乘，即可推导出保护率的后验概率分布：

$ Posterior \ beliefs \propto \theta^{\alpha+k-1}(1-\theta)^{\beta+n-k-1} $

上述结果可见，后验概率分布也符合 $ Beta(\alpha+k, \beta+n-k )$的Beta分布。

更概括一点说，如果一个似然函数是二项分布，先验分布是Beta分布，那么该模型称之为Beta二项分布模型，简写如下：

$ k \sim Binomial(n,\theta) $

$ \theta \sim Beta(\alpha,\beta) $

根据观察，我们认为人群中感染新冠的期望平均先验概率为0.01，根据Beta分布均值的求算公式，我们有

$ \frac{\alpha}{\alpha+\beta} = 0.01 $

对于α 和 β 我们同样至少给其中一个先验值。这种情况下，我们可以给一个较弱的先验值，因为一个较弱的先验值能够被我们的观测很容易地更改，以减少先验值对后验概率的影响程度。我们知道，对于和，值越大，先验信息越强，对后验概率影响越大。所以，我们这儿给一个较弱的先验值β =1。根据上式此时可以推算 α= 0.010101

此时，我们的先验概率服从如下的Beta分布：

curve(dbeta(x,0.010101,1), xlab="Incidence Rates")

你也可能会问，问什么不是α =1？如果α =1，那么β =99，不是一个弱信息先验概率。

六、似然函数

有了上述先验概率分布，下面我们加入我们的实验数据，根据我们上面的推导公式，后验概率会同样服从Beta分布，

$ Vaccine\ Incidence\ Rate \sim Beta(0.010101+8, 1+17411-8) $

$ Placebo\ Incidence\ Rate \sim Beta(0.010101+162, 1+17511-162) $

此时，实验组和对照组的患病率分布如下：

curve(dbeta(x,0.010101+8, 1+17411-8), col = "red", xlim=c(0,0.02), xlab="Incidence Rates", ylab="Density")
curve(dbeta(x,0.010101+162, 1+17511-162), col = "black", add=TRUE)
legend(0.015, 2000, legend=c("Vaccine","Placebo"), col=c("red","black"), lty=1)

【注：图中横坐标非0-1】

上图后验患病概率分布可以看出，安慰剂组的患病率明显要高于疫苗组的患病率。

蒙特卡洛抽样

更多内容，参见： https://zhuanlan.zhihu.com/p/396483260

参考资料

药企，独角兽，苏州。团队长期招人，感兴趣的都可以发邮件聊聊：tiehan@sina.cn
个人公众号，比较懒，很少更新，可以在上面提问题，如果回复不及时，可发邮件给我： tiehan@sina.cn