Hips model 3: MC estimates for
                           each strata, allowing
                        for parameter uncertainty in
                        revision hazard, h - gives
                        results for Table 3

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; monitor C, BL, BQ to get posterior mean and sd for each subgroup for results in top part of Table 3.

For results in bottom part of Table 3:

Approach 1: p(s, theta) = p(s | theta) p(theta)
m_theta and v_theta are the mean and variance across subgroups for a given value of theta
=> within BUGS code, calculate mean and variance of C[k], BL[k], BQ[k] across subgroups at each iteration, then take Monte Carlo expectation at end of run
=> monitor mean.C, mean.BL, mean.BQ, var.C, var.BL, var.BQ

Approach 2: p(s, theta) = p(theta | s) p(s)
overall mean, m = weighted mean of posterior means of C[k], BL[k], BQ[k] => calculate after BUGS run var due to uncertainty, vP2 = weighted mean of posterior variances of C[k], BL[k], BQ[k] => calculate after BUGS run var due to heterogeneity = vH2 = weighted variance of posterior means of C[k], BL[k], BQ[k] => calculate after BUGS run

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

    model {

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

       # Cost and benefit equations
       #######################

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

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

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


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

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

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

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

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

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

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

          # Marginal probability of being in each state at time 1
          pi[k,1,1] <- 1 - lambda.op pi[k,1, 2] <- 0 pi[k,1, 3] <- 0 ;
         pi[k,1, 4] <- 0 pi[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[k, t, s] <- inprod(pi[k, t - 1, ], Lambda[k, t, , s])
             }
          }
       }

       # age-sex specific revision hazard
       for(k in 1 : K) {
           logh[k] ~ dnorm(logh0[k], tau)
           h[k] <- exp(logh[k])
       }

       # Calculate mean and variance across strata at each iteration
       # (Gives overall mean and variance using approach 1)

       mean.C <- inprod(p.strata[], C[])
       mean.BL <- inprod(p.strata[], BL[])
       mean.BQ <- inprod(p.strata[], BQ[])

       for(k in 1:12) {
           C.dev[k] <- pow(C[k]-mean.C , 2)
           BL.dev[k] <- pow(BL[k]-mean.BL , 2)
           BQ.dev[k] <- pow(BQ[k]-mean.BQ , 2)
       }
       var.C <- inprod(p.strata[], C.dev[])
       var.BL <- inprod(p.strata[], BL.dev[])
       var.BQ <- inprod(p.strata[], BQ.dev[])

    }


Data ( click to open )

Results