利用SAS軟件PROC MCMC過程步實(shí)現(xiàn)診斷性試驗(yàn)的貝葉斯Meta分析

許佩華 鄭建清

診斷性試驗(yàn)最常報(bào)告的結(jié)局指標(biāo)是敏感度和特異度，其中敏感度反映的是正確診斷疾病或狀態(tài)的概率，而特異度反映正確識(shí)別其誤診的概率。診斷性試驗(yàn)的Meta分析可以定量地總結(jié)相似的原始研究，獲得各個(gè)研究的匯總效應(yīng)值[1]。研究結(jié)局的Meta分析可以使用固定效應(yīng)模型或隨機(jī)效應(yīng)模型進(jìn)行，前者假設(shè)研究間差異僅源于抽樣誤差，后者假設(shè)研究間差異不僅源于抽樣誤差，而且還來自研究之間未知的真實(shí)差異。對(duì)于隨機(jī)效應(yīng)薈萃分析，最常用的統(tǒng)計(jì)方法是雙變量隨機(jī)效應(yīng)模型(BREM)和層次綜合受試者工作特征(HSROC)模型[2-5]。

Meta分析可以基于頻率論(經(jīng)典)方法或貝葉斯方法實(shí)現(xiàn)。頻率論認(rèn)為概率是重復(fù)事件(假設(shè))發(fā)生的頻率；貝葉斯方法認(rèn)為概率是數(shù)值的可信程度[6]。換句話說，頻率論認(rèn)為模型參數(shù)是固定的，數(shù)據(jù)是隨機(jī)的；貝葉斯認(rèn)為模型參數(shù)是隨機(jī)的，數(shù)據(jù)是固定的。貝葉斯統(tǒng)計(jì)學(xué)與經(jīng)典統(tǒng)計(jì)學(xué)在統(tǒng)計(jì)推斷的理念和方法上存在本質(zhì)上的區(qū)別，包括對(duì)信息的利用，對(duì)未知參數(shù)θ的解釋以及統(tǒng)計(jì)推斷理念上的差異等[7,8]。Rutter和Gatsonis最早提出了用于診斷性試驗(yàn)Meta分析的貝葉斯HSROC方法，并且可以使用免費(fèi)統(tǒng)計(jì)軟件WinBUGS實(shí)現(xiàn)[5, 9]。然而，利用SAS軟件實(shí)現(xiàn)貝葉斯分析的文獻(xiàn)較少被報(bào)道。

Menke編制了一套基于SAS PROC MCMC過程步的Meta分析代碼，可以實(shí)現(xiàn)診斷性試驗(yàn)的貝葉斯Meta分析[10]，本文對(duì)其進(jìn)行了介紹，并結(jié)合實(shí)例介紹SAS代碼的實(shí)際應(yīng)用。有關(guān)貝葉斯統(tǒng)計(jì)方法的理論基礎(chǔ)，張?zhí)灬缘萚8]在《高級(jí)Meta分析-基于Stata實(shí)現(xiàn)》一書中已有詳細(xì)介紹。

1 Meta分析模型簡(jiǎn)介

1.1 診斷性試驗(yàn)BREM模型簡(jiǎn)介 BREM是廣義線性混合模型的特例[3]。在診斷性試驗(yàn)的Meta分析中，原始研究提供了2×2表格數(shù)據(jù)，從而可以計(jì)算出相應(yīng)的敏感度和特異度。對(duì)于研究i，TPi、FPi、FNi、TNi分別代表試驗(yàn)i內(nèi)的真陽(yáng)性、假陽(yáng)性、假陰性和真陰性。RPi=TNi+FPi表示研究i中金標(biāo)準(zhǔn)試驗(yàn)檢測(cè)出的陽(yáng)性結(jié)果(RP)。RNi=TNi+FPi表示研究i中金標(biāo)準(zhǔn)試驗(yàn)檢測(cè)出的陰性結(jié)果。

則TPi和RPi，TNi和RNi構(gòu)成二項(xiàng)式分布[10]：

TPi～Binomial(RPi,πA,i)；TNi～Binomial(RNi,πB,i)(公式1)

