Storyboard for GetTBinR 2017 data

storyboard-gettbinr-2017

## Get required packages - managed using pacman
if (!require(pacman)) install.packages("pacman"); library(pacman)
p_load("getTBinR")
p_load("ggplot2")
p_load("ggrepel")
p_load("scales")
p_load("viridis")
p_load("dplyr")
p_load("tidyr")
p_load("forcats")
p_load("ggridges")
p_load_gh("thomasp85/patchwork")

##Pull TB data 
tb <- get_tb_burden(verbose = FALSE) 


## Summarise global changes
global_tb <- summarise_tb_burden(compare_to_world = TRUE,
                                 annual_change = TRUE,
                                 stat = "rate",
                                 verbose = FALSE) %>% 
  filter(area == "Global")


## TB in 2017
tb_2017 <- global_tb %>% 
  filter(year == 2017)

## Global annual change
global_annual_change <- ggplot(global_tb, aes(year, e_inc_num)) +
  geom_smooth(se = FALSE, col = "black", size = 1.2, alpha = 0.7) +
  geom_point(size = 1.2, alpha = 0.8, col = "black") +
  scale_y_continuous(label = scales::percent, minor_breaks = NULL, breaks = seq(-0.025, 0, 0.0025)) +
  theme_minimal() +
  labs(
    y = "Annual Percentage Change",
    x = "Year",
    title = "Global Annual Percentage Change in Tuberculosis Incidence Rates",
    caption = ""
  )

## Remove countries with incidence below 1000 or incidence rates below 10 per 100,000 to reduce noise and cal country level annual change.
countries_with_tb_burden <- tb_burden %>% 
  filter(year == 2017,
         e_inc_100k > 10,
         e_inc_num > 1000) %>% 
  pull(country)

tb_annual_change <- summarise_tb_burden(countries = countries_with_tb_burden, 
                                        compare_to_region = FALSE, 
                                        compare_to_world = FALSE,
                                        metric = "e_inc_100k",
                                        annual_change = TRUE,
                                        verbose = FALSE) %>% 
  mutate(annual_change = e_inc_100k) %>% 
  left_join(tb_burden %>% 
              select(country, g_whoregion), 
                     by = c("area" = "country")) %>% 
  drop_na(g_whoregion)

## Function to plot annual change
plot_annual_change <- function(df, strat = NULL, subtitle = NULL, years = 2000:2017) {
  dist <- df %>% 
    filter(year %in% years) %>% 
    rename(Region = g_whoregion) %>% 
    mutate(year = year %>% 
             factor(ordered = TRUE) %>% 
             fct_rev) %>% 
    ggplot(aes_string(x = "annual_change", y = "year", col = strat, fill = strat)) +
    geom_density_ridges(quantile_lines = TRUE, quantiles = 2, alpha = 0.6) +
    scale_color_viridis(discrete = TRUE, end = 0.9) +
    scale_fill_viridis(discrete = TRUE, end = 0.9) +
    geom_vline(xintercept = 0, linetype = 2, alpha = 0.6) +
    scale_x_continuous(labels = scales::percent, breaks = seq(-0.4, 0.4, 0.1),
                       limits = c(-0.4, 0.4), minor_breaks = NULL) +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(x = paste0("Annual Change in ", search_data_dict("e_inc_100k")$definition),
         y = "Year",
         title = "Annual Percentage Change in Tuberculosis Incidence Rates",
         subtitle = subtitle,
         caption = "")
  
  return(dist)
}

## Overall country level annual change
overall <- plot_annual_change(tb_annual_change, NULL,
                              years = seq(2001, 2017, 2), subtitle = "By Country") 


## Regional country level annual change
region <-  plot_annual_change(tb_annual_change, "Region",
                              subtitle = "By Region", 
                              years = seq(2001, 2017, 2)) + 
  facet_wrap(~Region) +
  labs(caption = "")


## Regional and Global TB incidence rates over time
regional_incidence <- plot_tb_burden_summary(conf = NULL)

## Map global TB incidence rates for 2017 using getTBinR
map <- map_tb_burden(year = c(2005, 2009, 2013, 2017), facet = "year") +
  theme(strip.background = element_blank()) +
  labs(caption = "",
       title = "Tuberculosis Incidence Rates",
       subtitle = "By Country, per 100,000 population")

## Compose storyboard
storyboard <- (map + regional_incidence + plot_layout(widths = c(2, 1))) /
  (region + (global_annual_change /
               overall + labs(caption = "For country level annual percentages change countries with incidence above 1000 and an incidence rate above 10 per 100,000 are shown.
                    The global annual percentage change is shown with a LOESS fit. 
                    By @seabbs | Made with getTBinR | Source: World Health Organisation")) + plot_layout(widths = c(2, 1))) +
  plot_layout(widths = c(1, 1))

## Save storyboard
ggsave("../../static/img/getTBinR/storyboard-6-0.png",
       storyboard, width = 20, height = 15, dpi = 330)