Which statistical model is being used in the Pfizer study design for vaccine efficacy?

I know there's a similar question here:

How to calculate 95% CI of vaccine with 90% efficacy?

but it doesn't have an answer, at the moment. Also, my question is different: the other question asks how to compute VE, using functions from a R package. I want to know why vaccine efficacy is defined as illustrated at the bottom of this page:

$$\text{VE} = 1 - \text{IRR}$$

where

$$\text{IRR} = \frac{\text{illness rate in vaccine group}}{\text{illness rate in placebo group}}$$

and which is the statistical model behind it.

My attempts: I thought the researches would fit a logistic regression model a single binary predictor $$X$$, identifying subjects who got the vaccine ($$X=1$$) or not ($$X=0$$):

$$p(Y|X) = \frac{1}{1+\exp{-(\beta_0 +\beta_1 X)}}$$

However, this is clearly not the case, because for the Moderna vaccine we know that there were 5 cases in the vaccine arm, and 90 in the placebo arm, which corresponds to a $$\text{VE}$$ of $$94.\bar{4}\%$$. These data alone are enough to determine $$\text{VE}$$, but surely they're not enough to fit a LR model, and thus to determine $$\beta_1$$.

Also, by looking at page 111-113 of the Pfizer document, it looks like a different (Bayesian?) analysis is performed. Again, the point estimate seems to be $$\text{VE} = 1 - \text{IRR}$$, but the power of a test is mentioned, and two tables 7 and 8 are presented which show probability of success and failure. Can you show me how to obtain the results in such tables?

• skranz.github.io//r/2020/11/11/… might be of interest Commented Nov 17, 2020 at 12:31
• Andrew Gelman's blog had some discussion related to this here and here. Commented Nov 17, 2020 at 23:38
• @DeltaIV, yeah, fine (although I hope the reason is not that you aim to chat with someone with expertise in biostatistics...:-))! Some time in the afternoon (Berlin time)? Commented Nov 18, 2020 at 7:53
• @eipi10 wow, Andrew too? The stuff to read on this topic is growing rapidly :-) Commented Nov 18, 2020 at 7:58
• @DeltaIV, sounds good! Commented Nov 18, 2020 at 7:59

The relation between efficiency and illness risk ratio

I want to know why vaccine efficacy is defined as illustrated at the bottom of this page:

$$\text{VE} = 1 - \text{IRR}$$

where

$$\text{IRR} = \frac{\text{illness rate in vaccine group}}{\text{illness rate in placebo group}}$$

This is just a definition. Possibly the following expression may help you to get a different intuition about it

$$\begin{array}{} VE &=& \text{relative illness rate reduction}\\ &=& \frac{\text{change (reduction) in illness rate}}{\text{illness rate}}\\ &=& \frac{\text{illness rate in placebo group} -\text{illness rate in vaccine group}}{\text{illness rate in placebo group}}\\ &=& 1-IRR \end{array}$$

Modelling with logistic regression

These data alone are enough to determine $$\text{VE}$$, but surely they're not enough to fit a LR model, and thus to determine $$\beta_1$$.

Note that

$$\text{logit}(p(Y|X)) = \log \left( \frac{p(Y|X)}{1-p(Y|X)} \right) = \beta_0 + \beta_1 X$$

and given the two observations $$\text{logit}(p(Y|X=0))$$ and $$\text{logit}(p(Y|X=1))$$ the two parameters $$\beta_0$$ and $$\beta_1$$ can be computed

R-code example:

Note the below code uses cbind in the glm function. For more about entering this see this answer here.

vaccindata <- data.frame(sick    = c(5,90),
healthy = c(15000-5,15000-90),
X       = c(1,0)
)
mod <- glm(cbind(sick,healthy) ~ X, family = binomial, data = vaccindata)
summary(mod)


This gives the result:

Call:
glm(formula = cbind(sick, healthy) ~ X, family = binomial, data = vaccindata)

Deviance Residuals:
[1]  0  0

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -5.1100     0.1057 -48.332  < 2e-16 ***
X            -2.8961     0.4596  -6.301 2.96e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 9.2763e+01  on 1  degrees of freedom
Residual deviance: 2.3825e-12  on 0  degrees of freedom
AIC: 13.814

Number of Fisher Scoring iterations: 3


So the parameter $$\beta_1$$ is estimated as $$-2.8961$$ with standard deviation $$0.4596$$

From this, you can compute (estimate) the odds, the efficiency, and their confidence intervals. See also: How exactly is the "effectiveness" in the Moderna and Pfizer vaccine trials estimated?

The Bayesian model (Table 6)

Also, by looking at page 111-113 of the Pfizer document, it looks like a different (Bayesian?) analysis is performed. Again, the point estimate seems to be $$\text{VE} = 1 - \text{IRR}$$, but the power of a test is mentioned, and two tables 7 and 8 are presented which show probability of success and failure. Can you show me how to obtain the results in such tables?

