最近看到一篇文献,利用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:
- power model即指数模型,当时把这个模型应用于癌症预测的时候(Moller et al.),认为指数为5的时候最好用,nordpred包里也采用了指数为5的固定指数模型;
- 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年,中国(这里以男性的发病率为例)
注意:
- nordpred默认是18个年龄组,不符合的话会报错。由于食管癌15岁以下的没数据或都是为0,所以整成一组了。
- 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的人口结构整理
注意:
- 1990-2042年的数据按5年一组(1990-1992是3年一组),求了平均数;
- 源数据的单位是每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(