Home page

This script was developed by Lander Willem for illustration and teaching purposes only. Please make contact if you would like to use this model in any public way.

Goal

The goal is to illustrate the do’s and don’ts related to sampling for a probabilistic uncertainty analysis applied to the cost-effectiveness of a vaccination program for the elderly against influenza.

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

Input

We specify the target population and program details in which the vaccine is for free (in the current perspective) and there are no adverse events. Hence we expect only health and monetary gains. We aim to include uncertainty for the force of infection (mean = 4, sample size 60) and for the cost per case (mean = 10, se = 5).

# set number of samples to include in the uncertainty analysis
nsamples = 7

# population details
targetgroup     = 10000         # people>65 years

# force of infection
foi_mean        = 4   # cases
foi_sample_size = 60  # trial population

# cost per case
cost_mean = 10 # euro
cost_se   = 5  # euro

# program details
uptake_program          = 0.80
vaccine_efficacy        = 0.6
vaccine_price_per_dose  = 0 # euro

# QALY
QALYloss_case = 0.082

Inluding uncertainty: Don’t

Calculate the QALY loss and cost for the comparator and intervention:

# set the random-number-generator seed to get exactly the same results each time you run the code
set.seed(345698) 

# comparator
force_of_infection   = sample_beta_event(nsamples,foi_mean,foi_sample_size)  
cost_per_case        = sample_gamma(nsamples, cost_mean, cost_se)     # euro
comp_infections      = force_of_infection * targetgroup
comp_QALY_loss       = comp_infections * QALYloss_case
comp_total_cost      = comp_infections * cost_per_case

# intervention
force_of_infection   = sample_beta_event(nsamples,foi_mean,foi_sample_size)  
cost_per_case        = sample_gamma(nsamples, cost_mean, cost_se)     # euro
interv_infections    = force_of_infection * targetgroup * (1 - vaccine_efficacy * uptake_program)
interv_QALY_loss     = interv_infections * QALYloss_case
interv_total_cost    = interv_infections * cost_per_case + (targetgroup * vaccine_price_per_dose * uptake_program)

# calculate incremental cost and effect
incr_QALY_gain = comp_QALY_loss - interv_QALY_loss
incr_cost =  interv_total_cost - comp_total_cost

# plot CE-plane
plot(x = incr_QALY_gain, y = incr_cost, xlim = c(-10,80), ylim = c(-12e3,5e3), 
        xlab = 'Incremental effect', ylab = 'Incremental cost', pch = 19)
abline(h=0,lty=2)
abline(v=0,lty=2)

# show table
data_out = data.frame(sample = 1:nsamples,
                      comp_total_cost = comp_total_cost,
                      interv_total_cost = interv_total_cost,
                      incr_cost = incr_cost,
                      comp_QALY_loss = comp_QALY_loss,
                      interv_QALY_loss = interv_QALY_loss,
                      incr_QALY_gain = incr_QALY_gain)

# print table
knitr::kable(round(data_out), caption = "Summary Table of Costs and QALY Values")
Summary Table of Costs and QALY Values
sample comp_total_cost interv_total_cost incr_cost comp_QALY_loss interv_QALY_loss incr_QALY_gain
1 5290 7958 2668 50 32 18
2 11765 2709 -9056 44 32 12
3 3035 2356 -680 20 23 -3
4 7605 3541 -4064 89 35 54
5 12790 1686 -11104 87 34 53
6 2105 1693 -412 47 21 26
7 14608 3146 -11462 94 16 78

How is it possible that we have less health and/or more costs with a free program and no adverse events??!

# check "baseline" infections

# Calculate baseline number of infections in program
###interv_infections = force_of_infection * targetgroup * (1 - vaccine_efficacy * uptake_program)
interv_infections_baseline = interv_infections / (1 - vaccine_efficacy * uptake_program)

# calculate difference for "current practice"
diff_current_practice = interv_infections_baseline - comp_infections 

# calculate difference with program in place
incr_infections_program =  interv_infections - comp_infections

