本教程将介绍如何使用R的gemtc包对分类变量进行贝叶斯网状Meta分析。

前提条件：

需要下载R软件，本教程使用的R版本为3.5.3，URL：

https://cloud.r-project.org/bin/windows/base/old/3.5.1/R-3.5.1-win.exe

以及RStudio，一个R的友好交互界面，URL：

https://download1.rstudio.org/desktop/windows/RStudio-1.2.1335.exe

因为gemtc包现在只支持调用JAGS，所以还需要下载JAGS的exe文件

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

以上文件下载后直接双击安装即可。

准备工作完成后，打开RStudio后，在console里输入:

install.packages("gemtc")

第一次安装的时候，安装速度可能比较慢，耐心等待一下。

等安装完成后，设置我们的工作目录，后面可以把待分析的数据放入这个工作目录，敲入：

setwd("G:\\R网状教程\\分类变量")  #你可以更改为自己的工作目录，在后面环境变量增多的时候，在这一步前可以加一个rm(list=ls())清除环境变量。

然后敲入：

library(gemtc)     #装载gemtc包

将我们的数据准备成如下格式(与stata做网状的格式相同)：

其中每一行代表研究的一个臂，study代表研究，treatment代表干预名称，responders代表每一个臂的反应人数，sampleSize代表每一个研究中每一个臂的样本量。注意，以上英文单词必须准确无误，大小写都要一致，不然等会分析无法进行。

大家可以在excel里准备要分析的数据，然后另存为csv文件(以逗号为分隔符的表格文件)，在这里，我把数据文件存为data.csv(你可以更改为自己想要的名字)。然后把这个文件放到我们的工作目录之下。

把刚刚准备好的数据导入到R里，敲入：

data "data.csv",sep=",", header=T)   #sep代表逗号为分隔符，header=T代表把首行的数据作为标题，

把数据弄成gemtc的网络格式，敲入：

network #这里

输入：

plot(network)

可以看到网状图如下：

当然了，这些大小、颜色、字体等等的参数都是可以调整的，但是跟stata做网状图相比还是太麻烦了。Stata几乎不用调整一步到位，大家可以参考我之前制作的stata教程做网状图。

然后开始建立模型：

modeltype="consistency", n.chain=4,likelihood="binom",link="logit",linearModel="random")

其中，圆括号里的第一个参数为我们刚刚新建的gemtc网络类型的变量network，第二个参数type="consistency"，是指我们的模型为一致性模型，n.chain指的是马尔科夫蒙特卡洛MCMC的链条数目，一般为2-4均可。因为我们进行的是分类变量(OR)的网状分析，所以：likelihood="binom",link="logit"。如果我们用RR作为效应量，那么把link="log"。linearModel="random"指的是我们使用的模型为随机效应，当然也可以指定fixed作为固定效应模型。这一步可能会跳出一些红色的warning，强迫症的同学看到可能会感到不安，但不用紧张，不影响后续的分析，不管它即可。

现在就可以进行迭代了，敲入：

results <- mtc.run(model, n.adapt =50000, n.iter = 100000, thin = 1)

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

敲入：

summary(results)

可以看到如下结果：

其中结果包括DIC，DIC可以进一步与不一致性模型的DIC作比较，如果相差在5以内，就说明数据基本符合一致性的前提。DIC的绝对值没有任何意义，只用于相对比较。以及I2的总体异质性参数。

敲入：

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

可以看到每个干预相对于NT(no treatment)的森林图：

敲入：

plot(results)  #可以看到每个指标的收敛情况以及密度图

敲入：

gelman.plot(results)   # 收敛性诊断

进行排序，敲入：

ranks <- rank.probability(results,preferredDirection=1)  #若效应量越小越好，则preferredDirection=-1

print(ranks)

制作联赛表：

a <- relative.effect.table(results)

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

以上步骤可以在工作目录中导出这个表格。

不一致模型

很多同学不知道如何做不一致模型，想当然的以为将type="inconsistency"就OK了，其实并非如此。而应是在模型那一步中将type=”ume”或者”use”。至于ume和use具体是什么，本文不再赘述，感兴趣的同学可以查看交流群群文件中TSD4(inconsistency)的文件。如下：

modelume type="ume", n.chain=4,likelihood="binom",link="logit",linearModel="random")

resultsume <- mtc.run(modelume, n.adapt= 20000, n.iter = 50000, thin = 1)

summary(resultsume)    #查看不一致模型的DIC，然后与一致性模型进行比较。

用节点劈裂法探索局部不一致：

resultnodesplit <-mtc.nodesplit(network)

这里当然也可以写成复杂版的形式：

resultnodesplit <-mtc.nodesplit(network, n.adapt = 20000, n.iter= 50000, thin = 1, n.chain=4,likelihood="binom",link="logit",linearModel="random")

里面的参数都可以进行设定。在R的gemtc包里，我们无需像在WINBUGS软件里对初始值进行繁琐的设定，这一点是不是很棒。

print(summary(resultnodesplit))   #查看节点劈裂法的结果

plot(summary(resultnodesplit))

异质性分析：

resultanohe <- mtc.anohe(network)

这里当然也可以写成复杂版的形式：

resultanohe <- mtc.anohe(network, n.adapt = 20000, n.iter =50000, thin = 1, n.chain=4,likelihood="binom",link="logit",linearModel="random")

print(summary(resultanohe))

plot(summary(resultanohe))

至此，主分析的步骤基本已经完成。

如何把在RStudio中产生的图片导出来，可参考：手把手教你用R的gemtc包对连续变量进行贝叶斯网状Meta分析

其他教程：

贝叶斯网状Meta常用模型整理集合

手把手系列教程：将OpenBugs结果导入Stata作图SUCRA

(图片高清版)手把手教你用Stata的Network包实现连续性变量网状Meta分析

用贝叶斯方法实现网状meta中的亚组分析

手把手教你用R的gemtc包对连续变量进行贝叶斯网状Meta分析

我们的QQ交流群：519864737

更多教程，敬请关注：