simulating distribution of bootstrapped risk ratios https://x.com/nntaleb/status/1770385630163312902?s=20

risk-ratios-sim.md

library(tidyverse)

sim_risk_ratios <- function(x){
  events <- map2(rep(c(TRUE, FALSE), 5), c(31, 414 - 31, 82, 1492 - 82, 252, 4832 - 252, 423, 11831 - 423, 52, 1509-52), rep) %>% unlist()
  outcomes <- tibble(
    group = map2(c("<8 h", "8<10 h", "10<12 h", "12-16 h", ">16 h"), c(414, 1492, 4832, 11831, 1509), rep) %>% unlist()) %>% 
    mutate(event_sim = sample(events, n(), replace = TRUE)) %>% 
    group_by(group) %>% 
    summarise(risk = mean(event_sim))
  
  ref_risk <- filter(outcomes, group == "12-16 h") %>% pluck("risk")
  outcomes %>% 
    mutate(risk_ratio = risk / ref_risk)
}

simulations <- tibble(samp = 1:10000) %>% 
  mutate(sim_risk_ratios = map(samp, sim_risk_ratios)) %>% 
  unnest(sim_risk_ratios)

simulations %>% 
  arrange(desc(risk_ratio)) %>% 
  filter(group != "12-16 h") %>% 
  ggplot(aes(x = risk_ratio))+
  facet_wrap(~group)+
  geom_histogram(bins = 100)+
  theme_bw()

^{Created on 2024-03-25 with reprex v2.0.2}