COVID-19 Hospitalization Trends

This analysis compares the United States National and Regional trends of child and adult COVID-19 hospitalizations.

Data was acquired through the United States Department of Health and Human Services (HHS)

Helpful FAQ Resource

library(tidyverse)
library(ggplot2)
library(formattable)
library(data.table)
library(DT)

Data Import

Import HHS Hospitalization Counts and Census Data

hhs_data <- read_csv('data/COVID-19_Reported_Patient_Impact_and_Hospital_Capacity_by_State_Timeseries_20210417.csv')

census_data <- read.csv('data/sc-est2020-18+pop-res.csv')

state_abv <- read.delim2('data/state_abv.txt', sep = "|")

Data Pre-Processing

Add full state names to HHS Data

state_abv <- state_abv %>% 
  mutate(State = tolower(STATE_NAME)) %>%
  rename("state_abv" = STUSAB) %>%
  select(state_abv, State)

hhs_data <- hhs_data %>%
  rename("state_abv" = state) %>%
  left_join(., state_abv, by = "state_abv")

Remove US.Territories

hhs_data <- hhs_data %>%
  filter(!state_abv %in% c("GU", "MP", "PR", "VI"))

Clean US Census Data

Calculate new variable POPEST17uNDER2020 for estimated population size of ages \(\leq\) 17.

census_data <- census_data %>% 
  mutate(NAME = tolower(NAME)) %>%
  rename('State' = NAME) %>% 
  mutate(POPEST17Under2020 = POPESTIMATE2020 - POPEST18PLUS2020) %>%
  rename('Region' = 'REGION')

Recode Regions according to the U.S Census

census_data$Region <- recode(census_data$Region, 
                             '0' = 'US', 
                             '1' = 'Northeast', 
                             '2' = 'Midwest', 
                             '3' = 'South', 
                             '4' = 'West')

Add Census data to HHS data

hhs_data <- left_join(hhs_data, census_data, by = 'State')

Convert to date format

hhs_data$date <- as.Date(hhs_data$date)

Remove dates with no data collection

We will remove dates when there was noprevious_day_admission_pediatric_covid_confirmedcases reported for the whole week. Date collection appears to begin on July 15, 2020. However, previous structuring of the HHS had the data collection_week beginning on July 30 (see more note below). To stay consistent, we will remove data prior to July 30th

#View(hhs_data %>% select(date, previous_day_admission_pediatric_covid_confirmed) %>% group_by(date) %>% mutate(sum_date = sum(previous_day_admission_pediatric_covid_confirmed, na.rm = TRUE)) %>% distinct(date, sum_date))

hhs_data <- hhs_data %>% filter(!date < '2020-07-31')

Sum by month

Similar to the original HHS data release, we will sum by month. As described by the original data dictionary:

“For a given entry, the term collection_week signifies the start of the period that is aggregated. For example, a “collection_week” of 2020-11-20 means the average/sum/coverage of the elements captured from that given facility starting and including Friday, November 20, 2020, and ending and including reports for Thursday, November 26, 2020."

To develop a new collection_week variable, we will download the facility level HHS data and map each collection_week to a date_id. We will then group each unique data in the newer version of the data into bins of 7 days that will match the previous old collection_week bins.

facility_level <- read_csv("data/COVID-19_Reported_Patient_Impact_and_Hospital_Capacity_by_Facility_20210329.csv")

facility_level <- facility_level %>% 
  distinct(collection_week) %>%
  mutate(collection_week = as.Date(collection_week)) %>%
  data.frame() %>%
  arrange(collection_week)

# identify unique dates
date_week_id <- facility_level %>% distinct(collection_week) %>% arrange(collection_week)
date_week_id$date_id <- seq.int(nrow(date_week_id))
date_week_id <- date_week_id %>%
  rename("date" = collection_week)

# function to bin all dates every 7 days
all_dates <- hhs_data %>% distinct(date) %>% arrange(date)
all_dates$date_id <- c(0, rep(1:(nrow(all_dates)-1)%/%7))+1

