利用SAS軟件中的PROC GLIMMIX過程步實(shí)現(xiàn)二分類數(shù)據(jù)的網(wǎng)絡(luò)meta分析*

陳掌珠 黃碧芬 鄭建清△

【提 要】 目的 介紹利用SAS軟件中的PROC GLIMMIX過程步實(shí)現(xiàn)二分類數(shù)據(jù)的網(wǎng)絡(luò)Meta分析。方法 以《高級(jí)Meta分析方法——基于Stata實(shí)現(xiàn)》一書第十七章第三節(jié)的戒煙率數(shù)據(jù)作為實(shí)例，利用SAS軟件PROC GLIMMIX過程步實(shí)現(xiàn)基于廣義線性混合效應(yīng)模型的Meta分析，并提供編程代碼。結(jié)果 PROC GLIMMIX過程步可以分別擬合固定效應(yīng)模型NMA和隨機(jī)效應(yīng)模型NMA。實(shí)例數(shù)據(jù)的固定效應(yīng)NMA分析顯示，干預(yù)措施2、3、4與干預(yù)措施1相比的OR點(diǎn)估計(jì)值和95%CI分別為1.2531(0.9804，1.6026)、2.1231(1.8939，2.381)、2.2936(1.6287，3.2258)；而隨機(jī)效應(yīng)分別為1.5898(0.7918，3.1949)、2.2026(1.4641，3.3113)、2.8011(1.3106，5.988)。結(jié)論 PROC GLIMMIX可以非常方便地?cái)M合二分類數(shù)據(jù)的網(wǎng)絡(luò)Meta分析，并通過均數(shù)差分圖顯示成對(duì)比較，具有直觀的統(tǒng)計(jì)顯示效果。

在臨床醫(yī)學(xué)領(lǐng)域，傳統(tǒng)的Meta分析方法涉及的是干預(yù)和對(duì)照之間的比較，通常對(duì)照為安慰劑或其他常規(guī)療法，而干預(yù)則特指某種新療法，這種Meta分析被稱之為成對(duì)Meta分析(pairwise meta-analysis)或經(jīng)典Meta分析[1]。在存在足夠數(shù)量的針對(duì)特定臨床問題的原始研究的情況下，針對(duì)雙臂研究的傳統(tǒng)Meta分析可以獲得兩種不同治療方法的效果比較結(jié)果，提高檢驗(yàn)效能，比較容易獲得具有統(tǒng)計(jì)學(xué)顯著性的結(jié)果[2]。但在臨床實(shí)踐中，某些疾病存在多種新型治療方法，臨床醫(yī)生、患者或政府部門的決策者需要橫向比較這些治療方法的差別。在這種情況下，尤其缺乏頭對(duì)頭比較的原始研究時(shí)，傳統(tǒng)Meta分析方法變得無能為力[2]。為了克服傳統(tǒng)Meta分析的缺陷，在過去的十年里，針對(duì)隨機(jī)對(duì)照試驗(yàn)的網(wǎng)絡(luò)Meta分析(network meta-analysis，NMA)作為傳統(tǒng)meta分析的延伸而被引入循證醫(yī)學(xué)，其優(yōu)點(diǎn)是在缺乏直接頭對(duì)頭比較的情況下，實(shí)現(xiàn)多種干預(yù)措施的間接比較[3]。網(wǎng)絡(luò)Meta分析的核心思維是多重治療比較(multiple treatment comparison，MTC)，這種新方法對(duì)臨床研究人員很有吸引力，因?yàn)檫@種Meta分析方法可以確定最佳可行的干預(yù)措施。然而，盡管成對(duì)Meta分析的假設(shè)條件可以被系統(tǒng)評(píng)價(jià)員很好地理解，但是與NMA相關(guān)的假設(shè)則更復(fù)雜，而且容易被誤解。與成對(duì)Meta分析相比，網(wǎng)絡(luò)Meta分析可以獲得更多的可視化證據(jù)，從而有效評(píng)價(jià)不同干預(yù)措施之間的相對(duì)有效性，并對(duì)干預(yù)措施進(jìn)行排序。然而，在現(xiàn)實(shí)實(shí)踐中，NMA的實(shí)施仍然存在多種挑戰(zhàn)，尤其對(duì)于模型的構(gòu)建更是NMA最困難的任務(wù)。在本文中，筆者介紹了一種基于廣義線性混合模型(generalized linear mixed models，GLMMs)的NMA模型，該模型已被Jones和Madden等學(xué)者證明是一種非常高效的NMA方法[4-5]，但目前國內(nèi)尚無學(xué)者介紹該模型的具體實(shí)踐。在SAS軟件中，擬合廣義線性混合模型的過程步是PROC GLIMMIX，因此本文的編程基礎(chǔ)是PROC GLIMMIX。

