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

# nart1 is an categorical variable taking values 1 for normal and 2 for borderline intelligence for arm 1
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]<-mu.c1+beta.c1*(eff1[i]-phi.e1[i])+gc0.nart[nart1[i]]
    phi.e1[i]<-mu.e1+ge0.nart[nart1[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])
mu.c1~dunif(mu.c1.lim[1],mu.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])
mu.e1~dunif(mu.e1.lim[1],mu.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]<-mu.c2+beta.c2*(eff2[i]-phi.e2[i])+gc0.nart[nart2[i]]+dc0.nart[nart2[i]]
    phi.e2[i]<-mu.e2+ge0.nart[nart2[i]]+de0.nart[nart2[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])
mu.c2~dunif(mu.c2.lim[1],mu.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])
mu.e2~dunif(mu.e2.lim[1],mu.e2.lim[2])

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

# covariate identifiability constraints
# nart.n[i] gives how many patients are in nart sub group i
gc0.nart[1]<- -gc0.nart[2]*nart.n[2]/nart.n[1]
ge0.nart[1]<- -ge0.nart[2]*nart.n[2]/nart.n[1]
dc0.nart[1]<- -dc0.nart[2]*nart.n[2]/nart.n[1]
de0.nart[1]<- -de0.nart[2]*nart.n[2]/nart.n[1]

# rescale covarites so relative to the mean
for(j in 1:2){gc.nart[j]<-gc0.nart[j]}#-inprod(nart.n[1:2],gc0.nart[1:2])/sum(nart.n[1:2])}
for(j in 1:2){ge.nart[j]<-ge0.nart[j]}#-inprod(nart.n[1:2],ge0.nart[1:2])/sum(nart.n[1:2])}
for(j in 1:2){dc.nart[j]<-dc0.nart[j]}#-inprod(nart.n[1:2],dc0.nart[1:2])/sum(nart.n[1:2])}
for(j in 1:2){de.nart[j]<-de0.nart[j]}#-inprod(nart.n[1:2],de0.nart[1:2])/sum(nart.n[1: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]}

# covariate rescaled output
for(j in 1:2){
    gc0.nart.rs[j]<-gc0.nart[j]*cost.rs
    gc.nart.rs[j]<-gc.nart[j]*cost.rs
}
for(j in 1:2){
    dc0.nart.rs[j]<-dc0.nart[j]*cost.rs
    dc.nart.rs[j]<-dc.nart[j]*cost.rs
}
for(j in 1:2){
    ge0.nart.rs[j]<-ge0.nart[j]*eff.rs
    ge.nart.rs[j]<-ge.nart[j]*eff.rs
}
for(j in 1:2){
    de0.nart.rs[j]<-de0.nart[j]*eff.rs
    de.nart.rs[j]<-de.nart[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
}