簡(jiǎn)介一份用于全球疾病負(fù)擔(dān)GBD數(shù)據(jù)庫貝葉斯年齡-時(shí)期-隊(duì)列BAPC分析的R包資源主要面向公共衛(wèi)生研究人員和數(shù)據(jù)愛好者幫助完成疾病發(fā)病率和死亡率的趨勢(shì)分解與短期預(yù)測(cè)。壓縮包內(nèi)含41個(gè)文件大小約57KB其中13個(gè)R源代碼文件承擔(dān)核心統(tǒng)計(jì)函數(shù)13個(gè)Rd文檔提供函數(shù)接口說明3個(gè)Rda數(shù)據(jù)集用于實(shí)例演示另有RMarkdown、Markdown等文件輔助理解使用流程整體體量小巧但結(jié)構(gòu)清晰。資源以Nordpred包及其相關(guān)擴(kuò)展代碼為主體實(shí)現(xiàn)年齡、時(shí)期、隊(duì)列效應(yīng)估計(jì)、預(yù)測(cè)及可視化繪圖并附帶挪威男性癌癥數(shù)據(jù)作為教學(xué)實(shí)例同時(shí)提供數(shù)據(jù)生成腳本與原始文本數(shù)據(jù)便于研究者遷移至GBD各類疾病負(fù)擔(dān)分析場(chǎng)景。已有2089人學(xué)習(xí)適合需要快速上手BAPC方法、減少環(huán)境配置成本并希望結(jié)合示例數(shù)據(jù)理解建模思路的讀者。1. 用 R 包跑 GBD 數(shù)據(jù)庫的 BAPC 分析從下載到出圖先把這四件事搞清楚這幾年做疾病負(fù)擔(dān)測(cè)算的人越來越多GBD 數(shù)據(jù)庫的現(xiàn)成數(shù)據(jù)好拿但真正想從“描述發(fā)病率和死亡率”往前再走一步做未來趨勢(shì)預(yù)測(cè)繞不開 BAPC 這個(gè)模型。BAPC 全稱是 Bayesian Age-Period-Cohort analysis也就是貝葉斯年齡-時(shí)期-隊(duì)列分析它能把年齡結(jié)構(gòu)、時(shí)期趨勢(shì)和出生隊(duì)列效應(yīng)拆開再外推未來幾年的疾病負(fù)擔(dān)。這個(gè)分析在 R 里不是一個(gè)大而全的官方包而是由若干個(gè) R 包配合完成有的人卡在數(shù)據(jù)格式有的人卡在 INLA 安裝還有的人跑出來結(jié)果圖長(zhǎng)得很怪但不知道怎么排查。我這次分享的是我自己整理過的 R 包清單和配套腳本覆蓋從 GBD 數(shù)據(jù)提取、年齡分組、APC 模型擬合到 BAPC 預(yù)測(cè)出圖的完整流程。適合兩類人一類是剛接觸 GBD 的醫(yī)學(xué)生或公衛(wèi)從業(yè)者另一類是把 Stata 或 SPSS 用熟、但想在 R 里跑貝葉斯預(yù)測(cè)的進(jìn)階用戶。下面直接聊每個(gè)環(huán)節(jié)怎么落地參數(shù)怎么設(shè)以及我踩過的幾個(gè)實(shí)實(shí)在在的坑。2. BAPC 分析到底依賴哪些 R 包選型理由與版本邊界2.1 核心包拆解INLA、BAPC、AgePeriodCohort 各自的作用BAPC 在 R 里的實(shí)現(xiàn)并不像 lm() 那樣一個(gè)函數(shù)搞定它實(shí)際是一條鏈條。鏈條最底層是 INLA負(fù)責(zé)貝葉斯計(jì)算的近似推斷中間層是 BAPC 包把年齡-時(shí)期-隊(duì)列模型寫成 INLA 能吃的公式最上層是 AgePeriodCohort 包提供繪圖和結(jié)果的整理函數(shù)。這三個(gè)包缺一個(gè)都不行而且版本之間還有隱藏的兼容性問題。INLA 是整條鏈的引擎。它不是 CRAN 上的普通包而是需要從 R-INLA 官方倉(cāng)庫安裝因?yàn)樗挠?jì)算內(nèi)核依賴外部 C 庫和系統(tǒng)環(huán)境。我最早在 R 4.0 上直接 install.packages(INLA) 當(dāng)然失敗后來才知道必須用 inlainstall 的倉(cāng)庫地址。BAPC 包本身依賴 INLA 的結(jié)果對(duì)象它不自己算 MCMC而是把數(shù)據(jù)整理成 INLA 需要的格式然后調(diào)用 inla()。AgePeriodCohort 則負(fù)責(zé)可視化和一些描述性統(tǒng)計(jì)它的預(yù)測(cè)區(qū)間作圖函數(shù)是很多人最后出圖時(shí)最依賴的。選型理由一句話如果只是做簡(jiǎn)單的年齡標(biāo)化率趨勢(shì)用不了 BAPC但如果要分解年齡-時(shí)期-隊(duì)列效應(yīng)并外推預(yù)測(cè)BAPC 是目前 R 社區(qū)里相對(duì)成熟、文檔多、可復(fù)現(xiàn)的方案。相比之下手動(dòng)寫 JAGS 或 Stan 模型的靈活度更高但你要自己處理年齡分組矩陣、隊(duì)列共線性、預(yù)測(cè)期設(shè)計(jì)矩陣這些細(xì)節(jié)對(duì)大多數(shù)人來說門檻偏高。2.2 安裝順序與系統(tǒng)環(huán)境先裝底層再裝上層我一般建議按 INLA → BAPC → AgePeriodCohort 的順序裝先解決底層依賴再往上裝。INLA 的安裝用下面這段代碼注意不要用 install.packages而是用官方倉(cāng)庫的 URL# 安裝 INLA選擇測(cè)試版或穩(wěn)定版一般用 stable 即可 install.packages(INLA, reposc(getOption(repos), INLAhttps://inla.r-inla-download.org/R/stable), depTRUE) # 安裝 BAPC 和 AgePeriodCohort這兩個(gè)在 CRAN 上 install.packages(BAPC) install.packages(AgePeriodCohort)這段代碼里repos 參數(shù)把 INLA 的鏡像源追加到當(dāng)前源列表這樣 install.packages 才會(huì)從 R-INLA 官方源拉包。depTRUE 意思是一并安裝依賴包否則可能會(huì)漏掉 Matrix、sp 等底層包。安裝完后別急著直接用先跑一句 packageVersion(INLA) 確認(rèn)版本號(hào)常見能正常用的是 22.x 以后版本太老的版本在 R 4.2 以上容易出現(xiàn)矩陣維度不匹配的報(bào)錯(cuò)。裝完三個(gè)包后還有一個(gè)容易忽略的環(huán)節(jié)INLA 在運(yùn)行時(shí)會(huì)調(diào)用系統(tǒng)級(jí) BLAS/LAPACK 庫如果你的 R 是 Windows 下 32 位和 64 位共存的環(huán)境建議統(tǒng)一用 64 位 R 并且關(guān)閉多線程沖突。我見過有人跑 inla() 的時(shí)候 RStudio 直接死掉最后發(fā)現(xiàn)是系統(tǒng) PATH 里同時(shí)存在多個(gè) Rtools 版本導(dǎo)致編譯環(huán)境混亂解決辦法是只保留一個(gè) Rtools 并重啟 R 會(huì)話。3. 把 GBD 數(shù)據(jù)庫數(shù)據(jù)變成 BAPC 能吃的格式年齡分組與矩陣構(gòu)造3.1 GBD 數(shù)據(jù)下載后的結(jié)構(gòu)特征先看懂三個(gè)維度GBD 數(shù)據(jù)庫導(dǎo)出的數(shù)據(jù)通常是一個(gè)長(zhǎng)格式表格每一行代表某個(gè)國(guó)家、某個(gè)年份、某個(gè)年齡組、某個(gè)性別、某種疾病的患病率或死亡率。以 cause of death 或 incidence 為例列里常見的有 location_name、year、age_name、sex、val、upper、lower 這些字段。BAPC 分析前你必須把長(zhǎng)表拆成三維矩陣年齡組作為行年份作為列隊(duì)列作為對(duì)角線方向的衍生結(jié)構(gòu)。這里有一個(gè)關(guān)鍵點(diǎn)GBD 的年齡組并不是均勻分布的。GBD 2019 里年齡組通常按 5 歲一組劃分從 0-4 歲到 95 歲以上但有的疾病負(fù)擔(dān)研究會(huì)用 0-14 歲、15-49 歲、50 歲以上這種粗略分組。BAPC 的年齡-時(shí)期-隊(duì)列模型對(duì)年齡組的寬度和一致性很敏感如果你在數(shù)據(jù)準(zhǔn)備階段混用了不同寬度的年齡組預(yù)測(cè)結(jié)果會(huì)出現(xiàn)明顯的鋸齒。我的建議是初始整理盡量不動(dòng) GBD 原始的年齡組保持 5 歲一組再用 merge 或 pivot 操作把性別合并或分層。3.2 數(shù)據(jù)清洗與寬表轉(zhuǎn)換pivot_wider 與年齡按序排列下面這段代碼是我常用的清洗流程假設(shè)你已經(jīng)從 GBD 官網(wǎng)下載了 CSV讀進(jìn)來之后先過濾出你要的 location、sex 和 measure然后轉(zhuǎn)成寬表library(dplyr) library(tidyr) # 假設(shè) dat 是 GBD 導(dǎo)出數(shù)據(jù)列名按照下載時(shí)的標(biāo)準(zhǔn)列名 dat - read.csv(gbd_incidence.csv, stringsAsFactors FALSE) # 篩選只保留中國(guó)、全部年齡、兩性合計(jì)、發(fā)病數(shù) 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))) # 轉(zhuǎn)成寬表每行是一個(gè)年齡組每列是一個(gè)年份 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)這段代碼里有兩個(gè)地方容易翻車。第一是 age_name 的 levels 必須按時(shí)間順序排列R 的 factor 默認(rèn)按字母排序如果你不顯式指定 levels轉(zhuǎn)成矩陣后 10-14 歲會(huì)排在 0-4 歲前面整個(gè)模型的年齡效應(yīng)曲線就會(huì)對(duì)不上號(hào)。第二是 pivot_wider 之后可能出現(xiàn) NA 值原因是某些年份 GBD 沒有對(duì)應(yīng)年齡組的數(shù)據(jù)比如最早年份之前沒有 95 歲以上記錄這時(shí)候需要明確用 0 填充還是插補(bǔ)我一般傾向用前向填充或干脆刪掉缺失對(duì)應(yīng)的年份列避免 INLA 在擬合階段因 NA 報(bào)錯(cuò)。3.3 構(gòu)造年齡-時(shí)期-隊(duì)列索引矩陣別用手工循環(huán)BAPC 包內(nèi)部需要三個(gè)索引向量年齡組的序號(hào)、年份的序號(hào)、出生隊(duì)列的序號(hào)。出生隊(duì)列的計(jì)算規(guī)則是 period(年份) 減去 age(年齡組起點(diǎn))再映射到一個(gè)連續(xù)的整數(shù)編號(hào)。這里最容易出錯(cuò)的是年齡組用區(qū)間的起始年齡還是中點(diǎn)年齡。BAPC 的默認(rèn)邏輯是使用年齡組的編號(hào)本身作為年齡索引但隊(duì)列編號(hào)需要你預(yù)先算好。# 寬表 wide 已經(jīng)準(zhǔn)備好行數(shù)為年齡組數(shù)列數(shù)為年份數(shù) n_age - nrow(wide) n_year - ncol(wide) year_start - as.numeric(colnames(wide)[1]) year_end - as.numeric(colnames(wide)[ncol(wide)]) # 每個(gè)年齡組的起始年齡這里假設(shè)是 5 歲一組第一組 0-4 歲 age_start - seq(0, by 5, length.out n_age) # 構(gòu)造每個(gè)格的隊(duì)列編號(hào)隊(duì)列 年份 - 年齡組起始值 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] } } # 隊(duì)列編號(hào)需要從 1 開始連續(xù)編號(hào)不能有負(fù)數(shù)或跳號(hào) 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)])這段代碼本質(zhì)是把二維的年齡-年份格子影射成一個(gè)一維的隊(duì)列向量。為什么要手動(dòng)改編號(hào)因?yàn)?BAPC 包在擬合時(shí)會(huì)把隊(duì)列效應(yīng)當(dāng)作一個(gè)隨機(jī)效應(yīng)處理如果隊(duì)列編號(hào)是 1950、1951、1960 這種不連續(xù)的原始年份INLA 會(huì)默認(rèn)當(dāng)成連續(xù)型協(xié)變量而不是分組效應(yīng)導(dǎo)致結(jié)果完全錯(cuò)誤。映射成 1, 2, 3... 之后INLA 才能正確識(shí)別這是 46 個(gè)獨(dú)立的隊(duì)列水平。循環(huán)在這里并不慢因?yàn)槟挲g組一般不超過 20 個(gè)年份最多三四十年但如果年份跨度到 100 年可以改用 outer 函數(shù)加速。4. 跑通 BAPC 模型從 INLA 公式到未來預(yù)測(cè)4.1 使用 BAPC 包的推薦流程先讀官方示例再改數(shù)據(jù)BAPC 包自帶一個(gè)示例數(shù)據(jù)集是荷蘭的癌癥發(fā)病數(shù)據(jù)。我每次拿到新數(shù)據(jù)第一步不是直接替換而是跑一遍它的示例腳本確認(rèn)當(dāng)前 R 環(huán)境里的 INLA、BAPC、AgePeriodCohort 三個(gè)包都能協(xié)同工作。因?yàn)檫@三個(gè)包的版本更新節(jié)奏不一致有時(shí)候 BAPC 更新后調(diào)用 INLA 的對(duì)象結(jié)構(gòu)方式變了你不跑示例根本看不出來。library(BAPC) data(ICD) # 這個(gè)內(nèi)置數(shù)據(jù)集是 ICD 編碼的荷蘭癌癥數(shù)據(jù)包含了 age, period, cohort 三個(gè)維度 # 先跑一個(gè)默認(rèn)模型驗(yàn)證環(huán)境沒問題 bapc_model - BAPC(ICD, model APC, predict TRUE, agegroup 5year, type inla)BAPC() 函數(shù)里最核心的參數(shù)是 model可選 APC、AP、AC、PC其中 APC 表示年齡-時(shí)期-隊(duì)列全模型這也是 BAPC 名字的由來。agegroup 參數(shù)用來指定數(shù)據(jù)里年齡組的跨度內(nèi)置數(shù)據(jù)集是 5 年一組所以寫 5year。predictTRUE 表示在擬合的同時(shí)對(duì)未來年份做預(yù)測(cè)預(yù)測(cè)的年限由另一個(gè)參數(shù) default 控制默認(rèn)會(huì)外推同等長(zhǎng)度的時(shí)間段。typeinla 是指定使用 INLA 作為計(jì)算引擎早期版本還有 typemcmc 的選項(xiàng)但現(xiàn)在已經(jīng)基本廢棄除非你安裝了額外的 MCMC 依賴包。這一步跑通后再把自己的數(shù)據(jù)組織成和 ICD 一樣的格式。ICD 數(shù)據(jù)本質(zhì)上就是一個(gè)矩陣行是年齡組列是年份每個(gè)單元格是發(fā)病人數(shù)或死亡率。你前面做的 wide 矩陣只要把行名修改成年份字符串列名修改成年齡組字符串就可以直接傳給 BAPC()。注意如果數(shù)據(jù)里存在零值INLA 的泊松或二項(xiàng)模型可能給出極小但非零的后驗(yàn)估計(jì)值不需要過度清洗。4.2 直接調(diào)用 inla() 手動(dòng)構(gòu)建的姿勢(shì)當(dāng) BAPC 包不夠靈活時(shí)BAPC 包封裝的函數(shù)雖然方便但有時(shí)候你想自定義先驗(yàn)分布、改變預(yù)測(cè)年限或者調(diào)整隨機(jī)效應(yīng)模型這時(shí)候就得繞開 BAPC手動(dòng)構(gòu)建 INLA 公式。這個(gè)做法更接近底層也更適合需要調(diào)整模型結(jié)構(gòu)的人。library(INLA) # 構(gòu)建數(shù)據(jù)幀每一行對(duì)應(yīng)一個(gè)年齡-年份格 data_df - expand.grid(age 1:n_age, year 1:n_year, cohort as.vector(cohort_id)) data_df$y - as.vector(wide) # 響應(yīng)變量 # 定義隨機(jī)游走先驗(yàn)常用于年齡和隊(duì)列效應(yīng) 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))) # 構(gòu)建公式年齡用 rw2 隨機(jī)效應(yīng)時(shí)期用 rw1隊(duì)列用 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 定義隨機(jī)效應(yīng)的函數(shù)。age 使用 rw2 表示二階隨機(jī)游走適合年齡效應(yīng)這種平滑但有局部變化的曲線year 使用 rw1 是一階隨機(jī)游走對(duì)時(shí)期效應(yīng)的波動(dòng)有更好的適應(yīng)性cohort 使用 iid 意思是隊(duì)列之間完全獨(dú)立不做平滑。為什么要這樣搭配因?yàn)槟挲g效應(yīng)通常隨年齡緩慢上升再下降需要二階懲罰讓曲線平滑時(shí)期效應(yīng)可能受政策、治療水平等突發(fā)因素影響一階游走更敏感隊(duì)列效應(yīng)在很多疾病中表現(xiàn)為相鄰隊(duì)列相似但很難確定具體函數(shù)形式所以用獨(dú)立的隨機(jī)效應(yīng)最保守。這段代碼里 control.predictor 的 computeTRUE 很關(guān)鍵它讓 INLA 返回每個(gè)數(shù)據(jù)點(diǎn)的后驗(yàn)預(yù)測(cè)值否則后續(xù)畫置信區(qū)間時(shí)拿不到數(shù)據(jù)。control.compute 里 configTRUE 是為了后面做樣本回放或進(jìn)一步的后驗(yàn)推斷dicTRUE 則會(huì)計(jì)算 DIC 值用于比較不同模型設(shè)定。運(yùn)行到最后你可以用 summary(res) 查看各效應(yīng)的后驗(yàn)均值和 95% 區(qū)間這是我每次判斷模型收斂的第一入口。4.3 預(yù)測(cè)未來年份prediction 數(shù)據(jù)幀的構(gòu)造要點(diǎn)BAPC 內(nèi)置的預(yù)測(cè)功能很方便但如果你想自己控制預(yù)測(cè)年限就得手動(dòng)構(gòu)造一個(gè)新的數(shù)據(jù)幀把 age、year、cohort 三個(gè)索引擴(kuò)展到未來年份。這一步的難點(diǎn)在于 cohort 編號(hào)的銜接新的預(yù)測(cè)年份減去年齡組起始年齡得到的是新的隊(duì)列編號(hào)這個(gè)編號(hào)必須按照之前的 cohort_map 進(jìn)行映射如果出現(xiàn)數(shù)據(jù)里沒有的隊(duì)列需要追加新的編號(hào)。# 假設(shè)預(yù)測(cè) 15 年 pred_years - (year_end 1):(year_end 15) n_pred_year - length(pred_years) # 構(gòu)造預(yù)測(cè)數(shù)據(jù)幀的所有組合 pred_grid - expand.grid(age 1:n_age, year (n_year 1):(n_year n_pred_year)) # 計(jì)算預(yù)測(cè)格的隊(duì)列原始編號(hào)年份實(shí)際年 - 年齡起始年 pred_cohort_raw - pred_grid$year # 這里不能用 pred_grid$year 直接做隊(duì)列計(jì)算因?yàn)?year 是編號(hào)不是實(shí)際年份這一行注釋說明了一個(gè)極易踩的坑如果之前把年份轉(zhuǎn)換成了 1 到 n_year 的編號(hào)那么在預(yù)測(cè)時(shí)也必須用實(shí)際年份先算出隊(duì)列原始值再映射到編號(hào)。很多人的數(shù)據(jù)在進(jìn)入模型前就已經(jīng)把年份改成了序號(hào)導(dǎo)致預(yù)測(cè)時(shí)的隊(duì)列編號(hào)全部錯(cuò)位。我的建議是始終保留一個(gè)獨(dú)立的年份實(shí)際值向量比如 year_actual 1990:2019所有預(yù)測(cè)構(gòu)造都在實(shí)際值上進(jìn)行最后再映射到編號(hào)。否則預(yù)測(cè)結(jié)果看起來平滑實(shí)際上隊(duì)列效應(yīng)整體平移了若干年外推的曲線會(huì)系統(tǒng)性偏移。5. 結(jié)果解讀與出圖避坑從模型輸出到發(fā)表級(jí)圖表5.1 年齡、時(shí)期、隊(duì)列效應(yīng)圖怎么看不要只看預(yù)測(cè)線BAPC 的最終輸出包括年齡效應(yīng)、時(shí)期效應(yīng)、隊(duì)列效應(yīng)三個(gè)后驗(yàn)分布以及未來的預(yù)測(cè)發(fā)病/死亡數(shù)。很多人拿到結(jié)果圖只看最后預(yù)測(cè)的曲線打了個(gè)箭頭這不對(duì)。三個(gè)效應(yīng)揭示了完全不同的信息年齡效應(yīng)代表疾病隨年齡變化的固有風(fēng)險(xiǎn)時(shí)期效應(yīng)代表所有年齡組同時(shí)受到的影響比如篩查推廣、診斷標(biāo)準(zhǔn)改變隊(duì)列效應(yīng)代表不同出生年份的人由于成長(zhǎng)環(huán)境差異帶來的風(fēng)險(xiǎn)差異。這里要特別提醒一個(gè)容易誤讀的點(diǎn)時(shí)期效應(yīng)和隊(duì)列效應(yīng)存在內(nèi)在共線性。年齡-時(shí)期-隊(duì)列三者的線性趨勢(shì)無法被唯一識(shí)別INLA 通過隨機(jī)效應(yīng)的先驗(yàn)施加了額外的約束所以后驗(yàn)估計(jì)有縮水。因此你在解讀時(shí)期效應(yīng)的下降趨勢(shì)時(shí)不能說是絕對(duì)的下降而應(yīng)該描述為在模型約束下的相對(duì)變化。5.2 預(yù)測(cè)區(qū)間的含義寬區(qū)間不代表模型差I(lǐng)NLA 給出的預(yù)測(cè)帶是后驗(yàn)預(yù)測(cè)區(qū)間它包含了參數(shù)的不確定性和未來數(shù)據(jù)本身的隨機(jī)性。有時(shí)候你看到預(yù)測(cè)區(qū)間在遠(yuǎn)年部分變得很寬甚至上下限翻倍這不是模型報(bào)錯(cuò)而是外推時(shí)間越遠(yuǎn)隊(duì)列和時(shí)期效應(yīng)的不確定性累積越明顯。發(fā)表論文時(shí)不能只畫一個(gè)點(diǎn)預(yù)測(cè)線必須把 95% 區(qū)間帶上。我在實(shí)際處理時(shí)常用 AgePeriodCohort 的繪圖函數(shù)但它的默認(rèn)配色和坐標(biāo)軸刻度不一定合適。下面的代碼演示了怎么提取預(yù)測(cè)結(jié)果并手動(dòng)繪制一條帶區(qū)間的曲線library(ggplot2) # 從 INLA 結(jié)果中提取每年的總體預(yù)測(cè)均值 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] # 預(yù)測(cè)年份的年齡標(biāo)準(zhǔn)化或者總數(shù)需要按年齡匯總 # 這里假設(shè)每個(gè)格子的預(yù)測(cè)值已經(jīng)求和歸一化 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) # 如果希望按年齡加權(quán)合并可先按年聚合 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)這段代碼里最關(guān)鍵的是按 year 聚合前必須確認(rèn)你自己的研究希望按年齡分別展示還是合并展示。如果研究關(guān)注特定年齡組比如 65 歲以上人群的疾病負(fù)擔(dān)就不能直接聚合而是應(yīng)該 filter(age 13) 后再畫這個(gè)年齡組的預(yù)測(cè)序列。我見過有人把青年和老年混在一起預(yù)測(cè)總數(shù)看著平穩(wěn)實(shí)際上年齡內(nèi)部結(jié)構(gòu)已經(jīng)完全變形。5.3 避坑章BAPC 分析中最常見的五個(gè)具體問題問題一安裝 BAPC 時(shí)提示 R 版本過低現(xiàn)象安裝 BAPC 包時(shí)報(bào)錯(cuò) package ‘BAPC’ is not available for this version of R。 原因這臺(tái)機(jī)器上的 R 版本太舊BAPC 的當(dāng)前版本要求 R 4.1 以上。 解決先升級(jí) R 到最新穩(wěn)定版再重新安裝。如果因?yàn)轫?xiàng)目原因不能升 R可以嘗試安裝 BAPC 的歷史版本在 CRAN 的 archive 目錄里找舊版壓縮包用 install.packages(BAPC_1.0.0.tar.gz, repos NULL) 的方式安裝但舊版本可能與新 INLA 不兼容。問題二inla() 報(bào)錯(cuò) Error in inla.model.construct.rw2現(xiàn)象調(diào)用 BAPC() 或 inla() 時(shí)出現(xiàn) rw2 模型構(gòu)造失敗錯(cuò)誤信息里有一長(zhǎng)串矩陣維度。 原因多數(shù)情況下是因?yàn)?age 或者 cohort 因子變量里存在只有單一水平的組INLA 在構(gòu)建隨機(jī)游走精度矩陣時(shí)無法計(jì)算二階差分的懲罰。 解決檢查你傳入 BAPC 的矩陣是否出現(xiàn)了全零行或全零列。比如某一年齡組在整個(gè)時(shí)間段內(nèi)都是 0 發(fā)病數(shù)INLA 會(huì)認(rèn)為該水平?jīng)]有信息量但仍然要參與矩陣運(yùn)算。處理辦法是刪除全零的年齡組或年份或者在數(shù)據(jù)準(zhǔn)備階段用很小的非零值比如 0.01替代純零但這樣做會(huì)影響結(jié)果的解釋我更建議直接刪除完全無數(shù)據(jù)的維度。問題三預(yù)測(cè)值出現(xiàn)負(fù)數(shù)現(xiàn)象預(yù)測(cè)曲線在遠(yuǎn)期出現(xiàn)負(fù)的發(fā)病數(shù)或死亡率明顯不符合生物學(xué)常識(shí)。 原因默認(rèn)情況下 INLA 使用高斯近似作為默認(rèn)似然當(dāng)你沒有指定 family poisson 時(shí)BAPC 包裝函數(shù)可能沒有正確傳遞分布家族導(dǎo)致模型按正態(tài)分布處理了計(jì)數(shù)數(shù)據(jù)。 解決如果你發(fā)現(xiàn)預(yù)測(cè)值有負(fù)值回到代碼檢查是否顯式設(shè)置了 family poisson 或者 family binomial。另外檢查響應(yīng)變量是否真的是整數(shù)計(jì)數(shù)而不是已經(jīng)除以人口數(shù)得到的率。率數(shù)據(jù)應(yīng)該用二項(xiàng)模型或者先進(jìn)行 log 變換不能讓模型把率當(dāng)作正態(tài)分布來擬合。問題四繪圖時(shí)置信區(qū)間突然收窄成一條線現(xiàn)象AgePeriodCohort 畫出的預(yù)測(cè)區(qū)間在某個(gè)年份節(jié)點(diǎn)后突然從寬變窄幾乎貼住點(diǎn)估計(jì)。 原因這是預(yù)測(cè)期和觀測(cè)期銜接處處理方式的問題。如果你的預(yù)測(cè)期數(shù)據(jù)幀里 cohort 編號(hào)沒有正確擴(kuò)展INLA 會(huì)用已有的隊(duì)列效應(yīng)覆蓋新隊(duì)列導(dǎo)致預(yù)測(cè)的變化自由度降低。 解決檢查 cohort_id 的最大值是否小于未來預(yù)測(cè)所需的隊(duì)列編號(hào)最大值。比如觀測(cè)期是 1990-2019 年年齡組 0-4 歲的隊(duì)列編號(hào)范圍是 1986-2019但預(yù)測(cè)到 2034 年時(shí)會(huì)出現(xiàn)新的隊(duì)列編號(hào) 2030 左右如果沒有把這些新編號(hào)加入模型繪圖時(shí)區(qū)間就會(huì)被壓縮。問題五R 會(huì)話崩潰或內(nèi)存耗盡現(xiàn)象跑 inla() 到一半 RStudio 直接閃退或者提示 cannot allocate vector of size。 原因INLA 本身對(duì)內(nèi)存的占用較高尤其是當(dāng)你用了很細(xì)的年齡組1 歲一組且年份跨度長(zhǎng)時(shí)精度矩陣的維度可能達(dá)到幾千乘幾千稀疏矩陣存儲(chǔ)也要幾百 MB。更常見的是同時(shí)打開了太多 R 進(jìn)程或者瀏覽器占用了大量?jī)?nèi)存。 解決先減少數(shù)據(jù)維度比如把 1 歲年齡組合并成 5 歲一組或者把預(yù)測(cè)年限縮短到 10 年。另外在調(diào)用 inla() 前設(shè)置 inla.setOption(num.threads 2) 限制線程數(shù)防止 INLA 默認(rèn)占用所有 CPU 核心導(dǎo)致系統(tǒng)假死。如果內(nèi)存實(shí)在不夠用 64 位版本的 R 并增加虛擬內(nèi)存。6. 驗(yàn)證預(yù)測(cè)效果的一個(gè)實(shí)用技巧回帶歷史數(shù)據(jù)做軌道驗(yàn)證模型跑通不是終點(diǎn)關(guān)鍵要驗(yàn)證預(yù)測(cè)可信度。我習(xí)慣做一個(gè)簡(jiǎn)單的后驗(yàn)驗(yàn)證把觀測(cè)數(shù)據(jù)截?cái)嗟侥骋荒曛氨热鐚?shí)際數(shù)據(jù)到 2019 年我只用 2000-2014 年擬合模型預(yù)測(cè) 2015-2019 年的結(jié)果然后跟真實(shí)觀測(cè)值對(duì)比。這樣可以直觀看到模型在歷史階段的預(yù)測(cè)誤差有多大也能側(cè)面判斷你的年齡分組和模型設(shè)定是否合理。# 截?cái)鄶?shù)據(jù)僅使用 2000-2014 年作為擬合期 train_wide - wide[, colnames(wide) %in% as.character(2000:2014)] # 注意保持年齡順序不變 train_cohort_id - cohort_id[, 1:ncol(train_wide)] # 重新跑一遍 BAPC預(yù)測(cè)到 2019 年 bapc_val - BAPC(train_wide, model APC, predict TRUE, prediction list(year 2015:2019)) # 提取預(yù)測(cè)均值和觀測(cè)真實(shí)值 val_pred - bapc_val$prediction$rate val_obs - wide[, colnames(wide) %in% as.character(2015:2019)] # 計(jì)算相對(duì)誤差 rel_error - (val_pred - val_obs) / val_obs * 100這段驗(yàn)證代碼里的 prediction 參數(shù)接受一個(gè)列表可以指定預(yù)測(cè)的年份范圍。我一般會(huì)計(jì)算分年齡組的相對(duì)誤差看是不是特定年齡組誤差特別大。如果某個(gè)年齡組的誤差超過 30%就要考慮數(shù)據(jù)本身在該年齡組是否存在報(bào)告不全或者模型對(duì)隊(duì)列效應(yīng)的假設(shè)在這個(gè)年齡組不適用。比如中老年組往往因?yàn)樵\斷率提高導(dǎo)致時(shí)期效應(yīng)很強(qiáng)可能需要把時(shí)期效應(yīng)的隨機(jī)游走從 rw1 改成 rw2 來獲得更平滑的趨勢(shì)。另外還有一種更輕量的驗(yàn)證方式比較 BAPC 的預(yù)測(cè)結(jié)果和用簡(jiǎn)單負(fù)二項(xiàng)回歸外推的結(jié)果。雖然沒有概率上的嚴(yán)謹(jǐn)對(duì)比但可以看出貝葉斯模型的預(yù)測(cè)區(qū)間是否明顯比樸素方法寬以及中期預(yù)測(cè)點(diǎn)估計(jì)是否落在樸素方法的置信區(qū)間內(nèi)。這不會(huì)寫進(jìn)論文但能幫你判斷模型有沒有系統(tǒng)性偏倚。如果驗(yàn)證發(fā)現(xiàn)誤差在可接受范圍內(nèi)接下來才用全量數(shù)據(jù)重新擬合并把預(yù)測(cè)年限設(shè)為研究需要的長(zhǎng)度。經(jīng)過這個(gè)流程你對(duì)自己模型輸出的解釋就更有底氣了。從那以后我每次做 BAPC 分析都強(qiáng)制把回帶驗(yàn)證這一步加在正式預(yù)測(cè)之前哪怕結(jié)果已經(jīng)出來也返工重跑這個(gè)習(xí)慣幫我避免過至少一次因?yàn)槟挲g組順序錯(cuò)亂導(dǎo)致的完全錯(cuò)誤的趨勢(shì)預(yù)測(cè)。希望這些經(jīng)驗(yàn)?zāi)軒湍阍?GBD 數(shù)據(jù)的 BAPC 分析上少走彎路。本文還有配套的精品資源點(diǎn)擊獲取