---
title: "Dorling Cartograms HIV Rate and Burden, 2015"
author: "Danny Sack<br><small>Division of Epidemiology<br>Vanderbilt University School of Medicine</small>"
date: "<small>`r Sys.Date()`</small>"
output:
  rmdformats::readthedown:
    code_folding: hide
    lightbox: true
    code_download: true
description: "Dorling Cartograms HIV Rate and Burden, 2015"
---

```{r setup, include=FALSE}
library(tidyverse)
library(readxl)
library(Hmisc)
library(sf)
library(tigris) # for county level data
library(cartogram)
library(cdlTools)
library(extrafont)
library(ggtext)
library(grid)
library(gtable)
library(ggspatial)
library(ggpubr)
loadfonts()
knitrSet(lang='markdown', fig.path='fin_figs2/', fig.align='center', cache=TRUE)
options(prType='html')
```

# Set Up

```{r cleaning, results = "hide"}
##### STATE LEVEL DATA SET UP #####

# read in data for this project

# read in VCCC data, variable of interest "VCCC_state"
vccc <- read_excel("HIV_VCCC_ByState_2015.xlsx")

# now pull in 2015 HIV prevalence data from CDC https://www.cdc.gov/nchhstp/atlas/index.htm
prev <- read_csv("2015_states_prev.csv")

# now join vccc with prevalence data
vccc <- left_join(vccc, prev, by = c("geography" = "Geography")) %>%
  select(-prevalent_hiv_cases, -rate_per_100000) %>%
  rename(prevalent_hiv_cases = Cases,
         rate_per_100000 = `Rate per 100000`)

# get VCCC cases included
sum(vccc$VCCC_State, na.rm = TRUE)

# read in NA-ACCORD data, variable of interest “state_pop”
naaccord <- read_excel("HIV_NAACCORD_ByState_2015.xlsx")

# get NA-ACCORD cases included
sum(naaccord$state_pop, na.rm = TRUE)

# combine key rows of VCCC and NA-ACCORD Data
state_data <- left_join(naaccord, vccc, by = "state") %>%
  select(state, cohort_id:region_pop, VCCC_State, geography, year:FIPS, prevalent_hiv_cases, rate_per_100000) # select key rows

# start by getting state polygon data
states <- states(cb = TRUE, year = 2015, class = "sf")

# now merge with vccc and na-accord data
states <- left_join(states, state_data, by = c("GEOID" = "FIPS"))

# transform spatial object for appropriate mapping
states_trans <- st_transform(states, crs = st_crs("+init=epsg:3395"))

# remove non-continental states
# "02" Alaska
# "15" Hawaii
# "60" American Samoa
# "66" Guam
# "69" Northern Mariana Islands
# "72" Puerto Rico
# "78" Virgin Islands
states_trans_fin <- states_trans %>% 
  filter(GEOID %nin% c("02", "15", "60", "66", "69", "72", "78"))

# now make categorical state numbers for plots

# VCCC state prevalence categories
states_trans_fin$VCCC_State_cat <- NA
states_trans_fin$VCCC_State_cat[states_trans_fin$VCCC_State < 5] <- "0 - 4"
states_trans_fin$VCCC_State_cat[states_trans_fin$VCCC_State > 4 & states_trans_fin$VCCC_State < 20] <- "4 - 19"
states_trans_fin$VCCC_State_cat[states_trans_fin$VCCC_State > 19 & states_trans_fin$VCCC_State < 206] <- "19 - 205"
states_trans_fin$VCCC_State_cat[states_trans_fin$VCCC_State > 205] <- "205 - 7,345"
states_trans_fin$VCCC_State_cat <- factor(states_trans_fin$VCCC_State_cat,
                                          levels = c("0 - 4", "4 - 19", "19 - 205", "205 - 7,345"))

# NA-ACCORD state prevalence categories
states_trans_fin$state_pop_cat <- NA
states_trans_fin$state_pop_cat[states_trans_fin$state_pop < 1323] <- "0 - 1,322"
states_trans_fin$state_pop_cat[states_trans_fin$state_pop > 1322 & states_trans_fin$state_pop < 4092] <- "1,322 - 4,091"
states_trans_fin$state_pop_cat[states_trans_fin$state_pop > 4091 & states_trans_fin$state_pop < 8737] <- "4,091 - 8,736"
states_trans_fin$state_pop_cat[states_trans_fin$state_pop > 8736] <- "8,736 - 31,403"
states_trans_fin$state_pop_cat <- factor(states_trans_fin$state_pop_cat,
                                          levels = c("0 - 1,322", "1,322 - 4,091", "4,091 - 8,736", "8,736 - 31,403"))

# CDC prevalence state categories
states_trans_fin$prev_cat <- NA
states_trans_fin$prev_cat[states_trans_fin$prevalent_hiv_cases < 7804] <- "0 - 7,803"
states_trans_fin$prev_cat[states_trans_fin$prevalent_hiv_cases > 7803 & states_trans_fin$prevalent_hiv_cases < 21608] <- "7,803 - 21,607"
states_trans_fin$prev_cat[states_trans_fin$prevalent_hiv_cases > 21607 & states_trans_fin$prevalent_hiv_cases < 49464] <- "21,607 - 49,463"
states_trans_fin$prev_cat[states_trans_fin$prevalent_hiv_cases > 49463] <- "49,463 - 128,681"
states_trans_fin$prev_cat <- factor(states_trans_fin$prev_cat,
                                          levels = c("0 - 7,803", "7,803 - 21,607", "21,607 - 49,463", "49,463 - 128,681"))

# CDC prelavence proportion state categories
states_trans_fin$prop_cat <- NA
states_trans_fin$prop_cat[states_trans_fin$rate_per_100000 < 216] <- "0 - 215"
states_trans_fin$prop_cat[states_trans_fin$rate_per_100000 > 215 & states_trans_fin$rate_per_100000 < 406] <- "215 - 405"
states_trans_fin$prop_cat[states_trans_fin$rate_per_100000 > 405 & states_trans_fin$rate_per_100000 < 770] <- "405 - 769"
states_trans_fin$prop_cat[states_trans_fin$rate_per_100000 > 769] <- "769 - 2,590"
states_trans_fin$prop_cat <- factor(states_trans_fin$prop_cat,
                                          levels = c("0 - 215", "215 - 405", "405 - 769", "769 - 2,590"))

# get total hiv cases
sum(states_trans_fin$prevalent_hiv_cases)
# get summary of proportion
summary(states_trans_fin$rate_per_100000)

# make Gastner-Newman Cartograms
gn_prev <- cartogram_cont(states_trans_fin, "prevalent_hiv_cases", itermax = 30)
gn_prop <- cartogram_cont(states_trans_fin, "rate_per_100000", itermax = 30)

# make Dorling Cartograms
d_prev <- cartogram_dorling(states_trans_fin, "prevalent_hiv_cases", k = 9) # can fiddle with k from 5-10, all look okay
d_prop <- cartogram_dorling(states_trans_fin, "rate_per_100000", k = 9) 

##### COUNTY LEVEL DATA SET UP #####
# read in county-level data, needed to delete the first few rows in excel to make it work
vccc_cty <- read_excel("HIV_VCCC_ByCounty_2015.xlsx")

# get county level map data of continguous states from 2015 using tigris package
usa_county <- counties(state = vccc$state, year = 2015, class = "sf")

# subset hiv_data by matching fips codes
vccc_cty <- vccc_cty %>% filter(FIPS %in% usa_county$GEOID)

# merge HIV info from hiv_data to maps
usa_county <- left_join(usa_county, vccc_cty, by = c("GEOID" = "FIPS"))

# now do it with all counties from CDC https://www.cdc.gov/nchhstp/atlas/index.htm
cty_2015 <- read_csv("2015_counties_prev.csv")

# if data suppressed or not available, make 0
cty_2015 <- cty_2015 %>% mutate(Cases = parse_number(Cases, na = c("Data suppressed")),
                                Cases = ifelse(!is.na(Cases), Cases, 0),
                                `Rate per 100000` = parse_number(`Rate per 100000`, na = c("Data suppressed")),
                                `Rate per 100000` = ifelse(!is.na(`Rate per 100000`), `Rate per 100000`, 0))

# now combine with usa_county
all_cty <- usa_county %>% mutate(GEOID = parse_number(GEOID)) %>%
  left_join(., cty_2015, by = c("GEOID" = "FIPS"))

# now make Cases and Rates Zero if missing and pt_count_byCounty zero if missing
all_cty <- all_cty %>% mutate(Cases = ifelse(is.na(Cases), 0, Cases),
                              `Rate per 100000` = ifelse(is.na(`Rate per 100000`), 0, `Rate per 100000`),
                              pt_count_byCounty = ifelse(is.na(pt_count_byCounty), 0, pt_count_byCounty))

# now change coordinate system
all_cty_trans <- st_transform(all_cty, crs = st_crs("+init=epsg:3395"))

# now remove non-contiguous states
all_cty_trans <- all_cty_trans %>% filter(STATEFP %nin% c("02", "15", "60", "66", "69", "72", "78"))

# add new state label 
all_cty_trans <- all_cty_trans %>% mutate(state_abb = fips(STATEFP, "Abbreviation"))

# now add labels
all_cty_trans <- all_cty_trans %>% 
  group_by(state_abb) %>% 
  summarise(max_prev = max(Cases)) %>% 
  mutate(prev_lab = state_abb) %>% st_drop_geometry() %>%
  left_join(all_cty_trans, ., by = c("state_abb" = "state_abb", "Cases" = "max_prev"))

# add labels on counties above the 97.5 percentile of prevalent HIV
all_cty_trans <- all_cty_trans %>%
  mutate(prev_lab = ifelse(Cases > quantile(Cases, 0.975, na.rm = TRUE), state_abb, NA))

# now adding back counties with missing/0 cases and rates make them all 1 to avoid incorrect scaling that distorts output
all_cty_trans <- all_cty_trans %>% mutate(Cases1 = ifelse(Cases == 0, 1, Cases),
                                          Rate = ifelse(`Rate per 100000` == 0, 1, `Rate per 100000`))

# make county categories for vccc cases
all_cty_trans$cat_vccc_pts <- NA
all_cty_trans$cat_vccc_pts[all_cty_trans$pt_count_byCounty < 5] <- "0 - 4"
all_cty_trans$cat_vccc_pts[all_cty_trans$pt_count_byCounty > 4 & all_cty_trans$pt_count_byCounty < 20] <- "4 - 19"
all_cty_trans$cat_vccc_pts[all_cty_trans$pt_count_byCounty > 19 & all_cty_trans$pt_count_byCounty < 206] <- "19 - 205"
all_cty_trans$cat_vccc_pts[all_cty_trans$pt_count_byCounty > 205] <- "205 - 2,981"
all_cty_trans$cat_vccc_pts <- factor(all_cty_trans$cat_vccc_pts, 
                                     levels = c("0 - 4", "4 - 19", "19 - 205", "205 - 2,981"))


# make county dorling map based on county-level HIV prevalence
burden_dorling_prev <- cartogram_dorling(all_cty_trans, "Cases1", k = 2)

# subset to just the state of TN
tn_cty_trans <- all_cty_trans %>%
  filter(STATEFP == "47") %>%
  mutate(tn_lab = ifelse(NAME %in% c("Davidson", "Shelby"), NAME, ""))

# dorling cartogram for tennessee based on HIV prevalence
dorlingtn <- cartogram_dorling(tn_cty_trans, "Cases1", k = 3)

# gastner-newman cartogram for tennessee based on HIV prevalence
gncartn <- cartogram_cont(tn_cty_trans, "Cases1", itermax = 30)

## previous contiguous cartograms are from Dougenik, J. A., Chrisman, N. R., & Niemeyer, D. R. (1985). An Algorithm To Construct Continuous Area Cartograms. In The Professional Geographer, 37(1), 75-81.
## previous Dorling cartograms are from Dorling, D. (1996). Area Cartograms: Their Use and Creation. In Concepts and Techniques in Modern Geography (CATMOG), 59.

# follow directions here for Gastner-Newman cartograms, another type of contiguous cartograms
# https://www.r-bloggers.com/2016/10/cartograms-with-r/
# does not work on M1 macs yet
# requires link between R installation and "fftw", which can be installed via homebrew
##Install the R implementation of Cart by Gastner and Newman (2004)
# devtools::install_github("omegahat/Rcartogram")
# devtools::install_github('chrisbrunsdon/getcartr',subdir='getcartr')
```

