应用R软件metamisc程序包及CopulaREMADA程序包实现诊断准确性试验的Meta分析

王权，杨廉洁，何倩，喻亚宇，许杨鹏，张超

· 循证理论与实践 · 论著 ·

应用R软件metamisc程序包及CopulaREMADA程序包实现诊断准确性试验的Meta分析

王权1，杨廉洁2，何倩3，喻亚宇3，许杨鹏3，张超4

诊断性准确性试验（diagnostic test accuracy，DTA）Meta分析是一种全面评价诊断试验证据准确性的研究方法，R软件metamisc程序包与CopulaREMADA程序包是基于经典频率学方法用于DTA Meta分析制作及图形绘制的程序包，与传统的双变量模型相比，其所建立的分析模型减少了组间差异，简化了其繁琐的迭代运算过程，使诊断试验评价指标更加准确与高效。

诊断准确性试验；Meta分析；metamisc程序包；CopulaREMADA程序包；R软件

诊断准确性试验（diagnostic test accuracy，DTA）的Meta分析［1］是通过准确度评价指标，获得检测结果与参考标准结果的符合程度。传统的双变量模型可以保留灵敏度和特异度这两个不同指标各自的性质并获得在负相关下的综合估计量，但忽略了研究组内灵敏度和特异度之间的相关性［2，3］。R软件metamisc程序包及CopulaREMADA程序包所建立的双变量模型，能够通过减少其组间差异和抽样误差，提高诊断试验评价指标的准确性［4］。本文将以Walusimbi等［5］发表的文章中的GeneXpert 组的数据为例，来演示R软件metamisc程序包及CopulaREMADA程序包实例操作。

1　metamisc程序包

1.1metamisc程序包简介及其安装加载 metamisc

程序包由Thomas Debray研发，目前，最新更新时间为：2013-05-30，其最新版本为：V-0.1.1。metamisc程序包可以进行诊断试验和干预实验的Meta分析，该程序包进行诊断试验Meta分析是基于经典频率学的双变量随机效应模型，运用Newton-Raphson法［6］计算最大似然估计值（maximum likelihood estimation，MLE）［7］。在当前版本中，贝叶斯方法仅用于单变量模型的计算，其双变量模型的应用会在后续版本中更新。

2　CopulaREMADA程序包

2.1CopulaREMADA程序包简介及安装加载CopulaREMADA程序包可以通过概率模型来执行诊断准确性试验的meta分析，最新版本为V-0.9，更新时间为2015-06-11。该程序包实现诊断试验Meta分析的基本思路是，首先建立诊断准确性试验样本的概率模型，以此模型为介体，通过高斯求积［10］的方法来计算该模型的MLE，即通过样本估计总体参数的最可能值，则可估算出该诊断试验的诊断价值。CopulaREMADA程序包的安装及加载命令如下：

install.packages（‘CopulaREMADA'）

library（CopulaREMADA）

2.2数据加载 CopulaREMADA程序包与metamisc程序包的数据排列格式一致，数据加载的命令也相同。

2.3数据分析

2.3.1产生概率矩阵 使用高斯求积方法选取适当的节点和权重，并通过拟牛顿法（quasi-Newton）产生稳定的概率矩阵，比较发现=15足以产生至少有四个小数位的高精度数值。具体命令如下：

nq=15

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

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

其中，gauss.quad.prob（ ）为产生节点和权重的命令，“uniform”表示产生均匀分布的概率，该命令提供的分布类型还有“normal”，“beta”和“gamma”。命令“mgrid＜-meshgrid（gl$n，gl$n）” 即产生一个行数和列数为gl$n概率矩阵。

2.3.2将样本数据导入模型并计算MEI

attach（data）

est.n=countermonotonicCopulaREMADA. norm（TP，FN， FP，TN，gl，mgrid）

est.n

est.b=countermonotonicCopulaREMADA.beta（TP，FN，FP，TN，gl，mgrid）

est.b

est.n表示边缘为normal分布，est.b表示边缘为beta分布，命令“est.n”和命令“est.b”可分别得到一个结果列表，此处仅展示命令“est.n”结果（表3）及讲解，minimum为负对数似然的最小值；estimate即最大似然估计值（maximum likelihood estimation，MEI）；gradient表示梯度；iterations表示迭代运算次数。

2.3.3绘制受试者工作特征曲线（summary receiver operating characteristic，SROC） 受试者工作特征曲线通过分量回归技术拟建模型，比较双变量随机效应不同分布类型的特征，具体命令如下：

SROC（est.b$e，est.n$e，TP，FN，FP，TN）

其中est.b$e和est.n$e分别表示normal边缘分布和beta边缘分布的MEI（图3），黑色表示normal边缘分布SROC，红色表示beta边缘分布。

3　小结

两个程序包都是基于经典频率学方法进行DTA meta分析，不同的是，metamisc程序包建立了双变量随机效应模型，在迭代运算过程中采用的是Nelder-Mead的方法，结果为0.675（0.589，0.751），由于考虑了研究间相关性并给出了置信区间，其结果更为准确；而CopulaREMADA程序包建立的是双变量介体混合模型［11］，运用了高斯求积的方法，结果为0.671，由于进行多次反复的迭代运算，其结果更为精确，总体而言两者结果相差不大。两个程序包都提供了绘图功能，metamisc程序包所绘制的图形比较简单而且直观，但它只给出了合并敏感度的结果，其图形所呈现的意义并不大；CopulaREMADA程序包所绘制的SROC曲线图显示了不同边缘分布的特点，敏感度与特异度的变化关系体现了双变量模型的优越性，更利于曲线下面积的计算。

