傳統 Meta 分析的統計基礎是頻率學理論，這種 Meta 分析通常采用 Q 統計量衡量不同研究之間的異質性大小，然后根據異質性檢驗結果來確定采用何種模型（固定效應模型或隨機效應模型）進行 Meta 分析。經典頻率學隨機效應模型可通過 Q 檢驗獲得研究間方差的矩估計，但是無法獲得研究間協方差的 95% 可信區間（confidence interval，CI），而且基于頻率學的 Meta 分析檢驗效能較低。在傳統的 Meta 分析模型中，效應量指標的分布一般假設服從正態分布。例如，對于二分類資料的 Meta 分析而言，最常用的效應量為比值比（odds ratio，OR）^[1]。經典 Meta 分析方法是將 OR 對數化，然后將 log（OR）設為效應值，并假設其服從正態分布^[2]。這種 Meta 分析模型是以正態分布假設為前提的，因此在存在小樣本資料時，如果不符合近似正態的條件，經典方法合并效應值存在一定誤差^[3]。此外，經典 Meta 分析模型很難勝任復雜的統計模型，并且無法識別和處理極端值^[4]。

貝葉斯 Meta 分析（Bayesian meta-analysis，BMA）是近年來以貝葉斯統計理論為基礎發展起來的一種新型的 Meta 分析方法，主要采用“馬爾科夫鏈-蒙特卡羅”（Markov chain Monte Carlo，MCMC）方法進行 Meta 分析^[5]。貝葉斯方法與頻率論（經典）方法的原理完全不同，頻率論認為概率反映的是某事件在重復試驗中出現的次數與試驗次數之比，也就是事件發生的頻率，具有隨機性和穩定性特點，這是頻率論的統計基礎；而貝葉斯方法認為概率是數值的可信程度。這意味著，頻率論將模型參數設置為固定，而將數據設置為隨機；貝葉斯則與之相反^[6]。貝葉斯統計學與傳統頻率學在統計推斷的理念和方法上存在本質上的區別，包括對信息的利用、對未知參數 θ 的解釋及統計推斷理念上的差異等^{[4, 7]}。

Meta 分析同樣可采用貝葉斯分析實現。貝葉斯統計將參數 θ 作為一個隨機變量，有一定的先驗分布，在獲得樣本之后（給定的樣本信息），θ 的后驗分布 π（θ|x）應包含 θ 的綜合信息，可從 θ 的后驗分布獲得參數 θ 的貝葉斯估計。貝葉斯分析能很好地解決經典 Meta 分析存在的缺陷，近年來，采用貝葉斯方法的 Meta 分析引起了廣泛重視，基于貝葉斯的網絡 Meta 分析是 Meta 分析領域的研究熱點。貝葉斯 Meta 分析最常用的統計軟件是 WinBUGS、R 軟件，但實際上 SAS 提供了兩種過程來實現貝葉斯 Meta 分析，分別是 PROC GENMOD 與 PROC MCMC^[8]。本研究團隊先前已介紹了利用 PROC MCMC 實現診斷性試驗的貝葉斯 Meta 分析方法^[6]。然而，目前尚無學者介紹利用 SAS 軟件實現二分類數據貝葉斯 Meta 分析的方法。本文通過實例數據介紹利用 SAS 軟件中的 PROC MCMC 實現二分類數據貝葉斯 Meta 分析，同時將結果與基于頻率學 Meta 分析方法的結果進行對比，相關的 SAS 代碼在文后附錄中給出。

1   貝葉斯 Meta 分析理論簡介

1.1   貝葉斯統計理論簡介

貝葉斯分析是基于貝葉斯定理發展起來用于系統闡述和解決統計問題的方法^[3]。其理論簡述如下：根據既往經驗，假設某事件 θ 發生的概率為 P（θ），我們這個概率稱為先驗概率。然后我們現有一組新的數據 y，則 y 在 y/θ 的前提下發生的條件概率記為 P（y/θ），此條件概率稱為似然。根據先驗概率和似然可計算出概率 P（θ/y），表示 θ 在 y 存在的前提下發生事件的可能性大小，稱為后驗概率。后驗概率與先驗概率與似然的乘積成正比，即 P（θ/y）=α×P（y/θ）×P（θ）。在貝葉斯框架下，參數 θ 是一個隨機變量，服從一定的先驗分布，在給定的樣本信息（抽樣）條件下，通過抽樣獲取樣本，則 θ 的后驗分布 應包含 θ 的綜合信息，因此可通過 θ 的后驗分布來獲得參數 θ 的貝葉斯估計。因此，貝葉斯統計是由模型、參數和似然 3 個要素構成的統計方法。

1.2   二分類數據貝葉斯 Meta 分析的建模方法

二分類數據 Meta 分析最常用的效應指標是 OR，經過對數轉換后的 OR 值服從正態分布。與頻數法相對應，本文貝葉斯 Meta 分析的建模方法同樣基于 OR 的對數轉換值。具體模型構建如下：

假設某 Meta 分析共納入 n 個研究，nt、rt、nc、rc 分別為治療組和對照組總人數及事件發生例數。 為第 i 個研究的效應量，本文為 OR 的對數轉換值， 為第 i 個研究的效應量的研究內方差，theta 為研究的真實效應值。上述參數建模如下：

固定效應模型假設所有納入 Meta 分析的研究來自同一個總體，因此效應值固定。對于固定效應模型，，其中 是固定效應模型背景下的真實效應值，是待估計的固定效應模型參數。因此，貝葉斯固定效應模型可建模如下：

，需要獲得的數據集包括：；。模型待估計的參數為 。似然函數為正態分布密度 。后驗分布為 。

但由于納入 Meta 分析的原始研究通常存在異質性，導致使用固定效應模型的假設條件往往不成立。這時采用隨機效應模型，即假設納入 Meta 分析的研究可能來自不同的研究總體，每一個研究都擁有不同的效應值，這些效應值服從同一分布的總體。對于隨機效應模型，，而，方差 反映了研究間變異的度量，代表隨機效應模型的隨機部分。貝葉斯隨機效應模型可建模如下：

，。需要獲得的數據集包括：；。模型待估計的參數包括 ， 和 。其中 表示每個研究的效應， 為真實效應值； 為研究的方差。似然函數與固定效應模型相似。后驗分布為 。

2   資料與方法

