rm(list=ls())
require(ggthemes)

## Loading required package: ggthemes

require(tidyverse)

## Loading required package: tidyverse

## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

require(ellipse)

## Loading required package: ellipse

main.results <- FALSE
pkg = list("tidyverse","deSolve","ellipse","readxl","lhs")
invisible(lapply(pkg, require, character.only = TRUE))

## Loading required package: deSolve

## 
## Attaching package: 'deSolve'

## The following object is masked from 'package:graphics':
## 
##     matplot

## Loading required package: readxl

## Loading required package: lhs

run.id.base <- "vogi-nber"

mainDir <- "/Volumes/ACCRE/projects/rise-simulation/run-data"
subDir <- run.id.base

if (!file.exists(file.path(mainDir,subDir))){
  dir.create(file.path(mainDir, subDir))
}

ss_death <- read.csv("../numerical/ss-death-2011.csv")
source("../numerical/model-n-simple.R")
source("../numerical/numerical-functions.R")
source("../sub-files/metamodeling-functions.R")
times    <- seq(0, 80, by = 2 / 365)

The Value of Geomic Information

lambda <- 125000

Visualizing CE of Genetic Testing in a Cost-Effectiveness Plane

Our estimate of the VOGI is based on the Net Monetary Benefit (NMB) as a summary measure of cost-effectiveness. NMB is a common metric summarizing (in dollar terms) the health benefits and costs of a given treatment alternative. For treatment option ($\alpha$) and model parameters $\mathbf{\pi} = \{ \pi_1 , ...., \pi_I \}$ π = define:

\[ (1) \quad \quad NMB(\alpha, \mathbf{\pi}) = \lambda \cdot Q(\alpha,\mathbf{\pi}) – C(\alpha, \mathbf{\pi}) \]

Where $Q(\alpha, \mathbf{\pi})$ is the QALY measure from the intervention, $C(\alpha, \mathbf{\pi})$ is the cost of the intervention, and $$ is the fixed value of the marginal societal willingness to pay for an incremental health improvement (i.e., $50,000/QALY). For a given $ $ we define the optimal treatment as that which provides the maximum NMB on average:

\[ (2) \quad \quad \textrm{max}_\alpha \quad E_{\mathbf{\pi}}[ \textrm{NMB}(\alpha, \mathbf{\pi}) ] \]

When estimating equation (2) for the baseline no-genotyping scenario, all model parameters $\mathbf{\pi}=\{ \pi_1,..., \pi_I}$ are held fixed at their mean value. In other words, we assume that clinical decisions on drug therapy are based on the average population-level risk of an adverse event and not on information on the patient’s individual risk, preferences or characteristics.

Conversely, to study the value of information we directly model heterogeneity in $\mathbf{\pi}$ by drawing values from the joint probability distribution $p(\mathbf{\pi})$ at each iteration of the model.97 To estimate the EVGT we will utilize genetic data from the simulated patient cohort and model an alternative scenario in which clinicians obtain and act on information on genotype for that patient. The underlying genetic results used for these analyses demonstrate reproducibility of greater than 99%, including greater than 99% concordance and greater than 95% call rate. Under this “perfect information” scenario the physician utilizes information on each patient’s genetic makeup and chooses the optimal drug at the time of initial prescription, thus minimizing the risk of an adverse event. The average NMB under this scenario is value of perfect information on genotype and is given by:

\[ (3) \quad \quad \textrm{E}_{\mathbf{\pi}} [ \textrm{max}_{\alpha} \quad \textrm{NMB}(\alpha, \mathbf{\pi})] \]

Where the expectation is taken over the joint probability distribution of $\mathbf{\pi}$ (which in this case is simply the distribution of patients with each drug phenotype).

We obtain an estimate of the EVGT by taking the difference between equations (3) and (2):

\[ (4) \quad \quad \textrm{EVGT} =\textrm{E}_{\mathbf{\pi}} [ \textrm{max}_{\alpha} \quad \textrm{NMB}(\alpha, \mathbf{\pi})]- \textrm{max}_\alpha \quad E_{\mathbf{\pi}}[ \textrm{NMB}(\alpha, \mathbf{\pi}) ] \]

The EVGT provides an estimate of the opportunity cost (in dollar terms) of failing to incorporate genetic information into clinical decisions at the point of prescribing.98 In other words, it tells us maximum society would be willing to pay (per patient) to implement genotype-tailored care. An estimate of the total amount society (or a payer) would be willing to pay can be obtained by multiplying this value by the total number of patients and by the rate of indication incidence in the population.98

run.dir <- "/Volumes/ACCRE/projects/rise-simulation/"
#list.files(file.path(run.dir,"run-data/vogi-nber2/"))

todrop <- c("possible","fatal_b","living","disutil_a","disutil_b","dCOST.test","dCOST.drug","dCOST.treat")

load(file.path(run.dir,"run-data/vogi-nber2/nber-psa-none.RData"))
vogi.none <- psa[,-which(colnames(psa) %in% todrop)]
load(file.path(run.dir,"run-data/vogi-nber-seed1234/nber-psa-none.RData"))
vogi.none <- rbind(vogi.none,psa[,-which(colnames(psa) %in% todrop)])
load(file.path(run.dir,"run-data/vogi-nber-seed4321/nber-psa-none.RData"))
vogi.none <- rbind(vogi.none,psa[,-which(colnames(psa) %in% todrop)])


load(file.path(run.dir,"run-data/vogi-nber2/nber-psa-reactive-single.RData"))
vogi.single <-  psa[,-which(colnames(psa) %in% todrop)]
load(file.path(run.dir,"run-data/vogi-nber-seed1234/nber-psa-reactive-single.RData"))
vogi.single <- rbind(vogi.single,psa[,-which(colnames(psa) %in% todrop)])
load(file.path(run.dir,"run-data/vogi-nber-seed4321/nber-psa-reactive-single.RData"))
vogi.single <- rbind(vogi.single,psa[,-which(colnames(psa) %in% todrop)])


load(file.path(run.dir,"run-data/vogi-nber2/nber-psa-reactive-panel.RData"))
vogi.panel <-  psa[,-which(colnames(psa) %in% todrop)]
load(file.path(run.dir,"run-data/vogi-nber-seed1234/nber-psa-reactive-panel.RData"))
vogi.panel <- rbind(vogi.panel,psa[,-which(colnames(psa) %in% todrop)])
load(file.path(run.dir,"run-data/vogi-nber-seed4321/nber-psa-reactive-panel.RData"))
vogi.panel <- rbind(vogi.panel,psa[,-which(colnames(psa) %in% todrop)])


load(file.path(run.dir,"run-data/vogi-nber2/nber-psa-preemptive-panel.RData"))
vogi.prepanel<-  psa[,-which(colnames(psa) %in% todrop)]
load(file.path(run.dir,"run-data/vogi-nber-seed1234/nber-psa-preemptive-panel.RData"))
vogi.prepanel <- rbind(vogi.prepanel,psa[,-which(colnames(psa) %in% todrop)])
load(file.path(run.dir,"run-data/vogi-nber-seed4321/nber-psa-preemptive-panel.RData"))
vogi.prepanel <- rbind(vogi.prepanel,psa[,-which(colnames(psa) %in% todrop)])

ce.results <-
  plyr::rbind.fill(vogi.none,vogi.single,vogi.panel,vogi.prepanel) %>% mutate(strategy = gsub("-", ".", strategy)) 
varying1 <-
  suppressWarnings(ce.results %>% map_dbl( ~ var(., na.rm = TRUE)))
varying <- names(varying1[which(varying1 > 0)])
id.vars <-
  c("psa_id", "strategy", varying[-grep("psa_id|strategy|QALY|COST", varying)])
