用R语言进行连续型变量的网状meta分析--详细教程

最近在后台有不少读者留言想要了解一下网状meta分析，而我恰好前两天帮一位医生朋友处理了一次连续型变量的网状meta分析，那么本文我将以其中的部分数据为例，向大家详细展示一下此需求是如何用R软件实现的。需要本教程的原始数据及完整程序命令的可在公，众号（全哥的学习生涯）内回，复“网状meta”。

1. 网状meta分析的理论基础

网状meta分析是传统meta分析的拓展，可同时比较三个或三个以上干预措施的疗效，其主要包括调整间接比较和混合治疗分析。间接比较是指通过干预措施A vs C和干预措施B vs C的结果，间接得出A vs B的相对效果的一种方法。目前认为，Meta分析中进行间接比较的原因有二:一是无直接比较的原始研究;二是有直接比较的原始研究但这些研究数量较少或质量较低;混合治疗分析是在直接比较的基础上合并间接比较的证据，从而提高分析结果的精确性，但是这种方法主要是应用于干预措施可以形成具有闭合环路( loop)时。

网状meta分析的最大优势就是可以对治疗同类疾病的不同干预措施进行量化比较，并按照某一结果指标效果好坏进行排序，进而选择最优治疗方案。网状meta分析的基本假设包括同质性假设、相似性假设及一致性假设。其统计方法主要包括频率学法和贝叶斯方法。

2. 准备工作

2.1 安装与加载

本教程使用的R软件版本为3.6.2，RStudio版本为1.2.1335，需要软件的读者可在公，众号（全哥的学习生涯）内回，复“R”,”RS”自取。同时，我们还需要安装一款基于马尔科夫链-蒙特卡罗（Markov Chain Monte Carlo，MCMC）法的软件：JAGS，以完成对数据的建模、迭代运算等，需要本软件的伙伴可以在以下地址中下载：

https://nchc.dl.sourceforge.net/project/mcmcjags/JAGS/4.x/Windows/JAGS-4.3.0.exe

若下载缓慢，可直接在公，众号（全哥的学习生涯）内回，，复“JAGS”获取本软件。

同时，我们需要安装“gemtc”包，及加载此包。代码如下：

install.packages(“gemtc”)

library(gemtc)

2.2 数据处理

在本教程中，我们的目的是运用网状meta分析的方法对联用B药物与其他西药（包括A，C，D，E，F，G，H，L，M等药物）的疗效进行系统评价，以确定最优干预措施为单用Flunarizine或是联用其他药物。输入以下命令将数据录入进R中：

setwd(dir="C:/Users/Desktop") #可更改为自己的工作目录

data <- read.csv("data.csv",header = T)

原数据的前20行如图1所示：

图 1

在本数据中，每一行代表研究的一个臂（在本教程中选取的都是双臂实验，故一项研究拆分成两行），study代表研究，treatment代表干预名称，mean代表连续变量的均值，std.dev代表连续变量的标准差，sampleSize代表每一个研究中每一个臂的样本量。特别需要注意的是，由于分析软件基于文本识别系统，以上英文单词必须按以上描述所示，大小写都要一致，否则将出现运行错误。

2.3 生成gemtc格式数据

R软件还具有另一项功能，就是可以为那些数据录入繁琐的软件生成相关格式数据存储文件，从而为这些软件的数据录入提供方便。输入命令：

network <- mtc.network(data)

3. 实现分析

gemtc 程序包自身并没有运算程序功能，需要model 进行相关设置及需要相关指令调用JAGS分析软件来进行相关的迭代运算。

3.1 生成网状图

键入命令：

plot(network,mode="circle",displaylabels=TRUE,boxed.label=FALSE)

生成的网状图如图2所示：

图 2

可使用如下命令查看其内部联系：

summary(network)

如图3所示，从图中可以看出，双臂实验有48个。

图 3

3.2设置network model

在选取具有网状meta 分析功能软件之前，我们需要对model 进行相关设置，需要“mtc.model()”命令，具体命令如下：

model <- mtc.model(network,type = "consistency",n.chain = 4,likelihood=”normal”,link=”identity”,linearModel=”random”)

其中，network 为network 数据，type 为是否选取一致性模型，n.chain 为迭代运算中链的条数(马尔科夫蒙特卡洛MCMC链条)，一般取2-4。非一致性模型（inconsistency model）可将type调整为“ume”或“use”。likelihood=”normal”,link=”identity”,定性为连续变量。linearModel=”random”指选取的模型为随机效应模型，若填fixed则为固定效应模型。

3.3软件选取及迭代

gemtc 程序包提供了3 种软件的调用方式，依次为JAGS、OpenBUGS 及WinBUGS 软件，通过“mtc.run()”命令选取并执行。对于调用软件的选取可通过设置参数“sampler”值来实现，当取值“rjags”为调用JAGS软件，计算机将自动完成建模、数据加载及初始化等过程。具体命令如下：

results <- mtc.run(model,sampler = "rjags",n.adapt=20000,n.iter = 5000,thin=1)

其中，model是我们上一步建立的模型名称，n.adapt为退火次数，n.iter为迭代次数。这个时候我们耐心等待，RStudio界面如图4所示：

图 4

查看结果的命令为：

summary(results)

可看到图5及图6所示结果：

图 5

图 6

4.图形绘制与结果展示

4.1 森林图绘制

敲入命令：

forest(relative.effect(results,"Flunarizine"))

可以看到每个干预相对于B药物的森林图（图7）：

图 7

4.2 收敛诊断图

收敛诊断（convergence diagnostics）图绘制的命令为“gelman.plot(results)”，即可绘制图形。

图 8

此处应该有12 个gelman plot（图8），可以采用如下命令进行全部显示：

par(mfrow=c(3,4))

gelman.plot(results,auto.layout =F)

4.3轨迹和密度图

轨迹（trace）和密度（density）图的绘制的命令为“plot(results)”，如图9、图10。此处需要说明的是，若要在此窗口上显示全部图形，可通过将字体调小的方式实现。更多meta分析及数据分析知识，请关注公，众号：全哥的学习生涯。

图 9

图 10

4.4 排序图

进行排序，如图11：

图 11

ranks <- rank.probability(results)

sucra(ranks)

堆积排序图，如图12：

plot(ranks)

图 12

单个排序等级图，如图13：

plot(ranks,beside=TRUE)

图 13

4.5制作联赛表

输入命令：

a <- relative.effect.table(results)

write.csv(a,"1111111.csv")

可在工作路径下生成1111111的联赛表，如图14：

图 14

4.6 异质性检验