爱之初体验——网状meta分析必备技能（二）
R 软件风靡全球，R 软件的 gemtc 程序包通过调用对应的 rjags 程序包来执行网状 Meta 分析是最新的流行趋势。同时，还可以为 GeMTC 软件生成相关的数据存储文件。其本质是调用基于MCMC 法的软件进行网状 Meta 分析。
本文简要介绍 R 软件 gemtc 程序包调用rjags软件进行网状Meta。
废话少说，代码如下：
rm(list=ls())
getwd()
# 猴哥是放在F盘的0108gemtc_sm文件夹下，你可以自己新建一个文件夹
setwd("F:\\0108gemtc_sm")
# F:\0108gemtc-master
# 如果没有安装gemtc，请执行如下命令，执行时去掉#，并在Rstudio中按 ctrl+enter
# install.packages("gemtc")
#也可以是如下路径
# setwd("F:\\0108gemtc-master\\gemtc")
# setwd("F:\\0108gemtc-master\\gemtc\\tests")
library(lattice)
library(coda)
library(codetools)
# library(testthat)
library(gemtc)
library(R2OpenBUGS)
library(igraph)
library(meta)
# library(XML)
# igraph (>= 0.6.4), meta (>= 2.1), XML (>= 3.6)
# 导入治疗措施treatments 并编号A B C D E F 。。。
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)
# 默认有row.name 1 2 3 4.。。
# write.csv(treatments, file = "F:\\0108gemtc_sm\\treatments0001.csv")
# 导入数据data，为 study treatment responders sampleSize的arm水平数据格式--单臂长数据格式
data <- read.table(textConnection(
'study treatment responders sampleSize
1 A 73 152
1 B 76 150
1 G 83 154
2 A 66 124
2 B 78 122
2 L 66 118
3 A 55 121
3 B 77 120
3 L 79 119
4 A 40 161
4 D 205 324
5 A 39 122
5 D 52 125
6 A 48 126
6 D 142 249
7 A 36 121
7 D 46 123
8 A 61 164
8 D 132 315
8 E 74 159
9 A 54 141
9 E 55 141
10 A 49 139
10 E 64 128
11 A 26 122
11 E 54 123
12 A 44 137
12 E 117 273
12 F 112 274
13 A 24 70
13 E 32 70
13 G 15 33
14 A 51 99
14 E 129 196
14 K 59 97
15 A 41 93
15 E 126 188
15 K 63 86
16 A 18 78
16 G 132 285
17 A 10 19
17 G 31 54
17 K 32 55
18 A 37 102
18 G 45 104
18 N 51 102
19 A 41 98
19 G 52 103
19 N 54 100
20 A 5 18
20 H 9 18
21 A 15 42
21 J 25 39
22 A 14 45
22 J 41 90
23 A 12 56
23 K 24 55
24 A 45 129
24 L 70 129
25 A 16 49
25 L 19 49
26 A 49 150
26 L 77 149
27 A 43 129
27 L 65 132
28 A 13 116
28 L 26 111
29 A 29 102
29 N 53 95
30 B 37 61
30 G 35 62
31 B 81 122
31 L 93 126
32 B 33 63
32 M 21 61
33 C 87 120
33 F 83 120
34 C 33 108
34 H 31 109
35 E 66 138
35 F 83 140
36 E 81 151
36 F 94 144
37 E 144 238
37 K 157 240
38 F 94 123
38 G 89 117
39 F 175 232
39 K 146 227
40 F 75 107
40 L 74 108
41 F 59 98
41 N 47 100
42 G 30 66
42 I 35 66
43 G 27 61
43 J 29 64
44 G 67 101
44 K 67 102
45 G 27 45
45 K 30 45
46 G 26 50
46 K 25 50
47 G 57 92
47 K 64 96
47 L 70 96
48 G 35 120
48 L 48 118
49 G 63 144
49 L 73 142
50 G 31 54
50 N 36 55
51 G 35 47
51 N 35 40
52 G 98 170
52 N 81 171
53 G 153 186
53 N 170 196
54 G 95 161
54 N 107 153
55 G 34 73
55 N 48 73
56 I 74 139
56 K 66 136
57 I 61 100
57 M 51 100
58 J 11 20
58 K 16 20
59 J 42 78
59 L 41 82
60 K 48 53
60 M 48 55
61 L 37 60
61 M 46 62
62 L 41 72
62 N 49 75
63 L 56 79
63 N 56 84
64 L 45 82
64 N 49 78'
),header=TRUE)
# 默认有row.name 1 2 3 4....
# write.csv(data, file = "F:\\0108gemtc_sm\\data0001.csv")
# 3.1中的数据输入方式较为繁琐，大家习惯了在excel中处理数据，可以先转化为 .csv格式
# 读取CSV文件
# header=T,表示文件存在，表示需要读取数据的数据头部
# skip=0,表示数据中没有需要跳过的行, ./treatments0001.csv表示本文件夹中的treatments0001.csv文件，row.names = 1是第一列作为行名
treatments<- read.csv("./treatments0001.csv",header=T,sep = ",",skip=0,row.names = 1)
# 我们看一下数据
head(treatments)
View(treatments)
# 我们看一下数据
data <- read.csv("./data0001.csv",header=T,sep = ",",row.names = 1)
head(data)
View(data)
建立成network 单臂长数据格式,description="Example"是对network的描述
network<- mtc.network(data, description="Example", treatments=treatments)
#network是单臂长数据
network
# 进一步了解 network的属性
head(network)
# 数据内容如下：为strn格式
# MTC dataset: Example
# Arm-level data:
# study treatment responders sampleSize
# 1 1 A 73 152
# 2 1 B 76 150
# 3 1 G 83 154
#5####
#写入数据文件到Rwork文件夹里，写入gemtc文件中，方便我们GeMTC的软件的读取
# write.mtc.network(network,"C:/Users/Administrator/Desktop/file.gemtc")
# 新版gemtc包不能用，如果需要用，请安装旧版的gemtc包
# write.mtc.network(network, "./data0001.gemtc")
##1、读取系统数据
##file <- system.file("extdata/file.gemtc", package="gemtc")
##以network形式读取数据
##network <- read.mtc.network(file)
##1、1直接读取Rwork文件夹里文件到network里
# 需要说明的是，本例至 2.4 已经完成。2.4 中介绍的数据保存，若再次仍采用 R 读取，还可以使用 “.txt”格式。
# network <- read.mtc.network("C:/Users/Administrator/Desktop/Rwork/file.gemtc")
# network <- read.mtc.network("./data0001.gemtc")
#创建model中运算的属性
#model <- mtc.model(network)
#一致性模块
# 其中， network 为 network 数据， type 为是否选取一致性模型，
# n.chain 为迭代运算中链的条数。需要注意的是，目前该程序包通过 R 软件只能使用一致性模型，
# model <- mtc.model(network,type = "Consistency", n.chain = 3) #例子
model <- mtc.model(network,type = "consistency", factor = 2.5, n.chain = 3) #ok
#构建模型成功，出现warning，别害怕
#网状图在可这里可以直接出图
#模型网状图
# plot(network)
# # points(network, cex = .5, col = "dark red") #改变颜色没有用
# #network图汇总
# summary(network)
# 网状图导出
#绘制tiff图片并保存在Rwork文件夹下
# tiff(file="network.tiff")
# plot(network)
# dev.off()
# help(plot)
#使用JAGS、OpenBUGS及WinBUGS进行运算，NA为依次全部运算#JAGS
result <-mtc.run(model, sampler ="rjags", n.adapt = 5000, n.iter = 20000, thin = 1)
#Setting the sampler is deprecated, only JAGS is supported. #sunmin dell PC
#拟合模型，需要等一段时间，拟合完毕后，结果就ok了，先吃块西瓜。。。。。。
# Setting the sampler is deprecated, only JAGS is supported.不支持其他程序OpenBUGS
#使用JAGS、OpenBUGS及WinBUGS进行运算，NA为依次全部运算#OpenBUGS
# results <-mtc.run(model, sampler ="BRugs", n.adapt = 5000, n.iter = 20000, thin = 1)
#results <-mtc.run(model, sampler ="BRugs", n.adapt = 5000, n.iter = 20000, thin = 1)
# results <-mtc.run(model, sampler ="R2WinBUGS", n.adapt = 5000, n.iter = 20000, thin = 1)
#三个软件都运行
#results <-mtc.run(model, sampler =NA, n.adapt = 5000, n.iter = 20000, thin = 1)
#网状图在可这里可以直接出图
#模型网状图
plot(network)
# # points(network, cex = .5, col = "dark red") #改变颜色没有用
# #network图汇总，右下角的窗口中出现了，我们要的图形
# 网状图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0101network.tiff")
plot(network)
dev.off()
#文件夹中出现了，我们要的森林图，其他同理可得
summary(network)
# $Description
# [1] "MTC dataset: Example"
#
# $`Studies per treatment` #每个治疗措施A B C 各有多少个study
# 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臂和3臂各有多少个研究study
# 2-arm 3-arm
# 52 12
#
# $`Studies per treatment comparison`
# t1 t2 nr
# 1 A B 3
# 2 A D 5
#诊断收敛性
gelman.plot(result)
# 收敛图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0102gelman.plot.tiff")
# 此处应该有 14个 gelman plot（图 2） ，可以采用如下命令进行全部显示：
#3行5列
par(mfrow=c(3,5))
gelman.plot(result,auto.layout =F)
dev.off()
#结果汇总
summary(result)
#密度图
plot(result)
# 密度图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0103密度图.tiff")
# 此处应该有 14个 gelman plot（图 2） ，可以采用如下命令进行全部显示：
#3行5列
# par(mfrow=c(3,5))
plot(result)
dev.off()
#森林图
forest(result)
# 森林图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0104森林图.tiff")
# 此处应该有 14个 gelman plot（图 2） ，可以采用如下命令进行全部显示：
#3行5列
# par(mfrow=c(3,5))
forest(result)
dev.off()
#等级排名
ranks <- rank.probability(result)
print(ranks)
#排序图，堆积排序图
tiff(file="0105堆积排序图.tiff")
plot(ranks) # plot a cumulative rank plot
dev.off()
#单个排序图，等级图
tiff(file="0106单个排序图.tiff")
barplot(t(ranks), beside=TRUE) # plot a 'rankogram'
dev.off()
#单个排序图，等级图，我们让色彩丰富一些
tiff(file="0107单个排序图.tiff")
barplot(t(ranks), beside = TRUE,
col = c("lightblue", "mistyrose", "lightcyan",
"lavender", "cornsilk"),
legend = rownames(VADeaths), ylim = c(0, 1))
title(main = "Death Rates in Virginia", font.main = 8)
dev.off()
#单个排序图，等级图，我们让色彩更丰富一些
# help("colors")
# help(barplot)
barplot(t(ranks), beside=TRUE, col = colors(1:14))
# 单个排序图导出
#绘制tiff图片并保存在所在目录文件夹下
tiff(file="0108单个排序图.tiff")
# 此处应该有 14个 gelman plot（图 2） ，可以采用如下命令进行全部显示：
#3行5列
# par(mfrow=c(3,5))
barplot(t(ranks), beside=TRUE, col = colors(1:14))
dev.off()
# 相对影响森林图导出vsB
tiff(file="0109相对影响森林图导出vsB.tiff")
forest(relative.effect(result, "B"))
dev.off()
summary(relative.effect(result, "B", c("A", "C", "D")))
print(result$deviance)
## residual deviance plot
# #the first method，暂不能实现，但需要运行
if (model$data$ns.r2 + model$data$ns.rm == 0) {
tpl <- gemtc:::arm.index.matrix(model[['network']])
study <- matrix(rep(1:nrow(tpl), times=ncol(tpl)), nrow=nrow(tpl), ncol=ncol(tpl))
study <- t(study)[t(!is.na(tpl))]
devbar <- t(results$deviance$dev.ab)[t(!is.na(tpl))]
title <- "Per-arm residual deviance"
xlab <- "Arm"
} else {
nd <- model$data$na
nd[-(1:model$data$ns.a)] <- nd[-(1:model$data$ns.a)] - 1
devbar <- c(apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE), result$deviance$dev.re) / nd
study <- 1:length(devbar)
title <- "Per-study mean per-datapoint residual deviance"
xlab <- "Study"
}
plot(devbar, ylim=c(0,max(devbar, na.rm=TRUE)),
ylab="Residual deviance", xlab=xlab,
main=title, pch=c(1, 22)[(study%%2)+1])
# #the second method
###########################
for (i in 1:length(devbar)) {
lines(c(i, i), c(0, devbar[i]))
}
# w <- sign(r - rfit) * sqrt(devbar)
# plot(w, leverage, xlim=c(-3,3), ylim=c(0, 4.5))
# x <- seq(from=-3, to=3, by=0.05)
# for (c in 1:4) { lines(x, c - x^2) }
fit.ab <- apply(result$deviance$fit.ab, 1, sum, na.rm=TRUE)
dev.ab <- apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE)
lev.ab <- dev.ab - fit.ab
fit.re <- result$deviance$fit.re
dev.re <- result$deviance$dev.re
lev.re <- dev.re - fit.re
nd <- model$data$na
nd[-(1:model$data$ns.a)] <- nd[-(1:model$data$ns.a)] - 1
w <- sqrt(c(dev.ab, dev.re) / nd)
lev <- c(lev.ab, lev.re) / nd
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
# 重新运行下画图程序，可以出图
tiff(file="0110相对影响森林图导出vsB.tiff")
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
dev.off()
# #the third method第3种方法，以后备用哦
###########################
for (i in 1:length(devbar)) {
lines(c(i, i), c(0, devbar[i]))
}
# w <- sign(r - rfit) * sqrt(devbar)
# plot(w, leverage, xlim=c(-3,3), ylim=c(0, 4.5))
# x <- seq(from=-3, to=3, by=0.05)
# for (c in 1:4) { lines(x, c - x^2) }
fit.ab <- apply(result$deviance$fit.ab, 0.5, sum, na.rm=TRUE)
dev.ab <- apply(result$deviance$dev.ab, 0.5, sum, na.rm=TRUE)
lev.ab <- dev.ab - fit.ab
fit.re <- result$deviance$fit.re
dev.re <- result$deviance$dev.re
lev.re <- dev.re - fit.re
nd <- model$data$na
nd[-(1:model$data$ns.a)] <- nd[-(1:model$data$ns.a)] - 1
w <- sqrt(c(dev.ab, dev.re) / nd)
lev <- c(lev.ab, lev.re) / nd
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
# 目前实现网状 Meta 分析的软件中，进行数据处理模型多数采取贝叶斯模型。
# gemtc 程序包就是基于贝叶斯模型来实现数据处理的，其省略了建立相关的贝叶斯 model 的过程，甚至操作者不会使用调用的软件都可以得到相关结果。
# 下一讲我们深入讲解间接比较和异质性。
猴哥：Freescience公众号meta分析栏目现任主编。副主任医师，武汉大学肿瘤学博士，专注于胃肠道肿瘤分子生物学机制、系统评价/Meta分析、数据挖掘、临床统计研究。
想动手试一下上述流程？分享本文到百人群截图发送后台，获取R软件gemtc包；分享本文到朋友圈截图发送后台，获取R软件、Rstudio（Mac+windows版）和Stata。另有猴哥小礼包，获得方式戳这里
科研路，不孤单！
Freescience医学科研联盟全国火热招募ing
50家高校及医院的小伙伴已经加入啦
具体点这里
FS科研软件库，集合50+医学科研必备神器，现在统统打包分享，点这里
还有freescience交流群
加小秘书后拉进去
点开并长按二维码加小秘书为好友
未经许可 不得转载长按二维码关注