王權(quán)，楊廉潔，何倩，喻亞宇，許楊鵬，張超

· 循證理論與實(shí)踐 · 論著 ·

應(yīng)用R軟件metamisc程序包及CopulaREMADA程序包實(shí)現(xiàn)診斷準(zhǔn)確性試驗(yàn)的Meta分析

王權(quán)1，楊廉潔2，何倩3，喻亞宇3，許楊鵬3，張超4

診斷性準(zhǔn)確性試驗(yàn)（diagnostic test accuracy，DTA）Meta分析是一種全面評(píng)價(jià)診斷試驗(yàn)證據(jù)準(zhǔn)確性的研究方法，R軟件metamisc程序包與CopulaREMADA程序包是基于經(jīng)典頻率學(xué)方法用于DTA Meta分析制作及圖形繪制的程序包，與傳統(tǒng)的雙變量模型相比，其所建立的分析模型減少了組間差異，簡(jiǎn)化了其繁瑣的迭代運(yùn)算過(guò)程，使診斷試驗(yàn)評(píng)價(jià)指標(biāo)更加準(zhǔn)確與高效。

診斷準(zhǔn)確性試驗(yàn)；Meta分析；metamisc程序包；CopulaREMADA程序包；R軟件

診斷準(zhǔn)確性試驗(yàn)（diagnostic test accuracy，DTA）的Meta分析［1］是通過(guò)準(zhǔn)確度評(píng)價(jià)指標(biāo)，獲得檢測(cè)結(jié)果與參考標(biāo)準(zhǔn)結(jié)果的符合程度。傳統(tǒng)的雙變量模型可以保留靈敏度和特異度這兩個(gè)不同指標(biāo)各自的性質(zhì)并獲得在負(fù)相關(guān)下的綜合估計(jì)量，但忽略了研究組內(nèi)靈敏度和特異度之間的相關(guān)性［2，3］。R軟件metamisc程序包及CopulaREMADA程序包所建立的雙變量模型，能夠通過(guò)減少其組間差異和抽樣誤差，提高診斷試驗(yàn)評(píng)價(jià)指標(biāo)的準(zhǔn)確性［4］。本文將以Walusimbi等［5］發(fā)表的文章中的GeneXpert 組的數(shù)據(jù)為例，來(lái)演示R軟件metamisc程序包及CopulaREMADA程序包實(shí)例操作。

1　metamisc程序包

1.1metamisc程序包簡(jiǎn)介及其安裝加載 metamisc

程序包由Thomas Debray研發(fā)，目前，最新更新時(shí)間為：2013-05-30，其最新版本為：V-0.1.1。metamisc程序包可以進(jìn)行診斷試驗(yàn)和干預(yù)實(shí)驗(yàn)的Meta分析，該程序包進(jìn)行診斷試驗(yàn)Meta分析是基于經(jīng)典頻率學(xué)的雙變量隨機(jī)效應(yīng)模型，運(yùn)用Newton-Raphson法［6］計(jì)算最大似然估計(jì)值（maximum likelihood estimation，MLE）［7］。在當(dāng)前版本中，貝葉斯方法僅用于單變量模型的計(jì)算，其雙變量模型的應(yīng)用會(huì)在后續(xù)版本中更新。

2　CopulaREMADA程序包

2.1CopulaREMADA程序包簡(jiǎn)介及安裝加載CopulaREMADA程序包可以通過(guò)概率模型來(lái)執(zhí)行診斷準(zhǔn)確性試驗(yàn)的meta分析，最新版本為V-0.9，更新時(shí)間為2015-06-11。該程序包實(shí)現(xiàn)診斷試驗(yàn)Meta分析的基本思路是，首先建立診斷準(zhǔn)確性試驗(yàn)樣本的概率模型，以此模型為介體，通過(guò)高斯求積［10］的方法來(lái)計(jì)算該模型的MLE，即通過(guò)樣本估計(jì)總體參數(shù)的最可能值，則可估算出該診斷試驗(yàn)的診斷價(jià)值。CopulaREMADA程序包的安裝及加載命令如下：

install.packages（‘CopulaREMADA'）

library（CopulaREMADA）

