Model.6.2.txt

model {
    for (i in 1:30) {
     # Likelihood:
     r[i] ~ dbin(p[i], n[i])
      
      # Logistic regression model:
       logit(p[i]) <- mu[s[i]] + delta[s[i]] * equals(t[i], 2)
    }
  for (j in 1:15) {
   # Hierarchical random 
   # effects model:
    delta[j] ~ dnorm(d, prec) 
   
 # Prior for baseline
 # random effect, µ:
  mu[j] ~ dnorm(0, 1.0e-6)
  }
 
 # Prior for Pop. 
 # treatment effect:
  d ~ dnorm(0, 1.0e-6)
 # Population OR:
  OR <- exp(d)
 # Prior for btw. studies
 # sd:
  tau ~ dnorm(0, 1.0e-6)T(0, )
 # Variance:
  tau.sq <- (tau * tau)
 # Precision:
  prec <- 1 / (tau.sq)
 
 # Replicate LOR for 
 # prediction:
  delta.new ~ dnorm(d, prec)
 # Mean random
 # effect:
  delta[19] <- exp(d)
 
 # Estimated no. deaths
 # in control group:
  a <- r[31]
 # Estimated no. survivors
 # in control group:
  b <- n[31] - r[31]
 # Draw on new control
 # probability:
  p[31] ~ dbeta(a, b)
 
 # Form new treatment
 # probability:
  logit(p[32]) <- logit(p[31]) + delta.new
  
 # Draw new no. of deaths
 # in treatment group:
  r.new ~ dbin(p[32], n[32])
 # Record whether predicted no.
 # of deaths exceeds observed:
 p.crossval <- step(r.new - r[32]) - .5 * equals(r.new, r[32])

}