网络meta分析教程
转载丁香园
本来想用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$t1r1<-REmodeldata$r1n1<-REmodeldata$n1t2<-REmodeldata$t2r2<-REmodeldata$r2n2<-REmodeldata$n2t3<-REmodeldata$t3r3<-REmodeldata$r3n3<-REmodeldata$n3na<-REmodeldata$nat<-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<-25nt<-5data<-list("t","r","n","na","ns","nt")上述中REmodeldata为数据框名字，t为治疗方法，r为每项研究组发生事件数，n每项研究总数，na为每项研究中的比较组数，ns为本分析中总的研究数，nt为本分析中所涉及的治疗措施数
6.构建初始值
#Set Initial Valuesinits<-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时，模型收敛度较好欢迎关注
网络meta分析教程
相关文章