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.x[1]*x1.IPD[i] + beta.x[2]*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]) } invR.IPD[1:5, 1:5] <- inverse(R.IPD[ , ]) for(j in 1:5){ for(k in 1:5){ R.IPD[j, k]<-tau.IPD*tau.IPD*pow(rho.IPD,step(abs(j-k)-0.5)) } } tau.IPD ~ dunif(0.01,10) rho.IPD ~ dunif(0,1) ## Fit AD const <- 10000 # this constant has to be large enough to ensure all phi > 0 for (i in 1:NAD) { meanAD[s[i],t[i]] <- mu.AD[t[i]] + v[s[i],t[i]] + beta.x[1]*x1.AD[i] + beta.x[2]*x2.AD[i] # log likelihood and use zeros trick sdsd[i] <- sdAD[i]/sqrt(n[i]) phi[i] <- a0[i]*( log(sdsd[i]) + 0.5*pow((yAD[i]- meanAD[s[i],t[i]])/sdsd[i] , 2) ) + const zeros[i] ~ dpois(phi[i]) } for (j in 1:NSAD) { v[j, 1:27] ~ dmnorm(zero.AB[1:27], invR.AD[1:27, 1:27]) } invR.AD[1:27, 1:27] <- inverse(R.AD[ , ]) for(j in 1:27){ for(k in 1:27){ R.AD[j, k]<-tau.AD[1]*tau.AD*pow(rho.AD,step(abs(j-k)-0.5)) } } tau.AD ~ dunif(0.01,10) rho.AD ~ dunif(0,1) ## Redefine mu parameters including interactions for(i in 1:27) { mu.AD[i] <- mu[i] } for(i in 1:4) { mu.IPD[i] <- mu[i] } mu.IPD[5] <- mu[28] ## Prior to beta.x for (i in 1:2) { beta.x[i] ~ dnorm(0,0.01) } for (k in 1:NT) { mu[k] ~ dnorm(0, 0.001) d[k] <- mu[k]-mu[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] + beta.x[1]*(57-57) + beta.x[2]*(temp.x2[j]-9) } } }