Hospitalization_Trends_Vaccine.Rmd

---
title: "COVID-19 Hospitalization Trends"
author: "Meghan Hutch, Molei Liu, Paul Avillach, Yuan Luo, Florence Bourgeois"
date: "4/26/2021"
output: html_document
---

**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)](https://healthdata.gov/Hospital/COVID-19-Reported-Patient-Impact-and-Hospital-Capa/g62h-syeh)

[Helpful FAQ Resource](https://github.com/CareSet/COVID_Hospital_PUF?src=hd)

Vaccination data was acquired through the [CDC](https://covid.cdc.gov/covid-data-tracker/#vaccination-demographics-trends)

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(ggplot2)
library(formattable)
library(data.table)
library(DT)
library(stringr)
#remotes::install_github("mattcowgill/ggannotate")
library(ggannotate) #ggannotate(plot) 
```

## **Data Import**

**Import HHS Hospitalization Counts and Census Data**

```{r message=FALSE, warning=FALSE}
hhs_data <- read_csv('data/COVID-19_Reported_Patient_Impact_and_Hospital_Capacity_by_State_Timeseries_20210607.csv')

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

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

load("processed_data/vaccine_age_table.rda")
```

## **Data Pre-Processing**

**Add full state names to HHS Data**

```{r}
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**

```{r}
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.

```{r}
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](https://www.census.gov/geographies/reference-maps/2010/geo/2010-census-regions-and-divisions-of-the-united-states.html)

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

**Add Census data to HHS data**

```{r}
hhs_data <- left_join(hhs_data, census_data, by = 'State')
```

**Convert to date format**

```{r}
hhs_data$date <- as.Date(hhs_data$date)
```

**Remove dates with no data collection**

We will remove dates when there was no```previous_day_admission_pediatric_covid_confirmed```cases 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

```{r}
#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 Week**

Similar to the original HHS data release, we will sum by week. 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.

```{r message=FALSE, warning=FALSE}
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)

# 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-06-04") # 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

```{r}
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**

```{r}
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 stratify into pediatric and adult datasets

```{r}
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**

```{r}
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

```{r}
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**

```{r}
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**

```{r}
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**

```{r}
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**

```{r}
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 Hospitalization Rates**

Rates are expressed as per 100,000 adults or children

```{r}
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**

```{r}
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)
```

```{r}
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

```{r}
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**

```{r}
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**

```{r}
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**

```{r}
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**

```{r}
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)))
  )
```

## **National Vaccination Rates **

```{r}
ggplot(vaccine_age_table %>% 
         distinct(week_id, total_pediatric_week_stan, total_adult_week_stan) %>% 
         mutate(total_pediatric_week_stan = total_pediatric_week_stan*5) %>% 
         pivot_longer(!week_id, names_to = "population", values_to = "counts"))  +
       aes(x = week_id, y = counts, group = population, color = population) + 
  geom_line() + 
  geom_point() + 
  ylab("Vaccinations per 100K Adults") + 
  xlab("") + 
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./5, name = "Vaccinations 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')) + 
  theme(strip.text.x = element_text(size = 18))  
```

## **Change Point Analysis**

```{r}
source("R/change_point_function.R")
```

Modify the dates of the HHS and vaccination dataframes

```{r}
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")

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

### **National Hospitalizations Change Point **

```{r}
set.seed(101)
```


```{r}
## 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**

```{r}
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**

```{r}
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**

```{r}
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**

```{r}
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)
```


### **Vaccination Change Point **

