Optimizing the current model of care for knee osteoarthritis by the implementation of guideline recommended non-surgical care. A model-based health economic evaluation focusing on knee replacement delay.
Thomas Vetsch, Jan Taeymans, Nathanael Lutz
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
setwd("C:/Users/schwa/Desktop/markov_ONC_02122022/markov_ONC")Description
This file contains all R-Scripts that were used for modelling of the input parameters and creating plots for the manuscript.
Initialisation
In the first script all inputs for the model are defined. Age-specific transition probabilities are retrieved from stored .csv files. Values for costs, utilities and initial model counts are added directly in the script.
## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------------------------
library(tidyverse)
library(knitr)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------------------------
setwd("C:/Users/schwa/Desktop/markov_ONC_02122022/markov_ONC/input_files")
bfs_mort <- read.csv("p_mort.csv", sep = ";")
prob_SUR_uc <- read.csv("p_SUR_uc.csv", sep = ";")
prob_SUR_oc <- read.csv("p_SUR_oc.csv", sep = ";")
prob_REV_uc <- read.csv("p_REV_uc.csv", sep = ";")
prob_REV_oc <- read.csv("p_REV_oc.csv", sep = ";")
prob_PREV_m <- read.csv("PREV_m.csv", sep = ";")
prob_PREV_f <- read.csv("PREV_f.csv", sep = ";")
prob_PREV_m$value <- 0.02 #annual incidence for individuals aged 20+
prob_PREV_f$value <- 0.02
prob_prev_m_se <- 0.0512
prop_prev_f_se <- 0.0512
age_init <- 40 #Starting age chosen at 40 years. CAVE: Darf nicht tiefer sein als Alterskategorien aus .csv
# age increases with cycles
#cost usual care
c_NOSUR_uc <- 222
c_NOSUR_uc_sd <- 50
c_SUCC_NOSUR_uc <- 200
c_SUCC_NOSUR_uc_sd <- 20
c_SUR_uc <- 18325
c_SUR_uc_sd <- 1606/sqrt(1222)
c_SUCC_SUR_uc <- 500
c_SUCC_SUR_uc_sd <- 50
c_REV_uc <- 28775
c_REV_uc_sd <- 6436/sqrt(144)
c_SUCC_REV_uc <- 500
c_SUCC_REV_uc_sd <- 20
#cost optimized care
c_NOSUR_oc <- 1209
c_NOSUR_oc_sd <- 150
c_SUCC_NOSUR_oc <- 200
c_SUCC_NOSUR_oc_sd <- 20
c_SUR_oc <- 18325
c_SUR_oc_sd <- 1606/sqrt(1222)
c_SUCC_SUR_oc <- 500
c_SUCC_SUR_oc_sd <- 20
c_REV_oc <- 28775
c_REV_oc_sd <- 6436/sqrt(144)
c_SUCC_REV_oc <- 500
c_SUCC_REV_oc_sd <- 20
#utilities usual care
u_NOSUR_uc <- 0.658
u_SUCC_NOSUR_uc <- 0.749
u_SUCC_NOSUR_uc_sd <- 0.158/sqrt(41)
u_SUR_uc <- 0.658
u_SUCC_SUR_uc <- 0.878
u_SUCC_SUR_uc_sd <- 0.151/sqrt(41)
u_REV_uc <- 0.658
u_SUCC_REV_uc <- 0.658
#utilities optimized care
u_NOSUR_oc <- 0.658
u_SUCC_NOSUR_oc <- 0.783
u_SUCC_NOSUR_oc_sd <- 0.108/sqrt(41)
u_SUR_oc <- 0.658
u_SUCC_SUR_oc <- 0.878
u_SUCC_SUR_oc_sd <- 0.151/sqrt(41)
u_REV_oc <- 0.658
u_SUCC_REV_oc <- 0.658
#start count
n_male <- 60164
n_female <- 59508Model calculations
Current model of care
Now the markov model is run for 70 cycles under current conditions. This means we use cost and utility values defined for “current non-surgical care”.
## ----setup, include=FALSE-----------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE----------------------------------------------------------
cycle <- 1:70
n_cycle <- length(cycle)
POP_fu <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fu[i+1] <- POP_fu[i] -
(POP_fu[i] * prob_PREV_f[prob_PREV_f$age == age,2]) -
(POP_fu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mu <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mu[i+1] <- POP_mu[i] -
(POP_mu[i] * prob_PREV_m[prob_PREV_m$age == age,2]) -
(POP_mu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURu[i+1] <- (POP_fu[i]*prob_PREV_f[prob_PREV_f$age == age,2]) +
(POP_mu[i]*prob_PREV_m[prob_PREV_m$age == age,2]) -
(NOSURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURu[i+1] <- (SUCC_NOSURu[i] +NOSURu[i]) -
(SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SUCC_NOSURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURu[i+1] <- (SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURu[i+1] <- (SUCC_SURu[i] +SURu[i]) -
(SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(SUCC_SURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVu[i+1] <- (SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(REVu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVu[i+1] <- (SUCC_REVu[i] +REVu[i]) -
(SUCC_REVu[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHu[i+1] <- DEATHu[i] +
POP_fu[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mu[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURu[i] *bfs_mort[bfs_mort$age == age,2]+
REVu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVu[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}Optimized model of care
Following the modelling under “current non-surgical care” conditions, the model is run for “optimized model of care” conditions. Meaning that we use the input parameters that we’re defined for the “optimized model of care” (Higher costs for non-surgical treatment, higher utilities for non-surgical treatment). The transition probabilities remain the same as under “current model of care”.
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo[i+1] <- POP_fo[i] -
(POP_fo[i] * prob_PREV_f[prob_PREV_f$age == age,2]) -
(POP_fo[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo[i+1] <- POP_mo[i] -
(POP_mo[i] * prob_PREV_m[prob_PREV_m$age == age,2]) -
(POP_mo[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo[i+1] <- (POP_fo[i]*prob_PREV_f[prob_PREV_f$age == age,2]) +
(POP_mo[i]*prob_PREV_m[prob_PREV_m$age == age,2]) -
(NOSURo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo[i+1] <- (SUCC_NOSURo[i] +NOSURo[i])-
(SUCC_NOSURo[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo[i+1] <- (SUCC_NOSURo[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(SURo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo[i+1] <- (SUCC_SURo[i] +SURo[i]) -
(SUCC_SURo[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo[i+1] <- (SUCC_SURo[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo[i+1] <- (SUCC_REVo[i] +REVo[i]) -
(SUCC_REVo[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo[i+1] <- DEATHo[i] +
POP_fo[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo[i] *bfs_mort[bfs_mort$age == age,2]+
SURo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo[i] *bfs_mort[bfs_mort$age == age,2]+
REVo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}Optimized model of care +2 year delay
For the scenario analysis, the model is again run for 70 cycles with values for “optimized model of care”. To simulate a 2 year delay of total knee replacement we add two “tunnel states” after “successful non-surgical care”. The transition probabilities remain the same.
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo2 <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo2[i+1] <- POP_fo2[i] -
(POP_fo2[i] * prob_PREV_f[prob_PREV_f$age == age,2]) -
(POP_fo2[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo2 <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo2[i+1] <- POP_mo2[i] -
(POP_mo2[i] * prob_PREV_m[prob_PREV_m$age == age,2]) -
(POP_mo2[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo2[i+1] <- (POP_fo2[i]*prob_PREV_f[prob_PREV_f$age == age,2]) +
(POP_mo2[i]*prob_PREV_m[prob_PREV_m$age == age,2]) -
(NOSURo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo2[i+1] <- (SUCC_NOSURo2[i] +NOSURo2[i])-
(SUCC_NOSURo2[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del1_2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del1_2[i+1] <- (SUCC_NOSURo2[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(del1_2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del2_2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del2_2[i+1] <- (del1_2[i]) -
(del2_2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo2[i+1] <- (del2_2[i]) -
(SURo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo2[i+1] <- (SUCC_SURo2[i] +SURo2[i]) -
(SUCC_SURo2[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo2[i+1] <- (SUCC_SURo2[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo2[i+1] <- (SUCC_REVo2[i] +REVo2[i]) -
(SUCC_REVo2[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo2[i+1] <- DEATHo2[i] +
POP_fo2[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo2[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo2[i] *bfs_mort[bfs_mort$age == age,2]+
del1_2[i] *bfs_mort[bfs_mort$age == age,2]+
del2_2[i] *bfs_mort[bfs_mort$age == age,2]+
SURo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo2[i] *bfs_mort[bfs_mort$age == age,2]+
REVo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo2[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}Optimized model of care +5 year delay
Same procedure as described above. But now we add 5 “tunnel states” to simualte a 5 year delay of total knee replacement. Again, all transition probabilities remain the same.
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo5 <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo5[i+1] <- POP_fo5[i] -
(POP_fo5[i] * prob_PREV_f[prob_PREV_f$age == age,2]) -
(POP_fo5[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo5 <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo5[i+1] <- POP_mo5[i] -
(POP_mo5[i] * prob_PREV_m[prob_PREV_m$age == age,2]) -
(POP_mo5[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo5[i+1] <- (POP_fo5[i]*prob_PREV_f[prob_PREV_f$age == age,2]) +
(POP_mo5[i]*prob_PREV_m[prob_PREV_m$age == age,2]) -
(NOSURo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo5[i+1] <- (SUCC_NOSURo5[i] +NOSURo5[i])-
(SUCC_NOSURo5[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del1_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del1_5[i+1] <- (SUCC_NOSURo5[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(del1_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del2_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del2_5[i+1] <- (del1_5[i]) -
(del2_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del3_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del3_5[i+1] <- (del2_5[i]) -
(del3_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del4_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del4_5[i+1] <- (del3_5[i]) -
(del4_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del5_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del5_5[i+1] <- (del4_5[i]) -
(del5_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo5[i+1] <- (del5_5[i]) -
(SURo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo5[i+1] <- (SUCC_SURo5[i] +SURo5[i]) -
(SUCC_SURo5[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo5[i+1] <- (SUCC_SURo5[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo5[i+1] <- (SUCC_REVo5[i] +REVo5[i]) -
(SUCC_REVo5[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo5[i+1] <- DEATHo5[i] +
POP_fo5[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo5[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo5[i] *bfs_mort[bfs_mort$age == age,2]+
del1_5[i] *bfs_mort[bfs_mort$age == age,2]+
del2_5[i] *bfs_mort[bfs_mort$age == age,2]+
del3_5[i] *bfs_mort[bfs_mort$age == age,2]+
del4_5[i] *bfs_mort[bfs_mort$age == age,2]+
del5_5[i] *bfs_mort[bfs_mort$age == age,2]+
SURo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo5[i] *bfs_mort[bfs_mort$age == age,2]+
REVo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo5[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}Results
QALY calculation for all scenarios
QALY per person is calculated for each model.
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
discount <- 0.03
qaly_uc <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
qaly_uc[i] <- POP_fu[i] *u_NOSUR_uc+
POP_mu[i] *u_NOSUR_uc+
NOSURu[i] *u_NOSUR_uc+
SUCC_NOSURu[i] *u_SUCC_NOSUR_uc+
SURu[i] *u_NOSUR_uc+
SUCC_SURu[i] *u_SUCC_SUR_uc+
REVu[i] *u_NOSUR_uc+
SUCC_REVu[i] *u_SUCC_REV_uc
}
qaly_disc_uc <- vector()
for(i in 1:n_cycle){
qaly_disc_uc[i] <- qaly_uc[i]/(1+discount)^cycle_years[i]
}
total_qaly_uc <- sum(qaly_disc_uc)
person_qaly_uc <- total_qaly_uc/119672
cat("Model = ","UNC",", total qaly = ", total_qaly_uc, ", person qaly", person_qaly_uc)## Model = UNC , total qaly = 2024442 , person qaly 16.91659## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
discount <- 0.03
qaly_oc <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
qaly_oc[i] <- POP_fo[i] *u_NOSUR_oc+
POP_mo[i] *u_NOSUR_oc+
NOSURo[i] *u_NOSUR_oc+
SUCC_NOSURo[i] *u_SUCC_NOSUR_oc+
SURo[i] *u_NOSUR_oc+
SUCC_SURo[i] *u_SUCC_SUR_oc+
REVo[i] *u_NOSUR_oc+
SUCC_REVo[i] *u_SUCC_REV_oc
}
qaly_disc_oc <- vector()
for(i in 1:n_cycle){
qaly_disc_oc[i] <- qaly_oc[i]/(1+discount)^cycle_years[i]
}
total_qaly_oc <- sum(qaly_disc_oc)
person_qaly_oc <- total_qaly_oc/119672
cat("Model = ","ONC no delay",", total qaly = ", total_qaly_oc, ", person qaly", person_qaly_oc)## Model = ONC no delay , total qaly = 2042947 , person qaly 17.07122## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
discount <- 0.03
qaly_oc2 <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
qaly_oc2[i] <-POP_fo2[i] *u_NOSUR_oc+
POP_mo2[i] *u_NOSUR_oc+
NOSURo2[i] *u_NOSUR_oc+
SUCC_NOSURo2[i] *u_SUCC_NOSUR_oc+
del1_2[i] *u_SUCC_NOSUR_oc+
del2_2[i] *u_SUCC_NOSUR_oc+
SURo2[i] *u_NOSUR_oc+
SUCC_SURo2[i] *u_SUCC_SUR_oc+
REVo2[i] *u_NOSUR_oc+
SUCC_REVo2[i] *u_SUCC_REV_oc
}
qaly_disc_oc2 <- vector()
for(i in 1:n_cycle){
qaly_disc_oc2[i] <- qaly_oc2[i]/(1+discount)^cycle_years[i]
}
total_qaly_oc2 <- sum(qaly_disc_oc2)
person_qaly_oc2 <- total_qaly_oc2/119672
cat("Model = ","ONC 2 year delay",", total qaly = ", total_qaly_oc2, ", person qaly", person_qaly_oc2)## Model = ONC 2 year delay , total qaly = 2041053 , person qaly 17.05539## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
discount <- 0.03
qaly_oc5 <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
qaly_oc5[i] <- POP_fo5[i] *u_NOSUR_oc+
POP_mo5[i] *u_NOSUR_oc+
NOSURo5[i] *u_NOSUR_oc+
SUCC_NOSURo5[i] *u_SUCC_NOSUR_oc+
del1_5[i] *u_SUCC_NOSUR_oc+
del2_5[i] *u_SUCC_NOSUR_oc+
del3_5[i] *u_SUCC_NOSUR_oc+
del4_5[i] *u_SUCC_NOSUR_oc+
del5_5[i] *u_SUCC_NOSUR_oc+
SURo5[i] *u_NOSUR_oc+
SUCC_SURo5[i] *u_SUCC_SUR_oc+
REVo5[i] *u_NOSUR_oc+
SUCC_REVo5[i] *u_SUCC_REV_oc
}
qaly_disc_oc5 <- vector()
for(i in 1:n_cycle){
qaly_disc_oc5[i] <- qaly_oc5[i]/(1+discount)^cycle_years[i]
}
total_qaly_oc5 <- sum(qaly_disc_oc5)
person_qaly_oc5 <- total_qaly_oc5/119672
cat("Model = ","ONC 5 year delay",", total qaly = ", total_qaly_oc5, ", person qaly", person_qaly_oc5)## Model = ONC 5 year delay , total qaly = 2038503 , person qaly 17.03409Cost calculation for all scenarios
Cost per person is calculated for each model.
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----fig.cap= "", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE-----------------
cost_uc <- vector()
cycle_years <- rep(0:70)
discount <- 0.03
for(i in 1:n_cycle){
cost_uc[i] <-
NOSURu[i] *c_NOSUR_uc+
SUCC_NOSURu[i] *c_SUCC_NOSUR_uc+
SURu[i] *c_SUR_uc+
SUCC_SURu[i] *c_SUCC_SUR_uc+
REVu[i] *c_REV_uc+
SUCC_REVu[i] *c_SUCC_REV_uc
}
cost_disc_uc <- vector()
for(i in 1:n_cycle){
cost_disc_uc[i] <- cost_uc[i]/(1+discount)^cycle_years[i]
}
total_cost_uc <- sum(cost_disc_uc)
person_cost_uc <- total_cost_uc/119672
cat("Model = ","UNC",", total cost = ", total_cost_uc, ", person cost", person_cost_uc)## Model = UNC , total cost = 534362192 , person cost 4465.223## ----fig.cap= "", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE-----------------
cost_oc <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
cost_oc[i] <-
NOSURo[i] *c_NOSUR_oc+
SUCC_NOSURo[i] *c_SUCC_NOSUR_oc+
SURo[i] *c_SUR_oc+
SUCC_SURo[i] *c_SUCC_SUR_oc+
REVo[i] *c_REV_oc+
SUCC_REVo[i] *c_SUCC_REV_oc
}
cost_disc_oc <- vector()
for(i in 1:n_cycle){
cost_disc_oc[i] <- cost_oc[i]/(1+discount)^cycle_years[i]
}
total_cost_oc <- sum(cost_disc_oc)
person_cost_oc <- total_cost_oc/119672
cat("Model = ","ONC no delay",", total cost = ", total_cost_oc, ", person cost", person_cost_oc)## Model = ONC no delay , total cost = 575123320 , person cost 4805.83## ----fig.cap= "", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE-----------------
cost_oc2 <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
cost_oc2[i] <-
NOSURo2[i] *c_NOSUR_oc+
SUCC_NOSURo2[i] *c_SUCC_NOSUR_oc+
del1_2[i] *c_SUCC_NOSUR_oc+
del2_2[i] *c_SUCC_NOSUR_oc+
SURo2[i] *c_SUR_oc+
SUCC_SURo2[i] *c_SUCC_SUR_oc+
REVo2[i] *c_REV_oc+
SUCC_REVo2[i] *c_SUCC_REV_oc
}
cost_disc_oc2 <- vector()
for(i in 1:n_cycle){
cost_disc_oc2[i] <- cost_oc2[i]/(1+discount)^cycle_years[i]
}
total_cost_oc2 <- sum(cost_disc_oc2)
person_cost_oc2 <- total_cost_oc2/119672
cat("Model = ","ONC 2 year delay",", total cost = ", total_cost_oc2, ", person cost", person_cost_oc2)## Model = ONC 2 year delay , total cost = 532684821 , person cost 4451.207## ----fig.cap= "", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE-----------------
cost_oc5 <- vector()
cycle_years <- rep(0:70)
for(i in 1:n_cycle){
cost_oc5[i] <-
NOSURo5[i] *c_NOSUR_oc+
SUCC_NOSURo5[i] *c_SUCC_NOSUR_oc+
del1_5[i] *c_SUCC_NOSUR_oc+
del2_5[i] *c_SUCC_NOSUR_oc+
del3_5[i] *c_SUCC_NOSUR_oc+
del4_5[i] *c_SUCC_NOSUR_oc+
del5_5[i] *c_SUCC_NOSUR_oc+
SURo5[i] *c_SUR_oc+
SUCC_SURo5[i] *c_SUCC_SUR_oc+
REVo5[i] *c_REV_oc+
SUCC_REVo5[i] *c_SUCC_REV_oc
}
cost_disc_oc5 <- vector()
for(i in 1:n_cycle){
cost_disc_oc5[i] <- cost_oc5[i]/(1+discount)^cycle_years[i]
}
total_cost_oc5 <- sum(cost_disc_oc5)
person_cost_oc5 <- total_cost_oc5/119672
cat("Model = ","ONC 5 year delay",", total cost = ", total_cost_oc5, ", person cost", person_cost_oc5)## Model = ONC 5 year delay , total cost = 474355741 , person cost 3963.799# o5_cost_nosur <- sum(NOSURo5) * c_NOSUR_oc
# o5_cost_nosursuc <-
# (sum(SUCC_NOSURo5) + sum(del1_5) + sum(del2_5) + sum(del3_5) + sum(del4_5) + sum(del5_5)) *
# c_SUCC_NOSUR_oc
# o5_cost_sur <- sum(SURo5) * c_SUR_oc
# o5_cost_sursuc <- sum(SUCC_SURo5) * c_SUCC_SUR_oc
# o5_cost_rev <- sum(REVo5) * c_REV_oc
# o5_cost_revsuc <- sum(SUCC_REVo5) * c_SUCC_REV_oc
#
# o2_cost_nosur <- sum(NOSURo2) * c_NOSUR_oc
# o2_cost_nosursuc <-
# (sum(SUCC_NOSURo2) + sum(del1_2) + sum(del2_2)) * c_SUCC_NOSUR_oc
# o2_cost_sur <- sum(SURo2) * c_SUR_oc
# o2_cost_sursuc <- sum(SUCC_SURo2) * c_SUCC_SUR_oc
# o2_cost_rev <- sum(REVo2) * c_REV_oc
# o2_cost_revsuc <- sum(SUCC_REVo2) * c_SUCC_REV_oc
#
# oc_cost_nosur <- sum(NOSURo) * c_NOSUR_oc
# oc_cost_nosursuc <- sum(SUCC_NOSURo) * c_SUCC_NOSUR_oc
# oc_cost_sur <- sum(SURo) * c_SUR_oc
# oc_cost_sursuc <- sum(SUCC_SURo) * c_SUCC_SUR_oc
# oc_cost_rev <- sum(REVo) * c_REV_oc
# oc_cost_revsuc <- sum(SUCC_REVo) * c_SUCC_REV_oc
#
# uc_cost_nosur <- sum(NOSURu) * c_NOSUR_uc
# uc_cost_nosursuc <- sum(SUCC_NOSURu) * c_SUCC_NOSUR_uc
# uc_cost_sur <- sum(SURu) * c_SUR_uc
# uc_cost_sursuc <- sum(SUCC_SURu) * c_SUCC_SUR_uc
# uc_cost_rev <- sum(REVu) * c_REV_uc
# uc_cost_revsuc <- sum(SUCC_REVu) * c_SUCC_REV_uc
#
#
#
# total_uc <-
# sum(
# uc_cost_nosur,
# uc_cost_nosursuc,
# uc_cost_sur,
# uc_cost_sursuc,
# uc_cost_rev,
# uc_cost_revsuc
# )
# total_oc <-
# sum(
# oc_cost_nosur,
# oc_cost_nosursuc,
# oc_cost_sur,
# oc_cost_sursuc,
# oc_cost_rev,
# oc_cost_revsuc
# )
# total_o2 <-
# sum(
# o2_cost_nosur,
# o2_cost_nosursuc,
# o2_cost_sur,
# o2_cost_sursuc,
# o2_cost_rev,
# o2_cost_revsuc
# )
# total_o5 <-
# sum(
# o5_cost_nosur,
# o5_cost_nosursuc,
# o5_cost_sur,
# o5_cost_sursuc,
# o5_cost_rev,
# o5_cost_revsuc
# )
#
#
#
# df_costs <-
# data.frame(
# name = c(
# "uc_cost_nosur",
# "uc_cost_nosursuc",
# "uc_cost_sur",
# "uc_cost_sursuc",
# "uc_cost_rev",
# "uc_cost_revsuc",
# "oc_cost_nosur",
# "oc_cost_nosursuc",
# "oc_cost_sur",
# "oc_cost_sursuc",
# "oc_cost_rev",
# "oc_cost_revsuc",
# "o2_cost_nosur",
# "o2_cost_nosursuc",
# "o2_cost_sur",
# "o2_cost_sursuc",
# "o2_cost_rev",
# "o2_cost_revsuc",
# "o5_cost_nosur",
# "o5_cost_nosursuc",
# "o5_cost_sur",
# "o5_cost_sursuc",
# "o5_cost_rev",
# "o5_cost_revsuc"
# ),
# cost = c(
# uc_cost_nosur,
# uc_cost_nosursuc,
# uc_cost_sur,
# uc_cost_sursuc,
# uc_cost_rev,
# uc_cost_revsuc,
# oc_cost_nosur,
# oc_cost_nosursuc,
# oc_cost_sur,
# oc_cost_sursuc,
# oc_cost_rev,
# oc_cost_revsuc,
# o2_cost_nosur,
# o2_cost_nosursuc,
# o2_cost_sur,
# o2_cost_sursuc,
# o2_cost_rev,
# o2_cost_revsuc,
# o5_cost_nosur,
# o5_cost_nosursuc,
# o5_cost_sur,
# o5_cost_sursuc,
# o5_cost_rev,
# o5_cost_revsuc
# )
# )
#
# write.csv(df_costs, "export_new.csv")
#
# df_cycle_cost <-
# cbind(cost_oc,
# cost_disc_oc,
# cost_oc2,
# cost_disc_oc2,
# cost_oc5,
# cost_disc_oc5)
#
# write.csv(data_test, "cycle_cost.csv")
#
# counts_per_state <-
# cbind(NOSURo, SUCC_NOSURo, SURo, SUCC_SURo, REVo, SUCC_REVo)
#
# write.csv(counts_per_state, "counts_opti.csv")ICER calculation for all scenarios
The ICER is calculated for each scenario and used to plot on the cost-effectiveness panel
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----fig.cap= "counts by state", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE----
delta_costs_base <- person_cost_oc - person_cost_uc
delta_qalys_base <- person_qaly_oc - person_qaly_uc
icer_base <- delta_costs_base/delta_qalys_base
cat("Scenario = ","basecase",", Cost difference = ", delta_costs_base, ", QALY difference", delta_qalys_base, ", ICER = ", icer_base)## Scenario = basecase , Cost difference = 340.6071 , QALY difference 0.154632 , ICER = 2202.694delta_costs_scen2 <- person_cost_oc2 - person_cost_uc
delta_qalys_scen2 <- person_qaly_oc2 - person_qaly_uc
icer_scen2 <- delta_costs_scen2/delta_qalys_scen2
cat("Scenario = ","2 year delay",", Cost difference = ", delta_costs_scen2, ", QALY difference", delta_qalys_scen2, ", ICER = ", icer_scen2)## Scenario = 2 year delay , Cost difference = -14.01641 , QALY difference 0.138808 , ICER = -100.977delta_costs_scen5 <- person_cost_oc5 - person_cost_uc
delta_qalys_scen5 <- person_qaly_oc5 - person_qaly_uc
icer_scen5 <- delta_costs_scen5/delta_qalys_scen5
cat("Scenario = ","5 year delay",", Cost difference = ", delta_costs_scen5, ", QALY difference", delta_qalys_scen5, ", ICER = ", icer_scen5)## Scenario = 5 year delay , Cost difference = -501.4243 , QALY difference 0.1175037 , ICER = -4267.305icer_table <- data.frame(scenario = c("base", "2 year delay", "5 year delay"),
cost_delta = c(delta_costs_base, delta_costs_scen2, delta_costs_scen5),
qaly_delta = c(delta_qalys_base, delta_qalys_scen2, delta_qalys_scen5),
ICER = c(icer_base, icer_scen2, icer_scen5))
per_person_table <- data.frame(scenario = c("base case UNC","base case ONC", "2 year delay", "5 year delay"),
cost_person = c(person_cost_uc, person_cost_oc, person_cost_oc2, person_cost_oc5),
qaly_person = c(person_qaly_uc, person_qaly_oc, person_qaly_oc2, person_qaly_oc5))Probabilistic Sensitivity Analysis
PSA base case scenario
Probabilistic sensitivity analysis with 1000 (10000 in manuscript) iterations on utility and cost values for the “basecase scenario”
#PSA TEST
## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------------------------
# Beta distribution
estBetaParams <- function(mu, sd) {
var <- sd^2
alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
beta <- alpha * (1 / mu - 1)
return(params = list(alpha = alpha, beta = beta))
}
#create beta distribution for utilities, for successfull nosur and sur
alpha_beta <- estBetaParams(u_SUCC_NOSUR_uc, u_SUCC_NOSUR_uc_sd)
dist_u_SUCC_NOSUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_NOSUR_oc, u_SUCC_NOSUR_oc_sd)
dist_u_SUCC_NOSUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_uc, u_SUCC_SUR_uc_sd)
dist_u_SUCC_SUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_oc, u_SUCC_SUR_oc_sd)
dist_u_SUCC_SUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#create gamma distributions for costs
dist_c_NOSUR_oc <- rgamma(100000, (c_NOSUR_oc_sd/c_NOSUR_oc)^-2, scale = c_NOSUR_oc * (c_NOSUR_oc_sd/c_NOSUR_oc)^2)
dist_c_NOSUR_uc <- rgamma(100000, (c_NOSUR_uc_sd/c_NOSUR_uc)^-2, scale = c_NOSUR_uc * (c_NOSUR_uc_sd/c_NOSUR_uc)^2)
dist_c_SUR_oc <- rgamma(100000, (c_SUR_oc_sd/c_SUR_oc)^-2, scale = c_SUR_oc * (c_SUR_oc_sd/c_SUR_oc)^2)
dist_c_SUR_uc <- rgamma(100000, (c_SUR_uc_sd/c_SUR_uc)^-2, scale = c_SUR_uc * (c_SUR_uc_sd/c_SUR_uc)^2)
dist_c_REV_oc <- rgamma(100000, (c_REV_oc_sd/c_REV_oc)^-2, scale = c_REV_oc * (c_REV_oc_sd/c_REV_oc)^2)
dist_c_REV_uc <- rgamma(100000, (c_REV_uc_sd/c_REV_uc)^-2, scale = c_REV_uc * (c_REV_oc_sd/c_REV_uc)^2)
#create gamma distributions for probabilities healthy to non-surgical treatment
#prob_PREV_m$value <- 0.02 #annual incidence for individuals aged 20+
#prob_PREV_f$value <- 0.02
#prob_prev_m_se <- 0.00512
#prop_prev_f_se <- 0.00512
alpha_beta <- estBetaParams(0.02, 0.00510)
dist_p_prev_m <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
dist_p_prev_f <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#define vectors for psa
psa_cost_disc_uc <- vector()
psa_cost_uc <- vector()
psa_cost_disc_oc <- vector()
psa_cost_oc <- vector()
psa_qaly_disc_oc <- vector()
psa_qaly_oc <- vector()
psa_qaly_disc_uc <- vector()
psa_qaly_uc <- vector()
res_cost_uc <- vector()
res_cost_oc <- vector()
res_qaly_uc <- vector()
res_qaly_oc <- vector()
discount <- 0.03
n_cycle <- 70
#run psa x times
for(j in 1:1000){
c_NOSUR_oc_psa <- sample(dist_c_NOSUR_oc,1)
c_NOSUR_uc_psa <- sample(dist_c_NOSUR_uc,1)
c_SUR_oc_psa <- sample(dist_c_SUR_oc,1)
c_SUR_uc_psa <- sample(dist_c_SUR_uc,1)
c_REV_oc_psa <- sample(dist_c_REV_oc,1)
c_REV_uc_psa <- sample(dist_c_REV_uc,1)
u_SUCC_NOSUR_oc_psa <- sample(dist_u_SUCC_NOSUR_oc,1)
u_SUCC_NOSUR_uc_psa <- sample(dist_u_SUCC_NOSUR_uc,1)
u_SUCC_SUR_oc_psa <- sample(dist_u_SUCC_SUR_oc,1)
u_SUCC_SUR_uc_psa <- sample(dist_u_SUCC_SUR_uc,1)
p_prev_m_psa <- sample(dist_p_prev_m,1)
p_prev_f_psa <- sample(dist_p_prev_f,1)
## RUN MODEL USUAL CARE
cycle <- 1:70
n_cycle <- length(cycle)
POP_fu <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fu[i+1] <- POP_fu[i] -
(POP_fu[i] * p_prev_f_psa) -
(POP_fu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mu <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mu[i+1] <- POP_mu[i] -
(POP_mu[i] * p_prev_m_psa) -
(POP_mu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURu[i+1] <- (POP_fu[i]*p_prev_f_psa) +
(POP_mu[i]*p_prev_m_psa) -
(NOSURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURu[i+1] <- (SUCC_NOSURu[i] +NOSURu[i]) -
(SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SUCC_NOSURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURu[i+1] <- (SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURu[i+1] <- (SUCC_SURu[i] +SURu[i]) -
(SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(SUCC_SURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVu[i+1] <- (SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(REVu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVu[i+1] <- (SUCC_REVu[i] +REVu[i]) -
(SUCC_REVu[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHu[i+1] <- DEATHu[i] +
POP_fu[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mu[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURu[i] *bfs_mort[bfs_mort$age == age,2]+
REVu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVu[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
## RUN MODEL OPTIMIZED CARE
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo[i+1] <- POP_fo[i] -
(POP_fo[i] * p_prev_f_psa) -
(POP_fo[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo[i+1] <- POP_mo[i] -
(POP_mo[i] * p_prev_m_psa ) -
(POP_mo[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo[i+1] <- (POP_fo[i]*p_prev_f_psa) +
(POP_mo[i]*p_prev_m_psa) -
(NOSURo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo[i+1] <- (SUCC_NOSURo[i] +NOSURo[i])-
(SUCC_NOSURo[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo[i+1] <- (SUCC_NOSURo[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(SURo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo[i+1] <- (SUCC_SURo[i] +SURo[i]) -
(SUCC_SURo[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo[i+1] <- (SUCC_SURo[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo[i+1] <- (SUCC_REVo[i] +REVo[i]) -
(SUCC_REVo[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo[i+1] <- DEATHo[i] +
POP_fo[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo[i] *bfs_mort[bfs_mort$age == age,2]+
SURo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo[i] *bfs_mort[bfs_mort$age == age,2]+
REVo[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
#cost calc for UNC
for(i in 1:n_cycle){
psa_cost_uc[i] <- NOSURu[i] *c_NOSUR_uc_psa+
SUCC_NOSURu[i] *c_SUCC_NOSUR_uc+
SURu[i] *c_SUR_uc_psa+
SUCC_SURu[i] *c_SUCC_SUR_uc+
REVu[i] *c_REV_uc_psa+
SUCC_REVu[i] *c_SUCC_REV_uc
}
#discouting costs
for(i in 1:n_cycle){
psa_cost_disc_uc[i] <- psa_cost_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_uc <- sum(psa_cost_disc_uc)
psa_person_cost_uc <- psa_total_cost_uc/119672
res_cost_uc[j] <- psa_person_cost_uc
# cost calc for ONC
for(i in 1:n_cycle){
psa_cost_oc[i] <- NOSURo[i] *c_NOSUR_oc_psa+
SUCC_NOSURo[i] *c_SUCC_NOSUR_oc+
SURo[i] *c_SUR_oc_psa+
SUCC_SURo[i] *c_SUCC_SUR_oc+
REVo[i] *c_REV_oc_psa+
SUCC_REVo[i] *c_SUCC_REV_oc
}
#discounting costs
for(i in 1:n_cycle){
psa_cost_disc_oc[i] <- psa_cost_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_oc <- sum(psa_cost_disc_oc)
psa_person_cost_oc <- psa_total_cost_oc/119672
res_cost_oc[j] <- psa_person_cost_oc
#utility calc for UNC
for(i in 1:n_cycle){
psa_qaly_uc[i] <- POP_fu[i] *u_NOSUR_uc+
POP_mu[i] *u_NOSUR_uc+
NOSURu[i] *u_NOSUR_uc+
SUCC_NOSURu[i] *u_SUCC_NOSUR_uc_psa+
SURu[i] *u_NOSUR_uc+
SUCC_SURu[i] *u_SUCC_SUR_uc_psa+
REVu[i] *u_NOSUR_uc+
SUCC_REVu[i] *u_SUCC_REV_uc
}
for(i in 1:n_cycle){
psa_qaly_disc_uc[i] <- psa_qaly_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_uc <- sum(psa_qaly_disc_uc)
psa_person_qaly_uc <- psa_total_qaly_uc/119672
res_qaly_uc[j] <- psa_person_qaly_uc
#utility calc for ONC
for(i in 1:n_cycle){
psa_qaly_oc[i] <- POP_fo[i] *u_NOSUR_oc+
POP_mo[i] *u_NOSUR_oc+
NOSURo[i] *u_NOSUR_oc+
SUCC_NOSURo[i] *u_SUCC_NOSUR_oc_psa+
SURo[i] *u_NOSUR_oc+
SUCC_SURo[i] *u_SUCC_SUR_oc_psa+
REVo[i] *u_NOSUR_oc+
SUCC_REVo[i] *u_SUCC_REV_oc
}
for(i in 1:n_cycle){
psa_qaly_disc_oc[i] <- psa_qaly_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_oc <- sum(psa_qaly_disc_oc)
psa_person_qaly_oc <- psa_total_qaly_oc/119672
res_qaly_oc[j] <- psa_person_qaly_oc
}
mean(res_qaly_oc)## [1] 17.06963mean(res_cost_oc)## [1] 4771.391psabc_delta_costs <- res_cost_oc - res_cost_uc
psabc_delta_qalys <- res_qaly_oc - res_qaly_uc
kable(per_person_table, digits =2)| scenario | cost_person | qaly_person |
|---|---|---|
| base case UNC | 4465.22 | 16.92 |
| base case ONC | 4805.83 | 17.07 |
| 2 year delay | 4451.21 | 17.06 |
| 5 year delay | 3963.80 | 17.03 |
PSA 2 year delay scenario
Probabilistic sensitivity analysis with 1000 (10000 in manuscript) iterations on utility and cost values for the “2-year delay scenario”
#PSA TEST
## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------------------------
# Beta distribution
estBetaParams <- function(mu, sd) {
var <- sd^2
alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
beta <- alpha * (1 / mu - 1)
return(params = list(alpha = alpha, beta = beta))
}
#create beta distribution for utilities, for successfull nosur and sur
alpha_beta <- estBetaParams(u_SUCC_NOSUR_uc, u_SUCC_NOSUR_uc_sd)
dist_u_SUCC_NOSUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_NOSUR_oc, u_SUCC_NOSUR_oc_sd)
dist_u_SUCC_NOSUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_uc, u_SUCC_SUR_uc_sd)
dist_u_SUCC_SUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_oc, u_SUCC_SUR_oc_sd)
dist_u_SUCC_SUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#create gamma distributions for costs
dist_c_NOSUR_oc <- rgamma(100000, (c_NOSUR_oc_sd/c_NOSUR_oc)^-2, scale = c_NOSUR_oc * (c_NOSUR_oc_sd/c_NOSUR_oc)^2)
dist_c_NOSUR_uc <- rgamma(100000, (c_NOSUR_uc_sd/c_NOSUR_uc)^-2, scale = c_NOSUR_uc * (c_NOSUR_uc_sd/c_NOSUR_uc)^2)
dist_c_SUR_oc <- rgamma(100000, (c_SUR_oc_sd/c_SUR_oc)^-2, scale = c_SUR_oc * (c_SUR_oc_sd/c_SUR_oc)^2)
dist_c_SUR_uc <- rgamma(100000, (c_SUR_uc_sd/c_SUR_uc)^-2, scale = c_SUR_uc * (c_SUR_uc_sd/c_SUR_uc)^2)
dist_c_REV_oc <- rgamma(100000, (c_REV_oc_sd/c_REV_oc)^-2, scale = c_REV_oc * (c_REV_oc_sd/c_REV_oc)^2)
dist_c_REV_uc <- rgamma(100000, (c_REV_uc_sd/c_REV_uc)^-2, scale = c_REV_uc * (c_REV_oc_sd/c_REV_uc)^2)
#create gamma distributions for probabilities healthy to non-surgical treatment
#prob_PREV_m$value <- 0.02 #annual incidence for individuals aged 20+
#prob_PREV_f$value <- 0.02
#prob_prev_m_se <- 0.00512
#prop_prev_f_se <- 0.00512
alpha_beta <- estBetaParams(0.02, 0.00510)
dist_p_prev_m <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
dist_p_prev_f <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#define vectors for psa
psa_cost_disc_uc <- vector()
psa_cost_uc <- vector()
psa_cost_disc_oc <- vector()
psa_cost_oc <- vector()
psa_qaly_disc_oc <- vector()
psa_qaly_oc <- vector()
psa_qaly_disc_uc <- vector()
psa_qaly_uc <- vector()
res_cost_uc <- vector()
res_cost_oc2 <- vector()
res_qaly_uc <- vector()
res_qaly_oc2 <- vector()
discount <- 0.03
n_cycle <- 70
#run psa x times
for(j in 1:1000){
c_NOSUR_oc_psa <- sample(dist_c_NOSUR_oc,1)
c_NOSUR_uc_psa <- sample(dist_c_NOSUR_uc,1)
c_SUR_oc_psa <- sample(dist_c_SUR_oc,1)
c_SUR_uc_psa <- sample(dist_c_SUR_uc,1)
c_REV_oc_psa <- sample(dist_c_REV_oc,1)
c_REV_uc_psa <- sample(dist_c_REV_uc,1)
u_SUCC_NOSUR_oc_psa <- sample(dist_u_SUCC_NOSUR_oc,1)
u_SUCC_NOSUR_uc_psa <- sample(dist_u_SUCC_NOSUR_uc,1)
u_SUCC_SUR_oc_psa <- sample(dist_u_SUCC_SUR_oc,1)
u_SUCC_SUR_uc_psa <- sample(dist_u_SUCC_SUR_uc,1)
p_prev_m_psa <- sample(dist_p_prev_m,1)
p_prev_f_psa <- sample(dist_p_prev_f,1)
## RUN MODEL USUAL CARE
cycle <- 1:70
n_cycle <- length(cycle)
POP_fu <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fu[i+1] <- POP_fu[i] -
(POP_fu[i] * p_prev_f_psa) -
(POP_fu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mu <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mu[i+1] <- POP_mu[i] -
(POP_mu[i] * p_prev_m_psa) -
(POP_mu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURu[i+1] <- (POP_fu[i]*p_prev_f_psa) +
(POP_mu[i]*p_prev_m_psa) -
(NOSURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURu[i+1] <- (SUCC_NOSURu[i] +NOSURu[i]) -
(SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SUCC_NOSURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURu[i+1] <- (SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURu[i+1] <- (SUCC_SURu[i] +SURu[i]) -
(SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(SUCC_SURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVu[i+1] <- (SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(REVu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVu[i+1] <- (SUCC_REVu[i] +REVu[i]) -
(SUCC_REVu[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHu[i+1] <- DEATHu[i] +
POP_fu[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mu[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURu[i] *bfs_mort[bfs_mort$age == age,2]+
REVu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVu[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
## RUN MODEL OPTIMIZED CARE 2 year delay
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo2 <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo2[i+1] <- POP_fo2[i] -
(POP_fo2[i] * p_prev_f_psa) -
(POP_fo2[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo2 <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo2[i+1] <- POP_mo2[i] -
(POP_mo2[i] * p_prev_m_psa) -
(POP_mo2[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo2[i+1] <- (POP_fo2[i]*p_prev_f_psa) +
(POP_mo2[i]*p_prev_m_psa) -
(NOSURo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo2[i+1] <- (SUCC_NOSURo2[i] +NOSURo2[i])-
(SUCC_NOSURo2[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del1_2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del1_2[i+1] <- (SUCC_NOSURo2[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(del1_2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del2_2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del2_2[i+1] <- (del1_2[i]) -
(del2_2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo2[i+1] <- (del2_2[i]) -
(SURo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo2[i+1] <- (SUCC_SURo2[i] +SURo2[i]) -
(SUCC_SURo2[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo2[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo2[i+1] <- (SUCC_SURo2[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo2[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo2[i+1] <- (SUCC_REVo2[i] +REVo2[i]) -
(SUCC_REVo2[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo2 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo2[i+1] <- DEATHo2[i] +
POP_fo2[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo2[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo2[i] *bfs_mort[bfs_mort$age == age,2]+
del1_2[i] *bfs_mort[bfs_mort$age == age,2]+
del2_2[i] *bfs_mort[bfs_mort$age == age,2]+
SURo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo2[i] *bfs_mort[bfs_mort$age == age,2]+
REVo2[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo2[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
#cost calc for UNC
for(i in 1:n_cycle){
psa_cost_uc[i] <- NOSURu[i] *c_NOSUR_uc_psa+
SUCC_NOSURu[i] *c_SUCC_NOSUR_uc+
SURu[i] *c_SUR_uc_psa+
SUCC_SURu[i] *c_SUCC_SUR_uc+
REVu[i] *c_REV_uc_psa+
SUCC_REVu[i] *c_SUCC_REV_uc
}
#discouting costs UNC
for(i in 1:n_cycle){
psa_cost_disc_uc[i] <- psa_cost_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_uc <- sum(psa_cost_disc_uc)
psa_person_cost_uc <- psa_total_cost_uc/119672
res_cost_uc[j] <- psa_person_cost_uc
# cost calc for ONC
for(i in 1:n_cycle){
psa_cost_oc[i] <-
NOSURo2[i] *c_NOSUR_oc_psa+
SUCC_NOSURo2[i] *c_SUCC_NOSUR_oc+
del1_2[i] *c_SUCC_NOSUR_oc+
del2_2[i] *c_SUCC_NOSUR_oc+
SURo2[i] *c_SUR_oc_psa+
SUCC_SURo2[i] *c_SUCC_SUR_oc+
REVo2[i] *c_REV_oc_psa+
SUCC_REVo2[i] *c_SUCC_REV_oc
}
#discounting costs ONC
for(i in 1:n_cycle){
psa_cost_disc_oc[i] <- psa_cost_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_oc <- sum(psa_cost_disc_oc)
psa_person_cost_oc <- psa_total_cost_oc/119672
res_cost_oc2[j] <- psa_person_cost_oc
#utility calc for UNC
for(i in 1:n_cycle){
psa_qaly_uc[i] <- POP_fu[i] *u_NOSUR_uc+
POP_mu[i] *u_NOSUR_uc+
NOSURu[i] *u_NOSUR_uc+
SUCC_NOSURu[i] *u_SUCC_NOSUR_uc_psa+
SURu[i] *u_NOSUR_uc+
SUCC_SURu[i] *u_SUCC_SUR_uc_psa+
REVu[i] *u_NOSUR_uc+
SUCC_REVu[i] *u_SUCC_REV_uc
}
#discounting utility UNC
for(i in 1:n_cycle){
psa_qaly_disc_uc[i] <- psa_qaly_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_uc <- sum(psa_qaly_disc_uc)
psa_person_qaly_uc <- psa_total_qaly_uc/119672
res_qaly_uc[j] <- psa_person_qaly_uc
#utility calc for ONC
for(i in 1:n_cycle){
psa_qaly_oc[i] <- POP_fo2[i] *u_NOSUR_oc+
POP_mo2[i] *u_NOSUR_oc+
NOSURo2[i] *u_NOSUR_oc+
SUCC_NOSURo2[i] *u_SUCC_NOSUR_oc_psa+
del1_2[i] *u_SUCC_NOSUR_oc+
del2_2[i] *u_SUCC_NOSUR_oc+
SURo2[i] *u_NOSUR_oc+
SUCC_SURo2[i] *u_SUCC_SUR_oc_psa+
REVo2[i] *u_NOSUR_oc+
SUCC_REVo2[i] *u_SUCC_REV_oc
}
#discounting utility ONC
for(i in 1:n_cycle){
psa_qaly_disc_oc[i] <- psa_qaly_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_oc <- sum(psa_qaly_disc_oc)
psa_person_qaly_oc <- psa_total_qaly_oc/119672
res_qaly_oc2[j] <- psa_person_qaly_oc
}
mean(res_qaly_oc2)## [1] 17.04876mean(res_cost_oc2)## [1] 4405.482psa2_delta_costs <- res_cost_oc2 - res_cost_uc
psa2_delta_qalys <- res_qaly_oc2 - res_qaly_uc
kable(per_person_table, digits = 2)| scenario | cost_person | qaly_person |
|---|---|---|
| base case UNC | 4465.22 | 16.92 |
| base case ONC | 4805.83 | 17.07 |
| 2 year delay | 4451.21 | 17.06 |
| 5 year delay | 3963.80 | 17.03 |
PSA 5 year delay scenario
Probabilistic sensitivity analysis with 1000 (10000 in manuscript) iterations on utility and cost values for the 5-year delay scenario"
#PSA TEST
## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------------------------
# Beta distribution
estBetaParams <- function(mu, sd) {
var <- sd^2
alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
beta <- alpha * (1 / mu - 1)
return(params = list(alpha = alpha, beta = beta))
}
#create beta distribution for utilities, for successfull nosur and sur
alpha_beta <- estBetaParams(u_SUCC_NOSUR_uc, u_SUCC_NOSUR_uc_sd)
dist_u_SUCC_NOSUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_NOSUR_oc, u_SUCC_NOSUR_oc_sd)
dist_u_SUCC_NOSUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_uc, u_SUCC_SUR_uc_sd)
dist_u_SUCC_SUR_uc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
alpha_beta <- estBetaParams(u_SUCC_SUR_oc, u_SUCC_SUR_oc_sd)
dist_u_SUCC_SUR_oc <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#create gamma distributions for costs
dist_c_NOSUR_oc <- rgamma(100000, (c_NOSUR_oc_sd/c_NOSUR_oc)^-2, scale = c_NOSUR_oc * (c_NOSUR_oc_sd/c_NOSUR_oc)^2)
dist_c_NOSUR_uc <- rgamma(100000, (c_NOSUR_uc_sd/c_NOSUR_uc)^-2, scale = c_NOSUR_uc * (c_NOSUR_uc_sd/c_NOSUR_uc)^2)
dist_c_SUR_oc <- rgamma(100000, (c_SUR_oc_sd/c_SUR_oc)^-2, scale = c_SUR_oc * (c_SUR_oc_sd/c_SUR_oc)^2)
dist_c_SUR_uc <- rgamma(100000, (c_SUR_uc_sd/c_SUR_uc)^-2, scale = c_SUR_uc * (c_SUR_uc_sd/c_SUR_uc)^2)
dist_c_REV_oc <- rgamma(100000, (c_REV_oc_sd/c_REV_oc)^-2, scale = c_REV_oc * (c_REV_oc_sd/c_REV_oc)^2)
dist_c_REV_uc <- rgamma(100000, (c_REV_uc_sd/c_REV_uc)^-2, scale = c_REV_uc * (c_REV_oc_sd/c_REV_uc)^2)
#create gamma distributions for probabilities healthy to non-surgical treatment
#prob_PREV_m$value <- 0.02 #annual incidence for individuals aged 20+
#prob_PREV_f$value <- 0.02
#prob_prev_m_se <- 0.00512
#prop_prev_f_se <- 0.00512
alpha_beta <- estBetaParams(0.02, 0.00510)
dist_p_prev_m <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
dist_p_prev_f <- rbeta(100000, alpha_beta$alpha, alpha_beta$beta)
#define vectors for psa
psa_cost_disc_uc <- vector()
psa_cost_uc <- vector()
psa_cost_disc_oc <- vector()
psa_cost_oc <- vector()
psa_qaly_disc_oc <- vector()
psa_qaly_oc <- vector()
psa_qaly_disc_uc <- vector()
psa_qaly_uc <- vector()
res_cost_uc <- vector()
res_cost_oc5 <- vector()
res_qaly_uc <- vector()
res_qaly_oc5 <- vector()
discount <- 0.03
n_cycle <- 70
#run psa x times
for(j in 1:1000){
c_NOSUR_oc_psa <- sample(dist_c_NOSUR_oc,1)
c_NOSUR_uc_psa <- sample(dist_c_NOSUR_uc,1)
c_SUR_oc_psa <- sample(dist_c_SUR_oc,1)
c_SUR_uc_psa <- sample(dist_c_SUR_uc,1)
c_REV_oc_psa <- sample(dist_c_REV_oc,1)
c_REV_uc_psa <- sample(dist_c_REV_uc,1)
u_SUCC_NOSUR_oc_psa <- sample(dist_u_SUCC_NOSUR_oc,1)
u_SUCC_NOSUR_uc_psa <- sample(dist_u_SUCC_NOSUR_uc,1)
u_SUCC_SUR_oc_psa <- sample(dist_u_SUCC_SUR_oc,1)
u_SUCC_SUR_uc_psa <- sample(dist_u_SUCC_SUR_uc,1)
p_prev_m_psa <- sample(dist_p_prev_m,1)
p_prev_f_psa <- sample(dist_p_prev_f,1)
## RUN MODEL USUAL CARE
cycle <- 1:70
n_cycle <- length(cycle)
POP_fu <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fu[i+1] <- POP_fu[i] -
(POP_fu[i] * p_prev_f_psa) -
(POP_fu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mu <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mu[i+1] <- POP_mu[i] -
(POP_mu[i] * p_prev_m_psa) -
(POP_mu[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURu[i+1] <- (POP_fu[i]*p_prev_f_psa) +
(POP_mu[i]*p_prev_m_psa) -
(NOSURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURu[i+1] <- (SUCC_NOSURu[i] +NOSURu[i]) -
(SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SUCC_NOSURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURu[i+1] <- (SUCC_NOSURu[i]*prob_SUR_uc[prob_SUR_uc$age == age,2]) -
(SURu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURu[i+1] <- (SUCC_SURu[i] +SURu[i]) -
(SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(SUCC_SURu[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVu[i+1] <- (SUCC_SURu[i]*prob_REV_uc[prob_REV_uc$age == age,2]) -
(REVu[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVu[i+1] <- (SUCC_REVu[i] +REVu[i]) -
(SUCC_REVu[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHu <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHu[i+1] <- DEATHu[i] +
POP_fu[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mu[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURu[i] *bfs_mort[bfs_mort$age == age,2]+
SURu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURu[i] *bfs_mort[bfs_mort$age == age,2]+
REVu[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVu[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
## RUN MODEL OPTIMIZED CARE 2 year delay
cycle <- 1:70
n_cycle <- length(cycle)
POP_fo5 <- c(n_female, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_fo5[i+1] <- POP_fo5[i] -
(POP_fo5[i] * p_prev_f_psa) -
(POP_fo5[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
POP_mo5 <- c(n_male, rep(NA, n_cycle-1))
age = age_init
for(i in 1:(n_cycle-1)){
POP_mo5[i+1] <- POP_mo5[i] -
(POP_mo5[i] * p_prev_m_psa) -
(POP_mo5[i] * bfs_mort[bfs_mort$age == age,2])
age = age+1
}
NOSURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
NOSURo5[i+1] <- (POP_fo5[i]*p_prev_f_psa) +
(POP_mo5[i]*p_prev_m_psa) -
(NOSURo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_NOSURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_NOSURo5[i+1] <- (SUCC_NOSURo5[i] +NOSURo5[i])-
(SUCC_NOSURo5[i]*prob_SUR_oc[prob_SUR_oc$age == age,2])-
(SUCC_NOSURo5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del1_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del1_5[i+1] <- (SUCC_NOSURo5[i]*prob_SUR_oc[prob_SUR_oc$age == age,2]) -
(del1_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del2_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del2_5[i+1] <- (del1_5[i]) -
(del2_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del3_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del3_5[i+1] <- (del2_5[i]) -
(del3_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del4_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del4_5[i+1] <- (del3_5[i]) -
(del4_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
del5_5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
del5_5[i+1] <- (del4_5[i]) -
(del5_5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
SURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SURo5[i+1] <- (del5_5[i]) -
(SURo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_SURo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_SURo5[i+1] <- (SUCC_SURo5[i] +SURo5[i]) -
(SUCC_SURo5[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(SUCC_SURo5[i]*bfs_mort[bfs_mort$age == age,2])
age = age+1
}
REVo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
REVo5[i+1] <- (SUCC_SURo5[i]*prob_REV_oc[prob_REV_oc$age == age,2]) -
(REVo5[i] *bfs_mort[bfs_mort$age == age, 2])
age = age+1
}
SUCC_REVo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
SUCC_REVo5[i+1] <- (SUCC_REVo5[i] +REVo5[i]) -
(SUCC_REVo5[i]*bfs_mort[bfs_mort$age == age,2])
age=age+1
}
DEATHo5 <- c(0, rep(NA, n_cycle-1))
age=age_init
for(i in 1:(n_cycle-1)){
DEATHo5[i+1] <- DEATHo5[i] +
POP_fo5[i] *bfs_mort[bfs_mort$age == age,2]+
POP_mo5[i] *bfs_mort[bfs_mort$age == age,2]+
NOSURo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_NOSURo5[i] *bfs_mort[bfs_mort$age == age,2]+
del1_5[i] *bfs_mort[bfs_mort$age == age,2]+
del2_5[i] *bfs_mort[bfs_mort$age == age,2]+
del3_5[i] *bfs_mort[bfs_mort$age == age,2]+
del4_5[i] *bfs_mort[bfs_mort$age == age,2]+
del5_5[i] *bfs_mort[bfs_mort$age == age,2]+
SURo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_SURo5[i] *bfs_mort[bfs_mort$age == age,2]+
REVo5[i] *bfs_mort[bfs_mort$age == age,2]+
SUCC_REVo5[i] *bfs_mort[bfs_mort$age == age,2]
age=age+1
}
#cost calc for UNC
for(i in 1:n_cycle){
psa_cost_uc[i] <- NOSURu[i] *c_NOSUR_uc_psa+
SUCC_NOSURu[i] *c_SUCC_NOSUR_uc+
SURu[i] *c_SUR_uc_psa+
SUCC_SURu[i] *c_SUCC_SUR_uc+
REVu[i] *c_REV_uc_psa+
SUCC_REVu[i] *c_SUCC_REV_uc
}
#discouting costs UNC
for(i in 1:n_cycle){
psa_cost_disc_uc[i] <- psa_cost_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_uc <- sum(psa_cost_disc_uc)
psa_person_cost_uc <- psa_total_cost_uc/119672
res_cost_uc[j] <- psa_person_cost_uc
# cost calc for ONC
for(i in 1:n_cycle){
psa_cost_oc[i] <-
NOSURo5[i] *c_NOSUR_oc_psa+
SUCC_NOSURo5[i] *c_SUCC_NOSUR_oc+
del1_5[i] *c_SUCC_NOSUR_oc+
del2_5[i] *c_SUCC_NOSUR_oc+
del3_5[i] *c_SUCC_NOSUR_oc+
del4_5[i] *c_SUCC_NOSUR_oc+
del5_5[i] *c_SUCC_NOSUR_oc+
SURo5[i] *c_SUR_oc_psa+
SUCC_SURo5[i] *c_SUCC_SUR_oc+
REVo5[i] *c_REV_oc_psa+
SUCC_REVo5[i] *c_SUCC_REV_oc
}
#discounting costs ONC
for(i in 1:n_cycle){
psa_cost_disc_oc[i] <- psa_cost_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_cost_oc <- sum(psa_cost_disc_oc)
psa_person_cost_oc <- psa_total_cost_oc/119672
res_cost_oc5[j] <- psa_person_cost_oc
#utility calc for UNC
for(i in 1:n_cycle){
psa_qaly_uc[i] <- POP_fu[i] *u_NOSUR_uc+
POP_mu[i] *u_NOSUR_uc+
NOSURu[i] *u_NOSUR_uc+
SUCC_NOSURu[i] *u_SUCC_NOSUR_uc_psa+
SURu[i] *u_NOSUR_uc+
SUCC_SURu[i] *u_SUCC_SUR_uc_psa+
REVu[i] *u_NOSUR_uc+
SUCC_REVu[i] *u_SUCC_REV_uc
}
#discounting utility UNC
for(i in 1:n_cycle){
psa_qaly_disc_uc[i] <- psa_qaly_uc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_uc <- sum(psa_qaly_disc_uc)
psa_person_qaly_uc <- psa_total_qaly_uc/119672
res_qaly_uc[j] <- psa_person_qaly_uc
#utility calc for ONC
for(i in 1:n_cycle){
psa_qaly_oc[i] <- POP_fo5[i] *u_NOSUR_oc+
POP_mo5[i] *u_NOSUR_oc+
NOSURo5[i] *u_NOSUR_oc+
SUCC_NOSURo5[i] *u_SUCC_NOSUR_oc_psa+
del1_5[i] *u_SUCC_NOSUR_oc+
del2_5[i] *u_SUCC_NOSUR_oc+
del3_5[i] *u_SUCC_NOSUR_oc+
del4_5[i] *u_SUCC_NOSUR_oc+
del5_5[i] *u_SUCC_NOSUR_oc+
SURo5[i] *u_NOSUR_oc+
SUCC_SURo5[i] *u_SUCC_SUR_oc_psa+
REVo5[i] *u_NOSUR_oc+
SUCC_REVo5[i] *u_SUCC_REV_oc
}
#discounting utility ONC
for(i in 1:n_cycle){
psa_qaly_disc_oc[i] <- psa_qaly_oc[i]/(1+discount)^cycle_years[i]
}
psa_total_qaly_oc <- sum(psa_qaly_disc_oc)
psa_person_qaly_oc <- psa_total_qaly_oc/119672
res_qaly_oc5[j] <- psa_person_qaly_oc
}
mean(res_qaly_oc5)## [1] 17.02416mean(res_cost_oc5)## [1] 3912.628psa5_delta_costs <- res_cost_oc5 - res_cost_uc
psa5_delta_qalys <- res_qaly_oc5 - res_qaly_uc
kable(per_person_table, digits = 2)| scenario | cost_person | qaly_person |
|---|---|---|
| base case UNC | 4465.22 | 16.92 |
| base case ONC | 4805.83 | 17.07 |
| 2 year delay | 4451.21 | 17.06 |
| 5 year delay | 3963.80 | 17.03 |
CEAC
Cost-Effectiveness Acceptability Curve
Calculating and plotting cost-effectiveness acceptability curve using PSA output
## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----echo=TRUE, message=FALSE, warning=FALSE, include=TRUE-------------------------------------------------------------------------
wtp <- seq(0,100000,100)
ce_bc <- vector()
ce_2y <- vector()
ce_5y <- vector()
psa_icer_bc <- (res_cost_oc - res_cost_uc)/(res_qaly_oc - res_qaly_uc)
psa_e_diff_bc <- res_qaly_oc - res_qaly_uc
psa_c_diff_bc <- res_cost_oc - res_cost_uc
psa_icer_2y <- (res_cost_oc2 - res_cost_uc)/(res_qaly_oc2 - res_qaly_uc)
psa_e_diff_2y <- res_qaly_oc2 - res_qaly_uc
psa_c_diff_2y <- res_cost_oc2 - res_cost_uc
psa_icer_5y <- (res_cost_oc5 - res_cost_uc)/(res_qaly_oc5 - res_qaly_uc)
psa_e_diff_5y <- res_qaly_oc5 - res_qaly_uc
psa_c_diff_5y <- res_cost_oc5 - res_cost_uc
psa_icer_table <- data.frame(psa_icer_bc, psa_e_diff_bc,
psa_icer_2y, psa_e_diff_2y,
psa_icer_5y, psa_e_diff_5y)
## CAVE WURDE ANGEPASST AUF PSA, DESWEGEN /1000. AM SCHLUSS WIEDER AUF /10000
for(i in 1:1001){
ce_bc[i] <- sum(ifelse(psa_icer_bc<wtp[i] & psa_e_diff_bc >0, 1, 0))/10000
ce_2y[i] <- sum(ifelse(psa_icer_2y<wtp[i] & psa_e_diff_2y >0, 1, 0))/10000
ce_5y[i] <- sum(ifelse(psa_icer_5y<wtp[i] & psa_e_diff_5y >0, 1, 0))/10000
}
df <- data.frame(wtp_tres = c(wtp,wtp,wtp),
icer = c(ce_bc, ce_2y, ce_5y),
Scenario = c(rep("base case", length(ce_bc)),
rep("two-year delay", length(ce_2y)),
rep("five-year delay", length(ce_5y)))
)
df$Scenario <- factor(df$Scenario, levels = c("base case", "two-year delay", "five-year delay"))
ceac <- ggplot(df, aes(wtp_tres, icer, color = Scenario)) + geom_line(size = 1) +
labs(x = "Willingness to pay for one additional QALY", y = "Probability scenario is cost-effective \n compared to usual care")+
ylim(0, 1)+
xlim(0, 100000)+
theme_bw()
ceac
nw_bc <- sum(ifelse(psa_e_diff_bc<0 & psa_c_diff_bc > 0, 1, 0))/10000 #% of ICERs in NW
sw_bc <- sum(ifelse(psa_e_diff_bc<0 & psa_c_diff_bc < 0, 1, 0))/10000
se_bc <- sum(ifelse(psa_e_diff_bc>0 & psa_c_diff_bc < 0, 1, 0))/10000
ne_bc <- sum(ifelse(psa_e_diff_bc>0 & psa_c_diff_bc > 0, 1, 0))/10000
nw_2y <- sum(ifelse(psa_e_diff_2y<0 & psa_c_diff_2y > 0, 1, 0))/10000 #% of ICERs in NW
sw_2y <- sum(ifelse(psa_e_diff_2y<0 & psa_c_diff_2y < 0, 1, 0))/10000
se_2y <- sum(ifelse(psa_e_diff_2y>0 & psa_c_diff_2y < 0, 1, 0))/10000
ne_2y <- sum(ifelse(psa_e_diff_2y>0 & psa_c_diff_2y > 0, 1, 0))/10000
nw_5y <- sum(ifelse(psa_e_diff_5y<0 & psa_c_diff_5y > 0, 1, 0))/10000 #% of ICERs in NW
sw_5y <- sum(ifelse(psa_e_diff_5y<0 & psa_c_diff_5y < 0, 1, 0))/10000
se_5y <- sum(ifelse(psa_e_diff_5y>0 & psa_c_diff_5y < 0, 1, 0))/10000
ne_5y <- sum(ifelse(psa_e_diff_5y>0 & psa_c_diff_5y > 0, 1, 0))/10000
nw_bc + se_bc + sw_bc + ne_bc # Check if 1## [1] 0.1nw_2y + se_2y + sw_2y + ne_2y## [1] 0.1nw_5y + se_5y + sw_5y + ne_5y## [1] 0.1se_bc## [1] 0.0174ne_bc## [1] 0.0609se_2y## [1] 0.0261ne_2y## [1] 0.048se_5y## [1] 0.0798ne_5y## [1] 0Plots
Plots used for publication
## ----setup, include=FALSE----------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
## ----fig.cap= "counts per state for base case scenario", fig.align='center', fig.height=6, fig.width=8, message= FALSE, warning= FALSE------------------
model_base <- data.frame( n= c(SUCC_NOSURu, SUCC_SURu, SUCC_REVu, POP_fu, POP_mu, DEATHu),
cycle = 1:n_cycle,
State = c(rep("SUCC_NOSUR", length(SUCC_NOSURu)),
rep("SUCC_SUR", length(SUCC_SURu)),
rep("SUCC_REV", length(SUCC_REVu)),
rep("POP_fu", length(POP_fu)),
rep("POP_mu", length(POP_mu)),
rep("DEATH", length(DEATHu))
))
plot_base <- ggplot(model_base, aes(x=cycle , y= n, color = State)) +
geom_line(size =1) +
ylim(0, 120000)+
ylab("count") +
# ggtitle("counts per state for base case Scenario")+
theme_bw()
plot_base
## ----fig.cap= "counts per state for base case scenario", fig.align='center', fig.height=6, fig.width=8, message= FALSE, warning= FALSE------------------
model_base_nodeath <- data.frame( n= c(SUCC_NOSURu, SUCC_SURu, SUCC_REVu),
cycle = 1:n_cycle,
State = c(rep("SUCC_NOSUR", length(SUCC_NOSURu)),
rep("SUCC_SUR", length(SUCC_SURu)),
rep("SUCC_REV", length(SUCC_REVu))
))
plot_base_nodeath <- ggplot(model_base_nodeath, aes(x=cycle , y= n, color = State)) +
geom_line(size =1) +
ylim(0, 40000)+
ylab("count") +
# ggtitle("counts per state for base case Scenario")+
theme_bw()
plot_base_nodeath
## ----fig.cap= "", fig.align='center', fig.height=6, fig.width=8, message= FALSE, warning= FALSE------------------
comp_SUR <- data.frame(n= c(SUCC_SURu, SUCC_SURo2, SUCC_SURo5),
cycle = 1:n_cycle,
Scenario = c(rep("base case", length(SUCC_SURu)),
rep("two-year delay", length(SUCC_SURo2)),
rep("five-year delay", length(SUCC_SURo5))
))
plot_comp_SUR <- ggplot(comp_SUR, aes(x=cycle , y= n, color = Scenario)) +
geom_line(size =1) +
ylim(0, 30000)+
ylab("count") +
ggtitle("TKR")+
theme_bw()
plot_comp_SUR
count_sur <- data.frame(Scenario = c("base case", "two-year delay", "five-year delay"),
number = c(sum(SURu),sum(SURo2),sum(SURo5)),
change = c((sum(SURu)/sum(SURu)-1), (sum(SURo2)/sum(SURu)-1), (sum(SURo5)/sum(SURu)-1))
)
kable(count_sur, digits = 2)| Scenario | number | change |
|---|---|---|
| base case | 34275.41 | 0.00 |
| two-year delay | 35887.02 | 0.05 |
| five-year delay | 28932.45 | -0.16 |
## ----fig.cap= "", fig.align='center', fig.height=6, fig.width=8, message= FALSE, warning= FALSE------------------
comp_REV <- data.frame(n= c(SUCC_REVu, SUCC_REVo2, SUCC_REVo5),
cycle = 1:n_cycle,
Scenario = c(rep("base case", length(SUCC_REVu)),
rep("two-year delay", length(SUCC_REVo2)),
rep("five-year delay", length(SUCC_REVo5))
))
comp_REV$Scenario <- factor(comp_REV$Scenario, levels = c("base case", "two-year delay", "five-year delay"))
plot_comp_REV <- ggplot(comp_REV, aes(x=cycle , y= n, color = Scenario)) +
geom_line(size =1) +
ylim(0, 2500)+
ylab("Revision Surgery") +
xlab("Model Cycle")+
#ggtitle("Revisons")+
theme_bw()
plot_comp_REV 
count_rev <- data.frame(Scenario = c("base case", "two-year delay", "five-year delay"),
number = c(sum(REVu),sum(REVo2),sum(REVo5)),
change = c((sum(REVu)/sum(REVu)-1), (sum(REVo2)/sum(REVu)-1), (sum(REVo5)/sum(REVu)-1))
)
kable(count_rev, digits = 2)| Scenario | number | change |
|---|---|---|
| base case | 2816.59 | 0.00 |
| two-year delay | 2672.45 | -0.05 |
| five-year delay | 1813.05 | -0.36 |
## ----fig.cap= "Model Outputs", fig.align='center', fig.height=8, fig.width=10, message= FALSE, warning= FALSE----
cepanel <- ggplot(icer_table, aes(x=qaly_delta , y= cost_delta, color = scenario)) +
geom_point(size =2) +
ylim(-2000, 2000)+
ylab("Incremental Cost (CHF)") +
xlim(-1.0,1.0)+
xlab("Incremental QALY") +
# ggtitle("ICER")+
geom_hline(yintercept=0)+
geom_vline(xintercept = 0)+
theme_bw()
cepanel 
## ----include= TRUE, echo= TRUE, message= FALSE, warning= FALSE---------------------------------------------------
# calculate delta costs and delta qalys - plot results
# psabc_delta_costs <- res_cost_oc - res_cost_uc
# psabc_delta_qalys <- res_qaly_oc - res_qaly_uc
# psa2_delta_costs <- res_cost_oc2 - res_cost_uc
# psa2_delta_qalys <- res_qaly_oc2 - res_qaly_uc
# psa5_delta_costs <- res_cost_oc5 - res_cost_uc
# psa5_delta_qalys <- res_qaly_oc5 - res_qaly_uc
psa_table <- data.frame(qaly_delta = c(psabc_delta_qalys, psa2_delta_qalys, psa5_delta_qalys),
cost_delta = c(psabc_delta_costs, psa2_delta_costs, psa5_delta_costs),
icer = c(round(psabc_delta_costs/psabc_delta_qalys,0),
round(psa2_delta_costs/psa2_delta_qalys,0),
round(psa5_delta_costs/psa5_delta_qalys,0)),
Scenario = c(rep("base case", length(psabc_delta_costs)),
rep("two-year delay", length(psabc_delta_costs)),
rep("five-year delay", length(psabc_delta_costs))
))
#amount of ICER in upper east quadrant and the lower east quadrant, therefore dominating the basecase
psa_ue <- subset(psa_table, qaly_delta > 0)
#calculate how many of the psa ICERS are in dominant
nrow(subset(psa_ue, Scenario == "base case")) / nrow(subset(psa_table, Scenario == "base case"))*100## [1] 86.5nrow(subset(psa_ue, Scenario == "two-year delay")) / nrow(subset(psa_table, Scenario =="two-year delay"))*100## [1] 82.8nrow(subset(psa_ue, Scenario == "five-year delay")) / nrow(subset(psa_table, Scenario =="five-year delay"))*100## [1] 79.8psa_table$Scenario <- factor(psa_table$Scenario, levels = c("base case", "two-year delay", "five-year delay"))
psapanel <- ggplot(psa_table, aes(x=qaly_delta , y= cost_delta, color= Scenario) ) +
geom_point(size =0.2) +
geom_point(aes(x= delta_qalys_base, y= delta_costs_base), color = "red")+
geom_point(aes(x= delta_qalys_scen5, y= delta_costs_scen5), color = "blue")+
geom_point(aes(x= delta_qalys_scen2, y= delta_costs_scen2), color = "green")+
ylim(-1000, 1000)+
ylab("Incremental Cost (CHF)") +
xlim(-1,1)+
xlab("Incremental QALY") +
#ggtitle("ICER PSA")+
geom_hline(yintercept=0)+
geom_vline(xintercept = 0)+
guides(colour = guide_legend(override.aes = list(size=2)))+
# geom_mark_ellipse(expand = 0,aes(fill=icer ), color = "black", show.legend = FALSE, tol = 2)+
# annotate("text", x=0.7, y=1800, label= paste(psa_in_ue , "% in ue quadrant"))+
theme_bw()
psapanel
#stacked bar plot for cost overview
data_stack <- read.csv("C:/Users/schwa/OneDrive/99_REVISION/revision/bar_chart.csv", header = TRUE, sep = ";", check.names = FALSE)
data_stack$scenario <- factor(data_stack$scenario, # Change ordering manually
levels = c("CMOC", "OMOC base case", "OMOC two-year delay", "OMOC five-year delay"))
data_stack$state <- factor(data_stack$state,
levels = c("Non-surgical treatment", "Successful non-surgical", "TKR surgery", "Successful TKR", "Revision surgery", "Successful revision"))
stack_plot <- ggplot(data_stack, aes(fill=state, y=value, x=scenario)) +
geom_bar(position="stack", stat="identity")+
ylab("per-person cost [CHF]")+
xlab("model of care")+
scale_fill_brewer(palette="Spectral")+
theme_bw()
stack_plot
ggsave(plot = stack_plot, filename = "C:/Users/schwa/OneDrive/99_REVISION/revision/cost_by_state_compr.tiff", width = 10, height = 6, compression = "lzw" )
#plot ceac curve
# ceac
#save all plots
# setwd("C:/Users/schwa/Desktop/Revision/SMW/revision")
#
ggsave(plot = psapanel, filename = "C:/Users/schwa/OneDrive/99_REVISION/revision/psa_28122022.tiff", width = 6, height = 4 )
ggsave(plot = ceac, filename = "C:/Users/schwa/OneDrive/99_REVISION/revision/ceac_28122022.tiff", width = 6, height = 4 )
# ggsave(plot = cepanel, filename = "cepanel.tiff", width = 6, height = 4 )
# ggsave(plot = plot_base, filename = "plot_base.tiff", width = 6, height = 4 )
# ggsave(plot = plot_comp_SUR, filename = "plot_comp_SUR.tiff", width = 6, height = 4 )
ggsave(plot = plot_comp_REV, filename = "C:/Users/schwa/Desktop/markov_ONC_02122022/markov_ONC/plot_comp_REV_new.tiff", width = 6, height = 4 )
# ggsave(plot = plot_base_nodeath, filename = "counts_ppp.tiff", width = 6, height = 4 )
#
# ggsave(plot = psapanel, filename = "psapanel.jpeg", width = 6, height = 4 )
# ggsave(plot = ceac, filename = "ceac.jpeg", width = 6, height = 4 )
# ggsave(plot = cepanel, filename = "cepanel.jpeg", width = 6, height = 4 )
# ggsave(plot = plot_base, filename = "plot_base.jpeg", width = 6, height = 4 )
# ggsave(plot = plot_comp_SUR, filename = "plot_comp_SUR.jpeg", width = 6, height = 4 )
# ggsave(plot = plot_comp_REV, filename = "plot_comp_REV.jpeg", width = 6, height = 4 )
# ggsave(plot = plot_base_nodeath, filename = "counts_ppp.jpeg", width = 6, height = 4 )