利用GBD数据库做未来疾病负担预测

最近看到一篇文献，利用GBD数据库中的数据对疾病负担的发展趋势进行预测。觉得很有趣，就抽空用R语言复盘了一下里面关键的方法。

参考文献：
Changing trends in the disease burden of esophageal cancer in China from 1990 to 2017 and its predicted level in 25 years
DOI: 10.1002/cam4.3775

1. 相关R包的安装

主要用到的R包有3个，一个一个讲：

1.1 Nordpred package

用于做 Nordpred age-period-cohort (APC) analysis
注意，这个包不在r-cran上，需要利用remotes包从github上下载

install.packages("remotes")
library(remotes)
remotes::install_github("haraldwf/nordpred") 
library(nordpred)

如果下载比较慢，可以下载我存在百度云盘上的包
链接：https://pan.baidu.com/s/1DmooS88GtWVRflcelJ71sQ
提取码：jm4q
接着本地安装一下

install.packages("F:/R_packages/haraldwf-nordpred-ef83b2c.tar.gz", repos = NULL, type = "source")

1.2 INLA package

为了利用Integrated nested Laplace approximation (INLA) approach

install.packages("INLA", repos = "https://inla.r-inla-download.org/R/stable", dependencies = TRUE)

这个直接下比较麻烦，我试了好几次都半路掉线失败了，用我分享的包吧

#安装INLA前需要预装一个包 foreach
install.packages('foreach')
install.packages("F:/R_packages/INLA_21.02.23.zip", repos = NULL, type = "source")

1.3 BAPC package

为了利用Bayesian age-period-cohort models，可以产生 age-specific and age-standardized projected rates

#安装BAPC前需要预装一些包
install.packages("cmprsk")
install.packages("fanplot")
install.packages("Epi")
#接着下载安装BAPC包
install.packages("BAPC", repos = "http://R-Forge.R-project.org")
library(BAPC)

提醒：常见报错有rlang包版本低，重新卸了再装即可：

remove.packages('rlang')
install.packages('rlang')

同样，如果BAPC包下载不顺利，可以用我云盘上存的包进行本地安装：

install.packages("F:/R_packages/BAPC_0.0.34.tar.gz", repos = NULL, type = "source")

2. 数据准备

（1） GBD数据下载：http://ghdx.healthdata.org/gbdresults-tool
这里的案例用的是cancer medicine上发的中国食管癌文章，选用GBD 2017的数据；目前更新到GBD 2019的版本；在官网上按需要的条件一选，然后下载csv即可；

（2） 人口预测数据：https://population.un.org/wpp/Download/Standard/CSV/

（3） 人口结构标化，WHO2000-2025标准：https://seer.cancer.gov/stdpopulations/world.who.html

3. 模型讲解（可略过）

以下是我个人的理解，因为我只是个临床大夫，对模型的理解可能不够准确，所以这一段可以略过。

nordpred包和BAPC包的模型基础都是 **Age-Period-Cohort (APC)**模型，这个模型的理论基础是，发病率或死亡率和年龄结构和人口规模有关联。挖掘这种关联的基础还是广义线性模型（GLM）。nordpred包中有两种GLM：

1. power model即指数模型，当时把这个模型应用于癌症预测的时候（Moller et al.），认为指数为5的时候最好用，nordpred包里也采用了指数为5的固定指数模型；
2. Poisson log-link model，这个也是比较经典的了，log变换后满足泊松分布；

在上述基础上，APC模型可以在传统概率论框架下搭建，也可以在贝叶斯框架下进行，后者即为Bayesian age-period-cohort models（BAPC）。
有文献提到了一种在贝叶斯框架下，以随机指数模型作为链接函数的APC模型，讲得很仔细，具体的方程和评价方法都有，建议参考：
Verena Jürgens et al. A Bayesian generalized age–period–cohort power model for cancer projections；
DOI: 10.1002/sim.6248

