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

· 循證理論與實踐 · 論著 ·

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

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

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

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

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

1　metamisc程序包

1.1metamisc程序包簡介及其安裝加載 metamisc

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

2　CopulaREMADA程序包

2.1CopulaREMADA程序包簡介及安裝加載CopulaREMADA程序包可以通過概率模型來執(zhí)行診斷準(zhǔn)確性試驗的meta分析，最新版本為V-0.9，更新時間為2015-06-11。該程序包實現(xiàn)診斷試驗Meta分析的基本思路是，首先建立診斷準(zhǔn)確性試驗樣本的概率模型，以此模型為介體，通過高斯求積［10］的方法來計算該模型的MLE，即通過樣本估計總體參數(shù)的最可能值，則可估算出該診斷試驗的診斷價值。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é)點和權(quán)重，并通過擬牛頓法（quasi-Newton）產(chǎn)生穩(wěn)定的概率矩陣，比較發(fā)現(xiàn)=15足以產(chǎn)生至少有四個小數(shù)位的高精度數(shù)值。具體命令如下：

nq=15

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

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

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

2.3.2將樣本數(shù)據(jù)導(dǎo)入模型并計算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”可分別得到一個結(jié)果列表，此處僅展示命令“est.n”結(jié)果（表3）及講解，minimum為負(fù)對數(shù)似然的最小值；estimate即最大似然估計值（maximum likelihood estimation，MEI）；gradient表示梯度；iterations表示迭代運(yùn)算次數(shù)。

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

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é)

兩個程序包都是基于經(jīng)典頻率學(xué)方法進(jìn)行DTA meta分析，不同的是，metamisc程序包建立了雙變量隨機(jī)效應(yīng)模型，在迭代運(yùn)算過程中采用的是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é)果更為精確，總體而言兩者結(jié)果相差不大。兩個程序包都提供了繪圖功能，metamisc程序包所繪制的圖形比較簡單而且直觀，但它只給出了合并敏感度的結(jié)果，其圖形所呈現(xiàn)的意義并不大；CopulaREMADA程序包所繪制的SROC曲線圖顯示了不同邊緣分布的特點，敏感度與特異度的變化關(guān)系體現(xiàn)了雙變量模型的優(yōu)越性，更利于曲線下面積的計算。

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

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

［3］ 張俊，徐志偉，李克. 診斷性試驗Meta分析的效應(yīng)指標(biāo)評價［J］. 中國循證醫(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］ 祝慧萍，方龍，夏欣，等. 雙變量模型在診斷試驗Meta分析中的應(yīng)用［J］. 華中科技大學(xué)學(xué)報（醫(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

湖北省教育廳重點項目（D20142102）
共同第一作者：楊廉潔

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

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

10.3969/j.issn.1674-4055.2016.04.03

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

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

注：TP：真陽性數(shù)；FN：假陰性數(shù)；FP：假陽性數(shù)；TN：真陰性數(shù)；Mark為檢測樣本，其中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：假陽性率；beta1：敏感度的logit值；beta2：假陽性率的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ù)據(jù)進(jìn)行l(wèi)ogit轉(zhuǎn)換以滿足該模型的正態(tài)分布需要。命令如下：

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

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

1.3.2亞組分析 以Mark為分組依據(jù)，data［which（d ata$Mark==1），］即表示data這個數(shù)據(jù)集里Mark標(biāo)記為1的數(shù)據(jù)，Mark為2時類推，同樣通過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í)行繪圖功能，圖形展示了靈敏度與假陽性率的相關(guān)變化。

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

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

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

圖1整體分析曲線圖

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è)置曲線的線條類型，pch設(shè)置合并結(jié)果的形狀，其可選參數(shù)與R軟件plot（ ）一致，此處對分組的合并結(jié)果形狀及線條類型分別設(shè)置以便區(qū)分（圖2），實線和實心點是Mark為1所繪制的曲線圖，虛線和空心點是Mark為2所繪制的曲線圖，該程序包不足之處是無法設(shè)置線條及點的顏色。

圖2亞組分析曲線圖