An analysis of ivermectin efficacy

ivermectin.R

library(tidyverse)
library(brms)
library(tidybayes)
# https://docs.google.com/spreadsheets/d/1vG0WdjZaYlS4_7_if-OaE3uMv3aX6zfvQMx5JvDREF0/edit?usp=sharing

# Chi.sq
prop.test(c(103, 135), c(5947, 5609),
          alternative = "less",
          conf.level = 0.99)

# Data
tibble(outcome = c(103, 135),
       n = c(5947, 5609),
       Group = factor(c("Treatment", "Control"))
) -> d

# Bayesian model to quantify predicted differences by group
brm(outcome | trials(n) ~ Group,
    data = d,
    iter = 2e4,
    family = "binomial") -> fit

d %>%
    select(Group) %>%
    mutate(n = 1e5) %>% # assume 100k treated in each group
    add_predicted_draws(fit) -> preds

preds %>%
    ungroup() %>%
    select(Group, .draw, .prediction) %>%
    pivot_wider(id_cols = .draw,
                names_from = Group,
                values_from = .prediction) %>%
    mutate(outcomes_prevented = Control - Treatment) %>%
    ggplot(aes(outcomes_prevented)) +
    geom_hline(yintercept = 0.5) +
    geom_vline(xintercept = 0) +
    stat_ecdf(size = 2) +
    scale_x_continuous(breaks = seq(-600, 1700, 1e2)) +
    ylab("Quantile") +
    theme_bw(14) +
    theme(panel.grid = element_blank()) +
    xlab("Outcomes prevented") +
    ggtitle("Predicted combined hospitalizations and deaths prevented by treating 100k w/ ivermectin")


preds %>%
    ungroup() %>%
    ggplot(aes(.prediction, fill = Group)) +
    geom_density(alpha = 0.4) +
    theme_bw(14) +
    theme(panel.grid = element_blank()) +
    geom_hline(yintercept = 0) +
    expand_limits(x = c(0)) +
    xlab("Predicted combined hospitalizations and deaths by treatment\n(assuming 100k treated in each group)") +
    ggtitle("Predicted combined hospitalizations and deaths by treatment\n(assuming 100k treated in each group)")

preds %>%
    ungroup() %>%
    group_by(Group) %>%
    summarise(mean(.prediction))