Assignment: Shootings.

library(leaflet)
library(tidyverse)
library(lubridate)
library(htmltools)
library(rtweet)
library(readxl)
library(broom)
library(plotly)
library(DT)
shootings <- read_csv("https://docs.google.com/spreadsheets/d/1b9o6uDO18sLxBqPwl_Gh9bnhW-ev_dABH83M5Vb5L8o/export?format=csv")

For your notebook, do some further analyses on the shootings data.

  1. Create a map of the shootings.
shootings
shootings %>% 
  leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers()

Here is a map of all the shootings from the Mother Jones dataset. Most of the shootings occur in high-population areas, as there are none (fortunately) in Montana, the Dakotas, Wyoming and Idaho.

There are columns called latitude and longitude in the data that can generate dots on a map of where the shootings took place. A. Because there was a shooting in Hawaii, the map is spread out a lot by default. Center on the continental US by looking up the center of the continental US, translating it to computer, finding a good zoom level, and using setView(lng = x, lat = y, zoom = z).

shootings %>% 
  leaflet() %>% 
  addTiles() %>% 
  setView(lat = 39.833333, lng = -98.583333, zoom = 4.4) %>% 
  addCircleMarkers()

The geographic center of the lower 48 is located just south of the Kansas-Nebraska border, at Lat. 39°50’, Long. -98°35’.

B. Change the size of the dots so that the more victims, the larger the dot on the map. You can do this by setting radius = ~fatalities or radius = ~total_victims inside of addCircleMarkers(). If the dots are too big, you can do something like radius = ~total_victims/10, or radius = ~log(total_victims). I like that last one because larger numbers are reduced more than smaller ones.

shootings %>% 
  leaflet() %>% 
  addTiles() %>% 
  setView(lat = 39.833333, lng = -98.583333, zoom = 4.4) %>% 
  addCircleMarkers(radius = ~fatalities/2)

Here is the same map of the lower 48 shootings with circles indicating the number of fatalities at each shooting. The larger the circle, the more fatalities.

  1. Report some additional statistics, including:
    A. Median number of total_victims and median number of fatalities.
shootings %>% 
  summarize(count = n(), median_fatalities = median(fatalities))
  

Out of the 118 shootings, the median number of fatalities is six.

shootings %>% 
  summarize(count = n(), median_total_victims = median(total_victims))
  

Out of the 118 shootings, the median number of total victims is 10.5.

B. A histogram of the number of shootings per year. You will want to set nbinsx = 40 inside add_histogram(), because that is the approx. number of years covered by the data.

shootings %>% 
  group_by(year) %>% 
  plot_ly(x = ~count, y = ~year) %>% 
  add_histogram(nbinsx = 40)

Here is a histogram illustrating the number of shootings per year. Unfortunately the graph bars get larger the longer we measure these statistics, indicating that the prevalence of these events is increasing.

C. Heatmaps with add_histogram2dcontour() of the gender and race of the shooter, and the gender and age of the shooter. You’ll notice some problems with the data that you’ll need to fix with fct_collapse().

shootings$race %>% 
  table()
shootings %>% 
   mutate(race = as_factor(race)) %>% 
  mutate(race = fct_collapse(race, Unknown = c("-", "Other", "unclear"))) %>% 
  count(race)
shootings$gender %>% 
  table()
shootings %>% 
   mutate(gender = as_factor(gender)) %>% 
  mutate(gender = fct_collapse(gender, Unknown = c("-"))) %>% 
  count(gender)
shootings %>% 
  group_by(gender) %>% 
  plot_ly(x = ~gender, y = ~race) %>% 
  add_histogram2dcontour()

Here is a heatmap of the race and gender of the shooters. From the dataset, most of the shooters are white males as indicated by the color saturation on the graphic.

shootings %>% 
  group_by(gender) %>% 
  plot_ly(x = ~gender, y = ~age_of_shooter) %>% 
  add_histogram2dcontour()

Here is a heatmap of the gender and age of shooters. It’s heavily skewed toward male shooters, but the age range is far wider.

