County analysis

Let’s look at how many pills have gone to each county in West Virginia.

# Uncomment and run the lines below to see if you have the packages required already installed
# packages <- c("tidyverse", "jsonlite", "knitr", "geofacet", "scales", "forcats")
# 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(arcos)
library(knitr)
library(tigris)
library(viridis)
library(dplyr)
library(ggplot2)
library(scales)
library(forcats)

First, let’s look at the pharmacies in Mingo, West Virginia.

The total pills per pharmacy in a county can be pulled in with total_pharmacies_county().

mingo <- total_pharmacies_county(county = "Mingo", state="WV", key="WaPo")

kable(head(mingo))
mingo <- readRDS("data/mingo.RDS")

kable(head(mingo))
buyer_state buyer_county buyer_dea_no buyer_name buyer_city total_dosage_unit total_records
WV MINGO BS7437064 STROSNIDER KERMIT 13168350 7691
WV MINGO BH6954401 HURLEY DRUG COMPANY INC WILLIAMSON 8890370 11138
WV MINGO FT0251227 TUG VALLEY PHARMACY, LLC WILLIAMSON 8827860 5390
WV MINGO FR0261684 RIVERSIDE PHARMACY GILBERT 1780680 2249
WV MINGO AA8151728 ADKINS PHARMACY INC GILBERT 1576200 2491
WV MINGO FG0153863 GILBERT PHARMACY GILBERT 1403720 3138

Which were the most prolific pharmacies in Mingo between 2006 and 2014?

Let’s chart them out.

ggplot(mingo,
       aes(x=total_dosage_unit, y=fct_reorder(buyer_name, total_dosage_unit))) +
  geom_segment(
       aes(x = 0,
           y=fct_reorder(buyer_name, total_dosage_unit),
           xend = total_dosage_unit,
           yend = fct_reorder(buyer_name, total_dosage_unit)),
           color = "gray50") +
           geom_point() +
    scale_x_continuous(label=comma) +
  labs(x="Dosage units", y="", 
       title = "Pills sold at Mingo, WV pharmacies",
       subtitle = "Between 2006 and 2014",
       caption = "Source: The Washington Post, ARCOS") +
  theme_minimal()

Okay, now we can look at all the counties in West Virginia.

Pull that data with summarized_county_annual().

wv <- summarized_county_annual(state="WV", key="WaPo")

kable(head(wv))
BUYER_COUNTY BUYER_STATE year count DOSAGE_UNIT countyfips
BARBOUR WV 2006 1044 508100 54001
BARBOUR WV 2007 1310 625150 54001
BARBOUR WV 2008 1632 766200 54001
BARBOUR WV 2009 1823 869860 54001
BARBOUR WV 2010 2137 945110 54001
BARBOUR WV 2011 2174 969850 54001

For easy mapping, we’ll use Census shape files pulled with the Tigris package.

## Set the option for shapefiles to load with sf
options(tigris_class = "sf")

## Function to download county shapefiles in West Virginia
wv_shape <- counties(state="WV", cb=T)

## Join the county dosage data we pulled
wv <- left_join(wv, wv_shape, by=c("countyfips"="GEOID"))
# Mapping with ggplot2, sf, and viridis

wv %>%
  ggplot(aes(geometry=geometry, fill = DOSAGE_UNIT, color = DOSAGE_UNIT)) +
  facet_wrap(~year, ncol=2) +
  geom_sf() +
  coord_sf(crs = 26915) + 
  scale_fill_viridis(direction=-1, label = comma) +
  scale_color_viridis(direction=-1, label = comma) +
  theme_void() +
  theme(panel.grid.major = element_line(colour = 'transparent')) +
  labs(title="Oxycodone and hydrocodone pills in West Virginia", caption="Source: The Washington Post, ARCOS")

Looks nice. You should probably adjust for population next.


population <- county_population(state="WV", key="WaPo") %>% 
# isolate the columns so it doesn't conflict in a join (there are doubles, that's why)
    select(countyfips, year, population)

left_join(wv, population) %>% 
  mutate(per_person=DOSAGE_UNIT/population) %>%
  ggplot(aes(geometry=geometry, fill = per_person, color = per_person)) +
  facet_wrap(~year, ncol=2) +
  geom_sf() +
  coord_sf(crs = 26915) + 
  scale_fill_viridis(direction=-1, label = comma) +
  scale_color_viridis(direction=-1, label = comma) +
  theme_void() +
  theme(panel.grid.major = element_line(colour = 'transparent')) +
  labs(title="Oxycodone and hydrocodone pills in West Virginia per person", caption="Source: The Washington Post, ARCOS")