Mapping Point Data to Polygons

Now what we’d like to do is get the shootings-per-tract data, which we can then combine with our research or clinical data to see if violence near home has any effect on our outcomes. To do this, we take the latitude and longitude of our shootings and transform them slightly so that they are understood as spatial coordinates, not just pairs of numbers. We’ll use the same map projection used in our original philadelphiaCensusTracts.

library(sp)
coordinates <- SpatialPoints(shootings[c("lng", "lat")])
proj4string(coordinates) <- proj4string(philadelphiaCensusTracts)

Let’s now apply what we know about our polygons (from philadelphiaCensusTracts) and apply that to our points. We’ll end up with a table that has one row for each shooting coordinate. Essentially, what we’re doing is taking each point, lining it up with a matching polygon, and then getting the data about that polygon, which came along with the geoJSON file we downloaded. We don’t want our NAME10 to be a factor variable, but a character field (who knows if next Census we might have something like 176.02.A, which is why we don’t use numeric).

shooting_tract_data <- over(coordinates, philadelphiaCensusTracts)
shooting_tract_data$NAME10 <- shooting_tract_data$NAME10
head(shooting_tract_data)

##   OBJECTID STATEFP10 COUNTYFP10 TRACTCE10     GEOID10 NAME10
## 1      198        42        101    003600 42101003600     36
## 2      112        42        101    029000 42101029000    290
## 3      356        42        101    028200 42101028200    282
## 4       57        42        101    010300 42101010300    103
## 5      141        42        101    017602 42101017602 176.02
## 6       98        42        101    024700 42101024700    247
##            NAMELSAD10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10  INTPTLAT10
## 1     Census Tract 36   G5020          S  964539        0 +39.9279255
## 2    Census Tract 290   G5020          S  818133    20001 +40.0296096
## 3    Census Tract 282   G5020          S  855812        0 +40.0348704
## 4    Census Tract 103   G5020          S  281357        0 +39.9678322
## 5 Census Tract 176.02   G5020          S  400458        0 +39.9967182
## 6    Census Tract 247   G5020          S  941266        0 +40.0417034
##     INTPTLON10 LOGRECNO
## 1 -075.1920206    10377
## 2 -075.1166058    10605
## 3 -075.1403327    10596
## 4 -075.2254383    10437
## 5 -075.1305528    10507
## 6 -075.1645311    10560

We see the first few lines of the Census data for each of our shootings. For example, the first shooting in our shooting data corresponds to Census tract 36, which is in State 42 (Pennsylvania) and County 101 (Philadelphia County). We can use this to find out how many shootings take place in each Census tract.

shootings_by_census_shortname <- shooting_tract_data %>% 
                                 group_by(NAME10) %>% 
                                 summarise(num_shootings = n()) %>% 
                                 ungroup() 
head(shootings_by_census_shortname)

## # A tibble: 6 x 2
##   NAME10 num_shootings
##   <fct>          <int>
## 1 1                  2
## 2 10.02              1
## 3 100                3
## 4 101                9
## 5 102                9
## 6 103               14

Handling Empty Data

Don’t forget that there are some Census tracts that aren’t represented at all in our shooting_tract_data data frame, so let’s make sure we enrich it with all the tracts that aren’t included in the shooting data. We can get those by taking the data frame of our tract data, selecting the list of all the Census tracts in Philadelphia, and making sure that if they weren’t mentioned above, we add them, but with num_shootings equal to 0.

non_shooting_tracts <- philadelphiaCensusTracts@data %>% 
                       select(NAME10) %>%
                       filter(!NAME10 %in% shootings_by_census_shortname$NAME10) %>%
                       mutate(num_shootings = 0)
head(non_shooting_tracts)

##   NAME10 num_shootings
## 1    143             0
## 2    206             0
## 3     14             0
## 4     19             0
## 5     29             0
## 6     50             0

We can now combine the tracts-with-shootings and the tracts-with-no-shootings to get an overall picture of violence by census tract:

shooting_by_tract <- rbind(shootings_by_census_shortname, non_shooting_tracts)

Aside: Data Fabrication