基于廣義線性混合模型的網(wǎng)絡(luò)Meta分析模型簡介

1.廣義線性混合模型理論及介紹

GLMMs是廣義線性模型(generalized linear models，GLMs)和線性混合模型(linear mixed model，LMM)的擴(kuò)展形式，于二十世紀(jì)九十年代被提出，是一種比較新型的統(tǒng)計(jì)模型。GLMMs借鑒了混合模型的思想，因此可以通過隨機(jī)效應(yīng)來解釋數(shù)據(jù)間的相關(guān)性，也可以解決數(shù)據(jù)存在的過度離散(overdispersion)、異質(zhì)性(heterogeneity)等問題。此外，GLMMs模型中，因變量不再要求滿足正態(tài)分布。

GLMMs是在線性混合效應(yīng)模型的基礎(chǔ)上進(jìn)行擴(kuò)展的，LMM的基本表達(dá)形式是[6]：

y=βX+γZ+ε

(1)

對(duì)于GLMMs，令線性預(yù)測因子(linear predictor)η表示固定效應(yīng)和隨機(jī)效應(yīng)的組合(也就是說隨機(jī)誤差不包含在內(nèi))，即：

η=βX+γZ

(2)

令g(·)表示鏈接函數(shù)(link function)，用于描述影響因子線性組合的值與因變量之間的關(guān)系，h(·)為g(·)的反函數(shù)，則有：

g(E(y))=η，E(y)=h(η)=μ，因此，y=h(η)+ε。

因此，GLMMs模型表述為：

μy=E(yi|ui，Xi)=ηi=X′iβ+Z′iγi

(3)

ui為隨機(jī)效應(yīng)項(xiàng)，是反應(yīng)變量yi的條件分布的期望；條件均數(shù)ui(考慮了隨機(jī)效應(yīng))通過連接函數(shù)g(·)與條件線性預(yù)測值ηi連接。

g(μi)=ηi=X′iβ+Z′iγi

(4)

式(4)為廣義線性混合模型的一般式，yi是n×1維的反應(yīng)變量，在隨機(jī)效應(yīng)ui的條件下獨(dú)立，服從指數(shù)分布族，可以是二項(xiàng)分布、Poisson分布、Gamma分布等[7-8]。

在公式(3)中，Xi為解釋變量，β為固定效應(yīng)參數(shù)向量，ui是隨機(jī)效應(yīng)，服從均數(shù)為0，方差協(xié)方差矩陣為γ的多變量正態(tài)分布，ui解釋了由不可測因子引起的類間異質(zhì)性和同一類內(nèi)觀測到的相關(guān)，不同類間的ui是相互獨(dú)立的，Zi是與隨機(jī)效應(yīng)相關(guān)的解釋變量[9]。設(shè)計(jì)矩陣分固定效應(yīng)X與隨機(jī)效應(yīng)Z兩部分。分析的數(shù)據(jù)不同可以選擇不同的聯(lián)結(jié)函數(shù)g(·)，可以擬合含隨機(jī)效應(yīng)的logisitc回歸等多種模型。

2.基于廣義線性混合模型的網(wǎng)絡(luò)Meta分析(GLMM-NMA)理論及介紹

與成對(duì)Meta分析一樣，網(wǎng)絡(luò)Meta分析同樣可以分為固定效應(yīng)模型分析和隨機(jī)效應(yīng)分析，但在實(shí)踐操作中，現(xiàn)有的統(tǒng)計(jì)軟件，如R、GeMTC和Stata等軟件則幾乎不提供固定效應(yīng)模型的計(jì)算方法[10]。本文則介紹了基于GLMMs的固定效應(yīng)模型和隨機(jī)效應(yīng)模型。需要注意的是，本文討論的Meta分析模型是以對(duì)數(shù)比值比(log-odds ratio)為效應(yīng)尺度的，其他服從正態(tài)分布的效應(yīng)尺度也適合本文提供的模型。

