The objective of this notebook is to explore probabilistic sensitivity analysis with a discrete event simulation model with the following parameters:

A population of 40-year old women are at risk for both secular death and an event (A) that happens at a 10% [8-12%] rate over a 5 year period.
All those who experience the event incur a cost of $10,000.
Among those who experience the event, there is a 5% [2-7%] case fatality rate.
Among those who survive, they experience a 0.25 [0.2-0.3] utility decrement for the remainder of their life.
Among the surviors, there is a second event (B) that occurs with probability 50% [30-70] over a 5 year period.
Event B has no case fatality but incurs a $25,000 cost and a 0.10 [0.7-0.13] disutility for 1 year
There is a new therapy available that reduces the rate of the first event by a relative risk of 0.5 [0.45-0.55], however the therapy costs $8,000.

This model was intially coded up as a DES and then a numerical solution was derived using a delay differential equation.

Numerical Solution to Model

Single Model Run

lambda <- 50000

# Baseline Run
params['r_a'] <- params['r_a']
params['c_t']  <- 0
A <- as.data.frame(expected(params))
NMB.A <- lambda*A["qaly",] - A["cost",]; NMB.A

## [1] 1786903

# With Treatment
params['r_a'] <- params['r_a']*0.85
params['c_t']  <- 3500
B <- as.data.frame(expected(params))
NMB.B <- lambda*B["qaly",] - B["cost",]; NMB.B

## [1] 1786037

(ICER = (B["cost",]- A["cost",]) / (B["qaly",]-A["qaly",]))

## [1] 68099.35

Probabilistic Sensitivity Analysis

# Parameter List
params.psa <- list(
  r_a = list( # Rate of healthy having event A [uinform(8,12)]
    value = inst_rate(0.1, 5),
    type = "uniform",
    parameter = "r_a",
    min = 0.05,
    max = 0.15
  ),
  r_b = list( # Rate of post-A  having event B [unifirm(30-70)]
    value = inst_rate(0.5, 5),
    type = "uniform",
    parameter = "r_a",
    min = 0.30,
    max = 0.70
    ),
  r_ad = list( # Rate of death as direct result of A [beta(50,950)]
    value = 0.05,
    type = "beta",
    parameter = "r_ad",
    shape1 = 5,
    shape2 = 95
    ),
  c_a = list( # Cost of event A
    value = 10000,
    type = "uniform",
    parameter = "c_a",
    min = 8000,
    max = 12000
    ),  
  c_b = list( # Cost of event B
    value = 15000,
    type = "uniform",
    parameter = "c_b",
    min = 13000,
    max = 17000
    ), 
  c_t = list( # Cost of therapy
    value = 3500,
    type = "uniform",
    parameter = "c_b",
    min = 1000,
    max = 5000
    ),
  d_a = list( # Permanent disutility for a [beta(25,75)]
    value = 0.05,
    type = "beta",
    parameter = "d_a",
    shape1 = 5,
    shape2 = 95  
    ),
  d_b = list(    # 1-year disutility for b  [beta(100,900)]
    value = 0.10,
    type = "beta",
    parameter = "d_b",
    shape1 = 10,
    shape2 = 90
  ),
  rr_t = list( #RR of reducing A = 0.5 [beta(50,50)]
    value = 1,
    type = "beta",
    parameter = "rr_t",
    shape1 = 85,
    shape2 = 15             
  ),
  disc_rate = list( # Discount Rate
    value = 1e-12,
    type = "uniform",
    parameter = "disc_rate",
    min = 1e-12,
    max = 1e-12
    )
)

params <- unlist(lapply(params.psa,function(x) x$value))

We’ll now draw a random hypercube sample based on the number of parameters listed in the list object. Once we draw the uniform draws from the hypercube sample, we’ll apply the specific quantile function for the model parameter to get its draw for that iteration of the model.

PSA.N = 1000
require(lhs)

## Loading required package: lhs

