equity.Rmd

---
title: "Explore Equity"
description: |
  Are vaccines being responsibly and fairly administered?  <ul class="list-inline title-nav"> 
  <li class="list-inline-item-nav"><div class = "go-to-nav" >Go to: </div></li>
  <li class="list-inline-item-nav">
    <a href="providers.html"> <button type="button" class="btn btn-outline-light mb-2"> Providers </button> 
    </a>
  </li>
  <li class="list-inline-item-nav">
    <a href="pace.html"> <button type="button" class="btn btn-outline-light mb-2"> Pace </button> 
    </a>
  </li>
  </ul>
output: 
  distill::distill_article:
    toc: true
    css: equity_header.css
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(gt)
library(tidyverse)
library(highcharter)
library(readxl)
library(sf)
library(ggtext)
library(janitor)
library(ggforce)
library(mapboxapi)
library(leaflet)
library(tidycensus)
options(tigris_use_cache = TRUE)

today <- Sys.Date()
src <- "https://dshs.texas.gov/immunize/covid19/COVID-19-Vaccine-Data-by-County.xls" # URL For Daily COVID Data
lcl <- basename(src)
filepath <- paste0("raw_data/", today,"_",lcl,"x")
download.file(url = src, destfile = filepath) # Downloads The Data
today_lbl <- format(Sys.Date(), format = "%b %d, %Y")
# today_lbl <- format(Sys.Date(), format = "%b %d, %Y")
time_lbl <- format(Sys.time(), format = "%I:%M %p")

vaccine_hc <- hc_theme_merge(
  hc_theme_smpl(),
  hc_theme(chart = list(backgroundColor = "#fff", 
                        style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', 
                                     fontSize = "28px",  fontWeight = "bold",
                                     color="#4d4d4d", textTransform="initial")),
           title = list(style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', 
                                     fontWeight="bold",
                                     color="#2d2d2d"),
                        align = "left"), 
           subtitle = list(style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif',
                                        color="#7d7d7d",
                                        textTransform="initial",
                                        fontWeight="400",
                                        fontSize = "14px"),
                           align = "left"), 
           legend = list(align = "right", 
                         style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', color="#6d6d6d"),
                         itemStyle = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', color = '#6d6d6d'),
                         itemHoverStyle = list(color = 'gray'),   
                         verticalAlign = "top"),
           credits = list(style = list(color = "#2d2d2d")),
           xAxis = list(labels =list(style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', 
                                                  color="#2d2d2d")), 
                        title = list(style = list(color = "#4d4d4d", fontSize = "12px", 
                                                  color="#4d4d4d",fontWeight="500")),
                        gridLineWidth = 0,
                        gridLineColor = "#9d9d9d", 
                        lineColor = 'rgba(0,0,0,0.7)', 
                        minorGridLineColor = 'rgba(0,0,0,0.7)', 
                        tickColor = "#9d9d9d", 
                        tickWidth = 1), 
           yAxis = list(labels =list(style = list(fontFamily = '"Libre Franklin", -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", Arial, sans-serif', color="#2d2d2d")), 
                        title = list(style = list(color = "#fff", fontSize = "12px", 
                                                  color="#fff",fontWeight="500")), 
                        gridLineWidth = .5,
                        gridLineColor = 'rgba(0,0,0,0.15)', 
                        lineColor = 'rgba(0,0,0,0.15)', 
                        minorGridLineColor = 'rgba(0,0,0,0.15)', 
                        tickColor = "#9d9d9d", 
                        tickWidth = 2)))
```


# Vaccines + Equity

## Statewide Demographic Summaries

Towards the end of January, the Texas Department of State Health Services (DSHS) began publishing breakdowns of vaccinations by across different demographics, specifically age, gender, and race/ethnicity. While there's a lot of unknowns in the data--particularly across entries for "Race/Ethnicity"+"Gender"--and also because we don't have the population estimates for phases 1A and 1B broken down by race/age/gender, the data presented by DSHS struggles to paint an incomplete picture of the profile of vaccinated individuals. That said, the charts below attempt to breakdown that data with as much clarity as possible across all known and unknown classifications. 

```{r covid_case_fatality_dem_chart, echo=FALSE, fig.height=4.75, fig.width=12.25, layout = "l-page"}

## PULL ACS CENSUS DEMOGRAPHIC DATA ##

demvars <-  c(White = "B03002_003",
              Black = "B03002_004",
              Asian = "B03002_006",
              Hispanic = "B03002_012",
              Other = "B03002_005", #AIAN
              Other = "B03002_007", #NHPI
              Other = "B03002_008", #Some Other Race
              Other = "B03002_009") #Two Races

dem_tx_data <- get_acs(geography = "state",
                       variables = demvars,
                       state="TX",
                       survey = "acs5",
                       show_call = TRUE,
                       summary_var = "B01003_001") %>%
  group_by(variable) %>% 
  summarise(across(where(is.numeric), sum)) %>% 
  mutate(measure = "State Population",
         data="ACS 2019",
         pct = round(estimate/summary_est, digits=4)) %>% 
  select(data, measure, dem_group = variable, value = estimate, pct)