(1)基于固定效應(yīng)模型的GLMM-NMA

在固定效應(yīng)NMA模型中，效應(yīng)值log-odds被假設(shè)具有線性效應(yīng)，固定效應(yīng)的組成部分包括研究效應(yīng)和治療效應(yīng)。該模型是以對(duì)數(shù)-線性鏈接函數(shù)為尺度的二分類數(shù)據(jù)的廣義線性模型的實(shí)際應(yīng)用。

假設(shè)第i(i=1，…，N)個(gè)研究中受試者在第j(j=1，…，K)個(gè)治療組中的事件發(fā)生概率為pij，對(duì)pij進(jìn)行l(wèi)ogit轉(zhuǎn)換，則有：

(5)

其中，i=1，…，N，j=1，…，K。

假設(shè)受試者是獨(dú)立的，第i研究的治療組j的樣本量為nij，發(fā)生的事件數(shù)為rij。假設(shè)rij服從二項(xiàng)式分布(binomial)，即rij～Bin(nij，pij)。需要注意的是，上述模型并沒有截距項(xiàng)。此外，這個(gè)模型的構(gòu)建并不是唯一的，研究人員可以向si添加另外組分，然后從tj中去除該組分，得到的模型完全相同。該模型的構(gòu)建有一個(gè)約束條件：默認(rèn)情況下，SAS將添加一個(gè)截距用于描述全局均值，而且約定si和tj最后一個(gè)亞變量作為對(duì)照，即sN=0和tK=0，但這里的模型省略了截距項(xiàng)，并將第一個(gè)啞變量設(shè)置為對(duì)照(即t1=0)[5]。

在SAS中，PROC GENMOD、PROC LOGISTIC和PROC GLIMMIX等三種程序步可以非常方便地?cái)M合該模型[4]。

(2)基于隨機(jī)效應(yīng)模型的GLMM-NMA

隨機(jī)效應(yīng)模型的結(jié)構(gòu)與固定效應(yīng)模型相似，多出來的部分是假設(shè)不同研究中的治療效應(yīng)因研究而異。在固定效應(yīng)模型的基礎(chǔ)上，線性預(yù)測變量增加了隨機(jī)效應(yīng)項(xiàng)：

(6)

其中N個(gè)研究νij值構(gòu)成的向量νi，具有獨(dú)立性質(zhì)的K維多變量正態(tài)分布N(0，Ω)；同樣要求t1=0。

該模型中，si仍然屬于固定效應(yīng)部分。將si鎖定為固定效應(yīng)部分，可以有效避免不同研究間的潛在偏倚問題。在這種限制條件下，可以確保不同試驗(yàn)的治療差異均來源于隨機(jī)化信息。基于PROC GLIMMIX的編程方式有很多種，附錄代碼4是GLIMMIX過程步擬合隨機(jī)效應(yīng)模型的默認(rèn)代碼。此模型假設(shè)矩陣Ω=σ2IK。也就是說，不同試驗(yàn)中治療隨機(jī)效應(yīng)是獨(dú)立的，且具有相等的方差σ2。一般不建議使用此模型，因?yàn)樵撃Ｐ蜁?huì)嚴(yán)重低估方差σ2。

Lu和Ades證明了該模型的缺陷，因此，Lu和Ades建議隨機(jī)效應(yīng)模型應(yīng)該直接關(guān)注不同試驗(yàn)的治療差異及其在不同研究間的變異[11]。Lu和Ades將pij定義為：

(7)

對(duì)于所有第i=1，…，N個(gè)研究，δi1的參數(shù)被限定為零，同時(shí)假設(shè)(δi1，…，δiK)～N((d2，…，dK)，∑)。

資料與方法

1.實(shí)例數(shù)據(jù)介紹

為了比較上述三種模型[公式(5)、公式(6)、公式(7)]的差別，本文結(jié)合實(shí)例數(shù)據(jù)進(jìn)行模型擬合。該實(shí)例數(shù)據(jù)來源于“AHCPR Smoking Cessation Guideline Panel”報(bào)告，是目前介紹NMA最廣泛引用的實(shí)例數(shù)據(jù)[11-13]，在《高級(jí)Meta分析方法——基于Stata實(shí)現(xiàn)》一書中也被引用[14]。由于Stata是目前最廣泛使用的NMA軟件，本文以文獻(xiàn)[14]中的Meta分析結(jié)果作為對(duì)照，比較本文模型的有效性。