require(tidyverse)
X <- randomLHS(PSA.N, length(params.psa))
colnames(X) = names(params.psa)

lhc.draws.transformed <- cbind.data.frame(lapply(params.psa,function(x) 
  {
    if (x[["type"]]=="beta")
    {
      qbeta(X[,x[["parameter"]]],shape1=x[["shape1"]],shape2=x[["shape2"]])
    } 
    else if (x[["type"]]=="uniform")
    {
      qunif(X[,x[["parameter"]]],min=x[["min"]],max=x[["max"]])
    }   
  }
))
lhc.draws.transformed %>% head()

##          r_a       r_b       r_ad       c_a      c_b      c_t        d_a
## 1 0.09371402 0.4748561 0.06605514 11620.550 15757.18 3757.175 0.05258879
## 2 0.13033334 0.6213334 0.04451099  8557.211 16546.86 4546.860 0.03509287
## 3 0.06718434 0.3687374 0.03572792  8022.422 15296.05 3296.047 0.04213306
## 4 0.09074419 0.4629768 0.01775439  9687.381 15595.44 3595.443 0.04120641
## 5 0.05628350 0.3251340 0.04050782  8156.099 15238.22 3238.224 0.04682592
## 6 0.07429463 0.3971785 0.04466170 11977.241 13694.31 1694.315 0.05429759
##          d_b      rr_t disc_rate
## 1 0.14438305 0.8466872     1e-12
## 2 0.11809823 0.9079447     1e-12
## 3 0.08105263 0.8688823     1e-12
## 4 0.08035576 0.8652979     1e-12
## 5 0.12501556 0.8588456     1e-12
## 6 0.10450661 0.8667640     1e-12

# Set up Parallel Processing
# require(doParallel)
# nworkers <- detectCores()
# cl <- makePSOCKcluster(nworkers)
# clusterSetRNGStream(cl, c(1,2,3,4,5,6,7,8))
# registerDoParallel(cl)
# 
# # Split Into Chunks
# d = 1:dim(lhc.draws.transformed)[1]
# chunks = split(d,ceiling(seq_along(d)/100))
require(doParallel)

## Loading required package: doParallel

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

## Loading required package: parallel

if (!file.exists("./psa-toy-v2.Rdata"))
{
  PSA <- foreach (i = 1:dim(lhc.draws.transformed)[1],.combine=rbind) %do%
  {
    cat(i)
    # Standard of Care
    params = lhc.draws.transformed[i,]
    params$r_a = inst_rate(params$r_a, 5)
    params$r_b = inst_rate(params$r_b, 5)
    params['rr_t'] = 1
    params['r_a'] <- params['r_a']*params['rr_t']
    params['c_t']  <- 0
    A <- as.data.frame(expected(params))
    params.A <- params
    names(params.A) = paste0(names(params.A),".A")
    
    # Treatment
    params = lhc.draws.transformed[i,]
    params$r_a = inst_rate(params$r_a, 5)
    params$r_b = inst_rate(params$r_b, 5)
    params['r_a'] <- params['r_a']*params['rr_t']
    B <- as.data.frame(expected(params))
    params.B <- params
    names(params.B) = paste0(names(params.B),".B")
    
    
    out <- unlist(c(params.A,
                    params.B, 
                    cost.A = A["cost",] , qaly.A = A["qaly",] , Strategy.A = lambda*A["qaly",] - A["cost",],
                    cost.B = B["cost",] , qaly.B = B["qaly",] , Strategy.B = lambda*B["qaly",] - B["cost",]))
    out
  }
  
  PSA <- PSA %>% tbl_df()
  save(PSA,file="./psa-toy-v2.Rdata")
} else
{
  load("./psa-toy-v2.Rdata")
}

Probalisitc Sensitivity Analysis

load("./psa-toy-v2.Rdata")
select <- dplyr::select
mutate <- dplyr::mutate
Sim <- PSA %>% tbl_df()
Sim <- PSA %>% mutate(iteration=row_number()) %>% select(iteration,cost.A,qaly.A,cost.B,qaly.B,contains(".B"))

Strategies<-c("Standard", "Therapy")

#Determine the number of strategies
ndep<-length(Strategies)

#Create vector of variable names
Names <- colnames(Sim)
Parms <- Sim %>% select(-iteration,-cost.A,-cost.B,-qaly.A,-qaly.B,-Strategy.B,-disc_rate.B) %>% data.frame() #Sim2[,(ndep*2+2):dim(Sim2)[2]]
#Get parameter names
paramNames<-colnames(Parms)
indep<-ncol(Parms)
Outcomes <- Sim %>%  select(Standard_cost=cost.A,Standard_eff=qaly.A,Therapy_cost=cost.B,Therapy_eff=qaly.B) %>% data.frame

In the simulation file we only have each parameter’s sample and the resulting cost and effectiveness for each of the strategies for each set of the parameters’ realization. We need to construct the Net Health Benefit (NHB) for each of the strategies. The NHB for strategy $k$ is defined as \[ \bar{NHB}_k = \bar{E}_k - \bar{C}_k/\lambda ,\] where $\bar{E}_k$ and $\bar{C}_k$ are the expected effectiveness and cost for strategy $k$, respectively; and $\lambda$ is the willingness-to-pay threshold.

  lambda <- 50000

Calculate NHB For Each Strategy

  Sim$A_NHB <- Sim$qaly.A- Sim$cost.A/lambda
  Sim$B_NHB <- Sim$qaly.B- Sim$cost.B/lambda
  #Save the NHB for each strategy into a new data frame
  NHB <- Sim %>% select(A_NHB,B_NHB)

Linear Metamodel

# Identify the parameters that vary
paramvar <- Sim[,paramNames] %>% summarise_all(function(x) var(x))
paramNames.varying <- names(paramvar)[which(paramvar!=0)]

## Cost ##
YY = "cost.A"
fmla.unstd <- as.formula(paste0(YY,"~",paste0(paramNames.varying,collapse="+")))
fmla.std <- as.formula(paste0(YY,"~",paste0(paste0("scale(",paramNames.varying,")"),collapse="+")))
A.cost.fit <- lm(fmla.unstd, data=Sim)
A.cost.fit.std <- lm(fmla.std,data=Sim)

## Effectiveness ##
YY = "qaly.A"
fmla.unstd <- as.formula(paste0(YY,"~",paste0(paramNames.varying,collapse="+")))
fmla.std <- as.formula(paste0(YY,"~",paste0(paste0("scale(",paramNames.varying,")"),collapse="+")))
A.eff.fit <- lm(fmla.unstd, data=Sim)
A.eff.fit.std <- lm(fmla.std,data=Sim)

## NHB ##
YY = "A_NHB"
fmla.unstd <- as.formula(paste0(YY,"~",paste0(paramNames.varying,collapse="+")))
fmla.std <- as.formula(paste0(YY,"~",paste0(paste0("scale(",paramNames.varying,")"),collapse="+")))
A.nhb.fit <- lm(fmla.unstd, data=Sim)
A.nhb.fit.std <- lm(fmla.std,data=Sim)

stargazer::stargazer(A.cost.fit,A.eff.fit,A.nhb.fit,
          title="Regression Coefficients from Metamodels on the Outcomes of Chemotherapy",
          type='html',
          dep.var.caption=" With **unstandardized** predictors ",
          dep.var.labels=c("Cost","Effectiveness","NHB"),
          digits=2,
          digits.extra=1,
          keep.stat=c("rsq","n"),
          notes = c('NHB: Net health benefit  '),
          notes.align = 'l',
          notes.append = FALSE,
          omit.table.layout="#",
          report=c("vc"))

**Regression Coefficients from Metamodels on the Outcomes of Chemotherapy**

	With unstandardized predictors

	Cost	Effectiveness	NHB

