# setwd('path/to/GBD_Prediction') ##???ù???·??  # 已注释：请按本地环境设置工作目录
library(BAPC)
library(INLA)
library(nordpred)
library(reshape)
library(data.table)
library(tidyr)
library(tidyverse)
library(epitools)
library(ggplot2)

EC <-  read.csv('EC_predict.csv')
age_stand <- read.csv('GBD2019 world population age standard.csv')[-c(1:3),]
#### 用于疾病数据
ages <- c("20 to 24", "25 to 29",
          "30 to 34", "35 to 39", "40 to 44", "45 to 49", "50 to 54", "55 to 59",
          "60 to 64", "65 to 69", "70 to 74", "75 to 79", "80 to 84", "85 to 89",
          "90 to 94", "95 plus")

#### 读取人口学数据用
ages_2 <- c("<1 year","1 to 4", "5 to 9","10 to 14", "15 to 19","20 to 24", "25 to 29",
            "30 to 34", "35 to 39", "40 to 44", "45 to 49", "50 to 54", "55 to 59",
            "60 to 64", "65 to 69", "70 to 74", "75 to 79", "80 to 84", "85 to 89",
            "90 to 94", "95 plus")

###### 数据运算分析用
ages_3 <- c("0 to 4", "5 to 9","10 to 14", "15 to 19","20 to 24", "25 to 29",
           "30 to 34", "35 to 39", "40 to 44", "45 to 49", "50 to 54", "55 to 59",
           "60 to 64", "65 to 69", "70 to 74", "75 to 79", "80 to 84", "85 to 89",
           "90 to 94", "95 plus")
#
wstand <- c(age_stand$std_population[1:2] %>% as.numeric() %>% sum(),
            age_stand$std_population[3:21] %>% as.numeric())/sum(age_stand$std_population[1:21])

### for incidence for Male and female
EC_Male_incidence <- subset(EC,age %in% ages &
                              sex == 'Male' &
                              metric == 'Number' &
                              measure == 'Incidence' &
                              location== 'Global')[,c(3,4,7,8)]
EC_Male_incidence_n <- reshape2::dcast(data = EC_Male_incidence, year~age, value.var = "val")
rownames(EC_Male_incidence_n) <- EC_Male_incidence_n$year
EC_Male_incidence_n <- EC_Male_incidence_n[,-1]
##### 补齐0-19年龄段数据
EC_Male_incidence_n[,'0 to 4'] <- 0
EC_Male_incidence_n[,'5 to 9'] <- 0
EC_Male_incidence_n[,'10 to 14'] <- 0
EC_Male_incidence_n[,'15 to 19'] <- 0
EC_Male_incidence_n <- EC_Male_incidence_n[,c(17:20,1:16)]  ####将年龄段从小到大排列

EC_Male_incidence_n <- apply(EC_Male_incidence_n, c(1,2), as.integer) %>% as.data.frame()
EC_Male_incidence_n <- apply(EC_Male_incidence_n, c(1,2), round) %>% as.data.frame()  #### 数值取整


EC_Female_incidence <- subset(EC,age %in% ages &
                                sex == 'Female' &
                                metric == 'Number' &
                                measure == 'Incidence' &
                                location == 'Global')[,c(3,4,7,8)]
EC_Female_incidence_n <-  reshape2::dcast(data = EC_Female_incidence, year~age, value.var = "val")
rownames(EC_Female_incidence_n) <- EC_Female_incidence_n$year
EC_Female_incidence_n <- EC_Female_incidence_n[,-1]
##### 补齐0-19年龄段数据
EC_Female_incidence_n[,'0 to 4'] <- 0
EC_Female_incidence_n[,'5 to 9'] <- 0
EC_Female_incidence_n[,'10 to 14'] <- 0
EC_Female_incidence_n[,'15 to 19'] <- 0
EC_Female_incidence_n <- EC_Female_incidence_n[,c(17:20,1:16)]

