This tutorial has been developed by 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. For further information, questions or comments please contact: joke.bilcke@uantwerp.be or lander.willem@uantwerp.be from the Centre for Health Economics Research and Modeling of Infectious Diseases (CHERMID), Vaccine and Infectious Disease Institute (Vaxinfectio), University of Antwerp, Universiteitsplein 1, 2610 Antwerp, Belgium.

Research question

In this session, we delve into a critical question: Is vaccinating the elderly against Respiratory Syncytial Virus (RSV) a cost-effective strategy when compared to not vaccinating them? RSV is known to cause lower respiratory tract infections in infants, the elderly, and individuals with underlying health conditions.

Amidst the backdrop of numerous emerging interventions aimed at preventing RSV, our focus today is to conduct a comprehensive evaluation of the cost-effectiveness of implementing a vaccination program for the elderly. To achieve this, we will employ both static and dynamic disease transmission models, allowing us to gain valuable insights into the potential impact of vaccination strategies.

While the ICER serves as a crucial measure of cost-effectiveness, we’ll also encourage you to extend your analysis further. Consider calculating the Incremental Net Monetary Benefit (INMB) and the Net Loss (NL) associated with the vaccination program. This broader perspective will enhance your grasp of the economic implications and provide a more comprehensive view of the intervention’s impact.

Objectives of the Tutorial:

Evaluate the cost-effectiveness of an elderly vaccination program against RSV in comparison to a scenario without such a program.
Utilize disease transmission models to analyze the dynamics of RSV spread and its interactions within the population.
Obtain and interpret different measures of cost-effectiveness: the Incremental Cost-Effectiveness Ratio (ICER), quantified as the incremental cost per QALY gained, the Incremental Net Monetary Benefit (INMB) and the Net Loss (NL).

As we delve into the topic, we encourage you to explore various scenarios by manipulating input parameters. By doing so, you can gain a nuanced understanding of how altering these factors influences the cost-effectiveness outcome measures. This hands-on exploration will empower you to learn about the sensitivity of our analysis to various input parameters and its implications for decision-making.

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)

MODEL INPUT

The disease incidence and economic parameters are loosely based on the work of Zeevat et al (2022), published in the Journal of Infectious Diseases.

# population and disease input
population_size           = 4557901 # NL population +60
rsv_incidence             = 0.03525 # new infections per year, in the absence of vaccination
rsv_recovery_time         = 7 # days

# health economic input
cost_per_episode          = 122 # euro
morbidity_per_episode     = 0.002196 # QALY loss
proportion_mortality      = 0.006    # assume age-invariant
average_life_expectancy   = 9        # years

# vaccine characteristics
vaccine_efficacy  = 0.60
cost_per_dose     = 10 # euro 
admin_cost        = 50 # euro

# program settings
uptake = 0.2 # proportion of the population
time_horizon = 365 # days

PRE-PROCESSING

Based on the input above, we can compute the program cost and the proportion of the population that will be effectively protected:

program_costs = population_size * uptake * (cost_per_dose + admin_cost)
proportion_protected =  uptake * vaccine_efficacy
qaly_loss_per_episode = morbidity_per_episode + (proportion_mortality * average_life_expectancy)

STATIC DISEASE MODEL

We commence by delving into the static disease model. The following code delineates the model’s intrinsic parameters and conducts ICER computation. The outcomes revolve around the baseline disease impact, portraying anticipated infections in the absence of a vaccination program. Moreover, we explore averted infections in the case of a vaccination program, subsequently estimating the associated QALY improvement. This culmination leads us to one of the measures of cost-effectiveness: the ICER.

# TRANSLATE TO STATIC MODEL SETTING
force_of_infection  = rsv_incidence

# INFECTIONS (prevented)
reference_infections = population_size * force_of_infection
program_infections = population_size * (1 - proportion_protected) * force_of_infection
prevented_cases = reference_infections - program_infections

# COSTS AND EFFECTS
static_incr_cost = program_costs - (prevented_cases * cost_per_episode)
static_qaly_gain = prevented_cases * qaly_loss_per_episode

# INCREMENTAL COST EFFECTIVENESS RATIO
static_icer = static_incr_cost  / static_qaly_gain

# PRINT RESULTS
round(c(reference_infections = reference_infections,
        program_infections = program_infections,
        static_qaly_gain = static_qaly_gain,
        static_icer = static_icer),digits=2)

## reference_infections   program_infections     static_qaly_gain 
##            160666.01            141386.09              1083.45 
##          static_icer 
##             48310.90

