Population Data from Infoshare (quarterly, as at, DPE). Deaths from StatsNZ open data API, but also in their Covid data portal

excess_mort_2022.R

legpos="right"
theme_davidhood <- function(){
  # theme_minimal(base_family="OpenSans") %+replace% 
  theme_minimal() %+replace% 
    theme(panel.grid = element_blank(),
          axis.line.x = element_line(size=0.1),
          axis.line.y = element_line(size=0.1),
          axis.ticks = element_line(size=0.2),
          strip.background = element_rect(fill= "#FFFFFF", colour="#EFEFEF"),
          strip.text = element_text(size = 13,
                                    margin = margin(t = 5, r = 5, b = 5, l = 5, unit = "pt")),
          strip.placement = "inside",
          panel.background = element_rect(fill = "#FFFFFF", colour = "#FFFFFF"),
          panel.spacing = unit(1.5, "lines"),
          plot.title = element_text(size = 14,
                                    lineheight = 1.23,
                                    margin=margin(t = 0, r = 0, b = 10, l = 10, unit = "pt"),
                                    hjust=0),
          plot.background = element_rect(fill = "#F3F3F3"),
          axis.title = element_text(size=13),
          plot.caption = element_text(margin=margin(t = 5, r = 5, b = 5, l = 5, unit = "pt"),
                                      size=11, hjust=1),
          plot.caption.position = "plot",
          plot.margin = margin(0.3, 0.4, 0.3, 0.3, "cm"),
          legend.position = legpos)
}


source("get-odata-fun.R")

library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(ggthemes)
library(readr)

six_cols <- (colorblind_pal()(6))


if(!file.exists("../misc_data/NZ_weekly_deaths.csv")){
theres_data <- get_odata(
  service = "https://api.stats.govt.nz/opendata/v1",
  endpoint = "Covid-19Indicators",
  entity = "Observations",
  query_option = "$filter=(ResourceID eq 'CPWEE1')",
  service_api_key = "You_will_need_your_own_API_key_from_StatsNZ")
write.csv(theres_data, "../misc_data/NZ_weekly_deaths.csv")
}

weekly_deaths <- read_csv("../misc_data/NZ_weekly_deaths.csv",
                          col_types=cols(
                            ...1 = col_double(),
                            id = col_character(),
                            ResourceID = col_character(),
                            GeoUnit = col_logical(),
                            Geo = col_logical(),
                            Period = col_date(format = ""),
                            Duration = col_character(),
                            Label1 = col_character(),
                            Label2 = col_logical(),
                            Label3 = col_logical(),
                            Label4 = col_logical(),
                            Label5 = col_logical(),
                            Label6 = col_logical(),
                            Value = col_double(),
                            Unit = col_character(),
                            Measure = col_character(),
                            Multiplier = col_double(),
                            NullReason = col_logical(),
                            Status = col_logical()
                          ), name_repair = "unique")
latest_week = max(weekly_deaths$Period)
age_populations <- read_csv("DPE403901_20220311_082941_21.csv", skip=3,
                            col_types = cols(
                              ...1 = col_character(),
                              `0-4 Years` = col_double(),
                              `5-9 Years` = col_double(),
                              `10-14 Years` = col_double(),
                              `15-19 Years` = col_double(),
                              `20-24 Years` = col_double(),
                              `25-29 Years` = col_double(),
                              `30-34 Years` = col_double(),
                              `35-39 Years` = col_double(),
                              `40-44 Years` = col_double(),
                              `45-49 Years` = col_double(),
                              `50-54 Years` = col_double(),
                              `55-59 Years` = col_double(),
                              `60 Years and Over` = col_double(),
                              `80 Years and Over` = col_double()
                            ), name_repair = "unique")

age_populations$Total = rowSums(age_populations[,2:14])
age_populations$`80 and over` = age_populations$`80 Years and Over`
age_populations$`60 to 79` = age_populations$`60 Years and Over` - 
  age_populations$`80 Years and Over`
age_populations$`30 to 59` = age_populations$`30-34 Years` +
  age_populations$`35-39 Years` +
  age_populations$`40-44 Years` + 
  age_populations$`45-49 Years` + 
  age_populations$`50-54 Years` + 
  age_populations$`55-59 Years`
age_populations$`Under 30` = age_populations$`0-4 Years` +
  age_populations$`5-9 Years` +
  age_populations$`10-14 Years` + 
  age_populations$`15-19 Years` + 
  age_populations$`20-24 Years` + 
  age_populations$`25-29 Years`

quartly_pop <- age_populations %>%
  filter(!is.na(Total)) %>%
  separate(`...1`, into=c("Yr","Qt"), sep="Q", convert = TRUE) %>%
  mutate(Quarterly_date = floor_date(ISOdate(Yr, Qt*3,30) + days(3),"quarter")) %>%
  select(Quarterly_date, Total:`Under 30`) %>%
  gather(key="Label1", value="Population", Total:`Under 30`)

nextq <- quartly_pop %>%
  filter(Quarterly_date > max(Quarterly_date) - days(200)) %>%
  arrange(Label1, Quarterly_date) %>%
  group_by(Label1) %>%
  mutate(change = Population - lag(Population),
         meanchange = mean(change, na.rm=TRUE),
         Population = Population + meanchange) %>%
  slice(n()) %>%
  ungroup() %>%
  mutate(Quarterly_date = floor_date(Quarterly_date + days(120), unit="quarter")) %>%
  select(Label1, Quarterly_date, Population)

quartly_pop <- bind_rows(quartly_pop, nextq)
mortality <- weekly_deaths %>%
  select(Period, Label1, Value) %>%
  mutate(Quarterly_date = floor_date(Period, "quarter")) %>%
  inner_join(quartly_pop, by = c("Label1", "Quarterly_date")) %>% 
  mutate(mort100K = 100000 * Value/Population,
         week_ending = ISOdate(2022, month(Period), day(Period)),
         Years = ifelse(year(Period) == 2022, "2022", "2011-21")) %>%
  filter(Label1 != "Total") %>%
  mutate(Label1 = factor(Label1, 
                         levels=c("Under 30", "30 to 59",
                                  "60 to 79", "80 and over")))

lineonly <- mortality %>% filter(year(Period) == 2022)

graph_mort <- ggplot(mortality,
                     aes(x=week_ending, y=mort100K, colour=Years, alpha=Years)) +
  geom_point(size=0.8) + facet_wrap(~ Label1, ncol=2, scales = "free_y")+
  geom_line(data=lineonly) +
  theme_davidhood() + 
  scale_alpha_manual(values=c(0.6,1)) +
  scale_colour_manual(values=six_cols[c(2,1)]) +
  scale_x_datetime(date_labels = "%b") +
  labs(title=paste("NZ Mortality by week. To week ending", latest_week),
       subtitle="2022 (black) compared to 2011-2021 aggregate pattern (orange)",
       x="\nWeek ending\n", y="\nDeaths per 100,000\n",
       caption="Source:StatsNZ Open Data weekly deaths, Infoshare population")
graph_mort
ggsave(filename="~/Desktop/mort_22.png", plot=graph_mort,dpi=72, units="in", 
       bg="white", height = 5.556 * 1.6, width=9.877* 1.6)