2.2數(shù)據(jù)加載 CopulaREMADA程序包與metamisc程序包的數(shù)據(jù)排列格式一致，數(shù)據(jù)加載的命令也相同。

2.3數(shù)據(jù)分析

2.3.1產(chǎn)生概率矩陣 使用高斯求積方法選取適當(dāng)?shù)墓?jié)點(diǎn)和權(quán)重，并通過(guò)擬牛頓法（quasi-Newton）產(chǎn)生穩(wěn)定的概率矩陣，比較發(fā)現(xiàn)=15足以產(chǎn)生至少有四個(gè)小數(shù)位的高精度數(shù)值。具體命令如下：

nq=15

gl=gauss.quad.prob（nq，"uniform"）

mgrid＜- meshgrid（gl$n，gl$n）

其中，gauss.quad.prob（ ）為產(chǎn)生節(jié)點(diǎn)和權(quán)重的命令，“uniform”表示產(chǎn)生均勻分布的概率，該命令提供的分布類(lèi)型還有“normal”，“beta”和“gamma”。命令“mgrid＜-meshgrid（gl$n，gl$n）” 即產(chǎn)生一個(gè)行數(shù)和列數(shù)為gl$n概率矩陣。

2.3.2將樣本數(shù)據(jù)導(dǎo)入模型并計(jì)算MEI

attach（data）

est.n=countermonotonicCopulaREMADA. norm（TP，F(xiàn)N， FP，TN，gl，mgrid）

est.n

est.b=countermonotonicCopulaREMADA.beta（TP，F(xiàn)N，F(xiàn)P，TN，gl，mgrid）

est.b

est.n表示邊緣為normal分布，est.b表示邊緣為beta分布，命令“est.n”和命令“est.b”可分別得到一個(gè)結(jié)果列表，此處僅展示命令“est.n”結(jié)果（表3）及講解，minimum為負(fù)對(duì)數(shù)似然的最小值；estimate即最大似然估計(jì)值（maximum likelihood estimation，MEI）；gradient表示梯度；iterations表示迭代運(yùn)算次數(shù)。

2.3.3繪制受試者工作特征曲線(xiàn)（summary receiver operating characteristic，SROC） 受試者工作特征曲線(xiàn)通過(guò)分量回歸技術(shù)擬建模型，比較雙變量隨機(jī)效應(yīng)不同分布類(lèi)型的特征，具體命令如下：

SROC（est.b$e，est.n$e，TP，F(xiàn)N，F(xiàn)P，TN）

其中est.b$e和est.n$e分別表示normal邊緣分布和beta邊緣分布的MEI（圖3），黑色表示normal邊緣分布SROC，紅色表示beta邊緣分布。

3　小結(jié)

兩個(gè)程序包都是基于經(jīng)典頻率學(xué)方法進(jìn)行DTA meta分析，不同的是，metamisc程序包建立了雙變量隨機(jī)效應(yīng)模型，在迭代運(yùn)算過(guò)程中采用的是Nelder-Mead的方法，結(jié)果為0.675（0.589，0.751），由于考慮了研究間相關(guān)性并給出了置信區(qū)間，其結(jié)果更為準(zhǔn)確；而CopulaREMADA程序包建立的是雙變量介體混合模型［11］，運(yùn)用了高斯求積的方法，結(jié)果為0.671，由于進(jìn)行多次反復(fù)的迭代運(yùn)算，其結(jié)果更為精確，總體而言?xún)烧呓Y(jié)果相差不大。兩個(gè)程序包都提供了繪圖功能，metamisc程序包所繪制的圖形比較簡(jiǎn)單而且直觀，但它只給出了合并敏感度的結(jié)果，其圖形所呈現(xiàn)的意義并不大；CopulaREMADA程序包所繪制的SROC曲線(xiàn)圖顯示了不同邊緣分布的特點(diǎn)，敏感度與特異度的變化關(guān)系體現(xiàn)了雙變量模型的優(yōu)越性，更利于曲線(xiàn)下面積的計(jì)算。

