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

其明技术专家 meta分析, 博客, 基础科研 2017年5月3日

2.21W 0 0

我们仅仅是代码的编辑者、整合者、搬运工，仅免费传授方法，下文数据和代码取自于网络和免费软件“R语言说明书”，如果您觉得我们侵犯了您的版权，请通知我们撤稿。请大家谅解，谢谢！

#第一部分感谢兰州大学循证医学中心的代码和说明 #

#1 代码来源和文献解读####

下载链接： doi: 10.1111/jcpt.12410

文献解读（摘要）：FOLFOX方案化疗是结直肠癌标准方案。本文确定了63符合条件的研究（4837例患者），涉及9 CHIS。配对META分析表明，与单独的FOLFOX化疗方案相比，FOLFOX联合艾迪注射液和复方苦参注射液能显著提高整体的总反应率和生活质量（肿瘤学疗效评价指标），降低的恶心和呕吐（Ⅲ-Ⅳ），腹泻的发病率（Ⅲ-Ⅳ），血小板减少（Ⅲ-Ⅳ），白细胞减少（III-IV）和外周神经毒性（III-IV）。根据间接结果相比较而言，无统计学差异。综合考虑了小样本量效应，爱迪+ FOLFOX，参芪扶正+ FOLFOX和复方苦参+ FOLFOX临床最佳疗效和安全性的概率最高。

#第二部分网状meta分析的实现#

联合stata和R2winbugs软件实现

1.1 软件安装及程序包加载

STATA 13.1的安装和视频请参考猴哥以前的推文。

R 软件版本为R-3.1.2，获取网址为http://www.r-project.org/；采用OpenBUGS软件执行数据分析，因此还需从http://www.openbugs.net/w/Downloads下载该软件。所有软件下载完成后，双击安装包进行安装。

本文附件代码如下：

#Binary outcome - random-effects model 二分类数据随机效应模型

# Binomial likelihood, logit link二项式分布，连接函数：logit（变换）

# Random effects model随机效应模型

#vague priors本文采用模糊先验，当然，还可以用精确先验

#适合2arm和3arm、多arm

#计算lor or SUCRA cumeffectiveness effectiveness等

#本次代码为完整代码

model{ # *** PROGRAM STARTS

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 distributions (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 treatment

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

# pair wise ORs and LORs for all possible 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 of treatments#

for(k in 1:nt) {

order[k]<- nt+1-rank(d[],k)

#因本次为有利指标（有效events），故代码为order[k]<- nt+1-rank(d[],k)

#若为有害指标（副作用events），代码为order[k]<- rank(d[],k)

#order[k]<- rank(d[],k)

most.effective[k]<-equals(order[k],1)

for(j in 1:nt) {

effectiveness[k,j]<- equals(order[k],j)

}

for(k in 1:nt) {

for(j in 1:nt) {

cumeffectiveness[k,j]<- sum(effectiveness[k,1:j])

}

#SUCRAS#

for(k in 1:nt) {

SUCRA[k]<- sum(cumeffectiveness[k,1:(nt-1)]) /(nt-1)

}}

# *** PROGRAM ENDS

## （猴哥）固定效应模型代码如下，没有sd和taud的赋值是固定效应模型与随机效应模型的区别，如上文标红加粗部分##

# 下文红色加粗部分为精确先验的代码，如果没有精确先验的值，请忽略#

# Binomial likelihood, logit link

# Fixed effects model

model{ # *** PROGRAM STARTS

for(i in 1:ns){ # LOOP THROUGH STUDIES

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

# model for linear predictor

logit(p[i,k]) <- mu[i] + d[t[i,k]] - d[t[i,1]]

# expected value of the numerators

rhat[i,k] <- p[i,k] * n[i,k]

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

}

totresdev <- sum(resdev[]) # Total Residual Deviance

d[1]<-0 # treatment effect is zero for reference treatment

# vague priors for treatment effects

for (k in 2:nt){ d[k] ~ dnorm(0,.0001) }

# Provide estimates of treatment effects T[k] on the natural (probability) scale

# Given a Mean Effect, meanA, for 'standard' treatment A,

# with precision (1/variance) precA

A ~ dnorm(meanA,precA)

for (k in 1:nt) { logit(T[k]) <- A + d[k] }

} # *** PROGRAM ENDS

###暂讲到这里，谢谢大家###

################### 猴哥与三军医大博士对话###################

三军医大博士：@猴哥能讲讲R2winbug的初始值设置问题吗？
猴哥：excel可以得到初始值

你用r随机生成 -5 到5 的数据或者 excel 随机生成 -7到7
三军医大博士：
#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))
}
这个实例中，初始值定义中d代表什么，d=c( NA, 0, 0, 0, 0), sd=1中括号中的5个数字代表什么，NA代表什么，mu中有的是25个0，有的是25个-3，这25个数从何而来，是因为纳入了25个实验么？我怎么根据自己的数据设定初始值呢？