2  Map Neighborhoods

When considering the health of persons, we have to also consider the neighborhood environment. Sometimes this is looking at neighborhood level health outcomes, like premature mortality at the census tract scale, or cumulative COVID rates by zip code. Sometimes we’re interested in neighborhood factors like poverty, access to affordable housing, or distance to nearest health provider, or pollution-emitting facility. These measurements of the “social determinants of health” at the neighborhood scale are increasingly urgent in modern public health thinking, and are thought to drive and/or reinforce racial, social, and spatial inequity

In this module, we’ll learn about the basics of thematic mapping – known as choropleth mapping – to visualize neighborhood level health phenomena. This will allow you to begin the process of exploratory spatial data analysis and hypothesis generation & refinement.

2.1 Clean Attribute Data

Let’s consider COVID-19 cases by zip code in Chicago. We’ll upload and inspect a summary of cases from the Chicago Data Portal first:

COVID <- read.csv("~/Documents/GitHub/Spatial-Health-Workshop/data/COVID-19_Cases__Tests__and_Deaths_by_ZIP_Code.csv")

head(COVID)

  ZIP.Code Week.Number Week.Start   Week.End Cases...Weekly Cases...Cumulative
1    60603          39 09/20/2020 09/26/2020              0                 13
2    60604          39 09/20/2020 09/26/2020              0                 31
3    60611          16 04/12/2020 04/18/2020              8                 72
4    60611          15 04/05/2020 04/11/2020              7                 64
5    60615          11 03/08/2020 03/14/2020             NA                 NA
6    60603          10 03/01/2020 03/07/2020             NA                 NA
  Case.Rate...Weekly Case.Rate...Cumulative Tests...Weekly Tests...Cumulative
1                  0                 1107.3             25                327
2                  0                 3964.2             12                339
3                 25                  222.0            101                450
4                 22                  197.4             59                349
5                 NA                     NA              6                  9
6                 NA                     NA              0                  0
  Test.Rate...Weekly Test.Rate...Cumulative Percent.Tested.Positive...Weekly
1               2130                27853.5                              0.0
2               1534                43350.4                              0.0
3                312                 1387.8                              0.1
4                182                 1076.3                              0.1
5                 14                   21.7                               NA
6                  0                    0.0                               NA
  Percent.Tested.Positive...Cumulative Deaths...Weekly Deaths...Cumulative
1                                  0.0               0                   0
2                                  0.1               0                   0
3                                  0.2               0                   0
4                                  0.2               0                   0
5                                   NA               0                   0
6                                   NA               0                   0
  Death.Rate...Weekly Death.Rate...Cumulative Population   Row.ID
1                   0                       0       1174 60603-39
2                   0                       0        782 60604-39
3                   0                       0      32426 60611-16
4                   0                       0      32426 60611-15
5                   0                       0      41563 60615-11
6                   0                       0       1174 60603-10
             ZIP.Code.Location
1 POINT (-87.625473 41.880112)
2 POINT (-87.629029 41.878153)
3 POINT (-87.620291 41.894734)
4 POINT (-87.620291 41.894734)
5 POINT (-87.602725 41.801993)
6 POINT (-87.625473 41.880112)

Each row corresponds to a zip code at a different week. This data thus exists as a “long” format, which doesn’t work for spatial analysis. We need to convert to “wide” format, or at the very least, ensure that each zip code corresponds to one row.

To simplify, let’s identify the last week of the dataset, and then subset the data frame to only show that week. We will be interested in the cumulative case rate. Following is one way of doing this – can you think of another way? Try out different approaches of reshaping data to test your R and “tidy” skills.

## How many weeks are in our dataset?
range(as.numeric(COVID$Week.Number))

[1] 10 40

## Subset & inspect to week 39.
COVID.39 <- subset(COVID, COVID$Week.Number == "39")
head(COVID.39)

   ZIP.Code Week.Number Week.Start   Week.End Cases...Weekly Cases...Cumulative
1     60603          39 09/20/2020 09/26/2020              0                 13
2     60604          39 09/20/2020 09/26/2020              0                 31
36    60601          39 09/20/2020 09/26/2020              8                213
37    60602          39 09/20/2020 09/26/2020              0                 21
41    60605          39 09/20/2020 09/26/2020             12                391
66    60610          39 09/20/2020 09/26/2020             35                666
   Case.Rate...Weekly Case.Rate...Cumulative Tests...Weekly Tests...Cumulative
