#source: Price & Briggs (2002), working.dir <- paste(getwd(),"/",sep="") setwd(working.dir) # Load all the relevant data S <- 5 # number of clinical states # s=1: STW # s=2: UTW # s=3: Pex # s=4: Hex # s=5: TF J <- 12 # number of weeks in the follow up # Observed data on transitions among the states for the two treatments r.0 <- (matrix(c( 66,32,0,0,2, 42,752,0,5,20, 0,4,0,1,0, 0,0,0,0,0, 0,0,0,0,156),c(S,S),byrow=TRUE)) r.1 <- (matrix(c( 210,60,0,1,1, 88,641,0,4,13, 1,0,0,0,1, 0,0,0,0,0, 0,0,0,0,81),c(S,S),byrow=TRUE)) # Includes only the first three rows of data (Hex has no data and TP is modelled deterministically) r.0 <- r.0[1:3,] r.1 <- r.1[1:3,] n.0 <- apply(r.0,1,sum) # total number of patients in each state for t=0 n.1 <- apply(r.1,1,sum) # total number of patients in each state for t=1 # Parameters of the Dirichlets prior distributions scale <- 10 # scale factor to use different version of minimally informative prior alpha.0 <- rep(scale,S) alpha.1 <- rep(scale,S) # Run JAGS library(R2jags) dataJags <- list("n.0","n.1","r.0","r.1","alpha.0","alpha.1","S") filein <- "MarkovModel.txt" params <- c("lambda.0","lambda.1") # lambda.0 and lambda.1 need: 4 random rows (even Hex has random parameters) # 1 row of NAs (because TP is deterministic) inits <- function(){ temp.0 <- matrix(rgamma(4*S,scale,1),4,S) sum.temp.0 <- apply(temp.0,1,sum) mat.0 <- temp.0/sum.temp.0 temp.1 <- matrix(rgamma(4*S,scale,1),4,S) sum.temp.1 <- apply(temp.1,1,sum) mat.1 <- temp.1/sum.temp.1 list(lambda.0=rbind(mat.0,rep(NA,S)),lambda.1=rbind(mat.1,rep(NA,S))) } n.iter <- 10000 n.burnin <- 5000 n.thin <- floor((n.iter-n.burnin)/500) mm1 <- jags(data=dataJags,inits=inits,parameters.to.save=params,model.file=filein, n.chains=2, n.iter, n.burnin, n.thin, DIC=TRUE, working.directory=working.dir, progress.bar="text") print(mm1,digits=3,intervals=c(0.025, 0.975)) attach.bugs(mm1\$BUGSoutput) # Markov model start <- c(1,0,0,0,0) # NB analysis for 1 single patient! m.0 <- m.1 <- array(NA,c(n.sims,S,(J+1))) for (s in 1:S){ m.0[,s,1] <- start[s] m.1[,s,1] <- start[s] } for (i in 1:n.sims) { for (j in 2:(J+1)){ for (s in 1:S){ m.0[i,s,j] <- sum(m.0[i,,j-1]*lambda.0[i,,s]) m.1[i,s,j] <- sum(m.1[i,,j-1]*lambda.1[i,,s]) } } } # Costs unit.cost.0 <- c(2.38,2.38,1815.58,95.21) unit.cost.1 <- c(7.96,7.96,1821.17,100.79) cost0 <- cost1 <- matrix(NA,n.sims,J) for (i in 1:n.sims) { for (j in 2:(J+1)) { cost0[i,j-1] <- m.0[i,S,j]*(unit.cost.0%*%m.0[i,1:(S-1),j])/sum(m.0[i,1:(S-1),j]) + unit.cost.0%*%m.0[i,1:(S-1),j] cost1[i,j-1] <- m.1[i,S,j]*(unit.cost.0%*%m.0[i,1:(S-1),j])/sum(m.0[i,1:(S-1),j]) + unit.cost.1%*%m.1[i,1:(S-1),j] } } ## General formulation to apply discount delta.b <- 0.035 # discount rate for benefits (3.5%) delta.c <- 0.035 # discount rate for costs (3.5%) # Defines the discount factors disc.b <- numeric(); disc.c <- numeric() disc.b[1] <- 1; disc.c[1] <- 1 for (j in 2:J) { disc.b[j] <- (1+delta.b)^(j-1) disc.c[j] <- (1+delta.c)^(j-1) } disc.cost0 <- disc.eff0 <- disc.cost1 <- disc.eff1 <- matrix(NA,n.sims,J) for (j in 1:J) { disc.cost0[,j] <- cost0[,j]/disc.c[j] disc.cost1[,j] <- cost1[,j]/disc.c[j] disc.eff0[,j] <- m.0[,1,j]/disc.b[j] disc.eff1[,j] <- m.1[,1,j]/disc.b[j] } # Sums the values across all time points and creates matrix of costs c <- matrix(NA,n.sims,2) c[,1] <- apply(cost0,1,sum) c[,2] <- apply(cost1,1,sum) # Effectiveness e <- matrix(NA,n.sims,2) e[,1] <- apply(m.0[,1,],1,sum) e[,2] <- apply(m.1[,1,],1,sum) # Cost-effectiveness analysis library(BCEA) ints <- c("FP","SFC") m <- bcea(e,c,ref=2,interventions=ints,Kmax=300) contour2(m,300)