兩個(gè)程序包都能將已發(fā)表的診斷試驗(yàn)Meta分析數(shù)據(jù)與新的診斷試驗(yàn)數(shù)據(jù)進(jìn)行合并，實(shí)現(xiàn)臨床診斷證據(jù)的更新，或者當(dāng)診斷試驗(yàn)文獻(xiàn)中沒(méi)有給出原始數(shù)據(jù)時(shí)，程序包可以將其效應(yīng)指標(biāo)直接合并，從而進(jìn)行評(píng)價(jià)。

［1］ Moses LE，Shapiro D，Littenberg B. Combining independent studies of a diagnostic test into a summary ROCcurve： data-analytic approaches and some additional considerations［J］. Stat Med，1993，12（14）：1293-316.

［2］ 劉鴻，周潔，馮巧靈，等. 基于檢驗(yàn)效能的診斷性試驗(yàn)Meta分析及系統(tǒng)評(píng)價(jià)方法［J］. 轉(zhuǎn)化醫(yī)學(xué)雜志，2015，4（1）：51-5.

［3］ 張俊，徐志偉，李克. 診斷性試驗(yàn)Meta分析的效應(yīng)指標(biāo)評(píng)價(jià)［J］. 中國(guó)循證醫(yī)學(xué)雜志，2013，13（7）：890-5.

［4］ Johannes B Reitsma，Afina S Glas，Anns WS Rutjes，et al. Bivariate Analysis of Sensitivity and SpecificityProduces Informative Summary Measures in Diagnostic reviews［J］. J Clin Epidemiol，2005，58，982-90.

［5］ Walusimbi S，Bwanga F， Costa A. Meta-analysis to compare the accuracy of GeneXpert， MODS and theWHO 2007 algorithm for diagnosis smear-negative pulmonary tuberculosis［J］. BMC Infect Dis，2013，13：507.

［6］ Maiti T，Pradhan V. A comparative study of the bias corrected estimates in logistic regression［J］. Stat MethodsMed Res，2008，17（6）：621-34.

［7］ Luo X，Willse JT. A Dual-Purpose Rasch Model with Joint Maximum Likelihood Estimation［J］. Stat Methods Med Res，2015，16（3）：298-314. ［8］ Riley RD，Abrams KR，Lambert PC，et al. An Evaluation of Bivariate Random-effects Meta-analys is for theJoint Synthesis of Two Correlated Outcomes［J］. Stat Med，2007，26（1）：78-97.

［9］ 祝慧萍，方龍，夏欣，等. 雙變量模型在診斷試驗(yàn)Meta分析中的應(yīng)用［J］. 華中科技大學(xué)學(xué)報(bào)（醫(yī)學(xué)版），2010，39（1）：78-81.

［10］ García-Martínez L，Rosete-Aguilar M，Gardu?o-Mejia J. Gauss-Legendre quadrature method used toevaluate the spatio-temporal intensity of ultrashort pulses in the focal region of lenses［J］. Appl Opt，2012，51（3）：306-15.

［11］ Nikoloulopoulos AK. A vine copula mixed effect model for trivariate meta-analysis of diagnostic test accuracystudies accounting for disease prevalence［J］. Stat Methods Med Res，2015，11.

本文編輯：姚雪莉

Implementation of diagnostic test accuracy with metamisc package and CopulaREMADA package in Rsoftware: a Meta-analysis

WANG Quan*， YANG Lian-jie， HE Qian， YU Ya-yu， XU Yang-peng， ZHANG Chao.*Department of Stomatology， Taihe Hospital， Hubei University of Medicine， Shiyan， 442000， China.

ZHANG Chao， E-mail: zhangchao0803@126.com

The Meta-analysis of diagnostic test accuracy （DTA） is a research method that comprehensive evaluates the accuracy of diagnostic test evidence. The metamisc package and CopulaREMADA package in R software are used for implementing Meta-analysis and graphic plotting of DTA based on classic frequency study approach. Compared with traditional bivariate model， the analysis models established by these two packages can reduce intergroup difference and simplify the tedious iterative operation process， which make DTA evaluation indicator more accurate and efficient.

Diagnostic test accuracy； Meta-analysis； Metamisc package； CopulaREMADA package； R software

R4

A

1674-4055（2016）04-0392-04

湖北省教育廳重點(diǎn)項(xiàng)目（D20142102）
共同第一作者：楊廉潔

