[hips20]    Hips model 2: MC estimates for
                     each strata - give   results for
                     "Monte Carlo" columns in Table 2


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 patient) are required for this model; monitor C, mean.C, BL, mean.BL, BQ, mean.BQ.

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(y[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(y[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(y[k, t, ], bq[]) / pow(1 + delta.b, t - 1)
          }
          BQ[k] <- sum(bqt[k, ])


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

          # Transition matrix
          for(t in 1 : 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

          # state of each individual in strata k at time t =1
           y[k,1,1 : S] ~ dmulti(pi[k,1, ], 1)

           # state of each individual in strata k at time t > 1
           for(t in 2 : N) {
              for(s in 1:S) {
               # sampling probabilities
                 pi[k, t, s] <- inprod(y[k, t - 1, ], Lambda[k, t, , s])
              }
              y[k, t, 1 : S] ~ dmulti(pi[k, t, ], 1)
           }

       }

       # Mean of costs and benefits over strata
       #################################

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

    }


Data ( click to open )

Results


[hips21]

Overall SD for Monte Carlo estimates at bottom of Table 2 is just the weighted SD of the strata-specific Monte Carlo means