```{r}
## adult vaccines
adult_vaccine <- vaccine_age_table %>% 
  select(collection_week, total_adult_week_stan) %>% 
  arrange(collection_week) %>%
  distinct(total_adult_week_stan)

adult_vaccine_mat <- t(data.matrix(adult_vaccine))
row.names(adult_vaccine_mat) <- NULL

## peds vaccines
peds_vaccine <- vaccine_age_table %>% 
  select(collection_week, total_pediatric_week_stan) %>% 
  arrange(collection_week) %>%
  distinct(total_pediatric_week_stan)

peds_vaccine_mat <- t(data.matrix(peds_vaccine))
row.names(peds_vaccine_mat) <- NULL


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


cp_nat_vac_ped <- change_point_test(peds_vaccine_mat, w = rep(1 / nrow(peds_vaccine_mat), nrow(peds_vaccine_mat)), 
                              t_range = 2:(ncol(peds_vaccine_mat) - 2), boot_num = 10000)
```

### **Change Point Results**

**National Level**

```{r}
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**

```{r}
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)
```

**Vaccines**

Note: t.min is not at the same scale as hospitalization trends, since vaccination data began 12/18/2020

```{r}
cp_nat_vac_adult <- unlist(cp_nat_vac_adult) %>% 
  data.frame()
colnames(cp_nat_vac_adult)[1] <- "National Adult Vaccinations"

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

cp_nat_vac_results <- cbind(cp_nat_vac_adult, cp_nat_vac_ped) %>% t() %>% data.frame()

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

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

datatable(cp_nat_vac_results)
```

**Annotate plots with change points**

**Identify significance level based on number of analyses**

We wil use bonferroni correction to adjust our p-value. We have 12 analyses, 2 for each region, 2 for the national level analysis, and 2 for the vaccination analysis.

```{r}
sig_level = 0.05/12
```


**National Hospitalization Change points**

Note: Parameters to add vertical lines indicating vaccine change points (in progress!)

```{r}
cp_peds <- cp_nat_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level) %>% 
  filter(Analysis == "National Pediatric Hospitalizations") %>% 
  select(Change.Point.Week) %>% 
  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")

if(nrow(cp_peds_xy) == 0){
  
  cp_peds_xy[1,] <- NA
} 

cp_adults <- cp_nat_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level) %>% 
  filter(Analysis == "National Adult Hospitalizations") %>% 
  select(Change.Point.Week) %>% 
  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")

if(nrow(cp_adults_xy) == 0){
  
  cp_adults_xy[1,] <- NA
} 

# indicate vaccine change point dates
#vertical.lines <- c(as.Date("2021-03-05"), as.Date("2021-04-30"))
#vertical.lines.cols <- c("#1B9E77", "#D95F02")


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", date_breaks = "1 month") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./10, name = "Hospitalizations per 100K Children")) + 
  #geom_vline(xintercept = as.numeric(vertical.lines), linetype = "dashed", color = vertical.lines.cols) +
  theme_bw() + 
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 19),
        legend.text = element_text(size=19),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 17, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 17, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 17, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 17, color='#D95F02')) +
     geom_text(data = data.frame(x = cp_peds_xy$collection_week,
                              y = cp_peds_xy$hospitalizations_standardized_us*10, 
                              label = if_else(nrow(cp_peds_xy) != 0,"*", NULL)),
            mapping = aes(x = x, y = y, label = label),
            size = 16, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) +
     geom_text(data = data.frame(x = cp_adults_xy$collection_week,
                              y = cp_adults_xy$hospitalizations_standardized_us, 
                              label = if_else(nrow(cp_adults_xy) != 0,"*", NULL)),
            mapping = aes(x = x, y = y, label = label),
            size = 16, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE)
```

**Regional Change points**

```{r}
cp_peds <- cp_reg_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level, 
         grepl("Pediatric", Analysis)) %>% 
  select(Analysis, Change.Point.Week) %>% 
  mutate(Region = word(Analysis, 1))

cp_adults <- cp_reg_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level, 
         grepl("Adult", Analysis)) %>% 
  select(Analysis, Change.Point.Week) %>% 
  mutate(Region = word(Analysis, 1))

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)

# pediatric label indicators 
cp_midwest_peds <- cp_peds_xy %>% 
  filter(Region == "Midwest",
         population == "pediatric",
         date_id == Change.Point.Week)