1442000 十堰，十堰市太和醫(yī)院（湖北醫(yī)藥學(xué)院附屬醫(yī)院）口腔科；2442000 十堰，十堰市太和醫(yī)院（湖北醫(yī)藥學(xué)院附屬醫(yī)院）院務(wù)辦公室；3442000 十堰，湖北醫(yī)藥學(xué)院口腔醫(yī)學(xué)院12級(jí)；4442000十堰，十堰市太和醫(yī)院（湖北醫(yī)藥學(xué)院附屬醫(yī)院）循證醫(yī)學(xué)中心

張超，E-mail：zhangchao0803@126.com

10.3969/j.issn.1674-4055.2016.04.03

值得注意的是，metamisc程序包的計(jì)算內(nèi)核是通過(guò)JAGS（Just Another Gibbs Sampler）軟件來(lái)實(shí)現(xiàn)的，因此在使用時(shí)還需事先安裝JAGS軟件，JAGS軟件當(dāng)前最新版本為4.0.0，下載地址為：http：//sourceforge.net/projects/mcmc-jags/files/。R軟件metamisc程序包的安裝和加載命令如下：install.packages（‘metamisc'）library（metadisc）

1.2數(shù)據(jù)加載 首先，在數(shù)據(jù)加載之前，需要對(duì)數(shù)據(jù)進(jìn)行排列格式，具體數(shù)據(jù)排列格式詳見(jiàn)表1。同樣，在上述數(shù)據(jù)排列完成后，儲(chǔ)存在桌面的Rwork文件中的data.txt文本中。隨后，開(kāi)始進(jìn)行數(shù)據(jù) 加載，具體命令如下：data＜-read.table（"C：\Users\Administrator\ Desktop\Rwork\data.txt"，header=TRUE，sep=""，na. strings="NA"， dec="."， strip.white=TRUE）

1.3數(shù)據(jù)分析 在完成數(shù)據(jù)加載之后，開(kāi)始實(shí)現(xiàn)數(shù)據(jù)的分析，該程序包提供多種數(shù)據(jù)模型分析功能，在進(jìn)行診斷準(zhǔn)確性試驗(yàn) Meta分析時(shí)只需要使用riley（ ）函數(shù)進(jìn)行數(shù)據(jù)合并，summary（ ）函數(shù)進(jìn)行結(jié)果匯總，以及plot（ ）函數(shù)進(jìn)行圖形繪制即可。riley（ ）函數(shù)構(gòu)建了由Riley提出的一種雙變量隨機(jī)效應(yīng)模型［8］，由于各種研究的診斷界點(diǎn)選取不同，加上診斷試驗(yàn)操作規(guī)范的差異等原因，使得診斷試驗(yàn)各個(gè)研究之間的差異較大，研究間的異質(zhì)性較強(qiáng)，不能采用常規(guī)的固定效應(yīng)模型，而應(yīng)采用隨機(jī)效應(yīng)模型進(jìn)行分析。雙變量模型分析就是通過(guò)對(duì)各研究的靈敏度和特異度進(jìn)行l(wèi)ogit變換，使之在服從正態(tài)分布的前提下進(jìn)行隨機(jī)效應(yīng)模型分析，并獲得靈敏度和特異度之間負(fù)向關(guān)聯(lián)性大小值的一種方法［9］。該模型考慮整體相關(guān)性參數(shù)而不是將研究?jī)?nèi)相關(guān)性與研究間相關(guān)性分開(kāi)，這個(gè)模型不是完全分層的，并且能夠估算出抽樣誤差，抽樣誤差并不等于一般模型的研究間變異。因此當(dāng)研究間變異很大時(shí)，可用的原始研究太少時(shí)，或者一般模型所估算出的研究間相關(guān)性為1或-1時(shí)，這個(gè)模型就非常適用。本文將對(duì)整體分析和亞組分析分別講解。

表1數(shù)據(jù)排列表

注：TP：真陽(yáng)性數(shù)；FN：假陰性數(shù)；FP：假陽(yáng)性數(shù)；TN：真陰性數(shù)；Mark為檢測(cè)樣本，其中1為Sputum，2為Various

