扫码联系客服
转载丁香园
本来想用WinBUGS进行介绍的，但张版主的《实用循证医学方法学》已经介绍的很清楚了，所以小弟也就不班门弄斧了。关于贝叶斯定理在张版主的书中已经有介绍，有兴趣的同志可以详细阅读一下，小弟就不在这里重复了。本研究中的数据采用张版主的的《实用循证医学方法学》第十九章 贝叶斯Meta分析 实例4中的数据，计算出来的数据可以与原数据进行比较。数据随后附上。
闲话不多说了，下面就开始吧！
1.下载R2WinBUGS包与加载R2WinBUGS包
| 2.模型构建及写入 2.1先在R软件中新建程序脚本，以便于程序写入，如下图
|
2.2 模型构建
本文采用的模型为用于二分类数据分析的完全随机效应模型，程序及模型如下：
library(R2WinBUGS)
REmodel<-function(){
# Binomial likelihood, logit link
# Random effects model for multi-arm trials
for(i in 1:ns){ # LOOP THROUGH STUDIES
w[i,1] <- 0 # adjustment for multi-arm trials is zero for control arm
delta[i,1] < - 0 # treatment effect is zero for control arm
mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines
for (k in 1:na[i]) { # LOOP THROUGH ARMS
r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood
logit(p[i,k]) < - mu[i] + delta[i,k] # model for linear predictor
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators
#Deviance contribution
dev[i,k] < - 2 * (r[i,k] * (log(r[i,k]) -log(rhat[i,k]))
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k] -rhat[i,k]))) }
# summed residual deviance contribution for this trial
resdev[i] < - sum(dev[i,1:na[i]])
for (k in 2:na[i]) { # LOOP THROUGH ARMS
# trial-specific LOR distributions
delta[i,k] ~ dnorm(md[i,k],taud[i,k])
# mean of LOR di stributions (with multi-arm trial correction)
md[i,k] <- d[t[i,k]] - d[t[i,1]] + sw[i,k]
# precision of LOR distributions (with multi-arm trial correction)
taud[i,k] < - tau *2*(k -1)/k
# adjustment for multi-arm RCTs
w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]])
# cumulative adjustment for multi-arm trials
sw[i,k] < - sum(w[i,1:k -1])/(k-1)
}
}
totresdev < - sum(resdev[]) # Total Residual Deviance
d[1]<-0 # treatment effect is zero for reference t reatment
# vague priors for treatment effects
for (k in 2:nt){ d[k] ~ dnorm(0,.0001) }
sd ~ dunif(0,5) # vague prior for between -trial SD
tau <- pow(sd,-2) # between-trial precision = (1/between -trial variance)
# pairwise ORs and LORs for all possi ble pair-wise comparisons, if nt>2
for (c in 1:(nt -1)) {
for (k in (c+1):nt) {
or[c,k] < - exp(d[k] - d[c])
lor[c,k] < - (d[k] -d[c])
}
}
# ranking on relative scale
for (k in 1:nt) {
rk[k] < - nt+1 -rank(d[],k) # assumes events are “good”
# rk[k] < - rank(d[],k) # assumes events are “bad”
best[k] < - equals(rk[k],1) #calculate probability that treat k is best
}
# *** PROGRAM ENDS
}# End of model
第一句library(R2WinBUGS)为加载R2WinBUGS包
第二句开始到最后一句为随机效应模型的构建，REmodel为模型名称，这可以随自己需要更换，funtion()内为整个用于二分类数据分析的完全随机效用模型。当然，模型的选择需根据数据来定，这里的模型不代表适用于一切二分类数据。
3.存储模型及查看是否存储成功
filename <- file.path("D://","REmodel.bug")
## write model file:
write.model(REmodel, filename)
## and let’s take a look:
file.show(filename)
若存储成功R软件会自动弹出信息窗口，如下图
4.导入所需数据
请将文中所需数据导入R软件中，导入方法请参阅小弟帖子http://www.dxy.cn/bbs/topic/20225778
导入数据后如下图
5.在程度中加载数据框中的数据
t1<-REmodeldata$t1
r1<-REmodeldata$r1
n1<-REmodeldata$n1
t2<-REmodeldata$t2
r2<-REmodeldata$r2
n2<-REmodeldata$n2
t3<-REmodeldata$t3
r3<-REmodeldata$r3
n3<-REmodeldata$n3
na<-REmodeldata$na
t<-c(t1,t2,t3)
r<-c(r1,r2,r3)
n<-c(n1,n2,n3)
dim(t)<-c(25,3)
dim(r)<-c(25,3)
dim(n)<-c(25,3)
ns<-25
nt<-5
data<-list("t","r","n","na","ns","nt")
上述中REmodeldata为数据框名字，t为治疗方法，r为每项研究组发生事件数，n每项研究总数，na为每项研究中的比较组数，ns为本分析中总的研究数，nt为本分析中所涉及的治疗措施数
6.构建初始值
#Set Initial Values
inits<-function(){
#chain 1
list(d=c( NA, 0, 0, 0, 0), sd=1,
mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0))
#chain 2
list(d=c( NA, -1, -1, -1, -1) , sd=4,
mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3,
-3, -3, -3, -3, -3, -3, -3, -3, -3, -3,
-3, -3, -3, -3, -3))
#chain 3
list(d=c( NA, 2, 2, 2, 2), sd=2,
mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0,
-3, -3, 0, 3, 5, -3, -3, -1, -3, -7,
-3, -3, 5, -1, 7))
}
7.进行网络meta分析
REmodel.sim<-bugs(data,inits,model.file="D:/REmodel.bug",
parameters=c("or","totresdev"),
n.chains=3,n.iter=40000,n.burnin=10000,
bugs.directory="D:/Program Files/WinBUGS14")
REmodel.sim为数据分析综合后的变量名，可以根据自己需要改变。data，inits上面已有介绍。model.file="D:/REmodel.bug"为我们刚构建的模型，parameters=c("or","totresdev")为本文所需要监测的指标，包括or和totresdev，n.chains=3表示使用3条MCMC链，n.iter=40000表示迭代40000次，n.burnin=10000表示10000次退火，bugs.directory="D:/Program Files/WinBUGS14"表示WinBUGS的绝对安装路径。
8.显示分析结果
输入命令print(REmodel.sim, digits.summary=4)，得出结果如下图
digits.summary=4表示显示出的数据保留4位小数
与张版主的所得出的数据对比，发现无明显差异。
另介绍一下totresdev这个参数，此参数为评价模型契合度（model fit）的参数，当此参数与原数据的总数相等时，模型契合度最好，现totresdev为48.0463，原数据的总数为53，基本相符，但仍有一定的差异。
| 9.评价模型的收敛度（convergence） 输入命令plot(Remodel.sim)，得出图如下
根据参考文献介绍，当所有的Rhat()参数都接近1.0时，模型收敛度较好 |
欢迎关注生信人