2.1   示例數據介紹

該示例數據來源于陳艷麗等^[9]發表的 Meta 分析，該 Meta 分析旨在評價腹腔鏡下根治性全子宮切除對比開腹手術對早期宮頸癌的療效和預后的影響。本研究提取該文中并發癥的發生率數據（表 1），采用貝葉斯 Meta 分析比較兩種手術方式對術后并發癥發生率的影響。

將表 1 的數據錄入 SAS 軟件，并存儲在名為 Dichotomous_data 的數據集中，數據集錄入的代碼見附錄。為擬合本文的 BMA 模型，需要的數據集變量包括 author、year、rt、Nt、rc、Nc；分別代表作者、發表年份、研究組事件數（在本示例中，事件代表術后發生并發癥）、研究組樣本量、對照組事件數和對照組樣本量。

2.2   SAS 代碼介紹

利用附錄代碼 1 完成原始數據集的錄入之后，再利用附錄代碼 2 計算效應值 di 及標準誤 sigmai2。需要注意的是，本文示例數據不存在零事件試驗，對于零事件試驗，需要進行連續性校正，校正方法是 2×2 四格表每個單元格增加 0.5；然后使用 PROC MCMC 程序，分別擬合固定效應模型和隨機效應模型的貝葉斯 Meta 分析。

在固定效應模型 BMA 中，模型將迭代 50 000 次，前 5 000 次退火用于消除初始值的影響。固定效應值 θ 的先驗分布被設置為均數為 0，方差為 1×10⁶的正態分布。不同研究觀測的效應值 di 被建模為模型均數為 θ，方差為 1×10⁶的正態分布。與固定效應模型不同的是，隨機效應模型引入一個反映研究間變異的度量 tau2，并且 tau2 的先驗分布被假設服從形狀參數（shape parameter）為 0.001、逆尺度參數為 0.001 的伽瑪分布（Gamma distribution）。在隨機效應模型背景下，通過“random deltai ～ normal（mean=theta，var=tau2）subject=study monitor=（deltai）；”語句設置隨機效應。

為保證基于 SAS 軟件的貝葉斯 Meta 分析具有實用性，本文給出了 SAS 軟件背景下的森林圖繪制代碼^[10]。森林圖繪制代碼將筆者先前發表的 SAS 代碼包裝成兩個宏程序，分別為&forestplot 和&forestdata。前者用于繪制森林圖，其宏參數是需要繪制森林圖的 SAS 數據集 Forest；后者用于構建森林圖數據集，并調用宏程序 forestplot 實現森林圖繪制；其宏參數是 PROC MCMC 過程步完成數據分析后輸出的后驗摘要數據集 PostSummaries_data 和指定的 Meta 分析模型 type（type=1 表示固定效應模型，type=2 表示隨機效應模型）。

3   結果

3.1   固定效應模型貝葉斯 Meta 分析結果

固定效應模型貝葉斯 Meta 分析結果見表 2，表中的 theta 和 OR 是對數轉換關系，即 OR=exp（theta）。表 3 是 PROC MCMC 提供的參數值的 95% 可信區間。圖 1 是 MCMC 采樣分析圖，分別包括了 Gibbs 抽樣動態蹤跡圖、自相關圖和后驗分布核密度估計圖。從圖 1 中的 Gibbs 動態圖可看出，參數 OR 經過 50 000 次迭代后基本收斂。在自相關圖中，OR 的驗后自相關非常小，模型可靠性高。核密度圖顯示，OR 的后驗分布類似于正態分布，與研究預期相似。

3.2   隨機效應模型貝葉斯 Meta 分析結果

表 4 給出了隨機效應模型貝葉斯 Meta 分析結果。除了給出 θ 的估計值外，隨機效應模型還給出了各個研究的隨機效應值。表 4 中，除了 OR 值外，其他需要經過擬對數轉換后應用于分析。隨機效應模型的 MCMC 采樣分析圖見圖 2，除了圖 2，隨機效應模型同時還會對其他參數進行 MCMC 采樣分析，讀者可自行運行本文的 SAS 代碼查看。

3.3   不同 Meta 分析方法的結果比較

表 5 給出了頻率學方法和貝葉斯方法的 Meta 分析結果，兩種方法給出的模型參數估計值極其相似。

3.4   貝葉斯 Meta 分析的森林圖結果

森林圖可非常直觀地觀察到 Meta 分析的結果。但 PROC MCMC 過程步無法直接給出森林圖，為便于讀者利用 SAS 軟件開展研究，本文開發了 SAS 軟件繪制森林圖的代碼（見附錄），運行代碼后，可獲得固定效應模型和隨機效應模型的森林圖。圖 3 顯示了固定效應模型的森林圖。

4   討論

貝葉斯統計和經典統計學派（頻率學派）是目前兩大主要的統計學派。貝葉斯統計比頻率學統計方法需要更強大的計算能力來完成推斷，因此貝葉斯統計是以計算機運算能力發展為前提的。隨著計算機技術的發展，貝葉斯統計學在醫學領域內的應用越來越廣泛，其中也包括 Meta 分析。貝葉斯統計是綜合未知參數的先驗信息和樣本信息，根據貝葉斯定理，求出后驗分布，然后根據后驗分布對未知參數進行統計推斷，最后獲得參數估計值的統計方法。在方法學上，頻率統計學方法以樣本為基礎，通過假設檢驗和統計計算的方法進行統計推斷，最后獲得結論；而貝葉斯方法則是基于貝葉斯原理，根據先驗概率推斷后驗概率^[11]。雖然傳統 Meta 分析通常采用頻率學方法實現，但越來越多的學者在探索利用貝葉斯方法實現 Meta 分析的方法。目前已有不少統計軟件可實現貝葉斯 Meta 分析，包括 R、Stata、WinBugs 等眾多軟件^{[4, 8, 12-16]}。Winbugs 是用于 Gibbs 抽樣的專用軟件包，是免費軟件，目前已廣泛應用于實施貝葉斯方法^[4]。既往學者研究表明，與頻率學 Meta 分析相比，貝葉斯 Meta 分析具有以下優點：① 建模方式更靈活，在多種貝葉斯分析軟件中可操作性強；② 貝葉斯 Meta 分析充分考慮了模型參數（如方差等）的不確定性，例如可估計研究間方差的不精確性，最終獲得的合并效應量點估計及 95% 可信區間來源于參數的后驗分布，結果可靠；③ 貝葉斯 Meta 分析提供的統計檢驗，反映在某種程度上可改變人們的觀念^[17]。此外，貝葉斯分析方法最大的優勢是不受樣本量限制，因為貝葉斯的基本原理就是通過先驗概率推斷后驗概率；而頻率學派反映的是抽樣概率，因此需要樣本量的支持。而實踐中往往納入 Meta 分析的原始研究數量極為有限，這為貝葉斯統計提供了更廣闊的應用范圍。貝葉斯 Meta 分析不足之處有以下兩點：① 頻率學派認為一個事件的概率可通過大量重復試驗下事件出現的頻率來解釋，這種概率不依賴于主觀性，被認為是客觀概率，而貝葉斯統計需要指定先驗分布，而這種分布帶有主觀性，因此被稱為主觀概率；因此，其參數的先驗分布受統計分析師建模思想的影響；② 盡管貝葉斯統計并不接受頻率觀點，但貝葉斯統計仍然需要依賴樣本分布，而樣本則來源于抽樣。在實際應用中，往往很難完全區分頻率學派和貝葉斯統計孰優孰劣。在 Meta 分析領域，初學者往往難以把握統計建模，應以頻率學統計結果為準；而具有統計建模基礎的系統評價員，則應更多嘗試貝葉斯統計方法。

