Difference between revisions of "Economic evaluation"

#http://fi.opasnet.org/fi/Special:Opasnet_Base?id=op_fi4433.pneumokokki_vaestossa
library(OpasnetUtils)
library(ggplot2)		

openv.setN(100)

if (length(vac) == 0) stop("Mitään skenaariota ei valittu")

vac <- c("No_vaccination",vac)

if(price10 == '') price10 <- 0
if(price13 == '') price13 <- 0
n_vac <- 1.8e5

vacprice <- data.frame(
	Vaccine = c("No_vaccination", "PCV10", "PCV13"),
	Result = as.numeric(c(0, price10, price13)) * n_vac
)
vacprice <- EvalOutput(Ovariable("vacprice", data = vacprice[vacprice$Vaccine %in% vac , ]))

temp <- opbase.data("Op_en6353", subset = "serotypes_in_typical_pneumococcal_vaccines")
temp$Obs <- NULL
colnames(temp)[colnames(temp) == "Result"] <- "Serotype"

serotypes <- temp[temp$Vaccine == "Existing serotypes" , "Serotype"]

userserotypes <- temp[temp$Vaccine %in% vac , ]

if(custom_vac) {
	userserotypes <- data.frame(
		Vaccine = c(rep("PCV10", length(vac_user10)), rep("PCV13", length(vac_user13))),
		Serotype = c(vac_user10, vac_user13)
	)
}

# Näyttää monimutkaiselta tuo servacin määrittely. Eikö voisi tehdä helpomminkin?
# -- Pointti on siis että kullekin käyttäjän valitsemalle rokotteelle tehdään merkintä 
# sen sisältämistä serotyypeistä 1 sisältyy 0 ei. Näin skenaariot saadaan tehtyä yksinkertaisella
# kertolaskulla (ovariable). Alla oleva koodi on täysin vektorisoitu ja kiertää siten kaksi 
# lyhyttä for looppia (R:n puolella), mikä on kieltämättä aika pieni voitto tässä tapauksessa... 

servac <- merge(
		data.frame(userserotypes, Result = 1), # Serotypes, either default or user-defined
		merge(data.frame(Vaccine = vac), data.frame(Serotype = serotypes)), # All combinations of vaccines and serotypes
		all.y = TRUE
)
servac$Result <- as.numeric(!is.na(servac$Result))
servac <- Ovariable(
		"servac", 
		data = servac
)

objects.latest("Op_en6358", code_name = "initiate") # [[:op_en:Economic evaluation]] ovariable ICER, function sumtable
objects.latest("Op_en6353", code_name = "initiate") # [[:op_en:Epidemiological modelling]] ovariables VacCar, VacIPD
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]

## Read the annual IPD and carriage incidence data.
## The 0 entries in IPD and carriage data are replaced by small values.

#IPD <- Ovariable("IPD", ddata = "Op_fi4305.pneumokokki_vaestossa")
IPD <- Ovariable("IPD", ddata = "Op_fi4433.pneumokokki_vaestossa") # [[Markunkoesivu]]

IPD@data <- IPD@data[IPD@data$Observation == "Incidence" , colnames(IPD@data) != "Observation"]

#Car <- Ovariable("Car", ddata = "Op_fi4305.pneumokokki_vaestossa")
Car <- Ovariable("Car", ddata = "Op_fi4433.pneumokokki_vaestossa") # [[Markunkoesivu]]
Car@data <- Car@data[Car@data$Observation == "Carrier" , colnames(Car@data) != "Observation"]


p_user <- q_user <- adultcarriers <- 1

p <- Ovariable("p", data = data.frame(Result = p_user))
q <- EvalOutput(Ovariable("q", data = data.frame(Result = q_user))) 
# EvalOutput must be used because q is mentioned twice in the code and there will otherwise be a merge mismatch.

## The true number of adult carriers may actually be larger than estimated. This adjusts for that.
#Car <- Car * Ovariable("adjust", data = data.frame(Age = c("Under 5", "Over 5"), Result = c(1, adultcarriers)))

#VacCar <- EvalOutput(VacCar)
VacIPD <- EvalOutput(VacIPD)

if (1==0) {
	cat("servac\n")
	oprint(summary(servac))
	
	
	cat("Number of carriers\n")
	oprint(summary(VacCar))
	cat("Incidence of invasive pneumococcal disease.\n")
	oprint(summary(VacIPD))  
}

#if("Iter" %in% colnames(VacCar@output)) N <- max(VacCar@output$Iter) else N <- 1
if("Iter" %in% colnames(VacIPD@output)) N <- max(VacIPD@output$Iter) else N <- 1