1                   0                 1107.3             25                327
2                   0                 3964.2             12                339
36                 54                 1451.4            202               4304
37                  0                 1688.1             27                460
41                 44                 1420.8            291               7160
66                 90                 1706.9            500              10680
   Test.Rate...Weekly Test.Rate...Cumulative Percent.Tested.Positive...Weekly
1                2130                27853.5                              0.0
2                1534                43350.4                              0.0
36               1376                29328.8                              0.0
37               2170                36977.5                              0.0
41               1058                26018.4                              0.0
66               1281                27371.3                              0.1
   Percent.Tested.Positive...Cumulative Deaths...Weekly Deaths...Cumulative
1                                   0.0               0                   0
2                                   0.1               0                   0
36                                  0.0               1                   6
37                                  0.0               0                   0
41                                  0.1               1                   3
66                                  0.1               0                  10
   Death.Rate...Weekly Death.Rate...Cumulative Population   Row.ID
1                  0.0                     0.0       1174 60603-39
2                  0.0                     0.0        782 60604-39
36                 6.8                    40.9      14675 60601-39
37                 0.0                     0.0       1244 60602-39
41                 3.6                    10.9      27519 60605-39
66                 0.0                    25.6      39019 60610-39
              ZIP.Code.Location
1  POINT (-87.625473 41.880112)
2  POINT (-87.629029 41.878153)
36 POINT (-87.622844 41.886262)
37 POINT (-87.628309 41.883136)
41 POINT (-87.623449 41.867824)
66   POINT (-87.63581 41.90455)

To clean our data a bit, we’ll just keep the zip code name, and cumulative case rate for the week of September 20th, 2020.

COVID.39f <- COVID.39[,c("ZIP.Code", "Case.Rate...Cumulative")]

head(COVID.39f)

   ZIP.Code Case.Rate...Cumulative
1     60603                 1107.3
2     60604                 3964.2
36    60601                 1451.4
37    60602                 1688.1
41    60605                 1420.8
66    60610                 1706.9

Let’s rename the column headings to make them a bit more friendly.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

COVID.39f <- COVID.39f %>% 
  rename(ZIP.Code = ZIP.Code, 
         CaseRtWk39 = Case.Rate...Cumulative)

head(COVID.39f)

   ZIP.Code CaseRtWk39
1     60603     1107.3
2     60604     3964.2
36    60601     1451.4
37    60602     1688.1
41    60605     1420.8
66    60610     1706.9

2.2 Merge Spatial Data

Next, let’s merge this data to our zip code master spatial file. Reload if necessary:

library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

Chi_Zips = st_read("~/Documents/GitHub/Spatial-Health-Workshop/data/chicago_zips1.shp")

Reading layer `chicago_zips1' from data source 
  `/Users/ppgroup/Documents/GitHub/Spatial-Health-Workshop/data/chicago_zips1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 59 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
Geodetic CRS:  +proj=longlat +ellps=WGS84 +no_defs

head(Chi_Zips)

Simple feature collection with 6 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.80649 ymin: 41.88747 xmax: -87.59852 ymax: 41.93228
Geodetic CRS:  +proj=longlat +ellps=WGS84 +no_defs
  objectid shape_area shape_len   zip                       geometry
