Visualizing State-Level Changes in Opioid Overdose Deaths

I’ve been working a lot lately with the CDC’s Opioid Prescribing Rate data I posted about in May. A fellow grad student at IU, Kevin Wiley, and myself, have been working on an RShiny app for visualizing maps of this data and the trends over time. We have a lot left that we want to do with it, including interactive components, allowing people to choose their own control data from CDC WONDER and the Area Health Resource File. Ultimately, we’d really like this tool to enable automated graphing of difference-in-difference trends to study pre/post trends for exogenous events (our BHAG). More posts on that to come in the next few weeks/months as it evolves.

This past week, we’ve wanted to combine this data at the state and county level with relevant data related to the opioid crisis. Naturally, the rate of overdose deaths is an important aspect to study, as it has been rising rapidly since 2013. We downloaded data from CDC WONDER using the same methodology as the CDC link above. That classifies “opioid overdose deaths” as those with T40.0-T40.4 and T40.6 as the ICD-10 codes for multiple cause of death. For reference, the screenshot below is what your multiple cause of death code search should look like in your data request. And here’s the data file from our search (state level opioid overdose death rate for 2006-2016).

Screen Shot 2018-07-28 at 10.46.41 AM

Because there are so few of these at the county level, most county data is suppressed. But the state level data is available for most states and most of the years we’re interested in. So we combined this with the opioid prescribing rate data to generate some cool animations using the gganimate package in R.

How To

I’ll walk through how to prep the data and create the gifs at the top of this post, which I think complement each other fairly well. They’re the same data, just using slightly different animation syntax. One obviously includes the state abbreviations and the other does not. The balance is that to have state abbreviations makes the trajectories shown in the second gif very hard to read, but adding abbreviations to the points makes that graphic cluttered. Since they are timed identically, I thought it would be good to just make both and present them side by side. So here’s how to make both.

Get the Data

Here’s the link again for the overdose death rate data from CDC WONDER and an .Rdata file for the opioid prescribing rates by state. To start, we’ll load up the tidyverse and gganimate and set up the data.

library(tidyverse)
library(gganimate)

Make sure the two files are in your working directory for these next data loading steps to work properly.

CDC WONDER DATA - Opioid Overdose Death Rates
od_death_s <- read_delim("od_deaths_state.txt","\t", escape_double = FALSE, trim_ws = TRUE) %>% 
  filter(is.na(State)==FALSE) %>% # remove the CDC text at the bottom
  select(c(2,4,6:8,12,16)) %>% # grab the columns we want
  `colnames<-`(c("state","year","deaths","pop","crude_rate","age_ad_rate","pcttotdeaths")) %>% # rename the columns
  mutate(crude_rate = as.numeric(crude_rate), # make sure everything is numeric for plotting; this will create NAs
         age_ad_rate = as.numeric(age_ad_rate),
         deaths = as.integer(deaths))

#### Prescribing Rate Data (also from CDC)
load("rx_rate_06_16s.Rdata")

plotdat <- left_join(opioid_rx_06_16s,od_death_s) # should merge on state and year
abreg <- data.frame(state.abb,state.region) # we need the census regions data for color-coding by region
colnames(abreg) <- c("st","reg") # I rename for merge purposes so that is easier
plotdat <- left_join(plotdat,abreg) # get the state regions for color coding

Create the GIFs

Now that we have the plotdat object, we can plot both of these. Most of the code here is identical between the two, but I’m going to post the full syntax for both just so it is discrete.

State Abbreviations

ggplot(plotdat, aes(crude_rate, rx_rate, size = deaths, color=reg, label=st)) +
  geom_text(check_overlap = TRUE, show.legend = FALSE) +
  scale_size(range = c(3,10)) +
  labs(title = "State Opioid Prescribing Rate and Overdose Death Rate, {closest_state}",
       subtitle = "Colored by region, sized by count of opioid overdose deaths",
       x = "Opioid OD Death Rate (per 100,000)",
       y = "Opioid Rx Rate (per 100)",
       caption = "OD Death Rate Data from CDC Wonder, Rx Rate Data from CDC bit.ly/cdc-opioid-data") +
  # here are the gganimate parts:
  transition_states(as.integer(year),transition_length = 2, state_length = 2) +
  enter_fade() +
  exit_shrink() +
  ease_aes('linear') +
  # now for some formatting
  theme_bw() +
  theme(plot.title = element_text(color="#666666", face="bold", size=14, hjust=0.5),
        plot.subtitle = element_text(color="#666666",face = "italic", size=12, hjust=0.5),
        axis.title = element_text(color="#666666",face = "bold", size=12, hjust=0.5),
        axis.text = element_text(color="#666666",face = "bold", size=12))

anim_save("output4.gif")

State Points with Wakes

ggplot(plotdat, aes(crude_rate, rx_rate, size = deaths, color=reg, label=st)) +
geom_point(alpha=0.8,show.legend = FALSE) +
scale_size(range = c(2,12)) +
labs(title = "State Opioid Prescribing Rate and Overdose Death Rate, {frame_time}",
subtitle = "Colored by region, sized by count of opioid overdose deaths",
x = "Opioid OD Death Rate (per 100,000)",
y = "Opioid Rx Rate (per 100)",
caption = "OD Death Rate Data from CDC Wonder, Rx Rate Data from CDC bit.ly/cdc-opioid-data") +
# here are the gganimate parts:
transition_time(as.integer(year)) + # because year is a number already, we don't have to define transition states
enter_fade() +
exit_shrink() +
shadow_wake(wake_length = 0.3) + # this creates the tails
ease_aes('linear') + # this smooths their movement through the frames of the gif
# now for some formatting
theme_bw() +
theme(plot.title = element_text(color="#666666", face="bold", size=14, hjust=0.5),
plot.subtitle = element_text(color="#666666",face = "italic", size=12, hjust=0.5),
axis.title = element_text(color="#666666",face = "bold", size=12, hjust=0.5),
axis.text = element_text(color="#666666",face = "bold", size=12))

# if you run the above code and then run the anim_save with your filename, it will save the latest rendering as a gif in your working directory
anim_save("output_wake.gif")

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s