> library(BRugs)

> library(R2OpenBUGS)

> library(foreign)

> REmodeldata<-read.dta("G:/normal document/statistics software and video/R/network meta/network/REmodel_data.dta")

> 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 ~ dnorm(0,.0001) # vague priors for all trial baselines

+ for (k in 1:na) { # LOOP THROUGH ARMS

+ r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood

+ logit(p[i,k]) < - mu + 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 < - sum(dev[i,1:na])

+ for (k in 2:na) { # 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

>

> filename <- file.path("D://","REmodel.bug")

>

> ## write model file:

> write.model(REmodel, filename)

> ## and let’s take a look:

> file.show(filename)

> #load data from dataset

>

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

>

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

+ }

>

> 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(x86)/OpenBUGS/OpenBUGS323")

Error in bugs(data, inits, model.file = "D:/REmodel.bug", parameters = c("or", :

unused argument (bugs.directory = "D:/Program Files(x86)/OpenBUGS/OpenBUGS323")