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 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)

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 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) 
case_fatality_ratio      = sample_beta_event(nsamples,30,1000) 

# 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
admin_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 + admin_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:

### Similar to the previous tutorial: "CEA: Influenza"

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

### Similar to the previous tutorial: "CEA: Influenza"

Calculate the difference between the program and the comparator:

### Similar to the previous tutorial: "CEA: Influenza"

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           = mean(incr_total_cost) / mean(lifeyearslost_averted)
incr_cost_per_lifeyear_gained_disc      = mean(incr_total_cost) / mean(lifeyearslost_averted_disc)
incr_cost_per_QALY_gained               = mean(incr_total_cost) / mean(QALY_gained)
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           = incr_cost_per_lifeyear_gained,
               incr_cost_per_lifeyear_gained_disc      = incr_cost_per_lifeyear_gained_disc,
               incr_cost_per_QALY_gained               = incr_cost_per_QALY_gained,
               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,]                         2210                                   22110
##      incr_cost_per_death_averted incr_cost_per_lifeyear_gained
## [1,]                       73784                          6157
##      incr_cost_per_lifeyear_gained_disc incr_cost_per_QALY_gained
## [1,]                               8098                      6004
##      incr_cost_per_QALY_gained_disc
## [1,]                           7834

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_cost/QALY_gained_disc)       # mean of ratios

## [1] 11650.1

median(incr_total_cost/QALY_gained_disc)     # median of ratios

## [1] 8858.122

mean(incr_total_cost)/mean(QALY_gained_disc) # ratio of means: ICER!!

## [1] 7834.368

Cost-effectiveness plane

# plot CE plane
plot(QALY_gained_disc, incr_total_cost, xlab = 'Incremental effect', 
     ylab = 'Incremental cost', main = 'CE plane')

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

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:

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           -846667222
## Vaccine program           -614317689
## 
##                                          EIB   CEAC   ICER
## Vaccine program vs No intervention 232349533 0.9334 7834.4
## 
## Optimal intervention (max expected net benefit) for k = 25000: Vaccine program
##             
## EVPI 2104072

How much evidence is there in favour of the decision, for a range of WTP values? Look at the CEAC and CEAF.

m_ceaf = multi.ce(m)
ceaf.plot(m_ceaf) # uses base plot

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))

However, if are now interested in the EVPPI 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

# specify WTP
opt_wtp = seq(0,50000,500)

# combine input parameters with uncertainty
model_param = data.frame(life_expectancy = life_expectancy,
                         force_of_infection= force_of_infection,
                         p_hosp= p_hosp,
                         case_fatality_ratio= 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,
                         admin_cost_per_dose= admin_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(range(opt_wtp),range(evppi_edit),col=0, xlab = "WTP",ylab="EVPPI (euro)")
for(i in 1:ncol(evppi_edit)){
  lines(opt_wtp,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)

Compare with EVPPI from BCEA

plot(evppi_bcea,col=1:2,pos = c(1,1))
lines(opt_wtp, evppi_edit[,2],col=3)
legend('right',
       c('manual: EVPPI FoI'),
       col=3,
       lwd=1,
       cex = 0.7)

Navigate

Top: Index

Cost-effectiveness analysis for influenza vaccination with uncertainty