公式中，πA,i是研究i待估計(jì)的真實(shí)敏感度(概率)。πB,i是研究i待估計(jì)的真實(shí)特異度(概率)。 靈敏度和特異度的BREM由下面的公式給出[11]：

其中μA,i=logit(πA,i)，表示研究i待估計(jì)的真實(shí)敏感度的logit轉(zhuǎn)換值。μB,i=logit(πB,i)表示研究i待估計(jì)的真實(shí)特異度的logit轉(zhuǎn)換值。logit(π)= log(π /(1-π))表示概率π的logit轉(zhuǎn)換。

1.2 診斷性試驗(yàn)HSROC模型簡(jiǎn)介 HSROC模型由Rutter等提出，是對(duì)綜合受試者工作特征曲線模型(ROC)的擴(kuò)展，該模型利用靈敏度和特異度這一常用配對(duì)指標(biāo)對(duì)多個(gè)診斷性試驗(yàn)進(jìn)行綜合評(píng)價(jià)[5]。該模型具體如下：

假設(shè)研究i中某一患者疾病狀態(tài)j的陽(yáng)性概率為πij，其中j=0表示無病，j=1表示有病，則靈敏度PAi=πi1，特異度PBi=1-πi0。基于靈敏度和特異度，HSROC模型為：

水平1(研究?jī)?nèi))：logit(πij)=θi+αiXijexp(-βXij)

2 資料與方法

2.1 示例數(shù)據(jù)介紹 以發(fā)表的關(guān)于高分辨率超聲診斷顳下頜關(guān)節(jié)盤前移位的準(zhǔn)確性研究為例，文章直接報(bào)道了基于HSROC 統(tǒng)計(jì)模型的結(jié)果[12]，采用的統(tǒng)計(jì)軟件是STATA 12.0。張?zhí)灬缘冉榻B了有關(guān)HSROC在診斷試驗(yàn)Meta分析中的應(yīng)用及Stata實(shí)現(xiàn)[13]。作為對(duì)比，本文提供了SAS軟件的計(jì)算結(jié)果。

采用文后附錄代碼1錄入示例數(shù)據(jù)(表1)，并將數(shù)據(jù)集存儲(chǔ)在“file_data”。數(shù)據(jù)集“file_bidata”將數(shù)據(jù)集將file_data分成每個(gè)研究的2個(gè)記錄，并添加變量“A_pos”和“B_neg”，分別識(shí)別敏感度和特異度，構(gòu)建雙變量數(shù)據(jù)(表2)。

需要注意的是，本文模型擬合假設(shè)診斷正確例數(shù)服從二項(xiàng)式分布，故無需對(duì)數(shù)據(jù)集中零事件數(shù)的研究進(jìn)行校正。

2.2 SAS代碼介紹 本文提供的附錄代碼包括數(shù)據(jù)集的錄入及貝葉斯雙變量隨機(jī)效應(yīng)Meta分析的實(shí)現(xiàn)。本文利用SAS軟件PROC MCMC過程步編制了BREM的程序代碼，并提供了HSROC模型相關(guān)的參數(shù)計(jì)算。BREM的公式2可以用SAS 9.4中的PROC MCMC實(shí)現(xiàn)(見附錄代碼)。SAS PROC MCMC的詳細(xì)信息在聯(lián)機(jī)手冊(cè)中提供?；陔S機(jī)效應(yīng)模型的SAS代碼(附錄代碼3)主要介紹如下。

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

表2 待分析的雙變量數(shù)據(jù)集

注 A_pos代表真陽(yáng)性例數(shù)，B_neg代表真陰性例數(shù)

