|
|
Line 212: |
Line 212: |
| ) | | ) |
| VacIPD@marginal <- c(VacIPD@marginal, FALSE) | | 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) + | | 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") | | 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 | | # Sum over Serotype |
Line 231: |
Line 223: |
| Costs <- EvalOutput(Costs) # Healthcare costs | | Costs <- EvalOutput(Costs) # Healthcare costs |
| Total_costs <- oapply(Costs, NULL, sum, c("Outcome", "Age")) | | 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) | | Total_costs <- oapply(Total_costs, Total_costs@output[colnames(Total_costs@output) %in% c("Vaccine", "Iter")], mean) |
| health_care_costs <- Total_costs | | health_care_costs <- Total_costs |
Line 239: |
Line 230: |
|
| |
|
| QALYs <- EvalOutput(QALYs) | | QALYs <- EvalOutput(QALYs) |
|
| |
|
| |
|
| |
|
| #### Tässä voi tehdä tapauskohtaista säätöä valitsemalla sopivat indeksit. | | #### Tässä voi tehdä tapauskohtaista säätöä valitsemalla sopivat indeksit. |
Line 246: |
Line 235: |
| qalyind <- "Vaccine" | | qalyind <- "Vaccine" |
| if("Iter" %in% colnames(QALYs@output)) qalyind <- c(qalyind, "Iter") | | 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 <- oapply(QALYs, INDEX = QALYs@output[qalyind], FUN = sum) |
Line 254: |
Line 240: |
| colnames(qalysum@output)[colnames(qalysum@output) == "QALYsResult"] <- "Result" | | colnames(qalysum@output)[colnames(qalysum@output) == "QALYsResult"] <- "Result" |
|
| |
|
| #costsum <- oapply(Total_costs, INDEX = Total_costs@output[costind], FUN = sum)
| |
| costsum <- Total_costs | | costsum <- Total_costs |
|
| |
| #oprint(costsum)
| |
| #oprint(qalysum)
| |
|
| |
|
| #### The actual model | | #### The actual model |
|
| |
|
| ICER <- EvalOutput(ICER) | | 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)) { | | if (!is.null(debug_plot)) { |
Line 343: |
Line 287: |
| ICER2[1] <- 0 | | ICER2[1] <- 0 |
|
| |
|
| if (1==2) {
| | ipdtable <- oapply(VacIPD, VacIPD@output["Vaccine"], sum)@output |
| oprint(
| | colnames(ipdtable)[colnames(ipdtable) == "VacIPDResult"] <- "N_of_IPD_cases" |
| 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( | | oprint( |
| ipdtable, | | ipdtable[order(match(ipdtable$Vaccine, qorder)),], |
| sortable = FALSE, | | sortable = FALSE, |
| include.rownames = FALSE, | | include.rownames = FALSE, |
| caption = "Table 1. Number of cases of invasive pneumococcal disease (IPD) per year (see also Figures 1-2 below).", | | caption = "Table 1. Number of cases of invasive pneumococcal disease (IPD) per year (see also Figures 1-2 below).", |
| caption.placement = "top" | | caption.placement = "top", |
| | digits = rep(0, ncol(ipdtable) + 1) |
| ) | | ) |
|
| |
|
| |
|
| |
|
| ############################## | | ############################## |
| ## print healt care costs table | | ## print health care costs table |
|
| |
|
| sum_table1A <- data.frame( | | sum_table1A <- data.frame( |
| Vaccine = qorder, | | Vaccine = qorder, |
| Medical_costs = 0.01*round((result(health_care_costs)/1E4)[match(qorder,health_care_costs@output$Vaccine)]), | | Medical_costs = result(health_care_costs)[match(qorder,health_care_costs@output$Vaccine)] * 1e-6, |
| Vaccine_programme_cost = 0.01*round(result(vacprice)/1E4), | | Vaccine_programme_cost = result(vacprice) * 1e-6, |
| Health_care_costs = 0.01*round((result(costsum)/1E4)[match(qorder,costsum@output$Vaccine)]) | | Health_care_costs = result(costsum)[match(qorder,costsum@output$Vaccine)] * 1e-6 |
| ) | | ) |
| oprint( | | oprint( |
Line 381: |
Line 313: |
| include.rownames = FALSE, | | include.rownames = FALSE, |
| caption = "Table 2. Health care costs (in MEUR)", | | caption = "Table 2. Health care costs (in MEUR)", |
| caption.placement = "top" | | caption.placement = "top", |
| | digits = c(0,0,2,2,2) |
| ) | | ) |
|
| |
|
Line 404: |
Line 337: |
| "(*) QALYs and health-care costs refer to the Finnish population of 5.4 million individuals")) | | "(*) QALYs and health-care costs refer to the Finnish population of 5.4 million individuals")) |
|
| |
|
| oprint(tekstia, include.rownames = FALSE, include.colnames = FALSE, | | oprint( |
| caption = "Columns appearing in Table 3 (below)", | | tekstia, |
| caption.placement = "top") | | include.rownames = FALSE, |
| | include.colnames = FALSE, |
| | caption = "Columns appearing in Table 3 (below)", |
| | caption.placement = "top" |
| | ) |
|
| |
|
|
| |
|
Line 412: |
Line 349: |
| sum_table2 <- data.frame( | | sum_table2 <- data.frame( |
| Vaccine = qorder, | | Vaccine = qorder, |
| QALYs_gained = round(QALYs_gained), | | QALYs_gained = QALYs_gained, |
| Incremental_effect = round(QALYs_incremental), | | Incremental_effect = QALYs_incremental, |
| Health_care_costs = 0.01*round(Cost_total/1E4), | | Health_care_costs = Cost_total * 1e-6, |
| Incremental_cost = 0.01*round(Cost_incremental/1E4), | | Incremental_cost = Cost_incremental * 1e-6, |
| ICER = ICER2 | | ICER = ICER2 |
| ) | | ) |
Line 424: |
Line 361: |
| include.rownames = FALSE, | | include.rownames = FALSE, |
| caption = "Table 3. Cost-effectiveness analysis summary table ", | | caption = "Table 3. Cost-effectiveness analysis summary table ", |
| caption.placement = "top" | | caption.placement = "top", |
| | digits = c(0,0,0,0,2,2,6) |
| ) | | ) |
| </rcode> | | </rcode> |
Progression class
In Opasnet many pages being worked on and are in different classes of progression. Thus the information on those pages should be regarded with consideration. The progression class of this page has been assessed:
- This page is a full draft
- This page has been written through once, so all important content is already where it should be. However, the content has not been thoroughly checked yet, and for example important references might still be missing.
|
This page needs a curator. Learn more about curating Opasnet pages.
|
Question
How to identify the most cost-effective pneumococcal conjugate vaccine to the national immunisation programme?
- The health benefit (effectiveness) of the pneumococcal infant immunisation programme is assessed by the expected gain in Quality-Adjusted Life Years (QALYs), corresponding to the expected reduction in the annual number of invasive pneumococcal disease in the whole Finnish population.
- The perspective of the analysis is that of the health care provider.
- The analysis is based on incremental cost effectiveness
Answer
The answer to the question is based on the concept of incremental costs. For example, if there are only two vaccines to be compared, the more effective (and more expensive vaccine) is said to be more cost-effective if the incremental cost effectiveness ratio (ICER), comparing the vaccine to the less effective vaccine, exceeds the ICER of the less effective vaccine as compared to the alternative 'no vaccination'. The principle in general is explained below (see 'Rationale').
Computation
The following programme can be used to calculate the incremental cost effectiveness ratios (ICERs) for
two alternative vaccination programmes. The input required is:
(a) the serotype compositions of the two vaccines to be compared (the defaults are PCV10 and PCV13), and
(b) the prices per dose for the two vaccine products.
The computation utilises the epidemiological model[1] to predict the annual number of invasive pneumococcal disease (IPD) under both vaccination programmes and, for comparison, for the scenario 'no vaccination'. The summary table presents the ICERs. The vaccine programme with the lower ICER is identified as the more cost-effective of the two alternatives.
+ 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)
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")
######################
# 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"))
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")
qalysum <- oapply(QALYs, INDEX = QALYs@output[qalyind], FUN = sum)
qalysum@name <- ""
colnames(qalysum@output)[colnames(qalysum@output) == "QALYsResult"] <- "Result"
costsum <- Total_costs
#### The actual model
ICER <- EvalOutput(ICER)
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
ipdtable <- oapply(VacIPD, VacIPD@output["Vaccine"], sum)@output
colnames(ipdtable)[colnames(ipdtable) == "VacIPDResult"] <- "N_of_IPD_cases"
oprint(
ipdtable[order(match(ipdtable$Vaccine, qorder)),],
sortable = FALSE,
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",
digits = rep(0, ncol(ipdtable) + 1)
)
##############################
## print health care costs table
sum_table1A <- data.frame(
Vaccine = qorder,
Medical_costs = result(health_care_costs)[match(qorder,health_care_costs@output$Vaccine)] * 1e-6,
Vaccine_programme_cost = result(vacprice) * 1e-6,
Health_care_costs = result(costsum)[match(qorder,costsum@output$Vaccine)] * 1e-6
)
oprint(
sum_table1A,
sortable = FALSE,
include.rownames = FALSE,
caption = "Table 2. Health care costs (in MEUR)",
caption.placement = "top",
digits = c(0,0,2,2,2)
)
##############################
## 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 = QALYs_gained,
Incremental_effect = QALYs_incremental,
Health_care_costs = Cost_total * 1e-6,
Incremental_cost = Cost_incremental * 1e-6,
ICER = ICER2
)
oprint(
sum_table2,
sortable = FALSE,
include.rownames = FALSE,
caption = "Table 3. Cost-effectiveness analysis summary table ",
caption.placement = "top",
digits = c(0,0,0,0,2,2,6)
)
| |
Variable initiation (Only for developers)
+ Show code- Hide code
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")
| |
Cost calculation (Only for developers)
+ Show code- Hide code
library(OpasnetUtils)
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"
)
| |
Rationale
Vaccination programmes are ranked in ascending order according to their effectiveness. The effectiveness is measured as the expected reduction in invasive pneumococcal disease, as predicted by the epidemiological model.
Alternatives for which there is at least one other alternative with lower cost and better effectiveness are first excluded.
Each programme ('A') is then compared to the next more effective programme ('B') by the incremental cost-effectiveness ratio (ICER):
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ICER = \frac{(C_B-S_B) - (C_A-S_A)}{E_B-E_A},}
where C is the price of the vaccination program, S is the savings in health care costs (as compared to strategy 'no vaccination') and E is the savings in QALYs (as compared to 'no vaccination'). Any programme that is followed by a (more effective) programme with a smaller ICER (i.e. one which produces an additional unit of effect with lower cost) is dropped off from further consideration. The ICERs are then re-calculated and the procedure repeated as many times as needed to eventually identify the most cost-effective alternative. For a tutorial on incremental cost effectiveness analysis, see Phillips (2009) [2].
Costs
Health care resource use in secondary health care, per IPD case and sequelae after meningitis, were estimated from the Hospital Discharge Register (2000-2006). For each meningitis and bacteremia case, an episode of care was constructed by linking the outpatient visits and inpatient hospitalizations, using the unique personal identity code. The case fatality ratio (CFR) for IPD was obtained from a Finnish study [3]. The unit costs for hospitalizations and outpatient visits were estimated based on individual-level cost accounting data from one hospital district. Other unit cost estimates were mainly taken from a widely used national price list for the unit costs of health care in Finland. The costs were presented in 2012 prices and were evaluated from the health care provider perspective. Future costs and benefits were discounted at 3% per annum.
Sensitivity
The effects of alternative vaccine compositions on the outcomes of the cost-benefit analysis were assessed. Five modifications for PCV10 and one for PCV13 were considered Conclusion: The assumption about serotype 3 in PCV13 is crucial. In addition, assumptions about the role of 6A in PCV10 is important. For results, see Cost_effectiveness_sensitivity.
Data
Summary table of the data applied in the cost-effectiveness analysis.
1. QALY_menin = QALY losses due to meningitis (in years, *)
2. QALY_bact = QALY losses due to bacteremia (in years, *)
3. CFR = Case fatality ratio for meningitis and bacteremia
4. Life_y_lost = Life years lost due to IPD (mengitis or bacteremia, *)
5. Cost_ menin = Medical costs attributed to meningitis (in euros *)
6. Cost_ bact = Medical costs attributed to bacteremia (in euros *)
7. Menin_proportion = Proportion of meningitis cases of all IPD cases
(*) a discount rate of 3%/year was applied in all calculations
Estimated medical costs and years lost due to a single bacteremia or meningitis episode
Age class |
QALY_men |
QALY_bac |
CFR |
Life_y_lost |
COST_men |
COST_bac |
Menin_proportion
|
<5 years |
0.22 |
0.0079 |
0.014 |
31.1 |
22 070 |
1 986 |
0.037
|
5-64 years |
0.16 |
0.0079 |
0.112 |
20.7 |
26 488 |
9 000 |
0.046
|
65+ years |
0.08 |
0.0079 |
0.196 |
9.4 |
21 529 |
6 823 |
0.019
|
- Note: The above table lists averages within each age class. Cost-effectiveness analysis is based on age year -specific values.
Estimated medical costs and years lost in Finland without vaccination (per year)
Age group |
QALY_meningitis |
QALY_bacteremia |
Life_years_lost |
Cost_meningitis |
Cost_bacteremia
|
0-4y |
0.83 |
0.75 |
43.64 |
81 591 |
189 444
|
5-64y |
2.89 |
2.90 |
895.01 |
470 949 |
3 308 515
|
65+y |
0.51 |
2.34 |
555.60 |
125 916 |
2 020 437
|
See also
References