r_a.B	91,924.05	-21.44	-23.28

r_b.B	-18.84	0.75	0.75

r_ad.B	-745.88	-3.10	-3.09

c_a.B	0.10	0.000	-0.000

c_b.B	0.05	0.000	0.000

c_t.B

d_a.B	-5.42	-3.19	-3.19

d_b.B	-13.59	-0.02	-0.02

rr_t.B	-1,889.35	0.44	0.48

Constant	2.28	35.98	35.98


Observations	1,000	1,000	1,000
R²	0.99	0.94	0.95

Note:	NHB: Net health benefit

stargazer::stargazer(A.cost.fit.std,A.eff.fit.std,A.nhb.fit.std,
          title="Regression Coefficients from Metamodels on the Outcomes of Chemotherapy",
          type='html',
          dep.var.caption=" With **standardized** predictors ",
          dep.var.labels=c("Cost","Effectiveness","NHB"),
          digits=2,
          digits.extra=1,
          keep.stat=c("rsq","n"),
          notes = c('NHB: Net health benefit  '),
          notes.align = 'l',
          notes.append = FALSE,
          omit.table.layout="#",
          report=c("vc"))

**Regression Coefficients from Metamodels on the Outcomes of Chemotherapy**

	With standardized predictors

	Cost	Effectiveness	NHB

scale(r_a.B)	510.87	-0.12	-0.13

scale(r_b.B)	-0.91	0.04	0.04

scale(r_ad.B)	-16.16	-0.07	-0.07

scale(c_a.B)	114.66	0.000	-0.002

scale(c_b.B)	53.83	0.001	0.000

scale(c_t.B)

scale(d_a.B)	-0.12	-0.07	-0.07

scale(d_b.B)	-0.41	-0.001	-0.001

scale(rr_t.B)	-67.10	0.02	0.02

Constant	1,701.93	35.78	35.74


Observations	1,000	1,000	1,000
R²	0.99	0.94	0.95

Note:	NHB: Net health benefit

One-Way Sensitivity Analysis

Sim.Long <- Sim %>% tbl_df() %>% select(one_of(c("iteration",paramNames.varying,"A_NHB","B_NHB"))) %>% data.frame() %>% 
  melt(id.vars=c("iteration",paramNames.varying),measure.vars= c("A_NHB","B_NHB"))

Next we define the parameter for which we desire to do the one-way sensitivity analysis and its range.

parm<-'d_a.B'

The graph over the parameter’s sample domain looks like

OneWaySA(Strategies,Parms,NHB,parm)

# Tornado Diagram

NHB <- NHB %>% data.frame()
Parms <- Parms %>% data.frame()
TornadoOpt(Parms[,paramNames.varying[-6]],NHB)

TornadoAll(Strategies,Parms[,paramNames.varying[-6]],NHB)

Two-Way Sensitivity Analysis

parm1<-'r_a.B'#5
parm2<-'d_a.B'#1
range1<-quantile(PSA$r_a.B,c(0.025,0.975))
range2<-quantile(PSA$d_a.B,c(0.025,0.975))

The two-way sensitivity analysis of the NHB on r_a.B and d_a.B looks like

  TwoWaySA(Strategies,Parms,Outcomes=NHB,parm1,parm2,range1,range2)

  #TwoWaySA(Strategies,Parms,Outcomes=NHB,parm1,parm2)#Verify why not having ranges doesn't work!

Threshold Analysis

To perfom a threshold analysis we need to do a linear metamodel over all the different strategies. We will do it only for the standardized parameters. We will also include the interactions of the variables to perform a two-way sensitivity analysis in the next section.

Cost Effectiveness Plane

PlaneCE(Strategies,Outcomes)

## No id variables; using all as measure variables
## No id variables; using all as measure variables

Cost-Effectiveness Scatterplot

ScatterCE(Strategies,Outcomes)