雖然 SAS 軟件在統計學領域具有廣泛的應用基礎，但基于 SAS 軟件的貝葉斯 Meta 分析卻起步較晚。目前介紹 SAS 軟件實現貝葉斯 Meta 分析的文獻并不多。在 SAS 9.2 版本之前，SAS 無法實現貝葉斯模型分析，但在 SAS 9.2 之后，SAS 軟件提供了多種過程步來實現貝葉斯 Meta 分析，最常用的過程步包括 PROC GENMOD 和 PROC MCMC。PROC GENMOD 可通過“BAYES”語句實現貝葉斯 Meta 分析；但這種貝葉斯分析與 PROC MCMC 的貝葉斯分析基礎不同，前者類似于頻率學方法，后者是基于 MCMC 的統計方法。因此，在 SAS 中，貝葉斯 Meta 分析首選 PROC MCMC^[18]。國內外已有不少學者介紹了基于 PROC MCMC 實現不同數據的貝葉斯 Meta 分析^{[6, 8, 19]}；但目前尚無介紹 PROC MCMC 實現二分類數據 Meta 分析的文獻。鑒于二分類數據在醫學領域極為常見，有必要介紹二分類數據的貝葉斯 Meta 分析。曾平等^[20]詳細介紹了二分類數據的貝葉斯 Meta 分析的建模方法，本文的側重介紹該模型在 SAS 軟件中的實現過程。

本研究詳細介紹了利用 PROC MCMC 實現二分類數據貝葉斯 Meta 分析的 SAS 程序。通過示例數據分析顯示，貝葉斯 Meta 分析幾乎提供了與頻率學 Meta 分析完全一致的結果。但需要注意幾點：① PROC MCMC 是以馬爾科夫鏈-蒙特卡羅方法為基礎的。MCMC 是在貝葉斯理論框架下，通過計算機進行模擬的蒙特卡洛方法（Monte Carlo）。該方法將馬爾科夫（Markov）過程引入到 Monte Carlo 模擬中，實現抽樣分布隨模擬的進行而改變的動態模擬，彌補了傳統的蒙特卡羅積分只能靜態模擬的缺陷。② 在多數情況下，Metropolis-Hastings 抽樣（M-H sampling）是 MCMC 常用的算法；但對于高維的情形，M-H 算法效率太低。Meta 分析正是一種特殊的高維形式，因此 Gibbs 抽樣應用更廣泛。③ PROC MCMC 提供了兩種不同的后驗可信區間（posterior intervals），分別為等尾可信區間（equal-tail intervals，ET）和最高后驗密度可信區間（highest posterior density intervals，HPD）。讀者在使用可信區間時，應該首選 ET 可信區間，因為 ET 可信區間更窄，而 HPD 可信區間更寬。④ 貝葉斯 Meta 分析不提供統計學概率值（P值），因此不能利用P值來判斷效應值是否具有統計學意義，一般來說，貝葉斯通過后驗分布獲得的可信區間來反映效應值的統計學意義。以 OR 指標為例，根據既往統計學理論，OR 效應值是否具有統計學意義以是否跨越“1”為標準，當可信區間包含了“1”，就認為效應值無統計學意義。因此，結合本文的示例數據，隨機效應模型貝葉斯 Meta 分析獲得的 OR 值為 0.53，95% 可信區間為（0.30，0.85），該可信區間并未包含 1，因此具有統計學意義。其臨床意義為：與開腹手術相比，腹腔鏡下宮頸癌根治術的術后并發癥發生率明顯降低，差異具有統計學意義。最后，與目前多數貝葉斯統計軟件相比，SAS 軟件編程功能強大。基于 SAS 的編程基礎，本研究開發了一段 SAS 宏程序，可非常方便地實現森林圖繪制，這有效解決了目前多數貝葉斯統計軟件無法獲得森林圖的問題。

總之，基于 SAS 軟件強大的編程能力，PROC MCMC 可非常方便地實現二分類數據的貝葉斯 Meta 分析，隨著貝葉斯統計理論的快速發展，貝葉斯 Meta 分析在 Meta 分析領域將發揮越來越重要的作用。

附錄代碼

/*附錄代碼 1：數據集建立及資料錄入*/

Data Dichotomous_data；

　length author $20.；

　study=_N_；

　input author year rtNtrcNc；

　cards；

Lee2002430930

……

；

Run；

/*Proc print data=Dichotomous_data noobs；run；*/

/*附錄代碼 2：計算效應值 di 及標準誤 sigmai2*/

Data Bayesian_Dichotomous_data；

　set Dichotomous_data；

　Di=log（（rt/（Nt-rt））/（rc/（Nc-rc）））；

　sigmai2=1/rt+1/（Nt-rt）+1/rc+1/（Nc-rc）；

run；

proc print data=Bayesian_Dichotomous_data noobs；run；

