简介：这是一套基于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个维度
GBD（Global 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") + # 用position="fill"看占比分布 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%不确定性区间（UI，Uncertainty 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输出图和官网不一致，先回官网确认“是不是我筛选维度选错了”，再检查代码。这个顺序能帮你省下大量排查时间。
本文还有配套的精品资源，点击获取