INLA：integrated nested Laplace approximation (INLA) approach
贝叶斯推断中有一种比较有名的方法：马尔可夫链，但是这种方法效率不高，比较慢。而ILNA方法作为马尔科夫链的低计算密度的替代，用于隐高斯模型中进行近似贝叶斯推断。BAPC包中提供了在INLA框架下搭建的APC模型，官网描述：Package implementing Bayesian age-period-cohort models with the focus on projections. BAPC uses integrated nested Laplace approximations (INLA) for full Bayesian inference. BAPC generates age-specific and age-standardized projected rates. When interest lies in the predictive distribution, Poisson noise is automatically added.

4. 模型建立与运算

4.1 nordpred构建APC模型

nordpred的包详解见：https://rdrr.io/github/haraldwf/nordpred/man/nordpred.html
其中还自带一个示例，可以粘贴代码直接运行

4.1.1 数据整理

(1) 食管癌的数据，各个年龄段，1990-2017年，中国（这里以男性的发病率为例）
注意：

1. nordpred默认是18个年龄组，不符合的话会报错。由于食管癌15岁以下的没数据或都是为0，所以整成一组了。
2. nordpred包的模型最多预测5个值，如果你输入1990-2017年每年的数据，就可以预测2018-2022的值；所以，为了预测2042年的发病情况，这里把1990-2017年的数据按5年一组（1990-1992是3年一组），求了平均数；

esoph <-  fread('F:/nordpred/esoph.csv')
unique(esoph$age_name)
ages <- c("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")
### for incidence for male
esoph_ages <- esoph[age_name %in% ages & 
                      sex_name == 'Male' &
                      metric_name == 'Number' &
                      measure_name == 'Incidence',
                    .(age_id, age_name, year, val)]
#esoph_ages <- esoph_ages[order(age_id),]
esoph_ages_n <- dcast(data = esoph_ages, age_id + age_name ~ year)
rownames(esoph_ages_n) <- c("15-19", "20-24", "25-29", "30-34", 
                            "35-39", "40-44", "45-49", "50-54", 
                            "55-59", "60-64", "65-69", "70-74", 
                            "75-79", "80-84", "85-89", "90-94","95+")
esoph_ages_n["0-14",] <- 0
esoph_ages_n <- esoph_ages_n[order(esoph_ages_n$age_id),]

#这里需要按5年一组整理发病数， 我取了平均；因为有很多重复工作，所以写了个函数来完成，函数代码见文末
#参数分别是要整理的表格名，开始年份，结束年份，作为观察数据终点的年份
esoph_ages_g <- function_year5(esoph_ages_n , 1990, 2017, 2017)

整理后效果：

(2) 中国1990-2042的人口结构整理
注意：

1. 1990-2042年的数据按5年一组（1990-1992是3年一组），求了平均数；
2. 源数据的单位是每1000人，所以提取的时候*1000

population <- fread('F:/nordpred/WPP2019_PopulationByAgeSex_Medium.csv')
china_population <- population[Location == 'China',]
unique(china_population$AgeGrp)
china_population_1990_2017 <- china_population[Time %in% 1990:2042,
                                               .(Age_id = 1:21, AgeGrp, Time, PopMale = PopMale*1000)]
china_population_1990_2017_n <- dcast(data = china_population_1990_2017, Age_id +  AgeGrp ~ Time)
#calculation of 95+
china_population_1990_2017_n[22, 3:55] <- china_population_1990_2017_n[20, 3:48] + china_population_1990_2017_n[21, 3:55]
china_population_1990_2017_n[1, 3:55] <- china_population_1990_2017_n[1, 3:48] + china_population_1990_2017_n[2, 3:55] + china_population_1990_2017_n[3, 3:55]
china_population_1990_2017_n <- china_population_1990_2017_n[-c(2,3,20,21),]
rownames(china_population_1990_2017_n) <- c("0-14", "15-19", "20-24", "25-29", "30-34", 
                                            "35-39", "40-44", "45-49", "50-54", 
                                            "55-59", "60-64", "65-69", "70-74", 
                                            "75-79", "80-84", "85-89", "90-94","95+")
                                                                                       
china_population_1990_2017_g <- function_year5(china_population_1990_2017_n, 1990, 2042, 2017)

整理后的效果：

(3) 世界标准化人口

## data for age-standardization
age_stand <- read_file('F:/nordpred/AS.txt')
age_stand <- gsub(pattern = '\n', replacement = '\t', x = age_stand)
age_stand <- gsub(