# Tennessee Plots

```{r fig1, fig.width=7, fig.height=7.5}
# original width 14, height 15, text size 14, title size 18

# VCCC coordinates
# lat = 36.10868177979694, long = -86.76333513103516
vccclat <- 36.10868177979694
vccclong <- -86.76333513103516

# theme for figure 1
theme_map_fig1 <- theme_pubr() + 
  theme(text = element_text(family = "Arial", size = 9),
        axis.text.x = element_blank(), legend.text = element_text(face = "bold"),
        legend.title = element_text(face = "bold"), 
        legend.background = element_rect(fill = NA),
        plot.title = element_markdown(lineheight = 1.1, size = 12),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(), 
        axis.ticks.y = element_blank(), axis.title.x = element_blank(), 
        axis.title.y = element_blank(), axis.line = element_blank(),
        panel.background = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        plot.caption = element_text(hjust = 0.05, color = "#993D1B", 
                                    face = "bold", family = "Arial"))

# names of counties to be labeled on all plots
c_names <- c("Davidson", "Shelby", "Hamilton", "Knox") 

# add columns to for End the Epidemic County
tn_cty_trans <- tn_cty_trans %>%
  mutate(Name = ifelse(NAME == "Shelby", "EndEpi", NA))

# chloropleth with color by VCCC
#### will annotate this figure in detail, all other figures follow similar pattern ####
tn_chlor_vccc <- ggplot(tn_cty_trans) + # initializes plot with tn_cty_trans data
  geom_sf(aes(fill = cat_vccc_pts, color = Name)) + # adds shapefile geometry, with fill by VCCC patient quanity, colored by county name
  geom_point(x = vccclong, y = vccclat, 
             size = 1.5, color = "#E7BC51") + # add point for VCCC clinic
  geom_sf_text(aes(label = ifelse(NAME %in% c_names, NAME, "")), 
               size = 2.5, fontface = "bold", color = "white", family = "Arial") + #add layer with county names
  annotate(geom = "text", x = vccclong, y = vccclat - 0.12, 
           label = "VCCC", size = 2, fontface = "bold", 
           color = "#E7BC51", family = "Arial") + # adds VCCC clinic label on plot
  labs(x=NULL, y=NULL) + # removes coordinate axese
  ggtitle("<b>A) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort\n(VCCC)</b><br>
          <b style='color:#C7D9EE'>A) </b><b style='color:#080F8C'>Deformed by:</b>
          <b>Not Applicable</b>") + # adds title using ggtext package to allow element_markdown in theme
  scale_fill_manual(name = "VCCC\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.6,
                                         keyheight = 0.6)) + # manually scale county fill
  scale_color_manual(values = c("#993D1B", "#000000"),
                     guide = "none") + # manually scale county color
  annotation_scale(location = "tl", height = unit(0.05, "cm"), 
                   width_hint = 0.15,
                   text_family = "Arial") + # add map distance scale for choropleth plot
  annotation_north_arrow(location = "tl", pad_y = unit(0.4, "cm"),
                         height = unit(0.5, "cm"), width = unit(0.5, "cm"),
                         style = north_arrow_orienteering(text_family = "Arial",
                                                     text_size = 8)) + # add north arrow to orient reader
  labs(caption = "Priority Jurisdiction\nEnding the HIV Epidemic") + # add a caption to label Ending the Epidemic County
  annotation_custom(rectGrob(gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE")),
                    xmin = -81.25, xmax = -90.8,
                    ymin = 36.85, ymax = 37.35) + # add panel header
  coord_sf(clip = 'off', default_crs = st_crs(4326)) + # set coordinate system for plot
  theme_map_fig1 + # add overall map theme
  theme(legend.position = c(0.90, 0.18)) # make slight adjustment to legend position

# county names for GN cartogram
# names in white text
gn_names_w <- c("Rutherford", "Davidson", "Sumner", 
              "Maury", "Putnam",
              "Madison", "Shelby", "Hamilton", "Knox") 
# names in black text              
gn_names_b <- c("Washington")

# add columns for End the Epidemic County
gncartn <- gncartn %>%
  mutate(Name = ifelse(NAME == "Shelby", "EndEpi", NA))

# GN deformed by prevlent HIV cases, colored by VCCC
tn_gn_vccc <- ggplot(gncartn) +
  geom_sf(aes(fill = cat_vccc_pts, color = Name)) +
  geom_sf_text(aes(label = ifelse(NAME %in% gn_names_w, NAME, "")), 
               size = 3, fontface = "bold", color = "white", family = "Arial") +
  geom_sf_text(aes(label = ifelse(NAME %in% gn_names_b, NAME, "")), 
               size = 3, fontface = "bold", color = "black", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>B) </b><b style='color:#2A3389'>Colored by:</b>
          <b>Cumulative Cohort (VCCC)</b><br>
          <b style='color:#C7D9EE'>B) </b><b style='color:#2A3389'>Deformed by:</b>
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "VCCC\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C")) +
  scale_color_manual(values = c("#993D1B", "#000000"),
                     guide = "none") +
  guides(fill = guide_legend(keywidth = 0.6,
                             keyheight = 0.6)) +
  labs(caption = "Priority Jurisdiction\nEnding the HIV Epidemic") +
  annotation_custom(rectGrob(y = 1.174, height = 0.25,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig1 +
  theme(legend.position = c(0.88, 0.20))

# add columns for End the Epidemic County
dorlingtn <- dorlingtn %>%
  mutate(Name = ifelse(NAME == "Shelby", "EndEpi", NA))

# Dorling deformed by prevlent HIV cases, colored by VCCC
tn_dor_vccc <- ggplot(dorlingtn) +
  geom_sf(aes(fill = cat_vccc_pts, color = Name)) +
  geom_sf_text(aes(label = ifelse(NAME %in% c("Davidson", "Shelby", "Knox"), NAME, "")), 
               size = 2, fontface = "bold", color = "white", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>C) </b><b style='color:#2A3389'>Colored by:</b>
          <b>Cumulative Cohort (VCCC)</b><br>
          <b style='color:#C7D9EE'>C) </b><b style='color:#2A3389'>Deformed by:</b>
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "VCCC\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C")) +
  scale_color_manual(values = c("#993D1B", "#000000"),
                     guide = "none") +
  guides(fill = guide_legend(keywidth = 0.6,
                             keyheight = 0.6)) +
  labs(caption = "Priority Jurisdiction\nEnding the HIV Epidemic") +
  annotation_custom(rectGrob(y = 1.177, height = 0.27,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig1 +
  theme(legend.position = c(0.88, 0.20))

# combine plots to make the same width
# save plots as grobs
tn1 <- ggplotGrob(tn_chlor_vccc)
tn2 <- ggplotGrob(tn_gn_vccc)
tn3 <- ggplotGrob(tn_dor_vccc)

# get widths to be the same
# bind as row
tn4 <- rbind(tn1, tn2, tn3, size = "first")
tn4 $widths <- unit.pmax(tn1$widths, tn2$widths, tn3$widths)
# add space between rows
tn4 <- gtable_add_row_space(tn4, unit(0.01, "cm"))
# output
pdf("Final Figure PDFs/fig1.pdf", height = 7.5, width = 7)
grid.newpage()
grid.draw(tn4)
dev.off()
# embed the font
embed_fonts("Final Figure PDFs/fig1.pdf")
```

