# Binomial likelihood, logit link

# Random effects model for multi-arm trials

#

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  <- sqrt(1/tau)

tau ~ dgamma(0.001,0.001)   



# ranking 

for (k in 1:nt) {

   rk[k] <- rank(d[],k)                           # assumes events are “bad”

   for (m in 1:nt){

   	best[k,m] <- equals(rk[k],m)                   

	}

}



#All pairwise log odds ratios and odds ratios

for (c in 1:(nt-1)) { 

	for (k in (c+1):nt) { 

		lor[c,k] <- d[k] - d[c]

		log(or[c,k]) <- lor[c,k]

		# predictive distribution

		lornew[c,k] ~ dnorm(lor[c,k], tau)

		log(ornew[c,k]) <- lornew[c,k]

		}

	 }



# Non Inf given an OR non-inf margin

for (c in 1:(nt-1)) { 

	for (k in (c+1):nt) { 

		noninftemp[c,k] <- lormargin - lor[c,k] 

		noninf[c,k] <- step(noninftemp[c,k])

		}

	 }





}                                    # *** PROGRAM ENDS