JAGS(Just Another Gibbs Sampler)是為彌補BUGS軟件無法在Unix或Linux等非Windows系統上運行的缺陷而研發的編程軟件,其自身獨立擁有一套貝葉斯理論運算函數及公式,具有操作界面簡單、與系統兼容好、運行流暢、且與其他編程軟件具有良好交互性等特點。然而,由于JAGS軟件在結果數據讀取、解讀以及圖形繪制方面存在功能性缺失,限制了其推廣應用;通過R2jags、rjags及runjags程序包從R中調用JAGS軟件則克服了這些缺陷。R2jags、rjags及runjags這3款程序包的工作原理均是在R語言的框架下通過調用JAGS軟件來完成相關迭代運算,它們的功能結構基本相似,都具有操作簡單、命令簡潔、數據讀取和結果匯總以及圖形繪制功能較完善的特點,但在實踐運用中仍存在差異。本文介紹了如何使用R軟件通過R2jags、rjags及runjags程序包調用JAGS軟件來實現網狀Meta分析。
引用本文: 張超, 孫鳳, 曾憲濤. R軟件調用JAGS軟件實現網狀Meta分析. 中國循證醫學雜志, 2014, 14(2): 241-248. doi: 10.7507/1672-2531.20140042 復制
JAGS(Just Another Gibbs Sampler)是一款基于C++語言編寫的、自身獨立擁有一套完整的貝葉斯理論運算函數及公式的編程軟件。其具有操作界面簡單、與系統兼容好、運行流暢、且與其他編程軟件具有良好交互性等特點;但由于該軟件自身缺乏對結果數據讀取、解讀以及圖形繪制等功能,使得其應用及推廣較差[1]。因JAGS軟件具有良好的外接性,研究者將其與數據讀取能力強大及繪圖功能完善的編程軟件R相結合,使其能夠較完美地實現網狀Meta分析[2]。
當前,通過R軟件調用JAGS軟件的程序包有R2jags、rjags、runjags和bayesmix,這4個程序包的功能結構基本相似,都具有操作簡單、命令簡潔、數據讀取和匯總以及結果圖形繪制功能較完善的特點。通過這4個程序包,可以實現多種類型的Meta分析[3, 4]。本文主要介紹如何通過R2jags、rjags和runjags這3個程序包實現網狀Meta分析,實例及數據詳見《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中表 1。