附錄代碼第1～4行調(diào)用PROC MCMC，其中包含有關(guān)采樣和調(diào)整的代碼選項(xiàng)(第2行)，也包含報(bào)告輸出的選項(xiàng)(第3～4行)。丟棄的開動(dòng)采樣數(shù)量(nbi=100 000)被設(shè)定為MCMC采樣(nmc=1000 000)的10%。MCMC采樣間隔采用10來抽樣以減少自相關(guān)。因此，后驗(yàn)分布由100 000個(gè)樣本組成，對(duì)應(yīng)于每百分位數(shù)1 000個(gè)樣本。第6～8行3個(gè)ARRAY語句聲明了模型變量，在“[]”括號(hào)提供了數(shù)組下標(biāo)。“random eta”語句行聲明指定了納入分析的單個(gè)研究的雙變量eta數(shù)組。雙變量數(shù)組mu和協(xié)方差矩陣Sigma分別代表了公式2的BREM的5個(gè)參數(shù)。在第8行中，cov_uAB2實(shí)際上與cov_uAB是相同的。第10～11行PARMS語句將匯總估計(jì)“mu”和隨機(jī)效應(yīng)協(xié)方差矩陣“Sigma”聲明為模型參數(shù)，第12～13行prior語句指定先驗(yàn)分布。在第12行中，無信息先驗(yàn)(noninformative prior)的均勻先驗(yàn)分布(uniform priors)被用于mu_A和mu_B，范圍從-10到10。在靈敏度/特異度尺度上，對(duì)應(yīng)于0.005%至99.995%的足夠?qū)挼膮?shù)范圍。具有較大方差(例如100)的多變量正態(tài)先驗(yàn)也可以用于代替均勻先驗(yàn)。第13行將Wishart分布指定為隨機(jī)效應(yīng)協(xié)方差矩陣的多變量先驗(yàn)。

第9行的ARRAY語句定義了Wishart分布的同一性標(biāo)度矩陣。第14～16行針對(duì)file_ bidata數(shù)據(jù)集的每個(gè)研究的特定數(shù)據(jù)記錄進(jìn)行評(píng)估。第17行實(shí)現(xiàn)公式1，第15～16行提供了靈敏度和特異度的對(duì)數(shù)轉(zhuǎn)換，第14行實(shí)現(xiàn)了公式2的雙變量模型。在第18～25行中，提供了進(jìn)一步的模型計(jì)算結(jié)果，包括診斷比值比及其對(duì)數(shù)值、敏感度、特異度、陽(yáng)性似然比、陰性似然比及其對(duì)數(shù)轉(zhuǎn)換值等。第26～30行計(jì)算了與HSROC模型相關(guān)的5個(gè)參數(shù)。第32行調(diào)用內(nèi)置繪圖函數(shù)繪制Caterpillar圖。

固定效應(yīng)模型的SAS代碼省略了與隨機(jī)效應(yīng)相關(guān)的所有代碼(見附錄代碼2)。

運(yùn)行程序時(shí)，每個(gè)受監(jiān)控的模型參數(shù)都會(huì)生成一個(gè)MCMC樣本鏈(Chain)，這通過第4行代碼指定。這樣的馬爾可夫鏈(the Markov chains)是一個(gè)隨機(jī)過程，其未來狀態(tài)僅取決于當(dāng)前狀態(tài)，且獨(dú)立于過去狀態(tài)。作為穩(wěn)健結(jié)果的先決條件，馬爾可夫鏈應(yīng)收斂到每個(gè)模型參數(shù)的獨(dú)特后驗(yàn)分布。這可以通過程序的診斷代碼(即第3行)進(jìn)行評(píng)估。

2.3 評(píng)估鏈?zhǔn)諗?通過SAS PROC MCMC實(shí)現(xiàn)的馬爾可夫鏈的收斂，可以采用蹤跡圖和自相關(guān)圖來可視化評(píng)估，特別是通過評(píng)估鏈的混合性和平穩(wěn)性。除了圖形給出收斂性評(píng)估外，PROC MCMC同時(shí)給出不同收斂性診斷的統(tǒng)計(jì)值：蒙特卡羅標(biāo)準(zhǔn)誤、自相關(guān)參數(shù)、Geweke診斷、Raftery-Lewis診斷、Heidelberger-Welch診斷、有效樣本量等。需要注意的是，SAS PROC MCMC的所有收斂診斷都是使用默認(rèn)參數(shù)執(zhí)行的。顯著性水平設(shè)定為P<0.05[10]。

