GBD疾病负担数据分析：R语言从指标解读到可视化实战

简介这是一套基于R语言与GlobalBurdenR工具包分析全球疾病负担GBD数据的代码包面向公共卫生研究人员、流行病学分析者及有一定R基础的开发者用于解决GBD数据加载、合并、清洗、筛选及APC模型构建中的实际问题。压缩包共15个文件核心为6个R脚本apc_analysis.R、global_burden_utils.R、run_analysis.R等并包含2个CSV示例数据集、3张PNG结果图、1个HTML报告及说明文档体积仅75KB轻量易用目录中源代码、数据、输出分模块存放便于按需查阅和二次开发。目前已有361人学习下载。借助这份代码使用者既能掌握GlobalBurdenR工具包的操作流程和APC年龄-时期-队列模型的实现方法也能直接运行示例数据快速生成分析结果再替换自有数据完成扩展分析从而识别年龄、时期和队列效应对死亡率变化的影响辅助定位高风险群体为疾病负担评估与公共卫生政策制定提供可复现的研究工具R脚本中附有详细注释便于初学者理解每一步处理逻辑。1. GBD数据长什么样分析前必须先搞懂的5个指标和4个维度GBDGlobal Burden of Disease全球疾病负担研究数据由华盛顿大学健康指标与评估研究所IHME主导发布是目前全球范围内最系统的疾病负担数据库。它的核心价值不在于“某个病有多少人得”而在于把不同疾病的健康损失放在同一把尺子下比较——这就引出了几个绕不开的核心指标。做GBD数据分析前你必须先弄清楚五个关键指标的含义不然跑出来的代码再漂亮结论也可能是错的指标全称含义典型用途Incidence发病率某时间段内新发病例数评估疾病新发风险Prevalence患病率某时间点上现存病例数评估疾病负担存量Mortality死亡率因该病死亡的人数评估致命程度DALY伤残调整生命年因早死损失的生命年伤残损失的生命年综合评估总体负担YLD伤残损失健康生命年因患病或伤残损失的生命年评估非致命负担我第一次做GBD分析时上来就只取了死亡率数据结果做出来的疾病负担排名和官方报告差距很大。后来才明白DALY才是衡量疾病负担最核心的指标因为它同时覆盖了“死得太早”和“带病生存”两部分损失。如果你的分析目标是回答“哪个病对人群健康影响最大”首选DALY而不是单一死亡率。除了指标GBD数据还有四个重要的维度字段每个都会影响你的数据筛选策略location地区包含全球、国家、州省等多级区域划分。注意同一国家的国家级和省份级数据是分层的做国与国对比时不要混入子区域数据。year年份GBD 2019覆盖1990-2019年GBD 2021已扩展到1990-2021年。年份跨度不同版本有差异下载时注意版本说明。age年龄组从早产新生儿到80岁以上细分20多个年龄组也有年龄标化age-standardized后的汇总数据。sex性别男性、女性、两性合计。还有一个容易被忽略的字段是measure——区分你是取发病率、患病率还是DALY率。很多新手把“发病率”和“患病率”混成一个指标来跑结果趋势图出来两条线方向完全相反当场懵掉。R语言处理GBD数据的第一步就是把这些维度筛准确。有了这个基础认知下面我直接给出我实际在用的完整R语言分析流程从数据导入到最后出图。2. 从下载到清洗R语言处理GBD原始数据的标准流程2.1 数据导入与初始筛选IHME官网ghdx.healthdata.org下载的GBD数据通常是CSV压缩包格式文件较大我建议下载前先确认好自己要的分析范围尽量减少单次读取的数据量。library(tidyverse) library(data.table) # ---------- 数据导入 ---------- # 假设你下载的是gbd_data.csv放在工作目录下 gbd_raw - fread(gbd_data.csv, stringsAsFactors FALSE) # 查看数据结构 str(gbd_raw) # 实际字段一般为measure, location, sex, age, cause, year, val, upper, lower # ---------- 筛选关键维度 ---------- # 以两性合计 全年龄层 1990-2021年 国家层面为例 gbd_filtered - gbd_raw %% filter(sex Both, # 两性合计 age Age-standardized, # 年龄标化后的数据 year %in% 1990:2021, grepl(^[A-Z], location)) # 粗略过滤具体看你的location字段结构这里重点提醒筛选国家级数据时不要用location ! Global这种方式直接排除。GBD数据集中除了Global还有High-income这类收入分组、Latin America and Caribbean这类地区分组以及American Samoa这类非主权地区的细分区域。最稳妥的做法是把GBD官方提供的location hierarchy对照表下载下来筛选层级为国家或特定区域后再做分析。2.2 长表转宽表让数据结构匹配绘图逻辑GBD原数据是长表格式每个观测是一行但ggplot2绘图时不同场景需要长短表两种结构。比如你要画某个疾病从1990到2021年的DALY率变化趋势长表直接可以画但你要对比不同地区、不同疾病的负担时宽表反而方便。# ---------- 按指标拆分为宽表 ---------- gbd_wide - gbd_filtered %% select(measure, location, cause, year, val) %% pivot_wider(names_from measure, values_from val) # 保存清洗后的数据避免每次重复处理 saveRDS(gbd_wide, gbd_wide.rds)我习惯把清洗好的数据存成RDS格式因为GBD原始CSV动辄几百MB每次重新读取fread筛选非常浪费时间。RDS文件读取速度快一个量级而且保留因子类型、日期格式等R对象属性数据分析效率会高很多。2.3 计算变化率与年度变化百分比只是看绝对数值往往看不出趋势GBD分析里最常算的是相对于1990年的变化百分比和平均年度百分比变化AAPC。# ---------- 计算相对于1990年的变化百分比 ---------- gbd_full - gbd_wide %% group_by(location, cause, sex, age) %% mutate(year_1990 val[year 1990], pct_change (val - year_1990) / year_1990 * 100) %% ungroup()这里有个小坑如果你的数据是从2015年才开始有记录year 1990这一行会返回numeric(0)导致year_1990为NA后面的pct_change全部变成NA。稳妥的做法是先检查数据范围range(gbd_wide$year) # [1] 1990 2021确保你的时间覆盖范围完整再计算变化率。3. 时间趋势分析用ggplot2绘制疾病负担变化图的核心代码3.1 单疾病多地区趋势图GBD最常见的可视化需求之一是看某个疾病在几个代表性国家/地区的疾病负担随时间的变化趋势。我用的是这张图library(ggplot2) library(ggthemes) # 以缺血性心脏病Ischemic heart disease为例 ihd_trend - gbd_full %% filter(cause Ischemic heart disease, location %in% c(China, United States of America, India, Global)) # 设定地区顺序保证图例顺序可控 ihd_trend$location - factor(ihd_trend$location, levels c(China, United States of America, India, Global)) ggplot(ihd_trend, aes(x year, y val, color location)) geom_line(linewidth 1.1) scale_y_continuous(limits c(0, NA), expand expansion(mult c(0, 0.05))) labs(title Age-standardized DALY rate of Ischemic Heart Disease, 1990-2021, x Year, y DALY rate per 100,000, color NULL) theme_minimal(base_size 14) theme(legend.position bottom, plot.title element_text(face bold, hjust 0.5))这里注意geom_line()替代了旧版的geom_line(aes(group location))——新版ggplot2只要你映射了color或linetype就会自动按该变量分组不必再手动指定group。另外y轴建议用limits c(0, NA)强制从0开始否则疾病负担变化的视觉对比会被压缩得很夸张。3.2 按年龄组拆分的堆叠面积图如果要看疾病负担的年龄分布迁移比如这个病从老年人为主变为中年人为主堆叠面积图很直观# 取具体年龄组数据注意此处age不是Age-standardized gbd_age - gbd_raw %% filter(sex Both, location China, cause Ischemic heart disease, measure DALYs, age ! Age-standardized, year %in% c(1990, 2000, 2010, 2019)) %% mutate(age factor(age, levels unique(gbd_raw$age))) ggplot(gbd_age, aes(x year, y val, fill age)) geom_area(position fill) # 用positionfill看占比分布 scale_y_continuous(labels scales::percent) labs(title Contribution of age groups to IHD DALYs, China, x NULL, y Proportion, fill Age group) theme_minimal()position fill是把堆叠面积图归一化成百分比用来观察年龄结构变化非常方便。如果想看绝对值的年代变化用默认position stack。3.3 多个疾病放在同一张对比图我平时还会用facet_wrap把几种主要疾病并列展示这样一张图就能看出哪些疾病在下行、哪些在上行top_causes - c(Ischemic heart disease, Stroke, COPD, Diabetes, Low back pain) ggplot(gbd_full %% filter(cause %in% top_causes, location China, sex Both), aes(x year, y val)) geom_line(color #2E86AB, linewidth 1) facet_wrap(~ cause, scales free_y, ncol 2) labs(x Year, y DALY rate per 100,000) theme_minimal(base_size 13)这里scales free_y很关键——不同疾病的负担量级差异巨大比如下背痛的DALY率可能是缺血性心脏病的数倍共用一个y轴会把低值疾病压成一条直线。4. 森林图实战用forestploter包做疾病负担分点估计与区间展示4.1 数据准备怎么从GBD里提取点估计置信区间森林图是GBD数据分析出镜率极高的图尤其在展示多重对比不同地区、不同性别、不同年份的相对风险时。GBD数据自带upper和lower字段就是每个估计值的95%不确定性区间UIUncertainty Interval画森林图的原料本身就是现成的。# ---------- 准备森林图数据 ---------- forest_data - gbd_full %% filter(measure DALYs | measure Deaths, location %in% c(China, Japan, United States of America), cause %in% c(Ischemic heart disease, Stroke, COPD), sex %in% c(Male, Female)) %% select(measure, location, cause, sex, year, val, upper, lower) %% filter(year 2019) %% mutate(est_text sprintf(%.1f (%.1f-%.1f), val, lower, upper))需要注意森林图一般展示的是某一时点的对比不是全时段数据。截取一个年份比如2019或最新一期会让图面清晰得多。如果想把时间趋势也塞进森林图更高级的变体是把不同年份作为单独行分行排列。4.2 用forestploter绘制分组森林图forestploter是我用过圈子比较小但功能很顺手的R包它的优势是可直接用数据框控制左侧指标列右侧自动生成森林图排版控制力强。先安装# install.packages(forestploter) library(forestploter)然后整理成forestploter需要的格式# 准备左侧文字列可包含多个列 forest_plot_df - forest_data %% mutate(Region paste0(location, , sex), Cause cause, DALY rate est_text) %% select(Region, Cause, DALY rate, val, lower, upper) # 创建森林图对象 p - forest(forest_plot_df, est c(forest_plot_df$val, forest_plot_df$lower, forest_plot_df$upper), # 这个est参数在最新版本中更推荐用list方式传递见下方 ci_column 3, # 从第几列开始画森林图区域 ref_line NA, xlim c(0, 500), ticks_at seq(0, 500, 100), theme forest_theme(base_size 12))上面代码里的est参数在不同版本的forestploter中写法稍有差异。最新的forestploter推荐直接传入包含点估计、下限、上限的三列数据框。我更常用的写法是把val/lower/upper这三列单独抽出来传入# 更稳健的写法 forest_plot_df - forest_data %% mutate(Region paste0(location, , sex), Cause cause, DALY rate (95% UI) sprintf(%.1f (%.1f-%.1f), val, lower, upper)) %% select(Region, Cause, DALY rate (95% UI)) p - forest(forest_plot_df, est list(forest_data$val, forest_data$lower, forest_data$upper), ci_column 3, ref_line 0, xlim c(0, 400), ticks_at seq(0, 400, 50), theme forest_theme(base_size 12, footnote Data source: GBD 2019))注意ci_column指的是森林图从第几列开始占据的垂直区域。如果你左侧有三列文字ci_column通常设置为4或3取决于你希望在图右侧还是文字下方显示点估计。我一般把数值列放在左侧森林图只展示区间范围阅读更直观。4.3 一张图表多指标组合如果你想让森林图同时展示DALY率和死亡率可以在行标签加上分组变量forest_plot_df - forest_data %% mutate(Group paste(cause, |, measure), Region paste0(location, , sex), Disease burden sprintf(%.1f (%.1f-%.1f), val, lower, upper)) %% select(Group, Region, Disease burden) p - forest(forest_plot_df, est list(forest_data$val, forest_data$lower, forest_data$upper), ci_column 3, ref_line 0, xlim c(0, 500), ticks_at seq(0, 500, 100), theme forest_theme(base_size 12))画出来后保存图片我用的是ggsave(gbd_forest.png, p, width 10, height 8, dpi 300)实际跑下来forestploter的输出结果是ggplot对象所以ggsave()可以用。但如果用plot(p)预览会发现它在绘图设备上直接渲染了森林图的静态版本两者行为有差异——建议统一用ggsave()保存格式和分辨率都可控。5. 选代表性疾病标准化从GBD海量数据中提炼一篇论文级图表5.1 多疾病top10排序Bar图GBD覆盖几百种疾病和伤害分析时最常见的目标是“找出某地区负担最重的疾病并给出排序”。这一步我用一个极简的代码块完成# 选出2019年中国DALY率最高的10个三级病因 top10 - gbd_full %% filter(location China, measure DALYs, age Age-standardized, sex Both, year 2019) %% arrange(desc(val)) %% slice_head(n 10) top10$cause - factor(top10$cause, levels rev(top10$cause)) # 反向排序使y轴图上从大到小 ggplot(top10, aes(x cause, y val, fill cause)) geom_col(show.legend FALSE) coord_flip() labs(title Top 10 causes of DALYs in China, 2019, x NULL, y Age-standardized DALY rate per 100,000) theme_minimal()注意这里需要先按val排序再用factor(..., levels rev(...))反转因子顺序这样在coord_flip()之后最大的值出现在图的最上方。不反转的话ggplot2默认将第一个因子放在底部要么反过来图会很难看。5.2 性别差异分析男性和女性的疾病负担对比GBD分析中性别差异是高频分析点。我常做的是“男vs女差异百分比”用来筛选出性别差异最大的疾病gbd_sex - gbd_full %% filter(age Age-standardized, year 2019, cause %in% top10$cause) %% select(location, cause, sex, val) %% pivot_wider(names_from sex, values_from val) %% mutate(diff_pct (Male - Female) / Female * 100) ggplot(gbd_sex, aes(x reorder(cause, diff_pct), y diff_pct, fill diff_pct 0)) geom_col(show.legend FALSE) coord_flip() geom_hline(yintercept 0, linetype dashed, color grey40) scale_fill_manual(values c(#E07A5F, #3D405B)) labs(title Sex difference in age-standardized DALY rate, 2019, x NULL, y Male-Female difference (%)) theme_minimal()diff_pct (Male - Female) / Female * 100直接给出“男性比女性高百分之多少”的直观数值。这里的reorder(cause, diff_pct)是按差值排序省掉手动改因子顺序的功夫。图中正值男性负担更高用深蓝色负值女性更高用红色一眼就能看出性别差异最大的前几个病种。5.3 年龄-时期-队列的初步探索用热图看年龄和时间的交叉最后补充一个我很常用的交叉表热图用来快速探索“年龄组×年份”两个维度下疾病负担的分布模式gbd_year_age - gbd_raw %% filter(location China, cause Ischemic heart disease, measure DALYs, sex Both, age ! Age-standardized, year %% 5 0) %% # 每5年取一次降低密度 select(year, age, val) ggplot(gbd_year_age, aes(x factor(year), y age, fill val)) geom_tile() scale_fill_viridis_c(option C) labs(x Year, y Age group, fill DALY rate) theme_minimal(base_size 13) theme(axis.text.x element_text(angle 45, hjust 1))year %% 5 0可以快速降低横轴密度否则1990到2021每年一个格子标签会挤成一团。这种热图做出来非常像发表级别的Figure 1适合放在论文或报告开篇做全局概览。这几种图我组合起来用先热图看全局分布再映射趋势线看时间演化然后用森林图横向对比不同病种/地区/性别最后用top10条形图锁定重点疾病。一套流程下来一篇文章的核心图表基本够用了。6. 我踩过几次坑之后总结出来的GBD分析避坑清单6.1 版本选择与单位换算GBD官网每轮更新会同时发布多个版本如GBD 2019、GBD 2021不同版本中的疾病分类编码可能略有调整。分析时必须下载单一版本并记录版本号否则不同版本数据混用会导致同一疾病的发病率数值差异过大。另一点是率值单位GBD默认rate是“每10万人”的率有些场景你可能想换算成“每100万人”或“每千人”的率换算因子不要搞错。我在文中提供的示例代码全部按官方默认的每10万人率值来计算。6.2 不要直接用All causes做内部比较GBD数据里的cause字段包含层级关系All causes是顶层然后按病因层级分level 1/2/3/4。做疾病负担对比时如果直接混用不同层级的cause字段会重复计算负担。比如Cardiovascular diseases下包含Ischemic heart disease和Stroke如果你既选了父级又选子级加总后会大于真实负担。建议统一用level 3或level 4的病因编码并检查是否有重复行。我自己的检查习惯是# 找出重复组合 duplicates - gbd_filtered %% count(measure, location, age, sex, cause, year) %% filter(n 1) # 如果有重复需要去重或回看原始数据 nrow(duplicates)如果这个数量不为0说明筛选逻辑有漏洞得回头检查。6.3 处理age-standardized和“All ages”时的差异GBD数据中age字段常见取值包括Age-standardized年龄标化率、All ages全年龄层汇总值、Under 5、5-9 years一直到80 plus等具体年龄组。使用age Age-standardized得到的是经过世界标准人口加权后的标化率适合跨国家/跨年份比较而使用age All ages得到的是该地区所有年龄人口的原始负担总量表述为“总例数”而非“率”。两者经常被搞混错误出现最多的场景就是在标题里把“标化率”写成“实际发病例数”。6.4 保存中间结果别把原始文件反复读取在我刚开始用R处理GBD的几百兆CSV时每次改筛选条件都要重新读一遍文件RStudio被卡死是常事。后面养成习惯数据预处理一次后立即存RDS后续所有做图、建模都从RDS读取运行速度快非常多。这个习惯后来也延伸到其他大数据的分析工作中算是做数据项目最基础也最实用的工程化经验。saveRDS(gbd_data_clean, gbd_clean_2021.rds) # 下次直接用 gbd - readRDS(gbd_clean_2021.rds)最后再分享一个小技巧GBD的GBD Compare工具可以直接在线预览结果如果你的R输出图和官网不一致先回官网确认“是不是我筛选维度选错了”再检查代码。这个顺序能帮你省下大量排查时间。本文还有配套的精品资源点击获取