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

## 73.1 引言

group volunteers got_covid
placebo 22000 162
vaccinated 22000 8

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

## 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}

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

## 73.3 计算

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 %>%

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

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")