1 軟件安裝及程序包加載
JAGS軟件是一款自由獲取的免費軟件,可從http://mcmc-jags.sourceforge.net/下載。依據R調用原理[
R及JAGS安裝完成后,需要加載R2jags、rjags和runjags程序包。此3款程序包下載與加載的代碼見表 1。
2 軟件操作
2.1 模塊加載
本文的模塊與《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中所采取的模塊一致,但兩者的唯一區別在于:JAGS軟件自身沒有排列函數。因此,JAGS軟件需要在原始代碼中將排列函數刪除。JAGS軟件模塊加載的命令見框?1。
框1
model {for (i in 1:ns) ?{w[i, 1] < -0 ??delta[i, 1]<-0 ??mu[i]~dnorm(0, 0.0001) ???for(k in 1:na[i]) ????{r[i, k]~dbin(p[i, k], n[i, k]) ????logit(p[i, k]) < -mu[i]+delta[i, k] ????rhat[i, k] < -p[i, k]*n[i, k] dev[i, k] < -2*(r[i, k]*(log(r[i, k])-log(rhat[i, k]))+(n[i, k]-r[i, k])*(log(n[i, k]-r[i, k])-log(n[i, k]-rhat[i, k]))) ????} ??resdev[i] < -sum(dev[i, 1:na[i]]) ???for(k in 2:na[i]) ????{delta[i, k]~dnorm(md[i, k], taud[i, k]) ????md[i, k] < -d[t[i, k]]-d[t[i, 1]]+sw[i, k] ????taud[i, k] < -tau*2*(k-1)/k ????w[i, k] < -(delta[i, k]-d[t[i, k]]+d[t[i, 1]]) ????sw[i, k] < -sum(w[i, 1:k-1])/(k-1) ????} } totresdev < -sum(resdev[]) d[1] < -0 for(k in 2:nt) ?{d[k]~dnorm(0, 0.0001) ?} sd~dunif(0, 5) tau < -pow(sd, -2) for (c in 1:(nt-1)) ?{for(k in (c+1):nt) ????{or[c, k] < -exp(d[k]-d[c]) ????lor[c, k] < -(d[k]-d[c]) ????} }
處理完畢后,將其保存為:“.bug”或“.jags”格式。本例以“networkmodel”為名、以“.jags”格式保存在桌面的Rwork文件夾下。
根據程序包的工作原理,其最終迭代抽樣運算仍是使用JAGS軟件來實踐完成的。因此本處還需在R中對模塊進行加載,具體命令如下:
bugsname < -file.path("C:/Users/Administrator/Desktop/Rwork", "networkmodel.jags")
2.2 數據及初始值加載
本例3個程序包雖均采取《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中的數據和初始值,但相關的格式卻有所不同:
對于R2jags及rjags程序包而言,本例數據及初始值的排列格式及加載方式與上文中步驟2.2~2.5 [5]一致,本處不再講解。
對于runjags程序包,本例數據及初始值的排列格式及加載方式與上文中步驟2.2~2.4 [5]一致;但還需對數據及初始值格式作相關整合,具體命令分別如下:
data < -dump.format(list(t=t, r=r, n=n, na=na, ns=ns, nt=nt))
inits < -dump.format(list(d=c(NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), sd=1, mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)))
2.3 觀察指標設定
由于在上述模塊編譯中,JAGS軟件對lor和or數據集無法進行識別與監視,因此,此處需要將lor和or數據集中每個節點一一寫入監視。本例為了節省計算機內存資源以及方便后續結果展示,僅監視lor[1, 2],or[1, 2],lor[1, 3],or[1, 3]節點作為示范,監視代碼設定為:
parameters < -c("lor[1, 2]", "or[1, 2]", "lor[1, 3]", "or[1, 3]")。
2.4 迭代運算
此3個程序包中,雖在迭代運算基本思路與原理是相同的,但其迭代命令代碼仍有所不同。表 2及圖 1展示了這3個程序包迭代命令及迭代進程。


2.5 結果顯示及匯總
在上述執行迭代命令運行完成后,我們可以對結果進行顯示及匯總了。具體命令代碼見表 3,相關匯總結果見圖 2~圖 5。





3 結果圖形繪制
如前所述,圖形的繪制是JAGS軟件的短板,R則可在這方面給予其很好的補充,且實現方法多樣。本例將提供3種繪制結果圖的方法。
3.1 方法一
具體操作命令見表 4,執行命令后即可繪制相關的圖形(本處僅展示R2jags程序包繪制的森林圖,其余圖略)。

值得注意的是:盡管3個程序包雖均能繪制結果圖,但仍需相互結合使用方能繪制出森林圖、軌跡圖及密度圖。

3.2 方法二
首先,提取各監視節點的軌跡數據,再依據該數據進行軌跡圖與密度圖的繪制。提取各監視節點的軌跡數據命令見表 5。

提取監視節點迭代軌跡數據后,即可依據該數據繪制軌跡圖與密度圖,具體操作命令見表 6,結果圖略。