3 結(jié)果

3.1 貝葉斯固定效應(yīng)模型Meta分析結(jié)果 運(yùn)行附錄代碼2后，可以獲得貝葉斯固定效應(yīng)模型Meta分析結(jié)果。包括主要統(tǒng)計(jì)結(jié)果部分和鏈?zhǔn)諗啃栽u(píng)估部分。主要統(tǒng)計(jì)結(jié)果部分包含4個(gè)表格，分別為后驗(yàn)摘要表、后驗(yàn)可信區(qū)間表、后驗(yàn)相關(guān)系數(shù)矩陣和后驗(yàn)協(xié)方差矩陣。系統(tǒng)評(píng)價(jià)員主要關(guān)心后驗(yàn)摘要表和后驗(yàn)可信區(qū)間表。鏈?zhǔn)諗啃栽u(píng)估部分包含7個(gè)表格和針對(duì)本(Highest posterior density，HPD)文待研究變量的自相關(guān)圖。由于篇幅關(guān)系，這里僅提供后驗(yàn)摘要表和后驗(yàn)可信區(qū)間表的結(jié)果。后驗(yàn)可信區(qū)間表比較的是95%等尾(Equal-Tail)可信區(qū)間和95%最大后驗(yàn)密度可信區(qū)間的差異。資料見表3、表4。讀者可以運(yùn)行本文的附錄代碼，獲取其他表格數(shù)據(jù)及自相關(guān)圖。

表3 固定效應(yīng)模型后驗(yàn)摘要表

注 mu_A：敏感度的對(duì)數(shù)值，mu_B：特異度的對(duì)數(shù)值，LOR：診斷比值比的對(duì)數(shù)值，Sens：敏感度，spec：特異度，DOR：診斷比值比，LRP：陽(yáng)性似然比，LRN：陰性似然比，log_LRP：LRP的對(duì)數(shù)值，log_LRN：LRN的對(duì)數(shù)值，ALPHA：HSROC模型α參數(shù)，THETA：HSROC模型θ參數(shù)

表4 固定效應(yīng)模型后驗(yàn)可信區(qū)間表

注 mu_A：敏感度的對(duì)數(shù)值，mu_B：特異度的對(duì)數(shù)值，LOR：診斷比值比的對(duì)數(shù)值，Sens：敏感度，spec：特異度，DOR：診斷比值比，LRP：陽(yáng)性似然比，LRN：陰性似然比，log_LRP：LRP的對(duì)數(shù)值，log_LRN：LRN的對(duì)數(shù)值，ALPHA：HSROC模型α參數(shù)，THETA：HSROC模型θ參數(shù)

表3、4數(shù)據(jù)顯示，高分辨率超聲診斷可復(fù)性盤前移位(ADDWR)的合并敏感度(Sen)為0.823(95%CI：0.796～0.851)，合并特異度(Spe)為0.834(95%CI：0.801～0.865)，診斷比值比(DOR)為23.94(95%CI：17.568～32.119)，陽(yáng)性似然比(LRP)為5.007(95%CI：4.120～6.098)，陰性似然比(LRN)為0.2113(95%CI：0.178～0.247)；HSROC模型相關(guān)參數(shù)：α=3.164(95%CI：2.866～3.469)；θ=-0.036(95%CI：-0.187～0.114)。

3.2 貝葉斯隨機(jī)效應(yīng)模型Meta分析結(jié)果 運(yùn)行附錄代碼3后，可以獲得貝葉斯隨機(jī)效應(yīng)模型Meta分析結(jié)果。也包括主要統(tǒng)計(jì)結(jié)果部分和鏈?zhǔn)諗啃栽u(píng)估部分(表5和6)。需要注意的是，表5和6沒有提供單個(gè)試驗(yàn)的隨機(jī)效應(yīng)估計(jì)值。讀者可以運(yùn)行本文的附錄代碼，獲取表5和6的單個(gè)試驗(yàn)的隨機(jī)效應(yīng)估計(jì)值、其他表格數(shù)據(jù)及自相關(guān)圖。