D. A scatterplot with plotly of the number injured by the number of fatalities in the shooting.

shootings %>%                                       
  group_by(case) %>%
  plot_ly(x = ~injured, y = ~fatalities) %>% 
  add_markers()

Here is a scatterplot of the shootings by incident (known as “case” in the dataset) with number of fatalities and injuries. It seems as though most of the shootings share similar number of fatalities and injuries, with a few exceptions.

  1. Conduct a regression analysis testing the hypothesis that the number of shootings has increased over the years.

To create the regression model, y = the number of shootings and x = year, you need to set up the data like this:

num_per_year <- shootings %>% 
  filter(fatalities > 3) %>% 
  count(year) %>% 
  filter(year < 2019)

num_per_year

In your regression model, y = n and x = year.

shootings_reg <- lm(n ~ year, data = num_per_year)
tidy(shootings_reg)
glance(shootings_reg)

The regression analysis shows a slight increase in shootings per year, as characterized by the r squared coefficient of 0.4.

num_per_year %>% 
  plot_ly(x = ~year, 
          y = ~n,
          hoverinfo = "text", text = ~paste("Shootings", n, "Year: ", year)) %>% 
  add_markers(showlegend = F) %>% 
  add_lines(y = ~fitted(shootings_reg)) %>% 
            layout(title = "Shootings Per Year",
   xaxis = list(title = "Year"),
   yaxis = list(title = "Number of Shootings"))