# add the date_ids to main dataset and create new column "collection_week" with the first date for each date_id bin
hhs_data <- hhs_data %>%
  left_join(., all_dates, by = "date") %>% 
  group_by(date_id) %>% 
  mutate(collection_week = min(date)) %>%
  ungroup() %>%
  filter(!collection_week >= "2021-04-16") # remove the latest week due to incomplete data collection at the initial time of running this script

Data quality

In rare cases, if the previous_day_admission column is a negative number, we will replace with a 0

hhs_data <- hhs_data %>%
  mutate(previous_day_admission_pediatric_covid_confirmed = if_else(previous_day_admission_pediatric_covid_confirmed < 0, 0,
                                                                    previous_day_admission_pediatric_covid_confirmed),
         previous_day_admission_adult_covid_confirmed = if_else(previous_day_admission_adult_covid_confirmed < 0, 0,
                                                                previous_day_admission_adult_covid_confirmed))

Select columns of interest

hhs_data <- hhs_data %>% 
  select(collection_week, state_abv, State, Region, POPEST18PLUS2020, POPEST17Under2020,
         previous_day_admission_pediatric_covid_confirmed, previous_day_admission_adult_covid_confirmed)

Calculate Hospitalization Rates

First, we will straify into pediatric and adult datasets

hhs_pediatric_processed <- hhs_data %>%
  group_by(collection_week) %>%
  mutate(hospitalizations = sum(previous_day_admission_pediatric_covid_confirmed, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(collection_week, Region) %>%
  mutate(hospitalizations_region = sum(previous_day_admission_pediatric_covid_confirmed, na.rm = TRUE)) %>%
  ungroup() %>%
  select(collection_week, state_abv, State, Region, POPEST18PLUS2020, POPEST17Under2020, 
         hospitalizations, hospitalizations_region) %>%
  mutate(population = "pediatric")

hhs_adult_processed <- hhs_data %>%
  group_by(collection_week) %>%
  mutate(hospitalizations = sum(previous_day_admission_adult_covid_confirmed, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(collection_week, Region) %>%
  mutate(hospitalizations_region = sum(previous_day_admission_adult_covid_confirmed, na.rm = TRUE)) %>%
  ungroup() %>%
  select(collection_week, state_abv, State, Region, POPEST18PLUS2020, POPEST17Under2020,
         hospitalizations, hospitalizations_region) %>%
  mutate(population = "adult")

Combine adult and pediatric hospitalizations

hhs_data_proc <- rbind(hhs_pediatric_processed,
                       hhs_adult_processed)

Standardize Hospitalization Rates

Here we will standardize hospitalizations per 100,000 adults/children throughout the Nation and by Region

pop_adult_nation <- census_data %>% 
  filter(State == "united states") %>% 
  select(POPEST18PLUS2020) %>% 
  as.numeric()

pop_pediatric_nation <- census_data %>% 
  filter(State == "united states") %>% 
  select(POPEST17Under2020) %>% 
  as.numeric()

region_estimates <- hhs_data_proc %>% 
  ungroup() %>%
  distinct(Region, POPEST17Under2020, POPEST18PLUS2020) %>%
  group_by(Region) %>%
  mutate(Region_census_pediatric = sum(POPEST17Under2020),
         Region_census_adult = sum(POPEST18PLUS2020)) %>%
  select(Region, Region_census_pediatric, Region_census_adult) %>%
  distinct()


hhs_data_proc <- hhs_data_proc %>%
  left_join(., region_estimates, by = "Region") %>%
  mutate(hospitalizations_standardized_us = if_else(population == "pediatric", 
                                                    hospitalizations/pop_pediatric_nation*100000,
         hospitalizations/pop_adult_nation*100000)) %>%
  group_by(collection_week, population, Region) %>%
  mutate(hospitalizations_standardized_region = if_else(population == "pediatric", 
                                                    hospitalizations_region/Region_census_pediatric*100000,
         hospitalizations_region/Region_census_adult*100000)) %>%
  ungroup()

Statistics

Total Number of Pediatric and Adult Hospitalizations over data collection

hosp_counts <- hhs_data_proc %>% 
  distinct(population, collection_week, hospitalizations) %>%
  group_by(population) %>%
  mutate(`Total Hospitalizations` = sum(hospitalizations, na.rm = TRUE)) %>%
  ungroup() %>%
  distinct(population, `Total Hospitalizations`)

datatable(hosp_counts)

National Median Weekly Hospitalization Rate

hosp_counts_med <- hhs_data_proc %>% 
  group_by(population) %>%
  distinct(population, collection_week, hospitalizations) %>%
  mutate(Q1 = round(quantile(hospitalizations, 0.25), 1),
         `Median Weekly Hospitalizations` = median(hospitalizations, na.rm = TRUE),
         Q3 = round(quantile(hospitalizations, 0.75), 1),) %>%
  distinct(population, Q1, `Median Weekly Hospitalizations`, Q3)

datatable(hosp_counts_med)

National Median Hospitalization Rate per 100,000 Children/Adults

national_stats <- hhs_data_proc %>% 
  distinct(population, hospitalizations_standardized_us) %>%
  group_by(population) %>%
  mutate(Q1 = round(quantile(hospitalizations_standardized_us, 0.25), 1),
         Median = round(median(hospitalizations_standardized_us), 1),
         Q3 = round(quantile(hospitalizations_standardized_us, 0.75), 1)) %>%
  select(-hospitalizations_standardized_us) %>%
  distinct() %>%
  ungroup() %>%
  as.data.frame()

datatable(national_stats)

National Hospitalization Trends

ggplot(hhs_data_proc %>% 
         mutate(hospitalizations_standardized_us = if_else(population == "pediatric", hospitalizations_standardized_us*10, hospitalizations_standardized_us)) %>% 
       distinct(collection_week, hospitalizations_standardized_us, population),
       aes(x = collection_week, y = hospitalizations_standardized_us, color = population)) + 
  geom_line() + 
  geom_point() +
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  ylab("Hospitalizations per 100K Adults") + 
  xlab("") + 
  scale_x_date(date_labels = "%b-%d", date_breaks = "1 month") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./10, name = "Hospitalizations per 100K Children")) + 
  theme_bw() + 
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 15),
        legend.text = element_text(size=15),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 10),
        axis.title.y = element_text(face = "bold", size = 13, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 13, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 13, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 13, color='#D95F02'))