EC_Female_incidence_n <- apply(EC_Female_incidence_n, c(1,2), as.integer) %>% as.data.frame()
EC_Female_incidence_n <- apply(EC_Female_incidence_n, c(1,2), round) %>% as.data.frame() #### 数值取整

##### 读取人口学数据
dirname <- dir("path/to/GBD_Population")
file <- paste0("path/to/GBD_Population",dirname)
var_name <- c("location_name","sex_name","year_id","age_group_name","val")

GBD_population  <-  as.data.frame(matrix(nrow=0,ncol=length(var_name)))
names(GBD_population)=var_name
for (a in file) {
  data <- fread(a) %>% select(var_name) %>%
    filter(age_group_name %in% ages_2)
  GBD_population <- rbind(GBD_population,data)
}
GBD_population$sex_name[GBD_population$sex_name=='both'] <- 'Both'
GBD_population$sex_name[GBD_population$sex_name=='male'] <- 'Male'
GBD_population$sex_name[GBD_population$sex_name=='female'] <- 'Female'
GBD_population <- GBD_population[!duplicated(GBD_population),]

#### 读取预测人口数据
prediction_var_name <- c("location_name","sex","year_id","age_group_name","val")
GBD_population_prediction <- fread('IHME_POP_2017_2100_POP_REFERENCE_Y2020M05D01.csv') %>%
  select(prediction_var_name) %>%
  filter(year_id %in% 2020:2030)
names(GBD_population_prediction) <- var_name
GBD_1year <- GBD_population_prediction %>% filter(age_group_name %in% c("Early Neonatal","Late Neonatal", "Post Neonatal")) %>%
             group_by(location_name,sex_name,year_id) %>% summarise(val=sum(val)) %>%
             mutate(age_group_name="<1 year") %>% .[,c(1:3,5,4)]
GBD_population_prediction <- GBD_population_prediction %>% filter(!(age_group_name %in% c("Early Neonatal","Late Neonatal", "Post Neonatal"))) %>%
                             rbind(GBD_1year)
##### 合并现人口数+预测人口数
GBD <- rbind(GBD_population,GBD_population_prediction)
#### ?ϲ?0-19??????Ⱥ????
GBD_age4 <- GBD %>% subset(age_group_name %in% c("<1 year","1 to 4")) %>%
  group_by(location_name,sex_name,year_id) %>%
  summarize(val=sum(val)) %>%  mutate(age_group_name='0 to 4') %>% .[,c(1:3,5,4)]
GBD <- rbind(GBD,GBD_age4)
GBD <- subset(GBD, age_group_name %in% ages_3) %>%
  mutate(age_group_name=fct_relevel(age_group_name,ages_3)) %>%
  arrange(age_group_name)

#####################  读取男性和女性数据
GBD_Global_Male <- subset(GBD,location_name=='Global' & sex_name == 'Male')
GBD_Global_Female <- subset(GBD,location_name=='Global' & sex_name == 'Female')

GBD_Global_Male_n <- dcast(data = GBD_Global_Male, year_id ~ age_group_name,value.var = c("val"))
GBD_Global_Female_n <- dcast(data = GBD_Global_Female, year_id ~ age_group_name,value.var = c("val"))

rownames(GBD_Global_Male_n) <- GBD_Global_Male_n$year_id
rownames(GBD_Global_Female_n) <- GBD_Global_Female_n$year_id

GBD_Global_Male_n <- GBD_Global_Male_n %>% .[,-1] %>%
                       apply(c(1,2), as.numeric) %>%
                       apply(c(1,2), round) %>% as.data.frame()


GBD_Global_Female_n <- GBD_Global_Female_n %>% .[,-1] %>%
                     apply(c(1,2), as.numeric) %>%
                     apply(c(1,2), round) %>% as.data.frame()

GBD_Global_Both_n <- GBD_Global_Female_n + GBD_Global_Male_n
#

# 补全疾病数据
EC_pro <- matrix(data = NA, nrow = 2030-2019, ncol = ncol(GBD_Global_Male_n)) %>% as.data.frame()
rownames(EC_pro) <- seq(2020,2030,1)
colnames(EC_pro) <-  names(EC_Male_incidence_n)

