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