3.3 方法三
通過R軟件的coda程序包中的codamenu()菜單型命令進行結果圖形的繪制,首先需要加載coda程序包,命令如下:
install.packages(“coda”)
library(“coda”)
加載完成,需要進行結果的讀取,需要使用coda程序包中codamenu()菜單型命令,具體命令如下:
codamenu()
在執行上述命令后,將會出現以下3個菜單提供選擇:
1: Read BUGS output files
2: Use an mcmc object
3: Quit
Selection:
上述3個選擇的含義分別為:1,讀取BUGS軟件生成文件;2,使用mcmc函數對象;3,退出。在下方自動出現的“Selection:”后填寫2,來讀取“.mcmc”后綴結果數據。執行“Selection: 2”命令后,出現以下提示字樣:
Enter name of saved object (or type "exit" to quit)
1:
在“1:”的后方鍵入“.mcmc”后綴格式的結果數據變量集(即分別鍵入jagsfit.mcmc、coda.mcmc和run.jags.mcmc)并執行上述命令后,軟件會自動提示以下操作主菜單(標記為1),具體如下:
CODA Main Menu
1: Output Analysis
2: Diagnostics
3: List/Change Options
4: Quit
Selection:
在主菜單的目錄中,含義分別為:1,為結果分析菜單(標記為1.1);2,為結果診斷菜單(標記為1.2);3,對應的為該節點的列表信息,具體使用價值不大(有興趣讀者可自行嘗試);4,退出。此時,在“Selection:”后填寫1,將會出現子級菜單[即Output Analysis為結果分析菜單(標記為1.1)]供選擇,具體如下:
CODA Output Analysis menu
1: Plots
2: Statistics
3: List/Change Options
4: Return to Main Menu
Selection:
其中,選擇不同選項產生的結果不同:1,表示該節點軌跡圖與密度圖的繪制;2,表示該節點的相關信息匯總,其包含有Mean、SD、Na?ve SE、Time-series SE以及2.5%、25%、50%、75%、97.5%等分位數(見圖 7);3,對應的為該節點的列表信息;4,為返回主菜單。在“Selection:”后填寫2,將會出現下一子級菜單[即Diagnostics為結果診斷菜單(標記為1.2)]供選擇,具體如下:

CODA Diagnostics Menu
1: Geweke
2: Gelman and Rubin
3: Raftery and Lewis
4: Heidelberger and Welch
5: Autocorrelations
6: Cross-Correlations
7: List/Change Options
8: Return to Main Menu
Selection:
該菜單中:1~4均為各種收斂診斷方法;5為自相關性診斷;6為交互相關性;7為自行嘗試;8為返回主菜單。在“Selection:”處填寫5,即可顯示其自相關性的結果(見圖 8)并出現以下子菜單:

Autocorrelation Plots Menu
1: Plot autocorrelations
2: Return to Diagnostics Menu
Selection:
在其子菜單中(標記為1.2.5)下,若在“Selection:”處繼續選擇1時,即可繪制自相關性圖形。上述菜單中,還可繼續嵌套相關的子菜單,在此就不一一介紹,有興趣的讀者可自行嘗試。
4 結語
通過上述介紹,我們發現JAGS軟件與R的結合,使JAGS軟件生成的結果文件能在R軟件中得到完美的解讀,無論是結果準確性的解讀還是圖形精美繪制都可以以菜單命令形式展現,這也是在編程軟件中少見的,對于編程基礎語言欠缺的操作者來說,十分便捷。
然而,由于此3個程序包均是基于JAGS軟件來實現的,所以其共同的缺陷在于只能監視單一節點,而不能監視節點變量。與調用WinBUGS和OpenBUGS的程序相比[5, 7, 8],在節點較多的分析中,使用繁瑣;此外,這3個程序包以及JAGS軟件本身均缺乏rank函數功能[1],故無法直接進行排序,需要通過其他方式或運行其他代碼來實現。
從繪圖方面而言,本文展示的3種結果圖繪制方法中前兩種的方法操作相對較為簡單,但圖形的繪制需要相互補充;第三種方法功能較為全面,其不僅具有繪制多種圖形的能力,同時還能對其結果進行匯總并給與相關分析。故筆者建議使用第三種方法來完成相關操作,具體選取可依使用者需求而定。
細心的讀者會發現,在引言中筆者還提及了bayesmix程序包。然而,bayesmix程序包適用于單變量的混合模型數據且該程序包的model依據數據及參數的簡單設定在后臺自動生成,無法加載網狀model;再者,bayesmix程序包的實質還是運用rjags程序包來實現對JAGS軟件的調用。因此,本文未對其進行介紹。
總之,R2jags、rjags及runjags程序包均依附于R語言框架下,通過調用JAGS軟件的方式來實現JAGS軟件在網狀Meta分析中的運用。其不僅減少了JAGS軟件的一些不必要的繁瑣操作,而且彌補了該軟件在結果讀取、分析、匯總以及圖形繪制等方面的功能缺失;其為操作者在基于貝葉斯理論軟件的選取上又提供了新的友好選擇。
JAGS(Just Another Gibbs Sampler)是一款基于C++語言編寫的、自身獨立擁有一套完整的貝葉斯理論運算函數及公式的編程軟件。其具有操作界面簡單、與系統兼容好、運行流暢、且與其他編程軟件具有良好交互性等特點;但由于該軟件自身缺乏對結果數據讀取、解讀以及圖形繪制等功能,使得其應用及推廣較差[1]。因JAGS軟件具有良好的外接性,研究者將其與數據讀取能力強大及繪圖功能完善的編程軟件R相結合,使其能夠較完美地實現網狀Meta分析[2]。
當前,通過R軟件調用JAGS軟件的程序包有R2jags、rjags、runjags和bayesmix,這4個程序包的功能結構基本相似,都具有操作簡單、命令簡潔、數據讀取和匯總以及結果圖形繪制功能較完善的特點。通過這4個程序包,可以實現多種類型的Meta分析[3, 4]。本文主要介紹如何通過R2jags、rjags和runjags這3個程序包實現網狀Meta分析,實例及數據詳見《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中表 1。

