Pfizer's vaccine effectiveness simulation

Just for fun 😄. I saw this post about Pfizer's Vaccine Effectiveness Simulation. So I simply translate the Bayesian model (implemented in Stan) into my favorite Julia library Turing.jl. For details, please read the link.

Very briefly, from Vaccine Effectiveness Simulation

NYT reports a 44 thousand person trial with half of the people going to treatment and half to control. They further report that 162 people developed COVID in the control group and 8 where in the vaccine group. What is the probability that the vaccine is effective and what is the uncertainty in that probability? The Pfizer protocol defines vaccine effectiveness as follows:

using Random, Distributions, Turing
using RCall ; @rlibrary ggplot2 ; @rlibrary ggdist ; @rlibrary patchwork

Random.seed!(1)

# Model ------------------------------------------

# Numbers of volunteers, events in control, events in vaccine group, respectively
n, r_c, r_t = 4.4*10^4, 162, 8
n_c = n_t = n/2

@model PfizerVaccineEffectiveness(r_c, r_t) = begin
    # No prior beliefs for effectiveness
    p_c ~ Beta(1, 1)
    p_t ~ Beta(1, 1)

    r_c ~ Binomial(n_c, p_c)
    r_t ~ Binomial(n_t, p_t)

    # Generated quantities
    effect   = (p_t-p_c)*10^4  # Treatment effect
    ve       = 1 - p_t/p_c     # Vaccine effectiveness
    log_odds = log(p_t / (1-p_t)) - log(p_c / (1-p_c))

    return (effect, ve, log_odds)
end

N, n_chn, sampler_ = 10^4, 4, NUTS()

model = PfizerVaccineEffectiveness(r_c, r_t)
c = sample(model, sampler_, N)


# Plot -------------------------------------------

vals = generated_quantities(model, c)
effect, ve, log_odds = [getfield.(vals, i)[:] for i in 1:length(vals[1])]

p_effect = ggplot() + stat_halfeye(aes(effect)) +
  labs(title = "Reduction in infections on treatment per 10,000 people",
       x = "Size of effect", y = "Density")

p_log_odds = ggplot() + stat_halfeye(aes(log_odds)) +
  labs(title = """Log odds of treatment effect.
                  Negative means less likely to get infected on treatment""",
       x = "Log odds", y = "Density")

q_025, q_50, q_975 = round.(quantile(ve, [.025, .5, .975]) * 100, digits = 1)

p_ve = ggplot() + stat_halfeye(aes(ve)) +
  geom_vline(xintercept = q_50/100) +
  geom_text(aes(q_50/100, .5, label = "Median VE: $q_50%"), hjust = -.1) +
  labs(title = """Pfizer study protocol defines VE = 1 - Pt/Pc
               95% interval of effectiveness between $q_025% and $q_975%, compatible with model & data""",
       x = "Vaccine effectiveness", y = "Density")

wrap_plots(p_effect, p_log_odds, p_ve, ncol = 1)

kar9222/pfizer_vaccine_effectiveness_simulation.md

Pfizer's vaccine effectiveness simulation