# Ex 9.4  Hips 2.  Evidence synthesis       
 
#  Comparative analysis of Stanmore & Charnley incorporating all evidence  

#n = 10000 updates (1 per simulated set of parameter values) are required for this model; 
#For hazard ratio estimates, monitor HR.
# monitor C.incr, BQ.incr, BL.incr ICER,  INB and P.CEA


model {

  
 # Cost-effectiveness model
 ###################### 

  for(n in 1:2) {     # loop over protheses

 # Costs
   for(t in 1:N) {
     ct[n,t] <- inprod(pi[n,t,], c[n,])/pow((1+delta.c), (t-1))
   } 
   C[n] <- C0[n] + sum(ct[n,])

   # Benefits - life expectancy
   for(t in 1:N) {
     blt[n,t] <- inprod(pi[n,t,], bl[])/pow((1+delta.b), (t-1))
   } 
   BL[n] <- sum(blt[n,])

   # Benefits - QALYs
   for(t in 1:N) {
     bqt[n,t] <- inprod(pi[n,t,], bq[])/pow((1+delta.b), (t-1))
   } 
   BQ[n] <- sum(bqt[n,])

   # Markov model probabilities:
   #######################

   # Transition matrix
   for(t in 1:N) {
     Lambda[n,t,1,1] <- 1 -  gamma[n,t] - lambda[t]
     Lambda[n,t,1,2] <- gamma[n,t] * lambda.op
     Lambda[n,t,1,3] <- gamma[n,t] *(1-lambda.op)
     Lambda[n,t,1,4] <- 0
     Lambda[n,t,1,5] <- lambda[t] 

     Lambda[n,t,2,1] <- 0
     Lambda[n,t,2,2] <- 0 
     Lambda[n,t,2,3] <- 0 
     Lambda[n,t,2,4] <- 0 
     Lambda[n,t,2,5] <- 1 

     Lambda[n,t,3,1] <- 0
     Lambda[n,t,3,2] <- 0 
     Lambda[n,t,3,3] <- 0
     Lambda[n,t,3,4] <- 1 -  lambda[t]
     Lambda[n,t,3,5] <- lambda[t]

     Lambda[n,t,4,1] <- 0
     Lambda[n,t,4,2] <- rho * lambda.op
     Lambda[n,t,4,3] <- rho * (1- lambda.op)
     Lambda[n,t,4,4] <- 1 - rho - lambda[t]
     Lambda[n,t,4,5] <- lambda[t]

     Lambda[n,t,5,1] <- 0
     Lambda[n,t,5,2] <- 0 
     Lambda[n,t,5,3] <- 0
     Lambda[n,t,5,4] <- 0
     Lambda[n,t,5,5] <- 1

     gamma[n,t] <- h[n] * (t-1)
  }

 
  # Marginal probability of being in each state at time 1
    pi[n,1,1] <- 1-lambda.op ;  pi[n,1,2]<-0  ;    pi[n,1,3] <- 0 ;  pi[n,1,4]<-0;  pi[n,1,5] <- lambda.op

  # Marginal probability of being in each state at time t>1
    for(s in 1:S) {
       for(t in 2:N) {
       pi[n,t,s] <- inprod(pi[n,(t-1),], Lambda[n,t,,s])
     }
   }
}

  # Evidence
  #########

  for (i in 1:M){                                      # loop over studies
    rC[i] ~ dbin(pC[i],nC[i])                      # number of revisions on Charnley
    rS[i] ~ dbin(pS[i],nS[i])                      # number of revisions on Stanmore
#    logit(pC[i]) <- base[i] -logHR[i]/2      # 
#    logit(pS[i]) <- base[i]+logHR[i]/2 
    cloglog(pC[i]) <- base[i] -logHR[i]/2    # cloglog does not seem to work.
   cloglog(pS[i]) <- base[i]+logHR[i]/2 
    base[i] ~ dunif(-5,0)                   # need sensible prior restriction
    # log hazard ration for ith study
    logHR[i] ~ dnorm(LHR,tauHR[i]) 
    tauHR[i]<-qualweights[i]*tauh           # precision for ith study weighted by quality weights
 } 
 LHR~ dunif(-100,100) 
 log(HR)<-LHR 
 tauh<-1/ (sigmah*sigmah) 
 sigmah~dnorm( .2, 400)I(0,)                # between-trial sd = 0.05 (prior constrained to be positive)

# age-sex specific revision hazard

 
    logh  ~ dnorm(logh0, tau)
    logh0 <- log(h0)  
       h[1]      <- exp(logh)  			# revision hazard for Charnley
      h[2] <- HR * h[1]                           # revision hazard for Stanmore


 
 # Incremental costs and benefits
 ##########################

   C.incr    <- C[2] - C[1]
BL.incr <- BL[2] - BL[1]
   BQ.incr <-  BQ[2] - BQ[1]
   ICER      <- C.incr / BQ.incr

 # Probability of cost effectiveness @ KK pounds per QALY
  # (values of KK considered range from 500 to 20000 in 500 pound increments)  (m=12 for 6000, m=20 for 10000)
  for(m in 1:40) {
      KK[m]   <-  500*m         #   Amount health care provider is willing to pay for each additional QALY  
      INB[m] <- KK[m] * BQ.incr - C.incr
      P.CEA[m] <- step(INB[m])
  }
  }

 }