From 59f3b38a50b7030606d67aaa5024e1ec0794db00 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Thu, 4 Dec 2025 13:16:29 -0700 Subject: [PATCH] Minor improvments --- 1B_Run_Full_Simulation.r | 37 ++++++++++++++++++++----------------- 2B_Result_Analysis.r | 21 +++++++++++++++------ 2 files changed, 35 insertions(+), 23 deletions(-) diff --git a/1B_Run_Full_Simulation.r b/1B_Run_Full_Simulation.r index 160642c..7e1a63d 100644 --- a/1B_Run_Full_Simulation.r +++ b/1B_Run_Full_Simulation.r @@ -15,21 +15,20 @@ source("Scripts/Load_Custom_Functions/Increment_Data_Year.r") source("Scripts/Load_Custom_Functions/Single_Age_Mortality_Trend_Simulation.r") #######Preliminary Model Inputs -YEARS_AHEAD <- 10 -NUM_SIMULATIONS <- 10^4 -ST_YEAR <- 2017 +YEARS_AHEAD <- 43 +ST_YEAR <- 2023 ################################Load Data -DEMO <- readRDS("Data/Intermediate_Inputs/Starting_Demographic_Data_Sets_of_Monte_Carlo/2016_Starting_Kemmerer_Diamondville_Demographics_Matrix.Rds") +DEMO <- readRDS("Data/Intermediate_Inputs/Starting_Demographic_Data_Sets_of_Monte_Carlo/2023_Starting_Kemmerer_Diamondville_Demographics_Matrix.Rds") -BIRTH_MOD <- readRDS("Data/Intermediate_Inputs/Birth_Regressions/Birth_Regression_2016.Rds") +BIRTH_MOD <- readRDS("Data/Intermediate_Inputs/Birth_Regressions/Birth_Regression.Rds") #Must add region as a factor with multiple levels for predict to work. Seems to check for multiple levels although that is not needed econometrics. -BIRTH_DATA <- readRDS("Data/Intermediate_Inputs/Birth_Regressions/Regression_Data/Birth_Simulation_Key_Starting_Points.Rds") %>% mutate(Region=factor(Region)) %>% filter(KEM==1,Year==2016) +BIRTH_DATA <- readRDS("Data/Intermediate_Inputs/Birth_Regressions/Regression_Data/Birth_Simulation_Key_Starting_Points.Rds") %>% mutate(Region=factor(Region)) %>% filter(KEM==1,Year==2023) -MIGRATION_ARIMA <- readRDS("Data/Intermediate_Inputs/Migration_ARIMA_Models/Kemmerer_Diamondville_Adjusted_to_Lincoln_Model_Net_Migration_ARIMA_2016.Rds") +MIGRATION_ARIMA <- readRDS("Data/Intermediate_Inputs/Migration_ARIMA_Models/Kemmerer_Diamondville_Net_Migration_ARIMA.Rds") MIGRATION_ODDS <- readRDS("Data/Intermediate_Inputs/Migration_Trends/Migration_Age_Probability_Zero_to_85.Rds") ############## #Data for death rate trends - SINGLE_AGE_MODS <- readRDS("Data/Intermediate_Inputs/Mortality_Regression_Data/Single_Sex_Age_Time_Series_Regression_2016.Rds") + SINGLE_AGE_MODS <- readRDS("Data/Intermediate_Inputs/Mortality_Regression_Data/Single_Sex_Age_Time_Series_Regression.Rds") BOUNDS <- readRDS("Data/Intermediate_Inputs/Mortality_Regression_Data/Single_Sex_Age_Bounds_for_Predictions.Rds") MAX_MALE <- BOUNDS %>% filter(Sex=='Male') %>% pull(MAX_RATE) MIN_MALE <- BOUNDS %>% filter(Sex=='Male') %>% pull(MIN_RATE) @@ -37,16 +36,15 @@ MIGRATION_ODDS <- readRDS("Data/Intermediate_Inputs/Migration_Trends/Migration_A MIN_FEMALE <- BOUNDS %>% filter(Sex=='Female') %>% pull(MIN_RATE) MIN_GAP <- BOUNDS %>% filter(Sex=='Male') %>% pull(MIN_MALE_FEMALE_GAP) MAX_GAP <- BOUNDS %>% filter(Sex=='Male') %>% pull(MAX_MALE_FEMALE_GAP) - #Adjusted for 2016 by filtering REMOVE LATER - BASELINE_AGE_ADJUST_MEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Single_Sex_Age_Population_in_2000.Rds") %>% filter(Sex=='Male',Year<=2016) %>% pull(Percent_of_Population) - BASELINE_AGE_ADJUST_WOMEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Single_Sex_Age_Population_in_2000.Rds") %>% filter(Sex=='Female',Year<=2016) %>% pull(Percent_of_Population) + BASELINE_AGE_ADJUST_MEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Single_Sex_Age_Population_in_2000.Rds") %>% filter(Sex=='Male') %>% pull(Percent_of_Population) + BASELINE_AGE_ADJUST_WOMEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Single_Sex_Age_Population_in_2000.Rds") %>% filter(Sex=='Female') %>% pull(Percent_of_Population) #Adjust to just women popualtion (Not all population percent BASELINE_AGE_ADJUST_WOMEN <- BASELINE_AGE_ADJUST_WOMEN/sum(BASELINE_AGE_ADJUST_WOMEN ) BASELINE_AGE_ADJUST_MEN <- BASELINE_AGE_ADJUST_MEN/ sum(BASELINE_AGE_ADJUST_MEN ) - MOD_MEN_ALL <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_US_Men_Mortality_by_Age_2016.Rds") - MOD_WOMEN_ALL <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_US_Women_Mortality_by_Age_2016.Rds") - MOD_LIN_MEN <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_Lincoln_Men_Mortality_by_Age_2016.Rds") - MOD_LIN_WOMEN <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_Lincoln_Women_Mortality_by_Age_2016.Rds") + MOD_MEN_ALL <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_US_Men_Mortality_by_Age.Rds") + MOD_WOMEN_ALL <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_US_Women_Mortality_by_Age.Rds") + MOD_LIN_MEN <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_Lincoln_Men_Mortality_by_Age.Rds") + MOD_LIN_WOMEN <- readRDS("Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/ARIMA_Lincoln_Women_Mortality_by_Age.Rds") XREG <- cbind(rep(0.0001,YEARS_AHEAD+1),rep(0.0001,YEARS_AHEAD+1)) #Empty data set to simulate in the future XREG <- ts(XREG,start=ST_YEAR,frequency=1) SIMULATE_MORTALITY_RATE_TRENDS <- function(){ @@ -86,7 +84,12 @@ return(list(DEMO,BIRTH_DATA,c(TOTAL_POP,TOTAL_BIRTHS,TOTAL_DEATHS,TOTAL_MIGRATIO } MIGRATION_ARIMA_MODEL <- MIGRATION_ARIMA SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL){ - MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL) )) + TERRA_POWER_EFFECT <- rep(0,YEARS_AHEAD) + POP_WORK_RATIO <-3716/1920.54 #Total population of Kemmerer in 2024 divided total employment both are found in IMPLAN region details for zip code 83101 + TERRA_POWER_EFFECT[3:7] <- POP_WORK_RATIO*310.75/5 #Total IMPLAN job estimate times adjusted for families and spread over five years + + MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL)+runif(1,-55,0))+TERRA_POWER_EFFECT) + #The runif applies a downshift ranging from the historic decline rate all the way to the Lincoln rate applied in the model FINAL_REPORT_VALUES <- matrix(NA,ncol=6,nrow=YEARS_AHEAD) colnames(FINAL_REPORT_VALUES ) <- c("Sim_UUID","Year","Population","Births","Deaths","Net_Migration") @@ -111,7 +114,7 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL TOTAL_SIMULATIONS <- 10^6 N_RUNS <-ceiling(TOTAL_SIMULATIONS/BATCH_SIZE ) - SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2016_Simulation.csv") + SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2023_Simulation.csv") NEW_RES_FILE <- !file.exists(SIM_RES_FILE) for(i in 1:N_RUNS){ # MIGRATION_MATRIX <- simulate(nsim=BATCH_SIZE,MIGRATION_ARIMA,n=YEARS_AHEAD) diff --git a/2B_Result_Analysis.r b/2B_Result_Analysis.r index e5fa984..7b7836c 100644 --- a/2B_Result_Analysis.r +++ b/2B_Result_Analysis.r @@ -2,6 +2,8 @@ library(tidyverse) ###Process the simulations and save the main percentile results by year RES <- read_csv("Results/Simulations/Kemmerer_2023_Simulation.csv") RES[,"Year"] <- RES[,"Year"]+1 +RES<- RES %>% filter(!is.na(Year)) +RES <- RES %>% filter(!is.na(Population)) HIST <- readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940) YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) ######Population @@ -48,16 +50,23 @@ library(tidyverse) ALPHA=0.2 COLOR <- 'black' #GRAPH <- - nrow(RES)/10 - ggplot(data=GRAPH_DATA)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Population,group=Percentile,color=Percentile))+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED,aes(x=Year,y=Population),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 35000, by = 250))+theme_bw()+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0) + nrow(RES)/40 +POP_PLOT <- ggplot(data=GRAPH_DATA)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Population,group=Percentile,color=Percentile))+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED,aes(x=Year,y=Population),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 35000, by = 250))+theme_bw()+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0) +POP_PLOT +ggsave("Kemmerer_Population_Fan_Chart.png",POP_PLOT ) - ggplot(data=GRAPH_DATA_MIGRATION)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Migration,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Migration),color='black',linewidth=0.75) +geom_line(data=MEDIAN_PRED_MIGRATION,aes(x=Year,y=Migration),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 35000, by = 250))+theme_bw()+ggtitle("Kemmerer & Diamondville, Migration Forecast")+ expand_limits( y = 0) +MIGRATION_PLOT <- ggplot(data=GRAPH_DATA_MIGRATION)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Migration,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Migration),color='black',linewidth=0.75) +geom_line(data=MEDIAN_PRED_MIGRATION,aes(x=Year,y=Migration),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2065, by = 5)))+ scale_y_continuous(breaks = seq(-350, 350, by = 10))+theme_bw()+ggtitle("Kemmerer & Diamondville, Migration Forecast")+ expand_limits( y = 0) +MIGRATION_PLOT +ggsave("Kemmerer_Net_Migration_Fan_Chart.png",MIGRATION_PLOT) - - ggplot(data=GRAPH_DATA_MORTALITY)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Deaths,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Deaths),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_MORTALITY,aes(x=Year,y=Deaths),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 50, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Death Forecast")+ expand_limits( y = 0) +DEATHS_PLOT <- ggplot(data=GRAPH_DATA_MORTALITY)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Deaths,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Deaths),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_MORTALITY,aes(x=Year,y=Deaths),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 50, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Death Forecast")+ expand_limits( y = 0) +DEATHS_PLOT +ggsave("Kemmerer_Mortality_Fan_Chart.png",DEATHS_PLOT ) +2025-1940 -ggplot(data=GRAPH_DATA_BIRTHS)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Births,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Births),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_BIRTHS,aes(x=Year,y=Births),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 200, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Birth Forecast")+ expand_limits( y = 0) +BIRTH_PLOT <- ggplot(data=GRAPH_DATA_BIRTHS)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Births,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Births),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_BIRTHS,aes(x=Year,y=Births),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2065, by = 5)))+ scale_y_continuous(breaks = seq(0, 200, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Birth Forecast")+ expand_limits( y = 0) +ggsave("Kemmerer_Births_Fan_Chart.png",BIRTH_PLOT )