library(readstata13) library(ggplot2) masterdata.smr = NULL masterdata.mortprob = NULL masterdata.mortprobdiff = NULL for(scen in c("mortval_cohort_mortadj2016","mortval_cohort_minimal_mortadj2016","mortval_cohort_nomortadj","mortval_cohort_minimal_nomortadj","mortval_cohort_mortadj2004","mortval_cohort_minimal_mortadj2004")) { cdc.smr = read.dta13(paste0("./data/",scen,"_cdc_smr.dta")) hrs.smr = read.dta13(paste0("./data/",scen,"_hrs_smr.dta")) cdc.smr$ref.pop = "CDC" cdc.smr$scenario = scen hrs.smr$ref.pop = "HRS" hrs.smr$scenario = scen cdc.smr$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(cdc.smr)) cdc.smr$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(cdc.smr)) hrs.smr$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(hrs.smr)) hrs.smr$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(hrs.smr)) masterdata.smr = rbind(masterdata.smr, cdc.smr, hrs.smr) cdc.mortprob = read.dta13(paste0("./data/",scen,"_cdc_mortbyage.dta")) cdc.mortprob$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(cdc.mortprob)) cdc.mortprob$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(cdc.mortprob)) cdc.mortprob$logpdied2yr_cdc = log(cdc.mortprob$pdied2yr_cdc,10) cdc.mortprob$logpdied2yr = log(cdc.mortprob$pdied2yr,10) # scaling excess_died by units of 1,000, but also dividing by 2.0 to annualize the number of deaths over two years cdc.mortprob$excess_died1000 = cdc.mortprob$excess_died / 2000.0 masterdata.mortprob = rbind(masterdata.mortprob, cdc.mortprob) mortprobdiff = read.dta13(paste0("./data/",scen,"_cdc_bellmiller_mortprobdiff.dta")) mortprobdiff$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(mortprobdiff)) mortprobdiff$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(mortprobdiff)) masterdata.mortprobdiff = rbind(masterdata.mortprobdiff, mortprobdiff) } masterdata.smr$agegrplab = sapply(masterdata.smr$agebin, function(x) {return(ifelse(x==85,"Age 85+",paste("Age ",x,"-",x+9,sep="")))} ) for(male in 0:1) { sexlab = ifelse(male==1,"male","female") for(model in c("Minimal Model","Full Model")) { # SMR plots p <- ggplot(masterdata.smr[masterdata.smr$male==male & masterdata.smr$Model==model,], aes(year, smr)) + geom_line(aes(color=factor(ref.pop))) + facet_grid(agegrplab ~ Mortality.Adjustment) + scale_x_continuous(breaks=seq(2006, 2014, 2.0), limits=c(2006,2014)) + geom_hline(yintercept = 1.0) + labs(color="Reference Pop.") + ylab("Standardized Mortality Ratio (observed/expected)") + xlab("Year") ggsave(filename=paste0("./plots/fem_smr_",sexlab,"_",ifelse(model=="Full Model", "full", "minimal"),".pdf"), plot=p, device="pdf", width=11, height=8.5) #print(masterdata.mortprob[masterdata.mortprob$male==male & masterdata.smr$Model==model,]) p = ggplot(masterdata.mortprob[masterdata.mortprob$male==male & masterdata.mortprob$Model==model,], aes(x=logpdied2yr_cdc, y=logpdied2yr)) + geom_abline(slope=1, alpha=0.1) + geom_abline(slope=1, intercept=log(1.1,10), alpha=0.1) + geom_abline(slope=1, intercept=log(0.9,10), alpha=0.1) + geom_text(aes(label=agegrp), size=1.25) + facet_grid(year ~ Model + Mortality.Adjustment) p = p + xlab("CDC log(q_x)") + ylab("FEM log(q_x)") + labs(title="Age-specific two-year mortality probability in base 10 logs", subtitle=paste("FEM vs. CDC:",sexlab)) # print(p) ggsave(filename=paste0("./plots/qx_compare_",sexlab,"_",ifelse(model=="Full Model", "full", "minimal"),".pdf"),plot=p, device="pdf", height=8.5, width=11) p = ggplot(masterdata.mortprob[masterdata.mortprob$male==male & masterdata.mortprob$Model==model,], aes(x=agegrp,y=excess_died1000)) + geom_point() + geom_hline(yintercept=0.0) + facet_grid(year ~ Model + Mortality.Adjustment) p = p + xlab("Age") + ylab("Number of deaths (1000s)") + labs(title=paste0("Excess number of deaths in FEM for ",sexlab,"s"), subtitle="Annualized difference in FEM simulated mortality and FEM mortality expected at CDC rates") # print(p) ggsave(filename=paste0("./plots/excess_deaths_",sexlab,"_",ifelse(model=="Full Model", "full", "minimal"),".pdf"),plot=p, device="pdf", height=8.5, width=11) p = ggplot(masterdata.mortprob[masterdata.mortprob$male==male & masterdata.mortprob$Model==model,], aes(x=agegrp,y=excess_rate)) + geom_point() + geom_hline(yintercept=1.0) + facet_grid(year ~ Model + Mortality.Adjustment) p = p + xlab("Age") + ylab("Mortality Probability Ratio") + labs(title=paste0("Excess mortality probability in FEM for ",sexlab,"s"), subtitle="Ratio of two-year mortality probability in FEM and CDC") # print(p) ggsave(filename=paste0("./plots/excess_rate_",sexlab,"_",ifelse(model=="Full Model", "full", "minimal"),".pdf"),plot=p, device="pdf", height=8.5, width=11) p = ggplot(masterdata.mortprobdiff[masterdata.mortprobdiff$male==male & masterdata.mortprobdiff$Model==model & masterdata.mortprobdiff$Mortality.Adjustment != "Mortality Adjustment 2016",], aes(x=agegrp,y=mortprobdiff)) + geom_point(aes(color=factor(series)), alpha=0.6) + coord_cartesian(ylim = c(-0.081, 0.081)) + geom_hline(yintercept=0.0) + labs(color="Forecast model") + facet_grid(year ~ Model + Mortality.Adjustment) p = p + xlab("Age") + ylab("Mortality probability difference") + labs(title=paste0("Excess mortality probability: forecasted vs. actual for ",sexlab,"s"), subtitle="Difference in two-year mortality probability between FEM, Bell & Miller (2005) and CDC") # print(p) ggsave(filename=paste0("./plots/mortprob_diff_",sexlab,"_",ifelse(model=="Full Model", "full", "minimal"),".pdf"),plot=p, device="pdf", height=8.5, width=11) } } warnings()