Last active
September 29, 2023 20:03
-
-
Save CharlesFainLehman/5776b60a0e807138e8a8ad551aba1ba7 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment