软件程序包在网状分析中的应用

1、中国循证医学杂志 2014, 14(3): 355360R软件nlme程序包在网状Meta分析中的应用张1.超121,2*十堰 442000）；2.学院附属太和医院口腔医学中心（十堰 442000）学院附属太和医院循证医学中心（摘要 nlme 程序包是基于广义最小二乘法和线性混合效应模型研发的、可通过 R 软件实现广义线性和非线性混合效应模型下的 Meta 分析。该程序包实现 Meta 分析时，需要对数据先行转化为效应量的对数值本文介绍了使用 R 软件 nlme 程序包实现网状 Meta 分析的过程，详细呈现了如何转化数据这一步骤。网状 Meta 分析；广义线性混合效应模型；R 语言；nlme

2、 程序包进行。Implementation of Network Meta-Analysis with nlme Package in R SoftwareZHANG Chao1, NIU Yu-ming2, ZENG Xian-tao1,2*1. Center for Evidence-Based Medicine and Clinical Research, Taihe Hospital, Hubei University of Medicine, Shiyan 442000, China;2. Department of Stomatology, Taihe Hospital, Hube

3、i University of Medicine, Shiyan 442000, ChinaAbstract The nlme package is developed based on thralized least squares (gls) and linear mixed-effects mod(lme). It can perform meta-analysis based on linear and nonlinear mixed effects moin R language. When conductingmeta-analysis using nlme package in

4、R language, the first step is to translate the data into its logarithm estimation. In this article, we introduce how to perform network meta-analysis using R language nlme package and show the core step of data translation in detail.Key words Network meta-analysis; Generalized linear mixed m; R lang

5、uage; nlme package随着网状 Meta 分析的发展及各种软件的研发，运算模型的选择也日益受到关注，广义线性包在 R 语言中实现5,6。nlme 程序包亦可以实现网状 Meta 分析6,7，本文仍以R 软件 R2WinBUGS 程序包在网状 Meta 分析中的应用一文中实例8,9 为例进行展示。1软件及程序包的安装/加载R 软件的安装已在前文讲述，本文使用的是R-3.0.1 8,10,11。此外，需要安装 nlme 程序包，具体命令为：in-模型（generalized linear m，GLM）就是其中之一。GLM 是线性模型的扩展，其特点是不强行改变数据的自然度量，数据可以具

6、有非线性和非恒定方差结构；与线性模型相比，GLM 模型中 Y 的分布可以是的指数分布（如高斯分布、泊松分布、二项式分布），联结函数可以是任何单调可微函数（如对数函数 logarithm 或逻辑函数 logit），这些优点使得 GLM 模型可处理多种变量，日益受到关注（欲详细了解 GLM 的读者建议阅读参考文献1,2）。nlme 是一款基于 S 语言在 S-PLUS 软件中实现stall.packages(“nlme”)。在弹出的框中选择某个镜像（CRAN）安装，安装完成后再由 library(“nlme”)命令完成加载。2 数据的加载与预处理2.1 数据加载R 软件 nlme 程序包实现网状

7、Meta 分析的基本思路是将因变量 y 设置为效应量，自变量 x 设置为干预措施，结合相关的干预措施之间的比较来计算得出线性关系7。本文将 y 设为logOR，x 设为 13混合效应模型（mixed-effects mo）分析的软件。因此，其可以同时进行线性和非线性混合效应模型（linear and nonlinear mixed effects ms，NLME）分析3,4。当前，混合效应模型也可通过 nlme 程序DOI: 10.7507/1672-2531.201400602012 年度重点课题（编号：基金项目：省教育科学“十二五2012A050），学院 2011 年度优秀中青年科技创新团

8、队项目（编号：2011CZX01），学院附属太和医院 2013 年循证医学种以及安慰剂，并对相关之间的比较作一“苗圃基金”项目（编号：EBM2013004）。作者简介：，男（1988 年），本科，住院医师，以循证医学、系统评价与 Meta 分析为主要研究方向。: zhangchao0803* 通讯作者: zengxiantao1128定标记。因实例中含有 3 臂试验，所以，首先将数据中CJEBM 355 © 2014 中国循证医学杂志编辑部MethodologyChin J Evid-based Med 2014, 14(3): 355-360所有的 3 臂试验拆分成 2 臂试验数据