表5和6數(shù)據(jù)顯示，高分辨率超聲診斷ADDWR的合并Sen 0.852(95%CI：0.766～0.924)，合并Spe 0.865(95%CI：0.760～0.944)，DOR 53.413(95%CI：11.855～170.4)，LRP 7.306(95%CI：3.325～15.946)，LRN 0.174(95%CI：0.083～0.296)；HSROC模型相關(guān)參數(shù)：α=3.739(95%CI：2.526～5.162)；θ=0.153(95%CI：-0.404～0.732)；BETA=0.233(95%CI：-0.384～0.866)；var_uAlpha=4.105(95%CI：1.171～11.813)；var_uTheta=0.134(95%CI：0.043～0.356)。

3.3 貝葉斯隨機(jī)效應(yīng)模型Meta分析的鏈?zhǔn)諗啃栽u(píng)估結(jié)果 在本文的示例數(shù)據(jù)中，運(yùn)行附錄代碼后，可以獲得MCMC采樣評(píng)估圖，包括核密度圖、蹤跡圖和自相關(guān)圖。篇幅原因只提供了mu_A的結(jié)果(圖1)。圖1說明了針對(duì)SAS PROC MCMC參數(shù)mu_A的MCMC采樣。在蹤跡圖中(圖1A)，馬爾可夫鏈很好地混合并通過遍歷密度較低區(qū)域來探索其分布。在自相關(guān)圖中(圖1B)，在LAG=5的地方，mu_A的驗(yàn)后自相關(guān)是0.0296，因此很小。核密度圖顯示(圖1C)，mu_A的后驗(yàn)分布類似于正態(tài)分布，其參數(shù)Mean=1.7850，標(biāo)準(zhǔn)差SD=0.3319。

3.4 貝葉斯隨機(jī)效應(yīng)模型Meta分析的Caterpillar森林圖 SAS PROC MCMC可以調(diào)用內(nèi)置宏命令“%CATER”來繪制森林圖。本文示例數(shù)據(jù)繪制的森林圖見圖2。

表5 隨機(jī)效應(yīng)模型后驗(yàn)摘要表

注 mu_A：敏感度的對(duì)數(shù)值，mu_B：特異度的對(duì)數(shù)值，LOR：診斷比值比的對(duì)數(shù)值，Sens：敏感度，spec：特異度，DOR：診斷比值比，LRP：陽(yáng)性似然比，LRN：陰性似然比，log_LRP：LRP的對(duì)數(shù)值，log_LRN：LRN的對(duì)數(shù)值，ALPHA：HSROC模型α參數(shù)，THETA：HSROC模型θ參數(shù)

表6 隨機(jī)效應(yīng)模型后驗(yàn)可信區(qū)間表

4 討論

SAS 軟件使用廣泛，其強(qiáng)大的編程思維可以靈活地處理各種數(shù)據(jù)。相對(duì)于經(jīng)典的統(tǒng)計(jì)推論方法，貝葉斯方法長(zhǎng)短處兼有[14]。隨著計(jì)算方法的進(jìn)步，貝葉斯推論方法的發(fā)展加快。在SAS9.2版本之前，SAS無法實(shí)現(xiàn)貝葉斯模型分析。貝葉斯分析模型在SAS軟件里對(duì)口的過程是PROC MCMC，作為試驗(yàn)版本迄始于SAS 9.2，在SAS 9.3中有所增強(qiáng)，并提供了對(duì)應(yīng)于隨機(jī)效應(yīng)模型的Random Statement等。本文介紹了通過SAS 9.4版本的PROC MCMC 過程、隨機(jī)效應(yīng)模型來實(shí)現(xiàn)診斷性試驗(yàn)的貝葉斯Meta 分析。

基于特異度和靈敏度的BREM目前是最受推薦的診斷準(zhǔn)確性薈萃分析模型[6, 15, 16]。通過BREM的隨機(jī)效應(yīng)，研究人員可以獲得未知的靈敏度和特異度的真實(shí)差異。這種真實(shí)差異已被證明在原始研究之間是存在的。此外，該模型認(rèn)為這些敏感度和特異度的隨機(jī)效應(yīng)可能是相關(guān)的。Harbord等已經(jīng)證明BREM在數(shù)學(xué)上等同于HSROC模型[4]。本文SAS代碼獲取HSROC參數(shù)，便是直接從雙變量的結(jié)果獲得。