National Hospitalizations Rates

Rates are expressed as per 100,000 adults or children

nat_table <- hhs_data_proc %>% 
  distinct(collection_week, hospitalizations_standardized_us, population) %>% 
  mutate(hospitalizations_standardized_us = round(hospitalizations_standardized_us, 2)) %>%
  arrange(collection_week) %>%
  pivot_wider(names_from = population, values_from = hospitalizations_standardized_us) %>%
  rename("Week" = collection_week, 
         `Child Hospitalizations` = pediatric, 
         `Adult Hospitalizations` = adult)

lowValue = "#ffffff"
customOrange = "#D95F02"
customGreen = "#1B9E77"

as.datatable(
  formattable(nat_table, align = c('c', 'c','c'),  list(
            `Indicator Name` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")),
            `Child Hospitalizations`= color_tile(lowValue, customOrange),
            `Adult Hospitalizations`= color_tile(lowValue, customGreen)))
  )

Regional Hospitalization Trends

Regional Median Hospitalization Rates per 100,000 Children/Adults

regional_stats <- hhs_data_proc %>% 
  distinct(population, Region, hospitalizations_standardized_region) %>%
  group_by(population, Region) %>%
  mutate(Q1 = round(quantile(hospitalizations_standardized_region, 0.25), 1),
         Median = round(median(hospitalizations_standardized_region), 1),
         Q3 = round(quantile(hospitalizations_standardized_region, 0.75), 1)) %>%
  select(-hospitalizations_standardized_region) %>%
  distinct() %>%
  arrange(Region) %>%
  ungroup() %>%
  as.data.frame()

datatable(regional_stats)

