Cost-effectiveness analysis for influenza vaccination with uncertainty
This model was developed by Philippe Beutels, Joke Bilcke and Lander Willem for illustration and teaching purposes only. Please make contact if you would like to use this model in any public way.
In the field of economic evaluation in general, it could be applicable to any intervention with effectiveness realised within one year (like many curative interventions). In the field of infectious disease, it would only be applicable to influenza, and only if vaccination targets a small proportion of the population that does not at the same time form a core transmitter group of the virus (e.g. it is unlikely to be suitable to model childhood influenza vaccination). Note that for most applications, practical use of this model would likely entail having separate, though similarly simple decision trees to obtain the unit cost estimates (implying that these cells would be intermediary outcomes, and that there would be additional input cells for each relevant stage of disease (for instance for the number of consultations, the unit cost per consultation, medication use, hospital days etc. for cases of pneumonia, and all other relevant disease stages)).
Note that the final output is independent of coverage and target group size, if fixed delivery costs are set to 0. This is due to the exclusion of herd immunity effects in this model.
For further information, questions or comments please contact us.
Set working directory
The first step when creating a new R file is to specify a working
directory. You need to tell the R session where you want it to work. To
do this you use the command to “set the working directory” or
setwd().
You can label this working directory as ‘home’, and have other locations stored for plotting or for getting data. This location should be an existing folder on your computer. For example
Specific uncertainty distributions
For some parameters, the sample size is rather small or they should be strictly positive, hence other distribution might be preferred. Let’s first define some help functions to facilitate sampling using a given mean, se, proportion, samples size or confidence interval:
sample_gamma = function(nsamples, mean, se){
rgamma(nsamples, shape = mean^2 / se^2, rate = mean / se^2)
}
sample_beta_event = function(nsamples, events, sample_size){
return(rbeta(nsamples, shape1 = events, shape2 = sample_size - events))
}
sample_beta_mean = function(nsamples, mean, se){
alpha = mean^2 * (1 - mean) / se^2
beta = alpha * (1 - mean) / mean
return(rbeta(nsamples, shape1 = alpha, shape2 = beta))
}Goal and input
The goal is to calculate the cost-effectiveness of a vaccination program for the elderly against flu, taking into account disease burden, discounting, QALYs and parameter uncertainty as specified below:
# set number of samples to include in the uncertainty analysis
nsamples = 5000
# set the random-number-generator seed to get exactly the same results each time you run the code
set.seed(5698)
# population details
targetgroup = 2500000 # people>65 years
life_expectancy = sample_gamma(nsamples, 12, 2) # in years
discount_rate = 0.015
# disease burden
force_of_infection = sample_beta_event(nsamples, 4, 100)
p_hosp = sample_beta_event(nsamples, 10, 100)
hosp_case_fatality_ratio = sample_beta_event(nsamples, 30, 100)
# disease-related costs
cost_per_hosp_case = sample_gamma(nsamples, 2000, 100) # euro
cost_per_non_hosp_case = sample_gamma(nsamples, 100, 5) # euro
# program details
uptake_program = 0.80
vaccine_efficacy = sample_beta_mean(nsamples, mean = 0.6, se = (0.7 - 0.5) / (2 * 1.96))
vaccine_price_per_dose = 10 # euro
delivery_cost_per_dose = sample_gamma(nsamples, 50, 4) # euro
fixed_program_cost = 0 # euro
# QALY
QALYloss_hosp_case = sample_beta_mean(nsamples, 0.018, 0.0018)
QALYloss_non_hosp_case = sample_beta_mean(nsamples, 0.0082, 0.0018)Health Technology Assessment: step-by-step
Now, propagate the uncertainty from the input parameters of the flu model into the outcome of your flu model. To start the health technology assessment, we calculate the vaccine cost per dose and the treatment cost per case, taking into account the proportion of cases that lead to hospitalization, the QALY lost per case, and the discounted life expectancy:
vaccine_cost_per_dose = vaccine_price_per_dose + delivery_cost_per_dose
treatment_cost_per_case = cost_per_hosp_case * p_hosp + (1 - p_hosp) * cost_per_non_hosp_case
QALYloss_per_case = QALYloss_hosp_case * p_hosp + QALYloss_non_hosp_case * (1 - p_hosp)
life_expectancy_disc = rep(NA, nsamples)
for (i in 1:nsamples){
life_expectancy_disc[i] = sum(1 / ((1 + discount_rate)^(0:(life_expectancy[i] - 1))))
}Calculate the burden of disease with the program in place, along with the associated medical and program costs:
Calculate the burden of disease for the comparator without the program in place, including the medical costs:
Calculate the difference between the program and the comparator:
Calculate the incremental effectiveness ratios for infections, hospitalizations, and deaths averted. Warning, the expected ICER equals the expected incremental cost/average incremental effect. this is NOT the same as the average all ICER\(_i\) values!! Finally, calculate the incremental effectiveness ratios for life years saved:
# Calculate incremental effectiveness ratios
incr_cost_per_case_prevented = mean(incr_total_cost) / mean(infections_averted)
incr_cost_per_hospitalisation_prevented = mean(incr_total_cost) / mean(hospitalisations_averted)
incr_cost_per_death_averted = mean(incr_total_cost) / mean(deaths_averted)
incr_cost_per_lifeyear_gained_disc = mean(incr_total_cost) / mean(lifeyearslost_averted_disc)
incr_cost_per_QALY_gained_disc = mean(incr_total_cost) / mean(QALY_gained_disc)Show the results:
# combine results
# note: specify column names at once
ICER_values = cbind(incr_cost_per_case_prevented = incr_cost_per_case_prevented,
incr_cost_per_hospitalisation_prevented = incr_cost_per_hospitalisation_prevented,
incr_cost_per_death_averted = incr_cost_per_death_averted,
incr_cost_per_lifeyear_gained_disc = incr_cost_per_lifeyear_gained_disc,
incr_cost_per_QALY_gained_disc = incr_cost_per_QALY_gained_disc)
# print rounded results
round(ICER_values)## incr_cost_per_case_prevented incr_cost_per_hospitalisation_prevented
## [1,] 2212 22120
## incr_cost_per_death_averted incr_cost_per_lifeyear_gained_disc
## [1,] 73767 6957
## incr_cost_per_QALY_gained_disc
## [1,] 6761Please note that there may be slight variations due to platform and R-version specific differences in random number generation.
Second thoughts on the distributions for costs and cases
Kindly review the following plots and estimates:
## [1] 10893.83## [1] 6761.249# Try WTP values bracketing the ICER to find where the expected INMB changes sign.
# At the ICER, the expected INMB is exactly 0. Values above the ICER yield positive
# expected INMB (intervention is cost-effective); values below yield negative INMB.
k = 8000 # above the ICER: expected INMB should be positive
mean(QALY_gained_disc * k - incr_total_cost)## [1] 19431759# the INMB is not affected by the order of the operations
mean(QALY_gained_disc) * k - mean(incr_total_cost)## [1] 19431759k = 6000 # below the ICER: expected INMB should be negative
mean(QALY_gained_disc * k - incr_total_cost)## [1] -11941391Cost-effectiveness plane
# plot CE plane
plot(x = QALY_gained_disc,
y = incr_total_cost,
xlab = 'Incremental effect (QALY gain)',
ylab = 'Incremental cost (euro)',
main = 'CE plane')
Present CEA (w/o R package).
# set WTP range, and step size
opt_wtp = seq(0, 40000, 500)
# initialise matrix
inmb_wtp_summary = data.frame(wtp = opt_wtp,
inmb = NA,
ceac = NA,
ceaf = NA,
inmb_perfect_info = NA)
# loop over all WTP levels
for(i in 1:length(opt_wtp)){
# calculate the inmb per sample with the current WTP threshold
inmb_sample = (QALY_gained_disc * opt_wtp[i]) - incr_total_cost
# get the expected INMB
inmb_wtp_summary$inmb[i] = mean(inmb_sample)
# inspect the proportion of samples with a positive INMB
inmb_wtp_summary$ceac[i] = sum(inmb_sample >= 0) / nsamples
# to calculate INMB with perfect information, set all non-preferred options to zero
# note: if there are more than 2 options, other code is needed
if(inmb_wtp_summary$inmb[i] >= 0){
inmb_sample[inmb_sample < 0] = 0
} else{
inmb_sample[inmb_sample > 0] = 0
}
# get the expected INMB with perfect information
inmb_wtp_summary$inmb_perfect_info[i] = mean(inmb_sample)
}
# get preferred strategy by wtp, in this case the intervention when INMB is positive
bool_inmb_pos = inmb_wtp_summary$inmb > 0
# a boolean (false/true) + 1, gives id 1 and 2
inmb_wtp_summary$ceaf_intervention = bool_inmb_pos + 1
inmb_wtp_summary$ceaf_proportion[bool_inmb_pos] = inmb_wtp_summary$ceac[bool_inmb_pos]
inmb_wtp_summary$ceaf_proportion[!bool_inmb_pos] = 1 - inmb_wtp_summary$ceac[!bool_inmb_pos]
# net-loss: difference between INMB under perfect information and expected INMB.
# This equals the per-patient EVPI: the maximum a decision maker should pay to
# eliminate all parameter uncertainty before choosing between strategies.
inmb_wtp_summary$net_loss = abs(inmb_wtp_summary$inmb_perfect_info - inmb_wtp_summary$inmb)Inspect the CEA outcomes if the WTP threshold is 10000.
## wtp inmb ceac ceaf inmb_perfect_info ceaf_intervention
## 21 10000 50804909 0.6192 NA 66210741 2
## ceaf_proportion net_loss
## 21 0.6192 15405832Inspect the included WTP for which the intervention becomes cost-effective. Compare this with the ICER.
# summary statistics
ce_threshold = min(inmb_wtp_summary$wtp[inmb_wtp_summary$inmb > 0])
icer = mean(incr_total_cost) / mean(QALY_gained_disc)
data.frame(ce_threshold = ce_threshold,
icer = icer)## ce_threshold icer
## 1 7000 6761.249Cost-effectiveness acceptability curve (CEAC)
# CEAC
plot(x = inmb_wtp_summary$wtp,
y = inmb_wtp_summary$ceac,
type = 'l',
lwd = 2,
col = 2,
xlab = 'Willingness To Pay (euro)',
ylab = 'Probability to be cost-effective',
main = 'Cost-Effectiveness Acceptability Curve',
ylim = 0:1)
lines(x = inmb_wtp_summary$wtp,
y = 1 - inmb_wtp_summary$ceac,
lwd = 2,
col = 1)
legend('right',
c('current state','vaccination program'),
col = 1:2,
lwd = 2)
Cost-effectiveness acceptability frontier (CEAF)
plot(x = inmb_wtp_summary$wtp,
y = inmb_wtp_summary$ceaf_proportion,
col = inmb_wtp_summary$ceaf_intervention,
type = 'p',
lwd = 2,
xlab = 'Willingness to pay (euro)',
ylab = 'Probability to be cost-effective',
main = 'Cost-Effectiveness Acceptability Frontier',
ylim = 0:1,
xlim = range(opt_wtp))
legend('right',
c('current state','vaccination program'),
col = 1:2,
lwd = 2)
Or you can combine the CEAC and CEAF
# Combination of CEAC & CEAF
plot(x = inmb_wtp_summary$wtp,
y = inmb_wtp_summary$ceac,
type = 'l',
lwd = 2,
col = 2,
xlab = 'Willingness to pay (euro)',
ylab = 'Probability to be cost-effective',
main = '',
ylim = 0:1)
lines(x = inmb_wtp_summary$wtp,
y = 1 - inmb_wtp_summary$ceac,
type = 'l',
col = 1,
lwd = 2)
points(x = inmb_wtp_summary$wtp,
y = inmb_wtp_summary$ceaf_proportion,
col = inmb_wtp_summary$ceaf_intervention)
legend('right',
c('current state', 'vaccination program', 'CEAF'),
col = c(1, 2, 1),
lwd = c(2, 2, NA),
pch = c(NA, NA, 1))
EVPI
plot(x = inmb_wtp_summary$wtp,
y = inmb_wtp_summary$net_loss,
type = 'b',
xlab = 'Willingness to pay (euro)',
ylab = 'Expected net-loss (euro)',
main = "Expected Value of Perfect Information")
EVPPI
If we are interested in the EVPPI, hence per individual parameter of the influenza model above, we can adopt the GAM approximation using the calculated incremental cost and benefit with all sampled parameter values.
# load package for GAM
library(mgcv)
# incremental QALY gain and costs
# QALY_gained_disc
# incr_total_cost
# combine input parameters with uncertainty
model_param = data.frame(force_of_infection = force_of_infection,
p_hosp = p_hosp,
life_expectancy = life_expectancy,
hosp_case_fatality_ratio = hosp_case_fatality_ratio,
cost_per_hosp_case = cost_per_hosp_case,
cost_per_non_hosp_case = cost_per_non_hosp_case,
vaccine_efficacy = vaccine_efficacy,
delivery_cost_per_dose = delivery_cost_per_dose,
QALYloss_hosp_case = QALYloss_hosp_case,
QALYloss_non_hosp_case = QALYloss_non_hosp_case)
# specify the number of WTP and parameter values
num_param = ncol(model_param)
num_wtp = length(opt_wtp)
# initiate matrix to store EVPPI values
evppi = matrix(0, num_wtp, num_param)
# run over all WTP levels
for(j in 1:length(opt_wtp)){
# set WTP threshold
k = opt_wtp[j]
# calculate INMB
inmb_sample = (QALY_gained_disc * k) - incr_total_cost
# approximate NL for each parameter value
for(i in 1:num_param){
model = gam(inmb_sample ~ model_param[, i])
g.hat_new = model$fitted
perfect.info = mean(pmax(g.hat_new, 0))
baseline = pmax(mean(g.hat_new), 0)
## estimate EVPPI
evppi[j, i] = round(perfect.info - baseline, digits = 4)
}
}
# add column names
colnames(evppi) = colnames(model_param)Explore EVPPI
# if all EVPPI values are 0, exclude parameter
evppi_edit = evppi[, colSums(evppi, na.rm = T) > 0]
# start with empty plot, and add EVPPI one by one
plot(x = range(opt_wtp),
y = range(evppi_edit),
col = 0,
xlab = "WTP",
ylab = "EVPPI (euro)")
for(i in 1:ncol(evppi_edit)){
lines(x = opt_wtp,
y = evppi_edit[, i],
type = 'b',
lwd = 2,
col = i,
pch = i)
}
legend('topright',
colnames(model_param),
col = 1:num_param,
pch = 1:num_param,
lwd = 2,
cex = 0.7)
Present CEA (using the BCEA package).
If BCEA package is not installed yet, you can do so
by
Once installed, you have to load the package(s) using
- Run the
bcea()function to obtain the cost-effectiveness plane, the expected incremental net monetary benefit (INMB), the cost-effectiveness acceptability curve (CEAC) and the expected value for perfect information (EVPI):
# inspect the documentation of the 'bcea()' function: ?bcea
# Specify a matrix containing the clinical effectiveness for each intervention being considered.
# note: this needs to be specified as health gain (i.e. negative health loss)
ce_effects = cbind(-comparator_QALYloss_disc, -program_QALYloss_disc)
# Specify a matrix containing the cost for each intervention being considered.
# note: this needs to be specified as health gain (i.e. negative health loss)
ce_cost = cbind(comparator_total_cost, program_total_cost)
# Define the labels to be associated with each intervention.
interventions = c("No intervention", "Vaccine program")
# Run the 'bcea()' function
m = bcea(eff = ce_effects,
cost = ce_cost,
ref = 2, # the column of eff and cost considered to be the reference strategy.
interventions = interventions,
plot = T) # ref specifies to compare 'Vaccine program' to 'No intervention'
- Which decision should we take, given current information and for a
range of WTP values? Print the summary of the
bcea()function:
##
## Cost-effectiveness analysis summary
##
## Reference intervention: Vaccine program
## Comparator intervention: No intervention
##
## Optimal decision: choose No intervention for k < 6800 and Vaccine program for k >= 6800
##
##
## Analysis for willingness to pay parameter k = 10000
##
## Expected net benefit
## No intervention -356233083
## Vaccine program -305428174
##
## EIB CEAC ICER
## Vaccine program vs No intervention 50804909 0.6192 6761.2
##
## Optimal intervention (max expected net benefit) for k = 10000: Vaccine program
##
## EVPI 15405832- How much evidence is there in favour of the decision, for a range of WTP values? Look at the CEAC and CEAF.
- Is there value in further research? Look at the EVPI graph and summary above.
EVPPI
The BCEA package allows to calculate the expected value of partially perfect information for parameter subsets.
# combine input parameters with uncertainty
model_param = data.frame(force_of_infection = force_of_infection,
QALYloss_non_hosp_case = QALYloss_non_hosp_case)
inp = createInputs(model_param)
# explicitly use BCEA package namespace to avoid package conflict
evppi_bcea = BCEA::evppi(m, inp$parameters, inp$mat)
plot(evppi_bcea, col = 1:2, pos = c(1, 1))