> 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")
R2OpenBUGS的代码问题求大神