Let’s say we have some data about how much exercise children get per day. Note that I am going to completely fabricate this data below and fabricate a finding in order to have an interesting data product. You will have your own data to include, obviously, which isn’t made up. So if you don’t understand the next code chunk, don’t worry – I’m just making fake data.

set.seed(123)
shootings_q1 <- subset(shooting_by_tract, 
                       num_shootings <= quantile(num_shootings, 0.25))
shootings_q2 <- subset(shooting_by_tract, 
                       num_shootings > quantile(num_shootings, 0.25) & 
                         num_shootings <= quantile(num_shootings, 0.5))
shootings_q3 <- subset(shooting_by_tract, 
                       num_shootings > quantile(num_shootings, 0.5) & 
                         num_shootings <= quantile(num_shootings, 0.75))
shootings_q4 <- subset(shooting_by_tract, 
                       num_shootings > quantile(num_shootings, 0.75))

fake_mrns <- sample(c(1234567:9999999), 500, replace=FALSE)
fake_census_tract <- sample(shooting_by_tract$NAME10, 500, replace=TRUE)  
# this isn't really true, we don't have equal representation among our patients in each tract!

fake_exercise_data <- data.frame(mrn = fake_mrns, census_tract = fake_census_tract, 
                                 daily_exercise_minutes = NA_integer_, stringsAsFactors = FALSE)
fake_exercise_data$daily_exercise_minutes[which(fake_exercise_data$census_tract %in% shootings_q1$NAME10)] <- 
  sample(c(45:120), nrow(shootings_q1), replace = TRUE)
fake_exercise_data$daily_exercise_minutes[which(fake_exercise_data$census_tract %in% shootings_q2$NAME10)] <- 
  sample(c(35:90), nrow(shootings_q2), replace = TRUE)
fake_exercise_data$daily_exercise_minutes[which(fake_exercise_data$census_tract %in% shootings_q3$NAME10)] <- 
  sample(c(25:60), nrow(shootings_q3), replace = TRUE)
fake_exercise_data$daily_exercise_minutes[which(fake_exercise_data$census_tract %in% shootings_q4$NAME10)] <- 
  sample(c(15:30), nrow(shootings_q4), replace = TRUE)
fake_exercise_data$daily_exercise_minutes[which(is.na(fake_exercise_data$daily_exercise_minutes))] <- 
  sample(c(15:150), nrow(shootings_q4), replace = TRUE)
# add some noise!
fake_exercise_data$daily_exercise_minutes <- 
  round(fake_exercise_data$daily_exercise_minutes*sample(seq(.7,1.3, by=0.1), 
                                                         nrow(fake_exercise_data), replace = TRUE))

Let’s take a peek at our fake data:

head(fake_exercise_data)

##       mrn census_tract daily_exercise_minutes
## 1 3755308       274.02                     17
## 2 8144402       279.01                     25
## 3 4819425          247                     32
## 4 8974594          141                     32
## 5 9478166       279.01                     73
## 6 1633889       177.01                     31

We now have two data frames that interest us:

• One has a row per subject, and includes their Census tract and exercise amount
• One has a row per Census tract, and includes the number of shootings there

Adding Clinical/Research Data

Let’s combine them and take a quick look:

exercise_and_violence <- merge(x=fake_exercise_data, y = shooting_by_tract, 
                               by.x = "census_tract", by.y = "NAME10", all = TRUE)
head(exercise_and_violence)

##   census_tract     mrn daily_exercise_minutes num_shootings
## 1            1 8276686                     27             2
## 2        10.01      NA                     NA             0
## 3        10.02 4012634                     41             1
## 4        10.02 4819049                     86             1
## 5        10.02 6644697                     77             1
## 6          100 4318947                     40             3

Do we see any trends? Let’s do a quick plot.

ggplot(exercise_and_violence, aes(x=num_shootings, y=daily_exercise_minutes)) +
  geom_point()

Wow, it looks like there’s a trend here! (Which makes sense, since I baked it into my fake data…)

At this point, statistically, we could do a two-sample test on the top and bottom quartile of our subjects to see if near-home violence leads to decreased exercise, or create a predictive model using regression.