/*附錄代碼 3：固定效應模型貝葉斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_fixed；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=1234

monitor=（theta OR）stats（percentage=（2.5 97.5））=all dic；

parms theta 0；

prior theta ～ normal（mean = 0，var = 1.0e6）；

model di ～ normal（theta，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_fixed；run；

/*附錄代碼 4：森林圖繪制代碼*/

%macro forestplot（forestdata）；

data forest；

　　set & forestdata.；

　　format OddsRatio 6.2；

　　format Lowercl 6.2；format Uppercl 6.2；

　　format lcl2 6.2；format ucl2 6.2；

　　format SE 6.2；format weight percent5.；

　　　format Q1 Q3 4.2；

　if studyname="" then do；grp=2；study2="Overall"；end；

　else do；grp=1；end；

　　ObsId=_N_；

　　if mod（obsid，2）= 0 then studyref=studyname；

　　sum_weight=& sum_weight.；

　　weight=weight0/sum_weight；

　　lcl2=lowercl；

　　ucl2=uppercl；

　if grp=1 then do；

　　weightQ=weight0/sum_weight*0.5；

　　Q1=OddsRatio-OddsRatio*weightQ；

　　Q3=OddsRatio+OddsRatio*weightQ；

　　study2=""；

　　end；

　　OR='OR'；LCL='可信區間低限'；UCL='可信區間上限'；WT='權重'；

run；

data _null_；

　pct=0.75/nobs；height=3；dpi=100；

　call symputx（"pct"，pct）；

　call symputx（"pct2"，2*pct）；

　call symputx（"dpi"，dpi）；

　call symputx（"height"，height）；

　call symputx（"heightin"，height || "in"）；

　call symputx（"thickness"，floor（height*dpi*pct））；

　set forest nobs=nobs；

run；

proc template；

　define style styles.ForestColor93；

　　parent = Styles.htmlBlue；

　　　style GraphFonts from GraphFonts /

　　　　'GraphDataFont' =（"<sans-serif>，<MTsans-serif>"，7pt）

　　　　'GraphValueFont' =（"<sans-serif>，<MTsans-serif>"，7pt）；

　　　style GraphData1 from GraphData1 /

　　　　contrastcolor = GraphColors（'gcdata2'）

　　　　color = GraphColors（'gdata2'）；

　　　style GraphData2 from GraphData2 /

　　　　contrastcolor = GraphColors（'gcdata1'）

　　　　color = GraphColors（'gdata1'）；

　　end；

run；

/*設置森林圖樣式，DPI，圖像參數和標題等*/

ods graphics / reset width=6in height=& heightin imagename="ForestPlotColor" IMAGEFMT=EMF；

ods listing sge=off style=Styles.ForestColor93 image_dpi=& dpi；

title h=12pt "二分類數據的貝葉斯 Meta 分析的森林圖"；

title2 h=8pt '比值比 OR 及其 95% 可信區間'；

proc sgplot data=forest noautolegend nocycleattrs；

　refline studyref / lineattrs=（thickness=& thickness）transparency=0.85；

　scatter y=study2 x=oddsratio / markerattrs=（symbol=diamondfilled size=10）；

　highlow y=studyname low=lcl2 high=ucl2 / type=line；

　highlow y=studyname low=q1 high=q3 / type=bar；

　scatter y=studyname x=or / markerchar=oddsratio x2axis；

　scatter y=studyname x=lcl / markerchar=lowercl x2axis；

　scatter y=studyname x=ucl / markerchar=uppercl x2axis；

　scatter y=studyname x=wt / markerchar=weight x2axis；

　scatter y=study2 x=or / markerchar=oddsratio x2axis；

　scatter y=study2 x=lcl / markerchar=lowercl x2axis；

　scatter y=study2 x=ucl / markerchar=uppercl x2axis；

　scatter y=study2 x=wt / markerchar=weight x2axis；

　refline 1 100 / axis=x；

　refline 0.01 0.1 10 / axis=x lineattrs=（pattern=shortdash）transparency=0.5；

　inset '開腹手術組' / position=bottom；

　inset '腹腔鏡下手術組'/ position=bottomleft；

　xaxis type=log offsetmin=0 offsetmax=0.35 min=0.01 max=100 minor display=（nolabel）；

　x2axis offsetmin=0.7 display=（noticks nolabel）；

　yaxis display=（noticks nolabel）offsetmin=& pct offsetmax=& pct2 reverse；

run；

Title；

%mend；

%macro forestdata（PostSummaries_data，type）；

%if & type=1 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 2；

　run；

　%let tau=0；

%end；

%if & type=2 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 3；

　run；

　data temp_tau；

　　set & PostSummaries_data.；

　　if _N_= 2；

　　keep mean；

　run；

　proc sql noprint；

　　select mean into：tau0 from temp_tau；

　quit；

　%put & tau0.；

　%let tau=& tau0.；

%end；

data tem1；

　　set temp0；

　　keep mean stddev P2_5 P97_5；

　　rename mean=OddsRatio；rename StdDev=SE；

　　rename P2_5=Lowercl；rename P97_5=Uppercl；

run；

data Forest_plot_data；

　set Bayesian_Dichotomous_data；

/*計算效應值 OR 及標準誤*/

　OddsRatio=exp（di）；

　Lowercl=round（exp（di-1.96*sqrt（sigmai2）），.01）；

　Uppercl=round（exp（di+1.96*sqrt（sigmai2）），.01）；

　SE=（Uppercl-Lowercl）/（2*1.96）；

　studyname=author；

　weight0=1/（sigmai2+& tau.）；

　keep studyname OddsRatio SE Lowercl Uppercl weight0；

run；

proc sql noprint；

　select sum（weight0）into：sum_weight from Forest_plot_data；

quit；

　%put & sum_weight.；

data forest；

　set Forest_plot_data tem1；

run；

%forestplot（forest）；

%mend；

%forestdata（mcmc_Summaries_fixed，1）；