## PULL DSHS COVID CASE FATALITY DEMOGRAPHIC DATA ##

today <- Sys.Date()
src_cf_data <- "https://dshs.texas.gov/coronavirus/TexasCOVID19Demographics.xlsx.asp" # URL For Daily CF Data
src_cf_file <- "https://dshs.texas.gov/coronavirus/TexasCOVID19Demographics.xlsx" # URL For Daily CF Data
lcl_cf_data <- basename(src_cf_file)
# filepath_cf_data <- paste0("raw_data/",today,"_",lcl_cf_data)
# download.file(url = src_cf_data, destfile = filepath_cf_data) # Downloads The Data

filepath_cf_data <- "raw_data/2021-04-05_TexasCOVID19Demographics.xlsx"

sheets <- excel_sheets(filepath_cf_data) %>% set_names()

dem_cf_data <- map_df(sheets,
                    ~read_excel(filepath_cf_data, sheet=.x) %>% 
                      rename(dem_group=1,value=2,pct=3),
                    .id="Sheet") %>% 
  separate(Sheet, into=c("measure","dem_type"), sep = " by ") %>% 
  filter(str_detect(dem_type, regex("Race", ignore_case = TRUE)),
         !str_detect(dem_group, regex("Total", ignore_case = TRUE))) %>% 
  select(-dem_type) %>% 
  mutate(data="Case + Fatality")

## PULL DSHS COVID VACCINE DEMOGRAPHIC DATA ##

dem_vax_data_summ <- read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
  clean_names() %>%
  filter(!str_detect(gender, "Texas"),
         !str_detect(race_ethnicity, "Total")) %>% 
  group_by(dem_group=race_ethnicity) %>% 
  summarise(across(where(is.numeric), sum)) %>% 
  ungroup() %>% 
  pivot_longer(cols=2:4,
               names_to = "measure",
               values_to = "value") %>% 
  group_by(measure) %>% 
  summarise(sum = sum(value))

dem_vax_data <- read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
  clean_names() %>%
  filter(!str_detect(gender, "Texas"),
         !str_detect(race_ethnicity, "Total")) %>% 
  group_by(dem_group=race_ethnicity) %>% 
  summarise(across(where(is.numeric), sum)) %>% 
  ungroup() %>% 
  pivot_longer(cols=2:4,
               names_to = "measure",
               values_to = "value") %>%
  left_join(dem_vax_data_summ, by = "measure") %>% 
  mutate(pct=round(value/sum, digits=4)) %>% 
  select(-sum) %>% 
  mutate(measure = str_replace(measure,"doses_administered","Doses Administered"),
         measure = str_replace(measure,"people_vaccinated","Partially Vaccinated"),
         measure = str_replace(measure,"people_fully_vaccinated","Fully Vaccinated"),
         data="Vaccines")

## CREATE UNIFIED DEMOGRAPHIC DATASET ##

dem_data_full <- bind_rows(dem_cf_data,
                           dem_tx_data,
                           dem_vax_data) %>% 
  select(data, measure,dem_group,value,pct) %>% 
  mutate(measure_f = factor(measure, levels=c('Doses Administered','Partially Vaccinated',
                                              'Fully Vaccinated','Cases','Fatalities', 
                                              'State Population'))) %>% 
  # left_join(dem_tx_data, by = "dem_group") %>% 
  # filter(measure!='Doses Administered',
  #        dem_group!="Other") %>% 
  mutate(measure = paste0("% of ",measure),
         measure_f = paste0("% of ",measure_f),
         pct_lbl = scales::percent(pct, accuracy = 0.1)) %>% 
  filter(!str_detect(measure,'Partially Vaccinated|Fully Vaccinated'),
         dem_group!="Other") %>% 
  mutate(measure_f = factor(measure, levels=c('% of Doses Administered','% of Partially Vaccinated',
                                              '% of Fully Vaccinated','% of Cases','% of Fatalities', 
                                              '% of State Population'))) %>% 
  arrange((measure_f))

