Ore_Wash_Opioid_Deaths.R

#You'll need to get the data for deaths and population from CDC WONDER. 
#you can get my deaths file here: https://wonder.cdc.gov/controller/saved/D176/D358F098
#and you can get my population file here: https://wonder.cdc.gov/controller/saved/D176/D358F099
#Rename them appropriately for importation

library(tidyverse)
library(tidysynth)
library(ggpubr)

#don't forget to rename this file
deaths <- read.csv("Downloads/Provisional Mortality Statistics, 2018 through Last Week(1).txt", sep = "\t") %>%
  filter(Notes == "") %>%
  select(-Notes, -Occurrence.State.Code, -Month, -Population, -Crude.Rate) %>%
  separate(Month.Code, c("Year", "Month"), sep = "/", convert = T) %>%
  #I fill supprressed cells with values between 0 and 9
  complete(Occurrence.State = unique(Occurrence.State), Year = 2018:2022, Month = 1:12, fill = list(Deaths = floor(runif(1, 0, 10)))) %>%
  filter(!Year == 2023, !(Year == 2022 & Month > 3)) %>%
  rename(State = Occurrence.State)

#don't forget to rename this file
population <- read.csv("Downloads/Provisional Mortality Statistics, 2018 through Last Week(3).txt", sep = "\t") %>%
  filter(Notes == "") %>%
  select(Residence.State, Year.Code, Population) %>%
  rename(Year = Year.Code,
         #WONDER won't give you population estimates unless you use residence state, but it shouldn't matter
         State = Residence.State)

merge(deaths, population, by = c("State", "Year")) %>%
  mutate(Rate = Deaths/Population * 100000) %>%
  group_by(State) %>%
  arrange(Year, Month) %>%
  mutate(Period = 1:n()) -> opioid.death.rates

opioid.death.rates %>%
  filter(State != "Washington") %>%
  synthetic_control(outcome = Rate,
                    unit = State,
                    time = Period,
                    i_unit = "Oregon",
                    i_time = 38,
                    generate_placebos = T) %>%
  generate_predictor(time_window = 1:37,
                     Rate = mean(Rate, na.rm = T)
                    ) %>%
  generate_weights(optimization_window = 1:37, margin_ipop = .2) %>%
  generate_control() -> oregon

opioid.death.rates %>%
  filter(State != "Oregon") %>%
  synthetic_control(outcome = Rate,
                    unit = State,
                    time = Period,
                    i_unit = "Washington",
                    i_time = 39,
                    generate_placebos = T) %>%
  generate_predictor(time_window = 1:38,
                     Rate = mean(Rate, na.rm = T)
  ) %>%
  generate_weights(optimization_window = 1:38) %>%
  generate_control() -> washington

ggarrange(plot_trends(oregon) + labs(title = "", caption = ""),
          plot_trends(washington) + labs(title = "", caption = ""),
          labels = c("Oregon", "Washington"))

ggsave("Desktop/plot1.png", width = 8, height = 5)

#diagnostics
plot_mspe_ratio(oregon)
grab_outcome(oregon) %>%
  summarise(mean(Oregon))
grab_significance(oregon)

plot_mspe_ratio(washington)
plot_placebos(washington)
grab_outcome(washington) %>%
  summarise(mean(Washington))
grab_significance(washington)