data_out = data.frame(sample = 1:nsamples,
                      comp_infections = comp_infections,
                      interv_infections_baseline = interv_infections_baseline,
                      diff_current_practice = diff_current_practice,
                      interv_infections = interv_infections,
                      incr_infections_program = incr_infections_program
                       )

# print table
knitr::kable(round(data_out), caption = "Summary Table of Infections and QALY Values")
Summary Table of Infections and QALY Values
sample comp_infections interv_infections_baseline diff_current_practice interv_infections incr_infections_program
1 610 743 133 386 -223
2 534 748 213 389 -146
3 245 532 288 277 32
4 1085 810 -275 421 -664
5 1066 797 -269 414 -651
6 574 495 -79 257 -317
7 1143 379 -763 197 -945

Inluding uncertainty: Do

Include uncertainty from the input parameters of the influenza model:

# set the random-number-generator seed to get exactly the same results each time you run the code
set.seed(345698) 

# comparator
force_of_infection   = sample_beta_event(nsamples,foi_mean,foi_sample_size)  
cost_per_case        = sample_gamma(nsamples, cost_mean, cost_se)     # euro
comp_infections      = force_of_infection * targetgroup
comp_QALY_loss       = comp_infections * QALYloss_case
comp_total_cost      = comp_infections * cost_per_case

# intervention
####force_of_infection   = sample_beta_event(nsamples,foi_mean,foi_sample_size)  
####cost_per_case        = sample_gamma(nsamples, cost_mean, cost_se)     # euro
interv_infections    = comp_infections * (1 - vaccine_efficacy * uptake_program)
interv_QALY_loss     = interv_infections * QALYloss_case
interv_total_cost    = interv_infections * cost_per_case + (targetgroup * vaccine_price_per_dose * uptake_program)

# calculate incremental cost and effect
incr_QALY_gain = comp_QALY_loss - interv_QALY_loss
incr_cost =  interv_total_cost - comp_total_cost

# plot CE-plane
plot(incr_QALY_gain,incr_cost, xlim = c(-10,80), ylim = c(-12e3,5e3), 
        xlab = 'Incremental effect', ylab = 'Incremental cost', pch = 19)
abline(h=0,lty=2)
abline(v=0,lty=2)

# show table
data_out = data.frame(sample = 1:nsamples,
                      comp_total_cost = comp_total_cost,
                      interv_total_cost = interv_total_cost,
                      incr_cost = incr_cost,
                      comp_QALY_loss = comp_QALY_loss,
                      interv_QALY_loss = interv_QALY_loss,
                      incr_QALY_gain = incr_QALY_gain)

# print table
knitr::kable(round(data_out), caption = "Summary Table of Costs and QALY Values")
Summary Table of Costs and QALY Values
sample comp_total_cost interv_total_cost incr_cost comp_QALY_loss interv_QALY_loss incr_QALY_gain
1 5290 2751 -2539 50 26 24
2 11765 6118 -5647 44 23 21
3 3035 1578 -1457 20 10 10
4 7605 3955 -3651 89 46 43
5 12790 6651 -6139 87 45 42
6 2105 1094 -1010 47 24 23
7 14608 7596 -7012 94 49 45

Calculate INMB per sample

For illustrative purposes, assume a fixed program cost of 100,000 euro and WTP of 5000 euro per QALY gain. Calculate the incremental net monetary benefit per sample and the average. In addition, compare the INMB with the ICER. How certain are you about the cost-effective option?

# adjust program cost
incr_cost_adj = incr_cost + 100000

# set WTP threshold
k = 5000

inmb_sample = (incr_QALY_gain * k) - incr_cost_adj
inmb_pos    = as.numeric(inmb_sample > 0)

# show table
data_out = data.frame(sample = 1:nsamples,
                      incr_cost_adj = incr_cost_adj,
                      incr_QALY_gain = incr_QALY_gain,
                      inmb_sample = inmb_sample,
                      inmb_pos = inmb_pos)