equity_chart <-
  ggplot() +
    geom_col(data=dem_data_full %>% 
               # filter(dem_group!="Unknown") %>% 
               select(-measure),
             aes(y=reorder(dem_group, pct), x=pct, fill=dem_group), 
             stat="identity", na.rm=TRUE, width=.85) +
    facet_wrap(~measure_f, ncol=4) +
    geom_text(data=dem_data_full %>% 
                # filter(dem_group!="Unknown") %>% 
                select(-measure),
              aes(y=reorder(dem_group, pct), x=pct*1.05, label =pct_lbl), hjust=0, color="#3d3d3d", family = "Graphik-Bold") +
  hrbrthemes::scale_x_percent(limits = c(0,1)) +
  ggthemes::scale_fill_tableau(palette = "Tableau 10") +
  guides(fill=guide_legend(title = md("**Legend** | Dashed Box Indicates 2019 ACS Population"), 
                           title.position = "top",
                           keyheight = .8)) +
  theme_minimal(base_family = "Graphik") +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =20),
        plot.subtitle = element_text(family = "Graphik-Regular", color = "#333f48", margin = margin(b=15)),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.margin = margin(t=15, r=15, b=10, l = 15),
        plot.caption = element_markdown(lineheight = 1),
        # axis.text.x = element_text(angle =),
        legend.title = element_markdown(), 
        strip.text = element_text(family = "Graphik-Bold", color = "#9cadb7", size =14),
        legend.position = "none") +
  labs(title = "Comparing Cases, Fatalities, & Vaccines to Population Trends",
       subtitle = paste0("This chart shows break downs of known COVID Cases, COVID Fatalities, Vaccine Data for Texas, and the Statewide Population, by race and ethnicity.\nExcluded are the overall race and ethnic groups classified as 'Unknown' or 'Other' by DSHS. Administered doses is a summary figure of persons who have\nreceived a single dose and persons who are fully vaccinated.\nData pulled from sources on ", today_lbl, " at ", time_lbl),
       x = NULL,
       y = NULL,
       color = NULL,
       # fill = "Legend",
       caption = md("**Population Data**: 2019 ACS 5-Yr Series, Pulled with Tidycensus.<br>**Case Fatality Data**: 'Accessible Dashboard Data', Texas DSHS<br>**Vaccine Data**: 'Accessible Vaccine Dashboard Data', Texas DSHS"))

equity_chart
    
ggsave(paste0("figures/", today,"_","tx_covid_equity_analysis.png"),equity_chart,  dpi = 300, width = 12, height = 4.75)
  
```

The next chart is broken up into six different panels--one for each classification of a race or ethnicity. Within each of those panels, there's two sets of five-columned bar charts which represent the reported gender. The dark blue colored bars (facing upwards) represent the female population within each race and ethnicity classification while the gray colored bars (facing downwards) represent the male population within each race and ethnicity classification. Each column within a set of blue or gray colored bars represents a different age group. 

```{r dem_chart_hc, eval=FALSE, include=FALSE, layout="l-page"}

re_groups <- read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
  clean_names() %>% 
  select(race_ethnicity) %>% 
  distinct() %>% 
  simplify()

df <- read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
      clean_names() %>%
      filter(!str_detect(gender, "Unknown|Texas"),
         !str_detect(race_ethnicity, "Total")) %>% 
      mutate(value = ifelse(gender == "Male", -doses_administered, doses_administered),
             age_group = str_remove(age_group, " years")) %>% 
  arrange(race_ethnicity,age_group)

categories_grouped <- df %>%
  select(race_ethnicity, age_group) %>%
  distinct() %>% 
  group_by(name = race_ethnicity) %>% 
  summarise(categories = list(age_group)) %>% 
  list_parse()

df %>% 
  hchart(type="column", 
         hcaes(x = age_group, y = value, group=gender)) %>% 
  hc_title(text="Texas Distribution of All Administered Doses\nby Race, Age, and Reported Gender") %>% 
  hc_subtitle(text = paste0("This chart shows break downs by the number of doses administered in Texas. Excluded are the overall Texas number and gender\nthat the state reported as unknown. Administered doses is a summary figure of persons who have received a single dose\nand persons who are fully vaccinated. Source: Texas Department of State Health Services. As of ", today_lbl, " at ", time_lbl)) %>% 
  hc_yAxis(opposite = TRUE,
           labels = list(formatter = JS(paste0("function () {return Math.abs(this.value)/1000 + 'K';}")))) %>% 
  hc_xAxis(categories = categories_grouped,
           labels = list(groupedOptions = JS("[{
                         style: {color: '#bf5700'},
                         rotation: 0
                         }]"),
                         rotation = -90),
           title = NA,
           opposite=FALSE) %>%
  hc_add_dependency("plugins/grouped-categories.js") %>% 
  # hc_legend(enabled=FALSE) %>%
  # hc_tooltip(table = TRUE,
  #            sort = TRUE,
  #            shared=TRUE) %>%
  hc_plotOptions(series = list(stacking = "normal")) %>%
  hc_colors(colors = list('#005f86', '#9cadb7')) %>% 
  # hc_colors(colors = list('rgba(0, 95, 134,0.3)','rgba(0, 95, 134,1)')) %>%
  hc_caption("Chart by: Matt Worthington") %>%
  hc_credits(
    enabled = TRUE,
    text = "Data From: Texas DSHS 'Accessible Vaccine Dashboard Data'",
    href = "https://dshs.texas.gov/immunize/covid19/COVID-19-Vaccine-Data-by-County.xls") %>%
  hc_add_theme(vaccine_hc)
```


```{r dem_chart_ggplot, echo=FALSE, fig.height=6.5, fig.width=10, layout="l-page"}