在SAS軟件中，可以通過不同的過程步實(shí)現(xiàn)靈敏度和特異度的BREM，包括SASPROC GLIMMIX和SAS PROC MIXED等[3,17]。在SAS PROC GLIMMIX中，公式2的廣義線性混合模型由最大似然擬合，其中邊際分布在迭代中以數(shù)值方式近似，直到滿足收斂準(zhǔn)則。 SAS PROC GLIMMIX假定其模型參數(shù)的預(yù)定義固定分布，并且計(jì)算速度非常快。van Houwelingen等介紹了SAS PROC MIXED實(shí)現(xiàn)雙變量隨機(jī)效應(yīng)模型的原理，并提供了相關(guān)的SAS代碼[17]。

圖1 MCMC采樣分析圖(變量=mu_A)

圖2 各個(gè)研究的匯總Caterpillar森林圖

與SAS PROC GLIMMIX和SAS PROC MIXED相比，貝葉斯SAS PROC MCMC通過數(shù)值模擬獲得其結(jié)果，而不是通過分析或迭代近似獲得其結(jié)果。但需要更多的計(jì)算時(shí)間。由此產(chǎn)生的后驗(yàn)分布取決于先驗(yàn)和數(shù)據(jù)，不需要正態(tài)分布。貝葉斯分析可以直接估計(jì)參數(shù)的任何函數(shù)。例如，合并的陽(yáng)性似然比和陰性似然比及其95%可信區(qū)間均可以從合并的靈敏度和特異度估計(jì)，同理，可以獲得陽(yáng)性驗(yàn)后概率和陰性驗(yàn)后概率。

SAS PROC MCMC實(shí)現(xiàn)了馬爾可夫鏈蒙特卡羅(MCMC)方法，該方法結(jié)合了蒙特卡羅積分和馬爾可夫鏈采樣。蒙特卡羅積分是一種數(shù)值積分方法，通過采用離散樣本簡(jiǎn)化連續(xù)分布。樣本量較大時(shí)，樣本的直方圖近似于連續(xù)分布，可以通過樣本進(jìn)行有關(guān)分布的任何推斷，例如估計(jì)平均值和95%可信區(qū)間，提供了一種從目標(biāo)密度中提取樣本的方法。馬爾可夫鏈?zhǔn)且粋€(gè)隨機(jī)過程，其未來狀態(tài)僅依賴于當(dāng)前狀態(tài)并且與過去無關(guān)。這種屬性允許馬爾可夫鏈簡(jiǎn)化復(fù)雜問題，因?yàn)殒溨械南乱粋€(gè)樣本僅取決于前一個(gè)樣本[18]。馬爾可夫鏈需要一個(gè)在目標(biāo)分布范圍內(nèi)的初始值。樣本通過SAS PROC MCMC中的隨機(jī)漫步Metropolis算法和WinBUGS中的Gibbs采樣獲得[19]。

應(yīng)用MCMC方法，需要評(píng)估Markov鏈的收斂診斷，而這種收斂診斷可以幫助解決兩個(gè)問題。首先，明確馬爾可夫鏈?zhǔn)欠褡兊闷椒€(wěn)(最初的非平穩(wěn)樣本被丟棄為“退火(Burn-in)”)。其次，明確后驗(yàn)分布需要多少樣本。理論上，樣本越大，越容易準(zhǔn)確地估計(jì)模型參數(shù)，但計(jì)算時(shí)間會(huì)成比例增加。正常情況下，計(jì)算時(shí)間是限制MCMC樣本的主要原因。

