projection.R

library(classInt)
library(data.table)
library(dplyr)
library(lfe)
library(multcomp)
library(plotrix)
library(RColorBrewer)
library(readr)
library(reshape2)
library(sampling)
library(sp)
library(splines)



data <- read.csv('~/Desktop/GAM_final/GAM_monthly_dataframes/monthly_GAM_all_groups.csv')
est_us <- read.csv("~/Desktop/GAM_final/GAM_monthly_dataframes/est_us_new_updated.csv")
est_us <- est_us[, -which(names(est_us) == 'X')]

data$yr <- as.numeric(as.character(data$year))
yrs <- 2000:2019
ind <- data$yr%in%yrs
data <- data[ind,]
bus <- weighted.mean(data$suicide_rate,data$pop)

Tus <- read.csv('~/Downloads/NCC2018-master/inputs/projections/TemperatureProjections_us.csv') %>% as_data_frame()

Tus$deltaT <- Tus$deltaT 

quantile(Tus$deltaT,probs=c(0.025,0.5,0.975))  #statistics reported in paper
#quantile(Tmex$deltaT,probs=c(0.025,0.5,0.975))  #statistics reported in paper

dUS <- as.matrix(est_us)%*%matrix(Tus$deltaT,nrow=1)/bus*100

bxu <- c(dUS)
qu <- quantile(bxu,probs=c(0.025,0.5,0.975))
bxu[bxu<qu[1]] <- qu[1]
bxu[bxu>qu[3]] <- qu[3]

print(qu)

pp <- read_csv('~/Downloads/NCC2018-master/inputs/projections/PopulationProjections.csv')  #projections, medium variant
popc <- list()  #list to fill

## 2020 to 2100
p1 <- pp[pp[,3]=="United States of America",paste0(2020:2050)]
popc[[1]] <- as.numeric(gsub(" ","",unlist(p1)))*1000


yr=2020:2050

# excess suicides in US
dUS <- as.matrix(est_us)%*%matrix(Tus$deltaT,nrow=1)*12  #additional deaths per 100,000 per year by 2100 in US

# assume this rate increases linearly between 2000 and 2100 due to linear temperature increase
# looping over years, calculating effect for each bootstrap x model
ex <- array(dim=c(dim(dUS),length(yr)))
for (y in 1:length(yr)) {
  z <- (yr[y]-2020)/(2050-2020)
  ex[,,y] <- dUS*z*popc[[1]][y]/100000
}
exs <- apply(ex,c(1,2),sum)  #total cumulative deaths
exc_US <- exs

# define colors
alpha = 0.03
cll <- apply(sapply("orange", col2rgb)/255, 2, 
             function(x) rgb(x[1], x[2], x[3], alpha=alpha)) 

# define 2050 distributions for boxplotting
bxu <- c(exs)
qu <- quantile(bxu,probs=c(0.025,0.5,0.975))
bxu[bxu<qu[1]] <- qu[1]
bxu[bxu>qu[3]] <- qu[3]
print(qu)

#define median for each year
med <- apply(ex,3,median)
ub <- 30

#Panel (b)
plot(yr,cumsum(ex[1,1,])/1000,type="l",ylim=c(-8,ub),col=cll,las=1,ylab="excess suicides (thousands)",xlab="year",lwd=0.3,xlim=c(2020,2052),cex.axis=1.5,cex.lab=1.5,axes=F)
for (j in 1:dim(ex)[1]) {
  for (i in 1:dim(ex)[2]) {
    lines(yr,cumsum(ex[j,i,])/1000,col=cll,lwd=0.3)
  }}

lines(yr,cumsum(med)/1000,col="black",lwd=2)
abline(h=0,lty=2)
boxplot(bxu/1000,horizontal=F,range=0,at=2102,add=T,col="orange",axes=F,boxwex=3,lty=1)

mtext("Excess suicides by 2050, US", side=3, adj=0, line=0.5, cex=1.5, font=2)  #add title, left justified
axis(1,at=seq(2020,2050,10),labels=seq(2020,2050,10),cex.axis=1.5,las=1)
axis(2,at=seq(0,ub,10),labels=seq(0,ub,10),cex.axis=1.5,las=1)