# print table
knitr::kable(round(data_out), caption = "Summary Table per Sample (k = 5000)")
Summary Table per Sample (k = 5000)
sample incr_cost_adj incr_QALY_gain inmb_sample inmb_pos
1 97461 24 22522 1
2 94353 21 10831 1
3 98543 10 -50390 0
4 96349 43 117147 1
5 93861 42 115834 1
6 98990 23 13984 1
7 92988 45 131884 1 
# summary statistics
inmb_mean    = mean(incr_QALY_gain * k) - mean(incr_cost_adj)
inmb_mean2   = mean(inmb_sample)
icer         = mean(incr_cost_adj) / mean(incr_QALY_gain)

# show table
data_out = data.frame(inmb_mean = inmb_mean,
                      inmb_mean2 = inmb_mean2,
                      inmb_pos = sum(inmb_pos) / length(inmb_pos),
                      icer = icer)

# print table
knitr::kable(data_out, caption = "Summary Table (k = 5000)")
Summary Table (k = 5000)
inmb_mean inmb_mean2 inmb_pos icer
51687.46 51687.46 0.8571429 3251.028

What if WTP = 3500 euro/QALY

Repeat the calculation above, with a WTP of 4100 euro per QALY gain.

# adjust program cost
incr_cost_adj = incr_cost + 100000

# set WTP threshold
k = 3500

inmb_sample = (incr_QALY_gain * k) - incr_cost_adj
inmb_pos    = as.numeric(inmb_sample > 0)

# show table
data_out = data.frame(sample = 1:nsamples,
                      incr_cost_adj = incr_cost_adj,
                      incr_QALY_gain = incr_QALY_gain,
                      inmb_sample = inmb_sample,
                      inmb_pos = inmb_pos)


# print table
knitr::kable(round(data_out), caption = "Summary Table per Sample (k = 3500)")
Summary Table per Sample (k = 3500)
sample incr_cost_adj incr_QALY_gain inmb_sample inmb_pos
1 97461 24 -13473 0
2 94353 21 -20724 0
3 98543 10 -64836 0
4 96349 43 53098 1
5 93861 42 52925 1
6 98990 23 -19908 0
7 92988 45 64422 1 
# summary statistics
inmb_mean = mean(incr_QALY_gain * k) - mean(incr_cost_adj)
icer      = mean(incr_cost_adj) / mean(incr_QALY_gain)

# show table
data_out = data.frame(inmb_mean = inmb_mean,
                      inmb_pos = sum(inmb_pos) / length(inmb_pos),
                      icer = icer)

# print table
knitr::kable(data_out, caption = "Summary Table (k = 3500)")
Summary Table (k = 3500)
inmb_mean inmb_pos icer
7357.886 0.4285714 3251.028

Calculate Net Loss per sample

For illustrative purposes, assume a fixed program cost of 100,000 euro and WTP of 5000 euro per QALY gain. Calculate the Net Loss for WTP 0 till 5000.

# adjust program cost
incr_cost_adj = incr_cost + 100000

# set WTP threshold
k = seq(0,5000,50)

inmb_wtp = mean(incr_QALY_gain) * k - mean(incr_cost_adj)
inmb_pos = as.numeric(inmb_wtp > 0)

net_loss_comparator = (inmb_pos == TRUE) * inmb_wtp
net_loss_intervention = (inmb_pos == FALSE) * -inmb_wtp

# show table
data_out = data.frame(k = k,
                      inmb_wtp = inmb_wtp,
                      inmb_pos = inmb_pos,
                      net_loss_comparator = net_loss_comparator,
                      net_loss_intervention = net_loss_intervention)