两个程序包都能将已发表的诊断试验Meta分析数据与新的诊断试验数据进行合并，实现临床诊断证据的更新，或者当诊断试验文献中没有给出原始数据时，程序包可以将其效应指标直接合并，从而进行评价。

［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分析及系统评价方法［J］. 转化医学杂志，2015，4（1）：51-5.

［3］ 张俊，徐志伟，李克. 诊断性试验Meta分析的效应指标评价［J］. 中国循证医学杂志，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分析中的应用［J］. 华中科技大学学报（医学版），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 十堰，十堰市太和医院（湖北医药学院附属医院）口腔科；2442000 十堰，十堰市太和医院（湖北医药学院附属医院）院务办公室；3442000 十堰，湖北医药学院口腔医学院12级；4442000十堰，十堰市太和医院（湖北医药学院附属医院）循证医学中心

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

10.3969/j.issn.1674-4055.2016.04.03

值得注意的是，metamisc程序包的计算内核是通过JAGS（Just Another Gibbs Sampler）软件来实现的，因此在使用时还需事先安装JAGS软件，JAGS软件当前最新版本为4.0.0，下载地址为：http：//sourceforge.net/projects/mcmc-jags/files/。R软件metamisc程序包的安装和加载命令如下：install.packages（‘metamisc'）library（metadisc）

1.2数据加载 首先，在数据加载之前，需要对数据进行排列格式，具体数据排列格式详见表1。同样，在上述数据排列完成后，储存在桌面的Rwork文件中的data.txt文本中。随后，开始进行数据 加载，具体命令如下：data＜-read.table（"C：\Users\Administrator\ Desktop\Rwork\data.txt"，header=TRUE，sep=""，na. strings="NA"， dec="."， strip.white=TRUE）

1.3数据分析 在完成数据加载之后，开始实现数据的分析，该程序包提供多种数据模型分析功能，在进行诊断准确性试验 Meta分析时只需要使用riley（ ）函数进行数据合并，summary（ ）函数进行结果汇总，以及plot（ ）函数进行图形绘制即可。riley（ ）函数构建了由Riley提出的一种双变量随机效应模型［8］，由于各种研究的诊断界点选取不同，加上诊断试验操作规范的差异等原因，使得诊断试验各个研究之间的差异较大，研究间的异质性较强，不能采用常规的固定效应模型，而应采用随机效应模型进行分析。双变量模型分析就是通过对各研究的灵敏度和特异度进行logit变换，使之在服从正态分布的前提下进行随机效应模型分析，并获得灵敏度和特异度之间负向关联性大小值的一种方法［9］。该模型考虑整体相关性参数而不是将研究内相关性与研究间相关性分开，这个模型不是完全分层的，并且能够估算出抽样误差，抽样误差并不等于一般模型的研究间变异。因此当研究间变异很大时，可用的原始研究太少时，或者一般模型所估算出的研究间相关性为1或-1时，这个模型就非常适用。本文将对整体分析和亚组分析分别讲解。

表1数据排列表

注：TP：真阳性数；FN：假阴性数；FP：假阳性数；TN：真阴性数；Mark为检测样本，其中1为Sputum，2为Various

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

表2summary（fit）命令汇总结果

注：Sens：敏感度；FPR：假阳性率；beta1：敏感度的logit值；beta2：假阳性率的logit值；psi1：beta1的抽样误差；psi2：beta2的抽样误差；rho：psi1与psi2间的相关性

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（ ）函数进行数据合并时，将会对数据进行logit转换以满足该模型的正态分布需要。命令如下：

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

test.accuracy表示确定该数据为诊断数据，Nelder-Mead表示迭代运算方法，计算及汇总结果见表2。

1.3.2亚组分析 以Mark为分组依据，data［which（d ata$Mark==1），］即表示data这个数据集里Mark标记为1的数据，Mark为2时类推，同样通过summary（）函数得到汇总结果，此处不再赘述。

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图形绘制 该程序包采用函数plot（ ）执行绘图功能，图形展示了灵敏度与假阳性率的相关变化。

1.4.1整体分析图形绘制 具体命令如下：

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

其中，plotsumm为逻辑函数，表示是否显示合并结果，默认值为TRUE；plotnumerics表示是否显示灵敏度和假阳性率的汇总表，默认值为TRUE；cex.numerics设置汇总表的字体大小。在所绘制图形（图1）中，椭圆是置信区间曲线，中间的小方形为合并结果。

图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 设置曲线的线条类型，pch设置合并结果的形状，其可选参数与R软件plot（ ）一致，此处对分组的合并结果形状及线条类型分别设置以便区分（图2），实线和实心点是Mark为1所绘制的曲线图，虚线和空心点是Mark为2所绘制的曲线图，该程序包不足之处是无法设置线条及点的颜色。

图2亚组分析曲线图