library(tidyverse) library(fixest) #setwd("../") ##########################Model Population Trends ##Run Regression #Pull in Demographic data and create categories for key groups in the regression, male/female population with high fertility, children under one and two (but not zero). This data is broken down by each age group so aggregate to the county, year level for the final regression. #Fertility age bounds were informed by the regression found in the file ./Scripts/Other_Analysis/Select_Range_of_Male_Female_Fertility.r Which qualitatively supports that the number of people in these age ranges (18-28 Women, 18-30 Men) have the most significance in predicting birth rates. These two are combined into one variable which represent the minimum number of people in the key fertility window between the sexes, this is the binding fertility constraint and has more explanatory power than including either the number of men or women in the fertility window alone, providing a good trade off for including more variables or reducing variance. if(!exists("DEMOGRAPHIC_COUNTY_LOC")){DEMOGRAPHIC_COUNTY_LOC <- "./Data/Cleaned_Data/Demographic_Sex_Age_Data/RDS/All_Wyoming_Counties_Demographics.Rds"} if(!exists("DEMOGRAPHIC_KEM_LOC")){DEMOGRAPHIC_KEM_LOC <- "./Data/Cleaned_Data/Demographic_Sex_Age_Data/RDS/Kemmerer_Diamondville_Demographics.Rds"} if(!exists("DEMOGRAPHIC_OTHER_LIN_LOC")){DEMOGRAPHIC_OTHER_LIN_LOC <- "./Data/Cleaned_Data/Demographic_Sex_Age_Data/RDS/Other_Lincoln_Demographics.Rds"} if(!exists("POPULATION_COUNTY_LOC")){POPULATION_COUNTY_LOC <- "./Data/Cleaned_Data/Population_Data/RDS/All_Wyoming_County_Populations.Rds"} if(!exists("POPULATION_CITY_LOC")){POPULATION_CITY_LOC <- "./Data/Cleaned_Data/Population_Data/RDS/All_Wyoming_City_Populations.Rds"} #Function to make the data consistent for each data set used to run a birth simulation in the Monte Carlo #DEMO_DATA <- readRDS(DEMOGRAPHIC_KEM_LOC);POP_DATA <- readRDS(POPULATION_CITY_LOC) MAKE_REG_DATA <- function(DEMO_DATA){ return(DEMO_DATA %>% mutate(Male_Window=Age>=18 & Age<=30,Female_Window=Age>=18 & Age<=28) %>% group_by(County,Region,Year) %>% summarize(Female_Birth_Group=sum(Num_Female*Female_Window,na.rm=TRUE),Male_Birth_Group=sum(Num_Male*Male_Window,na.rm=TRUE),Min_Birth_Group=ifelse(Female_Birth_Group% ungroup) } DEMOGRAPHIC_COUNTY_DATA <- readRDS(DEMOGRAPHIC_COUNTY_LOC) COUNTY_POP <- readRDS(POPULATION_COUNTY_LOC) REG_DATA <- readRDS(POPULATION_COUNTY_LOC) %>% full_join(MAKE_REG_DATA(DEMOGRAPHIC_COUNTY_DATA)) REG_DATA <- REG_DATA %>% group_by(County,Region) %>% mutate(PREV_BIRTH=lag(Births),PREV_TWO_BIRTH=lag(Births,2)) %>% ungroup REG_DATA <- REG_DATA %>% select(-Female_Birth_Group,-Male_Birth_Group) #Store the data set of only the first year needing a birth forecast, to start the birth Monte Carlo simulations. ###Some of the years are missing births, previous births etc. Where missing fill this in by assuming all age zero children in the demographic data (DEMOGRAPHIC_LOC) were born in the last year. This makes a more complete data set. Some test find a near perfect 1 to 1 with this method #Data to fill in the missing records FILL_IN_DATA <- DEMOGRAPHIC_COUNTY_DATA %>% mutate(POP=Num_Male+Num_Female,BIRTHS=ifelse(Age==0,POP,0)) %>% group_by(County,Region,Year) %>% summarize(BIRTHS=sum(BIRTHS)) %>% arrange(County,Year) %>% mutate(ALT=lag(BIRTHS),ALT2=lag(BIRTHS,2)) %>% ungroup #Join and replace missing records REG_DATA <- REG_DATA %>% left_join(FILL_IN_DATA ) %>% mutate(Births=ifelse(is.na(Births),BIRTHS,Births),PREV_BIRTH=ifelse(is.na(PREV_BIRTH),ALT,PREV_BIRTH),PREV_TWO_BIRTH=ifelse(is.na(PREV_TWO_BIRTH),ALT2,PREV_TWO_BIRTH)) %>% select(-BIRTHS,-ALT,-ALT2) %>% select(Year,County,Region,everything()) %>% mutate(Region=County) ###Working on Kemmerer data DEMOGRAPHIC_KEM_DATA <- readRDS(DEMOGRAPHIC_KEM_LOC) readRDS(POPULATION_CITY_LOC) %>% filter(City %in% c("Kemmerer","Diamondville")) %>% group_by(Year) %>% mutate(Population=sum(Population,na.rm=TRUE)) %>% mutate(City='Kemmerer') %>% rename(Region=City) MAKE_REG_DATA(readRDS(DEMOGRAPHIC_KEM_LOC)) REG_DATA readRDS(DEMOGRAPHIC_KEM_LOC)%>% mutate(POP=Num_Male+Num_Female,Births=ifelse(Age==0,POP,0)) %>% group_by(County,Region,Year) %>% summarize(Births=sum(Births)) %>% arrange(County,Year) %>% mutate(PREV_BIRTH=lag(Births),PREV_TWO_BIRTH=lag(Births,2)) %>% ungroup readRDS(DEMOGRAPHIC_KEM_LOC)%>% mutate(Male_Window=Age>=18 & Age<=30,Female_Window=Age>=18 & Age<=28) %>% group_by(County,Year) %>% summarize(Female_Birth_Group=sum(Num_Female*Female_Window,na.rm=TRUE),Male_Birth_Group=sum(Num_Male*Male_Window,na.rm=TRUE),Min_Birth_Group=ifelse(Female_Birth_Group% ungroup DEMOGRAPHIC_DATA TEST <- readRDS(POPULATION_COUNTY_LOC) if(!("Births" %in% colnames(TEST))) "Deaths" %in% colnames(TEST) "Migration" %in% colnames(TEST) "Migration" %in% colnames(TEST) readRDS(DEMOGRAPHIC_COUNTY_LOC) readRDS(POPULATION_COUNTY_LOC) COUNTY_REG_DATA <- MAKE_REG_DATA(readRDS(DEMOGRAPHIC_COUNTY_LOC),readRDS(POPULATION_COUNTY_LOC) ) readRDS(DEMOGRAPHIC_KEM_LOC) readRDS(POPULATION_CITY_LOC) %>% readRDS(DEMOGRAPHIC_KEM_LOC) readRDS(DEMOGRAPHIC_KEM_LOC) MAKE_REG_DATA(readRDS(DEMOGRAPHIC_KEM_LOC),readRDS(POPULATION_CITY_LOC) ) %>% filter(!is.na(Region)) %>% pull(Region) %>% unique %>% pull(Region) %>% unique %>% filter(Region=='Kemmerer') readRDS(POPULATION_CITY_LOC) MAKE_REG_DATA(readRDS("Data/Cleaned_Data/Demographic_Sex_Age_Data/RDS/Kemmerer_Diamondville_Demographics.Rds"),readRDS("Data/Cleaned_Data/Population_Data/RDS/All_Wyoming_City_Populations.Rds")) #Extract the population trend data to connect with demographics (Population,births,deaths) POP_DATA <- readRDS(POPULATION_LOC) #Merger the two data sets and drop any records that cannot be used in the regression (this makes the "predict" function output the right number of records) REG_DATA <- POP_DATA %>% full_join(DEMOGRAPHIC_DATA) ###Predict the number of Births MOD_BIRTHS <- feols(log(Births)~log(PREV_BIRTH)+log(PREV_TWO_BIRTH)+log(Min_Birth_Group)+Year*County,cluster=~Year+County, data=REG_DATA ) #Higher AIC but worse acf #MOD_BIRTHS_ALT <- feols(log(Births)~log(PREV_BIRTH)+log(Min_Birth_Group)+Year*County,cluster=~Year+County, data=REG_DATA ) #AIC(MOD_BIRTHS)