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(vac, "No_vaccination")

if(price10 == '') price10 <- 0
if(price13 == '') price13 <- 0

vacprice <- data.frame(
	Vaccine = c("No_vaccination", "PCV10", "PCV13"),
	Result = c(0, price10, price13)
)
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 (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
)

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

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

qalyind <- "Vaccine"; if("Iter" %in% colnames(QALY@output)) qalyind <- c(qalyind, "Iter")
costind <- "Vaccine"; if("Iter" %in% colnames(cost@output)) costind <- c(costind, "Iter")

qalysum <- oapply(QALY, INDEX = QALY@output[qalyind], FUN = sum)
costsum <- oapply(cost, INDEX = cost@output[costind], FUN = sum)

#### The actual model

ICER <- EvalOutput(ICER)

cat("QALYs after different vaccination choices\n")

oprint(qalysum)

cat("Costs after different vaccination choices\n")

oprint(costsum)

cat("Cost-effectiveness after different vaccination choices\n")

oprint(ICER)

cat("Summary table\n")

oprint(sumtable())

library(OpasnetUtils)

#### Initiate ovariable ICER and function sumtable

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

		qalyorder <- oapply(QALY, INDEX = QALY@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) / (qalysum - qalysum2)

		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(ICER, sumtable)

cat("Initiated ovariable ICER and function sumtable\n")