Here is a graph depicting the shootings per year plotted with the regression analysis. It seems as though there is a positive correlation, indicating shootings are (unfortunately) increasing over time.

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCkFzc2lnbm1lbnQ6IFNob290aW5ncy4KCgpgYGB7cn0KbGlicmFyeShsZWFmbGV0KQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShsdWJyaWRhdGUpCmxpYnJhcnkoaHRtbHRvb2xzKQpsaWJyYXJ5KHJ0d2VldCkKbGlicmFyeShyZWFkeGwpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkocGxvdGx5KQpsaWJyYXJ5KERUKQpzaG9vdGluZ3MgPC0gcmVhZF9jc3YoImh0dHBzOi8vZG9jcy5nb29nbGUuY29tL3NwcmVhZHNoZWV0cy9kLzFiOW82dURPMThzTHhCcVB3bF9HaDlibmhXLWV2X2RBQkg4M001VmI1TDhvL2V4cG9ydD9mb3JtYXQ9Y3N2IikKYGBgCgoKRm9yIHlvdXIgbm90ZWJvb2ssIGRvIHNvbWUgZnVydGhlciBhbmFseXNlcyBvbiB0aGUgc2hvb3RpbmdzIGRhdGEuCgoxLiBDcmVhdGUgYSBtYXAgb2YgdGhlIHNob290aW5ncy4KCmBgYHtyfQpzaG9vdGluZ3MKYGBgCgoKYGBge3J9CnNob290aW5ncyAlPiUgCiAgbGVhZmxldCgpICU+JSAKICBhZGRUaWxlcygpICU+JSAKICBhZGRDaXJjbGVNYXJrZXJzKCkKYGBgCkhlcmUgaXMgYSBtYXAgb2YgYWxsIHRoZSBzaG9vdGluZ3MgZnJvbSB0aGUgTW90aGVyIEpvbmVzIGRhdGFzZXQuIE1vc3Qgb2YgdGhlIHNob290aW5ncyBvY2N1ciBpbiBoaWdoLXBvcHVsYXRpb24gYXJlYXMsIGFzIHRoZXJlIGFyZSBub25lIChmb3J0dW5hdGVseSkgaW4gTW9udGFuYSwgdGhlIERha290YXMsIFd5b21pbmcgYW5kIElkYWhvLiAKCgpUaGVyZSBhcmUgY29sdW1ucyBjYWxsZWQgbGF0aXR1ZGUgYW5kIGxvbmdpdHVkZSBpbiB0aGUgZGF0YSB0aGF0IGNhbiBnZW5lcmF0ZSBkb3RzIG9uIGEgbWFwIG9mIHdoZXJlIHRoZSBzaG9vdGluZ3MgdG9vayBwbGFjZS4KQS4gQmVjYXVzZSB0aGVyZSB3YXMgYSBzaG9vdGluZyBpbiBIYXdhaWksIHRoZSBtYXAgaXMgc3ByZWFkIG91dCBhIGxvdCBieSBkZWZhdWx0LiBDZW50ZXIgb24gdGhlIGNvbnRpbmVudGFsIFVTIGJ5IGxvb2tpbmcgdXAgdGhlIGNlbnRlciBvZiB0aGUgY29udGluZW50YWwgVVMsIHRyYW5zbGF0aW5nIGl0IHRvIGNvbXB1dGVyLCBmaW5kaW5nIGEgZ29vZCB6b29tIGxldmVsLCBhbmQgdXNpbmcgc2V0VmlldyhsbmcgPSB4LCBsYXQgPSB5LCB6b29tID0geikuICAKCgpgYGB7cn0Kc2hvb3RpbmdzICU+JSAKICBsZWFmbGV0KCkgJT4lIAogIGFkZFRpbGVzKCkgJT4lIAogIHNldFZpZXcobGF0ID0gMzkuODMzMzMzLCBsbmcgPSAtOTguNTgzMzMzLCB6b29tID0gNC40KSAlPiUgCiAgYWRkQ2lyY2xlTWFya2VycygpCmBgYApUaGUgZ2VvZ3JhcGhpYyBjZW50ZXIgb2YgdGhlIGxvd2VyIDQ4IGlzIGxvY2F0ZWQganVzdCBzb3V0aCBvZiB0aGUgS2Fuc2FzLU5lYnJhc2thIGJvcmRlciwgYXQgTGF0LiAzOcKwNTAnLCBMb25nLiAtOTjCsDM1Jy4KCgpCLiBDaGFuZ2UgdGhlIHNpemUgb2YgdGhlIGRvdHMgc28gdGhhdCB0aGUgbW9yZSB2aWN0aW1zLCB0aGUgbGFyZ2VyIHRoZSBkb3Qgb24gdGhlIG1hcC4gWW91IGNhbiBkbyB0aGlzIGJ5IHNldHRpbmcgcmFkaXVzID0gfmZhdGFsaXRpZXMgb3IgcmFkaXVzID0gfnRvdGFsX3ZpY3RpbXMgaW5zaWRlIG9mIGFkZENpcmNsZU1hcmtlcnMoKS4gSWYgdGhlIGRvdHMgYXJlIHRvbyBiaWcsIHlvdSBjYW4gZG8gc29tZXRoaW5nIGxpa2UgcmFkaXVzID0gfnRvdGFsX3ZpY3RpbXMvMTAsIG9yIHJhZGl1cyA9IH5sb2codG90YWxfdmljdGltcykuIEkgbGlrZSB0aGF0IGxhc3Qgb25lIGJlY2F1c2UgbGFyZ2VyIG51bWJlcnMgYXJlIHJlZHVjZWQgbW9yZSB0aGFuIHNtYWxsZXIgb25lcy4KCgpgYGB7cn0Kc2hvb3RpbmdzICU+JSAKICBsZWFmbGV0KCkgJT4lIAogIGFkZFRpbGVzKCkgJT4lIAogIHNldFZpZXcobGF0ID0gMzkuODMzMzMzLCBsbmcgPSAtOTguNTgzMzMzLCB6b29tID0gNC40KSAlPiUgCiAgYWRkQ2lyY2xlTWFya2VycyhyYWRpdXMgPSB+ZmF0YWxpdGllcy8yKQpgYGAKSGVyZSBpcyB0aGUgc2FtZSBtYXAgb2YgdGhlIGxvd2VyIDQ4IHNob290aW5ncyB3aXRoIGNpcmNsZXMgaW5kaWNhdGluZyB0aGUgbnVtYmVyIG9mIGZhdGFsaXRpZXMgYXQgZWFjaCBzaG9vdGluZy4gVGhlIGxhcmdlciB0aGUgY2lyY2xlLCB0aGUgbW9yZSBmYXRhbGl0aWVzLiAKCgoyLiBSZXBvcnQgc29tZSBhZGRpdGlvbmFsIHN0YXRpc3RpY3MsIGluY2x1ZGluZzogIApBLiBNZWRpYW4gbnVtYmVyIG9mIHRvdGFsX3ZpY3RpbXMgYW5kIG1lZGlhbiBudW1iZXIgb2YgZmF0YWxpdGllcy4gIAoKYGBge3J9CnNob290aW5ncyAlPiUgCiAgc3VtbWFyaXplKGNvdW50ID0gbigpLCBtZWRpYW5fZmF0YWxpdGllcyA9IG1lZGlhbihmYXRhbGl0aWVzKSkKICAKYGBgCk91dCBvZiB0aGUgMTE4IHNob290aW5ncywgdGhlIG1lZGlhbiBudW1iZXIgb2YgZmF0YWxpdGllcyBpcyBzaXguCgpgYGB7cn0Kc2hvb3RpbmdzICU+JSAKICBzdW1tYXJpemUoY291bnQgPSBuKCksIG1lZGlhbl90b3RhbF92aWN0aW1zID0gbWVkaWFuKHRvdGFsX3ZpY3RpbXMpKQogIApgYGAKT3V0IG9mIHRoZSAxMTggc2hvb3RpbmdzLCB0aGUgbWVkaWFuIG51bWJlciBvZiB0b3RhbCB2aWN0aW1zIGlzIDEwLjUuIAoKCkIuIEEgaGlzdG9ncmFtIG9mIHRoZSBudW1iZXIgb2Ygc2hvb3RpbmdzIHBlciB5ZWFyLiBZb3Ugd2lsbCB3YW50IHRvIHNldCBuYmluc3ggPSA0MCBpbnNpZGUgYWRkX2hpc3RvZ3JhbSgpLCBiZWNhdXNlIHRoYXQgaXMgdGhlIGFwcHJveC4gbnVtYmVyIG9mIHllYXJzIGNvdmVyZWQgYnkgdGhlIGRhdGEuICAKCmBgYHtyfQpzaG9vdGluZ3MgJT4lIAogIGdyb3VwX2J5KHllYXIpICU+JSAKICBwbG90X2x5KHggPSB+Y291bnQsIHkgPSB+eWVhcikgJT4lIAogIGFkZF9oaXN0b2dyYW0obmJpbnN4ID0gNDApCmBgYApIZXJlIGlzIGEgaGlzdG9ncmFtIGlsbHVzdHJhdGluZyB0aGUgbnVtYmVyIG9mIHNob290aW5ncyBwZXIgeWVhci4gVW5mb3J0dW5hdGVseSB0aGUgZ3JhcGggYmFycyBnZXQgbGFyZ2VyIHRoZSBsb25nZXIgd2UgbWVhc3VyZSB0aGVzZSBzdGF0aXN0aWNzLCBpbmRpY2F0aW5nIHRoYXQgdGhlIHByZXZhbGVuY2Ugb2YgdGhlc2UgZXZlbnRzIGlzIGluY3JlYXNpbmcuIAoKCgpDLiBIZWF0bWFwcyB3aXRoIGFkZF9oaXN0b2dyYW0yZGNvbnRvdXIoKSBvZiB0aGUgZ2VuZGVyIGFuZCByYWNlIG9mIHRoZSBzaG9vdGVyLCBhbmQgdGhlIGdlbmRlciBhbmQgYWdlIG9mIHRoZSBzaG9vdGVyLiBZb3UnbGwgbm90aWNlIHNvbWUgcHJvYmxlbXMgd2l0aCB0aGUgZGF0YSB0aGF0IHlvdSdsbCBuZWVkIHRvIGZpeCB3aXRoIGZjdF9jb2xsYXBzZSgpLiAgCgpgYGB7cn0Kc2hvb3RpbmdzJHJhY2UgJT4lIAogIHRhYmxlKCkKYGBgCgpgYGB7cn0Kc2hvb3RpbmdzICU+JSAKICAgbXV0YXRlKHJhY2UgPSBhc19mYWN0b3IocmFjZSkpICU+JSAKICBtdXRhdGUocmFjZSA9IGZjdF9jb2xsYXBzZShyYWNlLCBVbmtub3duID0gYygiLSIsICJPdGhlciIsICJ1bmNsZWFyIikpKSAlPiUgCiAgY291bnQocmFjZSkKYGBgCgpgYGB7cn0Kc2hvb3RpbmdzJGdlbmRlciAlPiUgCiAgdGFibGUoKQpgYGAKCmBgYHtyfQpzaG9vdGluZ3MgJT4lIAogICBtdXRhdGUoZ2VuZGVyID0gYXNfZmFjdG9yKGdlbmRlcikpICU+JSAKICBtdXRhdGUoZ2VuZGVyID0gZmN0X2NvbGxhcHNlKGdlbmRlciwgVW5rbm93biA9IGMoIi0iKSkpICU+JSAKICBjb3VudChnZW5kZXIpCmBgYAoKCgpgYGB7cn0Kc2hvb3RpbmdzICU+JSAKICBncm91cF9ieShnZW5kZXIpICU+JSAKICBwbG90X2x5KHggPSB+Z2VuZGVyLCB5ID0gfnJhY2UpICU+JSAKICBhZGRfaGlzdG9ncmFtMmRjb250b3VyKCkKYGBgCkhlcmUgaXMgYSBoZWF0bWFwIG9mIHRoZSByYWNlIGFuZCBnZW5kZXIgb2YgdGhlIHNob290ZXJzLiBGcm9tIHRoZSBkYXRhc2V0LCBtb3N0IG9mIHRoZSBzaG9vdGVycyBhcmUgd2hpdGUgbWFsZXMgYXMgaW5kaWNhdGVkIGJ5IHRoZSBjb2xvciBzYXR1cmF0aW9uIG9uIHRoZSBncmFwaGljLgoKCmBgYHtyfQpzaG9vdGluZ3MgJT4lIAogIGdyb3VwX2J5KGdlbmRlcikgJT4lIAogIHBsb3RfbHkoeCA9IH5nZW5kZXIsIHkgPSB+YWdlX29mX3Nob290ZXIpICU+JSAKICBhZGRfaGlzdG9ncmFtMmRjb250b3VyKCkKYGBgCkhlcmUgaXMgYSBoZWF0bWFwIG9mIHRoZSBnZW5kZXIgYW5kIGFnZSBvZiBzaG9vdGVycy4gSXQncyBoZWF2aWx5IHNrZXdlZCB0b3dhcmQgbWFsZSBzaG9vdGVycywgYnV0IHRoZSBhZ2UgcmFuZ2UgaXMgZmFyIHdpZGVyLiAKCgpELiBBIHNjYXR0ZXJwbG90IHdpdGggcGxvdGx5IG9mIHRoZSBudW1iZXIgaW5qdXJlZCBieSB0aGUgbnVtYmVyIG9mIGZhdGFsaXRpZXMgaW4gdGhlIHNob290aW5nLiAKCgoKYGBge3J9CnNob290aW5ncyAlPiUgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICBncm91cF9ieShjYXNlKSAlPiUKICBwbG90X2x5KHggPSB+aW5qdXJlZCwgeSA9IH5mYXRhbGl0aWVzKSAlPiUgCiAgYWRkX21hcmtlcnMoKQpgYGAKSGVyZSBpcyBhIHNjYXR0ZXJwbG90IG9mIHRoZSBzaG9vdGluZ3MgYnkgaW5jaWRlbnQgKGtub3duIGFzICJjYXNlIiBpbiB0aGUgZGF0YXNldCkgd2l0aCBudW1iZXIgb2YgZmF0YWxpdGllcyBhbmQgaW5qdXJpZXMuIEl0IHNlZW1zIGFzIHRob3VnaCBtb3N0IG9mIHRoZSBzaG9vdGluZ3Mgc2hhcmUgc2ltaWxhciBudW1iZXIgb2YgZmF0YWxpdGllcyBhbmQgaW5qdXJpZXMsIHdpdGggYSBmZXcgZXhjZXB0aW9ucy4gCgoKMy4gQ29uZHVjdCBhIHJlZ3Jlc3Npb24gYW5hbHlzaXMgdGVzdGluZyB0aGUgaHlwb3RoZXNpcyB0aGF0IHRoZSBudW1iZXIgb2Ygc2hvb3RpbmdzIGhhcyBpbmNyZWFzZWQgb3ZlciB0aGUgeWVhcnMuCgpUbyBjcmVhdGUgdGhlIHJlZ3Jlc3Npb24gbW9kZWwsIHkgPSB0aGUgbnVtYmVyIG9mIHNob290aW5ncyBhbmQgeCA9IHllYXIsIHlvdSBuZWVkIHRvIHNldCB1cCB0aGUgZGF0YSBsaWtlIHRoaXM6CgpgYGB7cn0KbnVtX3Blcl95ZWFyIDwtIHNob290aW5ncyAlPiUgCiAgZmlsdGVyKGZhdGFsaXRpZXMgPiAzKSAlPiUgCiAgY291bnQoeWVhcikgJT4lIAogIGZpbHRlcih5ZWFyIDwgMjAxOSkKCm51bV9wZXJfeWVhcgpgYGAKCkluIHlvdXIgcmVncmVzc2lvbiBtb2RlbCwgeSA9IG4gYW5kIHggPSB5ZWFyLgoKYGBge3J9CnNob290aW5nc19yZWcgPC0gbG0obiB+IHllYXIsIGRhdGEgPSBudW1fcGVyX3llYXIpCmBgYAoKYGBge3J9CnRpZHkoc2hvb3RpbmdzX3JlZykKZ2xhbmNlKHNob290aW5nc19yZWcpCmBgYApUaGUgcmVncmVzc2lvbiBhbmFseXNpcyBzaG93cyBhIHNsaWdodCBpbmNyZWFzZSBpbiBzaG9vdGluZ3MgcGVyIHllYXIsIGFzIGNoYXJhY3Rlcml6ZWQgYnkgdGhlIHIgc3F1YXJlZCBjb2VmZmljaWVudCBvZiAwLjQuCgoKYGBge3J9Cm51bV9wZXJfeWVhciAlPiUgCiAgcGxvdF9seSh4ID0gfnllYXIsIAogICAgICAgICAgeSA9IH5uLAogICAgICAgICAgaG92ZXJpbmZvID0gInRleHQiLCB0ZXh0ID0gfnBhc3RlKCJTaG9vdGluZ3MiLCBuLCAiWWVhcjogIiwgeWVhcikpICU+JSAKICBhZGRfbWFya2VycyhzaG93bGVnZW5kID0gRikgJT4lIAogIGFkZF9saW5lcyh5ID0gfmZpdHRlZChzaG9vdGluZ3NfcmVnKSkgJT4lIAogICAgICAgICAgICBsYXlvdXQodGl0bGUgPSAiU2hvb3RpbmdzIFBlciBZZWFyIiwKICAgeGF4aXMgPSBsaXN0KHRpdGxlID0gIlllYXIiKSwKICAgeWF4aXMgPSBsaXN0KHRpdGxlID0gIk51bWJlciBvZiBTaG9vdGluZ3MiKSkKYGBgCkhlcmUgaXMgYSBncmFwaCBkZXBpY3RpbmcgdGhlIHNob290aW5ncyBwZXIgeWVhciBwbG90dGVkIHdpdGggdGhlIHJlZ3Jlc3Npb24gYW5hbHlzaXMuIEl0IHNlZW1zIGFzIHRob3VnaCB0aGVyZSBpcyBhIHBvc2l0aXZlIGNvcnJlbGF0aW9uLCBpbmRpY2F0aW5nIHNob290aW5ncyBhcmUgKHVuZm9ydHVuYXRlbHkpIGluY3JlYXNpbmcgb3ZlciB0aW1lLiAKCgoKCgo=