StudyYear TPFNFPTNMark Helb201038150251 Malbruny2011600732 Bowles20112140231 Moure201161170201 Marlowe201143120471 Theron20112225193191 Rachow201111711021 Scott201111731041 Lawn2011233023201 Ioannidis20112932321 Miller2011322581 Teo20111362422 Nicol2011251801661 Rachow20121470221 Safianowska20124401812

表2summary（fit）命令匯總結(jié)果

注：Sens：敏感度；FPR：假陽(yáng)性率；beta1：敏感度的logit值；beta2：假陽(yáng)性率的logit值；psi1：beta1的抽樣誤差；psi2：beta2的抽樣誤差；rho：psi1與psi2間的相關(guān)性

Call：rileyES（X = NULL， Y1 = logit.sens， Y2 = logit.fpr， vars1 = var.logit.sens，vars2 = var.logit.fpr， optimization = optimization， control = control）Estimate2.50%97.50% Sens0.6748390.5885190.750721 FPR0.026350.0144930.047441 beta10.7301530.3578471.102459 beta2-3.60957-4.21946-2.99967 psi10.5376090.2184560.856763 psi20.625950.0744371.177463 rho0.255593-0.316890.691571

1.3.1整體分析 使用riley（ ）函數(shù)進(jìn)行數(shù)據(jù)合并時(shí)，將會(huì)對(duì)數(shù)據(jù)進(jìn)行l(wèi)ogit轉(zhuǎn)換以滿(mǎn)足該模型的正態(tài)分布需要。命令如下：

fit＜-riley（data，type="test. accuracy"，optimization="Nelder-Mead"）summary（fit）

test.accuracy表示確定該數(shù)據(jù)為診斷數(shù)據(jù)，Nelder-Mead表示迭代運(yùn)算方法，計(jì)算及匯總結(jié)果見(jiàn)表2。

1.3.2亞組分析 以Mark為分組依據(jù)，data［which（d ata$Mark==1），］即表示data這個(gè)數(shù)據(jù)集里Mark標(biāo)記為1的數(shù)據(jù)，Mark為2時(shí)類(lèi)推，同樣通過(guò)summary（）函數(shù)得到匯總結(jié)果，此處不再贅述。

fit1＜-riley（data［which（data$Mark==1），］，type="test. accuracy"，optimization="Nelder-Mead"） summary（fit1）

fit2＜-riley（data［which（data$Mark==2），］，type="test. accuracy"，optimization="Nelder-Mead"）summary（fit2）

1.4圖形繪制 該程序包采用函數(shù)plot（ ）執(zhí)行繪圖功能，圖形展示了靈敏度與假陽(yáng)性率的相關(guān)變化。

1.4.1整體分析圖形繪制 具體命令如下：

plot（fit， plotsumm=TRUE， plotnumerics=TRUE，cex. numerics=1）

其中，plotsumm為邏輯函數(shù)，表示是否顯示合并結(jié)果，默認(rèn)值為T(mén)RUE；plotnumerics表示是否顯示靈敏度和假陽(yáng)性率的匯總表，默認(rèn)值為T(mén)RUE；cex.numerics設(shè)置匯總表的字體大小。在所繪制圖形（圖1）中，橢圓是置信區(qū)間曲線(xiàn)，中間的小方形為合并結(jié)果。

圖1整體分析曲線(xiàn)圖

1.4.2亞組分析圖形繪制 具體命令如下：

plot（fit1，plotnumerics=FALSE， lty=1，pch=20）

plot（fit2，plotnumerics=FALSE，add=TRUE，lty=2，pch=1）

其中add=TRUE表示在原圖上繪制圖形，lty 設(shè)置曲線(xiàn)的線(xiàn)條類(lèi)型，pch設(shè)置合并結(jié)果的形狀，其可選參數(shù)與R軟件plot（ ）一致，此處對(duì)分組的合并結(jié)果形狀及線(xiàn)條類(lèi)型分別設(shè)置以便區(qū)分（圖2），實(shí)線(xiàn)和實(shí)心點(diǎn)是Mark為1所繪制的曲線(xiàn)圖，虛線(xiàn)和空心點(diǎn)是Mark為2所繪制的曲線(xiàn)圖，該程序包不足之處是無(wú)法設(shè)置線(xiàn)條及點(diǎn)的顏色。

圖2亞組分析曲線(xiàn)圖