第 73 章 贝叶斯分析案例-新冠疫苗有效率的计算

73.1 引言

纽约时报报道说,

美国制药公司辉瑞(Pfizer)和德国生物科技公司(BioNTech)11月9日率先宣布 ,根据在数国临床试验初步结果,其研发的新冠疫苗有效率达到90%以上,星期三,完整结果显示,参加疫苗实验的44000个志愿者中,共有170人确诊感染,其中安慰剂组162人,接种疫苗组仅8人,这证明了辉瑞开发的新冠疫苗有效率高达95%。

group volunteers got_covid
placebo 22000 162
vaccinated 22000 8

新冠疫苗是有效的,且有效率高达95%。 那么,这个95%是怎么计算出来的呢?它的概率是多少以及不确定性是多少呢? 回到这个问题,我们首先需要了解,辉瑞公司是如何定义疫苗有效率的

\[ \text{VE} = 1 - \frac{p_{t}}{p_{c}} \]

其中\(p_t\)疫苗组(vaccinated)的感染率,\(p_c\)安慰剂组(placebo)的感染率。

73.2 模型

library(tidyverse)
library(tidybayes)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

然后,我们建立如下数学模型:

\[ \begin{align} y_{c} \sim \textsf{binomial}(n_{c},p_{c}) \\ y_{t} \sim \textsf{binomial}(n_{t},p_{t}) \\ p_{c} \sim \textsf{beta}(1, 1) \\ p_{t} \sim \textsf{beta}(1, 1) \end{align} \]

通过模型可以直接计算干预效果\(\textsf{effect}\)和疫苗有效性\(VE\)

\[ \begin{align} \text{effect} = p_{t} - p_{c} \\ \text{VE} = 1 - \frac{p_{t}}{p_{c}} \end{align} \]

73.3 计算

具体Stan代码如下

stan_program <- "
data {
  int<lower=1> event_c;        // num events, control
  int<lower=1> event_t;        // num events, treatment
  int<lower=1> n_c;            // num of person trial, control
  int<lower=1> n_t;            // num of person trial, treatment
}
parameters {
  real<lower=0,upper=1> p_c;    
  real<lower=0,upper=1> p_t;    
}
model {
  event_c ~ binomial(n_c, p_c);
  event_t ~ binomial(n_t, p_t);
  p_c ~ beta(1, 1);
  p_t ~ beta(1, 1);
}
generated quantities {
  real effect   = p_t - p_c;
  real VE       = 1- p_t /p_c;
  real log_odds = log(p_t / (1- p_t)) - log(p_c / (1- p_c));
}
"


stan_data <- list(
  event_c = 162,
  event_t = 8,
  n_c     = 4.4e4 / 2,
  n_t     = 4.4e4 / 2
)

mod_vaccine <- stan(model_code = stan_program, data = stan_data)

73.4 结果

最后,我们后验概率抽样

draws <- mod_vaccine %>%
  tidybayes::spread_draws(effect, VE, log_odds)

draws %>% 
  head()
## # A tibble: 6 × 6
##   .chain .iteration .draw   effect    VE log_odds
##    <int>      <int> <int>    <dbl> <dbl>    <dbl>
## 1      1          1     1 -0.00674 0.944    -2.89
## 2      1          2     2 -0.00684 0.940    -2.82
## 3      1          3     3 -0.00777 0.961    -3.26
## 4      1          4     4 -0.00751 0.961    -3.25
## 5      1          5     5 -0.00697 0.947    -2.95
## 6      1          6     6 -0.00758 0.924    -2.59

73.4.1 干预效果

从结果中看到effect中很多负数。事实上,effect中越多的负值,即被感染的可能性越低,说明疫苗干预效果越好

mean(draws$effect < 0) %>% round(2)
## [1] 1

结果告诉我们,疫苗有明显的干预效果。比如,我们假定10000个人接受了疫苗,那么被感染的人数以及相应的可能性,如下图

draws %>%
  ggplot(aes(x = effect * 1e4)) +
  geom_density(fill = "blue", alpha = .2) +
  expand_limits(y = 0) +
  theme_minimal() +
  xlab("效应大小") +
  ggtitle("每10000个接种疫苗的人中被感染新冠的数量")

73.4.2 疫苗有效率

我们再看看疫苗有效率 VE 的结果

draws %>%
  select(VE) %>%
  ggdist::median_qi(.width = c(0.90))
## # A tibble: 1 × 6
##      VE .lower .upper .width .point .interval
##   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 0.947  0.908  0.971    0.9 median qi

通过数据看出,疫苗的有效性为0.95,在90%的可信赖水平, 中位数区间[0.91, 0.97].

当然,通过图可能理解的更清晰。

label_txt <- paste("median =", round(median(draws$VE), 2))

draws %>%
  ggplot(aes(x = VE)) +
  geom_density(fill = "blue", alpha = .2) +
  expand_limits(y = 0) +
  theme_minimal() +
  geom_vline(xintercept = median(draws$VE), size = 0.2) +
  annotate("text", x = 0.958, y = 10, label = label_txt, size = 3) +
  xlab("疫苗有效率") +
  ggtitle("辉瑞公司定义疫苗有效率为 VE = 1 - Pt/Pc")