R Notebook

library(tidyverse)
library(broom)
library(plotly)
library(tidycensus)      
library(sf)              
library(leaflet)         
library(trendyy)
library(usdata)          

states <- get_acs(geography = "state",               
                  variables = "B01003_001",          
                  geometry = TRUE,                   
                  shift_geo = T)                     

Getting data from the 2016-2020 5-year ACS
Warning: The `shift_geo` argument is deprecated and will be removed in a future release. We recommend using `tigris::shift_geometry()` instead.
Using feature geometry obtained from the albersusa package
Please note: Alaska and Hawaii are being shifted and are not to scale.
old-style crs object detected; please recreate object with a recent sf::st_crs()

states %>% 
  rename(location = NAME) %>% 
  inner_join(unemployment_states)

Joining, by = "location"
Simple feature collection with 51 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2100000 ymin: -2500000 xmax: 2516374 ymax: 732103.3
CRS:           +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
First 10 features:
   GEOID             location   variable estimate moe hits      keyword geo gprop                       geometry
1     04              Arizona B01003_001  7174064  NA   48 unemployment  US   web MULTIPOLYGON (((-1111066 -8...
2     05             Arkansas B01003_001  3011873  NA   35 unemployment  US   web MULTIPOLYGON (((557903.1 -1...
3     06           California B01003_001 39346023  NA   41 unemployment  US   web MULTIPOLYGON (((-1853480 -9...
4     08             Colorado B01003_001  5684926  NA   38 unemployment  US   web MULTIPOLYGON (((-613452.9 -...
5     09          Connecticut B01003_001  3570549  NA   82 unemployment  US   web MULTIPOLYGON (((2226838 519...
6     11 District of Columbia B01003_001   701974  NA   46 unemployment  US   web MULTIPOLYGON (((1960720 -41...
7     13              Georgia B01003_001 10516579  NA   41 unemployment  US   web MULTIPOLYGON (((1379893 -98...
8     17             Illinois B01003_001 12716164  NA   52 unemployment  US   web MULTIPOLYGON (((868942.5 -2...
9     18              Indiana B01003_001  6696893  NA   62 unemployment  US   web MULTIPOLYGON (((1279733 -39...
10    22            Louisiana B01003_001  4664616  NA   35 unemployment  US   web MULTIPOLYGON (((1080885 -16...

unemployment_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(fire)

Joining, by = "State"

depression_data <- depression_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(fire)

Joining, by = "State"

unemployment_data <- unemployment_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(fire)

Joining, by = "State"

states_leaflet <- get_acs(geography = "state",       
                  variables = "B19013_001",         
                  geometry = TRUE)                   

Getting data from the 2016-2020 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

fire <- read_csv("Firearm_death_rate.csv")

Rows: 150 Columns: 2
-- Column specification --------------------------------------------------------------------------------------------------------------
Delimiter: ","
chr (1): State
dbl (1): Rate

i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.

1. Here is a choropleth map of firearm death rates.

unemployment_colors <- colorNumeric(palette = "viridis", domain = unemployment_data$Rate)

states_leaflet %>% 
  rename(location = NAME) %>% 
  inner_join(unemployment_data) %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons(weight = 1,
              fillColor = ~unemployment_colors(Rate), 
              label = ~paste0(location, ", Firearm death rate = ", Rate),
              highlight = highlightOptions(weight = 2)) %>% 
  setView(-95, 40, zoom = 4) %>% 
  addLegend(pal = unemployment_colors, values = ~Rate)

Joining, by = "location"

Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'

Montana’s firearm death rate is 22.5

2. This is a choropleth of the search term ‘unemployment’.

unemployment <- trendy("unemployment", 
                   geo = "US", 
                   from = "2019-01-01", to = "2020-01-01")


unemployment_states <- unemployment %>%
  get_interest_region()

unemployment_states

NA

unemployment_colors <- colorNumeric(palette = "viridis", domain = unemployment_states$hits)

states_leaflet %>% 
  rename(location = NAME) %>% 
  inner_join(unemployment_states) %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons(weight = 1,
              fillColor = ~unemployment_colors(hits), 
              label = ~paste0(location, ", Search volume = ", hits),
              highlight = highlightOptions(weight = 2)) %>% 
  setView(-95, 40, zoom = 4) %>% 
  addLegend(pal = unemployment_colors, values = ~hits)

Joining, by = "location"

Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'

Montana’s search volume for the term ‘unemployment’ is 61.

3. Here is a linear model relating Firearm death rate data to the search term ‘unemployment’, and a scatterplot of the relationship.

unemployment_data %>% 
  drop_na() %>% 
  plot_ly(x = ~hits, 
          y = ~Rate,
          hoverinfo = "text", 
          text = ~paste("State: ", location, "<br>", "'unemployment' search rate: ", hits, "<br>", "Rate: ", Rate)) %>% 
  add_markers(showlegend = F) %>% 
  add_lines(y = ~fitted(unemployment_model)) %>% 
  layout(title = "Relationship between google searches for 'unemployment' and Firearm death rates, by state",
   xaxis = list(title = "Google search volume for 'unemployment'"),
   yaxis = list(title = "State Firearm death rate, per capita"))

NA
NA

There appears to be no correlation between the two items.

4. This is a choropleth of the search term ‘depression’.

depression <- trendy("depression", 
                   geo = "US", 
                   from = "2019-01-01", to = "2020-01-01")


depression_states <- depression %>%
  get_interest_region()

depression_states

NA

depression_colors <- colorNumeric(palette = "viridis", domain = depression_states$hits)

states_leaflet %>% 
  rename(location = NAME) %>% 
  inner_join(depression_states) %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons(weight = 1,
              fillColor = ~depression_colors(hits), 
              label = ~paste0(location, ", Search volume = ", hits),
              highlight = highlightOptions(weight = 2)) %>% 
  setView(-95, 40, zoom = 4) %>% 
  addLegend(pal = depression_colors, values = ~hits)

Joining, by = "location"

Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'

Montana’s search volume for the term ‘depression’ is 74.

5. Here is a linear model relating Firearm death rate data to the search term ‘depression’, and a scatterplot of the relationship.

depression_data %>% 
  drop_na() %>% 
  plot_ly(x = ~hits, 
          y = ~Rate,
          hoverinfo = "text", 
          text = ~paste("State: ", location, "<br>", "'depression' search rate: ", hits, "<br>", "Rate: ", Rate)) %>% 
  add_markers(showlegend = F) %>% 
  add_lines(y = ~fitted(depression_model)) %>% 
  layout(title = "Relationship between google searches for 'depression' and firearm death rates, by state",
   xaxis = list(title = "Google search volume for 'depression'"),
   yaxis = list(title = "State firearm death rate, per capita"))

NA
NA

There appears to be no correlation between the two.