# print table
knitr::kable(round(data_out), caption = "Summary Table Net Loss")
Summary Table Net Loss
k inmb_wtp inmb_pos net_loss_comparator net_loss_intervention
0 -96078 0 0 96078
50 -94600 0 0 94600
100 -93122 0 0 93122
150 -91645 0 0 91645
200 -90167 0 0 90167
250 -88690 0 0 88690
300 -87212 0 0 87212
350 -85734 0 0 85734
400 -84257 0 0 84257
450 -82779 0 0 82779
500 -81301 0 0 81301
550 -79824 0 0 79824
600 -78346 0 0 78346
650 -76868 0 0 76868
700 -75391 0 0 75391
750 -73913 0 0 73913
800 -72435 0 0 72435
850 -70958 0 0 70958
900 -69480 0 0 69480
950 -68002 0 0 68002
1000 -66525 0 0 66525
1050 -65047 0 0 65047
1100 -63569 0 0 63569
1150 -62092 0 0 62092
1200 -60614 0 0 60614
1250 -59136 0 0 59136
1300 -57659 0 0 57659
1350 -56181 0 0 56181
1400 -54704 0 0 54704
1450 -53226 0 0 53226
1500 -51748 0 0 51748
1550 -50271 0 0 50271
1600 -48793 0 0 48793
1650 -47315 0 0 47315
1700 -45838 0 0 45838
1750 -44360 0 0 44360
1800 -42882 0 0 42882
1850 -41405 0 0 41405
1900 -39927 0 0 39927
1950 -38449 0 0 38449
2000 -36972 0 0 36972
2050 -35494 0 0 35494
2100 -34016 0 0 34016
2150 -32539 0 0 32539
2200 -31061 0 0 31061
2250 -29583 0 0 29583
2300 -28106 0 0 28106
2350 -26628 0 0 26628
2400 -25150 0 0 25150
2450 -23673 0 0 23673
2500 -22195 0 0 22195
2550 -20718 0 0 20718
2600 -19240 0 0 19240
2650 -17762 0 0 17762
2700 -16285 0 0 16285
2750 -14807 0 0 14807
2800 -13329 0 0 13329
2850 -11852 0 0 11852
2900 -10374 0 0 10374
2950 -8896 0 0 8896
3000 -7419 0 0 7419
3050 -5941 0 0 5941
3100 -4463 0 0 4463
3150 -2986 0 0 2986
3200 -1508 0 0 1508
3250 -30 0 0 30
3300 1447 1 1447 0
3350 2925 1 2925 0
3400 4403 1 4403 0
3450 5880 1 5880 0
3500 7358 1 7358 0
3550 8836 1 8836 0
3600 10313 1 10313 0
3650 11791 1 11791 0
3700 13268 1 13268 0
3750 14746 1 14746 0
3800 16224 1 16224 0
3850 17701 1 17701 0
3900 19179 1 19179 0
3950 20657 1 20657 0
4000 22134 1 22134 0
4050 23612 1 23612 0
4100 25090 1 25090 0
4150 26567 1 26567 0
4200 28045 1 28045 0
4250 29523 1 29523 0
4300 31000 1 31000 0
4350 32478 1 32478 0
4400 33956 1 33956 0
4450 35433 1 35433 0
4500 36911 1 36911 0
4550 38389 1 38389 0
4600 39866 1 39866 0
4650 41344 1 41344 0
4700 42822 1 42822 0
4750 44299 1 44299 0
4800 45777 1 45777 0
4850 47254 1 47254 0
4900 48732 1 48732 0
4950 50210 1 50210 0
5000 51687 1 51687 0 
# Plot Net Loss
plot(k, net_loss_comparator, type = 'l', xlab = 'WTP threshold', ylab = 'Expected net loss', lwd=2)
lines(k, net_loss_intervention, xlab = 'WTP threshold', ylab = 'Expected net loss', col =  2, lty=2, lwd = 2)

# PLOT EVPI
evpi <- apply(cbind(net_loss_comparator,net_loss_intervention),1,min)

#plot subset
subs <- seq(1,length(evpi),length=15)
points(k[subs],evpi[subs],pch=20,col=3)

# add legend
legend('topleft',
       c('comparator',
       'intervention',
       'EVPI'),
       col=1:3,
       lwd=c(2,2,NA),
       pch=c(NA,NA,20),
       lty=1:2)