# vaccine_summ <- 
  read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
  clean_names() %>%
  filter(!str_detect(gender, "Unknown|Texas"),
         !str_detect(race_ethnicity, "Total")) %>% 
  mutate(value = ifelse(gender == "Male", -doses_administered, doses_administered)) %>% 
  ggplot(aes(x = age_group, y = value, fill = gender)) +
  geom_bar(stat = "identity", width = 1) +
  # coord_flip() +
  facet_wrap(~race_ethnicity, nrow=1) +
  scale_y_continuous(labels = function(y) paste0(abs(y / 1000), "k")) +
  scale_x_discrete(labels = function(x) gsub("Age | years", "", x)) +
  scale_fill_manual(values = c("#005f86","#9cadb7","#bf5700")) +
  guides(fill=guide_legend(title = "Legend", keyheight = .8)) +
  theme_minimal(base_family = "Graphik") +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =20),
        plot.subtitle = element_text(family = "Graphik-Regular", color = "#333f48", margin = margin(b=15)),
        panel.grid.minor.y = element_blank(),
        plot.margin = margin(t=15, r=15, b=10, l = 15),
        axis.text.x = element_text(angle = 90),
        legend.position = "top") +
  labs(title = "Texas Distribution of All Administered Doses\nby Race, Age, and Reported Gender",
       subtitle = paste0("This chart shows break downs by the number of doses administered in Texas. Excluded are the overall Texas number and gender\nthat the state reported as unknown. Administered doses is a summary figure of persons who have received a single dose\nand persons who are fully vaccinated. Source: Texas Department of State Health Services. As of ", today_lbl, " at ", time_lbl),
       x = NULL,
       y = NULL,
       color = NULL,
       caption = "Data File: 'Accessible Vaccine Dashboard Data'\nLink: https://dshs.texas.gov/immunize/covid19/COVID-19-Vaccine-Data-by-County.xls")
  # ggsave(paste0("figures/",today,"_","tx_shots_by_race_age_known_gender.png"), dpi = 300, width = 11, height = 8)

  
  read_excel(filepath, sheet = "By Age, Gender, Race", skip=0) %>% 
  clean_names() |> 
  filter(str_detect(gender, "Unknown"),
         !str_detect(race_ethnicity, "Total")) |> 
  ggplot(aes(x = age_group, y = doses_administered, fill = age_group)) +
  geom_bar(stat = "identity", width = 1) +
  facet_wrap(~race_ethnicity, nrow=1) +
  scale_y_continuous(labels = function(y) paste0(abs(y / 1000), "k")) +
  scale_x_discrete(labels = function(x) gsub("Age | years", "", x)) +
  scale_fill_manual(values = c("#f8971f", "#005f86","#9cadb7","#bf5700", "#333f48", "#579d42")) +
  guides(fill=guide_legend(title = "Legend", keyheight = .8)) +
  theme_minimal(base_family = "Graphik") +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =20),
        plot.subtitle = element_text(family = "Graphik-Regular", color = "#333f48", margin = margin(b=15)),
        panel.grid.minor.y = element_blank(),
        plot.margin = margin(t=15, r=15, b=10, l = 15),
        axis.text.x = element_text(angle = 90),
        legend.position = "top") +
  labs(title = "Texas Distribution of All Administered Doses\nby Race, Age, and Unknown Gender",
       subtitle = paste0("This chart shows break downs by the number of doses administered in Texas. Excluded are the overall Texas number and genders\nthat the state reported as known. Administered doses is a summary figure of persons who have received a single dose\nand persons who are fully vaccinated. Source: Texas Department of State Health Services. As of ", today_lbl, " at ", time_lbl),
       x = NULL,
       y = NULL,
       color = NULL,
       caption = "Data File: 'Accessible Vaccine Dashboard Data'\nLink: https://dshs.texas.gov/immunize/covid19/COVID-19-Vaccine-Data-by-County.xls")
  ggsave(paste0("figures/",today,"_","tx_shots_by_race_age_unknown_gender.png"), dpi = 300, width = 11, height = 8)
  

```

## Spatial Distributions, by County

One of the datasets that the Texas Department of State Health Services publishes is zip-code level administration data. While seeing the data at a statewide level with zip-codes is helpful, it also compares geographic areas that may have access to different levels of supply. Considering that, we've created a resource that allows users to select a county and visualize distribution trends at the county level so that you can spatial trends for zip codes in a county when compared only to other zip codes in a county. Once a user selects a county, the map shows you the rate of doses administered for every 1,000 persons in each zip code.

```{r vaccine-zip-cnty-map, echo=FALSE, message=FALSE, warning=FALSE, layout = "l-screen"}
knitr::include_app("https://mrworthington.shinyapps.io/cnty_zip_explorer/", height = "820px")
```


# Vaccines + Access

```{r vaccine-available-by-type, eval=FALSE, fig.width=10, include=FALSE, layout="l-page"}

