Home page

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 costs 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 administration 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

home <- "~/Documents/Modelling_intro/"     # on OS X
home <- "C:\\Documents\\Modelling_intro\\" # on windows
setwd(home)

Package(s)

If BCEA package is not installed yet, you can do so by

if(!"BCEA" %in% installed.packages()){
        install.packages('BCEA')
}

Once installed, you have to load the package(s) at the start of the R file. This is done using

library(BCEA)

Goal and input

The goal is to calculate the cost-effectiveness of a vaccination program for the elderly against flu, taking into account 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(698) 

# population details
targetgroup     <- 2500000                              #people>65 years
life_expectancy <- rnorm(Nsamples,mean=12,sd=2)                 #in years
discount_rate   <- 0.015

# disease burden
force_of_infection      <- rbeta(Nsamples,4,96)     #draw Nsamples randomly from a beta distribution with parameters 4 and 96
proportion_hospitalised <- rbeta(Nsamples,10,90)
casefatality_ratio      <- rbeta(Nsamples,30,970)

# disease-related costs
cost_per_hospitalised_case    <- rnorm(Nsamples,2000,100)  #dollars
cost_per_nonhospitalised_case <- rnorm(Nsamples,100,5)     #dollars

# program details
uptake_program                <- 0.80
vaccine_efficacy              <- rnorm(Nsamples,mean=0.6,sd=(0.7-0.5)/(2*1.96))
vaccine_price_per_dose        <- 10                #dollars
admin_cost_per_dose           <- rnorm(Nsamples,50,4)      #dollars
fixed_program_costs           <- 0             #dollars

# QALY
QALYslost_per_hospitalised_case    <- rnorm(Nsamples,mean=0.018,sd=0.0018)
alpha_qaly_nonhosp                 <- (0.0082^2*(1-0.0082)/0.0018^2)-0.0082
QALYslost_per_nonhospitalised_case <- rbeta(Nsamples,alpha_qaly_nonhosp,alpha_qaly_nonhosp*(1-0.0082)/0.0082)

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 + admin_cost_per_dose
treatment_cost_per_case   <- cost_per_hospitalised_case * proportion_hospitalised + 
                              (1-proportion_hospitalised) * cost_per_nonhospitalised_case
QALYslost_per_case        <- QALYslost_per_hospitalised_case * proportion_hospitalised +        
                               QALYslost_per_nonhospitalised_case * (1-proportion_hospitalised)

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:

### Similar to the previous tutorial: "HTA Flu (basic)"

Calculate the burden of disease for the comparator without the program in place, including the medical costs:

### Similar to the previous tutorial: "HTA Flu (basic)"

Calculate the difference between the program and the comparator:

### Similar to the previous tutorial: "HTA Flu (basic)"

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_costs) / mean(infections_averted)
incr_cost_per_hospitalisation_prevented <- mean(incr_total_costs) / mean(hospitalisations_averted)
incr_cost_per_death_averted             <- mean(incr_total_costs) / mean(deaths_averted)
incr_cost_per_lifeyear_gained           <- mean(incr_total_costs) / mean(lifeyearslost_averted)
incr_cost_per_lifeyear_gained_disc      <- mean(incr_total_costs) / mean(lifeyearslost_averted_disc)
incr_costs_per_QALY_gained              <- mean(incr_total_costs) / mean(QALY_gained)
incr_costs_per_QALY_gained_disc         <- mean(incr_total_costs) / 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           = incr_cost_per_lifeyear_gained,
               incr_cost_per_lifeyear_gained_disc      = incr_cost_per_lifeyear_gained_disc,
               incr_costs_per_QALY_gained              = incr_costs_per_QALY_gained,
               incr_costs_per_QALY_gained_disc         = incr_costs_per_QALY_gained_disc)

# print rounded results
round(ICER_values)
##      incr_cost_per_case_prevented incr_cost_per_hospitalisation_prevented
## [1,]                         2227                                   22212
##      incr_cost_per_death_averted incr_cost_per_lifeyear_gained
## [1,]                       73688                          6160
##      incr_cost_per_lifeyear_gained_disc incr_costs_per_QALY_gained
## [1,]                               8105                       6007
##      incr_costs_per_QALY_gained_disc
## [1,]                            7842

Please 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

Please note that the mean estimate for the incremental total cost per infection averted differs from the value we obtained when not accounting for uncertainty. This discrepancy could be due to extreme values in the ratio. Kindly review the following plots and estimates:

mean(incr_total_costs/infections_averted)       # mean of ratios
## [1] 3023.811
median(incr_total_costs/infections_averted)     # median of ratios
## [1] 2456.958
mean(incr_total_costs)/mean(infections_averted) # ratio of means: ICER!!
## [1] 2226.53

Cost-effectiveness analysis considering uncertainty (using BCEA package).

  1. 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_QALYslost_disc,-program_QALYslost_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_costs   <- cbind(comparator_totaldirectcosts,program_totaldirectcosts)

# 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_costs,
          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'

  1. Which decision should we take, given current information and for a range of WTP values? Print the summary of the bcea() function:
summary(m)
## 
## Cost-effectiveness analysis summary 
## 
## Reference intervention:  Vaccine program
## Comparator intervention: No intervention
## 
## Optimal decision: choose No intervention for k < 7900 and Vaccine program for k >= 7900
## 
## 
## Analysis for willingness to pay parameter k = 25000
## 
##                 Expected net benefit
## No intervention           -848207898
## Vaccine program           -615870901
## 
##                                          EIB   CEAC   ICER
## Vaccine program vs No intervention 232336997 0.9358 7841.8
## 
## Optimal intervention (max expected net benefit) for k = 25000: Vaccine program
##             
## EVPI 1824648
  1. How much evidence is there in favour of the decision, for a range of WTP values? Look also at the CEAC above.
m_ceaf <- multi.ce(m)
ceaf.plot(m_ceaf) # uses base plot

  1. Is there value in further research? Look at the EVPI graph and summary above.