# Plot INMB
plot(k, inmb_wtp, type = 'l', xlab = 'WTP threshold', ylab = 'Expected INMB', lwd=2)
abline(h=0)
abline(v=icer, lty=3)
text(x = icer, y = 3e4, labels = round(icer))

EVPPI

Assume the adjusted incremental costs, including the fixed program costs. What model parameters have the highest expected value of partial perfect information? hence, which parameters drive the decision uncertainty?

The EVPPI can be obtained using an iterative process (see below) by calculating INMB for different parameter subsets. Or, we can adopt a generalized additive model (GAM) approximation using the calculated incremental cost and benefit with all sampled parameter values.

# load package for GAM
library(mgcv)

# specify WTP
opt_wtp <- seq(0,5000,500)
num_wtp <- length(opt_wtp)

# combine input parameters with uncertainty
model_param = data.frame(force_of_infection= force_of_infection,
                         cost_per_case = cost_per_case)

# specify the number parameter values
num_param <- ncol(model_param)

# 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 = (incr_QALY_gain * k) - incr_cost_adj
  
  # 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)


data.frame(inmb_sample = inmb_sample,
          model_param = model_param[,i],
          g.hat_new = g.hat_new)
##   inmb_sample model_param g.hat_new
## 1    22522.29    8.677543  55500.61
## 2    10830.77   22.013042  35523.82
## 3   -50389.53   12.405731  49915.72
## 4   117147.44    7.010440  57997.96
## 5   115833.75   12.003076  50518.91
## 6    13983.61    3.666492  63007.24
## 7   131883.85   12.784757  49347.94

Explore EVPPI

# start with empty plot, and add EVPPI one by one
plot(range(opt_wtp),range(evppi),col=0, xlab = "WTP",ylab="EVPPI (euro)")
for(i in 1:ncol(evppi)){
  lines(opt_wtp,evppi[,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)

Iterative EVPPI for one WTP level

vaccine_price_per_dose <- 12.5

# select one WTP
k = 2500

evvpi_vector <- vector(length=nsamples)  # program vs doing nothing

for(i in 1:nsamples){

# let's say cost_per_case is known
cost_per_case_fixed <- cost_per_case[i]

# account for uncertainty in force of infection...

# comparator
comp_infections      = force_of_infection * targetgroup
comp_QALY_loss       = comp_infections * QALYloss_case
comp_total_cost      = comp_infections * cost_per_case_fixed

# intervention
interv_infections    = comp_infections * (1 - vaccine_efficacy * uptake_program)
interv_QALY_loss     = interv_infections * QALYloss_case
interv_total_cost    = interv_infections * cost_per_case_fixed + (targetgroup * vaccine_price_per_dose * uptake_program)

# calculate incremental cost and effect
incr_QALY_gain = comp_QALY_loss - interv_QALY_loss
incr_cost =  interv_total_cost - comp_total_cost

# calculate inmb
inmb_sample = (incr_QALY_gain * k) - incr_cost

perfect.info  <- mean(pmax(inmb_sample,0)) # select best options
baseline      <- pmax(mean(inmb_sample),0) # mean

## estimate EVPPI 
evvpi_vector[i]  <- round(perfect.info - baseline, digits=4) 
}

mean(evvpi_vector)
## [1] 5967.334
cbind(opt_wtp,evppi)
##       opt_wtp force_of_infection cost_per_case
##  [1,]       0              0.000        0.0000
##  [2,]     500              0.000        0.0000
##  [3,]    1000              0.000        0.0000
##  [4,]    1500              0.000        0.0000
##  [5,]    2000              0.000        0.0000
##  [6,]    2500           5818.053        0.0000
##  [7,]    3000          15076.103        0.0000
##  [8,]    3500          16976.268      453.3523
##  [9,]    4000          11457.795        0.0000
## [10,]    4500           8224.912        0.0000
## [11,]    5000           7188.797        0.0000