Appendix. WinBUGS and WBDiff code to implement the analysis using the Health Improvement (HI) model.
WinBUGS code
model {
for (m in 1:3) {
for (k in 1:2) {
#Estimation
# Note: The format of this part of the program will need to be altered if all trials in the network
# do not have the same comparator.
# Calculation of the transition Probabilities and
# specification of the likelihood
for (i in 1:3) {
r[m,k,i,1:4] ~ dmulti(p[m,k,i,1:4],n[m,k,i]) # likelihood
for(j in 1:4) {
# Vector containing rates to be turned into probabilities
theta[m,k,index[i,j]] <- Q[m,k,i,j]
#Transition probabilities set equal to solutions to KFE
p[m,k,i,j]<- solution[m,k,index[i,j]]
r.hat[m,k,i,j] <- p[m,k,i,j] * n[m,k,i] # Predicted data
}
}
# Function call to WBDIFF to solve Kolmogorov’s Forward Equations
solution[m,k,1:dim] <- four.state(init[1:dim],
grid[1:n.grid], theta[m,k,1:n.par], origin, tol)
# Calculation of the transition rates for a 6 transition model.
# Adjust accordingly for the 9 transition models.
Q[m,k,1,1] <- - q12[m,k]
Q[m,k,1,2] <- q12[m,k]
Q[m,k,1,3] <- 0
Q[m,k,1,4] <- 0
Q[m,k,2,1] <- q21[m,k]
Q[m,k,2,2] <- - q21[m,k] - q23[m,k] - q24[m,k]
Q[m,k,2,3] <- q23[m,k]
Q[m,k,2,4] <- q24[m,k]
Q[m,k,3,1] <- 0
Q[m,k,3,2] <- q32[m,k]
Q[m,k,3,3] <- - q32[m,k] - q34[m,k]
Q[m,k,3,4] <- q34[m,k]
# Loglinear models for the transition rates for the HI model.
# Adjust accordingly for the other treatment models.
log(q12[m,k]) <- alpha12[m]
log(q21[m,k]) <- alpha21[m] + t[m,k]
log(q23[m,k]) <- alpha23[m]
log(q24[m,k]) <- alpha24[m]
log(q32[m,k]) <- alpha32[m] + t[m,k]
log(q34[m,k]) <- alpha34[m]
}
}
# Reformat treatment parameters to coincide with notation in the paper
delta[1] <- 0
for (k in 2:4) {
delta[k] <- t[k-1,2]
}
# Priors for regression parameters
for (m in 1:3) {
alpha12[m] ~ dnorm(0,0.1)I(,10)
alpha21[m] ~ dnorm(0,0.1)I(,10)
alpha23[m] ~ dnorm(0,0.1)I(,10)
alpha24[m] ~ dnorm(0,0.1)I(,10)
alpha32[m] ~ dnorm(0,0.1)I(,10)
alpha34[m] ~ dnorm(0,0.1)I(,10)
t[m,1] <- 0
t[m,2] ~ dnorm(0,0.1)I(,10)
}
# Cost-effectiveness analysis
# Note that the code used to generate Figure 2 has been omitted as the
# charts are illustrative and not part of the primary analysis
# Cost and Utility proportions for the Exacerbation state
Pexprop ~ dbeta(115,19) # proportion of X from Pex
Hexprop <- 1 - Pexprop # proportion of X from Hex
HexERprop ~ dbeta(11,291) # proportion of HexER from Hex
HexIPprop <- 1 - HexERprop # proportion of HexIP from Hex
# Costs
# Fixed drug costs
TC[1] <- 9.20 + 0.14 # drug A
TC[2] <- 5.33 + 0.14 # drug B
TC[3] <- 10.36 + 0.14 # drug C
TC[4] <- 0.87 * 10.36 + 0.87 * 5.79 + 0.14 # drug D
# Exacerbation state costs
for (x in 1:3) {
SCX[x] ~ dnorm(msc[x],psc[x])
}
# cost of spending 1 week in each state
for (k in 1:4) {
C[k,1] <- TC[k]
C[k,2] <- TC[k]
C[k,3] <- Pexprop * SCX[1] + Hexprop * (HexERprop * SCX[2] +
HexIPprop * SCX[3]) + TC[k]
C[k,4] <- 13.66
}
# Utilities
SU[1] ~ dbeta(a[1],b[1])
SU[2] ~ dbeta(a[2],b[2])
x1u ~ dbeta(a[3],b[3])
x23u ~ dbeta(a[4],b[4])
SU[3] <- Pexprop * x1u + Hexprop * x23u
SU[4] <- 0.89
# Calculation of the calibrated transition rates
for (k in 1:4) {
Q.cal[k,1,1] <- - q12.cal[k]
Q.cal[k,1,2] <- q12.cal[k]
Q.cal[k,1,3] <- 0
Q.cal[k,1,4] <- 0
Q.cal[k,2,1] <- q21.cal[k]
Q.cal[k,2,2] <- - q21.cal[k] - q23.cal[k] - q24.cal[k]
Q.cal[k,2,3] <- q23.cal[k]
Q.cal[k,2,4] <- q24.cal[k]
Q.cal[k,3,1] <- 0
Q.cal[k,3,2] <- q32.cal[k]
Q.cal[k,3,3] <- - q32.cal[k] - q34.cal[k]
Q.cal[k,3,4] <- q34.cal[k]
log(q12.cal[k]) <- alpha12[3]
log(q21.cal[k]) <- alpha21[3] + delta[k]
log(q23.cal[k]) <- alpha23[3]
log(q24.cal[k]) <- alpha24[3]
log(q32.cal[k]) <- alpha32[3] + delta[k]
log(q34.cal[k]) <- alpha34[3]
}
# Calculation of calibrated transition probabilities for the CEA
for (k in 1:4) {
solution.cal[k,1:dim] <- four.state(init[1:dim], grid[1:n.grid],
theta.cal[k,1:n.par], origin, tol)
for (i in 1:3) {
for(j in 1:4) {
theta.cal[k,index[i,j]] <- Q.cal[k,i,j]
p.cal[k,i,j]<- solution.cal[k,index[i,j]]
}
}
}
# State occupancy times for a 52 week time horizon
for (k in 1:4) {
for (j in 1:3) {
pi.cal[k,j,1] <- start[1,1,j]
for (v in 2:52+1) { # No' of weeks for simulation
pi.cal[k,j,v] <- inprod(pi.cal[k,1:3, v-1], p.cal[k,1:3,j])
}
}
stime.cal[k,1] <- sum(pi.cal[k,1, ]) - start[1,1,1]
stime.cal[k,2] <- sum(pi.cal[k,2, ]) - start[1,1,2]
stime.cal[k,3] <- sum(pi.cal[k,3, ]) - start[1,1,3]
stime.cal[k,4] <- 52 - stime.cal[k,1] - stime.cal[k,2] - stime.cal[k,3]
# Total annual cost for treatment k
C.tot.cal[k] <- stime.cal[k,1] * C[k,1] + stime.cal[k,2] * C[k,2] +
stime.cal[k,3] * C[k,3] + stime.cal[k,4] * C[k,4]
# Total annual utility in QALY's for treatment k
U.cal[k] <- (stime.cal[k,1] * SU[1] + stime.cal[k,2] * SU[2] +
stime.cal[k,3] * SU[3] + stime.cal[k,4] * SU[4]) / 52
costdiffcal[k] <- C.tot.cal[k] - C.tot.cal[1]
bendiffcal[k] <- U.cal[k] - U.cal[1]
}
# Net Benefit and Incremental Net Benefit for different thresholds
K.space <- 0.02 # Increments of £1,000
for (o in 1:51) {
K[o] <- (o-1) * K.space * 50000
for (k in 1:4) {
INB.cal[k,o] <- K[o] * bendiffcal[k] - costdiffcal[k]
NB.cal[k,o] <- K[o] * U.cal[k] - C.tot.cal[k]
}
}
# variables to hold coda output used for model averaging
for (k in 1:4) {
NBcal0thresh[k] <- NB.cal[k,1]
NBcal5000thresh[k] <- NB.cal[k,6]
NBcal10000thresh[k] <- NB.cal[k,11]
NBcal15000thresh[k] <- NB.cal[k,16]
NBcal20000thresh[k] <- NB.cal[k,21]
NBcal22000thresh[k] <- NB.cal[k,23]
NBcal24000thresh[k] <- NB.cal[k,25]
NBcal25000thresh[k] <- NB.cal[k,26]
NBcal26000thresh[k] <- NB.cal[k,27]
NBcal28000thresh[k] <- NB.cal[k,29]
NBcal30000thresh[k] <- NB.cal[k,31]
NBcal35000thresh[k] <- NB.cal[k,36]
NBcal40000thresh[k] <- NB.cal[k,41]
NBcal45000thresh[k] <- NB.cal[k,46]
NBcal50000thresh[k] <- NB.cal[k,51]
}
# Calulation of the probabilities each of the treatments are the most cost effective (used for CEACs)
for (k in 1:4) {
for (o in 1:51) {
p.ce[k,o] <- equals(rank(INB.cal[1:4,o],k),4)
}
}
# EVPI
for (o in 1:4) {
# Maximum Expected Net Benefit for each threshold calculated in first MCMC #run
MaxINBcal[o] <- INB.cal[2,o]
}
for (o in 5:51) {
MaxINBcal[o] <- INB.cal[1,o]
}
# EVPI - only valid in the second MCMC run (which must be the same length as the first run and have # identical names as the first run for all stochastic nodes)
for (o in 1:51) {
EVPI[o] <- ranked(INB.cal[1:4,o],4) - MaxINBcal[o]
}
}
# Data (other than trial data)
List(
dim=12,origin=0,tol=1.0E-4, init=c(1,0,0,0, 0,1,0,0, 0,0,1,0),n.par=12,
n.grid=1,grid=c(1),
index=structure(.Data=c(1,2,3,4,
5,6,7,8,
9,10,11,12), .Dim=c(3,4)),
start = structure(.Data=c(
1,0,0,0, 1,0,0,0, 1,0,0,0, 1,0,0,0, 1,0,0,0, 1,0,0,0,
1,0,0,0, 1,0,0,0, 1,0,0,0, 1,0,0,0),.Dim=c(5,2,4))
msc = c(101.08,119.59,2047.24),
psc = c(0.20222,0.00218225,0.00001104),
a = c(99.01,18.19,36.37,73.74),
b = c(1,4.267,21.36,209.88))
# Initial Values
list(
alpha12 = c(-0.5,-0.5,-0.5),
alpha21 = c(-0.5,-0.5,-0.5),
alpha23 = c(-0.5,-0.5,-0.5),
alpha24 = c(-0.5,-0.5,-0.5),
alpha32 = c(-0.5,-0.5,-0.5),
alpha34 = c(-0.5,-0.5,-0.5),
t = structure(.Data=c(
NA,-5,NA,-5,NA,-5),.Dim=c(3,2)))
WBDiff code
(* code to be inserted into the main WBDiff template *)
(* code is suitable for all 4 state transition structures *)
CONST
nEq = 12;
nSt = 5; (* Number of states + 1 because vectors and Arrays rank from zero rather than 1*)
(* define index of parameters (look-up table) *)
index[1,1] := 0;
index[1,2] := 1;
index[1,3] := 2;
index[1,4] := 3;
index[2,1] := 4;
index[2,2] := 5;
index[2,3] := 6;
index[2,4] := 7;
index[3,1] := 8;
index[3,2] := 9;
index[3,3] := 10;
index[3,4] := 11;
(* define system of nEq Differential Equations *)
dCdt[index[1,1]]:= C[index[1,1]]*theta[index[1,1]] + C[index[1,2]]*theta[index[2,1]] +
C[index[1,3]]*theta[index[3,1]];
dCdt[index[1,2]]:= C[index[1,1]]*theta[index[1,2]] + C[index[1,2]]*theta[index[2,2]] +
C[index[1,3]]*theta[index[3,2]];
dCdt[index[1,3]]:= C[index[1,1]]*theta[index[1,3]] + C[index[1,2]]*theta[index[2,3]] +
C[index[1,3]]*theta[index[3,3]];
dCdt[index[1,4]]:= C[index[1,1]]*theta[index[1,4]] + C[index[1,2]]*theta[index[2,4]] +
C[index[1,3]]*theta[index[3,4]];
dCdt[index[2,1]]:= C[index[2,1]]*theta[index[1,1]] + C[index[2,2]]*theta[index[2,1]] +
C[index[2,3]]*theta[index[3,1]];
dCdt[index[2,2]]:= C[index[2,1]]*theta[index[1,2]] + C[index[2,2]]*theta[index[2,2]] +
C[index[2,3]]*theta[index[3,2]];
dCdt[index[2,3]]:= C[index[2,1]]*theta[index[1,3]] + C[index[2,2]]*theta[index[2,3]] +
C[index[2,3]]*theta[index[3,3]];
dCdt[index[2,4]]:= C[index[2,1]]*theta[index[1,4]] + C[index[2,2]]*theta[index[2,4]] +
C[index[2,3]]*theta[index[3,4]];
dCdt[index[3,1]]:= C[index[3,1]]*theta[index[1,1]] + C[index[3,2]]*theta[index[2,1]] +
C[index[3,3]]*theta[index[3,1]];
dCdt[index[3,2]]:= C[index[3,1]]*theta[index[1,2]] + C[index[3,2]]*theta[index[2,2]] +
C[index[3,3]]*theta[index[3,2]];
dCdt[index[3,3]]:= C[index[3,1]]*theta[index[1,3]] + C[index[3,2]]*theta[index[2,3]] +
C[index[3,3]]*theta[index[3,3]];
dCdt[index[3,4]]:= C[index[3,1]]*theta[index[1,4]] + C[index[3,2]]*theta[index[2,4]] +
C[index[3,3]]*theta[index[3,4]];