EC_Male_incidence_n  <- rbind(EC_Male_incidence_n , EC_pro)
EC_Female_incidence_n  <- rbind(EC_Female_incidence_n , EC_pro)

####### ģ??Ԥ??
Male_esoph <- APCList(EC_Male_incidence_n, GBD_Global_Male_n, gf = 5)
Male_bapc_result <- BAPC(Male_esoph, predict = list(npredict = 10, retro = T),
                         secondDiff = FALSE, stdweight = wstand, verbose = F)

Female_esoph <- APCList(EC_Female_incidence_n, GBD_Global_Female_n, gf = 5)
Female_bapc_result <- BAPC(Female_esoph, predict = list(npredict = 10, retro = T),
                           secondDiff = FALSE, stdweight = wstand, verbose = F)

#####
Male_proj <- agespec.proj(x = Male_bapc_result) %>% as.data.frame()   ### ????Ԥ????Ⱥ
Male_proj_mean <- Male_proj[,colnames(Male_proj) %like% 'mean']
colnames(Male_proj_mean) <- age_3

Female_proj <- agespec.proj(x = Female_bapc_result) %>% as.data.frame()   ### ????Ԥ????Ⱥ
Female_proj_mean <- Female_proj[,colnames(Female_proj) %like% 'mean']
colnames(Female_proj_mean) <- age_3

Both_proj_mean <- Female_proj_mean+ Male_proj_mean

##### ??ͬ?????㷢????
Male_rate <- agespec.rate(x = Male_bapc_result) %>% as.data.frame()
Male_rate_mean <- Male_rate[,colnames(Male_rate) %like% 'mean']*100000
colnames(Male_rate_mean) <- age_3

Female_rate <- agespec.rate(x = Female_bapc_result) %>% as.data.frame()
Female_rate_mean <- Female_rate[,colnames(Female_rate) %like% 'mean']*100000
colnames(Female_rate_mean) <- age_3

Both_rate_mean <- Both_proj_mean/GBD_Global_Both_n*100000


######  ????????У?????ķ?????
Male_ASR <- agestd.rate(x = Male_bapc_result) %>% as.data.frame()
Male_ASR$mean <- Male_ASR$mean*100000
Male_ASR$year <- rownames(Male_ASR)

Female_ASR <- agestd.rate(x =Female_bapc_result) %>% as.data.frame()
Female_ASR$mean <- Female_ASR$mean*100000
Female_ASR$year <- rownames(Female_ASR)

year_index <- 1990:2030
Both_ASR <- matrix(nrow = 0, ncol = 2) %>% as.data.frame()
names(Both_ASR) <- c("crude.rate","adj.rate")
#i=1
for (i in 1:(2030-1989)) {
  asr = ageadjust.direct(count = Both_proj_mean[i,], pop = GBD_Global_Both_n[i,],
                         stdpop = wstand)
  Both_ASR[i,1:2] <- round(asr[1:2]*10^5,2)
}
Both_ASR$year <- 1990:2030


Male_sum_year <- apply(Male_proj_mean, 1, sum) %>% as.data.frame()
colnames(Male_sum_year) <- 'number'
Male_sum_year$year <- rownames(Male_sum_year)


Female_sum_year <- apply(Female_proj_mean, 1, sum) %>% as.data.frame()
colnames(Female_sum_year) <- 'number'
Female_sum_year$year <- rownames(Female_sum_year)

Male_sum_year <- apply(Male_sum_year, 2, as.numeric) %>% as.data.frame()
Female_sum_year <- apply(Female_sum_year, 2, as.numeric) %>% as.data.frame()
Both_proj_mean <- Female_proj_mean + Male_proj_mean
Both_sum_year <- Male_sum_year + Female_sum_year
Both_sum_year$year <- 1990:2030


plotBAPC(Male_bapc_result, scale=10^5, type = 'ageStdProj', showdata = T)