本文采用的實(shí)例數(shù)據(jù)集的詳細(xì)信息如表1所示，該數(shù)據(jù)集文獻(xiàn)[14]中的橫向數(shù)據(jù)集經(jīng)過本文附錄代碼1轉(zhuǎn)換后的縱向數(shù)據(jù)，受篇幅所限，表1只納入前10個(gè)研究數(shù)據(jù)。

表1 待分析數(shù)據(jù)集

表1的數(shù)據(jù)存儲(chǔ)在名為NMA_data2的數(shù)據(jù)集中。相關(guān)的變量解釋如下：Design代表研究設(shè)計(jì)方法；Study代表研究編號(hào)；Narm代表某試驗(yàn)的臂數(shù)，也就是幾個(gè)干預(yù)組；Index用于指引某研究內(nèi)的干預(yù)組編號(hào)，根據(jù)治療變量t排序后分別標(biāo)示，最大值為Narm；Revindex是Index的反向標(biāo)示；Notfirst用于指引Index標(biāo)示下的非第一組臂，如果為第一組，則Notfirst取值為0，否則為1；NotLast意義與Notfirst相同，但指引最后一組。r代表某干預(yù)組內(nèi)成功戒煙患者數(shù)；n代表該干預(yù)組的樣本量；trt代表干預(yù)組(本例trt取值為1、2、3、4，代表4種干預(yù)措施，分別是no contact、self-help、individual counselling、group counselling，即無接觸戒煙、自助戒煙、個(gè)體輔導(dǎo)、團(tuán)體輔導(dǎo)等)。

2.SAS代碼簡介

本文提供了SAS軟件中實(shí)現(xiàn)了適合NMA的不同格式數(shù)據(jù)集的轉(zhuǎn)換代碼，鄢金柱等人提供了另一種轉(zhuǎn)換代碼，但本文代碼更簡潔[15]。附錄代碼同時(shí)提供了基于PROC GLIMMIX過程步擬合廣義線性混合效應(yīng)模型的代碼，為了便于讀者學(xué)習(xí)，我們同時(shí)給出了固定效應(yīng)模型和隨機(jī)效應(yīng)模型的擬合代碼。同時(shí)在PROC GLIMMIX中添加圖片輸出參數(shù)，讀者只需設(shè)置“圖片保存目錄”，即可便捷地輸出各種期刊格式的圖片到指定的文件夾。

結(jié) 果

1.基于標(biāo)準(zhǔn)策略的NMA分析結(jié)果

表2是文獻(xiàn)[14]第362頁的Meta分析結(jié)果，該結(jié)果來源于Salanti等人提出的基于標(biāo)準(zhǔn)對(duì)比模型(standard contrast-based model)，也是Stata軟件實(shí)現(xiàn)NMA的核心模型，也被稱為標(biāo)準(zhǔn)策略模型。此處提供相關(guān)Meta分析結(jié)果是為了方便比較本文構(gòu)建模型的有效性。

表2 基于Stata標(biāo)準(zhǔn)對(duì)比模型的NMA分析結(jié)果

2.基于GLMMs固定效應(yīng)模型的NMA分析結(jié)果

表3是PROC GLIMMIX擬合計(jì)算的單臂效應(yīng)值log-odds，通常不具有實(shí)際意義，但通過兩兩對(duì)比，獲得表4的效應(yīng)差值，則是系統(tǒng)評(píng)價(jià)員最關(guān)心的效應(yīng)值。如干預(yù)1和干預(yù)2對(duì)比，干預(yù)1的log-odds效應(yīng)值-干預(yù)2的log-odds效應(yīng)值即可得到效應(yīng)差值，該效應(yīng)差值為log-OR，通過逆函數(shù)exp即可獲得OR值。

表3 基于GLMMs固定效應(yīng)模型的單臂效應(yīng)值Meta分析結(jié)果