DYNAMIC DISEASE MODEL: general

Moving forward, our exploration advances to encompass the dynamic transmission model, while adhering to the identical model parameters and fundamental disease attributes as established in the static model.

For this example, we utilize the EpiDynamics package, which features an array of built-in functions tailored for such analyses. A noteworthy example is the incorporation of the ‘SIR’ function, adeptly employed for quantifying the disease burden. Through this approach, we effectively calculate the disease’s impact while considering distinct initial conditions and transmission parameters:

# check if the "EpiDynamics" package is installed, if not, install package
if(!'EpiDynamics' %in% installed.packages()){ install.packages('EpiDynamics') }

# load library
library(EpiDynamics)

# Parameters and initial conditions. 
parameters = list(beta = 5.147, gamma = 1 / rsv_recovery_time)
initials = c(S = 1 - 1e-06 - 0.2, I = 1e-06, R = 0.2) 

# Solve and plot. 
sir = SIR(pars = parameters, init = initials, time = 0:time_horizon) 

# Plot output
PlotMods(sir)

To ascertain the quantity of infections, we can deduce this by evaluating the alteration in the “Recovered” compartment:

sir$results$R[nrow(sir$results)] - sir$results$R[1]

## [1] 0.8

DYNAMIC DISEASE MODEL: RSV

In order to achieve the desired RSV incidence for this Cost-Effectiveness Analysis, we can posit a considerable existing level of immunity within the population:

# set initial immunity levels
initials['R'] = 0.95099
initials['S'] = 1 - sum(initials['I']+initials['R'])

# Solve 
dynamic_rsv_reference = SIR(pars = parameters,init = initials,time = 0:time_horizon)

# calculate total infections
dynamic_cases_reference = diff(dynamic_rsv_reference$results$R[c(1,nrow(dynamic_rsv_reference$results))]) * population_size
dynamic_cases_reference # print to terminal

## [1] 160638.9

To incorporate vaccination within the model, we can conceptualize the adoption of the vaccine as a transition from the “Susceptible” to the “Recovered” compartment:

initials['S'] = initials['S'] * (1-proportion_protected)
initials['R'] = 1 - sum(initials['S']+initials['I'])

# Solve 
dynamic_rsv = SIR(pars = parameters,init = initials,0:time_horizon)

#PlotMods(dynamic_rsv)
dynamic_cases_program = diff(dynamic_rsv$results$R[c(1,nrow(dynamic_rsv$results))]) * population_size

# prevented cases
dynamic_cases_prevented = dynamic_cases_reference - dynamic_cases_program

# CEA
dynamic_incr_cost = program_costs - (dynamic_cases_prevented * cost_per_episode)
dynamic_qaly_gain = dynamic_cases_prevented * qaly_loss_per_episode

# ICER
dynamic_icer = dynamic_incr_cost  / dynamic_qaly_gain

# PRINT RESULTS
round(c(dynamic_cases_reference=dynamic_cases_reference,
        dynamic_cases_program=dynamic_cases_program,
        dynamic_qaly_gain=dynamic_qaly_gain,
        dynamic_icer=dynamic_icer),digits=2)

## dynamic_cases_reference   dynamic_cases_program       dynamic_qaly_gain 
##               160638.94               121111.63                 2221.28 
##            dynamic_icer 
##                22452.17

Inspection of the disease prevalence:

SUMMARY

Inspect the cost-effectiveness plane:

plot(x = c(static_qaly_gain,dynamic_qaly_gain),
     y = c(static_incr_cost,dynamic_incr_cost),
     col = 1:2,
     xlab= 'QALY gain',
     ylab = 'incremental cost (euro)',
     pch=16)
legend('topright',
       c('static model',
         'dynamic model'),
       fill = 1:2,
       )

model	reference_infections	program_infections	qaly_gain	incr_cost	program_icer	infections_reduction
static	160666	141386	1083	52342662	48311	0.12
dynamic	160639	121112	2221	49872480	22452	0.25

Incremental Net Monetary Benefit (INMB)

Compute the Incremental Net Monetary Benefit for the RSV illustration provided, while adhering to willingness-to-pay thresholds of 20,000 and 40,0000 euro per Quality-Adjusted Life Year (QALY).

# define wtp thresholds
opt_wtp <- c(20e3,40e3)

# define function to calculate the inmb
get_inmb <- function(wtp,incr_qaly,incr_cost){
  return(incr_qaly*wtp-incr_cost)
}

# static model
static_inmb <- get_inmb(opt_wtp,static_qaly_gain,static_incr_cost)