本文提供的代碼程序是在Menke提供的SAS 軟件進(jìn)行貝葉斯診斷性試驗(yàn)Meta 分析的PROC MCMC 過程代碼上修改的。在Menke的文章中，作者針對(duì)敏感度和特異度的貝葉斯雙變量隨機(jī)效應(yīng)薈萃分析，提出了統(tǒng)計(jì)模型及其在SAS PROC MCMC的實(shí)施過程。提供的程序代碼在50個(gè)Meta分析中進(jìn)行了實(shí)證評(píng)估，結(jié)果與相應(yīng)的WinBUGS方法非常相似[9; 10]。本文在Menke提供的代碼上進(jìn)行實(shí)例應(yīng)用，并對(duì)代碼做了簡(jiǎn)單修改，可以直接獲得診斷比值比、陽(yáng)性似然比及陰性似然比，讀者可以比較兩種代碼的異同。鑒于學(xué)術(shù)性，讀者在學(xué)習(xí)本文提供的代碼時(shí)，請(qǐng)務(wù)必參考學(xué)習(xí)Menke的原文，并引用原文作為參考文獻(xiàn)[10]。

在本文引用的示例數(shù)據(jù)中，董小宇等利用STATA實(shí)現(xiàn)了HSROC模型的Meta分析[12]，這給本文的SAS代碼學(xué)習(xí)提供了便利的實(shí)證基礎(chǔ)。通過對(duì)比發(fā)現(xiàn)，董小宇等的研究結(jié)果中，高分辨率超聲診斷ADDWR的Sen=0.84(95%CI：0.76～0.90)，Spe=0.87(95%CI：0.77～0.93)，而本文Sen為0.852(95%CI：0.766～0.924)，Spe為0.865(95%CI：0.760～0.944)； 兩者統(tǒng)計(jì)結(jié)果極其接近，證明SAS PROC MCMC是可以實(shí)現(xiàn)HSROC模型計(jì)算的。需要注意的是，在貝葉斯統(tǒng)計(jì)理論中，參數(shù)θ是一個(gè)隨機(jī)變量，結(jié)合樣本信息和先驗(yàn)信息可以構(gòu)造一個(gè)區(qū)間，這個(gè)區(qū)間稱之為可信區(qū)間(Credit interval)，這個(gè)應(yīng)該區(qū)別于經(jīng)典統(tǒng)計(jì)學(xué)派的置信區(qū)間(Confidence interval)。貝葉斯理論可以簡(jiǎn)單描述為“有95%的概率使得參數(shù)θ落在可信區(qū)間內(nèi)”。董小宇等的文章提供的“CI”屬于經(jīng)典統(tǒng)計(jì)學(xué)的“置信區(qū)間”，在統(tǒng)計(jì)學(xué)概念上，略有不同。

5 結(jié)論

鑒于SAS軟件強(qiáng)大的編程能力，貝葉斯SAS PROC MCMC可以非常方便地實(shí)現(xiàn)基于敏感度和特異度的BREM的診斷性試驗(yàn)的Meta分析。

/* 附錄代碼1 建立數(shù)據(jù)集*/

ods html file='Data.html' style=statistical;

DATAfile_data;

length author $20;

input study author $ year tru_pos fal_pos fal_neg tru_neg;

ref_pos = tru_pos + fal_neg;

ref_neg = tru_neg + fal_pos;

datalines;

1 Yangjieping 2012 24 1 6 9

……

;run;

PROCPRINTdata=file_data;run;

DATAfile_bidata; set file_data;

status='A_pos'; true=tru_pos; total=ref_pos; output;

status='B_neg'; true=tru_neg; total=ref_neg; output;

keep study status true total;run;

PROCPRINTdata=file_bidata;run;

/* 附錄代碼2 固定效應(yīng)模型*/

ods html file='Fixed-effects_model.html' style=statistical;

ods graphics on;

PROCMCMCdata=file_bidata outpost=file_postout_fixed propcov=quanew

nbi=100000nmc=1000000ntu=1000thin=10seed=1simreport=100

diagnostics=all stats(percentage=(2.5510255075909597.5))=all

monitor=(mu_A mu_B LOR sens spec DOR LRP LRN log_LRP log_LRN ALPHA THETA) dic;

title 'Bayesian bivariate fixed-effects model';

