R Notebook

library(tidyverse)
library(broom)
library(plotly)
library(tidycensus)      # gets census data that we can use to create maps
library(sf)              # helper package for mapping
library(leaflet)         # interactive mapping package
library(trendyy)
library(usdata)          # this package has a conversion utility for state abbreviations to full names

census_api_key("edfe7dfb752b23fea5d7a7205c084760fa0f0707")

To install your API key for use in future sessions, run this function with `install = TRUE`.

1. Pick one of these and create a choropleth map with leaflet.

states_leaflet <- get_acs(geography = "state",       # gets state by state data
                  variables = "B19013_001",          # this is state income
                  geometry = TRUE)                   # gets geometry (the maps)

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

                  # shift_geo = T                    # shifts Hawaii and Alaska

states_leaflet %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons()

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

Here is a choropleth map using leaflet.

2. Choose a google search term or phrase that is relevant to your chosen CDC data, and create a choropleth of that as well.

drinking <- read_csv("Binge_drinking_percent.csv")

── Column specification ─────────────────────────────────────────────────────────────────────────────────────
cols(
  State = col_character(),
  Percentage = col_double()
)

states_leaflet <- get_acs(geography = "state",               # gets state by state data
                  variables = "B01003_001",          # this is state population
                  geometry = TRUE,                   # gets geometry (the maps)
                  shift_geo = T)                     # shifts Hawaii and Alaska

Getting data from the 2015-2019 5-year ACS
Using feature geometry obtained from the albersusa package
Please note: Alaska and Hawaii are being shifted and are not to scale.

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


alcohol_states <- alcohol %>%
  get_interest_region()

alcohol_states

NA