/*附錄代碼 5：隨機效應模型貝葉斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_random；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=246810

monitor=（theta tau2 OR）stats（percentage=（2.5 97.5））=all dic diagnostics=all；

parms theta=0 tau2 1；

prior theta ～ normal（mean = 0，var = 1e6）；

prior tau2 ～ igamma（0.001，scale =0.001）；

random deltai ～ normal（mean = theta，var = tau2）subject=study monitor=（deltai）；

model di ～ normal（deltai，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_random；run；

%forestdata（mcmc_Summaries_random，2）；

傳統 Meta 分析的統計基礎是頻率學理論，這種 Meta 分析通常采用 Q 統計量衡量不同研究之間的異質性大小，然后根據異質性檢驗結果來確定采用何種模型（固定效應模型或隨機效應模型）進行 Meta 分析。經典頻率學隨機效應模型可通過 Q 檢驗獲得研究間方差的矩估計，但是無法獲得研究間協方差的 95% 可信區間（confidence interval，CI），而且基于頻率學的 Meta 分析檢驗效能較低。在傳統的 Meta 分析模型中，效應量指標的分布一般假設服從正態分布。例如，對于二分類資料的 Meta 分析而言，最常用的效應量為比值比（odds ratio，OR）^[1]。經典 Meta 分析方法是將 OR 對數化，然后將 log（OR）設為效應值，并假設其服從正態分布^[2]。這種 Meta 分析模型是以正態分布假設為前提的，因此在存在小樣本資料時，如果不符合近似正態的條件，經典方法合并效應值存在一定誤差^[3]。此外，經典 Meta 分析模型很難勝任復雜的統計模型，并且無法識別和處理極端值^[4]。

貝葉斯 Meta 分析（Bayesian meta-analysis，BMA）是近年來以貝葉斯統計理論為基礎發展起來的一種新型的 Meta 分析方法，主要采用“馬爾科夫鏈-蒙特卡羅”（Markov chain Monte Carlo，MCMC）方法進行 Meta 分析^[5]。貝葉斯方法與頻率論（經典）方法的原理完全不同，頻率論認為概率反映的是某事件在重復試驗中出現的次數與試驗次數之比，也就是事件發生的頻率，具有隨機性和穩定性特點，這是頻率論的統計基礎；而貝葉斯方法認為概率是數值的可信程度。這意味著，頻率論將模型參數設置為固定，而將數據設置為隨機；貝葉斯則與之相反^[6]。貝葉斯統計學與傳統頻率學在統計推斷的理念和方法上存在本質上的區別，包括對信息的利用、對未知參數 θ 的解釋及統計推斷理念上的差異等^{[4, 7]}。

Meta 分析同樣可采用貝葉斯分析實現。貝葉斯統計將參數 θ 作為一個隨機變量，有一定的先驗分布，在獲得樣本之后（給定的樣本信息），θ 的后驗分布 π（θ|x）應包含 θ 的綜合信息，可從 θ 的后驗分布獲得參數 θ 的貝葉斯估計。貝葉斯分析能很好地解決經典 Meta 分析存在的缺陷，近年來，采用貝葉斯方法的 Meta 分析引起了廣泛重視，基于貝葉斯的網絡 Meta 分析是 Meta 分析領域的研究熱點。貝葉斯 Meta 分析最常用的統計軟件是 WinBUGS、R 軟件，但實際上 SAS 提供了兩種過程來實現貝葉斯 Meta 分析，分別是 PROC GENMOD 與 PROC MCMC^[8]。本研究團隊先前已介紹了利用 PROC MCMC 實現診斷性試驗的貝葉斯 Meta 分析方法^[6]。然而，目前尚無學者介紹利用 SAS 軟件實現二分類數據貝葉斯 Meta 分析的方法。本文通過實例數據介紹利用 SAS 軟件中的 PROC MCMC 實現二分類數據貝葉斯 Meta 分析，同時將結果與基于頻率學 Meta 分析方法的結果進行對比，相關的 SAS 代碼在文后附錄中給出。

1   貝葉斯 Meta 分析理論簡介

1.1   貝葉斯統計理論簡介

貝葉斯分析是基于貝葉斯定理發展起來用于系統闡述和解決統計問題的方法^[3]。其理論簡述如下：根據既往經驗，假設某事件 θ 發生的概率為 P（θ），我們這個概率稱為先驗概率。然后我們現有一組新的數據 y，則 y 在 y/θ 的前提下發生的條件概率記為 P（y/θ），此條件概率稱為似然。根據先驗概率和似然可計算出概率 P（θ/y），表示 θ 在 y 存在的前提下發生事件的可能性大小，稱為后驗概率。后驗概率與先驗概率與似然的乘積成正比，即 P（θ/y）=α×P（y/θ）×P（θ）。在貝葉斯框架下，參數 θ 是一個隨機變量，服從一定的先驗分布，在給定的樣本信息（抽樣）條件下，通過抽樣獲取樣本，則 θ 的后驗分布 應包含 θ 的綜合信息，因此可通過 θ 的后驗分布來獲得參數 θ 的貝葉斯估計。因此，貝葉斯統計是由模型、參數和似然 3 個要素構成的統計方法。

1.2   二分類數據貝葉斯 Meta 分析的建模方法

二分類數據 Meta 分析最常用的效應指標是 OR，經過對數轉換后的 OR 值服從正態分布。與頻數法相對應，本文貝葉斯 Meta 分析的建模方法同樣基于 OR 的對數轉換值。具體模型構建如下：

假設某 Meta 分析共納入 n 個研究，nt、rt、nc、rc 分別為治療組和對照組總人數及事件發生例數。 為第 i 個研究的效應量，本文為 OR 的對數轉換值， 為第 i 個研究的效應量的研究內方差，theta 為研究的真實效應值。上述參數建模如下：

固定效應模型假設所有納入 Meta 分析的研究來自同一個總體，因此效應值固定。對于固定效應模型，，其中 是固定效應模型背景下的真實效應值，是待估計的固定效應模型參數。因此，貝葉斯固定效應模型可建模如下：

，需要獲得的數據集包括：；。模型待估計的參數為 。似然函數為正態分布密度 。后驗分布為 。

但由于納入 Meta 分析的原始研究通常存在異質性，導致使用固定效應模型的假設條件往往不成立。這時采用隨機效應模型，即假設納入 Meta 分析的研究可能來自不同的研究總體，每一個研究都擁有不同的效應值，這些效應值服從同一分布的總體。對于隨機效應模型，，而，方差 反映了研究間變異的度量，代表隨機效應模型的隨機部分。貝葉斯隨機效應模型可建模如下：

，。需要獲得的數據集包括：；。模型待估計的參數包括 ， 和 。其中 表示每個研究的效應， 為真實效應值； 為研究的方差。似然函數與固定效應模型相似。后驗分布為 。

2   資料與方法

2.1   示例數據介紹

該示例數據來源于陳艷麗等^[9]發表的 Meta 分析，該 Meta 分析旨在評價腹腔鏡下根治性全子宮切除對比開腹手術對早期宮頸癌的療效和預后的影響。本研究提取該文中并發癥的發生率數據（表 1），采用貝葉斯 Meta 分析比較兩種手術方式對術后并發癥發生率的影響。

將表 1 的數據錄入 SAS 軟件，并存儲在名為 Dichotomous_data 的數據集中，數據集錄入的代碼見附錄。為擬合本文的 BMA 模型，需要的數據集變量包括 author、year、rt、Nt、rc、Nc；分別代表作者、發表年份、研究組事件數（在本示例中，事件代表術后發生并發癥）、研究組樣本量、對照組事件數和對照組樣本量。

2.2   SAS 代碼介紹

利用附錄代碼 1 完成原始數據集的錄入之后，再利用附錄代碼 2 計算效應值 di 及標準誤 sigmai2。需要注意的是，本文示例數據不存在零事件試驗，對于零事件試驗，需要進行連續性校正，校正方法是 2×2 四格表每個單元格增加 0.5；然后使用 PROC MCMC 程序，分別擬合固定效應模型和隨機效應模型的貝葉斯 Meta 分析。

在固定效應模型 BMA 中，模型將迭代 50 000 次，前 5 000 次退火用于消除初始值的影響。固定效應值 θ 的先驗分布被設置為均數為 0，方差為 1×10⁶的正態分布。不同研究觀測的效應值 di 被建模為模型均數為 θ，方差為 1×10⁶的正態分布。與固定效應模型不同的是，隨機效應模型引入一個反映研究間變異的度量 tau2，并且 tau2 的先驗分布被假設服從形狀參數（shape parameter）為 0.001、逆尺度參數為 0.001 的伽瑪分布（Gamma distribution）。在隨機效應模型背景下，通過“random deltai ～ normal（mean=theta，var=tau2）subject=study monitor=（deltai）；”語句設置隨機效應。

為保證基于 SAS 軟件的貝葉斯 Meta 分析具有實用性，本文給出了 SAS 軟件背景下的森林圖繪制代碼^[10]。森林圖繪制代碼將筆者先前發表的 SAS 代碼包裝成兩個宏程序，分別為&forestplot 和&forestdata。前者用于繪制森林圖，其宏參數是需要繪制森林圖的 SAS 數據集 Forest；后者用于構建森林圖數據集，并調用宏程序 forestplot 實現森林圖繪制；其宏參數是 PROC MCMC 過程步完成數據分析后輸出的后驗摘要數據集 PostSummaries_data 和指定的 Meta 分析模型 type（type=1 表示固定效應模型，type=2 表示隨機效應模型）。

3   結果

3.1   固定效應模型貝葉斯 Meta 分析結果

固定效應模型貝葉斯 Meta 分析結果見表 2，表中的 theta 和 OR 是對數轉換關系，即 OR=exp（theta）。表 3 是 PROC MCMC 提供的參數值的 95% 可信區間。圖 1 是 MCMC 采樣分析圖，分別包括了 Gibbs 抽樣動態蹤跡圖、自相關圖和后驗分布核密度估計圖。從圖 1 中的 Gibbs 動態圖可看出，參數 OR 經過 50 000 次迭代后基本收斂。在自相關圖中，OR 的驗后自相關非常小，模型可靠性高。核密度圖顯示，OR 的后驗分布類似于正態分布，與研究預期相似。

3.2   隨機效應模型貝葉斯 Meta 分析結果

表 4 給出了隨機效應模型貝葉斯 Meta 分析結果。除了給出 θ 的估計值外，隨機效應模型還給出了各個研究的隨機效應值。表 4 中，除了 OR 值外，其他需要經過擬對數轉換后應用于分析。隨機效應模型的 MCMC 采樣分析圖見圖 2，除了圖 2，隨機效應模型同時還會對其他參數進行 MCMC 采樣分析，讀者可自行運行本文的 SAS 代碼查看。

3.3   不同 Meta 分析方法的結果比較

表 5 給出了頻率學方法和貝葉斯方法的 Meta 分析結果，兩種方法給出的模型參數估計值極其相似。

3.4   貝葉斯 Meta 分析的森林圖結果

森林圖可非常直觀地觀察到 Meta 分析的結果。但 PROC MCMC 過程步無法直接給出森林圖，為便于讀者利用 SAS 軟件開展研究，本文開發了 SAS 軟件繪制森林圖的代碼（見附錄），運行代碼后，可獲得固定效應模型和隨機效應模型的森林圖。圖 3 顯示了固定效應模型的森林圖。

4   討論

貝葉斯統計和經典統計學派（頻率學派）是目前兩大主要的統計學派。貝葉斯統計比頻率學統計方法需要更強大的計算能力來完成推斷，因此貝葉斯統計是以計算機運算能力發展為前提的。隨著計算機技術的發展，貝葉斯統計學在醫學領域內的應用越來越廣泛，其中也包括 Meta 分析。貝葉斯統計是綜合未知參數的先驗信息和樣本信息，根據貝葉斯定理，求出后驗分布，然后根據后驗分布對未知參數進行統計推斷，最后獲得參數估計值的統計方法。在方法學上，頻率統計學方法以樣本為基礎，通過假設檢驗和統計計算的方法進行統計推斷，最后獲得結論；而貝葉斯方法則是基于貝葉斯原理，根據先驗概率推斷后驗概率^[11]。雖然傳統 Meta 分析通常采用頻率學方法實現，但越來越多的學者在探索利用貝葉斯方法實現 Meta 分析的方法。目前已有不少統計軟件可實現貝葉斯 Meta 分析，包括 R、Stata、WinBugs 等眾多軟件^{[4, 8, 12-16]}。Winbugs 是用于 Gibbs 抽樣的專用軟件包，是免費軟件，目前已廣泛應用于實施貝葉斯方法^[4]。既往學者研究表明，與頻率學 Meta 分析相比，貝葉斯 Meta 分析具有以下優點：① 建模方式更靈活，在多種貝葉斯分析軟件中可操作性強；② 貝葉斯 Meta 分析充分考慮了模型參數（如方差等）的不確定性，例如可估計研究間方差的不精確性，最終獲得的合并效應量點估計及 95% 可信區間來源于參數的后驗分布，結果可靠；③ 貝葉斯 Meta 分析提供的統計檢驗，反映在某種程度上可改變人們的觀念^[17]。此外，貝葉斯分析方法最大的優勢是不受樣本量限制，因為貝葉斯的基本原理就是通過先驗概率推斷后驗概率；而頻率學派反映的是抽樣概率，因此需要樣本量的支持。而實踐中往往納入 Meta 分析的原始研究數量極為有限，這為貝葉斯統計提供了更廣闊的應用范圍。貝葉斯 Meta 分析不足之處有以下兩點：① 頻率學派認為一個事件的概率可通過大量重復試驗下事件出現的頻率來解釋，這種概率不依賴于主觀性，被認為是客觀概率，而貝葉斯統計需要指定先驗分布，而這種分布帶有主觀性，因此被稱為主觀概率；因此，其參數的先驗分布受統計分析師建模思想的影響；② 盡管貝葉斯統計并不接受頻率觀點，但貝葉斯統計仍然需要依賴樣本分布，而樣本則來源于抽樣。在實際應用中，往往很難完全區分頻率學派和貝葉斯統計孰優孰劣。在 Meta 分析領域，初學者往往難以把握統計建模，應以頻率學統計結果為準；而具有統計建模基礎的系統評價員，則應更多嘗試貝葉斯統計方法。

雖然 SAS 軟件在統計學領域具有廣泛的應用基礎，但基于 SAS 軟件的貝葉斯 Meta 分析卻起步較晚。目前介紹 SAS 軟件實現貝葉斯 Meta 分析的文獻并不多。在 SAS 9.2 版本之前，SAS 無法實現貝葉斯模型分析，但在 SAS 9.2 之后，SAS 軟件提供了多種過程步來實現貝葉斯 Meta 分析，最常用的過程步包括 PROC GENMOD 和 PROC MCMC。PROC GENMOD 可通過“BAYES”語句實現貝葉斯 Meta 分析；但這種貝葉斯分析與 PROC MCMC 的貝葉斯分析基礎不同，前者類似于頻率學方法，后者是基于 MCMC 的統計方法。因此，在 SAS 中，貝葉斯 Meta 分析首選 PROC MCMC^[18]。國內外已有不少學者介紹了基于 PROC MCMC 實現不同數據的貝葉斯 Meta 分析^{[6, 8, 19]}；但目前尚無介紹 PROC MCMC 實現二分類數據 Meta 分析的文獻。鑒于二分類數據在醫學領域極為常見，有必要介紹二分類數據的貝葉斯 Meta 分析。曾平等^[20]詳細介紹了二分類數據的貝葉斯 Meta 分析的建模方法，本文的側重介紹該模型在 SAS 軟件中的實現過程。

本研究詳細介紹了利用 PROC MCMC 實現二分類數據貝葉斯 Meta 分析的 SAS 程序。通過示例數據分析顯示，貝葉斯 Meta 分析幾乎提供了與頻率學 Meta 分析完全一致的結果。但需要注意幾點：① PROC MCMC 是以馬爾科夫鏈-蒙特卡羅方法為基礎的。MCMC 是在貝葉斯理論框架下，通過計算機進行模擬的蒙特卡洛方法（Monte Carlo）。該方法將馬爾科夫（Markov）過程引入到 Monte Carlo 模擬中，實現抽樣分布隨模擬的進行而改變的動態模擬，彌補了傳統的蒙特卡羅積分只能靜態模擬的缺陷。② 在多數情況下，Metropolis-Hastings 抽樣（M-H sampling）是 MCMC 常用的算法；但對于高維的情形，M-H 算法效率太低。Meta 分析正是一種特殊的高維形式，因此 Gibbs 抽樣應用更廣泛。③ PROC MCMC 提供了兩種不同的后驗可信區間（posterior intervals），分別為等尾可信區間（equal-tail intervals，ET）和最高后驗密度可信區間（highest posterior density intervals，HPD）。讀者在使用可信區間時，應該首選 ET 可信區間，因為 ET 可信區間更窄，而 HPD 可信區間更寬。④ 貝葉斯 Meta 分析不提供統計學概率值（P值），因此不能利用P值來判斷效應值是否具有統計學意義，一般來說，貝葉斯通過后驗分布獲得的可信區間來反映效應值的統計學意義。以 OR 指標為例，根據既往統計學理論，OR 效應值是否具有統計學意義以是否跨越“1”為標準，當可信區間包含了“1”，就認為效應值無統計學意義。因此，結合本文的示例數據，隨機效應模型貝葉斯 Meta 分析獲得的 OR 值為 0.53，95% 可信區間為（0.30，0.85），該可信區間并未包含 1，因此具有統計學意義。其臨床意義為：與開腹手術相比，腹腔鏡下宮頸癌根治術的術后并發癥發生率明顯降低，差異具有統計學意義。最后，與目前多數貝葉斯統計軟件相比，SAS 軟件編程功能強大。基于 SAS 的編程基礎，本研究開發了一段 SAS 宏程序，可非常方便地實現森林圖繪制，這有效解決了目前多數貝葉斯統計軟件無法獲得森林圖的問題。

總之，基于 SAS 軟件強大的編程能力，PROC MCMC 可非常方便地實現二分類數據的貝葉斯 Meta 分析，隨著貝葉斯統計理論的快速發展，貝葉斯 Meta 分析在 Meta 分析領域將發揮越來越重要的作用。

附錄代碼

/*附錄代碼 1：數據集建立及資料錄入*/