/*變量聲明和初始化 */

array mu[2] mu_A mu_B; /*合并敏感度和特異度的Logit轉(zhuǎn)換*/

parms mu:0;

prior mu: ～ unif(-10,10);/*對(duì)每個(gè)數(shù)據(jù)記錄評(píng)估以下指令*/

if status='A_pos' then p = logistic(mu_A);

if status='B_neg' then p = logistic(mu_B);

model true ～ binomial(total, p);

/*計(jì)算其他雙變量參數(shù)*/

LOR = mu_A + mu_B; /*診斷比值比的對(duì)數(shù)值*/

sens = logistic(mu_A); /*合并敏感度*/

spec = logistic(mu_B); /*合并特異度*/

DOR=(sens/(1-spec))/((1-sens)/spec);

LRP = sens/(1-spec); /*陽(yáng)性似然比*/

LRN = (1-sens)/spec; /*陰性似然比*/

log_LRP = log10(sens/(1-spec)); /*陽(yáng)性似然比的log10轉(zhuǎn)換*/

log_LRN = log10((1-sens)/spec); /*陰性似然比的log10轉(zhuǎn)換*/

/*計(jì)算HSROC參數(shù)*/

THETA = (mu_A - mu_B)/2;

ALPHA = (mu_A + mu_B);

run;

/*附錄代碼3 隨機(jī)效應(yīng)模型*/

ods html file='Random-effects_model.html' style=statistical;

ods graphics on;

#1PROCMCMCdata=file_bidata outpost=file_postout_random propcov=quanew

#2nbi=100000nmc=1000000ntu=1000thin=10seed=1simreport=100

#3diagnostics=all stats(percentage=(2.5510255075909597.5))=all

#4monitor=(mu_A mu_B LOR sens spec DOR LRP LRN log_LRP log_LRN var_uA var_uB cov_uAB ALPHA THETA BETA var_uAlpha var_uTheta) dic;

#5title 'Bayesian bivariate random-effects model';

#6array eta[2] eta_A eta_B; /*納入研究特定的隨機(jī)效應(yīng)*/

#7array mu[2] mu_A mu_B;

#8array Sigma[2,2] var_uA cov_uAB cov_uAB2 var_uB;

#9array S[2,2] (1001);

#10parms mu:0;

#11parms Sigma {0.01000.01};

#12prior mu: ～ unif(-10,10);

#13prior Sigma ～ iwish(2, S);

#14random eta ～ mvn(mu, Sigma) subject=study monitor=(eta_A eta_B);

#15if status='A_pos' then p = logistic(eta_A);

#16if status='B_neg' then p = logistic(eta_B);

#17model true ～ binomial(total, p);

/*計(jì)算其他雙變量參數(shù)*/

#18LOR = mu_A + mu_B; /*診斷比值比的對(duì)數(shù)值*/

#19sens = logistic(mu_A); /*合并敏感度*/

#20spec = logistic(mu_B); /*合并特異度*/

#21DOR=(sens/(1-spec))/((1-sens)/spec);

#22LRP = sens/(1-spec); /*陽(yáng)性似然比*/

#23LRN = (1-sens)/spec; /*陰性似然比*/

#24log_LRP = log10(sens/(1-spec)); /*陽(yáng)性似然比的log10轉(zhuǎn)換*/

#25log_LRN = log10((1-sens)/spec); /*陰性似然比的log10轉(zhuǎn)換*/

/*計(jì)算HSROC參數(shù)*/

#26THETA = (mu_A*(var_uB/var_uA)**0.25- mu_B*(var_uA/var_uB)**0.25)/2;

#27ALPHA = (mu_A*(var_uB/var_uA)**0.25+ mu_B*(var_uA/var_uB)**0.25);

#28BETA = log(sqrt(var_uB/var_uA));

#29var_uTheta = (sqrt(var_uA*var_uB) - cov_uAB)/2;

#30var_uAlpha = (sqrt(var_uA*var_uB) + cov_uAB)*2;

#31run;

#32%CATER(data=file_postout_random, var=eta:);