從表4可以看出，兩兩對(duì)比中，只有1-3，1-4，3-2，4-2之間的對(duì)比具有統(tǒng)計(jì)學(xué)差異，說明這些干預(yù)效應(yīng)的差值具有統(tǒng)計(jì)學(xué)意義。圖1是模型提供的差分圖，以對(duì)角線為統(tǒng)計(jì)學(xué)檢驗(yàn)的界值，效應(yīng)線跨越對(duì)角線則提示差異無統(tǒng)計(jì)學(xué)意義。

表4 基于GLMMs固定效應(yīng)模型的NMA分析的兩兩對(duì)比結(jié)果

圖1 基于GLMMs固定效應(yīng)模型的NMA分析差分圖

3.基于GLMMs隨機(jī)效應(yīng)模型(1)的NMA分析結(jié)果

表5、表6是隨機(jī)效應(yīng)模型的擬合結(jié)果。模型的參數(shù)解釋與固定效應(yīng)模型相同，但根據(jù)表6可以看出，3-2，4-2之間的對(duì)比不具有統(tǒng)計(jì)學(xué)差異，與固定效應(yīng)模型具有差別。

表5 基于GLMMs隨機(jī)效應(yīng)模型(1)的單臂效應(yīng)值Meta分析結(jié)果

表6 基于GLMMs隨機(jī)效應(yīng)模型(1)的NMA分析的兩兩對(duì)比結(jié)果

圖2 基于GLMMs隨機(jī)效應(yīng)模型(1)的NMA分析差分圖

4.基于GLMMs隨機(jī)效應(yīng)模型(2)的NMA分析結(jié)果

表7、表8是第二種隨機(jī)效應(yīng)模型的擬合結(jié)果。

表7 基于GLMMs隨機(jī)效應(yīng)模型(2)的單臂效應(yīng)值Meta分析結(jié)果

表8 基于GLMMs隨機(jī)效應(yīng)模型(2)的NMA分析的兩兩對(duì)比結(jié)果

圖3 基于GLMMs隨機(jī)效應(yīng)模型(2)的NMA分析差分圖

5.兩種GLMMs隨機(jī)效應(yīng)模型的NMA分析結(jié)果比較

表9是本文兩種隨機(jī)效應(yīng)模型的協(xié)方差參數(shù)估計(jì)值，模型(1)的協(xié)方差參數(shù)明顯低于模型(2)。

表9 基于GLMMs隨機(jī)效應(yīng)模型協(xié)方差參數(shù)估計(jì)

討 論

在臨床實(shí)踐中，許多新的治療措施缺乏“頭對(duì)頭研究(head to head)”的證據(jù)，即直接比較證據(jù)。在這種情況下，傳統(tǒng)兩兩比較的Meta分析方法無法給出比較結(jié)果，臨床醫(yī)生無法獲得孰優(yōu)孰劣的準(zhǔn)確判斷[16]。網(wǎng)絡(luò)Meta分析是一種將多項(xiàng)研究(通常是隨機(jī)試驗(yàn))的證據(jù)結(jié)合在一起，進(jìn)行多重比較，從而確定最佳干預(yù)措施的Meta分析方法。它的主要特點(diǎn)是能夠同時(shí)通過直接比較和間接比較獲得兩種治療措施的比較證據(jù)。在缺乏頭對(duì)頭研究的情況下，NMA是獲得間接證據(jù)的最優(yōu)方法。

在所有統(tǒng)計(jì)軟件中，Stata是最常用的NMA軟件之一，其核心模塊是mvmeta程序包，可以執(zhí)行多變量Meta分析和Meta回歸[17]。利用mvmeta程序包可以實(shí)現(xiàn)網(wǎng)絡(luò)Meta分析模型的頻率估計(jì)[18]，但mvmeta包需要復(fù)雜的編程代碼。相比之下，基于mvmeta包的Network程序包則更為友好，可以通過簡潔的代碼實(shí)現(xiàn)NMA[13]。鑒于Network程序包的易用性，本文采用Stata的Network程序包的計(jì)算結(jié)果作為對(duì)照。從核心統(tǒng)計(jì)模型來看，無論是Stata軟件mvmeta程序包還是Network程序包，都是運(yùn)用多元線性回歸統(tǒng)計(jì)理論來實(shí)現(xiàn)多治療方法比較的NMA[19]。這種簡單的操作模式，可以解決NMA數(shù)據(jù)的內(nèi)在復(fù)雜性，已被證明是一種非常有效的多元分析方法。