if (1==0) {ggplot(VacCar@output, aes(x = Serotype, weight = result(VacCar) / N, fill = Vaccine)) + geom_bar(position = "dodge") + theme_gray(base_size = 24) + 
			labs(title = "Carriers", y = "Number of carriers in Finland") }

ggplot(VacIPD@output, aes(x = Serotype, weight = result(VacIPD) / N, fill = Vaccine)) + geom_bar(position = "dodge") + theme_gray(base_size = 24) +
		labs(title = "Incidence of invasive pneumococcal disease", y = "Number of cases per year")

VacIPD@output$Agegroup <- cut(
	as.numeric(levels(VacIPD@output$Age[VacIPD@output$Age])), 
	breaks = c(0, 3, 5, 15, 65, 80, 101),
	include.lowest = TRUE
)
VacIPD@marginal <- c(VacIPD@marginal, FALSE)
#oprint(VacIPD)
ggplot(VacIPD@output, aes(x = Vaccine, weight = result(VacIPD) / N, fill = Agegroup)) + geom_bar(position = "stack") + theme_gray(base_size = 24) +
		labs(title = "Incidence of invasive pneumococcal disease", y = "Number of cases per year")

######################

#QALYpercase <- Ovariable("QALYpc", ddata = "Op_en6358.qalys_lost") # [[Economic evaluation]] QALYs per case

#costpercase <- Ovariable("costpc", ddata = "Op_en6358.costs_incurred") # [[Economic evaluation]] QALYs per case

#QALY <- VacIPD * QALYpercase 

#cost <- VacIPD * costpercase + vacprice

# Sum over Serotype
VacIPD <- oapply(VacIPD, NULL, sum, c("Serotype"), na.rm = TRUE)

Costs <- EvalOutput(Costs) # Healthcare costs
Total_costs <- oapply(Costs, NULL, sum, c("Outcome", "Age"))
#oprint(Total_costs)
Total_costs <- oapply(Total_costs, Total_costs@output[colnames(Total_costs@output) %in% c("Vaccine", "Iter")], mean)
health_care_costs <- Total_costs
Total_costs <- Total_costs + vacprice
Total_costs@output <- Total_costs@output[c(colnames(Total_costs@output)[colnames(Total_costs@output) %in% c("Vaccine", "Iter")], "Result")]
Total_costs@marginal <- colnames(Total_costs@output) %in% c("Vaccine", "Iter")

QALYs <- EvalOutput(QALYs)

#### Tässä voi tehdä tapauskohtaista säätöä valitsemalla sopivat indeksit.

qalyind <- "Vaccine"
if("Iter" %in% colnames(QALYs@output)) qalyind <- c(qalyind, "Iter")

#costind <- "Vaccine"
#if("Iter" %in% colnames(Total_costs@output)) costind <- c(costind, "Iter")

qalysum <- oapply(QALYs, INDEX = QALYs@output[qalyind], FUN = sum)
qalysum@name <- ""
colnames(qalysum@output)[colnames(qalysum@output) == "QALYsResult"] <- "Result"

#costsum <- oapply(Total_costs, INDEX = Total_costs@output[costind], FUN = sum)
costsum <- Total_costs

#oprint(costsum)
#oprint(qalysum)

#### The actual model

ICER <- EvalOutput(ICER)

oprint(
	oapply(VacIPD, VacIPD@output["Vaccine"], sum), 
	include.rownames = FALSE, 
	caption = "IPD incidence (cases per year)", 
	caption.placement = "top"
)

oprint(
	qalysum, 
	include.rownames = FALSE, 
	caption = "QALYs lost due to IPD", 
	caption.placement = "top"
)

oprint(
	health_care_costs, 
	include.rownames = FALSE, 
	caption = "Health care costs due to IPD", 
	caption.placement = "top"
)

oprint(
	costsum,
	include.rownames = FALSE, 
	caption = "Total costs (health care + vaccination)", 
	caption.placement = "top"
)

oprint(
	ICER, 
	include.rownames = FALSE, 
	caption = "Cost-effectiveness of vaccination choices", 
	caption.placement = "top"
)

oprint(
	sumtable(), 
	include.rownames = FALSE, 
	caption = "Summary table", 
	caption.placement = "top"
)