Data Dichotomous_data；

　length author $20.；

　study=_N_；

　input author year rtNtrcNc；

　cards；

Lee2002430930

……

；

Run；

/*Proc print data=Dichotomous_data noobs；run；*/

/*附錄代碼 2：計算效應值 di 及標準誤 sigmai2*/

Data Bayesian_Dichotomous_data；

　set Dichotomous_data；

　Di=log（（rt/（Nt-rt））/（rc/（Nc-rc）））；

　sigmai2=1/rt+1/（Nt-rt）+1/rc+1/（Nc-rc）；

run；

proc print data=Bayesian_Dichotomous_data noobs；run；

/*附錄代碼 3：固定效應模型貝葉斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_fixed；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=1234

monitor=（theta OR）stats（percentage=（2.5 97.5））=all dic；

parms theta 0；

prior theta ～ normal（mean = 0，var = 1.0e6）；

model di ～ normal（theta，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_fixed；run；

/*附錄代碼 4：森林圖繪制代碼*/

%macro forestplot（forestdata）；

data forest；

　　set & forestdata.；

　　format OddsRatio 6.2；

　　format Lowercl 6.2；format Uppercl 6.2；

　　format lcl2 6.2；format ucl2 6.2；

　　format SE 6.2；format weight percent5.；

　　　format Q1 Q3 4.2；

