The following is from an email that Ian Fellows sent to ISOSTAT on 11/16/2020.
I wrote up two posts on this that may be of interest.
http://blog.fellstat.com/?p=468
http://blog.fellstat.com/?p=440
I think the methods in those posts could be great for an intro Bayesian course. For intro stats, this is a fantastic example. The hypothesis test here is just a simple binomial test!
Let’s do Moderna. It breaks down like this. 5 treated and 90 control cases were observed. Now, conditional on the fact that we’ve observed a case and there is no treatment effect, what is the probability that that case belongs to the treatment group? Answer: 0.5. Now scale it up. What is the distribution of the number of treatment cases given we’ve seen 95 total cases and there is no treatment effect? Answer: Binomial(0.5, 95)
Okay, so now all we need to do is run a test for H0: p = 0.5 and we have the answer as to whether the vaccine is more effective than placebo.
But… that is not really useful. Vaccines need to be A LOT more effective than placebo to make a difference, so we should really look at the confidence interval. Doing so gives us the interval for p, but we need to transform that to get vaccine efficacy. Efficacy is one minus the rate in the treatment group over the rate in the placebo group, or put another way one minus the odds ( VE = 1 - p / (1 - p) ).
That’s all there is to it.
(bt <- binom.test(5,95))##
## Exact binomial test
##
## data: 5 and 95
## number of successes = 5, number of trials = 95, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.01730788 0.11856279
## sample estimates:
## probability of success
## 0.05263158 p <- bt$estimate
p.ci <- bt$conf.int
print(c("Efficacy: ",Eff <- 1 - p/(1-p)))## probability of success
## "Efficacy: " "0.944444444444444" print(c("CI: ",Eff.CI <- 1 - (p.ci / (1 - p.ci))[c(2,1)]))## [1] "CI: " "0.865489240414379" "0.982387285477282"