9、进行排列（见表 1）。拆分之后，使用如下命令加载数据：表 88 个研究 组；c(“A”,“B”,“C”,“D”,“E”,“F”,“G”,“H”,“I”,“J”,“K”,“L”,“M”,“N”) 为列名（从 A 到 N），在此代表 13 种及安慰剂。linedata<-read.table("D:UsersAdministrator DesktopRworklinedata.txt", header=TRUE)。接着，设定配对变量矩阵，命令代码如下：for (aa in 1:88) for (bb in 6:7)matrixdataaa, as.numeric(char

10、ToRaw(as.cha racter(metadata$t1aa)-64<-1matrixdataaa, as.numeric(charToRaw(as.cha racter(metadata$t2aa)-64<-12.2数据转化在进行线性模拟之前，需用到 metafor 程序包10 来计算logOR 及方差。数据转化可分为三步。install.packages(“metafor”) 和第 一 步， 使 用 library(“metafor”) 命令安装和加载 Metafor 程序包10。第二步，计算因变量（即效应量）logOR。因实例为二分类数据，metafor 程序包使用的是

11、 escalc 命令：上述代码作用在于将两种治疗相互配对，metadata<- escalc(ai=r1, n1i=n1, ci=r2, n2i=n2, measure=“OR”, data=linedata)。将 metadata 试验组使用的配对的转换，相应的成 1 放置于 matrixdata 矩阵所对应的对照组设为 1。其中，as.character() 和 as.numeric() 是指将括号里面的元素分别转换成字符型和数字型，charToRaw() 为将字符转换成相应的ACS 码。2.4标识配对命令中：r1、n1、r2 和 n2 分别表示试验组的发生数与总数及对照组的发生数与

12、总数；measure 表示得出的效应量的对数如本例选“OR” 则结果自动转换为 log(OR)，但其结果是以变量名yi 来储存 log(OR) 数据集，而其方差集以变量名 vi 表示。需要注意的是，得出结果数据集将直接写入到原始数据中。第三步，转换变量名称。为方便记忆，本例将Meta 数据集中的 yi 及 vi 变量名称改为 logOR 及var，具体代码为：接着，需将 metadata 数据中的比较做相关的配对标识，以便后面 lme() 函数的使用。实现的命令为：mark<-paste(as.character(metadata$t1),as. character(metadata$t

13、2)。2.5预处理数据汇总在进行线性模型数据处理之前，我们需要搜集相关的预处理数据信息，具体命令为：totaldata<- cbind(metadata, matrixdata, mark)。3完成网状Meta分析在计算之前，可展示 totaldata 数据集（见表 1）， 命令为：totaldata。表 1 中，变量 logOR、A N 及mark 为最终 nlme 程序包中 lme() 函数所需的有效数据。接着，进行线性混合模型运算，命令如下：names(metadata)8:9<-c(“logOR”,“var”)。其中，8:9 代表第 8 行和第 9 行对应的数据集。2.3创

14、建哑变量矩阵本例我们需要将 88 个研究组中的 14 个自变量设定成 88×14 的哑变量矩阵，具体命令如下：matrixdata<-matrix(0,nrow=88,ncol=14,byrow= TRUE,dimnames=list(1:88,c(“A”,“B”,“C”,“D”,“E”,“F”,“G”,“H”,“I”,“J”,“K”,“L”,“M”,“N”)其中，0 为矩阵中的数值（在此全部设置为 0）； nrow 为矩阵行数；ncol 为矩阵列数；byrow 为依据行来对矩阵进行赋值；dimnames 为对矩阵的行和lme<-lme(logOR0+B+C+D+E+F+

15、G+H+I+J+K+L+M+N,random=1|mark,data= totaldata)。列进行命名；1:88 为行名（从 1 至 88），在此代命令中：logOR0+B+C+D+E+F+G+H+I+J+K+CJEBM 356 © 2014 Editorial Board of Chin J Evid-based Med中国循证医学杂志 2014, 14(3): 355360表 1 本文示例及预处理之后的数据123456789101112131415161718192021222324252627282930313233343536373839404142434445464748

16、49A A B A A B A A B A A A A A A D A A A A A E A A E A A E A A E A A A G A A G A A G A A A A A A AAB G G B L L B L L D D D D D E E E E E E F F E G G E K K E K K G G K K G N N G N N H J J K L L LL7373766666785555774039483661611521521501241241221211211201611221261211641647683837866667779792055214246132

17、741501541541221181181201191193241252491233151590.105 656 70.235 149 10.129 492 50.443 307 50.109 199 30.334 108 20.764 926 90.862 890 00.097 963 11.650 797 60.416 063 20.768 506 00.343 968 30.197 170 90.385 269 00.188 098 10.029 910 00.607 989 41.061 129 20.460 727 80.379 312 40.081 415 40.478 737 3

18、0.468 266 00.010 471 30.594 495 20.379 326 70.215 168 50.946 819 21.245 312 20.298 493 01.056 336 80.193 132 50.224 881 20.031 748 70.292 594 40.563 469 40.270 875 00.348 897 30.489 821 90.140 924 60.955 511 51.167 605 20.616 681 61.043 349 60.795 112 10.267 160 40.790 439 50.662 841 80.053 028 30.0

19、52 489 60.052 804 10.067 940 70.066 775 20.069 930 10.069 576 20.070 991 60.073 901 10.046 545 90.070 618 60.050 041 90.074 268 60.039 142 40.051 380 40.038 318 50.059 822 50.062 769 30.081 889 50.048 437 20.048 581 40.030 058 70.120 971 60.0.0.063 118 50.083 706 10.065 942 30.067 686 60.102 972 30.

20、083 416 80.086 333 90.286 847 40.285 839 40.0.081 583 00.081 627 30.078 387 10.080 772 70.0820.079 096 30.4990.2150.0.0.065 361 90.0.0570.065110110110111111011111011011011010000000000000000000000000000000000000000000000111110100000000000000000000011111101101101101000000000000000000000000000011000000

21、000000000000110000000000000000000011000000110110110000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011000000000000000000000000000000011011001100000000010000000011011000000000000000000000000000000000000111100000

22、000000000000000000000000000000000000000000000000000000000000000000000000000000001101100000000AB AG BG AB AL BL AB AL BL AD AD AD AD AD AE DE AE AE AE AE AF EF AE AG EG AE AK EK AE AK EK AG AG AK GK AG AN GN AG AN GN AH AJ AJ AK AL AL ALAL54492644441172424325151129414112618101031373745414152515141245

23、164943141139122137137273707070999919693931887819195410210210498981031842455612949150129556454117112112321515129595912663631323132324551515254549254124701977651411281232732742747033331969797188868628554555510410210210310010018399055129491491321101101110000000000000000000000000001101000000000000000000

24、10000000000001000CJEBM 357 © 2014 中国循证医学杂志编辑部sortt1t2r1n1r2n2logORvarABCD EFG H IJKLM N markMethodologyChin J Evid-based Med 2014, 14(3): 355-360续表 1505152535455565758596061626364656667686970717273747576777879808182838485868788A A B B B C C E E E F F F F G G G G G G G K G G G G G G G G I I J J

25、K L L LLL N G L M F H F F K G K L N I J K K K K L L L L N N N N N N K M K L M M N NN132937813387336681144941757559302767272657576435633135981539534746111424837415645116102611226312010813815123812323210798666110145509292961201445447170186161731391002078536072798226533593218331839415789146744735296730

26、2564707048733635811701074866511641484649564911195621266112010914014424011722710810066641024550969696118142554017119615373136100208255627584780.885 224 91.155 785 90.173 352 90.355 214 80.739 667 20.161 477 90.101 741 10.462 800 70.485 317 90.210 886 70.019 567 10.532 577 20.074 047 60.534 120 10.303

27、 682 40.042 471 40.028 987 50.287 682 10.080 042 70.205 444 00.502 695 50.297 251 50.510 009 00.307 667 40.340 587 00.875 468 70.413 661 90.343 771 50.479 965 30.789 526 30.188 518 30.407 306 91.185 623 70.154 150 70.336 472 20.580 629 00.354 138 90.196 710 30.328 779 90.0.090 858 80.0.077 791 60.0.

28、080 872 50.088 714 90.058 632 50.057 269 70.036 000 40.092 071 30.042 453 10.087 508 60.082 734 70.121 940 60.129 503 00.087 833 90.0.0.092 990 30.098 862 50.099 622 30.075 455 20.056 410 10.0.340 476 20.047 549 80.0810.056 762 80.115 886 10.058 335 40.082 050 50.514 520 20.100 367 80.384 523 80.0.1

29、15 518 00.114 906 80.104000000011100000000000000000000001011011110000000000000001000000010001111111011111000000100000000000000000000000000000000000000000000001000000000000000110000000000000000000000100000000000000001100000000000000101000011110100000000101010000100100000000100000001111000000000101111

30、000010000000000000000000000000010011000010000000000010000000000111111000000111AL AN BG BL BM CF CH EF EF EK FG FK FL FN GI GJ GK GK GK GK GL KL GL GL GN GN GN GN GN GN IK IM JK JL KM LM LN LNLN10000100000000000000000000000000110000000000000000000000000000000000000000000000000000000000000000000000000

31、0000000000000000000000000110000000000000000 sort 排序编号；t1、r1、n1分别为干预组的治疗方案、阳性例数及总例数；t3、r2、n2 分别为对照组的治疗方案、阳性例数及总例注数；治疗方案 t1 与 t2 中 A N 分别表示 13 种及安慰剂，A，placebo；B，bupropion；C，citalopram；D，desvenlafaxine；E，duloxetine；F， escitalopram；G，fluoxetine；H，fluvoxamine；I，mirtazapine；J，nefazodone；K，paroxetine；L，ser

32、traline；M，trazodone；N，velafaxine；此均为原始数据部分。 logOR 为计算之后各产生的研究组的效应量 OR 的对数；var 为计算之后各产生的研究组间的方差；此均为 metafor 程 序包的escalc() 函数转换而来。 A N 分别表示 13 种及安慰剂，与治疗方案中 A N 所代表的是一致的；其矩阵是由 matrix() 矩阵函数及配对赋值而得。 mark 为配对标记结果，由 paste() 函数合并而得。 对于上述数据中，对于本例 lme() 函数最终所需要的数据集为变量 logOR、A N 及 mark。L+M+N 成为固定部分，其“”左边代表因变量

33、（效此处，将“A”替换为“0”来进行运算。mark 为配对比较 的标记，同时作为 random 1 的随机效应范围。lindata 为储存数据的变量名。由于各种不同数据精准计算的需要，有时还需应量），“”右边代表自变量（13 种及安慰剂）；性模型中，对于存在 n 个自变量时，通常需要设置n1 个哑变量，另一个则需要用数字“0”替代，CJEBM 358 © 2014 Editorial Board of Chin J Evid-based Medsortt1t2r1n1r2n2logORvarABCD EFG H IJKLM N mark中国循证医学杂志 2014, 14(3): 35

34、5360对以下 3 个方面作相应的设置，本处均为默认设置如下： 算法（method）默认：限制性极大对数似然法（REML）； 权重（weight）默认：NULL，组间同方差（值得注意的是，通常在设置的时候，采取方差的协变量）；相关性（correlation）默认：NULL，组间不存在自相关性。上述代码运行之后，再运行汇总结果命令。命令为：summary(lme)。命令执行后的结果见框 1。上述结果中，是以“A”为基线作参考；框 1 R 软件 nlme 程序包进行网状 Meta 分析令执行结果Linear mixed-effects mo Data: linedatafit by REMLAIC

35、BIClogLik119.2982154.0605Random effects: Formula: 1 | mark-44.64908(Intercept)ResidualStdDev: 1.430292e-05 0.3715504Fixed effects: logor 0 + B + C + D + E + F + G + H + I + J + K + L + M + NValueStd.Error DFt-valuep-value 23232323232323232323232323B C D E F G H I J K L M N0.58385931.01646910.6272708

36、0.57521690.83459960.45941650.93511980.82744540.51190540.67008390.63474130.45332810.76018510.13539300.31750940.15262260.10136010.14054520.09021920.30696560.23653100.17177950.10518550.09614120.20834580.11378254.312331 0.00033.201383 0.00404.109947 0.00045.674985 0.00005.938300 0.00005.092224 0.00003.0

37、46334 0.00573.498253 0.00192.980014 0.00676.370499 0.00006.602176 0.00002.175845 0.04016.681039 0.0000Correlation: BCDEFGHIJKLMC D E F G H I J K L M N0.0290.0080.0730.0990.2700.0150.1350.0860.1630.3330.2600.1790.0110.0950.2950.0720.5170.0280.0230.070.0580.0280.0650.1110.0360.0220.0050.0100.0080.0310

38、.0130.0100.0150.3230.1990.0490.0930.0690.2810.1210.0880.1310.2440.1530.0940.0770.2370.1970.0950.2190.0370.2470.1950.4010.3670.2070.5050.0140.0120.0360.0300.0150.0330.0770.2750.1560.3580.1390.1920.1810.0810.1190.2490.2600.2310.2450.3410.137Standardized Within-Group Residuals:MinQ1MedQ3Max0.613952401.

39、92283637-2.76529166 -0.60957692 -0.05338441Number of Observations: 88 Number of Groups: 36CJEBM 359 © 2014 中国循证医学杂志编辑部MethodologyChin J Evid-based Med 2014, 14(3): 355-360当然，还可以使用数字“0”依次替换“B”至“N”，替析实现中的最大缺陷在于未考虑一致性、异质性及 协方差参数的设定，亦不具备检测功能；且有研究 显示使用非线性模型实现网状 Meta 分析，方法学仍在 。因此，介绍 nlme 程序包是必要的，该程换后

40、，即可得到相应的网状计算结果。有者可自行尝试。4 绘制图形4.1 网状关系图绘制的读模型之外的途径完善网状 Meta 分序从析提供新的可能途径。随着今后 nlme 程序包的不断完善，其在网状 Meta 分析中的运用将越发广泛。参 考 文 献同 R2WinBUGS 程序包一样，nlme 程序包当前尚无绘制网状关系图的功能，需要通过 network 程序包实现。有关使用此程序包及其他相关方法绘制 网状 Meta 分析关系图的方法，详见网状 Meta 分析中网状关系图的绘制一文12。4.2 相关图形绘制nlme 程序包可绘制的图形及命令分别如下（本文对这 5 种图形不一一展示）：（1） 截距汇总图：

41、plot(ranef(lme)（2） 残差分布图：plot(lme)（3） 按 mark 分组，标准化残差与拟合值关系图：plot(lme,resid(.,type=“p”)fitted(.)|mark,abline=0)（4） 按 mark 分 组，残差的箱型 图：plot(lme, markresid(.)（5） 按 mark 分组，效应量与拟合值关系图：plo t(lme,logorfitted(.)|mark,abline=c(0,1)5结语实际上，GLM 在网状 Meta 分析中的运用，虽1Dunteman GH, Ho MHR. An introduction to generalized linearms. Thousand Oaks, CA: Sage Publications, 2006.23(译). 广义线性模型. 上海: 格致Pinheiro JC, Bates DM. Mixed-Effects M USA, Ney Work: Springer, 2000., 2011.s in S and S-PLUS.4Tornøe CW, Agersø H, Jonsson EN,. Non-linear mixed-effectspharmacokinetic/ pharmacody