model {

  ## Fit IPD

  for (i in 1:NIPD) {
    yIPD[i] ~ dnorm(meanIPD[study[i],person[i]], precIPD[study[i],drug[i]])
    meanIPD[study[i],person[i]] <- mu.IPD[drug[i]] + v.IPD[study[i],drug[i]] + beta.x1.IPD[study[i],drug[i]]*x1.IPD[i] + beta.x2.IPD[study[i],drug[i]]*x2.IPD[i]   
  }

  for (i in 1:NSIPD) {
    for(k in 1:5) {
      precIPD[i,k] <- 1/pow(sdIPD[i,k],2) 
      sdIPD[i,k] ~ dunif(0.01, 10)
    }
  }

  for (j in 1:NSIPD) { 
    v.IPD[j, 1:5] ~ dmnorm(zero.AB.IPD[1:5], invR.IPD[1:5, 1:5]) 
    beta.x1.IPD[j, 1:5] ~ dmnorm(b.x1.IPD[1:5], invR.x1.IPD[1:5, 1:5])
    beta.x2.IPD[j, 1:5] ~ dmnorm(b.x2.IPD[1:5], invR.x2.IPD[1:5, 1:5])
  }

  invR.IPD[1:5, 1:5] <- inverse(R.IPD[ , ])
  invR.x1.IPD[1:5, 1:5] <- inverse(R.x1.IPD[ , ])
  invR.x2.IPD[1:5, 1:5] <- inverse(R.x2.IPD[ , ])

  for(j in 1:5){    
    for(k in 1:5){ 
      R.IPD[j, k]<-tau.IPD[1]*tau.IPD[1]*pow(rho.IPD[1],step(abs(j-k)-0.5))
      R.x1.IPD[j, k]<-tau.IPD[2]*tau.IPD[2]*pow(rho.IPD[2],step(abs(j-k)-0.5))
      R.x2.IPD[j, k]<-tau.IPD[3]*tau.IPD[3]*pow(rho.IPD[3],step(abs(j-k)-0.5))
    }
  }

  for (j in 1:3) {
    tau.IPD[j] ~ dunif(0.01,10)
    rho.IPD[j] ~ dunif(0,1)
  }

  ## Fit AD

  for (i in 1:NAD) {
    yAD[i] ~ dnorm(meanAD[s[i],t[i]], precAD[i])
    precAD[i] <- n[i]/pow(sdAD[i],2)
    meanAD[s[i],t[i]] <- mu.AD[t[i]] + v[s[i],t[i]] + beta.x1.AD[s[i],t[i]]*x1.AD[i] + beta.x2.AD[s[i],t[i]]*x2.AD[i]
  }

  for (j in 1:NSAD) { 
    v[j, 1:27] ~ dmnorm(zero.AB[1:27], invR.AD[1:27, 1:27]) 
    beta.x1.AD[j, 1:27] ~ dmnorm(b.x1.AD[1:27], invR.x1.AD[1:27, 1:27])
    beta.x2.AD[j, 1:27] ~ dmnorm(b.x2.AD[1:27], invR.x2.AD[1:27, 1:27])
  }

  invR.AD[1:27, 1:27] <- inverse(R.AD[ , ])
  invR.x1.AD[1:27, 1:27] <- inverse(R.x1.AD[ , ])
  invR.x2.AD[1:27, 1:27] <- inverse(R.x2.AD[ , ])

  for(j in 1:27){    
    for(k in 1:27){ 
      R.AD[j, k]<-tau.AD[1]*tau.AD[1]*pow(rho.AD[1],step(abs(j-k)-0.5))
      R.x1.AD[j, k]<-tau.AD[2]*tau.AD[2]*pow(rho.AD[2],step(abs(j-k)-0.5))
      R.x2.AD[j, k]<-tau.AD[3]*tau.AD[3]*pow(rho.AD[3],step(abs(j-k)-0.5))
    }
  }

  for (j in 1:3) {
    tau.AD[j] ~ dunif(0.01,10)
    rho.AD[j] ~ dunif(0,1)
  }

  ## Redefine mu parameters including interactions
  for(i in 1:27) { 
    mu.AD[i] <- mu[i] 
    b.x1.AD[i] <- b.x1[i]
    b.x2.AD[i] <- b.x2[i]    
  }
  for(i in 1:4) { 
    mu.IPD[i] <- mu[i] 
    b.x1.IPD[i] <- b.x1[i]
    b.x2.IPD[i] <- b.x2[i]
  }
  mu.IPD[5] <- mu[28]
  b.x1.IPD[5] <- b.x1[28]
  b.x2.IPD[5] <- b.x2[28]

  ## Prior to b.x1[] b.x2[]
  for (k in 1:NT) { 
    b.x1[k] ~ dnorm(bb.x[1],precx[1])
    b.x2[k] ~ dnorm(bb.x[2],precx[2])
  }

  for (i in 1:2) {
    bb.x[i] ~ dnorm(0, 0.001)
    precx[i] <- 1/pow(taux[i],2)
    taux[i] ~ dunif(0.01,10)
  }

  for (k in 1:NT) { 
    mu[k] ~ dnorm(0, 0.001) 
    d[k] <- mu[k]-mu[1]
    d.x1[k] <- b.x1[k] - b.x1[1]
    d.x2[k] <- b.x2[k] - b.x2[1]
  }  

  ## rank
  for (k in 1:NT) { T[k] <- mu[k] }
  T.rank <- rank(T)
  for (k in 1:NT) {
    for(j in 1:NT) {
      best[k,j] <- equals(T.rank[k],j)
    }
  }
  for (k in 1:NT) { best1[k] <- equals(T.rank[k],1) }
 
  # Prediction
  # at age 57

  for (k in 1:NT) {
    for (j in 1:8) {
      pred[k,j] <- mu[k] + b.x1[k]*(57-57) + b.x2[k]*(temp.x2[j]-9)
    }
  }

}