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

#cen1 is a categorical variable taking the value 1,2,3 or 4 denoting the centre a patient is in
model{
for(i in 1:N1){
    cost1[i]~dgamma(shape.c1,rate.c1[i])
    rate.c1[i]<-shape.c1/phi.c1[i]
    eff1[i]~dgamma(shape.e1,rate.e1[i])
    rate.e1[i]<-shape.e1/phi.e1[i]
    phi.c1[i]<-nu.c1+beta.c1*(eff1[i]-phi.e1[i])+gc0.cen[cen1[i]]
    phi.e1[i]<-nu.e1+ge0.cen[cen1[i]]
    da1[i]<-log(1/exp(loggam(shape.c1))*pow(rate.c1[i],shape.c1)*pow(cost1[i],shape.c1-1)*exp(-cost1[i]*rate.c1[i]))-log(cost.rs)
    db1[i]<-log(1/exp(loggam(shape.e1))*pow(rate.e1[i],shape.e1)*pow(eff1[i],shape.e1-1)*exp(-eff1[i]*rate.e1[i]))-log(eff.rs)
}

# prior distributions
shape.c1~dunif(shape.c1.lim[1],shape.c1.lim[2])
nu.c1~dunif(nu.c1.lim[1],nu.c1.lim[2])
beta.c1~dunif(beta.c1.lim[1],beta.c1.lim[2])
shape.e1~dunif(shape.e1.lim[1],shape.e1.lim[2])
nu.e1~dunif(nu.e1.lim[1],nu.e1.lim[2])

for(i in 1:N2){
    cost2[i]~dgamma(shape.c2,rate.c2[i])
    rate.c2[i]<-shape.c2/phi.c2[i]
    eff2[i]~dgamma(shape.e2,rate.e2[i])
    rate.e2[i]<-shape.e2/phi.e2[i]
    phi.c2[i]<-nu.c2+beta.c2*(eff2[i]-phi.e2[i])+gc0.cen[cen2[i]]+dc0.cen[cen2[i]]
    phi.e2[i]<-nu.e2+ge0.cen[cen2[i]]+de0.cen[cen2[i]]
    da2[i]<-log(1/exp(loggam(shape.c2))*pow(rate.c2[i],shape.c2)*pow(cost2[i],shape.c2-1)*exp(-cost2[i]*rate.c2[i]))-log(cost.rs)
    db2[i]<-log(1/exp(loggam(shape.e2))*pow(rate.e2[i],shape.e2)*pow(eff2[i],shape.e2-1)*exp(-eff2[i]*rate.e2[i]))-log(eff.rs)
}


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

# random effects
for(j in 1:4){
    de0.cen[j]~dnorm(0,tau.de.cen)
    dc0.cen[j]~dnorm(0,tau.dc.cen)
}

tau.de.cen<-1/ss.de.cen
ss.de.cen<-s.de.cen*s.de.cen
ss.de.cen.rs<-ss.de.cen*eff.rs*eff.rs

tau.dc.cen<-1/ss.dc.cen
ss.dc.cen<-s.dc.cen*s.dc.cen
ss.dc.cen.rs<-ss.dc.cen*cost.rs*cost.rs

s.de.cen~dunif(s.de.cen.lim[1],s.de.cen.lim[2])
s.dc.cen~dunif(s.dc.cen.lim[1],s.dc.cen.lim[2])

# covariate prior distributions
for(j in 2:4){gc0.cen[j]~dunif(gc0.cen.lim[j,1],gc0.cen.lim[j,2])}
for(j in 2:4){ge0.cen[j]~dunif(ge0.cen.lim[j,1],ge0.cen.lim[j,2])}

# adjustment to the intercept to get population mean parameters
ed.adj<-inprod(cen.n[1:4],de0.cen[1:4])/sum(cen.n[1:4])
cd.adj<-inprod(cen.n[1:4],dc0.cen[1:4])/sum(cen.n[1:4])
mu.e1<-nu.e1
mu.c1<-nu.c1
mu.e2<-nu.e2+ed.adj
mu.c2<-nu.c2+cd.adj

# covariate identifiability constraints
# cen.n[i] gives how many patients are in centre i
gc0.cen[1]<--inprod(cen.n[2:4],gc0.cen[2:4])/cen.n[1]
ge0.cen[1]<--inprod(cen.n[2:4],ge0.cen[2:4])/cen.n[1]

# rescale covarites so relative to the mean
for(j in 1:4){gc.cen[j]<-gc0.cen[j]}#-inprod(cen.n[1:4],gc0.cen[1:4])/sum(cen.n[1:4])}
for(j in 1:4){ge.cen[j]<-ge0.cen[j]}#-inprod(cen.n[1:4],ge0.cen[1:4])/sum(cen.n[1:4])}
for(j in 1:4){dc.cen[j]<-dc0.cen[j]-inprod(cen.n[1:4],dc0.cen[1:4])/sum(cen.n[1:4])}
for(j in 1:4){de.cen[j]<-de0.cen[j]-inprod(cen.n[1:4],de0.cen[1:4])/sum(cen.n[1:4])}

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

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

# covariate rescaled output
for(j in 1:4){
    gc0.cen.rs[j]<-gc0.cen[j]*cost.rs
    gc.cen.rs[j]<-gc.cen[j]*cost.rs
}
for(j in 1:4){
    dc0.cen.rs[j]<-dc0.cen[j]*cost.rs
    dc.cen.rs[j]<-dc.cen[j]*cost.rs
}
for(j in 1:4){
    ge0.cen.rs[j]<-ge0.cen[j]*eff.rs
    ge.cen.rs[j]<-ge.cen[j]*eff.rs
}
for(j in 1:4){
    de0.cen.rs[j]<-de0.cen[j]*eff.rs
    de.cen.rs[j]<-de.cen[j]*eff.rs
}

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