1 軟件安裝及程序包加載
JAGS軟件是一款自由獲取的免費軟件,可從http://mcmc-jags.sourceforge.net/下載。依據R調用原理[
R及JAGS安裝完成后,需要加載R2jags、rjags和runjags程序包。此3款程序包下載與加載的代碼見表 1。
2 軟件操作
2.1 模塊加載
本文的模塊與《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中所采取的模塊一致,但兩者的唯一區別在于:JAGS軟件自身沒有排列函數。因此,JAGS軟件需要在原始代碼中將排列函數刪除。JAGS軟件模塊加載的命令見框?1。
框1
model {for (i in 1:ns) ?{w[i, 1] < -0 ??delta[i, 1]<-0 ??mu[i]~dnorm(0, 0.0001) ???for(k in 1:na[i]) ????{r[i, k]~dbin(p[i, k], n[i, k]) ????logit(p[i, k]) < -mu[i]+delta[i, k] ????rhat[i, k] < -p[i, k]*n[i, k] dev[i, k] < -2*(r[i, k]*(log(r[i, k])-log(rhat[i, k]))+(n[i, k]-r[i, k])*(log(n[i, k]-r[i, k])-log(n[i, k]-rhat[i, k]))) ????} ??resdev[i] < -sum(dev[i, 1:na[i]]) ???for(k in 2:na[i]) ????{delta[i, k]~dnorm(md[i, k], taud[i, k]) ????md[i, k] < -d[t[i, k]]-d[t[i, 1]]+sw[i, k] ????taud[i, k] < -tau*2*(k-1)/k ????w[i, k] < -(delta[i, k]-d[t[i, k]]+d[t[i, 1]]) ????sw[i, k] < -sum(w[i, 1:k-1])/(k-1) ????} } totresdev < -sum(resdev[]) d[1] < -0 for(k in 2:nt) ?{d[k]~dnorm(0, 0.0001) ?} sd~dunif(0, 5) tau < -pow(sd, -2) for (c in 1:(nt-1)) ?{for(k in (c+1):nt) ????{or[c, k] < -exp(d[k]-d[c]) ????lor[c, k] < -(d[k]-d[c]) ????} }
處理完畢后,將其保存為:“.bug”或“.jags”格式。本例以“networkmodel”為名、以“.jags”格式保存在桌面的Rwork文件夾下。
根據程序包的工作原理,其最終迭代抽樣運算仍是使用JAGS軟件來實踐完成的。因此本處還需在R中對模塊進行加載,具體命令如下:
bugsname < -file.path("C:/Users/Administrator/Desktop/Rwork", "networkmodel.jags")
2.2 數據及初始值加載
本例3個程序包雖均采取《R軟件R2WinBUGS程序包在網狀Meta分析中的應用》一文[5]中的數據和初始值,但相關的格式卻有所不同:
對于R2jags及rjags程序包而言,本例數據及初始值的排列格式及加載方式與上文中步驟2.2~2.5 [5]一致,本處不再講解。
對于runjags程序包,本例數據及初始值的排列格式及加載方式與上文中步驟2.2~2.4 [5]一致;但還需對數據及初始值格式作相關整合,具體命令分別如下:
data < -dump.format(list(t=t, r=r, n=n, na=na, ns=ns, nt=nt))
inits < -dump.format(list(d=c(NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), sd=1, mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)))
2.3 觀察指標設定
由于在上述模塊編譯中,JAGS軟件對lor和or數據集無法進行識別與監視,因此,此處需要將lor和or數據集中每個節點一一寫入監視。本例為了節省計算機內存資源以及方便后續結果展示,僅監視lor[1, 2],or[1, 2],lor[1, 3],or[1, 3]節點作為示范,監視代碼設定為:
parameters < -c("lor[1, 2]", "or[1, 2]", "lor[1, 3]", "or[1, 3]")。
2.4 迭代運算
此3個程序包中,雖在迭代運算基本思路與原理是相同的,但其迭代命令代碼仍有所不同。表 2及圖 1展示了這3個程序包迭代命令及迭代進程。


