#From Nixon RM and Thompson SG, Methods for incorporating covariate adjustment, subgroup analysis and between-centre differences into cost-effectiveness evaluations, Health Economics, 2005.

#cost.re and eff.rs and factors by which the costs and effects are rescaled. 
#having a similar range of values for both makes for a better model simulation. 

model{
for(i in 1:N1){
    cost1[i]~dnorm(phi.c1[i],tau.c1)
    eff1[i]~dnorm(phi.e1[i],tau.e1)
    phi.c1[i]<-mu.c1+beta.c1*(eff1[i]-mu.e1)
    phi.e1[i]<-mu.e1
    da1[i]<-log(exp(-pow(cost1[i]-phi.c1[i],2)/(2*ss.c1))/(sqrt(2*3.1415927*ss.c1)))-log(cost.rs)
    db1[i]<-log(exp(-pow(eff1[i]-phi.e1[i],2)/(2*ss.e1))/(sqrt(2*3.1415927*ss.e1)))-log(eff.rs)
}

# node transformations
tau.c1<-1/ss.c1
ss.c1<-s.c1*s.c1
s.c1<-exp(ls.c1)
ss.c1.rs<-ss.c1*cost.rs*cost.rs
tau.e1<-1/ss.e1
ss.e1<-s.e1*s.e1
s.e1<-exp(ls.e1)
ss.e1.rs<-ss.e1*eff.rs*eff.rs

# prior distributions
# the limits of the prior distributions are put in a separate file cnen-lim.txt
# these may be need to be changed for different data sets.
ls.c1~dunif(ls.c1.lim[1],ls.c1.lim[2])
mu.c1~dunif(mu.c1.lim[1],mu.c1.lim[2])
beta.c1~dunif(beta.c1.lim[1],beta.c1.lim[2])
ls.e1~dunif(ls.e1.lim[1],ls.e1.lim[2])
mu.e1~dunif(mu.e1.lim[1],mu.e1.lim[2])

for(i in 1:N2){
    cost2[i]~dnorm(phi.c2[i],tau.c2)
    eff2[i]~dnorm(phi.e2[i],tau.e2)
    phi.c2[i]<-mu.c2+beta.c2*(eff2[i]-mu.e2)
    phi.e2[i]<-mu.e2
    da2[i]<-log(exp(-pow(cost2[i]-phi.c2[i],2)/(2*ss.c2))/(sqrt(2*3.1415927*ss.c2)))-log(cost.rs)
    db2[i]<-log(exp(-pow(eff2[i]-phi.e2[i],2)/(2*ss.e2))/(sqrt(2*3.1415927*ss.e2)))-log(eff.rs)
}

# node transformations
tau.c2<-1/ss.c2
ss.c2<-s.c2*s.c2
s.c2<-exp(ls.c2)
ss.c2.rs<-ss.c2*cost.rs*cost.rs
tau.e2<-1/ss.e2
ss.e2<-s.e2*s.e2
s.e2<-exp(ls.e2)
ss.e2.rs<-ss.e2*eff.rs*eff.rs

# prior distributions
ls.c2~dunif(ls.c2.lim[1],ls.c2.lim[2])
mu.c2~dunif(mu.c2.lim[1],mu.c2.lim[2])
beta.c2~dunif(beta.c2.lim[1],beta.c2.lim[2])
ls.e2~dunif(ls.e2.lim[1],ls.e2.lim[2])
mu.e2~dunif(mu.e2.lim[1],mu.e2.lim[2])

dev<--2*(sum(da1[])+sum(db1[])+sum(da2[])+sum(db2[]))

#unused variables
for(i in 1:N1){
    hos1.dum[i]<-hos1[i]
    cen1.dum[i]<-cen1[i]
}
for(i in 1:N2){
    hos2.dum[i]<-hos2[i]
    cen2.dum[i]<-cen2[i]
}
hos.bar.dum<-hos.bar
for(j in 1:4){cen.n.dum[j]<-cen.n[j]}

# ce[1]=c1, ce[2]=c2, ce[3]=e1, ce[4]=e2
ce[1]<-mu.c1*cost.rs
ce[2]<-mu.c2*cost.rs
ce[3]<-mu.e1*eff.rs
ce[4]<-mu.e2*eff.rs
}