nonpoisson.R

library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)
library(countrycode)
world <- readr::read_csv("https://population.un.org/wpp/Download/Files/1_Indicators%20(Standard)/CSV_FILES/WPP2019_TotalPopulationBySex.csv")
pops <- world %>% filter(Variant == "Medium") %>%
  mutate(country_code = countrycode(Location, origin = 'country.name', destination = 'iso3c')) %>%
  select(year=Time, country_code, PopTotal)
mort <- readr::read_csv("https://www.mortality.org/Public/STMF/Outputs/stmf.csv", skip = 1) %>%
  janitor::clean_names() 
joined <- mort %>% filter(sex == "b") %>%
  select(country_code, year, week, d_total) %>% 
  inner_join(pops, by=c("country_code", "year"))
step1 <- joined %>%
  filter() %>%
  arrange(country_code, week, year) %>%
  group_by(country_code,week) %>%
  mutate(rolling_d = d_total + lag(d_total),
         rolling_p = PopTotal + lag(PopTotal),
         lamb = PopTotal * rolling_d/rolling_p) %>%
  ungroup() %>%
  filter(!is.na(rolling_d)) %>%
  mutate(vs_pois = ppois(q = d_total, lambda = lamb)) %>% 
  arrange(country_code, year, week)
step1$abpoisson <- NA_real_
step1$abpoisson[104:nrow(step1)] <- sapply(104:nrow(step1), 
                                           function(x){sum(step1$vs_pois[(x-103):x] >= 0.975 |
                                                           step1$vs_pois[(x-103):x] < 0.025)/104})
step1 %>% group_by(country_code) %>%
  slice(104:n()) %>%
  ungroup() %>%
  mutate(Date = ISOdate(year, 1,1) + days(week*7-1)) %>%
  ggplot(aes(x=Date,y=abpoisson)) + geom_line() +
  facet_wrap(~ country_code, ncol=3) +
  ggtitle("How much do a country's weekly mortality rates resemble a poisson distribution")