# County Level Plots

```{r fig2, fig.width=7, fig.height=5.5}
# create vector of end the epidemic counties
endepicty <- c("Maricopa", # AZ
               "Alameda", "Los Angeles", "Orange", "Riverside",
               "Sacramento", "San Bernardino", "San Diego", "San Francisco", #CA
               "Broward", "Duval", "Hillsborough", "Miami-Dade",
               "Orange", "Palm Beach", "Pinellas", # FL,
               "Cobb", "DeKalb", "Fulton", "Gwinnett", # GA,
               "Cook", # IL,
               "Marion", #IN,
               "East Baton Rouge", "Orleans", #LA
               "Baltimore", "Montgomery", "Prince George's", #MD
               "Suffolk", #MA
               "Wayne", #MI
               "Clark", #NV
               "Essex", "Hudson", #NJ
               "Bronx", "Kings", "New York", "Queens", #NY
               "Mecklenburg", #NC
               "Cuyahoga", "Franklin", "Hamilton", #OH
               "Philadelphia", #PA
               "Shelby", #TN
               "Bexar", "Dallas", "Harris", "Tarrant", "Travis", #TX
               "King") #WA

# states where specific counties are part of end the epidemic phase 1
endepisc <- c("AZ", rep("CA", 8), rep("FL", 7),
              rep("GA", 4), "IL", "IN", rep("LA", 2),
              rep("MD", 3), "MA", "MI", "NV",
              rep("NJ", 2), rep("NY", 4), "NC",
              rep("OH", 3), "PA", "TN",
              rep("TX", 5), "WA")

# combine into tibble
endepitib <- tibble(NAME = endepicty, state_abb = endepisc)

# now state_abb where all counties are in end the epidemic
state_endepi <- c("DC", "AL", "AR", "KY", "MS", "MO", "OK", "SC")

# now add an end the epidemic colum to county-level dorling cartogram
# start with full states
burden_dorling_prev <- burden_dorling_prev %>%
  mutate(end_epi = ifelse(state_abb %in% state_endepi, "Yes", "No"),
         end_epi = factor(end_epi, levels = c("Yes", "No")))

# now do each county
for(i in 1:nrow(endepitib)){
  burden_dorling_prev$end_epi[burden_dorling_prev$NAME == endepitib$NAME[i] & burden_dorling_prev$state_abb == endepitib$state_abb[i]] <- "Yes"
}

# theme for figure 2
theme_map_fig2 <- theme_pubr() + 
  theme(legend.position = c(0.25, 0.1), text = element_text(family = "Arial", size = 9),
        axis.text.x = element_blank(), legend.text = element_text(face = "bold"),
        legend.title = element_text(face = "bold"),
        legend.background = element_rect(fill = NA),
        plot.title = element_markdown(lineheight = 1.1, size = 12),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(), 
        axis.ticks.y = element_blank(), axis.title.x = element_blank(), 
        axis.title.y = element_blank(), axis.line = element_blank(),
        panel.background = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

# Dorling all counties deformed by prevlent HIV cases, colored by VCCC
count_prev <- ggplot(burden_dorling_prev) +
  geom_sf(aes(fill = cat_vccc_pts, color = state_abb, size = end_epi)) +
  geom_sf_text(aes(label = prev_lab), size = 2, fontface = "bold", family = "Arial") +
  geom_sf_text(aes(label = ifelse(cat_vccc_pts %in% c("205 - 2,981"), state, "")), 
               size = 2, fontface = "bold", color = "white", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b style='color:#2A3389'>Colored by:</b>
          <b>Cumulative Cohort (VCCC)</b><br>
          <b style='color:#2A3389'>Deformed by:</b>
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "VCCC\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C")) +
  scale_size_manual(name = "Priority Jurisdiction\nEnding the HIV Epidemic", 
                    values = c(1, 0.2),
                    guide = guide_legend(override.aes = list(fill = NA),
                                         keywidth = 0.8, keyheight = 0.8)) +
  scale_color_viridis_d(name = "State") +
  annotation_custom(rectGrob(y = 1.065, height = 0.1,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  guides(color = "none", fill = guide_legend(order = 1,
                                             keywidth = 0.8, 
                                             keyheight = 0.8)) +
  theme_map_fig2 +
  theme(legend.box = "horizontal",
        legend.box.just = "bottom")

# repeat dorling all counties deformed by prevlent HIV cases, colored by VCCC
# for this version, change orientation of legend, so state legend color can be horizontal
state_leg <- ggplot(burden_dorling_prev) +
  geom_sf(aes(fill = cat_vccc_pts, color = state_abb)) +
  geom_sf_text(aes(label = prev_lab), size = 2, fontface = "bold", family = "Arial") +
  geom_sf_text(aes(label = ifelse(cat_vccc_pts %in% c("205 - 2,981"), state, "")), 
               size = 2, fontface = "bold", color = "white", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort (VCCC)</b><br>
          <b style='color:#080F8C'>Deformed by:
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "VCCC\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C")) +
  scale_color_viridis_d(name = "State", guide = guide_legend(nrow = 4, title.position = "top",
                              title.hjust = 0.5, keywidth = 0.6,
                              keyheight = 0.6,
                              override.aes = list(size = 3, fill = NA))) +
  guides(fill = "none") +
  theme_map_fig2 +
  theme(legend.position = "bottom")

# pull horizontal state legend into it's own grob
leg <- get_legend(state_leg)

# write out
pdf("Final Figure PDFs/fig2.pdf", height = 5.5, width = 7)
# combine full map with horizontal state legend
gridExtra::grid.arrange(count_prev, leg, ncol = 1, heights = c(1, 0.2))
dev.off()
embed_fonts("Final Figure PDFs/fig2.pdf")
```