We could also create a map – something that is easily understandable for policymakers and philanathropists to understand, or something you could put into a publication or on your website. Since in this case you’re talking about statisics per Census tract, you’ll want to simplify your data and find something like the mean or median amount of exercise per tract. Let’s do that now:

Aggregating Clinical/Research Data

exercise_per_tract <- exercise_and_violence %>% 
                      group_by(census_tract) %>%
                      summarise(mean_exercise = mean(daily_exercise_minutes)) %>%
                      ungroup()
head(exercise_per_tract)

## # A tibble: 6 x 2
##   census_tract mean_exercise
##   <fct>                <dbl>
## 1 1                     27.0
## 2 10.01                 NA  
## 3 10.02                 68.0
## 4 100                   29.5
## 5 101                   66.0
## 6 102                   NA

And we can combine our exercise by tract with our shootings by tract, for an almost-map-ready dataset:

exercise_shootings_per_tract <- merge(x=exercise_per_tract, y=shooting_by_tract, 
                                      by.x="census_tract", by.y="NAME10",
                                      all = TRUE)
head(exercise_shootings_per_tract)

##   census_tract mean_exercise num_shootings
## 1            1          27.0             2
## 2        10.01            NA             0
## 3        10.02          68.0             1
## 4          100          29.5             3
## 5          101          66.0             9
## 6          102            NA             9

Create Static Maps

Why do I say almost-map-ready? Because we have to combine our exercise and shooting data back into the data that includes our polygons. We can again do this using merge!

final_data <- merge(x = philadelphiaCensusTracts_fortified, y = exercise_shootings_per_tract, by.x = "id", by.y="census_tract", all=TRUE)

You’ll have to be creative about how to display the data in a map, because you want to show two different things: exercise amount (which we might not have for every single tract), and number of shootings (which we do have for each tract). If you’re using static maps (non-interactive images for a poster or publication), you could consider doing side-by-side maps like this:

library(scales)
philly_violence_map <- ggplot() +  
  geom_polygon(data = final_data, aes(x=long, y=lat, group=group, fill=num_shootings),
               color= "black", size = 0.1)  + 
  coord_map() +
  scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=5))  +
  theme_nothing(legend=TRUE) +
  labs(title="Shootings in Philadelphia, 2015-2018", fill="")
philly_violence_map


philly_exercise_map <- ggplot() +  
  geom_polygon(data = final_data, aes(x=long, y=lat, group=group, fill=mean_exercise),
               color= "black", size = 0.1)  + 
  coord_map() +
  scale_fill_distiller(type="seq", trans="reverse", palette = "Blues", breaks=pretty_breaks(n=5))  +
  theme_nothing(legend=TRUE) +
  labs(title="Exercise in Average Minutes Per Day by Subjects", fill="")
philly_exercise_map

Looking at the two side by side is a bit hard, right? There’s lots of grey in our exercise map, after all… it’s hard to tell what’s going on and if any pattern exists.

What if we used fill to show exercise amount and then had the border color different for high-violence sectors? We’ll map the areas with 10 or more shootings as a separate layer on top of the rest.

high_shooting_areas <- subset(final_data, num_shootings >=10)

philly_combined_map <- ggplot() +  
  geom_polygon(data = final_data, aes(x=long, y=lat, group=group, fill=mean_exercise),
               color = "black", size=0.2 )  + 
  coord_map() +
  geom_polygon(data = high_shooting_areas, aes(x=long, y=lat, fill=mean_exercise, group=group),
               color = "red", size=0.4 )  + 
  scale_fill_distiller(type="seq", trans="reverse", palette = "Blues", breaks=pretty_breaks(n=5))  +

  theme_nothing(legend=TRUE) +
  labs(title="Exercise in Average Minutes Per Day by Subjects,\nAreas with 10+ Shootings Outlined in Red", fill="")
philly_combined_map

This way we can see the effect of shootings on exercise at a glance: the areas with more shootings have a lighter color of blue, and those areas are often lined with red. Still, we can’t see the total number of shootings – the gradations of color to do that would be hard to detect, so we opted to use a cutoff for high violence areas.

Alternatively, we could create an interactive map. This would be more appropriate for a web page or presentation, as it provides data as-needed, when you move your mouse over an area.