　if studyname="" then do；grp=2；study2="Overall"；end；

　else do；grp=1；end；

　　ObsId=_N_；

　　if mod（obsid，2）= 0 then studyref=studyname；

　　sum_weight=& sum_weight.；

　　weight=weight0/sum_weight；

　　lcl2=lowercl；

　　ucl2=uppercl；

　if grp=1 then do；

　　weightQ=weight0/sum_weight*0.5；

　　Q1=OddsRatio-OddsRatio*weightQ；

　　Q3=OddsRatio+OddsRatio*weightQ；

　　study2=""；

　　end；

　　OR='OR'；LCL='可信區間低限'；UCL='可信區間上限'；WT='權重'；

run；

data _null_；

　pct=0.75/nobs；height=3；dpi=100；

　call symputx（"pct"，pct）；

　call symputx（"pct2"，2*pct）；

　call symputx（"dpi"，dpi）；

　call symputx（"height"，height）；

　call symputx（"heightin"，height || "in"）；

　call symputx（"thickness"，floor（height*dpi*pct））；

　set forest nobs=nobs；

run；

proc template；

　define style styles.ForestColor93；

　　parent = Styles.htmlBlue；

　　　style GraphFonts from GraphFonts /

　　　　'GraphDataFont' =（"<sans-serif>，<MTsans-serif>"，7pt）

