简介一份用于全球疾病负担GBD数据库贝叶斯年龄-时期-队列BAPC分析的R包资源主要面向公共卫生研究人员和数据爱好者帮助完成疾病发病率和死亡率的趋势分解与短期预测。压缩包内含41个文件大小约57KB其中13个R源代码文件承担核心统计函数13个Rd文档提供函数接口说明3个Rda数据集用于实例演示另有RMarkdown、Markdown等文件辅助理解使用流程整体体量小巧但结构清晰。资源以Nordpred包及其相关扩展代码为主体实现年龄、时期、队列效应估计、预测及可视化绘图并附带挪威男性癌症数据作为教学实例同时提供数据生成脚本与原始文本数据便于研究者迁移至GBD各类疾病负担分析场景。已有2089人学习适合需要快速上手BAPC方法、减少环境配置成本并希望结合示例数据理解建模思路的读者。1. 用 R 包跑 GBD 数据库的 BAPC 分析从下载到出图先把这四件事搞清楚这几年做疾病负担测算的人越来越多GBD 数据库的现成数据好拿但真正想从“描述发病率和死亡率”往前再走一步做未来趋势预测绕不开 BAPC 这个模型。BAPC 全称是 Bayesian Age-Period-Cohort analysis也就是贝叶斯年龄-时期-队列分析它能把年龄结构、时期趋势和出生队列效应拆开再外推未来几年的疾病负担。这个分析在 R 里不是一个大而全的官方包而是由若干个 R 包配合完成有的人卡在数据格式有的人卡在 INLA 安装还有的人跑出来结果图长得很怪但不知道怎么排查。我这次分享的是我自己整理过的 R 包清单和配套脚本覆盖从 GBD 数据提取、年龄分组、APC 模型拟合到 BAPC 预测出图的完整流程。适合两类人一类是刚接触 GBD 的医学生或公卫从业者另一类是把 Stata 或 SPSS 用熟、但想在 R 里跑贝叶斯预测的进阶用户。下面直接聊每个环节怎么落地参数怎么设以及我踩过的几个实实在在的坑。2. BAPC 分析到底依赖哪些 R 包选型理由与版本边界2.1 核心包拆解INLA、BAPC、AgePeriodCohort 各自的作用BAPC 在 R 里的实现并不像 lm() 那样一个函数搞定它实际是一条链条。链条最底层是 INLA负责贝叶斯计算的近似推断中间层是 BAPC 包把年龄-时期-队列模型写成 INLA 能吃的公式最上层是 AgePeriodCohort 包提供绘图和结果的整理函数。这三个包缺一个都不行而且版本之间还有隐藏的兼容性问题。INLA 是整条链的引擎。它不是 CRAN 上的普通包而是需要从 R-INLA 官方仓库安装因为它的计算内核依赖外部 C 库和系统环境。我最早在 R 4.0 上直接 install.packages(INLA) 当然失败后来才知道必须用 inlainstall 的仓库地址。BAPC 包本身依赖 INLA 的结果对象它不自己算 MCMC而是把数据整理成 INLA 需要的格式然后调用 inla()。AgePeriodCohort 则负责可视化和一些描述性统计它的预测区间作图函数是很多人最后出图时最依赖的。选型理由一句话如果只是做简单的年龄标化率趋势用不了 BAPC但如果要分解年龄-时期-队列效应并外推预测BAPC 是目前 R 社区里相对成熟、文档多、可复现的方案。相比之下手动写 JAGS 或 Stan 模型的灵活度更高但你要自己处理年龄分组矩阵、队列共线性、预测期设计矩阵这些细节对大多数人来说门槛偏高。2.2 安装顺序与系统环境先装底层再装上层我一般建议按 INLA → BAPC → AgePeriodCohort 的顺序装先解决底层依赖再往上装。INLA 的安装用下面这段代码注意不要用 install.packages而是用官方仓库的 URL# 安装 INLA选择测试版或稳定版一般用 stable 即可 install.packages(INLA, reposc(getOption(repos), INLAhttps://inla.r-inla-download.org/R/stable), depTRUE) # 安装 BAPC 和 AgePeriodCohort这两个在 CRAN 上 install.packages(BAPC) install.packages(AgePeriodCohort)这段代码里repos 参数把 INLA 的镜像源追加到当前源列表这样 install.packages 才会从 R-INLA 官方源拉包。depTRUE 意思是一并安装依赖包否则可能会漏掉 Matrix、sp 等底层包。安装完后别急着直接用先跑一句 packageVersion(INLA) 确认版本号常见能正常用的是 22.x 以后版本太老的版本在 R 4.2 以上容易出现矩阵维度不匹配的报错。装完三个包后还有一个容易忽略的环节INLA 在运行时会调用系统级 BLAS/LAPACK 库如果你的 R 是 Windows 下 32 位和 64 位共存的环境建议统一用 64 位 R 并且关闭多线程冲突。我见过有人跑 inla() 的时候 RStudio 直接死掉最后发现是系统 PATH 里同时存在多个 Rtools 版本导致编译环境混乱解决办法是只保留一个 Rtools 并重启 R 会话。3. 把 GBD 数据库数据变成 BAPC 能吃的格式年龄分组与矩阵构造3.1 GBD 数据下载后的结构特征先看懂三个维度GBD 数据库导出的数据通常是一个长格式表格每一行代表某个国家、某个年份、某个年龄组、某个性别、某种疾病的患病率或死亡率。以 cause of death 或 incidence 为例列里常见的有 location_name、year、age_name、sex、val、upper、lower 这些字段。BAPC 分析前你必须把长表拆成三维矩阵年龄组作为行年份作为列队列作为对角线方向的衍生结构。这里有一个关键点GBD 的年龄组并不是均匀分布的。GBD 2019 里年龄组通常按 5 岁一组划分从 0-4 岁到 95 岁以上但有的疾病负担研究会用 0-14 岁、15-49 岁、50 岁以上这种粗略分组。BAPC 的年龄-时期-队列模型对年龄组的宽度和一致性很敏感如果你在数据准备阶段混用了不同宽度的年龄组预测结果会出现明显的锯齿。我的建议是初始整理尽量不动 GBD 原始的年龄组保持 5 岁一组再用 merge 或 pivot 操作把性别合并或分层。3.2 数据清洗与宽表转换pivot_wider 与年龄按序排列下面这段代码是我常用的清洗流程假设你已经从 GBD 官网下载了 CSV读进来之后先过滤出你要的 location、sex 和 measure然后转成宽表library(dplyr) library(tidyr) # 假设 dat 是 GBD 导出数据列名按照下载时的标准列名 dat - read.csv(gbd_incidence.csv, stringsAsFactors FALSE) # 筛选只保留中国、全部年龄、两性合计、发病数 dat_sub - dat %% filter(location_name China, sex_name Both, measure_name Incidence, metric_name Number) %% select(year, age_name, val) %% mutate(age_name factor(age_name, levels c(0-4 years, 5-9 years, 10-14 years, 15-19 years, 20-24 years, 25-29 years, 30-34 years, 35-39 years, 40-44 years, 45-49 years, 50-54 years, 55-59 years, 60-64 years, 65-69 years, 70-74 years, 75-79 years, 80-84 years, 85-89 years, 90-94 years, 95 years))) # 转成宽表每行是一个年龄组每列是一个年份 wide - dat_sub %% pivot_wider(names_from year, values_from val) %% arrange(age_name) %% select(-age_name) %% as.matrix() rownames(wide) - levels(dat_sub$age_name)这段代码里有两个地方容易翻车。第一是 age_name 的 levels 必须按时间顺序排列R 的 factor 默认按字母排序如果你不显式指定 levels转成矩阵后 10-14 岁会排在 0-4 岁前面整个模型的年龄效应曲线就会对不上号。第二是 pivot_wider 之后可能出现 NA 值原因是某些年份 GBD 没有对应年龄组的数据比如最早年份之前没有 95 岁以上记录这时候需要明确用 0 填充还是插补我一般倾向用前向填充或干脆删掉缺失对应的年份列避免 INLA 在拟合阶段因 NA 报错。3.3 构造年龄-时期-队列索引矩阵别用手工循环BAPC 包内部需要三个索引向量年龄组的序号、年份的序号、出生队列的序号。出生队列的计算规则是 period(年份) 减去 age(年龄组起点)再映射到一个连续的整数编号。这里最容易出错的是年龄组用区间的起始年龄还是中点年龄。BAPC 的默认逻辑是使用年龄组的编号本身作为年龄索引但队列编号需要你预先算好。# 宽表 wide 已经准备好行数为年龄组数列数为年份数 n_age - nrow(wide) n_year - ncol(wide) year_start - as.numeric(colnames(wide)[1]) year_end - as.numeric(colnames(wide)[ncol(wide)]) # 每个年龄组的起始年龄这里假设是 5 岁一组第一组 0-4 岁 age_start - seq(0, by 5, length.out n_age) # 构造每个格的队列编号队列 年份 - 年龄组起始值 cohort_index - matrix(NA, nrow n_age, ncol n_year) for (i in 1:n_age) { for (j in 1:n_year) { cohort_index[i, j] - as.numeric(colnames(wide)[j]) - age_start[i] } } # 队列编号需要从 1 开始连续编号不能有负数或跳号 cohort_unique - sort(unique(as.vector(cohort_index))) cohort_map - setNames(seq_along(cohort_unique), cohort_unique) cohort_id - apply(cohort_index, 2, function(x) cohort_map[as.character(x)])这段代码本质是把二维的年龄-年份格子影射成一个一维的队列向量。为什么要手动改编号因为 BAPC 包在拟合时会把队列效应当作一个随机效应处理如果队列编号是 1950、1951、1960 这种不连续的原始年份INLA 会默认当成连续型协变量而不是分组效应导致结果完全错误。映射成 1, 2, 3... 之后INLA 才能正确识别这是 46 个独立的队列水平。循环在这里并不慢因为年龄组一般不超过 20 个年份最多三四十年但如果年份跨度到 100 年可以改用 outer 函数加速。4. 跑通 BAPC 模型从 INLA 公式到未来预测4.1 使用 BAPC 包的推荐流程先读官方示例再改数据BAPC 包自带一个示例数据集是荷兰的癌症发病数据。我每次拿到新数据第一步不是直接替换而是跑一遍它的示例脚本确认当前 R 环境里的 INLA、BAPC、AgePeriodCohort 三个包都能协同工作。因为这三个包的版本更新节奏不一致有时候 BAPC 更新后调用 INLA 的对象结构方式变了你不跑示例根本看不出来。library(BAPC) data(ICD) # 这个内置数据集是 ICD 编码的荷兰癌症数据包含了 age, period, cohort 三个维度 # 先跑一个默认模型验证环境没问题 bapc_model - BAPC(ICD, model APC, predict TRUE, agegroup 5year, type inla)BAPC() 函数里最核心的参数是 model可选 APC、AP、AC、PC其中 APC 表示年龄-时期-队列全模型这也是 BAPC 名字的由来。agegroup 参数用来指定数据里年龄组的跨度内置数据集是 5 年一组所以写 5year。predictTRUE 表示在拟合的同时对未来年份做预测预测的年限由另一个参数 default 控制默认会外推同等长度的时间段。typeinla 是指定使用 INLA 作为计算引擎早期版本还有 typemcmc 的选项但现在已经基本废弃除非你安装了额外的 MCMC 依赖包。这一步跑通后再把自己的数据组织成和 ICD 一样的格式。ICD 数据本质上就是一个矩阵行是年龄组列是年份每个单元格是发病人数或死亡率。你前面做的 wide 矩阵只要把行名修改成年份字符串列名修改成年龄组字符串就可以直接传给 BAPC()。注意如果数据里存在零值INLA 的泊松或二项模型可能给出极小但非零的后验估计值不需要过度清洗。4.2 直接调用 inla() 手动构建的姿势当 BAPC 包不够灵活时BAPC 包封装的函数虽然方便但有时候你想自定义先验分布、改变预测年限或者调整随机效应模型这时候就得绕开 BAPC手动构建 INLA 公式。这个做法更接近底层也更适合需要调整模型结构的人。library(INLA) # 构建数据帧每一行对应一个年龄-年份格 data_df - expand.grid(age 1:n_age, year 1:n_year, cohort as.vector(cohort_id)) data_df$y - as.vector(wide) # 响应变量 # 定义随机游走先验常用于年龄和队列效应 prior_age - list(prec list(prior pc.prec, param c(0.5, 0.01))) prior_period - list(prec list(prior rw2, param c(0.5, 0.01))) # 构建公式年龄用 rw2 随机效应时期用 rw1队列用 iid formula - y ~ f(age, model rw2, prior prior_age) f(year, model rw1, prior prior_period) f(cohort, model iid) # 拟合 INLA res - inla(formula, family poisson, data data_df, control.predictor list(compute TRUE), control.compute list(config TRUE, dic TRUE))这里的 f() 是 INLA 定义随机效应的函数。age 使用 rw2 表示二阶随机游走适合年龄效应这种平滑但有局部变化的曲线year 使用 rw1 是一阶随机游走对时期效应的波动有更好的适应性cohort 使用 iid 意思是队列之间完全独立不做平滑。为什么要这样搭配因为年龄效应通常随年龄缓慢上升再下降需要二阶惩罚让曲线平滑时期效应可能受政策、治疗水平等突发因素影响一阶游走更敏感队列效应在很多疾病中表现为相邻队列相似但很难确定具体函数形式所以用独立的随机效应最保守。这段代码里 control.predictor 的 computeTRUE 很关键它让 INLA 返回每个数据点的后验预测值否则后续画置信区间时拿不到数据。control.compute 里 configTRUE 是为了后面做样本回放或进一步的后验推断dicTRUE 则会计算 DIC 值用于比较不同模型设定。运行到最后你可以用 summary(res) 查看各效应的后验均值和 95% 区间这是我每次判断模型收敛的第一入口。4.3 预测未来年份prediction 数据帧的构造要点BAPC 内置的预测功能很方便但如果你想自己控制预测年限就得手动构造一个新的数据帧把 age、year、cohort 三个索引扩展到未来年份。这一步的难点在于 cohort 编号的衔接新的预测年份减去年龄组起始年龄得到的是新的队列编号这个编号必须按照之前的 cohort_map 进行映射如果出现数据里没有的队列需要追加新的编号。# 假设预测 15 年 pred_years - (year_end 1):(year_end 15) n_pred_year - length(pred_years) # 构造预测数据帧的所有组合 pred_grid - expand.grid(age 1:n_age, year (n_year 1):(n_year n_pred_year)) # 计算预测格的队列原始编号年份实际年 - 年龄起始年 pred_cohort_raw - pred_grid$year # 这里不能用 pred_grid$year 直接做队列计算因为 year 是编号不是实际年份这一行注释说明了一个极易踩的坑如果之前把年份转换成了 1 到 n_year 的编号那么在预测时也必须用实际年份先算出队列原始值再映射到编号。很多人的数据在进入模型前就已经把年份改成了序号导致预测时的队列编号全部错位。我的建议是始终保留一个独立的年份实际值向量比如 year_actual 1990:2019所有预测构造都在实际值上进行最后再映射到编号。否则预测结果看起来平滑实际上队列效应整体平移了若干年外推的曲线会系统性偏移。5. 结果解读与出图避坑从模型输出到发表级图表5.1 年龄、时期、队列效应图怎么看不要只看预测线BAPC 的最终输出包括年龄效应、时期效应、队列效应三个后验分布以及未来的预测发病/死亡数。很多人拿到结果图只看最后预测的曲线打了个箭头这不对。三个效应揭示了完全不同的信息年龄效应代表疾病随年龄变化的固有风险时期效应代表所有年龄组同时受到的影响比如筛查推广、诊断标准改变队列效应代表不同出生年份的人由于成长环境差异带来的风险差异。这里要特别提醒一个容易误读的点时期效应和队列效应存在内在共线性。年龄-时期-队列三者的线性趋势无法被唯一识别INLA 通过随机效应的先验施加了额外的约束所以后验估计有缩水。因此你在解读时期效应的下降趋势时不能说是绝对的下降而应该描述为在模型约束下的相对变化。5.2 预测区间的含义宽区间不代表模型差INLA 给出的预测带是后验预测区间它包含了参数的不确定性和未来数据本身的随机性。有时候你看到预测区间在远年部分变得很宽甚至上下限翻倍这不是模型报错而是外推时间越远队列和时期效应的不确定性累积越明显。发表论文时不能只画一个点预测线必须把 95% 区间带上。我在实际处理时常用 AgePeriodCohort 的绘图函数但它的默认配色和坐标轴刻度不一定合适。下面的代码演示了怎么提取预测结果并手动绘制一条带区间的曲线library(ggplot2) # 从 INLA 结果中提取每年的总体预测均值 pred_mean - res$summary.fitted.values[data_df$year n_year, mean] pred_ll - res$summary.fitted.values[data_df$year n_year, 0.025quant] pred_ul - res$summary.fitted.values[data_df$year n_year, 0.975quant] # 预测年份的年龄标准化或者总数需要按年龄汇总 # 这里假设每个格子的预测值已经求和归一化 plot_df - data.frame(year rep(pred_years, each n_age), age rep(1:n_age, times n_pred_year), mean pred_mean, ll pred_ll, ul pred_ul) # 如果希望按年龄加权合并可先按年聚合 agg_df - aggregate(cbind(mean, ll, ul) ~ year, data plot_df, FUN sum) ggplot(agg_df, aes(x year, y mean)) geom_ribbon(aes(ymin ll, ymax ul), alpha 0.3) geom_line(size 1) labs(x Year, y Predicted count)这段代码里最关键的是按 year 聚合前必须确认你自己的研究希望按年龄分别展示还是合并展示。如果研究关注特定年龄组比如 65 岁以上人群的疾病负担就不能直接聚合而是应该 filter(age 13) 后再画这个年龄组的预测序列。我见过有人把青年和老年混在一起预测总数看着平稳实际上年龄内部结构已经完全变形。5.3 避坑章BAPC 分析中最常见的五个具体问题问题一安装 BAPC 时提示 R 版本过低现象安装 BAPC 包时报错 package ‘BAPC’ is not available for this version of R。 原因这台机器上的 R 版本太旧BAPC 的当前版本要求 R 4.1 以上。 解决先升级 R 到最新稳定版再重新安装。如果因为项目原因不能升 R可以尝试安装 BAPC 的历史版本在 CRAN 的 archive 目录里找旧版压缩包用 install.packages(BAPC_1.0.0.tar.gz, repos NULL) 的方式安装但旧版本可能与新 INLA 不兼容。问题二inla() 报错 Error in inla.model.construct.rw2现象调用 BAPC() 或 inla() 时出现 rw2 模型构造失败错误信息里有一长串矩阵维度。 原因多数情况下是因为 age 或者 cohort 因子变量里存在只有单一水平的组INLA 在构建随机游走精度矩阵时无法计算二阶差分的惩罚。 解决检查你传入 BAPC 的矩阵是否出现了全零行或全零列。比如某一年龄组在整个时间段内都是 0 发病数INLA 会认为该水平没有信息量但仍然要参与矩阵运算。处理办法是删除全零的年龄组或年份或者在数据准备阶段用很小的非零值比如 0.01替代纯零但这样做会影响结果的解释我更建议直接删除完全无数据的维度。问题三预测值出现负数现象预测曲线在远期出现负的发病数或死亡率明显不符合生物学常识。 原因默认情况下 INLA 使用高斯近似作为默认似然当你没有指定 family poisson 时BAPC 包装函数可能没有正确传递分布家族导致模型按正态分布处理了计数数据。 解决如果你发现预测值有负值回到代码检查是否显式设置了 family poisson 或者 family binomial。另外检查响应变量是否真的是整数计数而不是已经除以人口数得到的率。率数据应该用二项模型或者先进行 log 变换不能让模型把率当作正态分布来拟合。问题四绘图时置信区间突然收窄成一条线现象AgePeriodCohort 画出的预测区间在某个年份节点后突然从宽变窄几乎贴住点估计。 原因这是预测期和观测期衔接处处理方式的问题。如果你的预测期数据帧里 cohort 编号没有正确扩展INLA 会用已有的队列效应覆盖新队列导致预测的变化自由度降低。 解决检查 cohort_id 的最大值是否小于未来预测所需的队列编号最大值。比如观测期是 1990-2019 年年龄组 0-4 岁的队列编号范围是 1986-2019但预测到 2034 年时会出现新的队列编号 2030 左右如果没有把这些新编号加入模型绘图时区间就会被压缩。问题五R 会话崩溃或内存耗尽现象跑 inla() 到一半 RStudio 直接闪退或者提示 cannot allocate vector of size。 原因INLA 本身对内存的占用较高尤其是当你用了很细的年龄组1 岁一组且年份跨度长时精度矩阵的维度可能达到几千乘几千稀疏矩阵存储也要几百 MB。更常见的是同时打开了太多 R 进程或者浏览器占用了大量内存。 解决先减少数据维度比如把 1 岁年龄组合并成 5 岁一组或者把预测年限缩短到 10 年。另外在调用 inla() 前设置 inla.setOption(num.threads 2) 限制线程数防止 INLA 默认占用所有 CPU 核心导致系统假死。如果内存实在不够用 64 位版本的 R 并增加虚拟内存。6. 验证预测效果的一个实用技巧回带历史数据做轨道验证模型跑通不是终点关键要验证预测可信度。我习惯做一个简单的后验验证把观测数据截断到某一年之前比如实际数据到 2019 年我只用 2000-2014 年拟合模型预测 2015-2019 年的结果然后跟真实观测值对比。这样可以直观看到模型在历史阶段的预测误差有多大也能侧面判断你的年龄分组和模型设定是否合理。# 截断数据仅使用 2000-2014 年作为拟合期 train_wide - wide[, colnames(wide) %in% as.character(2000:2014)] # 注意保持年龄顺序不变 train_cohort_id - cohort_id[, 1:ncol(train_wide)] # 重新跑一遍 BAPC预测到 2019 年 bapc_val - BAPC(train_wide, model APC, predict TRUE, prediction list(year 2015:2019)) # 提取预测均值和观测真实值 val_pred - bapc_val$prediction$rate val_obs - wide[, colnames(wide) %in% as.character(2015:2019)] # 计算相对误差 rel_error - (val_pred - val_obs) / val_obs * 100这段验证代码里的 prediction 参数接受一个列表可以指定预测的年份范围。我一般会计算分年龄组的相对误差看是不是特定年龄组误差特别大。如果某个年龄组的误差超过 30%就要考虑数据本身在该年龄组是否存在报告不全或者模型对队列效应的假设在这个年龄组不适用。比如中老年组往往因为诊断率提高导致时期效应很强可能需要把时期效应的随机游走从 rw1 改成 rw2 来获得更平滑的趋势。另外还有一种更轻量的验证方式比较 BAPC 的预测结果和用简单负二项回归外推的结果。虽然没有概率上的严谨对比但可以看出贝叶斯模型的预测区间是否明显比朴素方法宽以及中期预测点估计是否落在朴素方法的置信区间内。这不会写进论文但能帮你判断模型有没有系统性偏倚。如果验证发现误差在可接受范围内接下来才用全量数据重新拟合并把预测年限设为研究需要的长度。经过这个流程你对自己模型输出的解释就更有底气了。从那以后我每次做 BAPC 分析都强制把回带验证这一步加在正式预测之前哪怕结果已经出来也返工重跑这个习惯帮我避免过至少一次因为年龄组顺序错乱导致的完全错误的趋势预测。希望这些经验能帮你在 GBD 数据的 BAPC 分析上少走弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/10/25 18:20:29 · 阅读 1.2k GBD数据库疾病负担预测：R语言BAPC模型完整流程与避坑指南
简介一份用于全球疾病负担GBD数据库贝叶斯年龄-时期-队列BAPC分析的R包资源主要面向公共卫生研究人员和数据爱好者帮助完成疾病发病率和死亡率的趋势分解与短期预测。压缩包内含41个文件大小约57KB其中13个R源代码文件承担核心统计函数13个Rd文档提供函数接口说明3个Rda数据集用于实例演示另有RMarkdown、Markdown等文件辅助理解使用流程整体体量小巧但结构清晰。资源以Nordpred包及其相关扩展代码为主体实现年龄、时期、队列效应估计、预测及可视化绘图并附带挪威男性癌症数据作为教学实例同时提供数据生成脚本与原始文本数据便于研究者迁移至GBD各类疾病负担分析场景。已有2089人学习适合需要快速上手BAPC方法、减少环境配置成本并希望结合示例数据理解建模思路的读者。1. 用 R 包跑 GBD 数据库的 BAPC 分析从下载到出图先把这四件事搞清楚这几年做疾病负担测算的人越来越多GBD 数据库的现成数据好拿但真正想从“描述发病率和死亡率”往前再走一步做未来趋势预测绕不开 BAPC 这个模型。BAPC 全称是 Bayesian Age-Period-Cohort analysis也就是贝叶斯年龄-时期-队列分析它能把年龄结构、时期趋势和出生队列效应拆开再外推未来几年的疾病负担。这个分析在 R 里不是一个大而全的官方包而是由若干个 R 包配合完成有的人卡在数据格式有的人卡在 INLA 安装还有的人跑出来结果图长得很怪但不知道怎么排查。我这次分享的是我自己整理过的 R 包清单和配套脚本覆盖从 GBD 数据提取、年龄分组、APC 模型拟合到 BAPC 预测出图的完整流程。适合两类人一类是刚接触 GBD 的医学生或公卫从业者另一类是把 Stata 或 SPSS 用熟、但想在 R 里跑贝叶斯预测的进阶用户。下面直接聊每个环节怎么落地参数怎么设以及我踩过的几个实实在在的坑。2. BAPC 分析到底依赖哪些 R 包选型理由与版本边界2.1 核心包拆解INLA、BAPC、AgePeriodCohort 各自的作用BAPC 在 R 里的实现并不像 lm() 那样一个函数搞定它实际是一条链条。链条最底层是 INLA负责贝叶斯计算的近似推断中间层是 BAPC 包把年龄-时期-队列模型写成 INLA 能吃的公式最上层是 AgePeriodCohort 包提供绘图和结果的整理函数。这三个包缺一个都不行而且版本之间还有隐藏的兼容性问题。INLA 是整条链的引擎。它不是 CRAN 上的普通包而是需要从 R-INLA 官方仓库安装因为它的计算内核依赖外部 C 库和系统环境。我最早在 R 4.0 上直接 install.packages(INLA) 当然失败后来才知道必须用 inlainstall 的仓库地址。BAPC 包本身依赖 INLA 的结果对象它不自己算 MCMC而是把数据整理成 INLA 需要的格式然后调用 inla()。AgePeriodCohort 则负责可视化和一些描述性统计它的预测区间作图函数是很多人最后出图时最依赖的。选型理由一句话如果只是做简单的年龄标化率趋势用不了 BAPC但如果要分解年龄-时期-队列效应并外推预测BAPC 是目前 R 社区里相对成熟、文档多、可复现的方案。相比之下手动写 JAGS 或 Stan 模型的灵活度更高但你要自己处理年龄分组矩阵、队列共线性、预测期设计矩阵这些细节对大多数人来说门槛偏高。2.2 安装顺序与系统环境先装底层再装上层我一般建议按 INLA → BAPC → AgePeriodCohort 的顺序装先解决底层依赖再往上装。INLA 的安装用下面这段代码注意不要用 install.packages而是用官方仓库的 URL# 安装 INLA选择测试版或稳定版一般用 stable 即可 install.packages(INLA, reposc(getOption(repos), INLAhttps://inla.r-inla-download.org/R/stable), depTRUE) # 安装 BAPC 和 AgePeriodCohort这两个在 CRAN 上 install.packages(BAPC) install.packages(AgePeriodCohort)这段代码里repos 参数把 INLA 的镜像源追加到当前源列表这样 install.packages 才会从 R-INLA 官方源拉包。depTRUE 意思是一并安装依赖包否则可能会漏掉 Matrix、sp 等底层包。安装完后别急着直接用先跑一句 packageVersion(INLA) 确认版本号常见能正常用的是 22.x 以后版本太老的版本在 R 4.2 以上容易出现矩阵维度不匹配的报错。装完三个包后还有一个容易忽略的环节INLA 在运行时会调用系统级 BLAS/LAPACK 库如果你的 R 是 Windows 下 32 位和 64 位共存的环境建议统一用 64 位 R 并且关闭多线程冲突。我见过有人跑 inla() 的时候 RStudio 直接死掉最后发现是系统 PATH 里同时存在多个 Rtools 版本导致编译环境混乱解决办法是只保留一个 Rtools 并重启 R 会话。3. 把 GBD 数据库数据变成 BAPC 能吃的格式年龄分组与矩阵构造3.1 GBD 数据下载后的结构特征先看懂三个维度GBD 数据库导出的数据通常是一个长格式表格每一行代表某个国家、某个年份、某个年龄组、某个性别、某种疾病的患病率或死亡率。以 cause of death 或 incidence 为例列里常见的有 location_name、year、age_name、sex、val、upper、lower 这些字段。BAPC 分析前你必须把长表拆成三维矩阵年龄组作为行年份作为列队列作为对角线方向的衍生结构。这里有一个关键点GBD 的年龄组并不是均匀分布的。GBD 2019 里年龄组通常按 5 岁一组划分从 0-4 岁到 95 岁以上但有的疾病负担研究会用 0-14 岁、15-49 岁、50 岁以上这种粗略分组。BAPC 的年龄-时期-队列模型对年龄组的宽度和一致性很敏感如果你在数据准备阶段混用了不同宽度的年龄组预测结果会出现明显的锯齿。我的建议是初始整理尽量不动 GBD 原始的年龄组保持 5 岁一组再用 merge 或 pivot 操作把性别合并或分层。3.2 数据清洗与宽表转换pivot_wider 与年龄按序排列下面这段代码是我常用的清洗流程假设你已经从 GBD 官网下载了 CSV读进来之后先过滤出你要的 location、sex 和 measure然后转成宽表library(dplyr) library(tidyr) # 假设 dat 是 GBD 导出数据列名按照下载时的标准列名 dat - read.csv(gbd_incidence.csv, stringsAsFactors FALSE) # 筛选只保留中国、全部年龄、两性合计、发病数 dat_sub - dat %% filter(location_name China, sex_name Both, measure_name Incidence, metric_name Number) %% select(year, age_name, val) %% mutate(age_name factor(age_name, levels c(0-4 years, 5-9 years, 10-14 years, 15-19 years, 20-24 years, 25-29 years, 30-34 years, 35-39 years, 40-44 years, 45-49 years, 50-54 years, 55-59 years, 60-64 years, 65-69 years, 70-74 years, 75-79 years, 80-84 years, 85-89 years, 90-94 years, 95 years))) # 转成宽表每行是一个年龄组每列是一个年份 wide - dat_sub %% pivot_wider(names_from year, values_from val) %% arrange(age_name) %% select(-age_name) %% as.matrix() rownames(wide) - levels(dat_sub$age_name)这段代码里有两个地方容易翻车。第一是 age_name 的 levels 必须按时间顺序排列R 的 factor 默认按字母排序如果你不显式指定 levels转成矩阵后 10-14 岁会排在 0-4 岁前面整个模型的年龄效应曲线就会对不上号。第二是 pivot_wider 之后可能出现 NA 值原因是某些年份 GBD 没有对应年龄组的数据比如最早年份之前没有 95 岁以上记录这时候需要明确用 0 填充还是插补我一般倾向用前向填充或干脆删掉缺失对应的年份列避免 INLA 在拟合阶段因 NA 报错。3.3 构造年龄-时期-队列索引矩阵别用手工循环BAPC 包内部需要三个索引向量年龄组的序号、年份的序号、出生队列的序号。出生队列的计算规则是 period(年份) 减去 age(年龄组起点)再映射到一个连续的整数编号。这里最容易出错的是年龄组用区间的起始年龄还是中点年龄。BAPC 的默认逻辑是使用年龄组的编号本身作为年龄索引但队列编号需要你预先算好。# 宽表 wide 已经准备好行数为年龄组数列数为年份数 n_age - nrow(wide) n_year - ncol(wide) year_start - as.numeric(colnames(wide)[1]) year_end - as.numeric(colnames(wide)[ncol(wide)]) # 每个年龄组的起始年龄这里假设是 5 岁一组第一组 0-4 岁 age_start - seq(0, by 5, length.out n_age) # 构造每个格的队列编号队列 年份 - 年龄组起始值 cohort_index - matrix(NA, nrow n_age, ncol n_year) for (i in 1:n_age) { for (j in 1:n_year) { cohort_index[i, j] - as.numeric(colnames(wide)[j]) - age_start[i] } } # 队列编号需要从 1 开始连续编号不能有负数或跳号 cohort_unique - sort(unique(as.vector(cohort_index))) cohort_map - setNames(seq_along(cohort_unique), cohort_unique) cohort_id - apply(cohort_index, 2, function(x) cohort_map[as.character(x)])这段代码本质是把二维的年龄-年份格子影射成一个一维的队列向量。为什么要手动改编号因为 BAPC 包在拟合时会把队列效应当作一个随机效应处理如果队列编号是 1950、1951、1960 这种不连续的原始年份INLA 会默认当成连续型协变量而不是分组效应导致结果完全错误。映射成 1, 2, 3... 之后INLA 才能正确识别这是 46 个独立的队列水平。循环在这里并不慢因为年龄组一般不超过 20 个年份最多三四十年但如果年份跨度到 100 年可以改用 outer 函数加速。4. 跑通 BAPC 模型从 INLA 公式到未来预测4.1 使用 BAPC 包的推荐流程先读官方示例再改数据BAPC 包自带一个示例数据集是荷兰的癌症发病数据。我每次拿到新数据第一步不是直接替换而是跑一遍它的示例脚本确认当前 R 环境里的 INLA、BAPC、AgePeriodCohort 三个包都能协同工作。因为这三个包的版本更新节奏不一致有时候 BAPC 更新后调用 INLA 的对象结构方式变了你不跑示例根本看不出来。library(BAPC) data(ICD) # 这个内置数据集是 ICD 编码的荷兰癌症数据包含了 age, period, cohort 三个维度 # 先跑一个默认模型验证环境没问题 bapc_model - BAPC(ICD, model APC, predict TRUE, agegroup 5year, type inla)BAPC() 函数里最核心的参数是 model可选 APC、AP、AC、PC其中 APC 表示年龄-时期-队列全模型这也是 BAPC 名字的由来。agegroup 参数用来指定数据里年龄组的跨度内置数据集是 5 年一组所以写 5year。predictTRUE 表示在拟合的同时对未来年份做预测预测的年限由另一个参数 default 控制默认会外推同等长度的时间段。typeinla 是指定使用 INLA 作为计算引擎早期版本还有 typemcmc 的选项但现在已经基本废弃除非你安装了额外的 MCMC 依赖包。这一步跑通后再把自己的数据组织成和 ICD 一样的格式。ICD 数据本质上就是一个矩阵行是年龄组列是年份每个单元格是发病人数或死亡率。你前面做的 wide 矩阵只要把行名修改成年份字符串列名修改成年龄组字符串就可以直接传给 BAPC()。注意如果数据里存在零值INLA 的泊松或二项模型可能给出极小但非零的后验估计值不需要过度清洗。4.2 直接调用 inla() 手动构建的姿势当 BAPC 包不够灵活时BAPC 包封装的函数虽然方便但有时候你想自定义先验分布、改变预测年限或者调整随机效应模型这时候就得绕开 BAPC手动构建 INLA 公式。这个做法更接近底层也更适合需要调整模型结构的人。library(INLA) # 构建数据帧每一行对应一个年龄-年份格 data_df - expand.grid(age 1:n_age, year 1:n_year, cohort as.vector(cohort_id)) data_df$y - as.vector(wide) # 响应变量 # 定义随机游走先验常用于年龄和队列效应 prior_age - list(prec list(prior pc.prec, param c(0.5, 0.01))) prior_period - list(prec list(prior rw2, param c(0.5, 0.01))) # 构建公式年龄用 rw2 随机效应时期用 rw1队列用 iid formula - y ~ f(age, model rw2, prior prior_age) f(year, model rw1, prior prior_period) f(cohort, model iid) # 拟合 INLA res - inla(formula, family poisson, data data_df, control.predictor list(compute TRUE), control.compute list(config TRUE, dic TRUE))这里的 f() 是 INLA 定义随机效应的函数。age 使用 rw2 表示二阶随机游走适合年龄效应这种平滑但有局部变化的曲线year 使用 rw1 是一阶随机游走对时期效应的波动有更好的适应性cohort 使用 iid 意思是队列之间完全独立不做平滑。为什么要这样搭配因为年龄效应通常随年龄缓慢上升再下降需要二阶惩罚让曲线平滑时期效应可能受政策、治疗水平等突发因素影响一阶游走更敏感队列效应在很多疾病中表现为相邻队列相似但很难确定具体函数形式所以用独立的随机效应最保守。这段代码里 control.predictor 的 computeTRUE 很关键它让 INLA 返回每个数据点的后验预测值否则后续画置信区间时拿不到数据。control.compute 里 configTRUE 是为了后面做样本回放或进一步的后验推断dicTRUE 则会计算 DIC 值用于比较不同模型设定。运行到最后你可以用 summary(res) 查看各效应的后验均值和 95% 区间这是我每次判断模型收敛的第一入口。4.3 预测未来年份prediction 数据帧的构造要点BAPC 内置的预测功能很方便但如果你想自己控制预测年限就得手动构造一个新的数据帧把 age、year、cohort 三个索引扩展到未来年份。这一步的难点在于 cohort 编号的衔接新的预测年份减去年龄组起始年龄得到的是新的队列编号这个编号必须按照之前的 cohort_map 进行映射如果出现数据里没有的队列需要追加新的编号。# 假设预测 15 年 pred_years - (year_end 1):(year_end 15) n_pred_year - length(pred_years) # 构造预测数据帧的所有组合 pred_grid - expand.grid(age 1:n_age, year (n_year 1):(n_year n_pred_year)) # 计算预测格的队列原始编号年份实际年 - 年龄起始年 pred_cohort_raw - pred_grid$year # 这里不能用 pred_grid$year 直接做队列计算因为 year 是编号不是实际年份这一行注释说明了一个极易踩的坑如果之前把年份转换成了 1 到 n_year 的编号那么在预测时也必须用实际年份先算出队列原始值再映射到编号。很多人的数据在进入模型前就已经把年份改成了序号导致预测时的队列编号全部错位。我的建议是始终保留一个独立的年份实际值向量比如 year_actual 1990:2019所有预测构造都在实际值上进行最后再映射到编号。否则预测结果看起来平滑实际上队列效应整体平移了若干年外推的曲线会系统性偏移。5. 结果解读与出图避坑从模型输出到发表级图表5.1 年龄、时期、队列效应图怎么看不要只看预测线BAPC 的最终输出包括年龄效应、时期效应、队列效应三个后验分布以及未来的预测发病/死亡数。很多人拿到结果图只看最后预测的曲线打了个箭头这不对。三个效应揭示了完全不同的信息年龄效应代表疾病随年龄变化的固有风险时期效应代表所有年龄组同时受到的影响比如筛查推广、诊断标准改变队列效应代表不同出生年份的人由于成长环境差异带来的风险差异。这里要特别提醒一个容易误读的点时期效应和队列效应存在内在共线性。年龄-时期-队列三者的线性趋势无法被唯一识别INLA 通过随机效应的先验施加了额外的约束所以后验估计有缩水。因此你在解读时期效应的下降趋势时不能说是绝对的下降而应该描述为在模型约束下的相对变化。5.2 预测区间的含义宽区间不代表模型差INLA 给出的预测带是后验预测区间它包含了参数的不确定性和未来数据本身的随机性。有时候你看到预测区间在远年部分变得很宽甚至上下限翻倍这不是模型报错而是外推时间越远队列和时期效应的不确定性累积越明显。发表论文时不能只画一个点预测线必须把 95% 区间带上。我在实际处理时常用 AgePeriodCohort 的绘图函数但它的默认配色和坐标轴刻度不一定合适。下面的代码演示了怎么提取预测结果并手动绘制一条带区间的曲线library(ggplot2) # 从 INLA 结果中提取每年的总体预测均值 pred_mean - res$summary.fitted.values[data_df$year n_year, mean] pred_ll - res$summary.fitted.values[data_df$year n_year, 0.025quant] pred_ul - res$summary.fitted.values[data_df$year n_year, 0.975quant] # 预测年份的年龄标准化或者总数需要按年龄汇总 # 这里假设每个格子的预测值已经求和归一化 plot_df - data.frame(year rep(pred_years, each n_age), age rep(1:n_age, times n_pred_year), mean pred_mean, ll pred_ll, ul pred_ul) # 如果希望按年龄加权合并可先按年聚合 agg_df - aggregate(cbind(mean, ll, ul) ~ year, data plot_df, FUN sum) ggplot(agg_df, aes(x year, y mean)) geom_ribbon(aes(ymin ll, ymax ul), alpha 0.3) geom_line(size 1) labs(x Year, y Predicted count)这段代码里最关键的是按 year 聚合前必须确认你自己的研究希望按年龄分别展示还是合并展示。如果研究关注特定年龄组比如 65 岁以上人群的疾病负担就不能直接聚合而是应该 filter(age 13) 后再画这个年龄组的预测序列。我见过有人把青年和老年混在一起预测总数看着平稳实际上年龄内部结构已经完全变形。5.3 避坑章BAPC 分析中最常见的五个具体问题问题一安装 BAPC 时提示 R 版本过低现象安装 BAPC 包时报错 package ‘BAPC’ is not available for this version of R。 原因这台机器上的 R 版本太旧BAPC 的当前版本要求 R 4.1 以上。 解决先升级 R 到最新稳定版再重新安装。如果因为项目原因不能升 R可以尝试安装 BAPC 的历史版本在 CRAN 的 archive 目录里找旧版压缩包用 install.packages(BAPC_1.0.0.tar.gz, repos NULL) 的方式安装但旧版本可能与新 INLA 不兼容。问题二inla() 报错 Error in inla.model.construct.rw2现象调用 BAPC() 或 inla() 时出现 rw2 模型构造失败错误信息里有一长串矩阵维度。 原因多数情况下是因为 age 或者 cohort 因子变量里存在只有单一水平的组INLA 在构建随机游走精度矩阵时无法计算二阶差分的惩罚。 解决检查你传入 BAPC 的矩阵是否出现了全零行或全零列。比如某一年龄组在整个时间段内都是 0 发病数INLA 会认为该水平没有信息量但仍然要参与矩阵运算。处理办法是删除全零的年龄组或年份或者在数据准备阶段用很小的非零值比如 0.01替代纯零但这样做会影响结果的解释我更建议直接删除完全无数据的维度。问题三预测值出现负数现象预测曲线在远期出现负的发病数或死亡率明显不符合生物学常识。 原因默认情况下 INLA 使用高斯近似作为默认似然当你没有指定 family poisson 时BAPC 包装函数可能没有正确传递分布家族导致模型按正态分布处理了计数数据。 解决如果你发现预测值有负值回到代码检查是否显式设置了 family poisson 或者 family binomial。另外检查响应变量是否真的是整数计数而不是已经除以人口数得到的率。率数据应该用二项模型或者先进行 log 变换不能让模型把率当作正态分布来拟合。问题四绘图时置信区间突然收窄成一条线现象AgePeriodCohort 画出的预测区间在某个年份节点后突然从宽变窄几乎贴住点估计。 原因这是预测期和观测期衔接处处理方式的问题。如果你的预测期数据帧里 cohort 编号没有正确扩展INLA 会用已有的队列效应覆盖新队列导致预测的变化自由度降低。 解决检查 cohort_id 的最大值是否小于未来预测所需的队列编号最大值。比如观测期是 1990-2019 年年龄组 0-4 岁的队列编号范围是 1986-2019但预测到 2034 年时会出现新的队列编号 2030 左右如果没有把这些新编号加入模型绘图时区间就会被压缩。问题五R 会话崩溃或内存耗尽现象跑 inla() 到一半 RStudio 直接闪退或者提示 cannot allocate vector of size。 原因INLA 本身对内存的占用较高尤其是当你用了很细的年龄组1 岁一组且年份跨度长时精度矩阵的维度可能达到几千乘几千稀疏矩阵存储也要几百 MB。更常见的是同时打开了太多 R 进程或者浏览器占用了大量内存。 解决先减少数据维度比如把 1 岁年龄组合并成 5 岁一组或者把预测年限缩短到 10 年。另外在调用 inla() 前设置 inla.setOption(num.threads 2) 限制线程数防止 INLA 默认占用所有 CPU 核心导致系统假死。如果内存实在不够用 64 位版本的 R 并增加虚拟内存。6. 验证预测效果的一个实用技巧回带历史数据做轨道验证模型跑通不是终点关键要验证预测可信度。我习惯做一个简单的后验验证把观测数据截断到某一年之前比如实际数据到 2019 年我只用 2000-2014 年拟合模型预测 2015-2019 年的结果然后跟真实观测值对比。这样可以直观看到模型在历史阶段的预测误差有多大也能侧面判断你的年龄分组和模型设定是否合理。# 截断数据仅使用 2000-2014 年作为拟合期 train_wide - wide[, colnames(wide) %in% as.character(2000:2014)] # 注意保持年龄顺序不变 train_cohort_id - cohort_id[, 1:ncol(train_wide)] # 重新跑一遍 BAPC预测到 2019 年 bapc_val - BAPC(train_wide, model APC, predict TRUE, prediction list(year 2015:2019)) # 提取预测均值和观测真实值 val_pred - bapc_val$prediction$rate val_obs - wide[, colnames(wide) %in% as.character(2015:2019)] # 计算相对误差 rel_error - (val_pred - val_obs) / val_obs * 100这段验证代码里的 prediction 参数接受一个列表可以指定预测的年份范围。我一般会计算分年龄组的相对误差看是不是特定年龄组误差特别大。如果某个年龄组的误差超过 30%就要考虑数据本身在该年龄组是否存在报告不全或者模型对队列效应的假设在这个年龄组不适用。比如中老年组往往因为诊断率提高导致时期效应很强可能需要把时期效应的随机游走从 rw1 改成 rw2 来获得更平滑的趋势。另外还有一种更轻量的验证方式比较 BAPC 的预测结果和用简单负二项回归外推的结果。虽然没有概率上的严谨对比但可以看出贝叶斯模型的预测区间是否明显比朴素方法宽以及中期预测点估计是否落在朴素方法的置信区间内。这不会写进论文但能帮你判断模型有没有系统性偏倚。如果验证发现误差在可接受范围内接下来才用全量数据重新拟合并把预测年限设为研究需要的长度。经过这个流程你对自己模型输出的解释就更有底气了。从那以后我每次做 BAPC 分析都强制把回带验证这一步加在正式预测之前哪怕结果已经出来也返工重跑这个习惯帮我避免过至少一次因为年龄组顺序错乱导致的完全错误的趋势预测。希望这些经验能帮你在 GBD 数据的 BAPC 分析上少走弯路。本文还有配套的精品资源点击获取