These analyses are performed in an early stage to verify whether, given the outcomes, the vaccine is effective. The tables give hypothetical observations for which they would reach the tipping point to declare either failure (posterior probability of success <5%) or great success (the probability that VE>30% is larger than 0.995).

These percentages for the tipping points are actually based on controlling Type I error (more about that below). They control the overall type I error, but it is not clear how this is distributed among the multiple go/no-go points.

The outcome considered is the ratio/count of vaccinated people among all infected people. Conditional on the total infected people this ratio follows a binomial distribution*. For more details about the computation of the posterior in this case see: How does the beta prior affect the posterior under a binomial likelihood

*There is probably a question here about that; I still have to find a link for this; but you can derive this based on the idea that both groups are approximately Poisson distributed (more precisely they are binomial distributed) and the probability to observe a specific combination of cases $$k$$ and $$n-k$$ conditional on reaching $$n$$ total cases is $$\frac{\lambda_1^k e^{-\lambda_1}/k! \cdot \lambda_2^{n-k}e^{-\lambda_2}/(n-k)! }{\lambda_2^ne^{-(\lambda_1\lambda_2)}/n! } = {n \choose k} \left(\frac{\lambda_1}{\lambda_1+\lambda_2}\right)^k \left(1- \frac{\lambda_1}{\lambda_1+\lambda_2}\right)^{n-l}$$

The graphic below shows a plot for the output for these type of computations

• Success boundary This is computed by the posterior distribution for the value $$\begin{array}{}\theta &=& (1-VE)/(2-VE)\\ &=& RR/(1-RR) \\&=& \text{vaccinated among infected}\end{array}$$ For instance the in case of 6 vaccinated and 26 placebo among the first 32 infected people the posterior is Beta distributed with parameters 0.7+6 and 1+26 and the cumulative distribution for $$\theta < (1-0.3)/(2-0.3)$$ will be $$\approx 0.996476$$ for 7 vaccinated and 25 placebo it will be 0.989 which is below the level. In R you would compute these figures as pbeta(7/17,0.700102+6,1+26)

• Futility boundary For this they compute the probability of success which is the power of the test. Say for a given hypothesis the test criterium can be to observe 53 or less cases in the vaccine group among the first 164 cases. Then as function of the true VE you can estimate how probable it is to pass the test.

In the table 6 they compute this not as a function of a single VE, but as an integral over the posterior distribution of the VE or $$\theta$$ (and this $$\theta$$ is beta distributed and the test result will be beta-binomial distributed). It seems like they used something like the following:

 ### predict the probability of success (observing 53 or less in 164 cases at the end)
### k is the number of infections from vaccine
### n is the total number of infections
### based on k and n the posterior distribution can be computed
### based on the posterior distribution (which is a beta distribution)
### we can compute the success probability

predictedPOS <- function(k,n) {
#### posterior alpha and beta
alpha = 0.7+k
beta = 1+n-k
### dispersion and mean
s = alpha + beta
m = alpha/(alpha+beta)
### probability to observe 53 or less out of 164 in final test
### given we allread have observed k out of n (so 53-k to go for the next 164-n infections)
POS <- rmutil::pbetabinom(53-k,164-n,m,s)
return(POS)
}

# 0.03114652
predictedPOS(15,32)
# 0.02486854
predictedPOS(26,62)
# 0.04704588
predictedPOS(35,92)

# 0.07194807
predictedPOS(14,32)
# 0.07194807
predictedPOS(25,62)
# 0.05228662
predictedPOS(34,92)


The values 14, 25, 34 are the highest values for which the posterior POS is still above 0.05. For the values 15, 26, 35 it is below.

Controlling type I error (Table 7 and 8)

Table 7 and 8 give an analysis for the probability to succeed given a certain VE (they display for 30, 50, 60, 70, 80%). It gives the probability that the analysis passes the criterium for success during one of the interim analyses or with the final analysis.

The first column is easy to compute. It is binomially distributed. E.g. The probabilities 0.006, 0.054, 0.150, 0.368, 0.722 in the first columns are the the probability to have 6 cases or less when $$p=(100-VE)/(200-VE)$$ and $$n = 32$$.

The other columns are not similar binomial distributions. They represent the probability of reaching the success criterium if there wasn't success during the earlier analysis. I am not sure how they computed this (they refer to a statistical analysis plan, SAP, but it is unclear where this can be found and if it is open access). However, we can simulate it with some R-code

### function to simulate succes for vaccine efficiency analysis
sim <- function(true_p = 0.3) {
p <- (1-true_p)/(2-true_p)
numbers <- c(32,62,92,120,164)
success <- c(6,15,25,35,53)
failure <- c(15,26,35)
n <- c()
### simulate whether the infection cases are from vaccine or placebo group
n[1] <- rbinom(1,numbers[1],p)
n[2] <- rbinom(1,numbers[2]-numbers[1],p)
n[3] <- rbinom(1,numbers[3]-numbers[2],p)
n[4] <- rbinom(1,numbers[4]-numbers[3],p)
n[5] <- rbinom(1,numbers[5]-numbers[4],p)

### days with succes or failure
s <- cumsum(n) <= success
f <- cumsum(n)[1:3] >= failure

### earliest day with success or failure
min_s <- min(which(s==TRUE),7)
min_f <- min(which(f==TRUE),6)

### check whether success occured before failure
### if no success occured then it has value 7 and will be highest
### if no failure occured then it will be 6 and be highest unless no success occured either
result <- (min_s<min_f)

return(result)
}

