Hips model 4: Comparative
                        analysis of Stanmore &
                        Charnley incorporating evidence

Spiegelhalter, D.J. and Best, N.G. “Bayesian approaches to multiple sources of evidence and uncertainty in complex cost-effectiveness modelling”. Statistics in Medicine 22 , (2003), 3687-3709.

n = 10000 updates (1 per simulated set of parameter values) are required for this model;
For hazard ratio estimates in bottom of table 4, monitor HR. For results in table 5, monitor C.incr, BQ.incr, ICER.strata, mean.C.incr, mean.BQ.incr, mean.ICER, P.CEA.strata[30,],
P.CEA.strata[50,], P.CEA[30] and P.CEA[50]. To produce plots in Fig 2, use coda option to save samples of C.incr, BQ.incr, mean.C.incr, mean.BQ.incr. To produce plots in Fig 3, set summary monitors on P.CEA.strata and P.CEA to get posterior means

Sections of the code that have changed from Model 1 are shown in bold


    model {

    # 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
          cloglog(pC[i]) <- base[i] - logHR[i]/2
          cloglog(pS[i]) <- base[i] + logHR[i]/2
          base[i] ~ dunif(-100,100)
          # log hazard ratio 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( 0.2, 400)C(0, ) # between-trial sd = 0.05 (prior constrained to be positive)

       for(k in 1 : K) {
          logh[k] ~ dnorm(logh0[k], tau)
          h[1, k] <- exp(logh[k]) # revision hazard for Charnley
          h[2, k] <- HR * h[1, k] # revision hazard for Stanmore
       }

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

       for(k in 1 : K) { # loop over strata

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

          # Cost and benefit equations in closed form:
          ####################################

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

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

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

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

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

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

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

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

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

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

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

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

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

       for(k in 1 : K) {
          C.incr[k] <- C[2, k] - C[1, k]
          BQ.incr[k] <-BQ[2, k] - BQ[1, k]
          ICER.strata[k] <- C.incr[k] / BQ.incr[k]
       }

       # Probability of cost effectiveness @ KK pounds per QALY
       # (values of KK considered range from 200 to 20000 in 200 pound increments)
       for(m in 1 : 100) {
          for(k in 1 : 12) {
             P.CEA.strata[m,k] <- step(KK[m] * BQ.incr[k] - C.incr[k])
          }
          P.CEA[m] <- step(KK[m] * mean.BQ.incr - mean.C.incr)
       }

       # overall incremental costs and benefit
       for(n in 1 : 2) {
          mean.C[n] <- inprod(p.strata[], C[n, ])
          mean.BQ[n] <- inprod(p.strata[], BQ[n, ])
       }
       mean.C.incr <- mean.C[2] - mean.C[1]
       mean.BQ.incr <- mean.BQ[2] - mean.BQ[1]
       mean.ICER <- mean.C.incr / mean.BQ.incr

    }


Data ( click to open )

Inits  for chain 1    Inits  for chain 2 ( click to open )



Results


(quality weights c(0.5, 1, 0.2), delta.c = 0.06, delta.b = 0.06)