+ Show code- Hide code
#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 = c(0, price10, price13)
)
vacprice <- EvalOutput(Ovariable("vacprice", data = vacprice[vacprice$Vaccine %in% vac , ])) * n_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 = "Figure 1. Number of IPD cases per year, by serotype.", 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 = "Figure 2. Number of IPD cases per year, by age group.", 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)
if (1==2) {
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 cases per year", 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 <- c(0,diff( Cost_total))
ICER2 <- Cost_incremental / QALYs_incremental
ICER2[1] <- 0
if (1==2) {
oprint(
oapply(VacIPD, VacIPD@output["Vaccine"], sum),
include.rownames = FALSE,
caption = "Table 1. Number of cases of invasive pneumococcal disease (IPD) per year (see also Figures 1-2 below).",
caption.placement = "top"
)
}
vaccres<-matrix(result(VacIPD),101,3)[,c(3,1,2)]
ipdsums<-apply(vaccres,2,sum)
ipdtable<-data.frame(Vaccination_____=c("No vaccination ","PCV10 ","PCV13 "),N_of_IPD_cases____=round(ipdsums))
oprint(ipdtable,
include.rownames = FALSE,
caption = "Table 1. Number of cases of invasive pneumococcal disease (IPD) per year (see also Figures 1-2 below).",
caption.placement = "top"
)
##############################
## print healt care costs table
sum_table1A <- data.frame(
Vaccine__ = qorder,
Medical_costs__ = 0.01*round((result(health_care_costs)/1E4)[match(qorder,health_care_costs@output$Vaccine)]),
Vaccine_programme_cost__ = 0.01*round(result(vacprice)/1E4),
Health_care_costs__ = 0.01*round((result(costsum)/1E4)[match(qorder,costsum@output$Vaccine)])
)
oprint(
sum_table1A,
include.rownames = FALSE,
caption = "Table 2. Health care costs (in MEUR)",
caption.placement = "top"
)
##############################
## print summary table
tekstia<-data.frame(Columns=c(" 1 Vaccine ",
" 2 QALYs gained ",
" 3 Incremental effect ",
" 4 Health-case costs ",
" 5 Incremental cost ",
" 6 ICER ",
" "),
Content=c("vaccination programme",
"QALYs gained in the Finnish population (*) as compared to 'no vaccination'",
"difference in QALYs gained",
"medical costs due to IPD in the Finnish population(*) plus the cost of vaccination (in MEUR, 180000 doses) ",
"health-care cost difference (in MEUR)",
"incremental cost-effectiveness ratio (in euros). The programme with the lower ICER is identified as the more cost-effective",
"(*) QALYs and health-care costs refer to the Finnish population of 5.4 million individuals"))
oprint(tekstia, include.rownames = FALSE, include.colnames = FALSE,
caption = "Columns appearing in Table 3 (below)",
caption.placement = "top")
sum_table2 <- data.frame(
Vaccine = qorder,
QALYs_gained__ = round(QALYs_gained),
Incremental_effect__ = round(QALYs_incremental),
Health_care_costs__ = 0.01*round(Cost_total/1E4),
Incremental_cost__ = 0.01*round(Cost_incremental/1E4),
ICER__ = ICER2
)
oprint(
sum_table2,
include.rownames = FALSE,
caption = "Table 3. Cost-effectiveness analysis summary table ",
caption.placement = "top"
)
| |