簡介一份可安裝的R包源碼資源專門用于在R環(huán)境中調(diào)用JAGS完成貝葉斯統(tǒng)計分析適合具備R基礎(chǔ)并需要處理生態(tài)學(xué)、野生動物種群或統(tǒng)計學(xué)問題的中高級用戶。該包在rjags基礎(chǔ)上提供了更簡潔的封裝接口涵蓋數(shù)據(jù)檢查、參數(shù)初始化、模型運(yùn)行、自動續(xù)跑、后驗預(yù)測檢查、收斂診斷和圖形輸出等完整流程同時支持多條馬爾可夫鏈并行計算可明顯縮短復(fù)雜模型的運(yùn)行時間。資源共包含47個文件其中31個R腳本是核心實現(xiàn)兼顧高層函數(shù)與內(nèi)部工具9個Rd文檔提供每個函數(shù)的規(guī)范說明便于查閱和二次開發(fā)其余包括命名空間、DESCRIPTION、NEWS、構(gòu)建忽略文件等保證了R包結(jié)構(gòu)的完整性和可安裝性整個壓縮包僅44KB。讀者既可以直接加載使用快速完成模型擬合與結(jié)果展示也可以通過源碼學(xué)習(xí)JAGS接口的封裝技巧、MCMC并行策略及R包組織結(jié)構(gòu)。當(dāng)前已有955人參與學(xué)習(xí)或下載。 做貝葉斯數(shù)據(jù)分析的人大概率都繞不開MCMC。而在R語言生態(tài)里想跑JAGSJust Another Gibbs Sampler的模型jagsUI幾乎是我見過最省心的接口包。你不需要在R和JAGS之間來回倒文件不需要記一堆底層命令只要把模型、數(shù)據(jù)、參數(shù)丟給一個函數(shù)它就能自動完成采樣、收斂診斷和結(jié)果匯總。這幾年我用R語言做數(shù)據(jù)分析在需要快速出貝葉斯結(jié)果的場景下jagsUI一直是我最習(xí)慣的選擇。這篇文章就把我實際使用的經(jīng)驗、踩過的坑和核心用法一次講清楚希望能幫到剛接觸這塊的朋友。1. 為什么用jagsUI而不是直接用JAGS或其他包1.1 JAGS是什么為什么要在R里調(diào)用JAGS的全稱是Just Another Gibbs Sampler是一個用BUGS語言寫模型、用MCMC方法做貝葉斯推斷的獨(dú)立軟件。你寫好一個模型文件它負(fù)責(zé)編譯模型、生成采樣器、跑出后驗分布。但問題是JAGS自己不帶R那種方便的數(shù)據(jù)處理和可視化能力你要手動寫腳本管理數(shù)據(jù)和輸出非常別扭。R的好處是數(shù)據(jù)清洗、畫圖、報告生成都在一個環(huán)境里所以把JAGS嵌到R工作流里是很多人的剛需。你可以在R里整理數(shù)據(jù)調(diào)用JAGS跑MCMC再把后驗結(jié)果拿出來畫圖或做假設(shè)檢驗整個流程不用切換軟件。jagsUI就是這個“嵌入口”的一種實現(xiàn)。1.2 常見R接口包橫向?qū)Ρ萊里能調(diào)JAGS的包不止一個我最早用的是rjags后來也試過R2jags和runjags每個都有自己的特點(diǎn)。包底層封裝突出優(yōu)點(diǎn)缺點(diǎn)rjags直接封裝JAGS C接口靈活、穩(wěn)定底層控制力強(qiáng)寫起來繁瑣要自己處理模型更新和收斂判斷R2jags基于rjags提供jags()函數(shù)用法簡單輸出對象整合一般功能相對有限r(nóng)unjags基于rjags功能非常全支持并行、自動收斂擴(kuò)展參數(shù)復(fù)雜新手容易繞暈jagsUI基于rjags語法簡潔自動輸出Rhat、有效樣本量等統(tǒng)計量高級定制不如底層包靈活對比下來你會發(fā)現(xiàn)jagsUI不是功能最多的那個但它是綜合體驗最“現(xiàn)代”的。它把rjags里需要手動做的很多步驟封裝成了默認(rèn)行為比如自動判定模型是否有離散節(jié)點(diǎn)、自動生成初值、自動判斷哪些參數(shù)需要監(jiān)控這些設(shè)計讓入門門檻低了一大截。1.3 我選jagsUI的理由我個人更看重的是“出結(jié)果的速度”。實際項目里貝葉斯模型只是分析鏈路中的一環(huán)我不希望把大量時間花在接口包的使用細(xì)節(jié)上。jagsUI的jags()函數(shù)一次調(diào)用就能完成模型編譯、預(yù)熱、采樣、統(tǒng)計匯總返回的對象里直接帶著mean、sd、q2.5、q97.5、Rhat和n.eff拿來就能寫報告。另外jagsUI支持并行跑多條MCMC鏈在多核CPU上能明顯縮短等待時間這在跑復(fù)雜模型時非常關(guān)鍵。對于團(tuán)隊協(xié)作項目用jagsUI的人不需要額外熟悉一整套rjags命令代碼可讀性也更好。當(dāng)然如果你要高度定制采樣器或轉(zhuǎn)化器rjags可能更合適但日常絕大多數(shù)數(shù)據(jù)分析場景jagsUI足夠了。2. 環(huán)境準(zhǔn)備與安裝從JAGS本體到R包2.1 安裝JAGS獨(dú)立程序jagsUI只是R語言層面的接口真正干活的還是JAGS本體所以第一步是安裝JAGS。這里最容易被新手忽略光在R里裝包是不夠的系統(tǒng)里沒有JAGS可執(zhí)行文件后面運(yùn)行必然報錯。JAGS支持Windows、macOS和Linux。Linux用戶一般可以直接用軟件源安裝比如Ubuntu下執(zhí)行sudo apt install jagsmacOS用戶可以用Homebrew執(zhí)行brew install jags。Windows用戶需要去JAGS官網(wǎng)或CRAN的鏈接里下載安裝包安裝時記住安裝路徑我一般建議默認(rèn)路徑后續(xù)省事。下載時注意選擇對應(yīng)R版本的64位版現(xiàn)在基本都用64位。安裝完成后可以在命令行里輸入jags或查看安裝目錄下的JAGS.exe來確認(rèn)。不需要手動配置環(huán)境變量jagsUI在Windows下通常能自動找到JAGS的安裝位置。但如果你用的是定制安裝路徑后文會講怎么通過JAGS_HOME環(huán)境變量處理。2.2 安裝jagsUI并驗證在R環(huán)境里安裝jagsUI非常簡單install.packages(jagsUI)它會自動依賴rjags和coda等包CRAN上的版本一般都很新。裝完后加載并驗證一下是否能找到JAGSlibrary(jagsUI) # 用一個小模型測試 mod - jags( model.file textConnection( model { y ~ dnorm(0, 1) } ), data list(y 1), parameters.to.save y, n.chains 1, n.iter 100, n.burnin 0, n.adapt 10, verbose FALSE )如果這段代碼能順利跑完說明JAGS和jagsUI都安裝成功了。你會在返回的mod對象里看到$mean等輸出。2.3 環(huán)境變量與常見安裝坑一開始不知道JAGS路徑時我踩過一個坑裝完jagsUI后運(yùn)行模型直接報Error in jags.model(...): JAGS not found。后來才發(fā)現(xiàn)沒有把JAGS的安裝目錄告訴R。Windows下可以把JAGS安裝目錄加入環(huán)境變量比如JAGS_HOME C:/Program Files/JAGS/JAGS-4.3.0/x64或者直接把JAGS_HOME指向包含JAGS.exe的目錄。macOS如果從源碼編譯安裝也可能出現(xiàn)找不到JAGS的情況這時候用Sys.setenv(JAGS_HOME /usr/local/bin)暫時設(shè)置即可。不過說到底大部分情況默認(rèn)安裝路徑就行真遇到再排查環(huán)境變量不必一開始就折騰。3. 核心用法模型、數(shù)據(jù)與jags()函數(shù)全解析3.1 模型文件BUGS語言快速上手jagsUI要求你寫一個BUGS風(fēng)格的模型文件。這個文件可以放在磁盤上也可以直接用一個R字符串傳入。我通常在項目里單獨(dú)維護(hù)一個model.txt方便復(fù)用。一個最簡單的均值估計模型長這樣model { for (i in 1:N) { y[i] ~ dnorm(mu, tau) } mu ~ dnorm(0, 0.001) tau - 1 / (sigma * sigma) sigma ~ dunif(0, 100) }注意JAGS里的正態(tài)分布參數(shù)是均值和精度方差的倒數(shù)不是標(biāo)準(zhǔn)差。我一開始寫的時候經(jīng)常會順手寫成dnorm(mu, sigma)然后就發(fā)現(xiàn)后驗方差被嚴(yán)重高估。這里tau是精度所以采樣器里通常要設(shè)一個sigma的先驗再用tau - 1 / sigma^2轉(zhuǎn)換。模型文件中最重要的是“給每個參數(shù)指定先驗分布”JAGS會檢查模型閉合。如果某個節(jié)點(diǎn)沒有先驗?zāi)Ｐ途幾g會報錯。你的先驗選擇也直接影響MCMC收斂后面我會專門講。3.2 數(shù)據(jù)列表、初值與參數(shù)監(jiān)控數(shù)據(jù)必須整理成一個list名稱要和模型里的變量名嚴(yán)格對應(yīng)。比如上面模型需要y和N所以R里要準(zhǔn)備data_list - list( y c(3.2, 3.8, 2.9, 4.1, 3.5), N 5 )parameters.to.save用來告訴jagsUI你關(guān)心哪些參數(shù)/節(jié)點(diǎn)。比如parameters.to.save c(mu, sigma)它就只監(jiān)控這兩個節(jié)點(diǎn)的后驗。注意像tau這種確定性節(jié)點(diǎn)也可以監(jiān)控但沒必要。如果你好奇預(yù)測值可以把缺失值設(shè)為NAJAGS會自動當(dāng)作缺失數(shù)據(jù)預(yù)測但JAGS實際是用NA作為參數(shù)采樣這個特性可以用來做后驗預(yù)測。初值在jagsUI里可以完全交給它自動生成。它默認(rèn)會為隨機(jī)節(jié)點(diǎn)生成合理的初值但有時候復(fù)雜模型還是要手動指定比如給離散參數(shù)一個合適的整數(shù)初值。手動指定時你可以傳一個包含與n.chains等長list的inits參數(shù)。3.3 jags()函數(shù)關(guān)鍵參數(shù)逐個說jags()是核心函數(shù)參數(shù)很多但日常最常用的就這幾個model.file模型文件路徑或連接對象。data命名列表。parameters.to.save要監(jiān)控的參數(shù)名。n.chainsMCMC鏈數(shù)我一般設(shè)3或4。n.iter總迭代次數(shù)包含預(yù)熱階段。n.burnin預(yù)熱的迭代次數(shù)一般占總迭代的20%~50%。n.thin采樣間隔用來降低自相關(guān)。一般n.thin 1即可如果自相關(guān)高再調(diào)大。parallel是否并行跑多鏈設(shè)為TRUE能大幅提速。seed隨機(jī)種子讓結(jié)果可復(fù)現(xiàn)。我經(jīng)常這么設(shè)mod - jags( model.file model.txt, data data_list, parameters.to.save c(mu, sigma), n.chains 3, n.iter 10000, n.burnin 2000, n.thin 1, parallel TRUE, seed 123 )這樣每個參數(shù)會得到3 × (10000 - 2000) 24000個有效迭代樣本。parallel TRUE會讓每條鏈跑在獨(dú)立核心上但注意它消費(fèi)內(nèi)存復(fù)雜模型時別把核心數(shù)開得太大否則容易卡死。4. 實操案例用jagsUI跑一個線性回歸4.1 模擬數(shù)據(jù)與模型設(shè)定理論講再多不如實際跑一遍。我模擬一組簡單線性回歸數(shù)據(jù)想估計截距、斜率和方差。首先在R里造數(shù)據(jù)set.seed(42) N - 100 x - rnorm(N, 0, 1) true_a - 1.5 true_b - 2.0 true_sigma - 1.2 y - rnorm(N, true_a true_b * x, true_sigma)模型文件lm_model.txt內(nèi)容model { for (i in 1:N) { y[i] ~ dnorm(a b * x[i], tau) } a ~ dnorm(0, 0.001) b ~ dnorm(0, 0.001) tau - 1 / (sigma * sigma) sigma ~ dunif(0, 50) }這里我故意給b一個比較寬的弱先驗dnorm(0, 0.001)相當(dāng)于方差1000的正態(tài)分布表示我并沒有對斜率有太強(qiáng)的主觀預(yù)設(shè)。4.2 運(yùn)行模型與輸出解讀數(shù)據(jù)準(zhǔn)備和運(yùn)行data_list - list(y y, x x, N N) mod - jags( model.file lm_model.txt, data data_list, parameters.to.save c(a, b, sigma), n.chains 3, n.iter 10000, n.burnin 2000, parallel TRUE, seed 1 )運(yùn)行結(jié)束后mod對象里已經(jīng)包含所有統(tǒng)計量。你可以用print(mod)查看完整匯總也可以用mod$summary直接取數(shù)據(jù)框。我經(jīng)常用的幾個字段mod$mean后驗均值相當(dāng)于點(diǎn)估計。mod$q2.5和mod$q97.5后驗95%可信區(qū)間。mod$Rhat收斂診斷值一般要求小于1.1。mod$n.eff有效樣本量太小說明自相關(guān)嚴(yán)重。我這個模擬里后驗均值大概在a1.4~1.6、b1.9~2.1、sigma1.1~1.3的范圍內(nèi)和真實值比較接近說明模型恢復(fù)參數(shù)的能力是OK的。4.3 收斂診斷與后驗可視化MCMC跑完不能直接信結(jié)果要先看收斂。我一般會看Rhat是否都小于1.1再看n.eff有沒有低于幾百的。jagsUI還提供了traceplot(mod)函數(shù)可以快速看鏈的軌跡圖。軌跡圖要像毛毛蟲一樣來回扭動而不是一條直線或分段的水平線后者說明鏈卡在某個區(qū)域要么模型寫錯了要么先驗和似然沖突。后驗分布可視化我習(xí)慣轉(zhuǎn)成數(shù)據(jù)框再畫library(ggplot2) df - as.data.frame(mod$samples) ggplot(df, aes(x b)) geom_density(fill steelblue, alpha 0.4) labs(x 斜率 b)mod$samples返回的是mcmc.list對象as.data.frame()會把多鏈樣本合并成一個大數(shù)據(jù)框方便ggplot直接畫。你還可以畫后驗密度曲線疊加真實值或者畫斜率和截距的二維等高線觀察參數(shù)之間的相關(guān)性這些都是貝葉斯分析里很有價值的內(nèi)容。5. 常見問題與排查技巧實錄5.1 安裝與路徑問題最常遇到的錯誤是Could not find JAGS。這基本是JAGS本體沒裝好或R找不到JAGS路徑。Windows下可以用Sys.setenv(JAGS_HOME C:/Program Files/JAGS/JAGS-4.3.0/x64)指定注意路徑要寫到包含JAGS.exe那一層。如果你的R是32位而JAGS是64位也會出現(xiàn)連接問題盡量保證一致。另一個容易踩的坑是更新R之后舊包的二進(jìn)制不兼容。遇到package or namespace load failed時先試試重啟R再不行就重裝jagsUI和rjags。5.2 模型運(yùn)行報錯與排查模型編譯時報Error in node通常是模型里用了未定義的數(shù)據(jù)變量或者數(shù)據(jù)列表里變量名寫錯。比如模型寫了x[i]但你數(shù)據(jù)列表里寫的是X哪怕是大小寫不一致JAGS都會直接報錯。養(yǎng)成習(xí)慣模型文件里的變量名和list里的命名必須逐字對齊。另一個常見錯誤是Unknown variable。這個多半是因為數(shù)據(jù)列表只給了模型需要的部分變量。初值報錯也比較多尤其是指定初值時給了非數(shù)值或超出先驗范圍的初值。比如sigma的先驗是dunif(0, 50)初值給了負(fù)值JAGS就不干了。自動初值時遇到離散節(jié)點(diǎn)有時候也會出問題解決辦法是手動給離散參數(shù)設(shè)置整數(shù)初值。5.3 輸出處理與項目建議輸出里Rhat是NaN時先別慌可能是某條鏈在采樣時全部退化了或者n.eff太低。我遇到過因為模型參數(shù)化不當(dāng)導(dǎo)致多條鏈一直發(fā)散的情況后來把斜率的先驗從dnorm(0, 0.001)改成dnorm(0, 0.01)收斂就好了很多。所以當(dāng)Rhat異常時優(yōu)先檢查模型設(shè)定和數(shù)據(jù)標(biāo)準(zhǔn)化而不是盲目增大迭代次數(shù)。處理大型項目時我建議把模型文件、數(shù)據(jù)準(zhǔn)備、運(yùn)行腳本和結(jié)果輸出分開存放。jagsUI跑完后的對象可能很大尤其是mod$samples全部保存在內(nèi)存里會占很多空間?？梢韵萻aveRDS(mod, mod.rds)存檔后續(xù)分析再從文件讀取這樣R會話能輕松不少。結(jié)尾我在實際使用jagsUI的過程中最深的體會是它把貝葉斯分析的門檻降得很低尤其是對R語言用戶。你不需要先學(xué)完整個JAGS語法才能動手只需要會寫簡單的BUGS模型然后通過一個jags()函數(shù)和R的數(shù)據(jù)框無縫銜接。當(dāng)然它也不是萬能鑰匙遇到特別復(fù)雜的模型或高度非正態(tài)的后驗時還是需要回到底層包甚至換用Stan這類基于HMC的引擎。但如果你現(xiàn)在的工作流是R 貝葉斯統(tǒng)計并且需要一個開箱即用的JAGS接口jagsUI值得放進(jìn)你的工具箱。最后再分享一個小技巧處理新數(shù)據(jù)時先用很小的n.iter跑一次確認(rèn)模型能編譯、Rhat不爆再加大迭代量正式運(yùn)行這樣能省下大量調(diào)試時間。本文還有配套的精品資源點(diǎn)擊獲取