if(nrow(cp_midwest_peds) == 0){
  
  cp_midwest_peds[1,] <- NA
} 

cp_northeast_peds <- cp_peds_xy %>% 
  filter(Region == "Northeast",
         population == "pediatric",
         date_id == Change.Point.Week)

if(nrow(cp_northeast_peds) == 0){
  
  cp_northeast_peds[1,] <- NA
} 

cp_west_peds <- cp_peds_xy %>% 
  filter(Region == "West",
         population == "pediatric",
         date_id == Change.Point.Week)

if(nrow(cp_west_peds) == 0){
  
  cp_west_peds[1,] <- NA
} 

cp_south_peds <- cp_peds_xy %>% 
  filter(Region == "South",
         population == "pediatric",
         date_id == Change.Point.Week)

if(nrow(cp_south_peds) == 0){
  
  cp_south_peds[1,] <- NA
}


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)

# adult label indicators
cp_midwest_adult <- cp_adults_xy %>% 
  filter(Region == "Midwest",
         population == "adult",
         date_id == Change.Point.Week)

if(nrow(cp_midwest_adult) == 0){
  
  cp_midwest_adult[1,] <- NA
} 

cp_northeast_adult <- cp_adults_xy %>% 
  filter(Region == "Northeast",
         population == "adult",
         date_id == Change.Point.Week)

if(nrow(cp_northeast_adult) == 0){
  
  cp_northeast_adult[1,] <- NA
} 

cp_west_adult <- cp_adults_xy %>% 
  filter(Region == "West",
         population == "adult",
         date_id == Change.Point.Week)

if(nrow(cp_west_adult) == 0){
  
  cp_west_adult[1,] <- NA
} 

cp_south_adult <- cp_adults_xy %>% 
  filter(Region == "South",
         population == "adult",
         date_id == Change.Point.Week)