與多元線性回歸統(tǒng)計(jì)理論不同，本文構(gòu)建的NMA方法是基于廣義線性混合模型來實(shí)現(xiàn)的。這種模型最大的優(yōu)勢是無需對(duì)因變量進(jìn)行正態(tài)分布假設(shè)，同時(shí)可以通過不同的鏈接函數(shù)將觀測的均數(shù)向量與模型參數(shù)聯(lián)系起來，極大地?cái)U(kuò)大了模型的應(yīng)用范圍。GLMMs的另一個(gè)優(yōu)點(diǎn)是引入了隨機(jī)效應(yīng)項(xiàng)，并且這種隨機(jī)效應(yīng)項(xiàng)是以非線性形式的方式進(jìn)入模型。在GLMMs模型中，引入的各項(xiàng)隨機(jī)效應(yīng)變量既可以相互獨(dú)立，又可以相關(guān)。對(duì)于GLMMS模型而言，相關(guān)的結(jié)局變量可以是連續(xù)性變量，也可以是分類變量(包括二分類變量、有序多分類、無序多分類等)。因此，GLMMs不僅可以處理獨(dú)立數(shù)據(jù)，也可以處理相關(guān)數(shù)據(jù)，甚至對(duì)縱向數(shù)據(jù)、過離散數(shù)據(jù)也能夠進(jìn)行統(tǒng)計(jì)建模[20]。因此，廣義線性混合模型被譽(yù)為是最強(qiáng)大的統(tǒng)計(jì)模型之一，并被越來越多的學(xué)者引入Meta分析領(lǐng)域[9，21]，并被引入NMA[22]。我們前期的研究表明，GLMMs同樣適用于多結(jié)局測量的多變量Meta分析[23]，在經(jīng)典二分類數(shù)據(jù)的Meta分析中也是一種有效的統(tǒng)計(jì)模型[24-25]。

在SAS軟件中，實(shí)現(xiàn)廣義線性混合模型對(duì)應(yīng)的過程步是PROC GLIMMIX，但本文的統(tǒng)計(jì)模型同樣可以采用PROC GENMOD實(shí)現(xiàn)，但僅限于固定效應(yīng)模型[22]。對(duì)于多干預(yù)試驗(yàn)的Meta分析，GLMMs采用簡單的數(shù)據(jù)建模，即可實(shí)現(xiàn)這類復(fù)雜數(shù)據(jù)的NMA。從本文的分析結(jié)果來看，GLMMs的結(jié)果與Stata的結(jié)果相當(dāng)接近[14]。與Stata的組間結(jié)果相比，PROC GLIMMIX則通過回歸系數(shù)的方式，給出了單臂的對(duì)數(shù)比值(log odds)作為組水平的效應(yīng)值，系統(tǒng)評(píng)價(jià)員可以通過不同干預(yù)的單臂效應(yīng)值的減法計(jì)算來評(píng)價(jià)不同干預(yù)的效應(yīng)差別，這對(duì)系統(tǒng)評(píng)價(jià)員來說極為方便。

與Stata軟件中的Network模塊給出的森林圖不同，PROC GLIMMIX給出的是均數(shù)差分圖，并通過線型和線色的差異來顯示具有統(tǒng)計(jì)學(xué)意義的成對(duì)Meta分析結(jié)果。以本文圖1為例，差分圖輸出的結(jié)果表明，干預(yù)1與干預(yù)4的LS均值差異在5%的水平上具有統(tǒng)計(jì)學(xué)意義，其LS均值差為(-0.8298±0.1747)，相應(yīng)的OR為0.436，95%CI為(0.31，0.614)，提示與干預(yù)4相比，干預(yù)1更能提高戒煙效果。

SAS軟件需要一定的編程基礎(chǔ)，這是該軟件的不足，但以PROC GLIMMIX實(shí)現(xiàn)基于廣義線性混合效應(yīng)模型的Meta分析，則為NMA的實(shí)施提供了新的思路與途徑[26]。PROC GLIMMIX可以通過“random”語句來設(shè)置隨機(jī)效應(yīng)模型，通過簡單的編程，可以同時(shí)給出固定效應(yīng)模型Meta分析結(jié)果和隨機(jī)效應(yīng)模型Meta分析結(jié)果。據(jù)筆者所知，目前常用的主流NMA軟件，如R、GeMTC軟件、OpenBUGS等，只有PROC GLIMMIX過程步是基于GLMMs實(shí)現(xiàn)二分類數(shù)據(jù)的NMA。PROC GLIMMIX無法輸出森林圖和網(wǎng)狀關(guān)系圖，這是本文模型的不足之處，但森林圖的繪制，可以借助SAS的其他制圖程序來實(shí)現(xiàn)。