if (!is.null(debug_plot)) {
	temp <- QALYs
	temp <- oapply(temp, NULL, sum, "Outcome")
	temp@output$Age <- as.numeric(as.character(temp@output$Age))
	plot1 <- ggplot(
		temp@output, 
		aes(x = Age, y = QALYsResult, colour = Vaccine, group = Vaccine)
	) + geom_line() + theme_gray(base_size = 24) + labs(title = "QALYs lost due to IPD", y = "QALYs lost per year")
	# + facet_wrap(~ Outcome)
	
	temp <- Costs
	temp <- oapply(temp, NULL, sum, "Outcome")
	temp@output$Age <- as.numeric(as.character(temp@output$Age))
	plot2 <- ggplot(
		temp@output, 
		aes(x = Age, y = CostsResult, colour = Vaccine, group = Vaccine)
	) + geom_line() + theme_gray(base_size = 24) + labs(title = "IPD health care cost (excl. vaccination)", y = "")
	# + facet_wrap(~ Outcome)
	
	temp <- VacIPD
	temp@output$Age <- as.numeric(as.character(temp@output$Age))
	plot3 <- ggplot(
		temp@output, 
		aes(x = Age, y = VacIPDResult, colour = Vaccine, group = Vaccine)
	) + geom_line() + theme_gray(base_size = 24) + labs(title = "IPD incidence", y = "Cases per year")
}
if (!is.null(debug_plot)) plot3
if (!is.null(debug_plot)) plot2
if (!is.null(debug_plot)) plot1

# Rigid implementation which doesnt allow uncertainty, for debugging purposes

qorder <- qalysum@output$Vaccine[order(result(qalysum), decreasing = TRUE)]

QALYs_incremental <- c(0, -diff(result(qalysum)[match(qorder, qalysum@output$Vaccine)]))
QALYs_gained <- cumsum(QALYs_incremental)
Cost_total <- result(Total_costs)[match(qorder, Total_costs@output$Vaccine)]
Cost_incremental <- diff(c(0, Cost_total))
ICER2 <- Cost_incremental / QALYs_incremental
ICER2[1] <- 0

sum_table2 <- data.frame(
	Vaccine = qorder,
	QALYs_gained__ = QALYs_gained,
	Incremental_effect__ = QALYs_incremental,
	Health_care_costs__ = Cost_total,
	Incremental_cost__ = Cost_incremental,
	ICER__ = ICER2
)

oprint(
	sum_table2,
	include.rownames = FALSE, 
	caption = "Cost-effectiveness analysis", 
	caption.placement = "top"
)

library(OpasnetUtils)

# Initiate model components

primary_outcomes <- Ovariable("primary_outcomes", ddata = "Op_en6358.primary_outcomes")
secondary_outcomes <- Ovariable("secondary_outcomes", ddata = "Op_en6358.secondary_outcomes")
costs_per_outcomes <- Ovariable("costs_per_outcomes", ddata = "Op_en6358.costs_per_outcomes")
QALYs_per_outcomes <- Ovariable("QALYs_per_outcomes", ddata = "Op_en6358.QALYs_per_outcomes")

Outcomes <- Ovariable(
	"Outcomes", 
	dependencies = data.frame(
		Name = c("primary_outcomes", "secondary_outcomes", "VacIPD"),
		Ident = c(rep("Op_en6358/initiate", 2), "Op_en6353/initiate")
	),
	formula = function(...) {
		# Primaries
		out <- VacIPD * primary_outcomes
		
		# Secondaries
		temp <- out * secondary_outcomes
		
		# Combine outcomes under single index
		temp@output <- temp@output[!colnames(temp@output) %in% "Outcome"]
		colnames(temp@output)[colnames(temp@output) == "Outcome_new"] <- "Outcome"
		temp@output <- temp@output[colnames(temp@output) %in% colnames(out@output)]
		out <- orbind(out, temp)
		return(out)
	}
)

# Healthcare costs
Costs <- Ovariable(
	"Costs", 
	dependencies = data.frame(
		Name = c("Outcomes", "costs_per_outcomes"),
		Ident = rep("Op_en6358/initiate", 2)
	),
	formula = function(...) {
		out <- Outcomes * costs_per_outcomes
		return(out)
	}
)

# QALYs lost
QALYs <- Ovariable(
	"QALYs", 
	dependencies = data.frame(
		Name = c("Outcomes", "QALYs_per_outcomes"),
		Ident = rep("Op_en6358/initiate", 2)
	),
	formula = function(...) {
		out <- Outcomes * QALYs_per_outcomes
		return(out)
	}
)


# Initiate analysis ovariable ICER and function sumtable

