PrEP mappping

# Load necessary libraries
library(readxl)

## Warning: package 'readxl' was built under R version 4.3.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sf)

## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

library(RColorBrewer)
library(tigris)

## Warning: package 'tigris' was built under R version 4.3.3

## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

# Set options for tigris to use the sf class
options(tigris_use_cache = TRUE, tigris_class = "sf")

# Load data from the "State totals" sheet
AHEAD_data <- read_excel("C:/Users/bex/Downloads/AHEAD_data.xlsx", sheet = "State totals")

# Step 1: Filter the data for PrEP Coverage percentage
filtered_data <- AHEAD_data %>%
  filter(indicator == "PrEP coverage", statistic == "percentage") %>%  # Filter for PrEP coverage percentage
  mutate(state = toupper(trimws(state)))  # Standardize state names
# Step 2: Download the U.S. states shapefile
states <- states(cb = TRUE) %>%
  rename(state = NAME) %>%  # Rename `NAME` to `state` for merging
  mutate(state = toupper(trimws(state)))  # Standardize state names in the shapefile

## Retrieving data for the year 2021

# Step 3: Merge the filtered data with the shapefile data
map_data <- left_join(states, filtered_data, by = "state")
# Step 5: Plot PrEP Coverage map with blue gradient and adjusted limits
# This will enlarge the map view by focusing on the continental U.S.

plot(map_data["value"],
     main = "PrEP Coverage Percentage by State",
     breaks = "quantile",          # Use quantile breaks for better color distribution
     nbreaks = 9,                  # Divide into 9 bins
     pal = brewer.pal(9, "Blues"),  # Use Blue gradient palette
     xlim = c(-130, -65),          # Adjust x-limits to zoom in on the continental US
     ylim = c(20, 55)              # Adjust y-limits for the continental US
)