# Empirical power considering Welton et al approach library(R2WinBUGS) path.working="C:/Users/Elsa/Desktop/Welton" # file where the results are saved path.winbugs="C:/WinBUGS14" # the code requires the statistical software Winbugs M<-45 # M the number of meta-analysis nPlus<-3 # nPlus the mean number of trials with the characteristic of interest nMinus<-5 # nMinus the mean number of trials without the characteristic of interest beta<-log(0.8) # overall impact term phi<-0.04 # phi^2 the mean between meta-analysis variance kappa<-0.16 # kappa^2 the mean between-trial variance sigma<-0.5 # sigma^2 mean within-study variance of the treatment effect estimates # define Welton's bayesian model and parameters welton_model <- function(){ for (i in 1:K){ p[i]<-1/(sigma[i]*sigma[i]) y[i] ~ dnorm(theta[i], p[i] ) theta_0[i]<-d[meta[i]]+ beta[i] * x[i] theta[i]~dnorm(theta_0[i], p.d[meta[i]] ) beta[i] ~ dnorm( b[meta[i]] , p.kappa ) } for (m in 1:J){ b[m] ~ dnorm(b0,p.phi) } # priors for (m in 1:J){ d[m] ~ dnorm(0,0.001) s.d[m] ~ dunif(0,2) p.d[m]<-1/(s.d[m]*s.d[m]) } b0 ~ dnorm(0,0.001) s.phi ~ dunif(0,2) p.phi<-1/(s.phi*s.phi) s.kappa ~ dunif(0,2) p.kappa<-1/(s.kappa*s.kappa) } parameters <- c("b0", "s.phi", "s.kappa") # Due to issues inherent to the Winbugs software, # scenarios are run by groups of 10, saved on the fly # and results are merged when computation is complete. # function teniterations() makes ten iterations teniterations<-function(iteration){ results_10<-data.frame(matrix(NA,ncol=2,nrow=0)) names(results_10)<-c("power","Rhat") nsim<-0 while (nsim <10){ nsim<-nsim+1 index<-0 K<-0 # generate meta-epidemiological datasets base <- data.frame(matrix(nrow=0, ncol=5)) names(base)<-c("ma", "trial", "x", "y.nm", "sigma.nm") for (m in 1:M){ check<-0 # condition to have more than 2 trials # at least 1 with the characteristic # at least 1 without the characteristic while (check==0){ # number of trials with the characteristic, median=3 n.plus<-round(rlnorm(1, meanlog = log(nPlus)-1.056/2, sdlog = sqrt(1.056)),0) # number of trials without the characteristic, median=5 n.minus<-round(rlnorm(1, meanlog = log(nMinus)-1.056/2, sdlog = sqrt(1.056)),0) n.m<-n.plus+n.minus # number of trials in meta-analysis m if (n.m>2 & n.plus!=0 & n.minus!=0) check=1 } # between-trial variance of the effect estimate tau2<-rlnorm(1,-2.56,sqrt(1.74)) phi2<-phi*rnorm(1,0,1) for (n in 1:n.m){ index<-index+1 temp1<-m temp2<-n temp3<-1*(temp2>n.minus) # n.minus trials without characteristics temp5<-sqrt(runif(1,min=0.25*sigma^2,max=1.75*sigma^2)) temp4<-(beta+phi2+kappa*rnorm(1,0,1) )*temp3+sqrt(tau2+temp5^2)*rnorm(1,0,1) base[index,]<-c(temp1,temp2,temp3,temp4,temp5) } K<-K+n.m } # define data and inits for Welton model data<-list(J=M,K=K, meta=base$ma,x=base$x,y=base$y.nm,sigma=base$sigma.nm) inits <- function(){ list(d = c(rnorm(J, 0, 1000)),s.d = c(runif(J, 0, 2)),b0 = c(rnorm(1, 0, 1000)), s.phi=1, s.kappa=1) } # welton model welton.sim <- bugs(data, inits,parameters, welton_model, n.chains=3, n.iter=30000, n.burnin = 15000, n.thin = 1, bugs.directory=path.winbugs, working.directory=path.working, DIC=FALSE) # results: significance on beta estimate and Rhat for convergence results_10[nsim,]<-c(ifelse(welton.sim$summary[1,3]<0 & 0