　　　　'GraphValueFont' =（"<sans-serif>，<MTsans-serif>"，7pt）；

　　　style GraphData1 from GraphData1 /

　　　　contrastcolor = GraphColors（'gcdata2'）

　　　　color = GraphColors（'gdata2'）；

　　　style GraphData2 from GraphData2 /

　　　　contrastcolor = GraphColors（'gcdata1'）

　　　　color = GraphColors（'gdata1'）；

　　end；

run；

/*設置森林圖樣式，DPI，圖像參數和標題等*/

ods graphics / reset width=6in height=& heightin imagename="ForestPlotColor" IMAGEFMT=EMF；

ods listing sge=off style=Styles.ForestColor93 image_dpi=& dpi；

title h=12pt "二分類數據的貝葉斯 Meta 分析的森林圖"；

title2 h=8pt '比值比 OR 及其 95% 可信區間'；

proc sgplot data=forest noautolegend nocycleattrs；

　refline studyref / lineattrs=（thickness=& thickness）transparency=0.85；

　scatter y=study2 x=oddsratio / markerattrs=（symbol=diamondfilled size=10）；

　highlow y=studyname low=lcl2 high=ucl2 / type=line；

　highlow y=studyname low=q1 high=q3 / type=bar；

　scatter y=studyname x=or / markerchar=oddsratio x2axis；

　scatter y=studyname x=lcl / markerchar=lowercl x2axis；

　scatter y=studyname x=ucl / markerchar=uppercl x2axis；

　scatter y=studyname x=wt / markerchar=weight x2axis；

　scatter y=study2 x=or / markerchar=oddsratio x2axis；

　scatter y=study2 x=lcl / markerchar=lowercl x2axis；

　scatter y=study2 x=ucl / markerchar=uppercl x2axis；

　scatter y=study2 x=wt / markerchar=weight x2axis；

　refline 1 100 / axis=x；

　refline 0.01 0.1 10 / axis=x lineattrs=（pattern=shortdash）transparency=0.5；

　inset '開腹手術組' / position=bottom；

　inset '腹腔鏡下手術組'/ position=bottomleft；

　xaxis type=log offsetmin=0 offsetmax=0.35 min=0.01 max=100 minor display=（nolabel）；

　x2axis offsetmin=0.7 display=（noticks nolabel）；

　yaxis display=（noticks nolabel）offsetmin=& pct offsetmax=& pct2 reverse；

run；

Title；

%mend；

%macro forestdata（PostSummaries_data，type）；

%if & type=1 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 2；

　run；

　%let tau=0；

%end；

%if & type=2 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 3；

　run；

　data temp_tau；

　　set & PostSummaries_data.；

　　if _N_= 2；

　　keep mean；

　run；

　proc sql noprint；

　　select mean into：tau0 from temp_tau；

　quit；

　%put & tau0.；

　%let tau=& tau0.；

%end；

data tem1；

　　set temp0；

　　keep mean stddev P2_5 P97_5；

　　rename mean=OddsRatio；rename StdDev=SE；

　　rename P2_5=Lowercl；rename P97_5=Uppercl；

run；

data Forest_plot_data；

　set Bayesian_Dichotomous_data；

/*計算效應值 OR 及標準誤*/

　OddsRatio=exp（di）；

　Lowercl=round（exp（di-1.96*sqrt（sigmai2）），.01）；

　Uppercl=round（exp（di+1.96*sqrt（sigmai2）），.01）；

　SE=（Uppercl-Lowercl）/（2*1.96）；

　studyname=author；

　weight0=1/（sigmai2+& tau.）；

　keep studyname OddsRatio SE Lowercl Uppercl weight0；

run；

proc sql noprint；

　select sum（weight0）into：sum_weight from Forest_plot_data；

quit；

　%put & sum_weight.；

data forest；

　set Forest_plot_data tem1；

run；

%forestplot（forest）；

%mend；

%forestdata（mcmc_Summaries_fixed，1）；

/*附錄代碼 5：隨機效應模型貝葉斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_random；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=246810

monitor=（theta tau2 OR）stats（percentage=（2.5 97.5））=all dic diagnostics=all；

parms theta=0 tau2 1；

prior theta ～ normal（mean = 0，var = 1e6）；

prior tau2 ～ igamma（0.001，scale =0.001）；

random deltai ～ normal（mean = theta，var = tau2）subject=study monitor=（deltai）；

model di ～ normal（deltai，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_random；run；

%forestdata（mcmc_Summaries_random，2）；