kar9222 · December 31, 2020 03:22 · Nov 20, 2020 · Nov 20, 2020 · Nov 20, 2020 · Nov 20, 2020
diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -68,5 +68,5 @@ p_ve = ggplot() + stat_halfeye(aes(ve)) +
 wrap_plots(p_effect, p_log_odds, p_ve, ncol = 1)
 ```
 
-![Imgur](https://imgur.com/4tzl2BZ.png)
+![Imgur](https://i.imgur.com/OJHAqPS.png)
 
diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -62,7 +62,7 @@ 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 between $q_025% and $q_975%, compatible with model & data""",
+               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)

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -68,5 +68,5 @@ p_ve = ggplot() + stat_halfeye(aes(ve)) +
 wrap_plots(p_effect, p_log_odds, p_ve, ncol = 1)
 ```
 
-![Imgur](https://imgur.com/HIx0bBU.png)
+![Imgur](https://imgur.com/4tzl2BZ.png)
 
diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -56,13 +56,13 @@ p_log_odds = ggplot() + stat_halfeye(aes(log_odds)) +
                   Negative means less likely to get infected on treatment""",
        x = "Log odds", y = "Density")
 
-q_025, q_50, q_95 = round.(quantile(ve, [.025, .5, .95]) * 100, digits = 1)
+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 between $q_025% and $q_95%, compatible with model & data""",
+               95% interval 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)

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -56,13 +56,13 @@ p_log_odds = ggplot() + stat_halfeye(aes(log_odds)) +
                   Negative means less likely to get infected on treatment""",
        x = "Log odds", y = "Density")
 
-q_5, q_50, q_90 = round.(quantile(ve, [.025, .5, .95]) * 100, digits = 1)
+q_025, q_50, q_95 = round.(quantile(ve, [.025, .5, .95]) * 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 between $q_5% and $q_90%, compatible with model & data""",
+               95% interval between $q_025% and $q_95%, compatible with model & data""",
        x = "Vaccine effectiveness", y = "Density")
 
 wrap_plots(p_effect, p_log_odds, p_ve, ncol = 1)

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -9,7 +9,7 @@ Very briefly, from [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovi
 
 
 ```jl
-using Random, Distributions, Turing, DataFrames
+using Random, Distributions, Turing
 using RCall ; @rlibrary ggplot2 ; @rlibrary ggdist ; @rlibrary patchwork
 
 Random.seed!(1)

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -1,6 +1,6 @@
 # Pfizer's vaccine effectiveness simulation
 
-Just for fun. I saw this post about Pfizer's [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460). So I simply translate the Bayesian model (implemented in [Stan](https://mc-stan.org/)) into my favorite Julia library [Turing.jl](https://turing.ml/dev/). For details, please read the link.
+Just for fun 😄. I saw this post about Pfizer's [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460). So I simply translate the Bayesian model (implemented in [Stan](https://mc-stan.org/)) into my favorite Julia library [Turing.jl](https://turing.ml/dev/). For details, please read the link.
 
 Very briefly, from [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460)
 > 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:

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -1,4 +1,4 @@
-# Pfizer vaccine effectiveness simulation
+# Pfizer's vaccine effectiveness simulation
 
 Just for fun. I saw this post about Pfizer's [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460). So I simply translate the Bayesian model (implemented in [Stan](https://mc-stan.org/)) into my favorite Julia library [Turing.jl](https://turing.ml/dev/). For details, please read the link.
 

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -1,4 +1,4 @@
-# Pfizer Vaccine effectiveness simulation
+# Pfizer vaccine effectiveness simulation
 
 Just for fun. I saw this post about Pfizer's [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460). So I simply translate the Bayesian model (implemented in [Stan](https://mc-stan.org/)) into my favorite Julia library [Turing.jl](https://turing.ml/dev/). For details, please read the link.
 

diff --git a/pfizer_vaccine_effectiveness_simulation.md b/pfizer_vaccine_effectiveness_simulation.md
@@ -0,0 +1,72 @@
+# Pfizer Vaccine effectiveness simulation
+
+Just for fun. I saw this post about Pfizer's [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460). So I simply translate the Bayesian model (implemented in [Stan](https://mc-stan.org/)) into my favorite Julia library [Turing.jl](https://turing.ml/dev/). For details, please read the link.
+
+Very briefly, from [Vaccine Effectiveness Simulation](https://rpubs.com/ericnovik/692460)
+> 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:
+
+![equation](https://latex.codecogs.com/gif.latex?%5Ctext%7BVE%7D%20%3D%201%20-%20%5Cfrac%7Bp_%7Bt%7D%7D%7Bp_%7Bc%7D%7D)
+
+
+```jl
+using Random, Distributions, Turing, DataFrames
+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_5, q_50, q_90 = round.(quantile(ve, [.025, .5, .95]) * 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 between $q_5% and $q_90%, compatible with model & data""",
+       x = "Vaccine effectiveness", y = "Density")
+
+wrap_plots(p_effect, p_log_odds, p_ve, ncol = 1)
+```
+
+![Imgur](https://imgur.com/HIx0bBU.png)
+