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)
  }

  # spike and slab now
  for(i in 1:4) { 
    mu.IPD[i] ~ dnorm(mu.AD[i], sstau[i]) 
    b.x1.IPD[i] ~ dnorm(b.x1.AD[i], sstau.x1[i])
    b.x2.IPD[i] ~ dnorm(b.x2.AD[i], sstau.x2[i])
  }

  # sstau <- 1000    # *forces* commensurability of AD and IPD
  # sstau <- 0.001   # essentially disconnects AD and IPD
  # sstau ~ dgamma(.01, .01)  # Standard WinBUGS vague hyperprior

  for (k in 1:4) {
    tee[k,1] ~ dnorm(20,1)     # R_tau is (essentially) 20
    tee[k,2] ~ dgamma(0.1, 0.1)I(0.1, 5)  # replacement for the true slab
    flip[k] ~ dbern(0.5)     # p_tau is 0.5
    pick[k] <- flip[k] + 1
    sstau[k] <- tee[k,pick[k]]
    sstau_rev[k] <- sqrt(1/sstau[k])

    tee.x1[k,1] ~ dnorm(20,1)     # R_tau is (essentially) 20
    tee.x1[k,2] ~ dgamma(0.1, 0.1)I(0.1, 5)  # replacement for the true slab
    flip.x1[k] ~ dbern(0.5)     # p_tau is 0.5
    pick.x1[k] <- flip.x1[k] + 1
    sstau.x1[k] <- tee[k,pick.x1[k]]
    sstau_rev.x1[k] <- sqrt(1/sstau.x1[k])

    tee.x2[k,1] ~ dnorm(20,1)     # R_tau is (essentially) 20
    tee.x2[k,2] ~ dgamma(0.1, 0.1)I(0.1, 5)  # replacement for the true slab
    flip.x2[k] ~ dbern(0.5)     # p_tau is 0.5
    pick.x2[k] <- flip.x2[k] + 1
    sstau.x2[k] <- tee[k,pick.x2[k]]
    sstau_rev.x2[k] <- sqrt(1/sstau.x2[k])
  }

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

  for (k in 1:NT) { mu.temp[k] ~ dnorm(0, 0.001) }
  ## Prior to b.x1.AD[] b.x2.AD[]
  for (k in 1:NT) { 
    b.x1.temp[k] ~ dnorm(bb.x[1],precx[1])
    b.x2.temp[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:4) { 
    mu[k] <- mu.IPD[k] 
    b.x1[k] <- b.x1.IPD[k]
    b.x2[k] <- b.x2.IPD[k]
  }
  for(k in 5:28) { 
    mu[k] <- mu.temp[k] 
    b.x1[k] <- b.x1.temp[k]
    b.x2[k] <- b.x2.temp[k]
  }

  for (k in 1:NT) { 
    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)
    }
  }

}