diff --git a/2B_Result_Analysis.r b/2B_Result_Analysis.r index d573ab7..3d59456 100644 --- a/2B_Result_Analysis.r +++ b/2B_Result_Analysis.r @@ -1,9 +1,9 @@ library(tidyverse) library(gt) #For nice color coded capacity limits table. +source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Functions created to make nice graphs if(!exists("SAVE_RES_LOC")){SAVE_RES_LOC <- "./Results/Primary_Simulation_Results/Main_Results/"} dir.create(SAVE_RES_LOC, recursive = TRUE, showWarnings = FALSE) - ###Process the simulations and save the main percentile results by year RES <- read_csv("Results/Simulations/Kemmerer_2024_Simulation.csv") RES[,"Year"] <- RES[,"Year"] @@ -12,33 +12,6 @@ library(gt) #For nice color coded capacity limits table. HIST <- readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940) ######Population ####Fan New - GET_DATA <- function(RES,COL_NUM){ - - YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) - FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),seq(0.01,0.99,by=0.01))})) %>% as_tibble %>% mutate(Year=YEARS) - FAN_DATA <- rbind(FAN_DATA[1,],FAN_DATA) - START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric - FAN_DATA <- FAN_DATA %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique - NUM_YEARS <- length(unique(FAN_DATA$Year) ) - FAN_DATA$Group <- rep(c(1:49,0,rev(1:49)),NUM_YEARS) - FAN_DATA <- FAN_DATA %>% group_by(Year,Group) %>% summarize(MIN=min(value),MAX=max(value)) - TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup - TEMP[,3:4] <- START_VALUE - FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble - return(FAN_DATA) - } - RES %>% pull(Sim_UUID) %>% unique %>% length() -MAKE_GRAPH <- function(GRAPH_DATA,ALPHA=0.03,COLOR='cadetblue',LINE_WIDTH=0.75){ - PLOT <- ggplot(data=GRAPH_DATA) - for(i in 1:49){ - C_DATA <- GRAPH_DATA%>% filter(Group==i) - PLOT <- PLOT +geom_ribbon(data=C_DATA,aes(x=Year,ymin=MIN,ymax=MAX),alpha=ALPHA,fill=COLOR) - -} - CI_90 <- rbind(GRAPH_DATA%>% filter(Group==20) %>% mutate(Interval='80%'),GRAPH_DATA%>% filter(Group==5) %>% mutate(Interval='95%'),GRAPH_DATA%>% filter(Group==0) %>% mutate(Interval='Median Prediction')) - PLOT <- PLOT+geom_line(aes(x=Year,y=MIN,linetype=Interval,color=Interval),linewidth=LINE_WIDTH, data=CI_90)+geom_line(aes(x=Year,y=MAX,group=Interval,linetype=Interval,color=Interval),linewidth=LINE_WIDTH ,data=CI_90)+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash")) - return(PLOT) -} POP_DATA <- GET_DATA(RES,3) POP_PLOT <- MAKE_GRAPH(POP_DATA) POP_PLOT <- POP_PLOT+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2060, by = 10),2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 500))+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) diff --git a/2C_Upper_Bound_Result_Analysis.r b/2C_Upper_Bound_Result_Analysis.r index 85026c0..2f9d240 100644 --- a/2C_Upper_Bound_Result_Analysis.r +++ b/2C_Upper_Bound_Result_Analysis.r @@ -1,5 +1,6 @@ library(tidyverse) library(gt) #For nice color coded capacity limits table. +source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Functions created to make nice graphs if(!exists("SAVE_RES_LOC")){SAVE_RES_LOC <- "./Results/Primary_Simulation_Results/Upper_Bound_Results/"} dir.create(SAVE_RES_LOC, recursive = TRUE, showWarnings = FALSE) @@ -12,56 +13,29 @@ library(gt) #For nice color coded capacity limits table. HIST <- readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940) ######Population ####Fan New - GET_DATA <- function(RES,COL_NUM){ - - YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) - FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),seq(0.01,0.99,by=0.01))})) %>% as_tibble %>% mutate(Year=YEARS) - FAN_DATA <- rbind(FAN_DATA[1,],FAN_DATA) - START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric - FAN_DATA <- FAN_DATA %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique - NUM_YEARS <- length(unique(FAN_DATA$Year) ) - FAN_DATA$Group <- rep(c(1:49,0,rev(1:49)),NUM_YEARS) - FAN_DATA <- FAN_DATA %>% group_by(Year,Group) %>% summarize(MIN=min(value),MAX=max(value)) - TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup - TEMP[,3:4] <- START_VALUE - FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble - return(FAN_DATA) - } - RES %>% pull(Sim_UUID) %>% unique %>% length() -MAKE_GRAPH <- function(GRAPH_DATA,ALPHA=0.03,COLOR='springgreen4',LINE_WIDTH=0.75){ - PLOT <- ggplot(data=GRAPH_DATA) - for(i in 1:49){ - C_DATA <- GRAPH_DATA%>% filter(Group==i) - PLOT <- PLOT +geom_ribbon(data=C_DATA,aes(x=Year,ymin=MIN,ymax=MAX),alpha=ALPHA,fill=COLOR) - -} - CI_90 <- rbind(GRAPH_DATA%>% filter(Group==20) %>% mutate(Interval='80%'),GRAPH_DATA%>% filter(Group==5) %>% mutate(Interval='95%'),GRAPH_DATA%>% filter(Group==0) %>% mutate(Interval='Median Prediction')) - PLOT <- PLOT+geom_line(aes(x=Year,y=MIN,linetype=Interval,color=Interval),linewidth=LINE_WIDTH, data=CI_90)+geom_line(aes(x=Year,y=MAX,group=Interval,linetype=Interval,color=Interval),linewidth=LINE_WIDTH ,data=CI_90)+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash")) - return(PLOT) -} POP_DATA <- GET_DATA(RES,3) -POP_PLOT <- MAKE_GRAPH(POP_DATA) +POP_PLOT <- MAKE_GRAPH(POP_DATA,COLOR='springgreen4') POP_PLOT <- POP_PLOT+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2060, by = 10),2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 500))+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Population_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) POP_PLOT dev.off() BIRTH_DATA <- GET_DATA(RES,4) -BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA) +BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA,COLOR='springgreen4') BIRTH_PLOT <- BIRTH_PLOT+geom_line(data=HIST,aes(x=Year,y=Births),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Birth Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Births")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Birth_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) BIRTH_PLOT dev.off() DEATH_DATA <- GET_DATA(RES,5) %>% filter(!is.na(MIN)) -DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA) +DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA,COLOR='springgreen4') DEATH_PLOT <- DEATH_PLOT+geom_line(data=HIST,aes(x=Year,y=Deaths),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Mortality Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Deaths")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Mortality_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) DEATH_PLOT dev.off() MIGRATION_DATA <- GET_DATA(RES,6) %>% filter(!is.na(MIN)) -MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA) +MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA,COLOR='springgreen4') MIGRATION_PLOT <- MIGRATION_PLOT+geom_line(data=HIST,aes(x=Year,y=Migration),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(-1000, 1000, by = 50))+ggtitle("Kemmerer & Diamondville, Net Migration Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Migration")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Migration_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) MIGRATION_PLOT diff --git a/2D_Lower_Bound_Result_Analysis.r b/2D_Lower_Bound_Result_Analysis.r index 3ea18cb..e9fb98d 100644 --- a/2D_Lower_Bound_Result_Analysis.r +++ b/2D_Lower_Bound_Result_Analysis.r @@ -1,5 +1,6 @@ library(tidyverse) library(gt) #For nice color coded capacity limits table. +source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Functions created to make nice graphs if(!exists("SAVE_RES_LOC")){SAVE_RES_LOC <- "./Results/Primary_Simulation_Results/Lower_Bound_Results/"} dir.create(SAVE_RES_LOC, recursive = TRUE, showWarnings = FALSE) @@ -12,56 +13,29 @@ library(gt) #For nice color coded capacity limits table. HIST <- readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940) ######Population ####Fan New - GET_DATA <- function(RES,COL_NUM){ - - YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) - FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),seq(0.01,0.99,by=0.01))})) %>% as_tibble %>% mutate(Year=YEARS) - FAN_DATA <- rbind(FAN_DATA[1,],FAN_DATA) - START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric - FAN_DATA <- FAN_DATA %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique - NUM_YEARS <- length(unique(FAN_DATA$Year) ) - FAN_DATA$Group <- rep(c(1:49,0,rev(1:49)),NUM_YEARS) - FAN_DATA <- FAN_DATA %>% group_by(Year,Group) %>% summarize(MIN=min(value),MAX=max(value)) - TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup - TEMP[,3:4] <- START_VALUE - FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble - return(FAN_DATA) - } - RES %>% pull(Sim_UUID) %>% unique %>% length() -MAKE_GRAPH <- function(GRAPH_DATA,ALPHA=0.03,COLOR='firebrick2',LINE_WIDTH=0.75){ - PLOT <- ggplot(data=GRAPH_DATA) - for(i in 1:49){ - C_DATA <- GRAPH_DATA%>% filter(Group==i) - PLOT <- PLOT +geom_ribbon(data=C_DATA,aes(x=Year,ymin=MIN,ymax=MAX),alpha=ALPHA,fill=COLOR) - -} - CI_90 <- rbind(GRAPH_DATA%>% filter(Group==20) %>% mutate(Interval='80%'),GRAPH_DATA%>% filter(Group==5) %>% mutate(Interval='95%'),GRAPH_DATA%>% filter(Group==0) %>% mutate(Interval='Median Prediction')) - PLOT <- PLOT+geom_line(aes(x=Year,y=MIN,linetype=Interval,color=Interval),linewidth=LINE_WIDTH, data=CI_90)+geom_line(aes(x=Year,y=MAX,group=Interval,linetype=Interval,color=Interval),linewidth=LINE_WIDTH ,data=CI_90)+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash")) - return(PLOT) -} POP_DATA <- GET_DATA(RES,3) -POP_PLOT <- MAKE_GRAPH(POP_DATA) +POP_PLOT <- MAKE_GRAPH(POP_DATA,COLOR='firebrick2') POP_PLOT <- POP_PLOT+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2060, by = 10),2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 500))+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Population_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) POP_PLOT dev.off() BIRTH_DATA <- GET_DATA(RES,4) -BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA) +BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA,COLOR='firebrick2') BIRTH_PLOT <- BIRTH_PLOT+geom_line(data=HIST,aes(x=Year,y=Births),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Birth Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Births")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Birth_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) BIRTH_PLOT dev.off() DEATH_DATA <- GET_DATA(RES,5) %>% filter(!is.na(MIN)) -DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA) +DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA,COLOR='firebrick2') DEATH_PLOT <- DEATH_PLOT+geom_line(data=HIST,aes(x=Year,y=Deaths),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Mortality Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Deaths")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Mortality_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) DEATH_PLOT dev.off() MIGRATION_DATA <- GET_DATA(RES,6) %>% filter(!is.na(MIN)) -MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA) +MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA,COLOR='firebrick2') MIGRATION_PLOT <- MIGRATION_PLOT+geom_line(data=HIST,aes(x=Year,y=Migration),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(-1000, 1000, by = 50))+ggtitle("Kemmerer & Diamondville, Net Migration Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Migration")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank()) png(paste0(SAVE_RES_LOC,"Migration_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600) MIGRATION_PLOT diff --git a/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r b/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r new file mode 100644 index 0000000..2fc1724 --- /dev/null +++ b/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r @@ -0,0 +1,32 @@ +#Functions to create consistent fan functions from simulations + ###Create a data set which pairs of the upper and lower confidence intervals. This allows the layers to be plotted +GET_DATA <- function(RES,COL_NUM){ + YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) + LEVELS <- seq(0.01,0.99,by=0.005) + GROUPS <- floor(length(LEVELS)/2) + FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),LEVELS)})) %>% as_tibble %>% mutate(Year=YEARS) + FAN_DATA <- rbind(FAN_DATA[1,],FAN_DATA) + START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric + FAN_DATA <- FAN_DATA %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique + NUM_YEARS <- length(unique(FAN_DATA$Year) ) + FAN_DATA$Group <- rep(c(1:GROUPS,0,rev(1:GROUPS)),NUM_YEARS) + FAN_DATA <- FAN_DATA %>% group_by(Year,Group) %>% summarize(MIN=min(value),MAX=max(value)) + TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup + TEMP[,3:4] <- START_VALUE + FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble + return(FAN_DATA) + } +###Create a data set which loops through many confidence bounds to create a layered fan chart. This output can be stacked on with other data later +MAKE_GRAPH <- function(GRAPH_DATA,ALPHA=0.03,COLOR='cadetblue',LINE_WIDTH=0.75){ + PLOT <- ggplot(data=GRAPH_DATA) + LENGTH <- floor(nrow(GRAPH_DATA)/2) + for(i in 1:LENGTH){ + C_DATA <- GRAPH_DATA%>% filter(Group==i) + PLOT <- PLOT +geom_ribbon(data=C_DATA,aes(x=Year,ymin=MIN,ymax=MAX),alpha=ALPHA,fill=COLOR) + +} + CI_90 <- rbind(GRAPH_DATA%>% filter(Group==20) %>% mutate(Interval='80%'),GRAPH_DATA%>% filter(Group==5) %>% mutate(Interval='95%'),GRAPH_DATA%>% filter(Group==0) %>% mutate(Interval='Median Prediction')) + PLOT <- PLOT+geom_line(aes(x=Year,y=MIN,linetype=Interval,color=Interval),linewidth=LINE_WIDTH, data=CI_90)+geom_line(aes(x=Year,y=MAX,group=Interval,linetype=Interval,color=Interval),linewidth=LINE_WIDTH ,data=CI_90)+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash")) + return(PLOT) +} +