today <- Sys.Date()
src <- "https://genesis.soc.texas.gov/files/accessibility/vaccineprovideraccessibilitydata.csv"
lcl <- basename(src)
filepath <- paste0(today,"_",lcl)
download.file(url = src, destfile = filepath)

# provider_locations <- read_rds("clean_data/vaccine_supply_by_site.rds") %>%
  # select(name, street, lat, lon)

supply_by_provider_type <- read_csv(filepath, skip=0) %>%
  janitor::clean_names() %>% 
  group_by(type) %>% 
  summarise(tot_shipped = sum(total_shipped),
            tot_avail = sum(total_available)) 

# supply_by_provider_type %>% 
#   drop_na() %>% 
#   ggplot(aes(x=reorder(type, desc(tot_shipped)), y = tot_shipped, fill=type)) +
#   geom_col() +
#   theme_minimal(base_family = "Graphik") +
#   geom_text(angle=90, color="#2d2d2d", size = 3.1, family = "Graphik-Semibold",
#             aes(y= 0, x=type, label = type), hjust = 0) +
#   scale_y_continuous(labels = scales::unit_format(unit = "K", accuracy = 1, scale = 1e-3)) +
#   ggthemes::scale_fill_tableau(palette = "Tableau 20") +
#   theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =20),
#         plot.subtitle = element_text(family = "Graphik-Regular", color = "#333f48", margin = margin(b=15)),
#         panel.grid.major.y = element_blank(),
#         plot.margin = margin(t=15, r=15, b=10, l = 15),
#         axis.text.x = element_blank(),
#         legend.position = "none") +
#   labs(title = "Reported Total Supply by Provider Type",
#        subtitle = paste0("Source: Texas Department of State Health Services. As of ", today_lbl, " at ", time_lbl),
#        x = NULL,
#        y = "Doses (Shipped)",
#        color = NULL,
#        caption = "Data File: 'Texas Vaccine Provider Data'\nLink: https://genesis.soc.texas.gov/files/accessibility/vaccineprovideraccessibilitydata.csv\nChart by: @mrworthington") #+
#   # ggsave(paste0("figures/",today,"_","supply_by_provider_type.png"), dpi = 300, width = 11, height = 4.5)

supply_by_provider_type %>% 
  drop_na() %>% 
  ggplot(aes(x=reorder(type, desc(tot_avail)), y = tot_avail, fill=type)) +
  geom_col() +
  theme_minimal(base_family = "Graphik") +
  geom_text(angle=90, color="#2d2d2d", size = 3.1, family = "Graphik-Semibold",
            aes(y= 0, x=type, label = type), hjust = 0) +
  scale_y_continuous(labels = scales::unit_format(unit = "K", accuracy = 1, scale = 1e-3)) +
  ggthemes::scale_fill_tableau(palette = "Tableau 20") +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =20),
        plot.subtitle = element_text(family = "Graphik-Regular", color = "#333f48", margin = margin(b=15)),
        panel.grid.major.y = element_blank(),
        plot.margin = margin(t=15, r=15, b=10, l = 15),
        axis.text.x = element_blank(),
        legend.position = "none") +
  labs(title = "Reported Availability by Provider Type",
       subtitle = paste0("Source: Texas Department of State Health Services. As of ", today_lbl, " at ", time_lbl),
       x = NULL,
       y = "Doses (Available)",
       color = NULL,
       caption = "Data File: 'Texas Vaccine Provider Data'\nLink: https://genesis.soc.texas.gov/files/accessibility/vaccineprovideraccessibilitydata.csv") #+
  # ggsave(paste0("figures/",today,"_","availability_by_provider_type.png"), dpi = 300, width = 11, height = 4.5)
```

<!-- ##  Where Accessible Populations Live in Texas -->

```{r walking_accessibility, eval=FALSE, include=FALSE}


vaccine_sites_sf <- read_csv("https://utexas-lbjp-data.github.io/txcovidvaccinedata/clean_data/vaccine_provider_raw.csv",col_types = cols(time_pulled = col_character())) %>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>% 
    select(-Address) %>% 
    mutate(Address = paste0(Street, ", ", City, ", TX ", Zip)) %>% 
    select(-Street, -City, -Zip) %>% 
    mutate(Time = format(as.POSIXct(`Time Updated`), format = "%I:%M %p")) %>% 
    select(1:2, `Total Available`, `Last Updated`,Time, Phone, Address, everything())

vaccine_sites_sf_travis <- vaccine_sites_sf %>% 
  filter(County=="Travis")

# walking_isos <- mb_isochrone(
#   vaccine_sites_sf_travis,
#   profile = "walking",
#   time = 30,
#   id = "Provider Name",
#   rate_limit = 2500,
# ) %>% 
#   st_transform(4326) %>% 
#   st_as_sf() 
# 
# walking_combined <- st_union(walking_isos) %>% 
#   st_cast(to="POLYGON") %>% 
#   st_transform(4326) %>% 
#   st_as_sf() 
# 
# write_rds(walking_isos, "clean_data/walking_isos_30.rds")
# write_rds(walking_combined, "clean_data/walking_combined_30.rds")