# State Level Plots

```{r fig3_prep, fig.width=7, fig.height=5}
# theme for figure 3
theme_map_fig3 <- theme_pubr() + 
  theme(legend.position = c(0.15, 0.15), text = element_text(family = "Arial", size = 8),
        axis.text.x = element_blank(), legend.text = element_text(face = "bold"),
        legend.title = element_text(face = "bold"),
        legend.background = element_rect(fill = NA),
        plot.title = element_markdown(lineheight = 1.1, size = 9),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(), 
        axis.ticks.y = element_blank(), axis.title.x = element_blank(), 
        axis.title.y = element_blank(), axis.line = element_blank(),
        panel.background = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

# Choropleth colored by prevlent cases
choro_prev <- ggplot(states_trans_fin) +
  geom_sf(aes(fill = prev_cat)) +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>A) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Prevalent HIV+ Population (CDC)</b><br>
          <b style='color:#C7D9EE'>A) </b><b style='color:#080F8C'>Deformed by:</b>
          <b>Not Applicable</b>") +
  scale_fill_manual(name = "CDC\nPrev. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.5, keyheight = 0.5)) +
  annotation_scale(location = "br", height = unit(0.05, "cm"), 
                   pad_y = unit(0, "cm"),
                   width_hint = 0.15,
                   text_family = "Arial") +
  annotation_north_arrow(location = "br", pad_y = unit(0.25, "cm"),
                         height = unit(0.5, "cm"), width = unit(0.5, "cm"),
                         style = north_arrow_orienteering(text_family = "Arial",
                                                     text_size = 8)) +
  annotation_custom(rectGrob(y = 1.13, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3

# Choropleth colored by prevlanet proportion
choro_prop <- ggplot(states_trans_fin) +
  geom_sf(aes(fill = prop_cat)) +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>B) </b><b style='color:#080F8C'>Colored by:</b>
          <b>HIV Prevalence Proportion (CDC)</b><br>
          <b style='color:#C7D9EE'>B) </b><b style='color:#080F8C'>Deformed by:</b>
          <b>Not Applicable</b>") +
  scale_fill_manual(name = "CDC\nHIV Cases per 100,000", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.5, keyheight = 0.5)) +
  annotation_scale(location = "br", height = unit(0.05, "cm"), 
                   pad_y = unit(0, "cm"),
                   width_hint = 0.15,
                   text_family = "Arial") +
  annotation_north_arrow(location = "br", pad_y = unit(0.25, "cm"),
                         height = unit(0.5, "cm"), width = unit(0.5, "cm"),
                         style = north_arrow_orienteering(text_family = "Arial",
                                                     text_size = 8)) +
  annotation_custom(rectGrob(y = 1.13, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3

# GN states deformed by prevlent HIV cases, colored by NA-ACCORD
na_accord_gn_prev <- ggplot(gn_prev) +
  geom_sf(aes(fill = state_pop_cat)) +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>C) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort (NA-ACCORD)</b><br>
          <b style='color:#C7D9EE'>C) </b><b style='color:#080F8C'>Deformed by:</b>
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "NA-ACCORD\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.5, keyheight = 0.5)) +
  annotation_custom(rectGrob(y = 1.12, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3

# GN states deformed by prevlent proportion, colored by NA-ACCORD
na_accord_gn_prop <- ggplot(gn_prop) +
  geom_sf(aes(fill = state_pop_cat)) +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>D) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort (NA-ACCORD)</b><br>
          <b style='color:#C7D9EE'>D) </b><b style='color:#080F8C'>Deformed by:</b>
          <b>HIV Prevalence Proportion (CDC)</b>") +
  scale_fill_manual(name = "NA-ACCORD\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.5, keyheight = 0.5)) +
  annotation_custom(rectGrob(y = 1.12, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3

# Dorling states deformed by prevlent HIV cases, colored by NA-ACCORD
na_accord_d_prev <- ggplot(d_prev) +
  geom_sf(aes(fill = state_pop_cat)) +
  geom_sf_text(aes(label = state), 
               size = 2, fontface = "bold", family = "Arial") +
  geom_sf_text(aes(label = ifelse(state_pop_cat %in% c("4,091 - 8,736", "8,736 - 31,403"), state, "")), 
               size = 2, fontface = "bold", color = "white", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>E) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort (NA-ACCORD)</b><br>
          <b style='color:#C7D9EE'>E) </b><b style='color:#080F8C'>Deformed by:</b> 
          <b>Prevalent HIV+ Population (CDC)</b>") +
  scale_fill_manual(name = "NA-ACCORD\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.6, keyheight = 0.6)) +
  annotation_custom(rectGrob(y = 1.125, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3

# Dorling states deformed by prevlent proportion, colored by NA-ACCORD
na_accord_d_prop <- ggplot(d_prop) +
  geom_sf(aes(fill = state_pop_cat)) +
  geom_sf_text(aes(label = state), 
               size = 2, fontface = "bold", family = "Arial") +
  geom_sf_text(aes(label = ifelse(state_pop_cat %in% c("4,091 - 8,736", "8,736 - 31,403"), state, "")), 
               size = 2, fontface = "bold", color = "white", family = "Arial") +
  labs(x=NULL, y=NULL) +
  ggtitle("<b>F) </b><b style='color:#080F8C'>Colored by:</b>
          <b>Cumulative Cohort (NA-ACCORD)</b><br>
          <b style='color:#C7D9EE'>F) </b><b style='color:#080F8C'>Deformed by:</b> 
          <b>HIV Prevalence Proportion (CDC)</b>") +
  scale_fill_manual(name = "NA-ACCORD\nCum. HIV Cases", 
                    values = c("#B3E0F2", "#7A98D4", "#3F4EA6", "#080F8C"),
                    guide = guide_legend(keywidth = 0.6, keyheight = 0.6)) +
  annotation_custom(rectGrob(y = 1.125, height = 0.15,
                             gp = gpar(fill = "#C7D9EE",
                                       col = "#C7D9EE"))) +
  coord_sf(clip = 'off') +
  theme_map_fig3
```