if(nrow(cp_south_adult) == 0){
  
  cp_south_adult[1,] <- NA
}

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 labels
  geom_text(data = data.frame(x = cp_midwest_adult$collection_week,
                              y = cp_midwest_adult$hospitalizations_standardized_region, 
                              label = if_else(nrow(cp_midwest_adult) != 0,"*", NULL),
                              Region = "Midwest"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_midwest_peds$collection_week,
                              y = cp_midwest_peds$hospitalizations_standardized_region*10, 
                              label = if_else(nrow(cp_midwest_peds) != 0,"*", NULL),
                              Region = "Midwest"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) + 
  # northeast labels
  geom_text(data = data.frame(x = cp_northeast_adult$collection_week,
                              y = cp_northeast_adult$hospitalizations_standardized_region, 
                              label = if_else(nrow(cp_northeast_adult) != 0,"*", NULL),
                              Region = "Northeast"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_northeast_peds$collection_week,
                              y = cp_northeast_peds$hospitalizations_standardized_region*10, 
                              label = if_else(nrow(cp_northeast_peds) != 0,"*", NULL),
                              Region = "Northeast"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) + 
   # west labels
  geom_text(data = data.frame(x = cp_west_adult$collection_week,
                              y = cp_west_adult$hospitalizations_standardized_region, 
                              label = if_else(nrow(cp_west_adult) != 0,"*", NULL),
                              Region = "West"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_west_peds$collection_week,
                              y = cp_west_peds$hospitalizations_standardized_region*10, 
                              label = if_else(nrow(cp_west_peds) != 0,"*", NULL),
                              Region = "West"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) + 
  # south labels
  geom_text(data = data.frame(x = cp_south_adult$collection_week,
                              y = cp_south_adult$hospitalizations_standardized_region, 
                              label = if_else(nrow(cp_south_adult) != 0,"*", NULL),
                              Region = "South"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE) + 
  geom_text(data = data.frame(x = cp_south_peds$collection_week,
                              y = cp_south_peds$hospitalizations_standardized_region*10, 
                              label = if_else(nrow(cp_south_peds) != 0,"*", NULL),
                              Region = "South"),
            mapping = aes(x = x, y = y, label = label),
            size = 13, hjust = 0.5, vjust = 0.75, colour = "#D95F02", 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))  
```

**Vaccination Change points**

```{r}
cp_vac_peds <- cp_nat_vac_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level, 
         grepl("Pediatric", Analysis)) %>% 
  distinct(Change.Point.Week) %>% 
  as.numeric()

cp_vac_adults <- cp_nat_vac_results %>% 
  t() %>% 
  data.frame() %>% 
  tibble::rownames_to_column(., "Analysis") %>% 
  filter(P.value < sig_level, 
         grepl("Adult", Analysis)) %>% 
  distinct(Change.Point.Week) %>% 
  as.numeric()

cp_vac_peds_xy <- vaccine_age_table %>% 
  data.frame() %>% 
  distinct(week_id, collection_week, total_pediatric_week_stan) %>% 
  filter(week_id == cp_vac_peds)

if(nrow(cp_vac_peds_xy) == 0){
  
  cp_vac_peds_xy[1,] <- NA
} 

cp_vac_adults_xy <- vaccine_age_table %>% 
  data.frame() %>% 
  distinct(week_id, collection_week, total_adult_week_stan) %>% 
  filter(week_id == cp_vac_adults)

if(nrow(cp_vac_adults_xy) == 0){
  
  cp_vac_adults_xy[1,] <- NA
} 


ggplot(vaccine_age_table %>% 
         distinct(collection_week, week_id, total_pediatric_week_stan, total_adult_week_stan) %>% 
         mutate(total_pediatric_week_stan = total_pediatric_week_stan*5) %>% 
         pivot_longer(!c(collection_week, week_id), names_to = "population", values_to = "counts"))  +
       aes(x = collection_week, y = counts, group = population, color = population) + 
  geom_line() + 
  geom_point() + 
  ylab("Vaccinations per 100K Adults") + 
  xlab("") + 
  scale_x_date(date_labels = "%b-%d", date_breaks = "1 month") + 
  scale_color_brewer(palette="Dark2", labels = c("Adult", "Child"), name = "") + 
  scale_y_continuous(sec.axis = sec_axis(~ ./5, name = "Vaccinations per 100k Children")) + 
  theme_bw() + 
  theme(plot.title = element_text(""),
        legend.title = element_text(face = "bold", size = 19),
        legend.text = element_text(size=19),
        legend.position = "top",
        legend.direction = "horizontal",
        axis.text.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 17, color = "#1B9E77", margin=margin(0,10,0,0)),
        axis.text.y.left = element_text(face = "bold",  size = 17, color = "#1B9E77"),
        axis.title.y.right = element_text(face = "bold", size = 17, color='#D95F02', margin=margin(0,0,0,10)),
        axis.text.y.right = element_text(face = "bold", size = 17, color='#D95F02')) + 
  geom_text(data = data.frame(x = cp_vac_peds_xy$collection_week,
                              y = cp_vac_peds_xy$total_pediatric_week_stan*5, 
                              label = if_else(nrow(cp_vac_peds_xy) != 0,"*", NULL)),
            mapping = aes(x = x, y = y, label = label),
            size = 16, hjust = 0.5, vjust = 0.75, colour = "#D95F02", inherit.aes = FALSE) +
     geom_text(data = data.frame(x = cp_vac_adults_xy$collection_week,
                              y = cp_vac_adults_xy$total_adult_week_stan, 
                              label = if_else(nrow(cp_vac_adults_xy) != 0,"*", NULL)),
            mapping = aes(x = x, y = y, label = label),
            size = 16, hjust = 0.5, vjust = 0.75, colour = "#1B9E77", inherit.aes = FALSE)
```