walking_isos <- read_rds("clean_data/walking_isos_30.rds")
walking_combined <- read_rds("clean_data/walking_combined_30.rds")

# vars_2019 <- load_variables(2019, "acs5") #%>% 
#   filter(str_detect(concept,"POPULATION 12 YEARS AND OVER"),
#          label=="Estimate!!Total:")

# demvars <-  c(White = "B03002_003",
#               Black = "B03002_004",
#               Asian = "B03002_006",
#               Hispanic = "B03002_012")
# 
# travis_dems <- get_acs(geography = "block group", 
#                        variables = demvars,
#                        state="TX",
#                        county = "Travis County",
#                        survey = "acs5",
#                        show_call = TRUE,
#                        geometry = TRUE,
#                        summary_var = "B01003_001")

# travis_pop_blox <- get_acs(geography = "block group", 
#                       variables = "B01003_001",
#                       state="TX",
#                       county = "Travis County",
#                       survey = "acs5",
#                       geometry = TRUE,
#                       show_call = TRUE) %>% 
#   filter(estimate!=0) %>% 
#   st_transform(4326) %>% 
#   st_as_sf()

inaccessible_walking <- st_difference(travis_pop_blox,st_union(walking_combined)) %>%
  st_cast(to="POLYGON") %>%
  st_transform(4326) %>%
  st_as_sf()  %>%
  mutate(est_lab = scales::comma(estimate, accuracy=1, suffix = " People"))
# 
# inaccessible_walk_summ <- inaccessible_walking %>% 
#   as_tibble() %>% 
#   distinct(GEOID, .keep_all = TRUE) %>% 
#   mutate(state="TX") %>% 
#   group_by(state) %>% 
#   summarise(sum_inaccessible = sum(estimate))

mapbox_map <- leaflet() %>%
  addMapboxTiles(style_id = "light-v9",
                 username = "mapbox") 

# bins <- c(0, 10, 20, 50, 100, 200, 500, 1000, Inf)
pal <- colorBin("Reds", domain = inaccessible_walking$estimate)

mapbox_map %>%
  addPolygons(data = inaccessible_walking,
              stroke = FALSE,
              color = "#03F",
              fill = TRUE,
              # fillColor = ~pal(estimate),
              fillOpacity = 0.7,
              dashArray = NULL,
              smoothFactor = 1,
              popup = ~est_lab) %>% 
  addLegend("bottomleft",
            data = inaccessible_walking,
            values = ~estimate,
            # labFormat = labelFormat(
            #   suffix = "%"),
            pal = pal,
            title = "Population Est.")
```

## Where Inaccessible Populations Live in Texas

### From Vaccine Hubs

```{r inaccessible_hub_df, include=FALSE, layout="l-page"}

vaccine_sites_sf <- read_csv("https://utexas-lbjp-data.github.io/txcovidvaccinedata/clean_data/vaccine_provider_raw.csv") %>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>% 
    select(-Address) %>% 
    mutate(Address = paste0(Street, ", ", City, ", TX ", Zip)) %>% 
    select(-Street, -City, -Zip) %>% 
    mutate(Time = format(as.POSIXct(`Time Updated`), format = "%I:%M %p")) %>% 
    select(1:2, `Total Available`, `Last Updated`,Time, Phone, Address, everything())

vaccine_sites_sf_hub <- vaccine_sites_sf %>% 
  filter(`Organization Type`=="Vaccine Hub")

# driving_isos_30 <- mb_isochrone(
#   vaccine_sites_sf_hub,
#   profile = "driving",
#   time = 30,
#   id = "Provider Name"
# ) %>%
#   st_transform(4326) %>%
#   st_as_sf()
# 
# driving_combined_30 <- st_union(driving_isos_30) %>%
#   st_cast(to="POLYGON") %>%
#   st_transform(4326) %>%
#   st_as_sf()
# 
# write_rds(driving_isos_30, "clean_data/driving_isos_30_hub.rds")
# write_rds(driving_combined_30, "clean_data/driving_combined_30_hub.rds")

driving_isos_30 <- read_rds("clean_data/driving_isos_30_hub.rds")
driving_combined_30 <- read_rds("clean_data/driving_combined_30_hub.rds")

# mapbox_map <- leaflet() %>%
#   addMapboxTiles(style_id = "light-v9",
#                  username = "mapbox")
# 
# mapbox_map %>%
#   addPolygons(data = driving_combined_30,
#               stroke = FALSE,
#               color = "red",
#               fill = TRUE,
#               # fillColor = ~pal(estimate),
#               fillOpacity = 0.5,
#               dashArray = NULL,
#               smoothFactor = 1)
```

```{r inaccessible_hub_map, echo=FALSE, fig.height=11, fig.width=11, layout="l-body-outset"}