wide.fmla <-
  as.formula(paste0(paste0(c("psa_id", varying[-grep("psa_id|strategy|QALY|COST", varying)]), collapse =
                             "+"), "~variable+strategy"))
ce.results <-
  ce.results %>% reshape2::melt(id.vars = id.vars) %>% mutate(psa_id=as.character(psa_id)) %>% reshape2::dcast(wide.fmla)

ce.results %>% select(contains("dCOST"),contains("dQALY")) %>% colMeans() %>% data.frame() -> ce.results.summ
ce.results.summ$variable  = rownames(ce.results.summ)

ce.results.summ <- ce.results.summ %>% tidyr::separate(variable,c("variable","strategy"))
names(ce.results.summ) <- c("value","variable","strategy") 
ce.results.summ$strategy = gsub("dCOST_|dQALY_","",rownames(ce.results.summ)) 
ce.results.summ %>% tbl_df() %>% reshape2::dcast(strategy~variable) %>% Get_ICER()

##           strategy    dCOST    dQALY     ICER dominated.strong
## 1 preemptive.panel 38599.88 22.46424 113734.9                0
## 2   reactive.panel 36931.60 22.44957 112060.1                0
## 3  reactive.single 36822.31 22.44844       NA                0
## 4             none 33211.77 22.41638       NA                0
##   dominated.extended
## 1                  0
## 2                  0
## 3                  1
## 4                  0

source("../sub-files/metamodeling-functions.R")
Strategies <- c("none","reactive.single","reactive.panel","preemptive.panel")
#Determine the number of strategies
ndep <- length(Strategies)

#Create vector of variable names
Names <- ce.results %>% data.frame() %>% dplyr::select(-psa_id) %>% colnames()
Parms <- ce.results[, -grep("psa_id|Iteration|QALY|COST", colnames(ce.results))]   
#Get parameter names
paramNames <- colnames(Parms)
indep <- ncol(Parms)
Outcomes <- ce.results[, grep("psa_id|strategy|QALY|COST", colnames(ce.results))] 

for (i in seq(length(Strategies)))
  ce.results[, paste0("NHB_", Strategies[i])] <-
  ce.results %>% dplyr::select(-dplyr::contains("NHB"),-dplyr::contains("NMB"),-contains("global.")) %>% dplyr::select_if(grepl(Strategies[i], names(.))) %>% mutate_at(vars(contains("COST")), funs(-1 * . / lambda)) %>% mutate(NHB = rowSums(.)) %>% pull(NHB)

for (i in seq(length(Strategies)))
  ce.results[, paste0("NMB_", Strategies[i])] <-
  ce.results %>% dplyr::select(-contains("NHB"),-contains("NMB"),-contains("global.")) %>% dplyr::select_if(grepl(Strategies[i], names(.))) %>% mutate_at(vars(contains("COST")), funs(-1 *.)) %>% mutate_at(vars(contains("QALY")), funs(lambda * .)) %>%  mutate(NMB =rowSums(.)) %>% pull(NMB)

ce.results <- ce.results %>% mutate(Iteration=row_number())
NHB <- NMB <- NULL
NHB <- ce.results %>% dplyr::select(dplyr::contains("NHB"))
NMB <- ce.results %>% dplyr::select(dplyr::contains("NMB"))

nmb.gg <- reshape2::melt(NMB,  
               variable.name = "Strategy", 
               value.name = "NMB")
colnames(nmb.gg) <- c("Strategy","NMB")

