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