2.5 結果顯示及匯總
在上述執行迭代命令運行完成后,我們可以對結果進行顯示及匯總了。具體命令代碼見表 3,相關匯總結果見圖 2~圖 5。





3 結果圖形繪制
如前所述,圖形的繪制是JAGS軟件的短板,R則可在這方面給予其很好的補充,且實現方法多樣。本例將提供3種繪制結果圖的方法。
3.1 方法一
具體操作命令見表 4,執行命令后即可繪制相關的圖形(本處僅展示R2jags程序包繪制的森林圖,其余圖略)。

值得注意的是:盡管3個程序包雖均能繪制結果圖,但仍需相互結合使用方能繪制出森林圖、軌跡圖及密度圖。

3.2 方法二
首先,提取各監視節點的軌跡數據,再依據該數據進行軌跡圖與密度圖的繪制。提取各監視節點的軌跡數據命令見表 5。

提取監視節點迭代軌跡數據后,即可依據該數據繪制軌跡圖與密度圖,具體操作命令見表 6,結果圖略。

3.3 方法三
通過R軟件的coda程序包中的codamenu()菜單型命令進行結果圖形的繪制,首先需要加載coda程序包,命令如下:
install.packages(“coda”)
library(“coda”)
加載完成,需要進行結果的讀取,需要使用coda程序包中codamenu()菜單型命令,具體命令如下:
codamenu()
在執行上述命令后,將會出現以下3個菜單提供選擇:
1: Read BUGS output files
2: Use an mcmc object
3: Quit
Selection:
上述3個選擇的含義分別為:1,讀取BUGS軟件生成文件;2,使用mcmc函數對象;3,退出。在下方自動出現的“Selection:”后填寫2,來讀取“.mcmc”后綴結果數據。執行“Selection: 2”命令后,出現以下提示字樣:
Enter name of saved object (or type "exit" to quit)
1:
在“1:”的后方鍵入“.mcmc”后綴格式的結果數據變量集(即分別鍵入jagsfit.mcmc、coda.mcmc和run.jags.mcmc)并執行上述命令后,軟件會自動提示以下操作主菜單(標記為1),具體如下:
CODA Main Menu
1: Output Analysis
2: Diagnostics
3: List/Change Options
4: Quit
Selection:
在主菜單的目錄中,含義分別為:1,為結果分析菜單(標記為1.1);2,為結果診斷菜單(標記為1.2);3,對應的為該節點的列表信息,具體使用價值不大(有興趣讀者可自行嘗試);4,退出。此時,在“Selection:”后填寫1,將會出現子級菜單[即Output Analysis為結果分析菜單(標記為1.1)]供選擇,具體如下:
CODA Output Analysis menu
1: Plots
2: Statistics
3: List/Change Options
4: Return to Main Menu
Selection:
其中,選擇不同選項產生的結果不同:1,表示該節點軌跡圖與密度圖的繪制;2,表示該節點的相關信息匯總,其包含有Mean、SD、Na?ve SE、Time-series SE以及2.5%、25%、50%、75%、97.5%等分位數(見圖 7);3,對應的為該節點的列表信息;4,為返回主菜單。在“Selection:”后填寫2,將會出現下一子級菜單[即Diagnostics為結果診斷菜單(標記為1.2)]供選擇,具體如下:

