Network meta-analysis in R：gemtc包

传统meta分析与网状meta分析

不断增多的药物及治疗方案使医生和患者的选择增多,进而增加临床决策的困难,这就要求筛选出最有效最安全的治疗措施,提高卫生保健的质量。设计和实施良好的随机对照试验是评估临床效果和安全性的金标准,可以为临床治疗提供最佳证据,但许多不同干预措施的疗效差异尚无直接比较的随机对照试验证据,或相关研究数量较少、质量较低。传统Meta分析可以合并相关数据,提高统计效能和精确性, 但它仅比较2个干预措施间的疗效和安全性,无法比较同一疾病所有干预措施的疗效,且当2个干预措施暂无直接比较,就无法评估两者的有效性和安全性,因此传统meta分析为临床决策服务的能力是有限的。

网状meta分析(NMA)是传统meta分析的扩展,同时将一系列多个不同处理因素进行相互分析比较的方法。NMA最大的优越性在于它能将治疗同类疾病的不同干预措施(三个或三个以上)汇总后进行定量化的统计分析比较。调整间接比较和混合治疗效应是网状Meta分析的两种类型,调整间接比较即在没有直接比较时通过一个共同对照比较两个治疗措施间差异的统计学方法;当同时存在直接比较和间接比较时,可以将直接和间接比较结果进行合并，进而提高结果的精确性和统计学效能(混合治疗效应)。NMA还可以按照某一结果指标好坏进行排序,优选方案。

网状meta分析基本假设

NMA常涉及3个基本假设:同质性假设、相似性假设、一致性假设。同质性假设即直接比较的研究结果间不存在异质性,相似性假设即所有研究间影响效应量的因素相似,NMA结果不会产生偏倚,一致性假设是直接比较证据和间接比较证据同时存在,两种证据具有一致性。

用R实现网状meta分析

目前有很多用于实现NMA的软件,如ITC、WinBUGS、ADDIS、GeMTCR、SAS、State和R等。R软件广泛应用于Meta分析、诊断试验Meta分析和网状Meta分析（NMA）等，因其开放性和灵活性使其拥有广泛的程序包,NMA常用的程序包有gemtc、netmeta、pcnetmeta、nmalNLA、nmathresh等。本文使用gemtc包实现网状Meta分析。gemtc包基于贝叶斯广义线性模型，调用BUGS或JAGS软件进行网状Meta分析，处理的数据类型包括二分类变量、计数资料、连续变量、生存资料；分析模型有同质性假设模型、一致性假设模型、MCMC抽样、后验抽样控制、收敛性评估模型；统计假定方法有异质性评价、不一致性评价、模型拟合试验；统计推断和结果报告效应量及效应量评估。

---

使用gemtc实现网状meta分析

1安装及加载相关包（genmc包基于coda包，故需要先安装coda包，还有建议安装加载的包:testthat包、Matrix包、 XML包）

library(coda)
library(gemtc)

2 数据预处理

2.1读取数据

数据录入可读取GeMTC和ADDIS等软件整理的数据；或者在R窗口直接录入。

2.1.1 方法一：R窗口中录入数据

2.1.2 治疗方案的标识