states %>% 
  rename(location = NAME) %>% 
  inner_join(alcohol_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  7050299  NA   83 alcohol  US   web MULTIPOLYGON (((-1111066 -8...
2     05             Arkansas B01003_001  2999370  NA   79 alcohol  US   web MULTIPOLYGON (((557903.1 -1...
3     06           California B01003_001 39283497  NA   80 alcohol  US   web MULTIPOLYGON (((-1853480 -9...
4     08             Colorado B01003_001  5610349  NA   88 alcohol  US   web MULTIPOLYGON (((-613452.9 -...
5     09          Connecticut B01003_001  3575074  NA   84 alcohol  US   web MULTIPOLYGON (((2226838 519...
6     11 District of Columbia B01003_001   692683  NA   74 alcohol  US   web MULTIPOLYGON (((1960720 -41...
7     13              Georgia B01003_001 10403847  NA   77 alcohol  US   web MULTIPOLYGON (((1379893 -98...
8     17             Illinois B01003_001 12770631  NA   81 alcohol  US   web MULTIPOLYGON (((868942.5 -2...
9     18              Indiana B01003_001  6665703  NA   92 alcohol  US   web MULTIPOLYGON (((1279733 -39...
10    22            Louisiana B01003_001  4664362  NA   77 alcohol  US   web MULTIPOLYGON (((1080885 -16...

alcohol_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(drinking)

Joining, by = "State"

alcohol_data <- alcohol_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(drinking)

Joining, by = "State"

alcohol_colors <- colorNumeric(palette = "viridis", domain = alcohol_data$Percentage)

states_leaflet %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons(weight = 1,
              fillColor = ~state_colors(estimate), 
              label = ~paste0(NAME, ", income = ", estimate),           
              highlight = highlightOptions(weight = 2)) %>% 
  setView(-95, 40, zoom = 4) %>% 
  addLegend(pal = state_colors, values = ~estimate)

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

Here is my map using the term “alcohol”.

3. Do a linear model relating the CDC data to the good search term, and create a plotly scatterplot of the relationship.

alcohol_model <- lm(Percentage ~ hits, data = alcohol_data)

glance(alcohol_model)

tidy(alcohol_model)

alcohol_data %>% 
  drop_na() %>% 
  plot_ly(x = ~hits, 
          y = ~Percentage,
          hoverinfo = "text", 
          text = ~paste("State: ", location, "<br>", "'alcohol' search rate: ", hits, "<br>", "Binge Drinking Percentage: ", Percentage)) %>% 
  add_markers(showlegend = F) %>% 
  add_lines(y = ~fitted(alcohol_model)) %>% 
  layout(title = "Relationship between google searches for 'alcohol' and binge drinking percentage, by state",
   xaxis = list(title = "Google search volume for 'alcohol'"),
   yaxis = list(title = "State binge drinking percentage, per capita"))

NA
NA

This is a plotly scatterplot of the relationship between searches in “alcohol” and binge drinking percentage per state.

4. Choose a second google search term or phrase that is relevant to your chosen CDC data and create a choropleth map.

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


beer_states <- beer %>%
  get_interest_region()

beer_states

NA

states %>% 
  rename(location = NAME) %>% 
  inner_join(beer_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  7050299  NA   61    beer  US   web MULTIPOLYGON (((-1111066 -8...
2     05             Arkansas B01003_001  2999370  NA   48    beer  US   web MULTIPOLYGON (((557903.1 -1...
3     06           California B01003_001 39283497  NA   60    beer  US   web MULTIPOLYGON (((-1853480 -9...
4     08             Colorado B01003_001  5610349  NA   78    beer  US   web MULTIPOLYGON (((-613452.9 -...
5     09          Connecticut B01003_001  3575074  NA   60    beer  US   web MULTIPOLYGON (((2226838 519...
6     11 District of Columbia B01003_001   692683  NA   63    beer  US   web MULTIPOLYGON (((1960720 -41...
7     13              Georgia B01003_001 10403847  NA   50    beer  US   web MULTIPOLYGON (((1379893 -98...
8     17             Illinois B01003_001 12770631  NA   65    beer  US   web MULTIPOLYGON (((868942.5 -2...
9     18              Indiana B01003_001  6665703  NA   56    beer  US   web MULTIPOLYGON (((1279733 -39...
10    22            Louisiana B01003_001  4664362  NA   48    beer  US   web MULTIPOLYGON (((1080885 -16...

beer_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(drinking)

Joining, by = "State"

beer_data <- beer_states %>% 
  mutate(State = state2abbr(location)) %>% 
  inner_join(drinking)

Joining, by = "State"

 beer_colors <- colorNumeric(palette = "viridis", domain = beer_data$hits)

states_leaflet %>% 
  leaflet() %>% 
  addTiles() %>%
  addPolygons(weight = 1,
              fillColor = ~beer_colors(estimate), 
              label = ~paste0(NAME, ", income = ", estimate),           
              highlight = highlightOptions(weight = 2)) %>% 
  setView(-95, 40, zoom = 4) %>% 
  addLegend(pal = beer_colors, values = ~estimate)

sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'Some values were outside the color scale and will be treated as NASome values were outside the color scale and will be treated as NA

This is a choropleth map now using the word beer.

5. Do the linear model and plotly graph between your second term and the CDC data.

beer_model <- lm(Percentage ~ hits, data = beer_data)

glance(beer_model)

tidy(beer_model)

beer_data %>% 
  drop_na() %>% 
  plot_ly(x = ~hits, 
          y = ~Percentage,
          hoverinfo = "text", 
          text = ~paste("State: ", location, "<br>", "'beer' search rate: ", hits, "<br>", "Binge Drinking Percentage: ", Percentage)) %>% 
  add_markers(showlegend = F) %>% 
  add_lines(y = ~fitted(naloxone_model)) %>% 
  layout(title = "Relationship between google searches for 'beer' and binge drinking percentage, by state",
   xaxis = list(title = "Google search volume for 'beer'"),
   yaxis = list(title = "State binge drinking percentage, per capita"))

NA
NA

This is a linear model and plotly graph between the term “beer” and the CDC data.