ggplot(hhs_data_proc %>%
         mutate(hospitalizations_standardized_region = if_else(population == "pediatric", 
                                                 hospitalizations_standardized_region*10, hospitalizations_standardized_region)) %>%
         distinct(collection_week, Region, population, hospitalizations_standardized_region),
       aes(x = collection_week, y = hospitalizations_standardized_region, color = population)) + 
  geom_point() +
  geom_line(alpha=0.9) + 
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  ylab("Hospitalizations per 100k Adults") + 
  xlab("") + 
  scale_x_date(date_labels = "%b", date_breaks = "1 month") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./10, name = "Hospitalizations per 100k Children")) + 
  theme_bw() +
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 15),
        legend.text = element_text(size=15),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 10),
        axis.title.y = element_text(face = "bold", size = 13, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 13, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 13, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 13, color='#D95F02')) + 
  facet_wrap(~Region, scales = "fixed") + 
  theme(strip.text.x = element_text(size = 18))

Hospitalizations Rates by Region

Rates are expressed as per 100,000 adults or children

lowValue = "#ffffff"
customOrange = "#D95F02"
customGreen = "#1B9E77"

region_table <- hhs_data_proc %>% 
  distinct(collection_week, Region, hospitalizations_standardized_region, population) %>% 
  mutate(hospitalizations_standardized_region = round(hospitalizations_standardized_region, 2)) %>%
  arrange(collection_week) %>%
  pivot_wider(names_from = population, values_from = hospitalizations_standardized_region) %>%
  rename("Week" = collection_week, 
         `Child Hospitalizations` = pediatric, 
         `Adult Hospitalizations` = adult)

Midwest

as.datatable(
  formattable(region_table %>% filter(Region == "Midwest") %>% select(-Region), align = c('c', 'c','c'),  list(
            `Indicator Name` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")),
            `Child Hospitalizations`= color_tile(lowValue, customOrange),
            `Adult Hospitalizations`= color_tile(lowValue, customGreen)))
  )

Northeast

as.datatable(
  formattable(region_table %>% filter(Region == "Northeast") %>% select(-Region), align = c('c', 'c','c'),  list(
            `Indicator Name` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")),
            `Child Hospitalizations`= color_tile(lowValue, customOrange),
            `Adult Hospitalizations`= color_tile(lowValue, customGreen)))
  )

South

as.datatable(
  formattable(region_table %>% filter(Region == "South") %>% select(-Region), align = c('c', 'c','c'),  list(
            `Indicator Name` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")),
            `Child Hospitalizations`= color_tile(lowValue, customOrange),
            `Adult Hospitalizations`= color_tile(lowValue, customGreen)))
  )

West

as.datatable(
  formattable(region_table %>% filter(Region == "West") %>% select(-Region), align = c('c', 'c','c'),  list(
            `Indicator Name` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")),
            `Child Hospitalizations`= color_tile(lowValue, customOrange),
            `Adult Hospitalizations`= color_tile(lowValue, customGreen)))
  )

Change Point Analysis

source("R/change_point_function.R")

## modify the HHS datast
date_week_id <- hhs_data_proc %>% distinct(collection_week) %>% arrange(collection_week)
date_week_id$date_id <- seq.int(nrow(date_week_id))

hhs_data_proc <- hhs_data_proc %>% left_join(., date_week_id, by = "collection_week")

National Hospitalizations Change Point

set.seed(101)

## adult hosps
adult_hosps <- hhs_data_proc %>% 
  ungroup() %>% 
  filter(population == "adult") %>%
  select(collection_week, hospitalizations_standardized_us) %>% 
  arrange(collection_week) %>%
  distinct(hospitalizations_standardized_us)

adult_hosps_mat <- t(data.matrix(adult_hosps))
row.names(adult_hosps_mat) <- NULL

## pediatric hosps
ped_hosps <- hhs_data_proc %>% 
  ungroup() %>% 
  filter(population == "pediatric") %>%
  select(collection_week, hospitalizations_standardized_us) %>% 
  arrange(collection_week) %>%
  distinct(hospitalizations_standardized_us)