# texas_pop_block_grps_16plus <- get_acs(geography = "block group",
#                       variables = "B23027_001",
#                       state="TX",
#                       survey = "acs5",
#                       geometry = TRUE,
#                       show_call = TRUE) %>%
#   # filter(estimate==0) %>%
#   st_transform(4326) %>%
#   st_as_sf()

# write_rds(texas_pop_block_grps_16plus, "clean_data/texas_pop_block_grps_16plus.rds")

texas_pop_block_grps <-read_rds("clean_data/texas_pop_block_grps_16plus.rds")

texas_pop_blox <- texas_pop_block_grps %>%
  filter(estimate!=0)

# inaccessible_vaccine_hub <- st_difference(texas_pop_blox,st_union(driving_combined_30)) %>%
#   st_cast(to="POLYGON") %>%
#   st_transform(4326) %>%
#   st_as_sf()  %>%
#   mutate(est_lab = scales::comma(estimate, accuracy=1, suffix = " People"))
# 
# beepr::beep(1)

# write_rds(inaccessible_vaccine_hub, paste0("clean_data/", today, "_inaccessible_vaccine_hub_16plus.rds"))
# write_rds(inaccessible_vaccine_hub, "clean_data/inaccessible_vaccine_hub_16plus.rds")

inaccessible_vaccine_hub <-read_rds("clean_data/inaccessible_vaccine_hub_16plus.rds")

inaccessible_hub_summ <- inaccessible_vaccine_hub %>%
  as_tibble() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  mutate(state="TX") %>%
  group_by(state) %>%
  summarise(sum_inaccessible = sum(estimate))

# write_rds(inaccessible_hub_summ, paste0("clean_data/", today, "_inaccessible_vaccine_hub_summ.rds"))

date_pulled <- unique(vaccine_sites_sf$date_pulled)

  inaccessible_vaccine_hub %>% 
  ggplot() +
  geom_sf(data = inaccessible_vaccine_hub, aes(fill = "Inaccessible Population (+30min Drive)"), color = NA) +
  geom_sf(data = tidycensus::county_laea %>% filter(str_detect(GEOID, "^48")), fill=NA, color = "#d6d2c4") +
  geom_sf(data = driving_combined_30, aes(fill = "Accessible Population"), color = NA) +
  scale_fill_manual(values = c("Inaccessible Population (+30min Drive)"="#9cadb7","Accessible Population"="#bf5700")) +
  # guides(guide = guide_legend(override.aes = list(linetype = "blank", shape = NA))) +
  theme_minimal(base_family = "Graphik") +
  coord_sf(datum = NA, crs = 3082) +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =24),
        plot.subtitle = element_text(color = "#333f48", margin = margin(b=15), size = 16),
        plot.caption = element_text(color = "#333f48", size =12),
        # legend.key = element_rect(fill = "white", colour = NA),
        legend.position = c(.25,.15),
        legend.title = element_text(family = "Graphik-Bold", size=18),
        legend.text = element_text(size=16),
        plot.margin = margin(t=15, r=15, b=10, l = 15)) +
  labs(title = "~5.6m Texans are a 30+min drive from a Vaccine Hub",
       subtitle = paste0("Shown below are census blocks where the population aged 12 years and older would\nhave to drive 30 minutes or more to one of the state's existing Vaccine Hubs. Areas\nmarked in orange are within a 30 min drive to a Vaccine Hub.\nHub Locations Accurate as of ", date_pulled, "."),
       fill = "Legend",
       caption = "Provider Data: 'Texas Vaccine Provider Data', Texas Department of State Health Services.\nLink: https://genesis.soc.texas.gov/files/accessibility/vaccineprovideraccessibilitydata.csv\nPopulation Data: U.S. Census Bureau, ACS 5-Year Series ('B23027_001' = 12+ Population Estimate)") 
  
```

### From Pharmacies

```{r inaccessible_pharm_df, include=FALSE, layout="l-page"}

vaccine_sites_sf <- read_csv("https://utexas-lbjp-data.github.io/txcovidvaccinedata/clean_data/vaccine_provider_raw.csv") %>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>% 
    select(-Address) %>% 
    mutate(Address = paste0(Street, ", ", City, ", TX ", Zip)) %>% 
    select(-Street, -City, -Zip) %>% 
    mutate(Time = format(as.POSIXct(`Time Updated`), format = "%I:%M %p")) %>% 
    select(1:2, `Total Available`, `Last Updated`,Time, Phone, Address, everything())

vaccine_sites_sf_pharm <- vaccine_sites_sf %>% 
  filter(`Organization Type`=="Pharmacy")

# driving_isos_30 <- mb_isochrone(
#   vaccine_sites_sf_pharm,
#   profile = "driving",
#   time = 30,
#   id = "Provider Name"
# ) %>%
#   st_transform(4326) %>%
#   st_as_sf()
# 
# driving_combined_30 <- st_union(driving_isos_30) %>%
#   st_cast(to="POLYGON") %>%
#   st_transform(4326) %>%
#   st_as_sf()
# 
# write_rds(driving_isos_30, "clean_data/driving_isos_30_pharm.rds")
# write_rds(driving_combined_30, "clean_data/driving_combined_30_pharm.rds")

