# R analysis from pfizer vaccine paper and slides

# multinational, placebo-controlled and
# observed-blinded trial

n <- 43448   # total sample size
n_vax <- 21720  # sample size of vaccinated
n_placebo <- 21728  # sample size of placebo
cases_vax <- 8  # cases with vaccine
cases_placebo <- 162  # cases with placebo
p_vax <- cases_vax/n_vax # proportion of confirmed covid-19 among vaccinated
p_placebo <- cases_placebo/n_placebo # proportion of confirmed covid-19 among placebo
IRR <- p_vax/p_placebo  # incidence rate ratio: treated divided by control incidence rates
VE <- 1-IRR # vaccine efficacy
theta <- (1-VE)/(2-VE)  # parameter: probability of being vaccinated given being infected

# use of a Bayesian beta-binomial model
# to compute vaccine efficacy

# we need a prior on theta: Beta(0.700102,1):
# slide 19.
# used to guarantee a threshold of prior vaccine efficacy
# not lower than 30%
alpha_prior <- 0.700102
beta_prior <- 1

# expected prior values 
expected_prior_theta <- alpha_prior/(alpha_prior+beta_prior)
expected_prior_VE <- (1-2*expected_prior_theta)/(1-expected_prior_theta)

# 95% credible prior interval of vaccine efficacy
interval_prior_theta <- qbeta(c(0.025, 0.975), alpha_prior, beta_prior)
interval_prior_VE <- (1-2*interval_prior_theta[c(2,1)])/(1-interval_prior_theta[c(2,1)])

# Then we may compute the posterior as:
# prior X likelihood: Beta(alpha_prior + cases_vax, beta_prior + cases_placebo)
# slides: 20-21
# we define the posterior parameters:

alpha_posterior <- alpha_prior + cases_vax
beta_posterior <- beta_prior + cases_placebo

# we compute the expected posterior values:
expected_posterior_theta <- alpha_posterior/(alpha_posterior+beta_posterior)
expected_posterior_VE <- (1-2*expected_posterior_theta)/(1-expected_posterior_theta)

# 95% credible prior interval of vaccine efficacy
interval_posterior_theta <- qbeta(c(0.025, 0.975), alpha_posterior, beta_posterior)
interval_posterior_VE <- (1-2*interval_posterior_theta[c(2,1)])/(1-interval_posterior_theta[c(2,1)])

# Then we may compute final probabilities: such as
# the probability that the VE is greater than the 30%,
# which is approximately 1

# posterior distribution for theta
curve(dbeta(x, alpha_posterior, beta_posterior),
      xlab = expression(theta), ylab = "beta distribution",
      main ="Probability of being vaccinated given infected",
      cex.main=0.6)
abline(v = expected_posterior_theta, col = "blue")
abline(v=interval_posterior_theta[1], col ="blue", lty =2)
abline(v=interval_posterior_theta[2], col = "blue", lty =2)
abline(v=expected_prior_theta, col="red")

# transform in VE
new_function <- function(x){
 dbeta((1-x)/(2-x), alpha_posterior, beta_posterior) *abs( (-(2-x)+(1-x))/  ((2-x)^2))
}

curve(new_function(x), xlab = "VE", 
      ylab = "Distribution",
      main = "Distribution of vaccine efficacy",
      cex.main =0.6)
abline(v = expected_posterior_VE, col = "blue")
abline(v=interval_posterior_VE[1], col ="blue", lty =2)
abline(v=interval_posterior_VE[2], col = "blue", lty =2)
abline(v = expected_prior_VE, col = "red")

# posterior distribution for VE greater than 30%
# slide: 25
# Use the threshold: 0.986

threshold <- 0.986
ve_eff_threshold <- 0.3

# compute this threshold for the original data:
# the probability that theta is lower than 0.4118 is 1,
# which is greater than 0.986
pbeta(expected_prior_theta, alpha_posterior, beta_posterior,
      lower.tail = TRUE) > threshold

# I accept the alternative hypothesis, success.

# what if I have other number of cases?
# For instance, suppose to have the new numbers:
new_cases_vax <- 53
new_cases_placebo <- 111

# I redefine the posterior parameters
new_alpha_posterior <- alpha_prior + new_cases_vax
new_beta_posterior <- beta_prior + new_cases_placebo

# I recompute:

pbeta(expected_prior_theta, new_alpha_posterior,
      new_beta_posterior,
      lower.tail = TRUE) > threshold

# which is 0.99 > 0.986: still accepting the
# alternative hypothesis of success!

# what happens if new_cases_vax = 54?
new_cases_vax <- 54
new_cases_placebo <- 110

# I redefine the posterior parameters
new_alpha_posterior <- alpha_prior + new_cases_vax
new_beta_posterior <- beta_prior + new_cases_placebo

# I recompute:

pbeta(expected_prior_theta, new_alpha_posterior,
      new_beta_posterior,
      lower.tail = TRUE) > threshold

# 0.9853 < 0.986, so in this case I accept 
# the null hypothesis, the non-success.


# So the experimenters fixed some cases threshold
# to verify the efficacy during the trials.


# more r code and comment here:
# http://skranz.github.io//r/2020/11/11/CovidVaccineBayesian.html?utm_source=feedburner&utm_medium=feed&utm_campaign=Feed%3A+skranz_R+%28Economics+and+R+%28R+Posts%29%29
# https://statmodeling.stat.columbia.edu/2020/11/13/pfizer-beta-prior-vaccine-effect/