From daf0587f7b5d60f532178aaefc334f5c9538d7d6 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Thu, 20 Nov 2025 17:25:35 -0700 Subject: [PATCH] Working on ARIMA of Wyoming mortality. --- Mortality_Rate_Analysis.r | 168 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 160 insertions(+), 8 deletions(-) diff --git a/Mortality_Rate_Analysis.r b/Mortality_Rate_Analysis.r index 7901640..9cef05e 100644 --- a/Mortality_Rate_Analysis.r +++ b/Mortality_Rate_Analysis.r @@ -1,4 +1,11 @@ library(tidyverse) +LIN_1979 <- read_csv("Data/Raw_Data/Mortality_Rates_New/Lincoln_Age_Adjusted_1979-1998.csv") %>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`)%>% mutate(Region='Lincoln') +WY_1979 <- read_csv("Data/Raw_Data/Mortality_Rates_New/Wyoming_Age_Adjusted_1979-1998.csv") %>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`)%>% mutate(Region='Wyoming') +US_1979 <- read_csv("Data/Raw_Data/Mortality_Rates_New/US_Age_Adjusted_1979-1998.csv")%>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`) %>% filter(!is.na(Sex),!is.na(Year)) %>% mutate(Region="US")%>% filter(Year<2018) + + + + LIN_1999 <- read_csv("Data/Raw_Data/Mortality_Rates_New/Lincoln_Age_Adjusted_1999-2020.csv") %>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`)%>% mutate(Region='Lincoln') WY_1999 <- read_csv("Data/Raw_Data/Mortality_Rates_New/Wyoming_Age_Adjusted_1999-2020.csv") %>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`)%>% mutate(Region='Wyoming') %>% filter(Year<2018) US_1999 <- read_csv("Data/Raw_Data/Mortality_Rates_New/US_Age_Adjusted_1999-2020.csv")%>% select(Year,Sex,Mort_Rate=`Age Adjusted Rate`) %>% filter(!is.na(Sex),!is.na(Year)) %>% mutate(Region="US")%>% filter(Year<2018) @@ -11,22 +18,167 @@ LIN_2018<- read_csv("Data/Raw_Data/Mortality_Rates_New/Lincoln_Not_Age_Adjusted_ ADJUST_TERM <- LIN_2018 %>% rename(UNADJUSTED=Mort_Rate) %>% inner_join(LIN_1999) %>% filter(!is.na(Year)) %>% mutate(Ratio=Mort_Rate/UNADJUSTED) %>% group_by(Sex) %>% summarize(Ratio=mean(Ratio)) LIN_2018 <- LIN_2018 %>% filter(Year>2020) %>% left_join(ADJUST_TERM) %>% mutate(Mort_Rate=Mort_Rate*Ratio) %>% select(-Ratio) -DF <- rbind(LIN_1999,LIN_2018,WY_1999,US_1999,US_2018,WY_2018) %>% filter(!is.na(Year),!is.na(Sex)) +DF <- rbind(LIN_1999,LIN_2018,WY_1999,US_1999,US_2018,WY_2018,WY_1979,US_1979,LIN_1979) %>% filter(!is.na(Year),!is.na(Sex)) ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+geom_smooth(method="lm")+ facet_grid(. ~ Sex) -ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+geom_smooth()+ facet_grid(. ~ Sex) -ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+ facet_grid(. ~ Sex) +ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+geom_smooth(span=0.4)+ facet_grid(. ~ Sex) +ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+geom_line()+ facet_grid(. ~ Sex) ggplot(DF,aes(x=Year,y=Mort_Rate,group=Region,color=Region,fill=Region))+geom_point()+geom_smooth() +########################################################ARIMA +PANDIMIC_INDEX <- read_csv("https://fred.stlouisfed.org/graph/fredgraph.csv?bgcolor=%23ebf3fb&chart_type=line&drp=0&fo=open%20sans&graph_bgcolor=%23ffffff&height=450&mode=fred&recession_bars=off&txtcolor=%23444444&ts=12&tts=12&width=1320&nt=0&thu=0&trc=0&show_legend=yes&show_axis_titles=yes&show_tooltip=yes&id=WUPI&scale=left&cosd=1996-01-01&coed=2025-07-01&line_color=%230073e6&link_values=false&line_style=solid&mark_type=none&mw=3&lw=3&ost=-99999&oet=99999&mma=0&fml=a&fq=Annual&fam=sum&fgst=lin&fgsnd=2020-02-01&line_index=1&transformation=lin&vintage_date=2025-11-20&revision_date=2025-11-20&nd=1996-01-01") %>% mutate(Year=year(observation_date)) %>% select(Year,WUPI)%>% mutate(L_WUPI=lag(WUPI),L_TWO_WUPI=lag(WUPI,2)) +DF <- DF %>% ungroup%>% group_by(Sex,Region) %>% arrange(Sex,Region,Year) %>% mutate(L_Mort_Rate=lag(Mort_Rate)) %>% ungroup +REG_DATA <- DF %>% pivot_wider(values_from=c(Mort_Rate,L_Mort_Rate),names_from=Region) +REG_DATA <- REG_DATA %>% left_join(PANDIMIC_INDEX) %>% mutate(WUPI=ifelse(is.na(WUPI),0,WUPI),L_WUPI=ifelse(is.na(L_WUPI),0,L_WUPI)) - -REG_DATA <- DF %>% pivot_wider(values_from=Mort_Rate,names_from=Region) library(fixest) +MOD_WOMEN <- feols(Mort_Rate_US~L_Mort_Rate_US+Year+WUPI,REG_DATA %>% filter(Sex=='Female')) +acf(resid(MOD_WOMEN)) +MOD_WOMEN +MOD_MEN <- feols(Mort_Rate_US~L_Mort_Rate_US+Year+WUPI,REG_DATA %>% filter(Sex=='Male')) +###Lincoln REG_DATA -feols(log(Lincoln)~log(US)+Sex+Year,REG_DATA) -feols((Lincoln)~Sex+Year,REG_DATA) +MOD_WOMEN <- feols(Mort_Rate_Lincoln~Mort_Rate_US,REG_DATA %>% filter(Sex=='Female')) +acf(resid(MOD_WOMEN)) +MOD_WOMEN +MOD_MEN <- feols(Mort_Rate_US~L_Mort_Rate_US+Year+WUPI,REG_DATA %>% filter(Sex=='Male')) + +acf(resid(MOD_MEN)) +plot(resid(MOD_MEN)) +plot(predict(MOD_MEN)) +library(forecast) +DATA_MEN <- REG_DATA %>% filter(Sex=='Male') +DATA_WOMEN <- REG_DATA %>% filter(Sex=='Female') + +ST_YEAR <- DATA_MEN %>% pull(Year) %>% min + +TS_MEN_US <- DATA_MEN %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1) +TS_MEN_LIN <- DATA_MEN %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1) +TS_PANDEMIC <- DATA_MEN %>% select(WUPI,L_WUPI) %>% ts(start=ST_YEAR,frequency=1) + +TS_WOMEN_US <- DATA_WOMEN %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1) +TS_WOMEN_LIN <- DATA_WOMEN %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1) + +TS_MEN_US_INV <- TS_MEN_US + +TS_WOMEN_US_INV <- TS_WOMEN_US +FORECAST_XREG <- TS_PANDEMIC +FORECAST_XREG[,] <- 0 +MOD_US_MEN <- auto.arima(TS_MEN_US,lambda=0,biasadj=TRUE,xreg=TS_PANDEMIC) +checkresiduals(MOD_US_MEN) + +plot(forecast(MOD_US_MEN,xreg=FORECAST_XREG)) +MOD_US_WOMEN <- auto.arima(TS_WOMEN_US,lambda=0,biasadj=TRUE,xreg=TS_PANDEMIC) +checkresiduals(MOD_US_WOMEN) +plot(forecast(MOD_US_WOMEN,xreg=FORECAST_XREG)) + +MOD_LIN <- auto.arima(TS_MEN_LIN,biasadj=TRUE,xreg=TS_MEN_US) +simulate(MOD_US_MEN,xreg=FORECAST_XREG) +plot(simulate(MOD_LIN,xreg=simulate(MOD_US_MEN,xreg=FORECAST_XREG))) +plot(forecast(MOD_LIN,xreg=simulate(MOD_US_MEN,xreg=FORECAST_XREG))) +################################Other work +SINGLE_DATA <- read_csv("Data/Raw_Data/Mortality_Rates_New/US_Single_Age_1999-2020.csv") %>% select(Year,Sex,Age=`Single-Year Ages Code`,Mortality_Rate=`Crude Rate`) %>% mutate(Mortality_Rate=parse_number(Mortality_Rate)) %>% filter(!is.na(Mortality_Rate)) + +OLDER <- read_csv("Data/Raw_Data/Mortality_Rates_New/US_10_Year_Age_Groups_1999-2020.csv")%>% rename(Age=`Ten-Year Age Groups Code`,Mortality_Rate=`Crude Rate`) %>% filter(Age=='85+')%>% mutate(Age=85,Year=as.numeric(Year),Mortality_Rate=parse_number(Mortality_Rate)) %>% select(Year,Sex,Age,Mortality_Rate)%>% select(colnames(SINGLE_DATA)) +SINGLE_DATA <- rbind(SINGLE_DATA,OLDER) %>% left_join(REG_DATA %>% select(Year,Sex,US_Rate=Mort_Rate_US)) +SINGLE_DATA_PLAIN <- SINGLE_DATA + +SINGLE_DATA <- SINGLE_DATA %>% group_by(Age,Sex) %>% mutate(INDEX=Mortality_Rate/sum(ifelse(Year==min(Year),Mortality_Rate,0)),US_INDEX=US_Rate/sum(ifelse(Year==min(Year),US_Rate,0))) %>% ungroup +ggplot(SINGLE_DATA,aes(x=Year,y=INDEX,group=Age,color=Age))+geom_point()+geom_smooth(aes(y=US_INDEX),se=FALSE,color="black",linetype=2,linewidth=2) + facet_grid(rows= vars(Sex)) + +ggplot(SINGLE_DATA,aes(x=Year,y=INDEX,group=Age,color=Age))+geom_point(size=0.5)+geom_line(aes(y=US_INDEX),color="black",linetype=2,linewidth=1) + facet_grid(rows= vars(Sex)) + +ggplot(SINGLE_DATA,aes(x=Year,y=Mortality_Rate,group=Age,color=Age))+geom_point() + scale_y_log10() + facet_grid(. ~ Sex) + +ggplot(SINGLE_DATA,aes(x=Year,y=INDEX,group=Age,color=Age))+geom_smooth(se=FALSE,method="lm")+ scale_y_log10() + facet_grid(. ~ Sex) +COR_DAT <- SINGLE_DATA %>% arrange(Sex,Age) %>% select(-US_Rate,-US_INDEX,-INDEX) %>% pivot_wider(values_from=c(Mortality_Rate),names_from=c(Sex,Age)) %>% arrange(Year) %>% select(-Year) %>% as.matrix +rownames(COR_DAT) <- 1999:2020 + +TEST <- matrix(as.numeric(COR_DAT[1,]),nrow(COR_DAT),nrow=nrow(COR_DAT),ncol=ncol(COR_DAT)) +COR_DAT <- COR_DAT/TEST +#COR_DAT <- t(COR_DAT) +#US_COR_DATA <- SINGLE_DATA %>% select(Year,Sex,US_Rate) %>% unique %>% pivot_wider(values_from=US_Rate,names_from=Sex) %>% arrange(Year) +#COR_DAT <- cbind(US_COR_DATA,COR_DAT) %>% as.matrix +library(factoextra) +COR_DAT +corrplot(cor(COR_DAT[,1:86]),type="lower") +corrplot(cor(COR_DAT[,87:172]),type="lower") + +COR_DAT <- t(COR_DAT) +COR_DAT <- COR_DAT[,-1] + +fviz_nbclust(COR_DAT,kmeans,"gap_stat") +fviz_nbclust(COR_DAT,kmeans,"wss",nboot=1000) +fviz_nbclust(COR_DAT,kmeans,nboot=1000) +km.res <- kmeans(COR_DAT, 3, nstart = 1) +fviz_cluster(km.res,COR_DAT) +print(km.res) + +COR_DAT +SINGLE_DATA %>% filter(!is.numeric(Mortality_Rate)) +SINGLE_DATA <- SINGLE_DATA %>% mutate(Age=as.numeric(Age),US_Rate=as.numeric(US_Rate),Year=as.numeric(Year),Mortality_Rate=as.numeric(Mortality_Rate)) +SINGLE_DATA %>% filter(is.numeric(Age),is.numeric(Mortality_Rate),is.numeric(US_Rate)) +DATA <- SINGLE_DATA %>% left_join(PANDIMIC_INDEX ) %>% ungroup %>% group_by(Sex,Year) %>% mutate(Rank=rank(Mortality_Rate)) %>% ungroup %>% group_by(Year,Sex) %>% mutate(PER_MORT_MAX=Mortality_Rate/max(Mortality_Rate)) +feols(Mortality_Rate~factor(Rank)+Sex+Year,DATA) +MOD <- feols(log(Mortality_Rate)~log(US_Rate)*Sex+factor(Rank)+Sex|Year,DATA) +MOD +plot(resid(MOD)) +DATA$RESID <- resid(MOD) +acf((DATA %>% filter(Age==85,Sex=='Female') %>% pull(RESID))) + DATA +ggplot(SINGLE_DATA,aes(x=Year,y=Mortality_Rate,group=Age,color=Age))+geom_point() + scale_y_log10() + facet_grid(. ~ Sex) + +TEST <- SINGLE_DATA %>% group_by(Year,Sex) %>% mutate(TEST=rank(Mortality_Rate)) %>%ungroup +ggplot(TEST,aes(x=Year,y=TEST,group=Age,color=Age)) +geom_point()+ geom_smooth() +TEMP <- SINGLE_DATA %>% group_by(Sex,Year) %>% summarize(GAP=max(Mortality_Rate)-min(Mortality_Rate)) +ggplot(TEMP,aes(x=Year,y=GAP,group=Sex,color=Sex))+geom_line() + +MOD <- feols(log(Mortality_Rate)~factor(Age)*(log(US_Rate)+WUPI+L_WUPI)+Year,data=as_tibble(SINGLE_DATA)) +SINGLE_DATA +plot(resid(MOD)) +TEMP <- SINGLE_DATA +TEMP$RESID <- resid(MOD) +acf(TEMP %>% filter(Age==80) %>% pull(RESID)) +etable(MOD,group=list("Single Age"="factor")) +SINGLE_DATA[87,] +resid(MOD) +predict + + +REG_SINGLE_DATA <- SINGLE_DATA_PLAIN %>% pivot_wider(values_from="Mortality_Rate",names_from=c("Age")) +MALE <- REG_SINGLE_DATA %>% filter(Sex=='Male') +FEMALE <- REG_SINGLE_DATA %>% filter(Sex=='Female') +library(corrplot) +corrplot(cor(MALE %>% select(-Sex))) +MALE +corrplot(cor(cbind(MALE[,1],MALE[,4:ncol(MALE)]/t(MALE[,3])))) +corrplot(cor(log(FEMALE %>% select(-Sex)))) +TEMP <- MALE %>% select(RATE=`36`,US_Rate,Year) %>% as.data.frame +TEST <- feols((RATE)~US_Rate+Year,TEMP) +TEST <- feols(RATE/US_Rate~Year,TEMP) +acf(as.numeric((TEMP[,'RATE']-predict(TEST)) %>% unlist)) +acf +TEST <- lm(log(`22`)~log(US_Rate)+Year,data=MALE) +predict(TEST +resid(TEST) +lm(`20`~Year+ +plot(as.vector(t(FEMALE %>% dplyr::select(`20`)))) +%>% pull(`22`) +?corrplot +plot(cor(FEMALE %>% select(-Sex))[2,]) +plot(cor(FEMALE %>% select(-Sex))[1,]) +MAT <- (cbind(abs(cor(FEMALE %>% select(-Sex))[2,]),abs(cor(FEMALE %>% select(-Sex))[1,]))) +plot(apply(MAT,1,min)) -%>% select(Year,Sex,County,Mort_Rate=`Age Adjusted Rate`) %>% filter(!is.na(County)) +MALE + +library(fixest) +MOD +MOD <- feols(Age_85~US_Rate,data=MALE) +acf(MALE[,"Age_85"]-predict(MOD)) +residuals(MOD) +MOD0 +resid(MOD0)