ped_hosps_mat <- t(data.matrix(ped_hosps))
row.names(ped_hosps_mat) <- NULL


cp_nat_adult <- change_point_test(adult_hosps_mat, w = rep(1 / nrow(adult_hosps_mat), nrow(adult_hosps_mat)), 
                              t_range = 2:(ncol(adult_hosps_mat) - 2), boot_num = 10000)


cp_nat_ped <- change_point_test(ped_hosps_mat, w = rep(1 / nrow(ped_hosps_mat), nrow(ped_hosps_mat)), 
                              t_range = 2:(ncol(ped_hosps_mat) - 2), boot_num = 10000)

Hospitalizations Change Point by Region

Midwest

midwest_adult_hosps <- hhs_data_proc %>% 
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "adult" & Region == "Midwest") %>%
  distinct(hospitalizations_standardized_region) 

midwest_adult_hosps_mat <- t(data.matrix(midwest_adult_hosps))
row.names(midwest_adult_hosps_mat) <- NULL

cp_reg_midwest_adult <- change_point_test(midwest_adult_hosps_mat, w = rep(1 / nrow(midwest_adult_hosps_mat), nrow(midwest_adult_hosps_mat)), 
                              t_range = 2:(ncol(midwest_adult_hosps_mat) - 2), boot_num = 10000)

midwest_ped_hosps <- hhs_data_proc %>%
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "pediatric" & Region == "Midwest") %>%
  distinct(hospitalizations_standardized_region) 

midwest_ped_hosps_mat <- t(data.matrix(midwest_ped_hosps))
row.names(midwest_ped_hosps_mat) <- NULL


cp_reg_midwest_ped <- change_point_test(midwest_ped_hosps_mat, w = rep(1 / nrow(midwest_ped_hosps_mat), nrow(midwest_ped_hosps_mat)), 
                              t_range = 2:(ncol(midwest_ped_hosps_mat) - 2), boot_num = 10000)

Northeast

northeast_adult_hosps <- hhs_data_proc %>% 
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "adult" & Region == "Northeast") %>%
  distinct(hospitalizations_standardized_region) 

northeast_adult_hosps_mat <- t(data.matrix(northeast_adult_hosps))
row.names(northeast_adult_hosps_mat) <- NULL

cp_reg_northeast_adult <-change_point_test(northeast_adult_hosps_mat, w = rep(1 / nrow(northeast_adult_hosps_mat), nrow(northeast_adult_hosps_mat)), 
                              t_range = 2:(ncol(northeast_adult_hosps_mat) - 2), boot_num = 10000)

northeast_ped_hosps <- hhs_data_proc %>%
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "pediatric" & Region == "Northeast") %>%
  distinct(hospitalizations_standardized_region) 

northeast_ped_hosps_mat <- t(data.matrix(northeast_ped_hosps))
row.names(northeast_ped_hosps_mat) <- NULL

cp_reg_northeast_ped <- change_point_test(northeast_ped_hosps_mat, w = rep(1 / nrow(northeast_ped_hosps_mat), nrow(northeast_ped_hosps_mat)), 
                              t_range = 2:(ncol(northeast_ped_hosps_mat) - 2), boot_num = 10000)

South

south_adult_hosps <- hhs_data_proc %>% 
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "adult" & Region == "South") %>%
  distinct(hospitalizations_standardized_region) 

south_adult_hosps_mat <- t(data.matrix(south_adult_hosps))
row.names(south_adult_hosps_mat) <- NULL

cp_reg_south_adult <-change_point_test(south_adult_hosps_mat, w = rep(1 / nrow(south_adult_hosps_mat), nrow(south_adult_hosps_mat)), 
                              t_range = 2:(ncol(south_adult_hosps_mat) - 2), boot_num = 10000)

south_ped_hosps <- hhs_data_proc %>%
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "pediatric" & Region == "South") %>%
  distinct(hospitalizations_standardized_region) 

south_ped_hosps_mat <- t(data.matrix(south_ped_hosps))
row.names(south_ped_hosps_mat) <- NULL