### compute power (probability of success)
### for different efficienc<y of vaccine
set.seed(1)
nt <- 10^5
x <- c(sum(replicate(nt,sim(0.3)))/nt,
sum(replicate(nt,sim(0.5)))/nt,
sum(replicate(nt,sim(0.6)))/nt,
sum(replicate(nt,sim(0.7)))/nt,
sum(replicate(nt,sim(0.8)))/nt)
x


This gives 0.02073 0.43670 0.86610 0.99465 0.99992 which is close to the overall probability of success in the final column.

Although they use a Bayesian analysis to compute values in table 6. They have chosen the boundaries, based on which they performed the Bayesian analysis, according to controlling the type I error (I think that they use the probability to have success given VE = 0.3, p=0.021, as the basis for the type I error. This means that if the true VE = 0.3 then they might, erroneously, still declare success with probability 0.021, and if the true VE<0.3 this type I error will be even less)

• When you say "we can estimate $\beta_1$", don't you think you'll get a biased estimated due to the fact the model may be misspecified due to omitted variables like propensity to wear mask, maintain social distance and others, or, in statistics parlor, don't you think having a shot is endogenous? Commented Aug 31, 2021 at 19:59
• @SergeyBushmanov the misspecified variables should be of equal influence in both the vaccine group and the placebo group. These variables do cause potentially some additional noise (not everyone is gonna be the same regarding wearing masks, social distance, etc.) but I do not see any strong reason to believe that this is gonna be a systematic error... Commented Sep 2, 2021 at 13:04
• ... there might be some criticism like the double-blind study being influenced by the difference in the effects from the vaccine and placebo (people recognizing the vaccine treatment), which would reduce the strength to the comparison of a placebo on the basis of health outcomes. But, I consider such effects as not so large. By now we also have a much larger amount of data available and we could safely say that the vaccine works very well in most of the population. The issue seems to be much more about heterogeneity now, discovering the part for which it not does not work (as well). Commented Sep 2, 2021 at 13:10
• the misspecified variables should be of equal influence in both the vaccine group and the placebo group It's not necessarily true. I'm asking this specific question having in mind countries like Russia, where (1) there is high opposition to getting a shot (2) it's more cautious (and more educated) people who get vaccinated more often. Surprisingly efficacy rates for Sputnik V reported abroad are even higher than in Russia. Commented Sep 2, 2021 at 13:20
• @SergeyBushmanov now you are talking about something different than the situation of a controlled study for which the $\beta_1$ had been estimated (without much bias). Sure, in other situations the effects of vaccination, or the risks among a vaccinated population versus not-vaccinated population, might be different. But for those other situations we talk about something different than the $\beta_1$... Commented Sep 3, 2021 at 11:39

All these results are consistent with using the conditional Maximum Likelihood Estimate as implemented in the base R implementation of the fisher's exact test:

splits <- matrix(c(6,26,15,47,25,67,35,85,53,111), ncol = 2, byrow = T)
total <- 43000

for(interim in 1:nrow(splits)) {
positive_vax <- splits[interim, 1]
positive_pla <- splits[interim, 2]
negative_vax <- (total / 2 ) - positive_vax
negative_pla <- (total / 2 ) - positive_pla

cont_tab <- matrix(c(positive_vax, positive_pla, negative_vax, negative_pla), nrow = 2)

test <- fisher.test(cont_tab)
VE <- 1 - test$estimate print(paste(VE, "% (", positive_vax, ":", positive_pla, ")")) }  Result: [1] "0.769425572629548 % ( 6 : 26 )" [1] "0.681342630733629 % ( 15 : 47 )" [1] "0.627606975573189 % ( 25 : 67 )" [1] "0.589208653283242 % ( 35 : 85 )" [1] "0.523803347975998 % ( 53 : 111 )"  • This is interesting, but could you add a few more details on the Fisher's exact test, or at least on why is it the right test to apply here? What hypothesis are we testing here? If you just substitute the 3 lines of code cont_tab <- matrix(c(positive_vax, positive_pla, negative_vax, negative_pla), nrow = 2), test <- fisher.test(cont_tab)and VE <- 1 - test$estimate with the single line VE <- 1 -positive_vax/positive_pla you get results that are even closer to those in the original table, so I'm not sure what using this Fisher's test is buying us here. Commented Nov 17, 2020 at 22:28
• I no longer think this is a good answer. The other answer is far better, perhaps I should delete this one?
– Dave
Commented Nov 18, 2020 at 9:34
• And note that Fisher's "exact" test is often not very accurate. Commented Nov 18, 2020 at 13:39