ICER <- Ovariable("ICER", 
	dependencies = data.frame(Name = c(
		"qalysum", 
		"costsum",
		"QALYs"
	)),
	formula = function(...) {

		qalyorder <- oapply(QALYs, INDEX = QALYs@output["Vaccine"], FUN = sum)
		qalyorder <- as.character(qalyorder@output$Vaccine[order(result(qalyorder), decreasing = TRUE)])

		qalysum2 <- qalysum
		costsum2 <- costsum

		# Take the Vaccine group from the previous group (based on reverse QALY order, i.e. worst first.
		levels(qalysum2@output$Vaccine) <- qalyorder[match(levels(qalysum2@output$Vaccine), qalyorder) + 1]
		levels(costsum2@output$Vaccine) <- qalyorder[match(levels(costsum2@output$Vaccine), qalyorder) + 1]

		# Remove NAs from the index or otherwise they will match anything.
		qalysum2@output <- qalysum2@output[!is.na(qalysum2@output$Vaccine) , ]
		costsum2@output <- costsum2@output[!is.na(costsum2@output$Vaccine) , ]

		out <- (costsum - costsum2) / (-1 * (qalysum - qalysum2)) # The formula calls for QALY _savings_, hence * -1

		return(out)
	}
)

sumtable <- function() {
	out <- merge(
		merge(
			merge(
				qalysum@output, 
				costsum@output, by = "Vaccine"
			),
			vacprice@output, all.x = TRUE
		),
		ICER@output, all.x = TRUE
	)

	out <- out[c("Vaccine", "Result.x", "Result.y", "vacpriceResult", "ICERResult")]
	colnames(out) <- c("Vaccine", "QALY", "Costs incl. price", "Vaccination price", "ICER")
	out <- out[ order(out$QALY, decreasing = TRUE) , ]

	return(out)
}

objects.store(primary_outcomes, secondary_outcomes, costs_per_outcomes, QALYs_per_outcomes, Outcomes, Costs, QALYs, ICER, sumtable)

cat("Initiated ovariables primary_outcomes, secondary_outcomes, costs_per_outcomes, QALYs_per_outcomes, Outcomes, Costs, QALYs, ICER and function sumtable\n")

library(OpasnetUtils)

##sou#rce("F://TenderModel_R_8Aug2014//IPD_19NOV2013b//Cost_calculation.R")

cost_table <- opasnet.csv("/0/0e/Pneumococcus_cost_table.csv", wiki = "opasnet_en")
#cost_table<-re#ad.table("Cost_Table.dat")
## 101*8 taulukko

## Title of cost_table:
##  QALY losses and medical costs per case, separately for meningitis and bacteremia. 
##  (Note: QALY losses and costs for meningitis cases include sequlae.)


##Columns of  cost_table :
#1# Age (years)
age<-cost_table[,1]
#2# QALYs lost due to one meningitis case (incl. sequlae)
QALY_men<-cost_table[,2]
#3# QALYs lost due to one bacteremia case
QALY_bac<-cost_table[,3]
#4# case-fatality ratio for a meningitis or bacteremia case (ie for an IPD case)
CFR<-cost_table[,4]
#5# life years lost per one fatal IPD case
LYL<-cost_table[,5]
#6# Medical costs due to one meningitis case (including sequlae)
COST_men<-cost_table[,6]
#7# Medical costs due to one bacteremia case
COST_bac<-cost_table[,7]
#8# Proportion of meningitis cases among all IPD cases (rest are bacteremia)
PROP_men<-cost_table[,8]

## Tässä  koodissa "Cost_calculation.R" luetaan taulukko "Cost_Table.dat" ja muunnetaan 
## se taukukoksi "Loss_per_IPDcase" vastaamaan yhtä IPD tapausta. 
##
## Tällöin kust.vaik.-mallin antamat tulokset saadaan funktiossa 
## "calc_qalys_and_med_costs" kun argumentiksi annetaan IPD tapausten määrät 
## Suomessa ikävuosittain (101 kpl). Nämä IPD tapausten määrät vastaavat joko 
## "ei rokoteta" tilannetta tai lasketaan epidemiologisen mallin avulla eri 
## rokotevaihtiehdoille. (opasnetissä IPD-vektorit saadaan siis ovariablien kautta).
##
## Funktio "calc_3_ouput_tables" tuottaa 3 tulostaulukkoa. 
## Nämä ovat kust.vaik.-mallin lopputulokset.

##                           Markku Nurhonen 15.8.2014
######################################################################################




## Adjust matrix "Loss_per_case"  to correspond to one ipd case
## (instead of just meningitis or bacterremia case)
onevec<-rep(1,101)
adjustment<-cbind(onevec,PROP_men,(onevec-PROP_men),onevec,CFR,PROP_men,(onevec-PROP_men),onevec)
Loss_per_case<-cbind(age,QALY_men,QALY_bac,CFR,LYL,COST_men,COST_bac,PROP_men)
Loss_per_IPDcase<-Loss_per_case*adjustment