```{r fig3, fig.width = 7, fig.height = 7.5}
# original figure height 21, width 30
# save plots as grobs for combining plots
pana <- ggplotGrob(choro_prev)
panb <- ggplotGrob(choro_prop)
panc <- ggplotGrob(na_accord_gn_prev)
pand <- ggplotGrob(na_accord_gn_prop)
pane <- ggplotGrob(na_accord_d_prev)
panf <- ggplotGrob(na_accord_d_prop)


# get heights to be the same
r1 <- cbind(pana, panb, size = "first")
r1$heights <- unit.pmax(pana$heights, panb$heights)
r2 <- cbind(panc, pand, size = "first")
r2$heights <- unit.pmax(panc$heights, pand$heights)
r3 <- cbind(pane, panf, size = "first")
r3$heights <- unit.pmax(pane$heights, panf$heights)
# bind as row
all3 <- rbind(r1, r2, r3)
all3$widths<- unit.pmax(r1$widths, r2$widths, r3$widths)
# add space between rows
fig3 <- gtable_add_row_space(all3, unit(0.01, "cm"))
# output
pdf("Final Figure PDFs/fig3.pdf", height = 7.5, width = 7)
grid.newpage()
grid.draw(fig3)
dev.off()
embed_fonts("Final Figure PDFs/fig3.pdf")
```

# Session Info

```{r sessioninfo}
sessionInfo()
```