Individual drug mapping

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

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", "ggplot2", "lubridate", "tidyr", "scales", "leaflet")
# if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
#   install.packages(setdiff(packages, rownames(installed.packages())), repos = "http://cran.us.r-project.org") # }

library(arcos)
library(dplyr)
library(ggplot2)
library(lubridate)
library(tidyr)
library(scales)
library(leaflet)

The Washington Post's series The Opioid Files focused extensively on Oxycodone and Hydrocodone but the ARCOS database includes 12 other drugs that did not get as much attention.

This is an exploratory analysis that highlights possibilities from drilling into ARCOS data inspired by this story from The Boston Globe about Methadone Mile, an area near the South End of Boston with many recovery clinics but plagued by the effects of drug use.

First, let's get a list of the types of drugs tracked in this ARCOS data.

# what drugs are there

drugs <- drug_list(key="WaPo")

drugs

# what drugs are there
drugs <- readRDS("data/drugs_list.RDS")

drugs

We have 14 options, but let's focus on methadone, which is used to treat narcotic drug addiction.

And Boston is in Suffolk County in Massachusetts. Let's pull that data.

# okay, we want METHADONE from Suffolk County

methadone <- drug_county_biz(drug = "METHADONE", county = "Suffolk", state = "MA", key = "WaPo")

# how big is this file

nrow(methadone)

methadone <- readRDS("data/methadone.RDS")

# how big is this file

nrow(methadone)

That's a lot of orders.

Who's been ordering methadone?

# What type of buyers

methadone %>% 
  select(BUYER_DEA_NO, BUYER_BUS_ACT) %>% 
  unique() %>% 
  count(BUYER_BUS_ACT)

There are many clinics and pharmacies in Boston.

Just curious, but which type of buyer orders the most methadone?

# Clean up and consolidate the categories

methadone <- methadone %>% 
  mutate(type=case_when(
    grepl("PHARMACY", BUYER_BUS_ACT) ~ "PHARMACY",
    grepl("DETOX", BUYER_BUS_ACT) ~ "DETOX/MAINTENANCE",
    grepl("MAINT", BUYER_BUS_ACT) ~ "DETOX/MAINTENANCE",
    grepl("-VA", BUYER_BUS_ACT) ~ "VA",
    grepl("CLINIC", BUYER_BUS_ACT) ~ "CLINIC",
    TRUE ~ BUYER_BUS_ACT
  ))

# totals
methadone_summary <- methadone %>% 
  # focus only on purchase orders, which is transaction code S
  filter(TRANSACTION_CODE=="S") %>% 
  group_by(type) %>% 
  # calculate amount
  mutate(amount=CALC_BASE_WT_IN_GM*MME_Conversion_Factor) %>% 
  summarize(total=sum(amount, na.rm=T)) %>% 
  arrange(desc(total))

# chart it  
ggplot(methadone_summary, aes(x=total, y=type)) +
  geom_col() +
  scale_x_continuous(label=comma) +
  theme_minimal() +
  labs(title="Methadone dosages ordered by buyer type")

So even though there are only 14 detox and maintenance clinics in Suffolk County, they out-order everyone else when it comes to Methadone.

Map

Alright, let's map out where these orders are happening.

We need to bring in a couple new data sets from ARCOS.

Buyer details such as name and address of business using buyer_addresses() as well as the available latitude and longitude data for pharmacies using pharm_latlon().

# Let's calculate the amount of methadone purchased by every buyer in Suffolk County, MA
methadone_by_buyer <- methadone %>% 
  filter(TRANSACTION_CODE=="S") %>% 
  group_by(BUYER_DEA_NO, type) %>% 
  mutate(amount=CALC_BASE_WT_IN_GM*MME_Conversion_Factor) %>% 
  summarize(total=sum(amount, na.rm=T)) %>% 
  arrange(desc(total))

# Saving the BUYER_DEA_NO of each purchaser from our data set so far
buyer_ids <- methadone %>% 
  select(BUYER_DEA_NO) %>% 
  unique()

# Download the addresses of buyers in Suffolk County, MA
buyer_addresses <- buyer_addresses(county = "Suffolk", state = "MA", key = "WaPo")

glimpse(buyer_addresses)

# We only care about the addresses of buyers who purchased methadone
buyers <- left_join(buyer_ids, buyer_addresses)

# Let's bring over the methadone order amounts
buyers <- left_join(buyers, methadone_by_buyer)

# Download the lat and lon of pharmacies in Suffolk County, MA
buyer_latlon <- pharm_latlon(county = "Suffolk", state = "MA", key = "WaPo")

glimpse(buyer_latlon)

# Let's just focus on pharmacies and detox/maintenence facilities

buyers_pharm <- buyers %>% 
  filter(type=="PHARMACY")
buyers_detox <- buyers %>% 
  filter(type=="DETOX/MAINTENANCE")

# We can easily join buyer_latlon because pharmacies are geolocated
buyers_pharm <- left_join(buyers_pharm, buyer_latlon)

# But non-pharmacies are not geolocated.
# So I'm going to geolocate these detox/maintenance locations by hand
# Here's the result

detox <- tribble(
  ~BUYER_DEA_NO, ~lat,  ~lon,
  "PV0130562",  42.366429,  -71.058752,
  "PA0203024",  42.327501,  -71.083288,
  "RB0192574",  42.363672,  -71.05972,
  "RC0304395",  42.319285,  -71.052785,
  "RH0102549",  42.333502,  -71.066541,
  "RH0345783",  42.333502,  -71.066541,
  "RC0252558",  42.333435,  -71.073025,
  "RB0307480",  42.334507,  -71.0741,
  "RC0441751",  42.300033,  -71.101911,
  "RD0284581",  42.319239,  -71.096921,
  "RR0198336",  42.300033,  -71.101911,
  "RC0463327",  42.333138,  -71.070542
)


# okay, let's join the detox dataframe to this locations dataframe above
buyers_detox <- left_join(buyers_detox, detox)

# Now that buyers_pharm and buyers_detox both have lat and lon data, we can join them
buyers <- rbind(buyers_pharm, buyers_detox) 

buyers <- buyers %>% 
  filter(!is.na(lon))

# Just ordering the data frame for mapping purposes
buyers <- buyers %>% 
  arrange(desc(total))

# Setting up some color options
cof <- colorFactor(c("#ffa500", "#13ED3F"), domain=c("PHARMACY", "DETOX/MAINTENANCE"))

# mapping with leaflet
m <- leaflet(buyers) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  setView(-71.091112, 42.338827, zoom = 12) %>% 
  addCircleMarkers(~lon, ~lat, 
                   popup=paste0(buyers$BUYER_NAME, "<br />",
                                buyers$total, " total methadone"),
                   weight = 3,
                   radius=sqrt(buyers$total)/30, 
                   color=~cof(type),
                   stroke = FALSE, 
                   fillOpacity = 0.3)%>% 
  addLegend("bottomright", 
            colors= c("#ffa500", "#13ED3F"), 
            labels=c("Detox", "Pharmacy"), 
            title="Buyer type") 
m

Zoom into central Boston and you can see the cluster of 4-5 orange circles representing the methadone clinics.

That's Methadone Mile.

arcos documentation built on Feb. 19, 2021, 1:06 a.m.