R 软件 gemtc 程序包在网状 Meta 分析中的应用 张超 1 董圣杰 2 曾宪涛 1,3* 1. 湖北医药学院附属太和医院循证医学中心 ( 湖北十堰 442000);2. 苏州大学附属第一医院骨科 ( 江苏苏州 215006); 3. 湖北医药学院附属太和医院口腔科 ( 湖北十堰 442000) 摘要 R 软件相关程序包对于网状 Meta 分析的理论原理大多是基于贝叶斯统计模型, 也有少部分使用广义线性模型 R 软件 gemtc 程序包通过调用对应的 rjags 程序包 BRugs 程序包或 R2WinBUGS 程序包 ( 软件依次为 JAGS OpenBUGS 及 WinBUGS) 来执行网状 Meta 分析 同时, 还可以为 GeMTC 软件生成相关的数据存储文件 其本质是调用基于 MCMC 法的软件进行网状 Meta 分析 本文简要介绍 R 软件 gemtc 程序包调用 OpenBUGS 软件进行网状 Meta 分析 关键词 R 软件 ; 网状 Meta 分析 ;gemtc 程序包 Application of GeMTC R Package of R Software for Network Meta-Analysis ZHANG Chao 1, DONG Sheng-jie 2, ZENG Xian-tao 1,3* 1. Center for Evidence-Based Medicine and Clinical Research, Taihe Hospital, Hubei University of Medicine, Shiyan 442000, China; 2. Department of Orthopaedic, The First Affiliated Hospital of Soochow University, Suzhou 215006, China; 3. Department of Stomatology, Taihe Hospital, Hubei University of Medicine, Shiyan 442000, China Abstract The theoretical foundation of relevant packages of R software for network meta-analysis is mainly based on Bayesian statistical model and a few of them use generalized linear model. Network meta-analysis is performed using GeMTC R package through calling the corresponding rjags package, BRugs package, or R2WinBUGS package (namely, JAGS, OpenBUGS, and WinBUGS software, respectively). Meanwhile, GeMTC R package can generate data storage files for GeMTC software. Techonically, network meta-analysis is performed through calling the software based on Markov Chain Monte Carlo method. In this article, we briefly introduce how to use GeMTC R package to perform network metaanalysis through calling the OpenBUGS software. Key words R software; Network meta-analysis; GeMTC package 随着网状 Meta 分析 [1] 的发展, 关于网状 Meta 分析的软件日益增多 R 软件因其开放性和灵活性使其拥有广泛的程序包, 如 nlme 程序包 R2Win- BUGS 程序包 gemtc 程序包等 gemtc 程序包是通过相关指令来调用基于马尔科夫链 - 蒙特卡罗 (Markov Chain Monte Carlo,MCMC) 法的相关软件, 如 JAGS 软件 OpenBUGS 软件及 WinBUGS 软件来实现网状 Meta 分析 现已有介绍 nlme 程序包 [2] 和 R2WinBUGS 程序包 [3] 的文献, 本文对 R 软件 gemtc 程序包 [4,5] 在网状 Meta 分析中的应用进行简介, 调用的是 OpenBUGS 软件 DOI: 10.7507/1672-2530130215 基金项目 : 湖北省教育科学 十二五 规划 2012 年度重点课题 ( 编号 : 2012A050), 湖北医药学院 2011 年度优秀中青年科技创新团队项目 ( 编号 :2011CZX01), 湖北医药学院附属太和医院 2013 年循证医学 苗圃基金 项目 ( 编号 :EBM2013004) 作者简介 : 张超, 男 (1988 年 ~), 本科, 住院医师, 以循证医学 系统评价与 Meta 分析为主要研究方向 Email: zhangchao0803@126.com * 通讯作者,Email: zengxiantao1128@163.com 1 安装及加载 1.1 软件安装 R 软件当前最新版本为 R-3.0.1, 用户可从官方网站 http://www.r.project.org 上获取最新的版本及相应的统计包 [6] OpenBUGS 当前最新版本为 3.2.2, 下载网址为 http://www.openbugs.info/w.cgi/downloads 从官方网站下载软件后, 安装到本地计算机 安装完毕后,R 软件需要进一步安装和加载所需的程序包 程序包加载软件安装完毕后, 双击桌面的 R 软件图标, 可启动 R 的交互式窗口 (R GUI) 于命令提示符 > 后输入命令 install.packages( gemtc ), 选择镜像安装 (CRAN), 安装完成后可由 library( gemtc ) 命令完成加载 注意在加载 gemtc 包之前, 还需安装 lattice 及 coda 程序包, 具体方法如上 1258
中国循证医学杂志 2013, 13(10): 1258~1264 2 数据处理 我们先前介绍了命令读取和手动输入两种数据 录入方式 [3] 但 gemtc 程序包的数据录入方式有所 不同 本文同样引用 R 软件 R2WinBUGS 程序包在 网状 Meta 分析中的应用 文中 [3] 的实例为例, 具 体操作方法如下 : 2.1 治疗方案标识首先, 我们需要对治疗方案进行相关标记, 此处我们需要用到 read.table() 命令和 textconnection() 选择项, 对数据进行读取并进行关联 我们所引用的实例中的 14 种药物的对应关系如下 : A-placebo,B-bupropion,C-citalopram,Ddesvenlafaxine,E-duloxetine,F-escitalopram,Gfluoxetine,H-fluvoxamine,I-mirtazapine,Jnefazodone,K-paroxetine,L-sertraline,M-trazodone,Nvelafaxine 其在 R 软件运用 gemtc 程序包的录入标记命令如下 : treatments <- read.table(textconnection( 'id description A placebo B bupropion C citalopram D desvenlafaxine E duloxetine F escitalopram G fluoxetine H fluvoxamine I mirtazapine J nefazodone K paroxetine L sertraline M trazodone N velafaxine' ),header=true) 需要注意的是 :textconnection( ) 里面为英文模式的单引号,id 标识列与 description 标识列之间均以空格键一一隔开, 如 A 与 placebo 间空一格 2.2 数据录入同样, 以上述相同方式进行录入 其所需录入的列依次为 Study Treatment Responders Sample size, 分别代表研究组 治疗方案 阳性组数及总例数 具体数据及格式见表 1 2.3 network 数据创建 2013 中国循证医学杂志编辑部 论著 方法学研究与常用的数据存储不同,gemtc 程序包其数据的存储方式较为严格, 其采用的方式为 GeMTC 软件 [4] 的数据识别存储模式 对于 network 数据创建需要运用 mtc.network() 命令来进行执行, 命令如下 : network<-mtc.network(data, description ="Example", treatments=treatments) data 为录入的研究数据,description 为描述研究的名称,treatments 为录入的治疗方案标识 如此, 一个 network 数据已成功创建 2.4 生成 gemtc 格式数据 R 软件还具有另一项功能, 就是可以为那些数据录入繁琐的软件生成相关格式数据存储文件, 从而为这些软件的数据录入提供方便 gemtc 程序包可以为 GeMTC 软件生成相关的.gemtc 格式存储文件, 具体命令如下 : write.mtc.network(network,"c:/users/administrator/desktop/rwork/file.gemtc") 其中 :network 为事先临时存放的数据集 因默认保存路径在软件深层, 难以找到, 建议以 C:/ Users/Administrator/Desktop/Rwork/file.gemtc 路径保存, file.gemtc 为储存数据文件 2.5 network 数据读取 network 数据读取与以往的读取不同, 其需要命令 read.mtc.network() 来读取, 具体命令代码如下 : network<-read.mtc.network("c:/users/administrator/desktop/rwork/file.gemtc") read.mtc.network() 为使用 network 数据格式来进行读取 需要说明的是, 本例至 2.4 已经完成 2.4 中介绍的数据保存, 若再次仍采用 R 读取, 还可以使用.txt 格式 2.5 是介绍再次使用 R 读取采用 2.4 保存的文件源的方法 3 实现分析 gemtc 程序包自身并没有运算程序功能, 需要 model 进行相关设置及需要相关指令调用其他网状 meta 分析软件来进行相关的迭代运算 3.1 设置 network model 在选取具有网状 meta 分析功能软件之前, 我们需要对 model 进行相关设置, 需要 mtc.model() 命令, 具体命令如下 : model <- mtc.model(network,type = "Consistency", n.chain = 3) 1259
表 1 本文使用的数据及排列方式 Study Treatment Responders Sample size Study Treatment Responders Sample size Study Treatment Responders Sample size 1 A 73 152 19 G 52 103 42 I 35 66 1 B 76 150 19 N 54 100 43 G 27 61 1 G 83 154 20 A 5 18 43 J 29 64 2 A 66 124 20 H 9 18 44 G 67 101 2 B 78 122 21 A 15 42 44 K 67 102 2 L 66 118 21 J 25 39 45 G 27 45 3 A 55 121 22 A 14 45 45 K 30 45 3 B 77 120 22 J 41 90 46 G 26 50 3 L 79 119 23 A 12 56 46 K 25 50 4 A 40 161 23 K 24 55 47 G 57 92 4 D 205 324 24 A 45 129 47 K 64 96 5 A 39 122 24 L 70 129 47 L 70 96 5 D 52 125 25 A 16 49 48 G 35 120 6 A 48 126 25 L 19 49 48 L 48 118 6 D 142 249 26 A 49 150 49 G 63 144 7 A 36 121 26 L 77 149 49 L 73 142 7 D 46 123 27 A 43 129 50 G 31 54 8 A 61 164 27 L 65 132 50 N 36 55 8 D 132 315 28 A 13 116 51 G 35 47 8 E 74 159 28 L 26 111 51 N 35 40 9 A 54 141 29 A 29 102 52 G 98 170 9 E 55 141 29 N 53 95 52 N 81 171 10 A 49 139 30 B 37 61 53 G 153 186 10 E 64 128 30 G 35 62 53 N 170 196 11 A 26 122 31 B 81 122 54 G 95 161 11 E 54 123 31 L 93 126 54 N 107 153 12 A 44 137 32 B 33 63 55 G 34 73 12 E 117 273 32 M 21 61 55 N 48 73 12 F 112 274 33 C 87 120 56 I 74 139 13 A 24 70 33 F 83 120 56 K 66 136 13 E 32 70 34 C 33 108 57 I 61 100 13 G 15 33 34 H 31 109 57 M 51 100 14 A 51 99 35 E 66 138 58 J 11 20 14 E 129 196 35 F 83 140 58 K 16 20 14 K 59 97 36 E 81 151 59 J 42 78 15 A 41 93 36 F 94 144 59 L 41 82 15 E 126 188 37 E 144 238 60 K 48 53 15 K 63 86 37 K 157 240 60 M 48 55 16 A 18 78 38 F 94 123 61 L 37 60 16 G 132 285 38 G 89 117 61 M 46 62 17 A 10 19 39 F 175 232 62 L 41 72 17 G 31 54 39 K 146 227 62 N 49 75 17 K 32 55 40 F 75 107 63 L 56 79 18 A 37 102 40 L 74 108 63 N 56 84 18 G 45 104 41 F 59 98 64 L 45 82 18 N 51 102 41 N 47 100 64 N 49 78 19 A 41 98 42 G 30 66 A:placebo;B:bupropion;C:citalopram,;D:desvenlafaxine;E:duloxetine;F:escitalopram;G:fluoxetine;H:fluvoxamine;I:mirtazapine; J:nefazodone;K:paroxetine;L:sertraline;M:trazodone;N:velafaxine 1260
中国循证医学杂志 2013, 13(10): 1258~1264 其中,network 为 network 数据,type 为是否选取一致性模型,n.chain 为迭代运算中链的条数 需要注意的是, 目前该程序包通过 R 软件只能使用一致性模型, 非一致性模型 (inconsistency model) 尚不能使用 3.2 软件选取 gemtc 程序包提供了 3 种软件的调用方式, 依次为 JAGS OpenBUGS 及 WinBUGS 软件, 通过 mtc.run() 命令选取并执行 对于调用软件的选取可通过设置参数 sampler 值来实现, 其取值 论著 方法学研究 rjags BRugs R2WinBUGS 依次为调用 JAGS OpenBUGS 及 WinBUGS 软件 当 sampler 设为 NA 时, 则为依次调用上述三种软件进行迭代运算 上述操作计算机将自动完成建模 数据加载及初始化等过程 需要注意的是 : 调用三种软件的前提是计算机系统盘上必须安装相应的软件及在 R 软件拥有相关的程序包 由于 WinBUGS 软件已停止更新, 对于有些复杂数据可能存在无法处理的情况, 相比之下, 笔者建议最好能够选用 OpenBUGS 软件 当然, >results <-mtc.run(model, sampler ="BRugs", n.adapt =, n.iter =, thin = 1) Welcome to BRugs connected to OpenBUGS version 3.2.2 model is syntactically correct data loaded model compiled [1] "D:\\Users\\ADMINI~1\\AppData\\Local\\Temp\\RtmpuCqg5W/inits1.txt" [2] "D:\\Users\\ADMINI~1\\AppData\\Local\\Temp\\RtmpuCqg5W/inits2.txt" [3] "D:\\Users\\ADMINI~1\\AppData\\Local\\Temp\\RtmpuCqg5W/inits3.txt" Initializing chain 1: initial values loaded and chain initialized but another chain contain uninitialized variables Initializing chain 2: initial values loaded and chain initialized but another chain contain uninitialized variables Initializing chain 3: model is initialized model is already initialized updates took 12 s deviance set monitor set for variable 'd.a.b' monitor set for variable 'd.a.d' monitor set for variable 'd.a.e' monitor set for variable 'd.a.f' monitor set for variable 'd.a.g' monitor set for variable 'd.a.h' monitor set for variable 'd.a.j' monitor set for variable 'd.a.k' monitor set for variable 'd.a.l' monitor set for variable 'd.a.n' monitor set for variable 'd.b.m' monitor set for variable 'd.f.c' monitor set for variable 'd.g.i' monitor set for variable 'sd.d' updates took 52 s 图 1A OpenBUGS 行网状 Meta 分析的迭代进程 2013 中国循证医学杂志编辑部 1261
选择 WinBUGS 软件仍是可以的 3.3 实现分析本例以调用 OpenBUGS 为例作为演示, 具体命令如下 : results <-mtc.run(model, sampler ="BRugs", n.adapt =, n.iter =, thin = 1) n.adapt 为预迭代次数, n.iter 为迭代运算次 数, thin 为步长 命令运行的过程及结果见图 1 4 绘制图形在完成上述命令操作后, 我们基于 R 软件绘图功能, 对数据结果进行相关图形的绘制 需要说明的是, 以下截取的图形均可采用 R 软件保存图形的方式选择对应的位置及格式进行保存 Iterations = 5001: Thinning interval = 1 Number of chains = 3 Sample size per chain = 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naïve SE Time-series SE d.a.b 0.5884 0.16065 0.0006558 0.0017701 d.a.d 0.6877 0.15163 0.000619 0.0012211 d.a.e 0.5802 0.11339 0.0004629 0.0011434 d.a.f 0.8756 0.14932 0.0006096 0.0015964 d.a.g 0.5003 0.11076 0.0004522 0.0015197 d.a.h 0.951 0.45808 0.0018701 0.0066882 d.a.j 0.6594 0.22344 0.0009122 0.0026155 d.a.k 0.6476 0.13466 0.0005497 0.001661 d.a.l 0.6771 0.10881 0.0004442 0.0012514 d.a.n 0.7617 0.13254 0.0005411 0.0015976 d.b.m -0.1263 0.27025 0.0011033 0.0032112 d.f.c 0.1689 0.35278 0.0014402 0.0043635 d.g.i 0.3428 0.25177 0.0010279 0.0024953 sd.d 0.2565 0.04943 0.0002018 0.0007479 2. Quantiles for each variable: 2.50% 25% 50% 75% 97.50% d.a.b 0.27374 0.48075 0.5876 0.6954 0.9082 d.a.d 0.38536 0.58868 0.6886 0.7877 0.9863 d.a.e 0.35769 0.50468 0.5806 0.6564 0.8033 d.a.f 0.58434 0.77492 0.8745 0.9752 1.1732 d.a.g 0.2822 0.42603 0.5009 0.5742 0.7173 d.a.h 0.05157 0.64703 0.9503 561 1.8546 d.a.j 0.22627 0.50856 0.6572 0.8091 981 d.a.k 0.38647 0.55701 0.6469 0.7359 0.9172 d.a.l 0.46328 0.60387 0.677 0.7503 0.8914 d.a.n 0.50241 0.67242 0.7618 0.8498 219 d.b.m -0.65977-0.3072-0.1262 0.0563 0.406 d.f.c -0.52389-0.06604 0.1682 0.4039 0.8602 d.g.i -0.14992 0.17446 0.3421 0.5108 0.8408 sd.d 0.16513 0.22291 0.2542 0.2884 0.359 图 1B 网状 Meta 分析结果 1262
1.30 0 1.10 0 2 4 6 8 1.10 1.12 0 1 d.a.n 1.6 1 1 1 sd.d 1.12 d.g.i 1 d.a.l 1.6 1 1.8 1 0 1.1 1.3 1.5 1.6 d.f.c 1.6 1.8 d.b.m 1.1 1.3 1.5 1.6 1 d.a.k 1 1.1 1.3 1.5 1.6 1.7 1.6 1.8 2.0 1 d.a.j 1.3 2.0 d.a.h 1.1 1 d.a.g 8 d.a.f 4 1 0 2 4 6 8 1.10 1.12 d.a.e 1.3 1.1 0 5 1.10 1.15 0 5 1.30 d.a.d 1.5 d.a.b 著 方法学研究 论 中国循证医学杂志 2013, 13(10): 1258 1264 1 图2 1 收敛诊断图 Density of d.g.i 0.0 0.5 Trace of d.g.i 1 0.5 Iterations 0.0 0.5 1.5 N = Bandwidth = 0.0297 Density of sd.d 0 0.1 4 0.4 8 Trace of sd.d 1 0.1 Iterations 1.7996 (1.3149, 2.4798) 1.9910 (701, 2.6813) 1.7871 (300, 2.2329) 2.3977 (1.7938, 3.2324) 1.6503 (1.3260, 2.0489) 2.5865 (529, 6.3891) 1.9294 (539, 2.9983) 1.9095 (718, 2.5024) 1.9679 (1.5893, 2.4386) 2.1422 (1.6527, 2.7785) d.b.m 0.8814 (0.5170, 1.5008) d.f.c 1.1831 (0.5922, 2.3636) d.g.i 079 (0.8608, 2.3183) 1图 4 7 网状 Meta 分析森林图 2013 中国循证医学杂志编辑部 0.4 0.5 轨迹和密度图 Odds Ratio (95% CrI) d.a.b d.a.d d.a.e d.a.f d.a.g d.a.h d.a.j d.a.k d.a.l d.a.n 0.5 0.3 N = Bandwidth = 0.005735 图3 Comparison 0.2 E D F C G B H A N I J M K 图5 L 网状关系图 1263
> summary(network) $Description [1] "MTC dataset: Example" $`Studies per treatment` A B C D E F G H I J K L M N 29 6 2 5 11 8 22 2 3 5 13 17 4 13 $`Number of n-arm studies` 2-arm 3-arm 52 12 图 6 网状数据结构分布及内涵 4.1 收敛诊断图收敛诊断 (convergence diagnostics) 图绘制的命令为 gelman.plot(results), 即可绘制图形 此处应该有 13 个 gelman plot( 图 2), 可以采用如下命令进行全部显示 : par(mfrow=c(3,5)) gelman.plot(results,auto.layout =F) 4.2 轨迹和密度图轨迹 (trace) 和密度 (density) 图的绘制的命令为 plot(results), 如图 3 此处需要说明的是, 若要在此窗口上显示全部图形, 可通过将字体调小的方式实现 4.3 森林图森林图 (forest) 绘制的命令为 forest(results), 如图 4 4.4 网状关系图在完成 2.3 的步骤后, 即可进行网状关系图 (network plot) 的绘制 绘制命令为 plot(network), 如图 5 为进一步了解此图的内涵, 使用 summary (network) 查看其内部联系, 如图 6 如, 从图中可以看出, 两臂研究有 52 个, 三臂研究有 12 个 5 结语目前实现网状 Meta 分析的软件中, 进行数据 处理模型多数采取贝叶斯模型 gemtc 程序包就是基于贝叶斯模型来实现数据处理的, 与 WinBUGS OpenBUGS JAGS 等相比, 其省略了建立相关的贝叶斯 model 的过程, 甚至操作者不会使用调用的软件都可以得到相关结果 对于那些操作 BUGS 等软件欠佳操作者而言, 该程序的便利及灵活操作, 为他们在处理网状 Meta 分析数据提供了可能, 同时, 其也可以生成.gemtc 格式数据存储文件, 再通过 GeMTC 软件来实现网状 Meta 分析, 这也避免了该软件繁琐的数据录入工作, 为操作者节约了宝贵的时间 此外,GeMTC 软件 [4] 亦可独立进行网状 Meta 分析 当前,gemtc 程序包不仅适合作二分类型数据的网状 Meta 分析, 还适合制作连续数据的网状 Meta 分析 但仅基于一致性模型, 且仍需要进一步优化, 对于非一致性模型尚待研发 此外, 绘制的网状关系图只能反映各研究组之间存在直接比较, 而无法反映研究组间的具体比较量 若需要显示比较量, 则需要通过其他程序包来实现 [7] 参考文献 1 曾宪涛, 曹世义, 孙凤, 等. Meta 分析系列之六 : 间接比较及网状 Meta 分析. 中国循证心血管医学杂志, 2012, 4(5):399-402. 2 张天嵩, 熊茜. 网络 Meta 分析在 R 软件中的实现. 循证医学, 2012, 12(3):185-188. 3 曾宪涛, 张超, 郭毅. R 软件 R2WinBUGS 程序包在网状 Meta 分析中的应用. 中国循证医学杂志, 2013, 13(9): 1137-1144. 4 van Valkenhoef G, Lu B, de Brock B, et al. Automating network meta-analysis. Research Synthesis Methods, 2012, 3(4): 285-299. 5 van Valkenhoef G, Tervonen T, de Brock B, et al. Algorithmic Parameterization of Mixed Treatment Comparisons. Statistics and Computing, 2012, 22(5): 1099-1111. 6 董圣杰, 曾宪涛, 郭毅. R 软件 Metafor 程序包在 Meta 分析中的应用. 中国循证医学杂志, 2012, 12(9): 1141-1147. 7 张超, 徐畅, 曾宪涛. 网状 Meta 分析中网状关系图的绘制. 中国循证医学杂志, 2013: In press. 收稿日期 :2013 06 28 修回日期 :2013 09 25 本文编辑 : 张永刚 杜亮 1264