綜上所述，PROC GLIMMIX可以非常方便地?cái)M合二分類數(shù)據(jù)的網(wǎng)絡(luò)Meta分析，并通過均數(shù)差分圖顯示成對(duì)比較，具有非常直觀的統(tǒng)計(jì)顯示效果。

附錄代碼1：

title “數(shù)據(jù)集建立及資料錄入”；

Data NMA_data1；

length design $20.；

input study design ar an br bn cr cn dr dn@@；

datalines；

1 ACD 9 140 ..23 140 10 138

……

；

run；

附錄代碼2：

title “橫向數(shù)據(jù)轉(zhuǎn)換成縱向數(shù)據(jù)”；

data NMA_data2；

set NMA_data1 end=myend；

drop ar an br bn cr cn dr dn；

Narm=((ar+an) > 0) +((br+bn) > 0) +((cr+cn) > 0)+((dr+dn) > 0)；

Index= 0；

if ar+an >0 then do；

Revindex=Narm-index；index=index+1；

Notfirst=(index>1)；NotLast=(index

r=ar；n=an；trt=1；output；

end；

if br+bn >0 then do；

Revindex=Narm-index；index=index+1；

Notfirst=(index>1)；NotLast=(index

r=br；n=bn；trt=2；output；

end；

if cr+cn >0 then do；

Revindex=Narm-index；index=index+1；

Notfirst=(index>1)；NotLast=(index

r=cr；n=cn；trt=3；output；

end；

if dr+dn >0 then do；

Revindex=Narm-index；index=index+1；

Notfirst=(index>1)；NotLast=(index

r=dr；n=dn；trt=4；output；

end；

run；

附錄代碼3：

title “模型1：基于固定效應(yīng)模型的NMA分析”；

proc glimmix data=NMA_data2 order=data plots=Diffogram；

ods listing gpath=“圖片保存目錄” style=default image_dpi=300；

ods graphics/reset

imagename=“基于固定效應(yīng)模型的NMA分析”

outputfmt=epsi；

title “基于固定效應(yīng)模型的NMA分析”；

class study Trt；

model R/N = Study Trt /link=logit dist=bin ddfm=none；

lsmeans Trt /diff cl oddsratios；

run；

附錄代碼4：

title “模型2：基于隨機(jī)效應(yīng)模型的NMA分析(1)”；

proc glimmix data=NMA_data2 order=data plots=Diffogram；

ods listing gpath=“圖片保存目錄” style=default image_dpi=300；

ods graphics/reset

imagename=“基于隨機(jī)效應(yīng)模型的NMA分析(1)”

outputfmt=epsi；

title “基于隨機(jī)效應(yīng)模型的NMA分析(1)”；

class study Trt；

model R/N = Study Trt /link=logit dist=bin ddfm=none；

random Trt / subject=Study；

lsmeans Trt /diff cl oddsratios；

run；

附錄代碼5：

title “模型3：基于隨機(jī)效應(yīng)模型的NMA分析(2)”；

data jr；

set NMA_data2；

array x[4]x1-x4；

do i=1 to 4；

if i<= narm then x[i]=sqrt(0.5)*((i=index)-1/narm)；

else x[i]=0；

end；

run；

proc glimmix data=jr method=QUAD order=data plots=Diffogram；

ods listing gpath=“圖片保存目錄” style=default image_dpi=300；

ods graphics/reset

imagename=“基于隨機(jī)效應(yīng)模型的NMA分析(2)”

outputfmt=epsi；

title “基于隨機(jī)效應(yīng)模型的NMA分析(2)”；

class study Trt；

model R/N = Trt Study/link=logit dist=bin ddfm=none；

random X1 X2 X3 X4/ subject=study type=TOEP(1)；

lsmeans Trt /diff cl oddsratios；

run；