## HEALTH ECONOMIC EVALUATION --- EXAMPLE THROUGHOUT CHAPTER 3 AND CHAPTER 4

## Loosely based on Fox-Rushby & Cairns (2005). Economic Evaluation. Open University Press
## Parameters of the model
# pi = probability of side effects (specific to each treatment)
# rho = reduction in probability of side effects for t=2
# phi = probability of ambulatory visit, given side effects (constant across treatments)
# c.amb = cost of ambulatory visit after side effects
# c.hosp = cost of hospitalisation after side effects
# c.drug = cost of each treatment (known constant)

# define the working directory assuming - the user is already in ~/WebMaterial/ch2
working.dir <- paste(getwd(),"/",sep="")
source("https://gianluca.statistica.it/BMHE/WebMaterial/Utils.R")		# loads utility functions

## Generates data for the number of side effects under t=0
N.studies <- 5
sample.size <- 32
prop.pi <- .24
n <- rpois(N.studies,sample.size)
se <- rbinom(N.studies,n,prop.pi)

## Generates the data for the number of patients with side effects requiring ambulatory care
prop.gamma <- .619
amb <- rbinom(N.studies,se,prop.gamma)

## Computes the hyperparameters for the informative prior distributions 
m.pi <- .5 
s.pi <- sqrt(.125)
a.pi <- betaPar(m.pi,s.pi)$a
b.pi <- betaPar(m.pi,s.pi)$b

m.gamma <- .5
s.gamma <- sqrt(.125)
a.gamma <- betaPar(m.gamma,s.gamma)$a
b.gamma <- betaPar(m.gamma,s.gamma)$b

m.rho <- 0.8
s.rho <- 0.2
tau.rho <- 1/s.rho^2

mu.amb <- 120
sd.amb <- 20
m.amb <- lognPar(mu.amb,sd.amb)$mulog
s.amb <- lognPar(mu.amb,sd.amb)$sigmalog
tau.amb <- 1/s.amb^2

mu.hosp <- 5500
sd.hosp <- 980
m.hosp <- lognPar(mu.hosp,sd.hosp)$mulog
s.hosp <- lognPar(mu.hosp,sd.hosp)$sigmalog
tau.hosp <- 1/s.hosp^2

c.drug <- c(110,520)

# Number of patients in the population
N <- 1000


## Run JAGS
library(R2jags)
dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","a.gamma",		# data list
		"b.gamma","m.amb","tau.amb","m.hosp",
		"tau.hosp","m.rho","tau.rho","N") 
filein <- "https://gianluca.statistica.it/BMHE/WebMaterial/ch3/model.txt"	# model file
params <- c("pi","gamma","c.amb","c.hosp","rho","SE","A","H")			# parameters list
inits <- function(){								# randomly initialise relevant variables
	SE=rbinom(2,N,.2)
	list(pi=c(runif(1),NA),gamma=runif(1),c.amb=rlnorm(1),
	c.hosp=rlnorm(1),rho=runif(1),SE=SE,A=rbinom(2,SE,.2))
}

n.iter <- 20000
n.burnin <- 9500
n.thin <- floor((n.iter-n.burnin)/250)

chemo <- jags(dataJags, inits, params, model.file=filein,
	n.chains=2, n.iter, n.burnin, n.thin,
	DIC=TRUE, working.directory=working.dir, progress.bar="text")
## NB: If initial values for SE are not compatible with the rest of the model, JAGS stops with an error.
##     The solution is to re-launch the whole model until it does work. 

print(chemo,digits=3,intervals=c(0.025, 0.975))
attach.bugs(chemo$BUGSoutput)


## Cost-effectiveness analysis
# Defines the variables of cost and effectiveness
e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
e <- N - SE
for (t in 1:2) {
	c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
}

library(BCEA)
treats <- c("Old Chemotherapy","New Chemotherapy")
m <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
summary(m)



## Computes EVPPI
# Analysis of Expected Value of Partial Information using Strong & Oakley's or Sadatsafavi et al.'s methods
inp <- CreateInputs(chemo)	# Create the inputs (parameters & matrix of MCMC simulations for all of them)
x.so <- evppi("rho",inp$mat,m,n.blocks=50)		# Uses S&O method
plot(x.so)
x.sal <- evppi("rho",inp$mat,m,n.seps=2)		# Uses Sadatsafavi et al method
plot(x.sal)
# Compare the two methods
plot(x.so$k,x.so$evi,t="l",lwd=2,xlab="Willingness to pay",ylab="",main="Expected value of partial information")
points(x.so$k,x.so$evppi,t="l",col="blue")
points(x.sal$k,x.sal$evppi,t="l",col="red")
legend("topleft",c("EVPI","EVPPI for rho (S&O)","EVPPI for rho (S et al)"),cex=.7,bty="n",lty=1,
	col=c("black","blue","red"))