## Plot NMB for different strategies
require(ggplot2)
require(scales)  # For dollar labels
require(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()

n.sim <- dim(NMB)[1]

Strategy1 <- "none"
Strategy2 <- "reactive.panel"

inmb <- NMB %>% dplyr::select_at(vars(contains(Strategy1),contains(Strategy2))) %>% mutate_at(vars(contains(Strategy2)),funs(-1*.)) %>%  mutate(Net_NMB=rowSums(.)) %>% 
  mutate(Simulation = row_number()) %>% dplyr::select(Simulation,Net_NMB)

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

## Format data frame suitably for plotting
inmb.gg <- reshape2::melt(inmb, id.vars = "Simulation", 
                variable.name = "Comparison", 
                value.name = "INMB")
colnames(inmb.gg) <- c("Simulation","Comparison","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)

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


# Or without iterating (much faster!)
loss <- as.matrix(NMB - NMB[, d.star])

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

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

require(tidyr)
require(broom)
# Unscaled #
Cost.Fit <- Strategies %>% purrr::map(~broom::tidy(lm(
  as.formula(paste0("dCOST_",.x,"~",paste0(colnames(Parms),collapse="+")))
  ,data=ce.results)))
Effectiveness.Fit <-Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("dQALY_",.x,"~",paste0(colnames(Parms),collapse="+")))
  ,data=ce.results)))
NHB.Fit <- Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("NHB_",.x,"~",paste0(colnames(Parms),collapse="+")))
  ,data=ce.results)))

# Scaled
Cost.Scaled.Fit <- Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("dCOST_",.x,"~",paste0(paste0("scale(",colnames(Parms),")"),collapse="+")))
  ,data=ce.results)))
Effectiveness.Scaled.Fit <- Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("dQALY_",.x,"~",paste0(paste0("scale(",colnames(Parms),")"),collapse="+")))
  ,data=ce.results)))
NHB.Scaled.Fit <- Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("NHB_",.x,"~",paste0(paste0("scale(",colnames(Parms),")"),collapse="+")))
  ,data=ce.results)))

One Way Sensitivity Analysis

Sim.long<- reshape2::melt(ce.results, id.vars=c("Iteration",grep("Probability",Names,value=TRUE)),measure.vars=paste0("NHB_",Strategies))

parm<-'risk.vGene_SC_B'
OneWaySA(Strategies,Parms,NHB,parm)

# Need to re-scale versions of each parameter
Parms %>% map_df(~mean(.x)) -> Parms.mean
Parms %>% map_df(~sd(.x)) -> Parms.sd

NHB.Scaled.Fit <- Strategies %>% purrr::map(~tidy(lm(
  as.formula(paste0("NHB_",.x,"~",paste0(paste0("scale(",colnames(Parms),")"),collapse="+")))
  ,data=ce.results)))
names(NHB.Scaled.Fit) = Strategies
Strategy.Combinations <- t(combn(Strategies,2))  %>% data.frame() %>% filter(as.character(X1)!=as.character(X2)) %>% 
  mutate_all(.funs=as.character)
delta <- map2(.x=as.character(Strategy.Combinations$X1),.y=as.character(Strategy.Combinations$X2),.f=~NHB.Scaled.Fit[[.x]]$estimate-NHB.Scaled.Fit[[.y]]$estimate) 
names(delta) = paste0(as.character(Strategy.Combinations$X1),"_vs_",as.character(Strategy.Combinations$X2))
delta <- delta %>% bind_rows()
zeta <- delta %>% purrr::map(~-.x[1]/.x[-1]) %>% bind_rows() 
zeta.rescaled <- t(Parms.mean)+t(Parms.sd)*zeta
rownames(zeta.rescaled) <- names(Parms)

Tornado Diagram

TornadoOpt(Parms=Parms,Outcomes=Outcomes[,-1])

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

# TornadoAll(Strategies,Parms[,-which(names(Parms)=="global.single_test")],NHB)
TornadoAll(Strategies,Parms,NHB)

Cost-Effectiveness Scatterplot

 ScatterCE(Strategies,Outcomes)

Two-Way Sensitivity Analysis

parm1<-'cost.alt_SC_A'
parm2<-'disutility.B_Survive_SC_A'
range1<-range(Parms[,parm1])
range2<-range(Parms[,parm2])
TwoWaySA(Strategies,Parms,Outcomes=NHB,parm1,parm2,range1,range2)

RISE Simulation v. 2.0