meta分析

网状meta分析必备技能2-R软件gemtc包

R 软件风靡全球。R 软件 gemtc 程序包通过调用对应的 rjags 程序包来执行网状 Meta 分析是最新的流行趋势。同时,还可以为 GeMTC 软件生成相关的数据存储文件。其本质是调用基于 MCMC 法的软件进行网状 Meta 分析。本文简要介绍 R 软件 gemtc 程序包调用rjags软件进行网状 Meta。

废话少说,代码如下:

#1清除环境变量,设置路径--------------------------------------------------

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")

#2载入程序包--------------------------------------------------

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)

#3导入数据-----------------------------------------------------

#3.1导入数据,第1中方法----------------------------------------

# 导入治疗措施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.2导入数据,第2中方法----------------------------------------

# 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)

#4 创建network对象,建立成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)

#6####也可以直接读取

##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)

#7构建模型####

#一致性模块

# 其中, 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,别害怕

#7网状证据图####

#网状图在可这里可以直接出图

#模型网状图

# plot(network)

# # points(network, cex = .5, col = "dark red")  #改变颜色没有用

# #network图汇总

# summary(network)

# 网状图导出

#绘制tiff图片并保存在Rwork文件夹下

# tiff(file="network.tiff")

# plot(network)

# dev.off()

# help(plot)

#8运行结果####

#使用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)

#10网状证据图####

#网状图在可这里可以直接出图

#模型网状图

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()

#11相对影响森林图

# 相对影响森林图导出vsB

tiff(file="0109相对影响森林图导出vsB.tiff")

forest(relative.effect(result, "B"))

dev.off()

summary(relative.effect(result, "B", c("A", "C", "D")))

#12残差

print(result$deviance)

#13杠杆效应####

## 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 的过程,甚至操作者不会使用调用的软

# 件都可以得到相关结果。

# 当前, gemtc 程序包不仅适合作二分类型数据

# 的网状 Meta 分析,还适合制作连续数据、生存分析数据的网状

# Meta 分析。但仅基于一致性模型,且仍需要进一步

# 优化,对于非一致性模型尚待研发。本例绘制的

# 网状关系图只能反映各研究组之间存在直接比较,

# 下一讲我们深入讲解间接比较和异质性。

 

拓展知识:

SNPmeta

SNPMETA1

SNPMETA2_检索式

SNPMETA3_索要原始文献+数据提取

 

教你如何做一篇meta分析1_

教你如何做一篇meta分析2_

教你如何做一篇meta分析3_哈温平衡的计算

教你如何做一篇meta分析4_NOS量表的制作

 

网状emta分析必备技能

网状emta分析必备技能1

网状meta分析必备技能2_R软件gemtc包

网状Meta分析必备技能3_数据初始值的设定

网状Meta分析必备技能4_Gemtc实现生存分析

网状meta分析必备技能5_HR风险比

网状meta分析必备技能6_R+Rstudio运用meta包做简单meta分析

网状meta分析必备技能7_使用R、GeMTC和STATA软件实现连续变量的网状Meta分析

 

网状meta分析stata简易教程

NMA(网状meta分析)stata简易教程(1)

NMA(网状meta分析)stata简易教程(2)网状图

NMA(网状meta分析)stata简易教程(3)贡献图

NMA(网状meta分析)stata简易教程(4)漏斗图的制作

NMA(网状meta分析)stata简易教程(5)排序图的制作

 

诊断性meta分析简单实现

诊断性meta分析简单实现(1)之revman

诊断性meta分析简单实现(2)

诊断性meta分析简单实现(3)

诊断性meta分析简单实现(4)_meta-disc

诊断性meta分析简单实现(5)_结果的展示和解读

 

其他

meta分析统计方法-随机对照试验的风险评估rob

生物标记物(基因)联合诊断模型的stata实现ROC(AUC)

生物标记物(基因)联合诊断模型的R实现ROC(AUC)

(0)

热评文章

发表评论