treatments <- read.table(textConnection('
                                        id description
                                        A "treatment1"
                                        B "treatment2"
                                        C "treatment3"
                                        D "treatment4"'),
                         header = T)

2.1.3 按下列格式录入从文献中提取的疗效数据,所需录入的列依次为“Study”、“Treatment”、“Responders”、“Samplesize”，分别代表研究序号、治疗方案标识、阳性例数及总例数。

data <- read.table(textConnection('
study treatment responders sampeSize
1 A 155 253
1 B 63 133
2 A 72 137
2 B 42 135
2 C 63 140
3 A 89 175
3 B 72 172
3 C 48 103
3 D 50 102
4 B 49 130
4 C 23 132
5 B 110 267
5 D 49 134
6 C 29 100
6 D 44 91
7 A 48 105
7 C 59 106 
8 A 46 77
8 D 60 78
9 B 53 161 
9 C 82 155
9 D 68 125'),
header = T)

2.1.4 创建network数据

data为录入的数据,description为“研究名称”,treatments为录入的干预措施标识

network <- mtc.network(data,description = "研究名称", treatments = treatments)
print(network)

2.2 方法二：读取gemtc文件，并保存为net（本例读取数据为该程序包的示例数据）

sysfile <- system.file('extdata/luades-smoking.gemtc', package='gemtc')
net <- read.mtc.network(sysfile)

3 网状Meta分析

3.1 证据网络图

网络图中线的粗细与研究数成正比，线段越粗，比较这两个干预措施的研究数目越多。

plot(net)

3.2 内部关系总结

summary(net)

3.3 设置network model

type为是否选择一致性模型

model <- mtc.model(net,type = "consistency")
summary(model)

3.4 运行结果

sampler为调用软件的选取,这里调用JAGS实现分析,n.adapt为预迭代次数,n.iter为迭代运算次数，thin为步长

NMA分析结果

summary(result)

Results on the Log Odds Ratio scale

Iterations = 5001:25000
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean     SD  Naive SE Time-series SE
d.A.B 0.4931 0.4047 0.0014309       0.002578
d.A.C 0.8383 0.2399 0.0008482       0.001988
d.A.D 1.0995 0.4422 0.0015632       0.003501
sd.d  0.8437 0.1855 0.0006557       0.002370

2. Quantiles for each variable:

         2.5%    25%    50%    75% 97.5%
d.A.B -0.2941 0.2270 0.4902 0.7517 1.312
d.A.C  0.3846 0.6793 0.8305 0.9898 1.337
d.A.D  0.2579 0.8045 1.0892 1.3797 2.009
sd.d   0.5510 0.7118 0.8198 0.9498 1.274

-- Model fit (residual deviance):

    Dbar       pD      DIC 
54.15844 43.96942 98.12786 

50 data points, ratio 1.083, I^2 = 10%

3.5 绘图

3.5.1诊断收敛性

Brooks-Gelman-Rubin诊断图

Brooks-Gelman-Rubin诊断图为一种结合基于区间的图形评估和规模缩减因子(PSRF)定量分析的诊断模型收敛程度。图形评估方法为每次迭代均计算出缩减因子（shrinkfactor）的中位值及97.5%值,经过n次迭代后形成曲线并观察曲线是否相互拟合并持续稳定。将所有MCMC链最后一次迭代的结果进行分析,将链中与链间的变量相互对比。缩减因子（shrink factor）的中位值及 97.5% 值迅速趋向于1并达到稳定,且PSRF值趋向于1, 提示模型收敛满意。

par(mfrow=c(2,2))
gelman.plot(result)

这里n.adapt设为5000,n.iter设为20000,结果显示经过5000次预迭代后缩减因子的中位值及97.5%值迅速趋向于 1 并达到稳定,且PSRF值为1,提示模型收敛程度满意。

轨迹密度图

轨迹图为一种在迭代计算过程中表现MCMC链（Markov Chain Monte Carlo chain，MCMCchain,马尔科夫链-蒙特卡罗链）波动过程的图,用于诊断模型的收敛程度（convergence）。通过轨迹图形可了解各条在计算过程中是否达到稳定及良好的重叠,以判断模型收敛程度是否满意。

密度图作用与轨迹图一致,均用于诊断模型的程度,为一种表现参数的后验数值的分布形态图。密度图的评估方法为在预先设定分布的基础上，经过数值模拟计算后，观察后验数值的分布形态是否与预设分布一致。Bandwidth值可作为定量评估,其值越小,表示参数后验数值的分布范围与预设分布范围的差距越小。

这里,轨迹图预迭代次数设定为5000,及迭代次数设定为20000,每条MCMC链从起始部分已达到稳定融合,重叠面积占链波动范围的大部分,肉眼不能识别单条链的波动;密度图分布形态呈正态分布,收敛程度满意。

3.5.2 forest图

forest(result)

由森林图结果,我们可以认为四种方案效果由优至劣的顺序为D>C>B>A。

3.5.3 排序概率表和排序概率图

排序概率表（rank probabilities）由行数与列数等于总干预数的表格组成，首行为排序位数,首列为干预措施,表格中的数据为排序概率,表示干预措施排列在第n位的概率。

排序概率图（rankogram）以柱状图或曲线图形式表示各干预的排序概率更直观评估各干预措施的排列顺序。

概率排序表

ranks <- rank.probability(result)
print(ranks)

Rank probability; preferred direction = 1
       [,1]      [,2]      [,3]      [,4]
A 0.0000375 0.0025250 0.1066500 0.8907875
B 0.0592875 0.1753375 0.6607125 0.1046625
C 0.2325375 0.5999125 0.1672750 0.0002750
D 0.7081375 0.2222250 0.0653625 0.0042750

这里,A方案排序为4的概率为89%,B方案排序为3的概率为66%,C方案排序为2的概率为60%,D方案排序为1的概率为71%,因此,我们可以认为四种方案效果由优至劣的顺序为D>C>B>A,下图,排序概率柱状图也直观地展示此排序结果。 概率排序结果与森林图。

概率排序图

barplot(t(ranks), beside = TRUE,
       col = c("lightgreen", "mistyrose", "lightcyan","lavender"))
title(main = "Smoking cessation rates", font.main = 8)

4 参考文献

[1]张超,董圣杰,曾宪涛. R软件gemtc程序包在网状Meta分析中的应用[J]. 中国循证医学杂志,2013,13(10):1258-1264.

[2]易跃雄,张蔚,刘小媛,张娟,朱定军,吕琼莹. 网状Meta分析图形结果解读[J]. 中国循证医学杂志,2015,15(01):103-109.

[3]李伦,田金徽,姚亮,陈耀龙,杨克虎. 网状Meta分析的统计学基础、假设和证据质量评估[J]. 循证医学,2015,15(03):180-183.

[4]Gert van Valkenhoef, Joel Kuiper.gemtc程序包使用手册.https://CRAN.R-project.org/package=gemtc