Annual maps
In arcos: Load ARCOS Prescription Data Prepared by the Washington Post

How to make some quick maps with the ARCOS api.

First, let's load some packages.

options(rmarkdown.html_vignette.check_title = FALSE)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# Uncomment and run the lines below to see if you have the packages required already installed
# packages <- c("dplyr", "ggplot", "jsonlite", "knitr", "geofacet", "scales")
# if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
#   install.packages(setdiff(packages, rownames(installed.packages())), repos = "http://cran.us.r-project.org") # }

# These are all the packages you'll need to run everything below 

library(dplyr)
library(ggplot2)
library(arcos)
library(jsonlite)
library(knitr)
library(geofacet)
library(scales)

# Uncomment and run the lines below to see if you have the packages required already installed
# packages <- c("dplyr", "ggplot2", "jsonlite", "knitr", "geofacet", "scales")
# if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
#   install.packages(setdiff(packages, rownames(installed.packages())), repos = "http://cran.us.r-project.org") # }

# These are all the packages you'll need to run everything below 

library(dplyr)
library(ggplot2)
library(arcos)
library(jsonlite)
library(knitr)
library(geofacet)
library(scales)

states <- combined_buyer_annual(key="WaPo")

kable(head(states))

We've pulled every pharmacy across the country and how many pills they sold per year.

Let's aggregate and combine by states and year.

annual_states <- states %>% 
  group_by(BUYER_STATE, year) %>% 
  summarize(pills=sum(DOSAGE_UNIT)) %>% 
  filter(!is.na(BUYER_STATE))

kable(head(annual_states))

Looks good. Let's use the geofacet package to map out the annual pill trends.

ggplot(annual_states, aes(year, pills)) +
  geom_col() +
  facet_geo(~ BUYER_STATE, grid = "us_state_grid2") +
    scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  scale_y_continuous(label=comma) +
  labs(title = "Annual oxycodone and hydrocodone pills by state",
    caption = "Source: The Washington Post, ARCOS",
    x = "",
    y = "Dosage units") +
  theme(strip.text.x = element_text(size = 6))

annual_states %>% 
  filter(!BUYER_STATE %in% c("AE", "GU", "MP", "PR", "PW", "VI")) %>% 
ggplot(aes(year, pills)) +
  geom_col() +
  facet_geo(~ BUYER_STATE, grid = "us_state_grid2") +
    scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  scale_y_continuous(label=comma) +
  labs(title = "Annual oxycodone and hydrocodone pills by state",
    caption = "Source: The Washington Post, ARCOS",
    x = "",
    y = "Dosage units") +
  theme(strip.text.x = element_text(size = 6))

This is nice but Florida, California, and Texas stand out because of their higher population.

We should normalize this with population.

Get the annual state population with state_population().

population <- state_population(key="WaPo")

kable(head(population))

Now, we join the two data sets.

Adjust for population.

annual_states_joined <- left_join(annual_states, population) %>% 
  filter(!is.na(population))

kable(head(annual_states_joined))

Do some math...

annual_states_joined <- annual_states_joined %>% 
  mutate(pills_per=pills/population)

kable(head(annual_states_joined))

Map it out again this time per capita.

ggplot(annual_states_joined, aes(year, pills_per)) +
  geom_col() +
  facet_geo(~ BUYER_STATE, grid = "us_state_grid2") +
    scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  scale_y_continuous(label=comma) +
  labs(title = "Annual oxycodone and hydrocodone pills per person by state",
    caption = "Source: The Washington Post, ARCOS",
    x = "",
    y = "Dosage units") +
  theme(strip.text.x = element_text(size = 6))