# Analysis of Expected Value of Partial Information using 2stage approach
# Main parameters: 
# rho = reduction in side effects 
# gamma = chance of ambulatory care
rho.sim <- rho
gamma.sim <- gamma

detach.bugs()
inits <- NULL

# EVPPI for rho
# Replicates the C/E analysis for a fixed value of rho at each iteration
Ustar.phi <- matrix(NA,chemo$BUGSoutput$n.sims,length(m$k))
Rhat <- matrix(NA,chemo$BUGSoutput$n.sims,length(chemo$BUGSoutput$summary[,"Rhat"])-1)
for (i in 1:chemo$BUGSoutput$n.sims) {
	rho <- rho.sim[i]	

	dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","a.gamma","b.gamma","m.amb",
		"tau.amb","m.hosp","tau.hosp","N","rho") 
	filein <- "modelEVPPI_rho.txt"
	params <- c("pi","gamma","c.amb","c.hosp","SE","A","H")
	n.iter <- 10000
	n.burnin <- 9750
	n.thin <- floor((n.iter-n.burnin)/250)
	chemo.evppi <- jags(dataJags, inits, params, model.file=filein,
		n.chains=2, n.iter, n.burnin, n.thin,
		DIC=TRUE, working.directory=working.dir, progress.bar="text")
	Rhat[i,] <- chemo.evppi$BUGSoutput$summary[,"Rhat"]
	attach.bugs(chemo.evppi$BUGSoutput)

	# Performs the economic analysis given the simulations for all the
	# other parameters, conditionally on rho
	e.temp <- c.temp <- matrix(NA,chemo$BUGSoutput$n.sims,2)
	e.temp <- N - SE
	for (t in 1:2) {
		c.temp[,t] <- c.drug[Rt]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
	}
	m.evppi <- bcea(e=e.temp,c=c.temp,ref=2,interventions=treats,Kmax=50000,Ktable=25000)

	# Computes the maximum expected utility for that iteration for each value of k
	Ustar.phi[i,] <- apply(m.evppi$Ustar,2,mean)
	rm(m.evppi); detach.bugs()
	print(i)
}

# Computes the average value of the maximum expected utility for each value of k
Vstar.phi <- apply(Ustar.phi,2,mean)
Umax <- apply(apply(m$U,c(2,3),mean),1,max)

# Computes the EVPPI for each value of k
EVPPI_rho <- Vstar.phi - Umax

# Plot the results
plot(m$k,m$evi,t="l",xlab="Willingness to pay",ylab="") 
points(m$k,EVPPI_rho,t="l",lty=2)
text <- c("EVPI",expression(paste("EVPPI for ",rho,sep="")))
legend("topleft",text,lty=1:2,cex=.9,box.lty=0)

# Checks for convergence in all (inner) simulations
par(mfrow=c(4,3))
for (i in 1:12) {
	plot(Rhat[,i],t="l",main=rownames(chemo.evppi$BUGSoutput$summary)[i],ylim=c(.99,1.15))
	abline(h=1.1)
}