## Matriisia Loss_per_IPDcase käytetään päivitettäessä
## kustannuksia ja QALY-arvoja IPD insidenssien muuttuessa
## rokotteiden vaihtuessa

calc_qalys_and_med_costs<-function(ipd_novacc,ipd,Loss_per_IPDcase)
	## for two given 101-long IPD vectors
	## ipd_novacc = ipd under NO vaccination
	## ipd        = ipd under vaccination
	## this function gives a list of 
	## non-fatal,fatal and total QALYs gained:  result[[1]]:(1,2,3)
	## and medical costs under novacc and vacc: result[[2]]:(1,2)
	## Loss_per_IPDcase is a 101*8 matrix
{
	Loss_total_novacc<-matrix(ipd_novacc,101,8)*Loss_per_IPDcase
	Loss_total<-matrix(ipd,101,8)*Loss_per_IPDcase
	Gain<-apply(Loss_total_novacc-Loss_total,2,sum)  ##koko populaatio
	## Now columns 2+3 are nonfatal, 5 is fatal QALYs
	## list Qalys gained: nonfatal, fatal and total
	QALYs<-c(Gain[2]+Gain[3], Gain[5], Gain[2]+Gain[3]+Gain[5])
	## Now columns 6+7 are medical costs
	## list med cost under novacc and vacc
	medical_cost0<-cbind(Loss_total_novacc[,6]+Loss_total_novacc[,7],Loss_total[,6]+Loss_total[,7])
	medical_cost<-apply(medical_cost0,2,sum)
	list(QALYs,medical_cost)
}


calc_3_output_tables<-function(ipd0,ipd1,ipd2,vaccine_cost1,vaccine_cost2,Loss_per_IPDcase)
	## for 3 given 101-long IPD vectors
	## ipd0 = ipd under NO vaccination
	## ipd1= ipd under vaccination 1
	## ipd1= ipd under vaccination 2
	## and
	## vaccine_cost1,vaccine_cost2=
	##   per dose costs of vaccines 1 and 2
	## Loss_per_IPDcase is a 101*8 matrix
	##
	## calculate a list of 3 output tables
	##  rows and columns as indicated below
	##
	## typical call of this function:
	## calc_3_ouput_tables(IPD_noVac,IPD_pcv10,IPD_pcv13,20,40,Loss_per_IPDcase)
{
	c1<-calc_qalys_and_med_costs(ipd0,ipd1,Loss_per_IPDcase)
	c2<-calc_qalys_and_med_costs(ipd0,ipd2,Loss_per_IPDcase)

	## output table 1
	## columns(3): vaccination, non fatal, fatal and total qalys gained
	## rows: no_vacc, vacc1, vacc2
	table1<-rbind(rep(0,3),c1[[1]],c2[[1]])
	qalys_gained<-table1[,3]

	## output table 2
	## columns(3): medical costs, vaccination programme costs, health care costs
	##rows: no_vacc, vacc1, vacc2
	vaccine_cost_tot<-180000*c(0,vaccine_cost1,vaccine_cost2)
	med_cost<-c(c1[[2]],c2[[2]][2])
	healthcare_cost<-med_cost+vaccine_cost_tot
	table2<-cbind(med_cost,vaccine_cost_tot,healthcare_cost)

	## ouput table3
	## columns(5):   1.QALYs gained compared to no_vacc
	##               2.incremental effects (=incremental QALYS gained)
	##               3.Health care costs  4.incremental costs
	##	 	  	     5.ICER=column4/column2
	##rows: no_vacc, vacc1, vacc2

	incr_qalys<-(c(qalys_gained,0)-c(0,qalys_gained))[seq(3)]
	incr_costs<-(c(healthcare_cost,0)-c(0,healthcare_cost))[seq(3)]
	table3<-cbind(qalys_gained,incr_qalys,healthcare_cost,incr_costs,c(0,incr_costs[-1]/incr_qalys[-1]))

	list(table1,table2,table3)
} 

objects.store(age, QALY_men, QALY_bac, CFR, LYL, COST_men, COST_bac, PROP_men, onevec, adjustment, Loss_per_case, 
	Loss_per_IPDcase, calc_qalys_and_med_costs, calc_3_output_tables
)

cat("Objects age, QALY_men, QALY_bac, CFR, LYL, COST_men, COST_bac, PROP_men, onevec, adjustment, Loss_per_case, 
	Loss_per_IPDcase, calc_qalys_and_med_costs, calc_3_output_tables successfully stored.\n"
)