第 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.00716 0.942    -2.86
## 2      1          2     2 -0.00699 0.945    -2.91
## 3      1          3     3 -0.00763 0.960    -3.22
## 4      1          4     4 -0.00661 0.944    -2.89
## 5      1          5     5 -0.00722 0.967    -3.43
## 6      1          6     6 -0.00720 0.969    -3.47

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.946  0.911  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")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.