# dynamic model
dynamic_inmb <- dynamic_inmb <- get_inmb(opt_wtp,dynamic_qaly_gain,dynamic_incr_cost)

# aggregate into table
db_out_inmb <- db_out[c(1,1,2,2),c(1,4,5,6)]
db_out_inmb$wtp <- rep(opt_wtp,2)
db_out_inmb$inmb <- c(static_inmb,dynamic_inmb)

kable(db_out_inmb)

	model	qaly_gain	incr_cost	program_icer	wtp	inmb
1	static	1083	52342662	48311	20000	-30673573
1.1	static	1083	52342662	48311	40000	-9004483
2	dynamic	2221	49872480	22452	20000	-5446948
2.1	dynamic	2221	49872480	22452	40000	38978585

Exploring model behaviour and sensitivity

# To conduct model explorations, it is convenient to create a function for 
# the static and dynamic approach

run_static_model <- function(uptake = 0.2,
                             admin_cost = 50,
                             wtp = 20e3,
                             rsv_incidence = 0.03525){
  # PRE-PROCESSING
  program_costs = population_size * uptake * (cost_per_dose + admin_cost)
  proportion_protected =  uptake * vaccine_efficacy

 # TRANSLATE TO STATIC MODEL SETTING
  force_of_infection  = rsv_incidence
  
  # INFECTIONS (prevented)
  reference_infections = population_size * force_of_infection
  program_infections   = population_size * (1 - proportion_protected) * force_of_infection
  prevented_cases      = reference_infections - program_infections
  
  # COSTS AND EFFECTS
  static_incr_cost = program_costs - (prevented_cases * cost_per_episode)
  static_qaly_gain = prevented_cases * qaly_loss_per_episode
  
  # INCREMENTAL COST EFFECTIVENESS RATIO
  static_icer = static_incr_cost  / static_qaly_gain
 
  # INMB
  static_inmb = wtp * static_qaly_gain - static_incr_cost
  
  # RETURN RESULTS
  return(c(
           uptake = uptake,
           round(c(
                 admin_cost = admin_cost,
                 wtp = wtp,
                 reference_infections = reference_infections,
                 program_infections = program_infections,
                 qaly_gain = static_qaly_gain,
                 incr_cost = static_incr_cost,
                 icer = static_icer,
                 inmb_k=static_inmb/1e3),digits=0)))
}

run_static_model(uptake = 0.1)

##               uptake           admin_cost                  wtp 
##                  0.1                 50.0              20000.0 
## reference_infections   program_infections            qaly_gain 
##             160666.0             151026.0                542.0 
##            incr_cost                 icer               inmb_k 
##           26171331.0              48311.0             -15337.0

run_dynamic_model <- function(uptake = 0.2,
                              admin_cost = 50,
                              wtp = 20e3,
                              initial_R = 0.95099){
  
  # pre-processing
  program_costs = population_size * uptake * (cost_per_dose + admin_cost)
  proportion_protected =  uptake * vaccine_efficacy

  # Reference
  initials['I'] = 1e-06
  initials['R'] = initial_R
  initials['S'] = 1 - sum(initials['I']+initials['R'])
  dynamic_rsv = SIR(pars = parameters,init = initials,0:time_horizon)
  dynamic_cases_reference = diff(dynamic_rsv$results$R[c(1,nrow(dynamic_rsv$results))]) * population_size
  
  # Program
  initials['S'] = initials['S'] * (1-proportion_protected)
  initials['R'] = 1 - sum(initials['S']+initials['I'])
  dynamic_rsv = SIR(pars = parameters,init = initials,0:time_horizon)
  dynamic_cases_program = diff(dynamic_rsv$results$R[c(1,nrow(dynamic_rsv$results))]) * population_size
  
  # prevented cases
  dynamic_cases_prevented = dynamic_cases_reference - dynamic_cases_program
  
  # CEA
  dynamic_incr_cost = program_costs - (dynamic_cases_prevented * cost_per_episode)
  dynamic_qaly_gain = dynamic_cases_prevented * qaly_loss_per_episode
  
  # ICER
  dynamic_icer = dynamic_incr_cost  / dynamic_qaly_gain
  
  # INMB
  dynamic_inmb = wtp * dynamic_qaly_gain - dynamic_incr_cost
  
  # RETURN RESULTS
  return(c(
           uptake = uptake,
           round(c(
                 admin_cost = admin_cost,
                 wtp = wtp,
                 reference_infections = dynamic_cases_reference,
                 program_infections = dynamic_cases_program,
                 qaly_gain = dynamic_qaly_gain,
                 incr_cost = dynamic_incr_cost,
                 icer = dynamic_icer,
                 inmb_k=dynamic_inmb/1e3),digits=0)))
}