cp_reg_south_ped <- change_point_test(south_ped_hosps_mat, w = rep(1 / nrow(south_ped_hosps_mat), nrow(south_ped_hosps_mat)), 
                              t_range = 2:(ncol(south_ped_hosps_mat) - 2), boot_num = 10000)

West

west_adult_hosps <- hhs_data_proc %>% 
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "adult" & Region == "West") %>%
  distinct(hospitalizations_standardized_region) 

west_adult_hosps_mat <- t(data.matrix(west_adult_hosps))
row.names(west_adult_hosps_mat) <- NULL

cp_reg_west_adult <- change_point_test(west_adult_hosps_mat, w = rep(1 / nrow(west_adult_hosps_mat), nrow(west_adult_hosps_mat)), 
                              t_range = 2:(ncol(west_adult_hosps_mat) - 2), boot_num = 10000)

west_ped_hosps <- hhs_data_proc %>%
  distinct(collection_week, Region, population, hospitalizations_standardized_region) %>% 
  arrange(collection_week) %>%
  filter(population == "pediatric" & Region == "West") %>%
  distinct(hospitalizations_standardized_region) 

west_ped_hosps_mat <- t(data.matrix(west_ped_hosps))
row.names(west_ped_hosps_mat) <- NULL

cp_reg_west_ped <- change_point_test(west_ped_hosps_mat, w = rep(1 / nrow(west_ped_hosps_mat), nrow(west_ped_hosps_mat)), 
                              t_range = 2:(ncol(west_ped_hosps_mat) - 2), boot_num = 10000)

Change Point Results

National Level

cp_nat_adult <- unlist(cp_nat_adult) %>% 
  data.frame()
colnames(cp_nat_adult)[1] <- "National Adult Hospitalizations"

cp_nat_ped <- unlist(cp_nat_ped) %>% 
  data.frame()
colnames(cp_nat_ped)[1] <- "National Pediatric Hospitalizations"

cp_nat_results <- cbind(cp_nat_adult, cp_nat_ped) %>% t() %>% data.frame()

cp_nat_results <- cp_nat_results %>%
  mutate(p.value = round(p.value,4),
         stat = round(stat,1)) %>%
  t()

rownames(cp_nat_results) <- c("Change Point Week", "P-value", "F-Stat")

datatable(cp_nat_results)

Regional Level

cp_reg_midwest_adult <- unlist(cp_reg_midwest_adult) %>% 
  data.frame()
colnames(cp_reg_midwest_adult)[1] <- "Midwest Adult"

cp_reg_midwest_ped <- unlist(cp_reg_midwest_ped) %>% 
  data.frame()
colnames(cp_reg_midwest_ped)[1] <- "Midwest Pediatric"

cp_reg_south_adult <- unlist(cp_reg_south_adult) %>% 
  data.frame()
colnames(cp_reg_south_adult)[1] <- "South Adult"

cp_reg_south_ped <- unlist(cp_reg_south_ped) %>% 
  data.frame()
colnames(cp_reg_south_ped)[1] <- "South Pediatric"

cp_reg_west_adult <- unlist(cp_reg_west_adult) %>% 
  data.frame()
colnames(cp_reg_west_adult)[1] <- "West Adult"

cp_reg_west_ped <- unlist(cp_reg_west_ped) %>% 
  data.frame()
colnames(cp_reg_west_ped)[1] <- "West Pediatric"

cp_reg_northeast_adult <- unlist(cp_reg_northeast_adult) %>% 
  data.frame()
colnames(cp_reg_northeast_adult)[1] <- "Northeast Adult"

cp_reg_northeast_ped <- unlist(cp_reg_northeast_ped) %>% 
  data.frame()
colnames(cp_reg_northeast_ped)[1] <- "Northeast Pediatric"

cp_reg_results <- cbind(cp_reg_midwest_adult, cp_reg_midwest_ped,
                        cp_reg_south_adult, cp_reg_south_ped,
                        cp_reg_west_adult, cp_reg_west_ped,
                        cp_reg_northeast_adult, cp_reg_northeast_ped) %>% 
  t() %>% 
  data.frame()

