test

I want to try and probe a question that has been raised since Las Vegas, can a state regulate firearms on it’s own, given the free trade between jurisdictions. Since there is not an open electronic federal database for firearm ownership and transactions, one ipportant source of information is the Bureau of Alcohol, Tobacco and Firearms, (ATF). They publish the trace of firearms that were recovered every year, and when possible trace the state where the firearm originated. This creates a matrix that is similar to what in economics is called Direction of Trade.

In this matrix there is the rows depict the source and the columns depict the destination. This lets one get an idea where firearms that are confiscated by the ATF orginated from. From this we can also infer which states are net importers of firearms and which states are net exporters.

I will explore this matrix in an attempt to better understand if firearms are more likely to flow between geographically adjacent states. In the end I will get to a shiny app that ties everything together, for those who want to stop here…

Here is the a short gif showing the app

Below is the script to run the app from R

pkgs <- c('reshape2','geojson','readxl','ggplot2',
'leaflet','httr','rgeolocate','shiny','sp','dplyr')

check <- sapply(pkgs,require,warn.conflicts = TRUE,character.only = TRUE)

## Loading required package: reshape2

## Loading required package: geojson

## 
## Attaching package: 'geojson'

## The following object is masked from 'package:graphics':
## 
##     polygon

## Loading required package: readxl

## Loading required package: ggplot2

## Loading required package: leaflet

## Loading required package: httr

## Loading required package: rgeolocate

## Loading required package: shiny

## Loading required package: sp

## Loading required package: dplyr

## 
## 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

if(any(!check)){
  pkgs.missing <- pkgs[!check]
  install.packages(pkgs.missing)
  check <- sapply(pkgs.missing,require,warn.conflicts = TRUE,character.only = TRUE)
  }

#shiny::runGitHub('yonicd/gunflow')

So we start by loading all the pacakges we need

pkgs <- c('reshape2','geojson','readxl','ggplot2','heatmaply',
'leaflet','httr','rgeolocate','shiny','sp','dplyr')

sapply(pkgs,require,character.only = TRUE)

## Loading required package: heatmaply

## Loading required package: plotly

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:httr':
## 
##     config

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

## Loading required package: viridis

## Loading required package: viridisLite

## 
## ---------------------
## Welcome to heatmaply version 0.10.1
## Type ?heatmaply for the main documentation.
## The github page is: https://github.com/talgalili/heatmaply/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/heatmaply/issues
## Or contact: <tal.galili@gmail.com>
## 
## 	To suppress this message use:  suppressPackageStartupMessages(library(heatmaply))
## ---------------------

##   reshape2    geojson     readxl    ggplot2  heatmaply    leaflet 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
##       httr rgeolocate      shiny         sp      dplyr 
##       TRUE       TRUE       TRUE       TRUE       TRUE

Read in the data

f0 <- tempfile(fileext = '.json')
download.file("https://raw.githubusercontent.com/rstudio/leaflet/gh-pages/json/us-states.geojson",destfile = f0)
states <- geojsonio::geojson_read(f0, what = "sp")

f1 <- tempfile(fileext = '.xlsx')
download.file('https://www.atf.gov/docs/undefined/sourcerecoverybystatecy2016xlsx/download',destfile = f1)

gun_mat <- readxl::read_xlsx(f1,col_names = TRUE,range = 'B2:BD56')%>%
data.frame()

A natural place to inspect a matrix, a heatmap. We clean up the data.frame to set it up to use the heatmaply package.

g <- gun_mat[,-1]
row.names(g) <- gun_mat[,1]

mg <- as.matrix(g)

#remove diagonal (source of firearm in the state itself)
diag(mg) <- NA

mg <- mg[!grepl('^GUAM|^US VIRGIN',rownames(mg)),!grepl('^GUAM|^US VIRGIN',colnames(mg))]

# heatmaply::heatmaply(heatmaply::percentize(mg),
#           main='Firearms Sourced and Recovered in the United States and Territories 2016',
#           xlab = 'Source',
#           ylab = 'Recovered',
#           colors = rev(heatmaply::RdYlBu(256)),
#           k_row = 5,
#           k_col = 5,
#           fontsize_row = 6,
#           fontsize_col = 6)

We can see that geographically adjacent states are placed in clusters (west coast, east coast, et. al.), and they exhibit high levels of similarity. But this still keeps the picture pretty complex.

We move on to geographic based visualizations, using the leaflet package. Using the maps we can set a state as the state of interest and check what states are recieving firearms from it or what states are supplying it with firearms. Below is an example of such a map where Ohio is the state of interest. This map shows the outflow of firearms which is transformed as a percentage of total firearms (excl Ohio itself).

#<iframe src="gunflow.html" height="600" width="100%"></iframe>

Here we see that firearms that originated in Ohio found there way to all their adjecent states (MI, IN, KY, PA, WV), but also non-adjecent ones (NY, NC, FL, CA, IL,TX).

leaflet(states) %>%
  setView(-96, 37.8, 4) %>%
  addProviderTiles("MapBox", options = providerTileOptions(
    id = "mapbox.light",
    accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN')))%>% 
  addPolygons(
  fillColor = ~pal(pct),
  weight = 2,
  smoothFactor = 0.2,
  stroke=FALSE,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 1,
    bringToFront = TRUE),
  label = labels,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto"))%>% 
  addLegend(pal = pal, values = ~pct, opacity = 0.7, title = 'Percent',
            position = "bottomleft",na.label = 'Selected State') 