-
Notifications
You must be signed in to change notification settings - Fork 0
/
04a_CEA_1yr_Sensitivity_DiscountRate.Rmd
240 lines (198 loc) · 11.1 KB
/
04a_CEA_1yr_Sensitivity_DiscountRate.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
---
title: "Sensitivity - discount rates"
author: "Emily O'Neill, Andrew Huang"
date: "2023-08-17"
output: html_document
---
Testing ICER sensitivity to discount rates: 3%, 5% (base), 7%
Setup
```{r}
##Cost utility analysis of Rivaroxaban in the prevention of VTE reoccurrence versus standard or care with warfarin: A markov Model-Based Simulation
##remove any variables in R's Memory
rm(list = ls())
##Install packages for Markov Model
library(here)
library(tidyverse)
# Load functions
source(here::here("functions", "function_markov.R"))
source(here::here("functions", "function_plot-trace.R"))
source(here::here("functions", "function_create-sa.R"))
source(here::here("functions", "function_make-psa-obj.R"))
source(here::here("functions", "function_summary-psa.R"))
source(here::here("functions", "function_plot-psa.R"))
```
Set Parameters with Base Case values
```{r}
### Hyperparameters
cycle_length <- 1/12 # cycle length equals 1 month
n.t <- 12 # number of cycles, time horizon equals 1 year
v.n <- c("NE", "VR", "MB", "NMB", "OT") # Model states
v.M_1 <- c(1, 0, 0, 0, 0) # Everyone begins in no event state
#' @param v.cost is a numeric vector containing cost values for each state
##Rivaroxaban Costs (per 6 months (year/2))
riv.c.ne <- 142 / 12 #Cost for staying healthy in rivaroxaban arm
riv.c.vr <- 215 / 12 #Cost for VTE recurrence in Rivaroxaban arm ($215 annually)
riv.c.mb <- 490 / 12 #Cost for major bleed in Rivaroxaban arm ($490 annually)
riv.c.nmb <- 304 / 12 #Cost for CRNMB in Rivaroxaban arm ($304 annually)
riv.c.ot <- 0 / 12 #Cost of off treatment ($0 annually)
##Warfarin Costs (per 6 months (year/2))
war.c.ne <- 165 / 12 #Cost for staying healthy in Warfarin arm ($165 annually)
war.c.vr <- 267 / 12 #Cost for VTE recurrence in Warfarin arm ($267 annually)
war.c.mb <- 435 / 12 #Cost for major bleed in Warfarin arm ($435 annually)
war.c.nmb <- 332 / 12 #Cost for CRNMB in Warfarin arm ($332 annually)
war.c.ot <- 0 / 12 #Cost of off treatment ($0 annually)
##Vectors for cost values
riv.v.cost <- c(riv.c.ne, riv.c.vr, riv.c.mb, riv.c.nmb, riv.c.ot)
war.v.cost <- c(war.c.ne, war.c.vr, war.c.mb, war.c.nmb, war.c.ot)
#' @param v.util is a numeric vector containing utility values for each state
##Rivaroxaban Values
riv.u.ne <- 0.825 #Utility of No Event
riv.u.vr <- 0.76 #Utility of recurrence
riv.u.mb <- 0.61 #Utility of major bleed
riv.u.nmb <- 0.65 #Utility of CRNMB
riv.u.ot <- 0.68 #Utility of off treatment
##Warfarin Values
war.u.ne <- 0.825 #Utility of No Event
war.u.vr <- 0.76 #Utility of recurrence
war.u.mb <- 0.61 #Utility of major bleed
war.u.nmb <- 0.65 #Utility of CRNMB
war.u.ot <- 0.68 #Utility of off treatment
##Vectors for utility values
riv.v.util <- c(riv.u.ne, riv.u.vr, riv.u.mb, riv.u.nmb, riv.u.ot)
war.v.util <- c(war.u.ne, war.u.vr, war.u.mb, war.u.nmb, war.u.ot)
#' @param pr._ represent parameter values for state transition probabilities
## Rivaroxaban Values
#Event rates per person-6months converted to person-year:
riv.r.ne_vr <- 0.000 * 2 # rate of no event to VTE recurrence
riv.r.ne_mb <- 0.019 * 2 # rate of no event to major bleed
riv.r.ne_nmb <- 0.125 * 2 # rate of no event to CRNMB
riv.r.nmb_ot <- 0.100 * 2 # rate of CNMB to come off treatment
riv.r.mb_ot <- 0.100 * 2 # rate of MB to come off treatment
#Transition probabilities (per month cycle):
riv.pr.ne_vr <- 1 - exp(-riv.r.ne_vr * cycle_length) # from no event to VTE recurrence
riv.pr.ne_mb <- 1 - exp(-riv.r.ne_mb * cycle_length) # from no event to major bleed
riv.pr.ne_nmb <- 1 - exp(-riv.r.ne_nmb * cycle_length) # from no event to CRNMB
riv.pr.ne_ne <- 1 - sum(riv.pr.ne_vr, riv.pr.ne_mb, riv.pr.ne_nmb) # from no event to no event
riv.pr.nmb_ot <- 1 - exp(-riv.r.nmb_ot * cycle_length) # from CNMB to come off treatment
riv.pr.mb_ot <- 1 - exp(-riv.r.mb_ot * cycle_length) # from MB to come off treatment
## Warfarin Values
#Event rates per person-6months converted to person-year:
war.r.ne_vr <- 0.026 * 2 # rate of no event to VTE recurrence
war.r.ne_mb <- 0.008 * 2 # rate of no event to major bleed
war.r.ne_nmb <- 0.308 * 2 # rate of no event to CRNMB
war.r.nmb_ot <- 0.100 * 2 # rate of CNMB to come off treatment
war.r.mb_ot <- 0.100 * 2 # rate of MB to come off treatment
#Transition probabilities (per month cycle):
war.pr.ne_vr <- 1 - exp(-war.r.ne_vr * cycle_length) # from no event to VTE recurrence
war.pr.ne_mb <- 1 - exp(-war.r.ne_mb * cycle_length) # from no event to major bleed
war.pr.ne_nmb <- 1 - exp(-war.r.ne_nmb * cycle_length) # from no event to CRNMB
war.pr.ne_ne <- 1 - sum(war.pr.ne_vr, war.pr.ne_mb, war.pr.ne_nmb) # from no event to no event
war.pr.nmb_ot <- 1 - exp(-war.r.nmb_ot * cycle_length) # from CNMB to come off treatment
war.pr.mb_ot <- 1 - exp(-war.r.mb_ot * cycle_length) # from MB to come off treatment
```
Set Parameter with Sensitivity values
```{r}
### Parameter with Sensitivity values
d.c <- d.e <- c(0.03, 0.05, 0.07)
d.c <- d.e <- d.e / (1/cycle_length)
n_samples <- length(d.c)
```
Generate params
```{r}
## Generate an input dataframe containing parameters for each simulation
riv.input_params <- as.data.frame(cbind(nsim = as.character(seq(1, n_samples)),
n.t = as.numeric(rep(n.t, n_samples)),
d.c = as.numeric(d.c), ### SENS
d.e = as.numeric(d.e), ### SENS
v.n = rep(list(v.n), n_samples),
v.cost = rep(list(riv.v.cost), n_samples),
v.util = rep(list(riv.v.util), n_samples),
v.M_1 = rep(list(v.M_1), n_samples),
pr.ne_ne = as.numeric(riv.pr.ne_ne),
pr.ne_vr = as.numeric(riv.pr.ne_vr),
pr.ne_mb = as.numeric(riv.pr.ne_mb),
pr.ne_nmb = as.numeric(riv.pr.ne_nmb),
pr.mb_ot = as.numeric(riv.pr.mb_ot),
pr.nmb_ot = as.numeric(riv.pr.nmb_ot)))
war.input_params <- as.data.frame(cbind(nsim = as.character(seq(1, n_samples)),
n.t = as.numeric(rep(n.t, n_samples)),
d.c = as.numeric(d.c), ### SENS
d.e = as.numeric(d.e), ### SENS
v.n = rep(list(v.n), n_samples),
v.cost = rep(list(war.v.cost), n_samples),
v.util = rep(list(war.v.util), n_samples),
v.M_1 = rep(list(v.M_1), n_samples),
pr.ne_ne = as.numeric(war.pr.ne_ne),
pr.ne_vr = as.numeric(war.pr.ne_vr),
pr.ne_mb = as.numeric(war.pr.ne_mb),
pr.ne_nmb = as.numeric(war.pr.ne_nmb),
pr.mb_ot = as.numeric(war.pr.mb_ot),
pr.nmb_ot = as.numeric(war.pr.nmb_ot)))
```
Run Markov models on sensitivity values
```{r}
### Riv
#Initialize dataframe to store results
sim_results_riv <- data.frame(nsim = as.character(seq(1:n_samples)))
for (i in 1:n_samples){
sim_markov_riv <- Markov(n.t = riv.input_params$n.t[[i]],
d.c = riv.input_params$d.c[[i]],
d.e = riv.input_params$d.e[[i]],
v.n = riv.input_params$v.n[[i]],
v.cost = riv.input_params$v.cost[[i]],
v.util = riv.input_params$v.util[[i]],
v.M_1 = riv.input_params$v.M_1[[i]],
pr.ne_ne = riv.input_params$pr.ne_ne[[i]],
pr.ne_vr = riv.input_params$pr.ne_vr[[i]],
pr.ne_mb = riv.input_params$pr.ne_mb[[i]],
pr.ne_nmb = riv.input_params$pr.ne_nmb[[i]],
pr.mb_ot = riv.input_params$pr.mb_ot[[i]],
pr.nmb_ot = riv.input_params$pr.nmb_ot[[i]])
sim_results_riv$trace[i] <- list(sim_markov_riv$m.TRr)
sim_results_riv$tc_hat[i] <- sim_markov_riv$tc_hat
sim_results_riv$te_hat[i] <- sim_markov_riv$te_hat
}
### War
#Initialize dataframe to store results
sim_results_war <- data.frame(nsim = as.character(seq(1:n_samples)))
for (i in 1:n_samples){
sim_markov_war <- Markov(n.t = war.input_params$n.t[[i]],
d.c = war.input_params$d.c[[i]],
d.e = war.input_params$d.e[[i]],
v.n = war.input_params$v.n[[i]],
v.cost = war.input_params$v.cost[[i]],
v.util = war.input_params$v.util[[i]],
v.M_1 = war.input_params$v.M_1[[i]],
pr.ne_ne = war.input_params$pr.ne_ne[[i]],
pr.ne_vr = war.input_params$pr.ne_vr[[i]],
pr.ne_mb = war.input_params$pr.ne_mb[[i]],
pr.ne_nmb = war.input_params$pr.ne_nmb[[i]],
pr.mb_ot = war.input_params$pr.mb_ot[[i]],
pr.nmb_ot = war.input_params$pr.nmb_ot[[i]])
sim_results_war$trace[i] <- list(sim_markov_war$m.TRr)
sim_results_war$tc_hat[i] <- sim_markov_war$tc_hat
sim_results_war$te_hat[i] <- sim_markov_war$te_hat
}
```
Calculate ICERs
```{r}
sim_results_cost <- data.frame(war = sim_results_war$tc_hat,
riv = sim_results_riv$tc_hat)
sim_results_effectiveness <- data.frame(war = sim_results_war$te_hat,
riv = sim_results_riv$te_hat)
psa <- make_psa_obj(cost = sim_results_cost, effectiveness = sim_results_effectiveness, parameters = NULL,
strategies = c("Warfarin", "Rivaroxaban"), currency = "$", other_outcome = NULL)
# Save psa object
saveRDS(psa, here::here("output","results_1yr_sens_04a_DiscountRate.rds"))
# See mean ICER
summary.psa(psa, calc_sds=TRUE)
```
Distribution of ICER estimates
```{r}
plot.psa(psa) +
ggthemes::scale_color_colorblind() +
ggthemes::scale_fill_colorblind() +
xlab("Effectiveness (QALYs)") +
guides(col = guide_legend(nrow = 2)) +
theme(legend.position = "right")
```