CODEWITHSUNDEEP
rzipcode
Preston Carey
Preston Carey

Reputation: 11

Trying to make a heatmap based on percentage values from zip code area

I've been trying to use the choroplethrZip package to do this based on a previous post, but my version of R will not support this package. I'm using the data set and code below, and believe I have it formatted correctly.

region value
19102 0.33
19103 0.0625
19106 0.25
19107 0.119
19111 0
19114 0
19115 0.476
19116 0.0909
19120 0.0556
19121 0.107
19122 0.158
19124 0.1
19125 0.133
19128 0.05
19129 0
19130 0.048
19131 0
19132 0.158
19133 0.0667
19134 0.108
19135 0.2
19136 0.08
19137 0
19138 0.037
19139 0
19140 0.048
19141 0
19142 0.0909
19143 0.027
19144 0.174
19145 0.108
19146 0.079
19147 0.048
19148 0.156
19149 0
19150 0.045
19152 0
19154 0 
library(choroplethrZip)

coag20 <- read_excel("comorbheatmap.xlsx", sheet = 1)

zip_choropleth(coag20, zip_zoom = (coag20$region),
               num_colors = 1, 
               title = "Comorbidity Documentation per Encounter by Zip Code: Coagulopathy",
               legend = "Percentage")

Upvotes: 0

Views: 49

Answers (1)

L Tyrone
L Tyrone

Reputation: 7075

For a more future-proof option, you can use tigris for the zipcode geometries, and ggplot2 for plotting:

library(tigris)
library(sf)
library(dplyr)
library(ggplot2)

# Your example data
coag20 <- read.table(text = "region     value
19102   0.33
19103   0.0625
19106   0.25
19107   0.119
19111   0
19114   0
19115   0.476
19116   0.0909
19120   0.0556
19121   0.107
19122   0.158
19124   0.1
19125   0.133
19128   0.05
19129   0
19130   0.048
19131   0
19132   0.158
19133   0.0667
19134   0.108
19135   0.2
19136   0.08
19137   0
19138   0.037
19139   0
19140   0.048
19141   0
19142   0.0909
19143   0.027
19144   0.174
19145   0.108
19146   0.079
19147   0.048
19148   0.156
19149   0
19150   0.045
19152   0
19154   0", header = TRUE)

# Load sf of all zipcodes starting with the same two digits of your area of interest
sf_zip <- zctas(starts_with = "19", class = "sf")

# Subset sf_zip to just the zipcodes in coag20, join values to sf
sf_coag20 <- filter(sf_zip, ZCTA5CE20 %in% coag20$region) |>
  mutate(ZCTA5CE20 = as.integer(ZCTA5CE20)) |>
  left_join(coag20, by = join_by(ZCTA5CE20 == region))

sf_coag20[,c(1,10)]
# Simple feature collection with 38 features and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: -75.26443 ymin: 39.89184 xmax: -74.95576 ymax: 40.13799
# Geodetic CRS:  NAD83
# First 10 features:
#    ZCTA5CE20  value                       geometry
# 1      19143 0.0270 MULTIPOLYGON (((-75.25165 3...
# 2      19137 0.0000 MULTIPOLYGON (((-75.09842 3...
# 3      19129 0.0000 MULTIPOLYGON (((-75.20801 4...
# 4      19115 0.4760 MULTIPOLYGON (((-75.07456 4...
# 5      19103 0.0625 MULTIPOLYGON (((-75.18325 3...
# 6      19142 0.0909 MULTIPOLYGON (((-75.24776 3...
# 7      19131 0.0000 MULTIPOLYGON (((-75.24977 3...
# 8      19107 0.1190 MULTIPOLYGON (((-75.1653 39...
# 9      19124 0.1000 MULTIPOLYGON (((-75.12714 4...
# 10     19133 0.0667 MULTIPOLYGON (((-75.1549 39...

# Plot
ggplot() +
  geom_sf(data = sf_coag20, aes(fill = value)) +
  scale_fill_continuous(name = "Percentage") +
  labs(title = "Comorbidity Documentation per Encounter by Zip Code: Coagulopathy") +
  theme(title = element_text(size = 8))

1

Upvotes: 1

Related Questions