# EVPI for gamma
# Replicates the C/E analysis for a fixed value of rho at each iteration
Ustar.phi <- matrix(NA,chemo$BUGSoutput$n.sims,length(m$k))
Rhat <- matrix(NA,chemo$BUGSoutput$n.sims,length(chemo$BUGSoutput$summary[,"Rhat"])-1)
for (i in 1:chemo$BUGSoutput$n.sims) {
	gamma <- gamma.sim[i]	

	dataJags <- list("se","amb","n","N.studies","a.pi","b.pi","m.rho","tau.rho","m.amb",
		"tau.amb","m.hosp","tau.hosp","N","gamma") 
	filein <- "modelEVPPI_gamma.txt"
	params <- c("pi","rho","c.amb","c.hosp","SE","A","H")
	n.iter <- 10000
	n.burnin <- 9750
	n.thin <- floor((n.iter-n.burnin)/250)
	chemo.evppi <- jags(dataJags, inits, params, model.file=filein,
		n.chains=2, n.iter, n.burnin, n.thin,
		DIC=TRUE, working.directory=working.dir, progress.bar="text")
	Rhat[i,] <- chemo.evppi$BUGSoutput$summary[,"Rhat"]
	attach.bugs(chemo.evppi$BUGSoutput)

	# Performs the economic analysis given the simulations for all the
	# other parameters, conditionally on rho
	e.temp <- c.temp <- matrix(NA,chemo$BUGSoutput$n.sims,2)
	e.temp <- N - SE
	for (t in 1:2) {
		c.temp[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
	}
	m.evppi <- bcea(e=e.temp,c=c.temp,ref=2,interventions=treats,Kmax=50000,Ktable=25000)

	# Computes the maximum expected utility for that iteration for each value of k
	Ustar.phi[i,] <- apply(m.evppi$Ustar,2,mean)
	rm(m.evppi); detach.bugs()
	print(i)
}

# Computes the average value of the maximum expected utility for each value of k
Vstar.phi <- apply(Ustar.phi,2,mean)
Umax <- apply(apply(m$U,c(2,3),mean),1,max)

# Computes the EVPPI for each value of k
EVPPI_gamma <- Vstar.phi - Umax

# Plot the results
plot(m$k,m$evi,t="l",xlab="Willingness to pay",ylab="") 
points(m$k,EVPPI_gamma,t="l",lty=2)
text <- c("EVPI",expression(paste("EVPPI for ",gamma,sep="")))
legend("topleft",text,lty=1:2,cex=.9,box.lty=0)


# Checks for convergence in all (inner) simulations
par(mfrow=c(4,3))
for (i in 1:12) {
	plot(Rhat[,i],t="l",main=rownames(chemo.evppi$BUGSoutput$summary)[i],ylim=c(.99,1.15))
	abline(h=1.1)
}



## Model averaging
# M1 = complete model
# M2 = rho =: 0.3
# M3 = rho =: 1

# analysis for M2
detach.bugs()
rho <- 1

dataJags <- list("se","n","amb","N.studies","a.pi","b.pi","a.gamma","b.gamma","m.amb","tau.amb","m.hosp","tau.hosp","N","rho") 
filein <- "modelEVPPI.txt"
params <- c("pi","gamma","c.amb","c.hosp","rho","SE","A","H")
inits <- function(){
	SE=rbinom(2,N,.2)
	list(pi=c(runif(1),NA),gamma=runif(1),c.amb=rlnorm(1),c.hosp=rlnorm(1),SE=SE,A=rbinom(2,SE,.2))
}

n.iter <- 20000
n.burnin <- 9500
n.thin <- floor((n.iter-n.burnin)/250)
chemo2 <- jags(dataJags, inits, params, model.file=filein,
	n.chains=2, n.iter, n.burnin, n.thin,
	DIC=TRUE, working.directory=working.dir, progress.bar="text")
print(chemo2,digits=3,intervals=c(0.025, 0.975))
attach.bugs(chemo2$BUGSoutput,overwrite=TRUE)

## Cost-effectiveness analysis
# Defines the variables of cost and effectiveness
e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
e <- N - SE
for (t in 1:2) {
	c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
}

m2 <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
fileout <- "Model_rho_fixed.ps"
postscript(fileout,paper="special",width=6,height=6,horizontal=F)
plot(m2)
dev.off()

## Compares the models based on the DIC
d <- w <- numeric()
d[1] <- chemo$BUGSoutput$DIC
d[2] <- chemo2$BUGSoutput$DIC
dmin <- min(d)

# Computes the model weights
w <- exp(-.5*(d-dmin))/sum(exp(-.5*(d-dmin)))

# Computes the model average health economic analysis
SE <- w[1]*chemo$BUGSoutput$sims.list$SE + w[2]*chemo2$BUGSoutput$sims.list$SE
c.hosp <- w[1]*chemo$BUGSoutput$sims.list$c.hosp + w[2]*chemo2$BUGSoutput$sims.list$c.hosp
c.amb <- w[1]*chemo$BUGSoutput$sims.list$c.amb + w[2]*chemo2$BUGSoutput$sims.list$c.amb 
A <- w[1]*chemo$BUGSoutput$sims.list$A + w[2]*chemo2$BUGSoutput$sims.list$A
H <- w[1]*chemo$BUGSoutput$sims.list$H + w[2]*chemo2$BUGSoutput$sims.list$H

e <- c <- matrix(NA,chemo$BUGSoutput$n.sims,2)
e <- N - SE
for (t in 1:2) {
	c[,t] <- c.drug[t]*(N-SE[,t]) + (c.amb+c.drug[t])*A[,t] + (c.hosp+c.drug[t])*H[,t]
}
m.avg <- bcea(e=e,c=c,ref=2,interventions=treats,Kmax=50000)
plot(m.avg)

par(mfrow=c(2,2))
contour(m)
contour(m2)
contour(m.avg)