driving_isos_30 <- read_rds("clean_data/driving_isos_30_pharm.rds")
driving_combined_30 <- read_rds("clean_data/driving_combined_30_pharm.rds")
```

```{r inaccessible_pharm_map, echo=FALSE, fig.height=11, fig.width=11, layout="l-body-outset"}

# texas_pop_block_grps_16plus <- get_acs(geography = "block group",
#                       variables = "B23027_001",
#                       state="TX",
#                       survey = "acs5",
#                       geometry = TRUE,
#                       show_call = TRUE) %>%
#   # filter(estimate==0) %>%
#   st_transform(4326) %>%
#   st_as_sf()

# write_rds(texas_pop_block_grps_16plus, "clean_data/texas_pop_block_grps_16plus.rds")

texas_pop_block_grps <-read_rds("clean_data/texas_pop_block_grps_16plus.rds")

texas_pop_blox <- texas_pop_block_grps %>%
  filter(estimate!=0)

# inaccessible_vaccine_pharm <- st_difference(texas_pop_blox,st_union(driving_combined_30)) %>%
#   st_cast(to="POLYGON") %>%
#   st_transform(4326) %>%
#   st_as_sf()  %>%
#   mutate(est_lab = scales::comma(estimate, accuracy=1, suffix = " People"))
# 
# write_rds(inaccessible_vaccine_pharm, paste0("clean_data/", today, "_inaccessible_pharm_16plus.rds"))
# write_rds(inaccessible_vaccine_pharm, "clean_data/inaccessible_pharm_16plus.rds")

inaccessible_vaccine_pharm <-read_rds("clean_data/inaccessible_pharm_16plus.rds")

# inaccessible_hub_summ <- inaccessible_vaccine_pharm %>%
#   as_tibble() %>%
#   distinct(GEOID, .keep_all = TRUE) %>%
#   mutate(state="TX") %>%
#   group_by(state) %>%
#   summarise(sum_inaccessible = sum(estimate))

# write_rds(inaccessible_hub_summ, paste0("clean_data/", today, "_inaccessible_pharm_summ.rds"))

date_pulled <- unique(vaccine_sites_sf$date_pulled)

  inaccessible_vaccine_pharm %>% 
  ggplot() +
  geom_sf(data = inaccessible_vaccine_pharm, aes(fill = "Inaccessible Population (+30min Drive)"), color = NA) +
  geom_sf(data = tidycensus::county_laea %>% filter(str_detect(GEOID, "^48")), fill=NA, color = "#d6d2c4") +
  # geom_sf(data = driving_combined_30, aes(fill = "Accessible Population"), color = NA) 
  scale_fill_manual(values = c("Inaccessible Population (+30min Drive)"="#9cadb7","Accessible Population"="#bf5700")) +
  # guides(guide = guide_legend(override.aes = list(linetype = "blank", shape = NA))) +
  theme_minimal(base_family = "Graphik") +
  coord_sf(datum = NA, crs = 3082) +
  theme(plot.title = element_text(family = "Graphik-Bold", color = "#333f48", size =24),
        plot.subtitle = element_text(color = "#333f48", margin = margin(b=15), size = 16),
        plot.caption = element_text(color = "#333f48", size =12),
        # legend.key = element_rect(fill = "white", colour = NA),
        legend.position = c(.25,.15),
        legend.title = element_text(family = "Graphik-Bold", size=18),
        legend.text = element_text(size=16),
        plot.margin = margin(t=15, r=15, b=15, l = 15)) +
  labs(title = "~2 million Texans are a 30+min drive from a pharmacy",
       subtitle = paste0("Shown below are census blocks where the population aged 12 years and older would\nhave to drive 30 minutes or more to one of the state's existing pharmacies providing\nCOVID-19 vaccines.Areas marked in orange are within a 30 min drive to a pharmacy.\nPharmacy Locations Accurate as of ", date_pulled, "."),
       fill = "Legend",
       caption = "Provider Data: 'Texas Vaccine Provider Data', Texas Department of State Health Services.\nLink: https://genesis.soc.texas.gov/files/accessibility/vaccineprovideraccessibilitydata.csv\nPopulation Data: U.S. Census Bureau, ACS 5-Year Series ('B23027_001' = 12+ Population Estimate)") 
  
  # +
  # ggsave(paste0("figures/",today,"_","inaccessible_pharmacies_map.png"), dpi = 300, width = 13.5, height = 11.5)
  
```


## Photography Credit {.appendix}

<span>Photo by <a href="https://unsplash.com/@rymagsino?utm_source=unsplash&amp;utm_medium=referral&amp;utm_content=creditCopyText">Ryan Magsino</a> on <a href="https://unsplash.com/s/photos/austin-texas-roads?utm_source=unsplash&amp;utm_medium=referral&amp;utm_content=creditCopyText">Unsplash</a></span>