CODA Diagnostics Menu
1: Geweke
2: Gelman and Rubin
3: Raftery and Lewis
4: Heidelberger and Welch
5: Autocorrelations
6: Cross-Correlations
7: List/Change Options
8: Return to Main Menu
Selection:
該菜單中:1~4均為各種收斂診斷方法;5為自相關性診斷;6為交互相關性;7為自行嘗試;8為返回主菜單。在“Selection:”處填寫5,即可顯示其自相關性的結果(見圖 8)并出現以下子菜單:

Autocorrelation Plots Menu
1: Plot autocorrelations
2: Return to Diagnostics Menu
Selection:
在其子菜單中(標記為1.2.5)下,若在“Selection:”處繼續選擇1時,即可繪制自相關性圖形。上述菜單中,還可繼續嵌套相關的子菜單,在此就不一一介紹,有興趣的讀者可自行嘗試。
4 結語
通過上述介紹,我們發現JAGS軟件與R的結合,使JAGS軟件生成的結果文件能在R軟件中得到完美的解讀,無論是結果準確性的解讀還是圖形精美繪制都可以以菜單命令形式展現,這也是在編程軟件中少見的,對于編程基礎語言欠缺的操作者來說,十分便捷。
然而,由于此3個程序包均是基于JAGS軟件來實現的,所以其共同的缺陷在于只能監視單一節點,而不能監視節點變量。與調用WinBUGS和OpenBUGS的程序相比[5, 7, 8],在節點較多的分析中,使用繁瑣;此外,這3個程序包以及JAGS軟件本身均缺乏rank函數功能[1],故無法直接進行排序,需要通過其他方式或運行其他代碼來實現。
從繪圖方面而言,本文展示的3種結果圖繪制方法中前兩種的方法操作相對較為簡單,但圖形的繪制需要相互補充;第三種方法功能較為全面,其不僅具有繪制多種圖形的能力,同時還能對其結果進行匯總并給與相關分析。故筆者建議使用第三種方法來完成相關操作,具體選取可依使用者需求而定。
細心的讀者會發現,在引言中筆者還提及了bayesmix程序包。然而,bayesmix程序包適用于單變量的混合模型數據且該程序包的model依據數據及參數的簡單設定在后臺自動生成,無法加載網狀model;再者,bayesmix程序包的實質還是運用rjags程序包來實現對JAGS軟件的調用。因此,本文未對其進行介紹。
總之,R2jags、rjags及runjags程序包均依附于R語言框架下,通過調用JAGS軟件的方式來實現JAGS軟件在網狀Meta分析中的運用。其不僅減少了JAGS軟件的一些不必要的繁瑣操作,而且彌補了該軟件在結果讀取、分析、匯總以及圖形繪制等方面的功能缺失;其為操作者在基于貝葉斯理論軟件的選取上又提供了新的友好選擇。