1       33  106052287  42720.04 60647 MULTIPOLYGON (((-87.67762 4...
2       34  127476051  48103.78 60639 MULTIPOLYGON (((-87.72683 4...
3       35   45069038  27288.61 60707 MULTIPOLYGON (((-87.785 41....
4       36   70853834  42527.99 60622 MULTIPOLYGON (((-87.66707 4...
5       37   99039621  47970.14 60651 MULTIPOLYGON (((-87.70656 4...
6       38   23506056  34689.35 60611 MULTIPOLYGON (((-87.61401 4...

Next, merge on zip code ID. The key in the Chi_Zips object is zip, whereas the key for the COVID data is ZIP.code. Always merge non-spatial to spatial data, not the other way around. Think of the spatial file as your master file that you will continue to add on to…

Chi_Zipsf <- merge(Chi_Zips, COVID.39f, by.x = "zip", by.y = "ZIP.Code", all = TRUE)

head(Chi_Zipsf)

Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.63999 ymin: 41.85317 xmax: -87.60246 ymax: 41.88913
Geodetic CRS:  +proj=longlat +ellps=WGS84 +no_defs
    zip objectid shape_area shape_len CaseRtWk39                       geometry
1 60601       27    9166246  19804.58     1451.4 MULTIPOLYGON (((-87.62271 4...
2 60602       26    4847125  14448.17     1688.1 MULTIPOLYGON (((-87.60997 4...
3 60603       19    4560229  13672.68     1107.3 MULTIPOLYGON (((-87.61633 4...
4 60604       48    4294902  12245.81     3964.2 MULTIPOLYGON (((-87.63376 4...
5 60605       20   36301276  37973.35     1420.8 MULTIPOLYGON (((-87.62064 4...
6 60606       31    6766411  12040.44     2289.6 MULTIPOLYGON (((-87.63397 4...

dim(Chi_Zips)

[1] 59  5

dim(COVID.39f)

[1] 60  2

dim(Chi_Zipsf)

[1] 60  6

2.3 Quantile Maps

Starting with a “classic epi” approach, let’s look at case rates as quantiles. We use the tmap library, and update the choropleth data classification using the style parameter. We use the Blue-Purple palette, or BuPu, from Colorbrewer.

Colorbrewer Tip: To display all Colorbrewer palette options, load the RColorBrewer library and run display.brewer.all() – or just Google “R Colorbrewer palettes.”

library(tmap)

tmap_mode("plot")

tm_shape(Chi_Zipsf) +
  tm_polygons("CaseRtWk39", 
              fill.scale  = tm_scale_intervals(style="quantile",
                                               values="BuPu"),
              fill.legend = tm_legend(title="COVID Case Rate"),
              fill.chart = tm_chart_histogram()
              )

Let’s try tertiles:

tm_shape(Chi_Zipsf) +
  tm_polygons("CaseRtWk39",
              fill.scale=tm_scale_intervals(style='quantile', 
                                            values="BuPu",
                                            n=3),
              fill.legend=tm_legend(title="COVID Case Rate"),
              fill.chart = tm_chart_histogram()
              )

2.4 Standard Deviation Maps

While quantiles are a nice start, let’s classify using a standard deviation map. Standard deviation is a statistical technique type of map based on how much the data differs from the mean.

tm_shape(Chi_Zipsf) +
  tm_polygons("CaseRtWk39", 
              fill.scale=tm_scale_intervals(style='sd',
                                            values="BuPu"),
              fill.legend=tm_legend("COVID Case Rate"),
              fill.chart = tm_chart_histogram()
              )

2.5 Jenks Maps

Another approach of data classification is natural breaks, or jenks. This approach looks for “natural breaks” in the data using a univariate clustering algorithm.

tm_shape(Chi_Zipsf) +
  tm_polygons("CaseRtWk39",
              fill.scale=tm_scale_intervals(style="jenks",
                                            values="BuPu"),
              fill.legend=tm_legend(title="COVID Case Rate"),
              fill.chart = tm_chart_histogram()
              )

The first bin doesn’t seem very intuitive. Let’s try 4 bins instead of 5 by changing the n parameter. In this version, we’ll also had a histogram and scale bar, and move the legend outside the frame to make it easier to view.

tm_shape(Chi_Zipsf) +
  tm_polygons("CaseRtWk39",
              fill.scale=tm_scale_intervals(style="jenks",
                                            values="BuPu",
                                            n=4),
              fill.legend=tm_legend(title="COVID Case Rate"),
              fill.chart = tm_chart_histogram()
              ) +
  tm_scalebar(position="left")

2.6 Integrate More Data

To explore potential disparities in COVID health outcomes, let’s bring in pre-cleaned demographic, racial, and ethnic data from the Opioid Environment Policy Scan database. This data is orginally sourced from the American Community Survey 2018 5-year estimate, which you could also pull using the tidycensus.

#CensusVar <- read.csv("data/DS01_Z.csv")
#head(CensusVar)

CensusVar <- read.csv("https://github.com/healthyregions/oeps/raw/refs/heads/main/backend/oeps/data/tables/zcta-2018.csv")
glimpse(CensusVar)

Rows: 32,989
Columns: 82
$ HEROP_ID     <chr> "860US01001", "860US01002", "860US01003", "860US01005", "…
$ ZCTA5        <dbl> 1001, 1002, 1003, 1005, 1007, 1008, 1009, 1010, 1011, 101…
$ AreaSqMi     <dbl> 11.50, 55.06, 0.71, 44.26, 52.60, 53.80, 0.81, 34.77, 31.…
$ TotPop       <int> 17621, 30066, 11238, 4991, 14967, 1197, 249, 3739, 1448, …
$ Age15_44     <dbl> 5735, 18634, 11192, 1518, 5109, 344, 119, 1042, 474, 182,…
$ Ovr65P       <dbl> 24.11, 11.87, 0.04, 15.41, 14.76, 23.98, 38.96, 18.91, 18…
$ TotPopHh     <int> 17077, 24798, 77, 4763, 14961, 1197, 249, 3739, 1437, 560…
$ NonRelFhhP   <dbl> 4.34, 2.97, 22.22, 0.71, 1.62, 1.35, 0.00, 3.25, 2.50, 5.…
$ NonRelNfhhP  <dbl> 16.15, 49.07, 45.76, 34.28, 27.19, 14.56, 62.35, 13.06, 2…
$ DisbP        <dbl> 14.3, 8.6, 2.0, 9.4, 9.5, 10.5, 0.0, 10.8, 18.1, 10.0, 14…
$ WhiteP       <dbl> 92.48, 75.77, 74.31, 95.15, 94.73, 98.83, 100.00, 95.16, …
$ BlackP       <dbl> 1.96, 6.24, 5.08, 2.30, 0.57, 0.00, 0.00, 3.61, 0.28, 0.0…
$ HispP        <dbl> 5.52, 7.34, 5.55, 1.24, 1.39, 1.67, 0.00, 2.62, 0.76, 0.5…
$ AmIndP       <dbl> 0.03, 0.41, 0.52, 0.00, 0.00, 0.33, 0.00, 0.00, 0.00, 0.0…
$ AsianP       <dbl> 3.05, 10.95, 16.22, 0.78, 3.31, 0.00, 0.00, 0.45, 0.14, 0…
$ PacIsP       <dbl> 0.00, 0.11, 0.00, 0.48, 0.00, 0.00, 0.00, 0.00, 0.00, 0.8…
$ OtherP       <dbl> 2.47, 6.52, 3.87, 1.28, 1.39, 0.84, 0.00, 0.78, 1.24, 0.0…
$ DsmBlk       <dbl> 0.16, 0.30, 0.18, 0.30, 0.27, 0.19, 0.00, 0.41, 0.16, 0.0…
$ IntrBlkWht   <dbl> 0.89, 0.69, 0.76, 0.95, 0.92, 0.96, 0.00, 0.92, 0.96, 0.0…
$ IsoBlk       <dbl> 0.02, 0.08, 0.05, 0.02, 0.01, 0.01, 0.00, 0.03, 0.01, 0.0…
$ DsmHsp       <dbl> 0.33, 0.22, 0.19, 0.13, 0.15, 0.01, 0.00, 0.19, 0.00, 0.0…
$ IntrHspWht   <dbl> 0.87, 0.71, 0.74, 0.95, 0.93, 0.96, 0.91, 0.92, 0.96, 0.9…
$ IsoHsp       <dbl> 0.08, 0.07, 0.06, 0.02, 0.02, 0.02, 0.06, 0.04, 0.02, 0.0…
$ DsmAs        <dbl> 0.30, 0.26, 0.22, 0.12, 0.07, 0.30, 0.00, 0.37, 0.39, 0.0…
$ IntrAsWht    <dbl> 0.88, 0.70, 0.74, 0.95, 0.93, 0.96, 0.91, 0.92, 0.96, 0.0…
$ IsoAs        <dbl> 0.04, 0.13, 0.14, 0.01, 0.03, 0.00, 0.03, 0.02, 0.01, 0.0…
$ TotVetPop    <dbl> 1350, 659, 0, 333, 1023, 104, 12, 232, 157, 34, 1263, 275…
$ VetP         <dbl> 9.44, 2.50, 0.00, 8.40, 8.80, 10.13, 4.82, 7.77, 13.25, 7…
$ AlcTot       <dbl> 6, 6, NA, 1, 4, NA, NA, NA, NA, NA, 5, 8, 1, NA, 6, 4, NA…
$ AlcDens      <dbl> 0.52, 0.11, NA, 0.02, 0.08, NA, NA, NA, NA, NA, 0.89, 0.6…
$ AlcPerCap    <dbl> 0, 0, NA, 0, 0, NA, NA, NA, NA, NA, 0, 0, 0, NA, 0, 0, NA…
$ TotUnits     <int> 7769, 10947, 60, 2074, 5925, 634, 206, 1710, 696, 333, 10…
$ VacantP      <dbl> 3.71, 9.77, 40.00, 13.31, 7.59, 19.56, 49.03, 16.32, 18.2…
$ MobileP      <dbl> 0.91, 0.14, 0.00, 1.30, 5.96, 0.47, 0.00, 8.95, 1.01, 0.0…
$ LngTermP     <dbl> 32.56, 24.29, 23.38, 44.15, 37.04, 40.18, 68.27, 43.81, 4…
$ RentalP      <dbl> 26.69, 53.65, 88.89, 16.24, 18.01, 4.90, 0.00, 11.53, 20.…
$ UnitDens     <dbl> 675.31, 198.80, 84.35, 46.86, 112.64, 11.78, 253.54, 49.1…
$ Ndvi         <dbl> 0.21, 0.27, 0.23, 0.31, 0.30, 0.34, 0.26, 0.31, 0.35, 0.3…
$ PovP         <dbl> 9.0, 28.7, 79.3, 6.9, 5.3, 5.0, 24.9, 2.9, 10.1, 9.3, 20.…
$ UnempP       <dbl> 3.4, 5.9, 18.4, 5.9, 4.0, 4.5, 12.4, 8.4, 4.4, 7.4, 7.3, …
$ MedInc       <dbl> 33505, 16070, 4200, 35901, 42413, 39563, 39972, 36754, 30…
$ PciE         <dbl> 35135, 29499, 4590, 35190, 42021, 36847, 32651, 38489, 30…
$ GiniCoeff    <dbl> 0.42, 0.52, 0.65, 0.41, 0.41, 0.34, 0.21, 0.43, 0.42, 0.4…
$ TotWrkE      <dbl> 8867, 15647, 4101, 2650, 8642, 655, 155, 1890, 762, 299, …
$ HghRskP      <dbl> 16.17, 5.76, 4.24, 20.08, 15.46, 26.56, 21.29, 22.49, 24.…
$ HltCrP       <dbl> 14.50, 11.18, 5.61, 12.98, 13.21, 9.31, 7.10, 17.78, 12.6…
$ RetailP      <dbl> 13.54, 11.45, 12.44, 8.60, 10.33, 11.30, 20.65, 9.84, 19.…
$ EssnWrkE     <dbl> 3906, 5110, 1675, 1242, 3443, 331, 56, 843, 421, 159, 567…
$ EssnWrkP     <dbl> 44.05, 32.66, 40.84, 46.87, 39.84, 50.53, 36.13, 44.60, 5…
$ SviTh1       <dbl> 0.37, 0.41, 0.50, 0.24, 0.21, 0.15, 0.21, 0.26, 0.15, 0.0…
$ SviTh2       <dbl> 0.32, 0.46, 0.93, 0.28, 0.18, 0.29, 0.17, 0.31, 0.28, 0.4…
$ SviTh3       <dbl> 0.50, 0.21, 0.01, 0.23, 0.22, 0.21, 0.75, 0.46, 0.22, 0.1…
$ SviTh4       <dbl> 0.29, 0.45, 0.34, 0.05, 0.29, 0.14, 0.08, 0.17, 0.14, 0.0…
$ SviSmryRnk   <dbl> 0.51, 0.57, 0.75, 0.60, 0.39, 0.24, 0.28, 0.35, 0.24, 0.0…
$ MaleP        <dbl> 47.06, 48.55, 51.33, 49.45, 50.18, 52.05, 53.82, 49.02, 4…
$ FemP         <dbl> 52.94, 51.45, 48.67, 50.55, 49.82, 47.95, 46.18, 50.98, 5…
$ Ovr16P       <dbl> 82.86, 89.19, 100.00, 83.59, 81.57, 87.05, 100.00, 82.51,…
$ Ovr18P       <dbl> 81.30, 87.64, 99.41, 79.40, 77.64, 85.80, 100.00, 79.81, …
$ Ovr21P       <dbl> 71.51, 62.36, 25.16, 69.16, 64.07, 74.02, 95.58, 67.69, 7…
$ SRatio       <dbl> 88.88, 94.36, 105.49, 97.82, 100.71, 108.54, 116.52, 96.1…
$ SRatio18     <dbl> 78.41, 95.03, 106.01, 98.15, 88.25, 105.74, 116.52, 94.29…
$ SRatio65     <dbl> 63.42, 84.11, 0.00, 130.24, 65.97, 100.70, 203.12, 85.08,…
$ WhiteE       <dbl> 16296, 22780, 8351, 4749, 14178, 1183, 249, 3558, 1424, 5…
$ AsianE       <dbl> 538, 3293, 1823, 39, 496, 0, 0, 17, 2, 0, 528, 610, 50, 5…
$ PacIsE       <dbl> 0, 34, 0, 24, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ HisE         <dbl> 973, 2206, 624, 62, 208, 20, 0, 98, 11, 3, 6376, 4924, 18…
$ AmIndE       <dbl> 6, 124, 58, 0, 0, 4, 0, 0, 0, 0, 8, 60, 0, 0, 3, 30, 0, 1…
$ TwoRaceE     <dbl> 316, 1498, 337, 62, 170, 10, 0, 29, 11, 0, 367, 731, 176,…
$ OtherE       <dbl> 119, 462, 98, 2, 38, 0, 0, 0, 7, 0, 1765, 865, 0, 0, 100,…
$ SomeCollegeP <dbl> 17.54, 8.18, 41.53, 11.63, 12.20, 14.54, 37.88, 14.76, 14…
$ EduNoHsP     <dbl> 8.52, 4.84, 22.88, 3.87, 5.27, 10.12, 5.56, 5.83, 11.29, …
$ EduHsP       <dbl> 28.25, 13.42, 11.86, 36.85, 22.32, 31.82, 6.06, 26.50, 39…
$ BachelorsP   <dbl> 19.26, 25.55, 11.02, 17.23, 24.43, 19.81, 33.84, 19.27, 1…
$ GradSclP     <dbl> 13.24, 41.46, 4.24, 15.13, 23.38, 12.33, 0.00, 16.67, 6.8…
$ ChildrenP    <dbl> 0.19, 0.12, 0.01, 0.21, 0.22, 0.14, 0.00, 0.20, 0.18, 0.1…
$ Ovr16        <dbl> 14601, 26816, 11238, 4172, 12209, 1042, 249, 3085, 1202, …
$ Ovr18        <int> 14325, 26349, 11172, 3963, 11620, 1027, 249, 2984, 1185, …
$ Ovr21        <int> 12601, 18748, 2827, 3452, 9590, 886, 238, 2531, 1037, 442…
$ Ovr65        <int> 4249, 3568, 4, 769, 2209, 287, 97, 707, 262, 135, 3407, 5…
$ Age15_44P    <dbl> 32.55, 61.98, 99.59, 30.41, 34.14, 28.74, 47.79, 27.87, 3…
$ BlackE       <int> 346, 1875, 571, 115, 85, 0, 0, 135, 4, 0, 1166, 1072, 464…
$ EngProf      <dbl> 95.65, 95.01, 96.55, 99.92, 96.64, 99.13, 95.58, 99.27, 9…

CensusVar2 <- CensusVar %>%
                  select(HEROP_ID,Ovr65P,EduNoHsP,PovP,Ovr65P,BlackP,WhiteP,HispP)
glimpse(CensusVar2)

Rows: 32,989
Columns: 7
$ HEROP_ID <chr> "860US01001", "860US01002", "860US01003", "860US01005", "860U…
$ Ovr65P   <dbl> 24.11, 11.87, 0.04, 15.41, 14.76, 23.98, 38.96, 18.91, 18.09,…
$ EduNoHsP <dbl> 8.52, 4.84, 22.88, 3.87, 5.27, 10.12, 5.56, 5.83, 11.29, 9.28…
$ PovP     <dbl> 9.0, 28.7, 79.3, 6.9, 5.3, 5.0, 24.9, 2.9, 10.1, 9.3, 20.5, 1…
$ BlackP   <dbl> 1.96, 6.24, 5.08, 2.30, 0.57, 0.00, 0.00, 3.61, 0.28, 0.00, 5…
$ WhiteP   <dbl> 92.48, 75.77, 74.31, 95.15, 94.73, 98.83, 100.00, 95.16, 98.3…
$ HispP    <dbl> 5.52, 7.34, 5.55, 1.24, 1.39, 1.67, 0.00, 2.62, 0.76, 0.54, 2…

CensusVar2$Zcta <- substr(CensusVar2$HEROP_ID, 6,12)
glimpse(CensusVar2)

Rows: 32,989
Columns: 8
$ HEROP_ID <chr> "860US01001", "860US01002", "860US01003", "860US01005", "860U…
$ Ovr65P   <dbl> 24.11, 11.87, 0.04, 15.41, 14.76, 23.98, 38.96, 18.91, 18.09,…
$ EduNoHsP <dbl> 8.52, 4.84, 22.88, 3.87, 5.27, 10.12, 5.56, 5.83, 11.29, 9.28…
$ PovP     <dbl> 9.0, 28.7, 79.3, 6.9, 5.3, 5.0, 24.9, 2.9, 10.1, 9.3, 20.5, 1…
$ BlackP   <dbl> 1.96, 6.24, 5.08, 2.30, 0.57, 0.00, 0.00, 3.61, 0.28, 0.00, 5…
$ WhiteP   <dbl> 92.48, 75.77, 74.31, 95.15, 94.73, 98.83, 100.00, 95.16, 98.3…
$ HispP    <dbl> 5.52, 7.34, 5.55, 1.24, 1.39, 1.67, 0.00, 2.62, 0.76, 0.54, 2…
$ Zcta     <chr> "01001", "01002", "01003", "01005", "01007", "01008", "01009"…

Merge to our master Zip Code dataset. There are many more zip codes than we have in our Chicago dataset, so we only keep the zip codes in the main dataset, Chi_Zipf. In this case, let’s test the merge to see if we retain 61 observations.

Chi_Zipsf.new <- merge(Chi_Zipsf, CensusVar2, by.x = "zip", by.y = "Zcta", all.x=TRUE)

dim(Chi_Zipsf.new)

[1] 60 13

glimpse(Chi_Zipsf.new)

Rows: 60
Columns: 13
$ zip        <chr> "60601", "60602", "60603", "60604", "60605", "60606", "6060…
$ objectid   <dbl> 27, 26, 19, 48, 20, 31, 29, 28, 22, 54, 38, 16, 53, 32, 8, …
$ shape_area <dbl> 9166246, 4847125, 4560229, 4294902, 36301276, 6766411, 6466…
$ shape_len  <dbl> 19804.58, 14448.17, 13672.68, 12245.81, 37973.35, 12040.44,…
$ CaseRtWk39 <dbl> 1451.4, 1688.1, 1107.3, 3964.2, 1420.8, 2289.6, 2504.1, 340…
$ HEROP_ID   <chr> "860US60601", "860US60602", "860US60603", "860US60604", "86…
$ Ovr65P     <dbl> 14.14, 0.40, 9.54, 11.89, 9.34, 13.90, 5.42, 9.91, 11.35, 1…
$ EduNoHsP   <dbl> 0.00, 0.00, 0.00, 0.00, 2.39, 0.73, 3.94, 29.85, 25.76, 3.8…
$ PovP       <dbl> 8.4, 2.4, 14.2, 22.3, 9.6, 7.6, 21.7, 22.5, 27.9, 11.7, 10.…
$ BlackP     <dbl> 5.57, 3.78, 3.24, 5.63, 17.18, 2.35, 14.51, 17.65, 24.71, 1…
$ WhiteP     <dbl> 74.17, 68.17, 63.46, 63.43, 61.20, 72.75, 58.17, 44.61, 44.…
$ HispP      <dbl> 8.68, 6.51, 9.80, 4.35, 5.84, 6.29, 8.30, 50.69, 53.44, 6.6…
$ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((-87.62271 4..., MULTIPOLYGON (…

2.7 Thematic Map Panel

To facilitate data discovery, we likely want to explore multiple maps at once. Here we’ll generate maps for multiple variables, and plot them as a map panel.

Can you think of more efficient ways to run this code? There are also other tmap tricks to optimize this further, so enjoy your journey!

tm_shape(Chi_Zipsf.new) + 
  tm_polygons("Ovr65P",
          fill.scale=tm_scale_intervals(style='jenks',
                                        values='BuPu',
                                        n=4),
          fill.chart = tm_chart_histogram()
          )

COVID <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("CaseRtWk39",
          fill.scale  = tm_scale_intervals(style='jenks', values='Reds', n=4),
          fill.legend = tm_legend("COVID Rt"),
          fill.chart = tm_chart_histogram())

Senior <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("Ovr65P", fill.scale = tm_scale_intervals(style="jenks", 
                                                   values="BuPu", n=4),
          fill.chart = tm_chart_histogram()) 

EduNoHsP <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("EduNoHsP", fill.scale = tm_scale_intervals(style='jenks', 
                                                   values='BuPu', n=4),
          fill.chart = tm_chart_histogram())

BlkP <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("BlackP", fill.scale = tm_scale_intervals(style='jenks',
                                                    values='BuPu', n=4),
          fill.chart = tm_chart_histogram()) 

Latnx <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("HispP", fill.scale = tm_scale_intervals(style='jenks', 
                                                   values='BuPu', n=4),
          fill.chart = tm_chart_histogram())

WhiP <- tm_shape(Chi_Zipsf.new) + 
  tm_polygons("WhiteP", fill.scale = tm_scale_intervals(style='jenks', 
                                                    values='BuPu', n=4),
          fill.chart = tm_chart_histogram()) 

tmap_arrange(COVID, Senior, EduNoHsP, BlkP, Latnx, WhiP)

From the results, we see that cumulative COVID outcomes for one week in September 2020 seemed to have some geographic correlation with the Latinx/Hispanic community in Chicago. At the same time, low high school diploma rates are also concentrated in these areas, and there is some intersection with other variables considered. What are additional variables you could bring in to refine your approach? Perhaps percentage of essential workers; a different age group; internet access? What about linking in health outcomes like Asthma, Hypertension, and more at a similar scale?

In modern spatial epidemiology, associations must never be taken at face value. For example, we know that it is not “race” but “racism” that drives multiple health disparities – simply looking at a specific racial/ethnic group is not enough. Thus exploring multiple variables and nurturing a curiosity to understand these complex intersections will support knowledge discovery.

2.8 Data

We’re done! Well… not so fast. Let’s save the data so we don’t have to run the codebook again to access the data. Here, we’ll save as a geojson file. This spatial format is more forgiving with long column names, which is a long-standing challenge with shapefiles. But sometimes it can be written oddly, so double-check the dimensions when reading in later.

#Chi_Zipsf_recast <- st_cast(Chi_Zipsf, "MULTIPOLYGON")

#st_write(Chi_Zipsf_recast, "data/ChiZipMaster.geojson", driver = "GeoJSON")

We could also just the data as a CSV file which may be easier for future linking. For this, we use the st_drop_geometry() function to remove the geometry data, so we’re just left with the attributes.

#write.csv(st_drop_geometry(Chi_Zipsf), "data/ChiZipMaster.csv")

Practice in Sweden

The data folder has a dataset entitled sweden_income_2022.csv. Bring it into the workspace. Then:

1. Bring it into the workspace, and inspect. Prepare for a merge.
2. Merge the income-level data to the DeSo dataset you worked with in the prior practice. Inspect.
3. Generate a series of choropleth maps of median income using different styles. Select the best, and explain why you identified it as such.

More Resources

For choropleth mapping in R:

• https://tmap.geocompx.org/charts
• https://spatialanalysis.github.io/lab_tutorials/4_R_Mapping.html