cp_reg_results <- cp_reg_results %>%
  mutate(p.value = round(p.value,4),
         stat = round(stat,1)) %>%
  t()

rownames(cp_reg_results) <- c("Change Point Week", "P-value", "F-Stat")

datatable(cp_reg_results)

Next, we will annotate plots with change points

#remotes::install_github("mattcowgill/ggannotate")
library(ggannotate) #ggannotate(plot) 

cp_peds <- cp_nat_results %>% 
  data.frame() %>%
  select(National.Pediatric.Hospitalizations) %>% 
  slice(1L) %>%
  as.numeric()

cp_peds_xy <- hhs_data_proc %>% 
  data.frame() %>% 
  distinct(date_id, collection_week, hospitalizations_standardized_us, population) %>% 
  filter(date_id == cp_peds, 
         population == "pediatric")


cp_adults <- cp_nat_results %>% 
  data.frame() %>%
  select(National.Adult.Hospitalizations) %>% 
  slice(1L) %>%
  as.numeric()

cp_adults_xy <- hhs_data_proc %>% 
  data.frame() %>% 
  distinct(date_id, collection_week, hospitalizations_standardized_us, population) %>% 
  filter(date_id == cp_adults, 
         population == "adult")

ggplot(hhs_data_proc %>% 
         mutate(hospitalizations_standardized_us = if_else(population == "pediatric", hospitalizations_standardized_us*10, hospitalizations_standardized_us)) %>% 
       distinct(collection_week, hospitalizations_standardized_us, population),
       aes(x = collection_week, y = hospitalizations_standardized_us, color = population)) + 
  geom_line() + 
  geom_point() +
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  ylab("Hospitalizations per 100K Adults") + 
  xlab("") + 
  scale_x_date(date_labels = "%b-%d", date_breaks = "1 month") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./10, name = "Hospitalizations per 100K Children")) + 
  geom_text(data = data.frame(x = cp_adults_xy$collection_week,
                              y = cp_adults_xy$hospitalizations_standardized_us,
                              label = "*"),
            mapping = aes(x = x, y = y, label = label),
            size = 12, hjust = 0.45, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_peds_xy$collection_week,
                              y = cp_peds_xy$hospitalizations_standardized_us*10, 
                              label = "*"),
            mapping = aes(x = x, y = y, label = label),
            size = 12, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) + 
  theme_bw() + 
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 15),
        legend.text = element_text(size=15),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 10),
        axis.title.y = element_text(face = "bold", size = 13, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 13, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 13, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 13, color='#D95F02'))

cp_peds <- cp_reg_results %>% 
  data.frame() %>%
  select(contains("Pediatric")) %>% 
  slice(1L) %>%
  t() 

region_col <- data.frame(Region = c("Midwest", "South", "West", "Northeast"))

cp_peds <- cbind(cp_peds, region_col)
rownames(cp_peds) <- NULL


cp_adults <- cp_reg_results %>% 
  data.frame() %>%
  select(contains("Adult")) %>% 
  slice(1L) %>%
  t() 

cp_adults <- cbind(cp_adults, region_col)
rownames(cp_adults) <- NULL

cp_peds_xy <- hhs_data_proc %>% 
  data.frame() %>% 
  left_join(., cp_peds, by = "Region") %>% 
  distinct(date_id, collection_week, hospitalizations_standardized_region, population, Region, `Change Point Week`)

cp_adults_xy <- hhs_data_proc %>% 
  data.frame() %>% 
  left_join(., cp_adults, by = "Region") %>% 
  distinct(date_id, collection_week, hospitalizations_standardized_region, population, Region, `Change Point Week`)