## No id variables; using all as measure variables
## No id variables; using all as measure variables

  lambda_range<-c(1000,150000,50)

  CEAC(lambda_range,Strategies,Outcomes)

Value of Information

Plot Net Monetary Benefit For Each Strategy

lambda <- 50000
n.sim <- nrow(PSA)
nmb <-  PSA %>% mutate(Standard = qaly.A*lambda-cost.A, Therapy = qaly.B*lambda - cost.B) %>% select(Standard, Therapy) %>% data.frame()

# Number of Strategies
n.strategies <- ncol(nmb)
n.strategies

## [1] 2

# Assign name of strategies
strategies <- c("Standard","Therapy")
colnames(nmb) <- strategies

## Format data frame suitably for plotting
require(reshape2)
nmb.gg <- melt(nmb,  
               variable.name = "Strategy", 
               value.name = "NMB")

## No id variables; using all as measure variables

## Plot NMB for different strategies
require(ggplot2)
require(scales)  # For dollar labels
require(grid)

## Loading required package: grid

number_ticks <- function(n) {function(limits) pretty(limits, n)} 
# Faceted plot by Strategy
ggplot(nmb.gg, aes(x = NMB/1000)) +
  geom_histogram(aes(y =..density..), col="black", fill = "gray") +
  geom_density(color = "red") +
  facet_wrap(~ Strategy, scales = "free_y") +
  xlab("Net Monetary Benefit (NMB) x10^3") +
  scale_x_continuous(breaks = number_ticks(5), labels = dollar) + 
  scale_y_continuous(breaks = number_ticks(5)) + 
  theme_bw()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#### Incremental NMB (INMB) ####
# Calculate INMB of B vs A
# Only B vs A but we could have plotted all combinations

inmb <- data.frame(Simulation = 1:n.sim,
                   `Therapy vs Standard` = nmb$`Therapy` - nmb$`Standard`) 

## Format data frame suitably for plotting
inmb.gg <- melt(inmb, id.vars = "Simulation", 
                variable.name = "Comparison", 
                value.name = "INMB")
txtsize<-16

## Plot INMB
ggplot(inmb.gg, aes(x = INMB/1000)) +
  geom_histogram(aes(y =..density..), col="black", fill = "gray") +
  geom_density(color = "red") +
  geom_vline(xintercept = 0, col = 4, size = 1.5, linetype = "dashed") +
  facet_wrap(~ Comparison, scales = "free_y") +
  xlab("Incremental Net Monetary Benefit (INMB) in thousand $") +
  scale_x_continuous(breaks = number_ticks(5), limits = c(-100, 100)) + 
  scale_y_continuous(breaks = number_ticks(5)) + 
  theme_bw(base_size = 14)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 1 rows containing missing values (geom_bar).

#### Loss Matrix ####  
# Find optimal strategy (d*) based on the highest expected NMB
d.star <- which.max(colMeans(nmb))
d.star

## Standard 
##        1

## Compute Loss matrix iterating over all strategies
# Initialize loss matrix of dimension: number of simulation by number of strategies
loss <- matrix(0, n.sim, n.strategies)
# Or without iterating (much faster!)
loss <- as.matrix(nmb - nmb[, d.star])
head(loss)

##      Standard    Therapy
## [1,]        0  -809.5403
## [2,]        0   327.5891
## [3,]        0  1821.2237
## [4,]        0  2900.0905
## [5,]        0 -1009.8404
## [6,]        0 -1752.3566

#### EVPI ####
## Find maximum loss overall strategies at each state of the world 
## (i.e., PSA sample)
require(matrixStats)
max.loss.i <- rowMaxs(loss)
head(max.loss.i)

## [1]    0.0000  327.5891 1821.2237 2900.0905    0.0000    0.0000

## Average across all states of the world
evpi <- mean(max.loss.i)
evpi

## [1] 469.6443

Simple Event Model Value of Information