1. 安裝gemtc包
使用的是Rstudio收恢，如果Rstudio安裝不成功俏脊，可以先用Rgui安裝：
install.packages('code')
2. 載入程序包
library(gemtc)
library(code)
library(codetools)
library(testthat)
3. 導(dǎo)入數(shù)據(jù)
這里是一個(gè)gemtc包中的例子柠辞，命名為file:
file <- system.file('extdata/luades-smoking.gemtc', package='gemtc')? ? #導(dǎo)入數(shù)據(jù)鲸阔，并命名為file
net = read.mtc.network(file)? ? #讀取gemtc文件订晌，并保存為net
4. 保存數(shù)據(jù)
write.csv(network$treatments, file = "D:\\Desktop\\network0001_treatments.csv")
write.csv(network$data.ab, file = "D:\\Desktop\\network0001_data.ab.csv")
5.讀取數(shù)據(jù)
treatments<- read.csv("D:\\Desktop\\network0001_treatments.csv",header=T,sep = ",",skip=0,row.names = 1)
#header=T,表示文件存在菲嘴，表示需要讀取數(shù)據(jù)的數(shù)據(jù)頭部饿自，skip=0,表示數(shù)據(jù)中沒(méi)有需要跳過(guò)的行，row.names = 1是第一列作為行名
#同樣龄坪，把另外一個(gè)表也讀取進(jìn)來(lái)
data <- read.csv("D:\\Desktop\\network0001_data.ab.csv",header=T,sep = ",",row.names = 1)
6. 創(chuàng)建network對(duì)象昭雌，建立成network 單臂長(zhǎng)數(shù)據(jù)格式
network<- mtc.network(data, description="Luades_smoking", treatments=treatments)
7. 構(gòu)建模型，一致性模型
model <- mtc.model(network, type = 'consistency')
#其中健田， network 為 network 數(shù)據(jù)烛卧， type 為是否選取一致性模型
#這里提示沒(méi)有安裝slam包
8.運(yùn)算結(jié)果
result <- mtc.run(model)? ? #簡(jiǎn)單命令；需要電腦上安裝jags抄课，安裝地址下載jags;
result <-mtc.run(model, sampler ="rjags", n.adapt = 5000, n.iter = 20000, thin = 1)? ?
#復(fù)雜命令唱星，通過(guò)sampler命令選取軟件的調(diào)用方式；“n.adapt”為預(yù)迭代次數(shù)跟磨，“n.iter”為迭代運(yùn)算次數(shù)间聊，“thin”為步長(zhǎng)。
9. 網(wǎng)狀證據(jù)圖等
plot（network）? #模型網(wǎng)狀圖
#繪制tiff圖片并保存
tiff(file="network.tiff")
plot(network)
dev.off()
#診斷收斂性
gelman.plot(result)
#收斂圖導(dǎo)出
tiff(file="gelman.plot.tiff")
gelman.plot(result,auto.layout =F)
dev.off()
plot(result)? ? #密度圖
forest(result)? ? #森林圖
ranks <- rank.probability(result)? ? #等級(jí)排名
print(ranks)
#堆積排序圖
tiff(file="堆積排序圖.tiff")
plot(ranks)
dev.off()
#單個(gè)排序圖
tiff(file="單個(gè)排序圖.tiff")
barplot(t(ranks), beside=TRUE)
dev.off()
# 相對(duì)影響森林圖導(dǎo)出vsB
tiff(file="0109相對(duì)影響森林圖導(dǎo)出vsB.tiff")
forest(relative.effect(result, "B"))
dev.off()
10. 殘差
print(result$deviance)
11.異質(zhì)性圖
result.anohe <- mtc.anohe(network, n.adapt=1000, n.iter=5000)
summary.anohe <- summary(result.anohe)
plot(summary.anohe, xlim=log(c(0.2, 5)))
summary.anohe
summary.anohe$consEffects
summary.anohe$studyEffects
summary.anohe$pairEffects
12.節(jié)點(diǎn)劈裂法,探討模型的一致性和不一致性
network
result <-mtc.nodesplit(network)
summary(result)
names(result)
# [1] "d.A.C" "d.A.D" "d.B.D" "d.C.D" "consistency"
summary(result$d.A.C)
# Overall summary and plot
summary.ns <- summary(result)
print(summary.ns)
plot(summary.ns)