ggplot(hhs_data_proc %>%
         mutate(hospitalizations_standardized_region = if_else(population == "pediatric", 
                                                 hospitalizations_standardized_region*10, hospitalizations_standardized_region)) %>%
         distinct(collection_week, Region, population, hospitalizations_standardized_region),
       aes(x = collection_week, y = hospitalizations_standardized_region, color = population)) + 
  geom_point() +
  geom_line(alpha=0.9) +
  # midwest change point identifier
  geom_text(data = data.frame(x = cp_adults_xy %>% 
                                filter(Region == "Midwest" & population == "adult",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_adults_xy %>% 
                                filter(Region == "Midwest" & population == "adult",
                                       date_id == `Change Point Week`) %>% select(hospitalizations_standardized_region) %>%
                                rename('y' = hospitalizations_standardized_region), 
                              label = "*", Region = "Midwest"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#1B9E77", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_peds_xy %>% 
                                filter(Region == "Midwest" & population == "pediatric",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_peds_xy %>% 
                                filter(Region == "Midwest" & population == "pediatric",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                mutate(hospitalizations_standardized_region = hospitalizations_standardized_region*10) %>%
                                rename('y' = hospitalizations_standardized_region),
                              label = "*", Region = "Midwest"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#D95F02", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  # Northeast change point identifier
  geom_text(data = data.frame(x = cp_adults_xy %>% 
                                filter(Region == "Northeast" & population == "adult",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_adults_xy %>% 
                                filter(Region == "Northeast" & population == "adult",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                rename('y' = hospitalizations_standardized_region), 
                              label = "*", Region = "Northeast"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#1B9E77", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_peds_xy %>% 
                                filter(Region == "Northeast" & population == "pediatric",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_peds_xy %>% 
                                filter(Region == "Northeast" & population == "pediatric",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                mutate(hospitalizations_standardized_region = hospitalizations_standardized_region*10) %>%
                                rename('y' = hospitalizations_standardized_region),
                              label = "*", Region = "Northeast"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#D95F02", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  # South change point identifier
  geom_text(data = data.frame(x = cp_adults_xy %>% 
                                filter(Region == "South" & population == "adult",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_adults_xy %>% 
                                filter(Region == "South" & population == "adult",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                rename('y' = hospitalizations_standardized_region), 
                              label = "*", Region = "South"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#1B9E77", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_peds_xy %>% 
                                filter(Region == "South" & population == "pediatric",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_peds_xy %>% 
                                filter(Region == "South" & population == "pediatric",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                mutate(hospitalizations_standardized_region = hospitalizations_standardized_region*10) %>%
                                rename('y' = hospitalizations_standardized_region),
                              label = "*", Region = "South"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#D95F02", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) +
  # West change point identifier
  geom_text(data = data.frame(x = cp_adults_xy %>% 
                                filter(Region == "West" & population == "adult",
                                       date_id == `Change Point Week`) %>%
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_adults_xy %>% 
                                filter(Region == "West" & population == "adult",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                rename('y' = hospitalizations_standardized_region), 
                              label = "*", Region = "West"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#1B9E77", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_peds_xy %>% 
                                filter(Region == "West" & population == "pediatric",
                                       date_id == `Change Point Week`) %>% 
                                select(collection_week) %>%
                                rename('x' = collection_week),
                              y = cp_peds_xy %>% 
                                filter(Region == "West" & population == "pediatric",
                                       date_id == `Change Point Week`) %>%
                                select(hospitalizations_standardized_region) %>%
                                mutate(hospitalizations_standardized_region = hospitalizations_standardized_region*10) %>%
                                rename('y' = hospitalizations_standardized_region),
                              label = "*", Region = "West"),
            mapping = aes(x = x, y = y, label = label),
            colour = "#D95F02", size = 11, hjust = 0.5, vjust = 0.8, inherit.aes = FALSE) +
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  ylab("Hospitalizations per 100k Adults") + 
  xlab("") + 
  scale_x_date(date_labels = "%b", date_breaks = "1 month") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./10, name = "Hospitalizations per 100k Children")) + 
  theme_bw() +
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 15),
        legend.text = element_text(size=15),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 10),
        axis.title.y = element_text(face = "bold", size = 13, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 13, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 13, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 13, color='#D95F02')) + 
  facet_wrap(~Region, scales = "fixed") + 
  theme(strip.text.x = element_text(size = 18))