Probabilistic Sensitivity Analysis: methods
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:
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 = 1000
# 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.082Inluding 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(head(round(data_out), 10), caption = "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 | 6309 | 1826 | -4483 | 50 | 29 | 21 |
| 2 | 5787 | 1174 | -4613 | 44 | 10 | 33 |
| 3 | 1299 | 5606 | 4307 | 20 | 38 | -18 |
| 4 | 10302 | 8096 | -2206 | 89 | 34 | 55 |
| 5 | 15214 | 9412 | -5802 | 87 | 43 | 44 |
| 6 | 5107 | 5230 | 122 | 47 | 21 | 26 |
| 7 | 12353 | 3146 | -9206 | 94 | 27 | 67 |
| 8 | 5925 | 1377 | -4548 | 54 | 10 | 44 |
| 9 | 3156 | 4377 | 1220 | 51 | 30 | 21 |
| 10 | 15648 | 2122 | -13526 | 71 | 16 | 55 |
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(head(round(data_out), 10), caption = "Summary Table of Infections and QALY Values")| sample | comp_infections | interv_infections_baseline | diff_current_practice | interv_infections | incr_infections_program |
|---|---|---|---|---|---|
| 1 | 610 | 684 | 74 | 356 | -254 |
| 2 | 534 | 246 | -288 | 128 | -407 |
| 3 | 245 | 883 | 639 | 459 | 215 |
| 4 | 1085 | 797 | -288 | 414 | -671 |
| 5 | 1066 | 1014 | -52 | 527 | -538 |
| 6 | 574 | 499 | -75 | 259 | -315 |
| 7 | 1143 | 630 | -513 | 327 | -815 |
| 8 | 661 | 244 | -417 | 127 | -534 |
| 9 | 621 | 699 | 78 | 363 | -257 |
| 10 | 870 | 375 | -496 | 195 | -676 |
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(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(head(round(data_out), 10), caption = "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 | 6309 | 3280 | -3028 | 50 | 26 | 24 |
| 2 | 5787 | 3009 | -2778 | 44 | 23 | 21 |
| 3 | 1299 | 675 | -623 | 20 | 10 | 10 |
| 4 | 10302 | 5357 | -4945 | 89 | 46 | 43 |
| 5 | 15214 | 7911 | -7303 | 87 | 45 | 42 |
| 6 | 5107 | 2656 | -2452 | 47 | 24 | 23 |
| 7 | 12353 | 6423 | -5929 | 94 | 49 | 45 |
| 8 | 5925 | 3081 | -2844 | 54 | 28 | 26 |
| 9 | 3156 | 1641 | -1515 | 51 | 26 | 24 |
| 10 | 15648 | 8137 | -7511 | 71 | 37 | 34 |
Calculate ICER per sample: don’t!
For illustrative purposes, assume a fixed program cost of 100,000 euro and WTP of 4000 euro per QALY gain. Calculate the ratio of means and mean of ratios with respect to the incremental costs and effects.
# program cost
program_cost = 1e5
# adjust program cost
incr_cost_adj = incr_cost + program_cost
# willingness to pay
k = 4000
# calculate ratio of means and mean of ratios
mean_of_ratios = mean(incr_cost_adj / incr_QALY_gain)
ratio_of_means = mean(incr_cost_adj) / mean(incr_QALY_gain)
# calcualte also INMB
inmb_mean = mean(incr_QALY_gain * k) - mean(incr_cost_adj)
# show table
data_out = data.frame(mean_of_ratios = mean_of_ratios,
ratio_of_means = ratio_of_means,
inmb_mean = inmb_mean)
# print table
knitr::kable(head(data_out, 10), caption = "Summary Table (k = 4000)")| mean_of_ratios | ratio_of_means | inmb_mean |
|---|---|---|
| 4943.034 | 3684.851 | 8279.476 |
Is the intervention cost-effective when the WTP is 4000 euro per QALY gain.
Calculate INMB per sample
Assume a 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?
# 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(head(round(data_out), 10), caption = "Summary Table per Sample (k = 5000)")| sample | incr_cost_adj | incr_QALY_gain | inmb_sample | inmb_pos |
|---|---|---|---|---|
| 1 | 96972 | 24 | 23011 | 1 |
| 2 | 97222 | 21 | 7961 | 1 |
| 3 | 99377 | 10 | -51223 | 0 |
| 4 | 95055 | 43 | 118442 | 1 |
| 5 | 92697 | 42 | 116998 | 1 |
| 6 | 97548 | 23 | 15425 | 1 |
| 7 | 94071 | 45 | 130801 | 1 |
| 8 | 97156 | 26 | 32876 | 1 |
| 9 | 98485 | 24 | 23649 | 1 |
| 10 | 92489 | 34 | 78806 | 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(head(data_out, 10), caption = "Summary Table (k = 5000)")| inmb_mean | inmb_mean2 | inmb_pos | icer |
|---|---|---|---|
| 34551.06 | 34551.06 | 0.686 | 3684.851 |
What if WTP = 3500 euro/QALY
Repeat the calculation above, with a WTP of 3500 euro per QALY gain.
# 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(head(round(data_out), 10), caption = "Summary Table per Sample (k = 3500)")| sample | incr_cost_adj | incr_QALY_gain | inmb_sample | inmb_pos |
|---|---|---|---|---|
| 1 | 96972 | 24 | -12984 | 0 |
| 2 | 97222 | 21 | -23594 | 0 |
| 3 | 99377 | 10 | -65669 | 0 |
| 4 | 95055 | 43 | 54393 | 1 |
| 5 | 92697 | 42 | 54089 | 1 |
| 6 | 97548 | 23 | -18467 | 0 |
| 7 | 94071 | 45 | 63340 | 1 |
| 8 | 97156 | 26 | -6133 | 0 |
| 9 | 98485 | 24 | -12991 | 0 |
| 10 | 92489 | 34 | 27418 | 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(head(data_out, 10), caption = "Summary Table (k = 3500)")| inmb_mean | inmb_pos | icer |
|---|---|---|
| -4856.318 | 0.414 | 3684.851 |
Calculate Net Loss per sample
Still assume a fixed program cost of 100,000 euro. Calculate the Net Loss for WTP 0 till 5000 euro per QALY.
# set WTP threshold
k = seq(0, 5000, 500)
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(head(round(data_out), 10), caption = "Summary Table Net Loss")| k | inmb_wtp | inmb_pos | net_loss_comparator | net_loss_intervention |
|---|---|---|---|---|
| 0 | -96807 | 0 | 0 | 96807 |
| 500 | -83671 | 0 | 0 | 83671 |
| 1000 | -70535 | 0 | 0 | 70535 |
| 1500 | -57399 | 0 | 0 | 57399 |
| 2000 | -44264 | 0 | 0 | 44264 |
| 2500 | -31128 | 0 | 0 | 31128 |
| 3000 | -17992 | 0 | 0 | 17992 |
| 3500 | -4856 | 0 | 0 | 4856 |
| 4000 | 8279 | 1 | 8279 | 0 |
| 4500 | 21415 | 1 | 21415 | 0 |
# Plot Net Loss
plot(x = k,
y = net_loss_comparator,
type = 'l',
xlab = 'WTP threshold',
ylab = 'Expected net loss',
lwd = 2)
lines(x = k,
y = 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(x = k,
y = 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)Explore EVPPI
# start with empty plot, and add EVPPI one by one
plot(x = range(opt_wtp),
y = range(evppi),
col = 0,
xlab = "WTP",
ylab = "EVPPI (euro)")
for(i in 1:ncol(evppi)){
lines(x = opt_wtp,
y = 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 and one parameter
# select one WTP
k = 2500
g_hat = vector(length = nsamples)
for(i in 1:nsamples){
foi_fixed = force_of_infection[i] # fix FOI
comp_infections = foi_fixed * targetgroup # scalar
incr_QALY_gain = comp_infections * QALYloss_case * vaccine_efficacy * uptake_program
incr_cost_adj = -vaccine_efficacy * uptake_program * comp_infections * cost_per_case + program_cost # vector over all cost samples
inmb_sample = (incr_QALY_gain * k) - incr_cost_adj
g_hat[i] = mean(inmb_sample) # E[NMB | FOI_i]
}
perfect.info = mean(pmax(g_hat, 0))
baseline = pmax(mean(g_hat), 0)
evppi_foi = perfect.info - baseline
evppi_foi## [1] 3779.916## opt_wtp force_of_infection cost_per_case
## [1,] 0 0.0000 0
## [2,] 500 0.0000 0
## [3,] 1000 0.0000 0
## [4,] 1500 107.2607 0
## [5,] 2000 1096.8666 0
## [6,] 2500 3758.6191 0
## [7,] 3000 8603.1494 0
## [8,] 3500 15489.2750 0
## [9,] 4000 15789.8221 0
## [10,] 4500 12396.8706 0
## [11,] 5000 9875.9433 0