Impact of Vaccine Uptake: Evaluate the effect of altering vaccine uptake to 0.1, 0.2, 0.7 or 0.9 on ICERs and INMB.

static_out <- rbind(run_static_model(uptake = 0.1),
                    run_static_model(uptake = 0.2),
                    run_static_model(uptake = 0.7),
                    run_static_model(uptake = 0.9))

dynamic_out <-rbind(run_dynamic_model(uptake = 0.1),
                    run_dynamic_model(uptake = 0.2),
                    run_dynamic_model(uptake = 0.7),
                    run_dynamic_model(uptake = 0.9))

kable(cbind(model='static',static_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
static	0.1	50	20000	160666	151026	542	26171331	48311	-15337
static	0.2	50	20000	160666	141386	1083	52342662	48311	-30674
static	0.7	50	20000	160666	93186	3792	183199316	48311	-107358
static	0.9	50	20000	160666	73906	4876	235541977	48311	-138031

kable(cbind(model='dynamic',dynamic_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
dynamic	0.1	50	20000	160639	141211	1092	24977179	22877	-3142
dynamic	0.2	50	20000	160639	121112	2221	49872480	22452	-5447
dynamic	0.7	50	20000	160639	461	9001	171890164	19096	8137
dynamic	0.9	50	20000	160639	24	9026	226531664	25098	-46014

# Dynamic disease models enable herd immunity effects. As such, the 
# number of prevented cases can outnumber the (effectively) protected 
# individuals. This reduces the ICER, since the health related gains 
# increase for the same investment. 

# When the uptake is already high, additional uptake can become less 
# rewarding.

# note on INMB:
inmb_example1 <- c(wtp = 20e3,
                   incr_cost = 4e3,
                   qaly_gain = 2,
                   icer = 4e3 / 2,
                   inmb = 20e3*2 - 4e3)
inmb_example2 <- c(wtp = 20e3,
                   incr_cost = 4e3*2,
                   qaly_gain = 2*2,
                   icer = (4e3*2) / (2*2) ,
                   inmb = 20e3*(2*2) - (4e3*2))
inmb_example3 <- c(wtp = 20e3,
                   incr_cost = 4e3*6,
                   qaly_gain = 2*6,
                   icer = (4e3*6) / (2*6) ,
                   inmb = 20e3*(2*6) - (4e3*6))
kable(rbind(inmb_example1,inmb_example2,inmb_example3))

	wtp	incr_cost	qaly_gain	icer	inmb
inmb_example1	20000	4000	2	2000	36000
inmb_example2	20000	8000	4	2000	72000
inmb_example3	20000	24000	12	2000	216000

Consistency Across Models: Examine whether changes in vaccine uptake lead to similar outcomes in both disease models. Provide an explanation for the observed results.

# Static disease models do not take herd immunity effects into account. Only 
# personal protection is considered in the economic evaluation. This holds 
# when the target group has almost no effect on the transmission dynamics. 
# The ICER does not change when the vaccine uptake varies.

Combining Vaccines: Investigate the consequences of combining RSV and Influenza vaccine uptake on ICERs. Consider the influence of admission costs in your analysis.

static_out <- rbind(run_static_model(admin_cost = admin_cost),
                    run_static_model(admin_cost = admin_cost/2))

dynamic_out <-rbind(run_dynamic_model(admin_cost = admin_cost),
                    run_dynamic_model(admin_cost = admin_cost/2))

kable(cbind(model='static',static_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
static	0.2	50	20000	160666	141386	1083	52342662	48311	-30674
static	0.2	25	20000	160666	141386	1083	29553157	27277	-7884

kable(cbind(model='dynamic',dynamic_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
dynamic	0.2	50	20000	160639	121112	2221	49872480	22452	-5447
dynamic	0.2	25	20000	160639	121112	2221	27082975	12193	17343

# Combining the administration of different vaccines, hence reducing the 
# administration cost, is beneficial for the cost-effectiveness.

Matching Static and Dynamic Models: Identify the parameter settings within the dynamic model that would yield results consistent with the static model. Discuss the significance of these settings.

# If the (basic) reproduction number is very high, vaccine-related immunization 
# has limited impact on the effective transmission rate, hence herd immunity 
# effects are limited. As such, the CEA results are (more) in line with the CEA 
# based on the static disease model.

# If the target group is mainly a sink for infections, and not driving the 
# transmission dynamics,vaccine-related immunization has limited impact on 
# the effective transmission rate, hence herd immunity effects are limited. 
# As such, the CEA results are (moreà in line with the CEA based on the static 
# disease model.

# If the (basic) reproduction number is close to one, vaccine-related immunization 
# can tip the balance quickly and result in large herd immunity effects, hence more 
# positive CEA outcomes. 

dynamic_out <-rbind(run_dynamic_model(uptake = 0.1, initial_R = 0),
                    run_dynamic_model(uptake = 0.2, initial_R = 0))

rsv_incidence_new <- as.numeric(dynamic_out[1,'reference_infections'] / population_size)

static_out <- rbind(run_static_model(uptake = 0.1, rsv_incidence = rsv_incidence_new),
                    run_static_model(uptake = 0.2, rsv_incidence = rsv_incidence_new))


kable(cbind(model='static',static_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
static	0.1	50	20000	4557901	4284427	15368	-6016429	-391	313379
static	0.2	50	20000	4557901	4010953	30736	-12032859	-391	626759

kable(cbind(model='dynamic',dynamic_out))

model	uptake	admin_cost	wtp	reference_infections	program_infections	qaly_gain	incr_cost	icer	inmb_k
dynamic	0.1	50	20000	4557901	4284427	15368	-6016397	-391	313379
dynamic	0.2	50	20000	4557901	4010953	30736	-12032792	-391	626758

The Role of Discounting: Reflect on the decision not to incorporate discounting. Offer your recommendation on whether discounting should be added and provide a rationale.

# Discounting does affect future costs and effects, in this case the life years lost due to disease-related mortality.

The Role of Uncertainty:

Reflect on the decision when

the cost per case is represented by rgamma(100, shape=(122/70)^2, rate = 122/70^2))
the vaccine efficacy by rnorm(100,vaccine_efficacy,0.1).

Try to reproduce the figure below for 100 samples:

# Note: it is also possible to use the summary functions above
# For illustration purposes, we explicitly include all code here.

num_samples <- 100
set.seed(1234)
unc_cost_per_episode <- rgamma(num_samples, shape=(122/70)^2, rate = 122/70^2)
unc_vaccine_efficacy <- rnorm(num_samples,vaccine_efficacy,0.1)

# static
unc_static_prevented_cases = population_size * (uptake * unc_vaccine_efficacy) * force_of_infection
unc_static_incr_qaly       = unc_static_prevented_cases * qaly_loss_per_episode
unc_static_incr_cost       = program_costs - (unc_static_prevented_cases * unc_cost_per_episode)
unc_static_icer            = mean(unc_static_incr_cost)  / mean(unc_static_incr_qaly)

# dynamic
unc_dynamic_prevented_cases <- vector(length=num_samples)
initial_prop_R = 0.95099
for(i in 1:num_samples){
 unc_proportion_protected = uptake * unc_vaccine_efficacy[i]
 initials['I'] = 1e-06
 initials['S'] = (1-initial_prop_R) * (1-unc_proportion_protected)
 initials['R'] = 1 - sum(initials['S']+initials['I'])
 dynamic_rsv = SIR(pars = parameters,init = initials,0:time_horizon)
 dynamic_cases_program = diff(dynamic_rsv$results$R[c(1,nrow(dynamic_rsv$results))]) * population_size
 unc_dynamic_prevented_cases[i] = dynamic_cases_reference - dynamic_cases_program
} 

# prevented cases
unc_dynamic_incr_qaly       = unc_dynamic_prevented_cases * qaly_loss_per_episode
unc_dynamic_incr_cost       = program_costs - (unc_dynamic_prevented_cases * unc_cost_per_episode)
unc_dynamic_icer            = mean(unc_dynamic_incr_cost)  / mean(unc_dynamic_incr_qaly)


plot(x=unc_static_incr_qaly,
     y=unc_static_incr_cost,
      xlim = range(unc_static_incr_qaly,unc_dynamic_incr_qaly),
     # ylim = range(unc_static_incr_cost,unc_dynamic_incr_cost),
     ylim = c(4,6)*1e7,
     xlab = 'incremental QALY gain',
     ylab = "incremental cost (euro)",
     pch=1)
points(x=unc_dynamic_incr_qaly,
       y=unc_dynamic_incr_cost,
     col = 2,
     pch=20)
legend('topright',
       c('static model',
         'dynamic model'),
       col = 1:2,
       pch = c(1,20)
       )

Navigate

Top: Index — Results

CEA for RSV prevention based on static and dynamic disease models