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)T(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 }
   Results

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