-
Notifications
You must be signed in to change notification settings - Fork 0
/
08_CEA_1yr_MC_CostsUtility.Rmd
328 lines (266 loc) · 14.2 KB
/
08_CEA_1yr_MC_CostsUtility.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
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
---
title: "CEA anticoag - MC all"
author: "Emily O'Neill, Andrew Huang"
date: "2024-01-17"
output: html_document
---
MC draws for event rates, costs and utilities
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 static values
```{r}
### Hyperparameters
cycle_length <- 1/12 # cycle length equals 1 month
n.t <- 12 # number of cycles, time horizon equals 1 year
d.c <- d.e <- 0.05 / (1/cycle_length) # discounting of costs and QALYs, annual rate = 5%
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
#Number of MC Samples
n_samples <- 1000
#Set the random number seed for reproducibility
set.seed(123)
```
Set MC draws for event rate parameters
```{r}
##Rivaroxaban
#Health States
#Event rates per person-6months converted to person-year:
sample_r.ne_vr <- runif(n_samples, min = 0, max = 0.000 + 0.50) * 2 # rate of no event to VTE recurrence
sample_r.ne_mb <- runif(n_samples, min = 0, max = 0.019 + 0.05) * 2 # rate of no event to major bleed
sample_r.ne_nmb <- runif(n_samples, min = 0, max = 0.125 + 0.05) * 2 # rate of no event to CRNMB
sample_r.nmb_ot <- runif(n_samples, min = 0, max = 0.20) * 2 # rate of CNMB to come off treatment
sample_r.mb_ot <- runif(n_samples, min = 0, max = 0.20) * 2 # rate of MB to come off treatment
#Transition probabilities (per month cycle):
sample_pr.ne_vr <- 1 - exp(-sample_r.ne_vr * cycle_length) # from no event to VTE recurrence
sample_pr.ne_mb <- 1 - exp(-sample_r.ne_mb * cycle_length) # from no event to major bleed
sample_pr.ne_nmb <- 1 - exp(-sample_r.ne_nmb * cycle_length) # from no event to CRNMB
sample_pr.nmb_ot <- 1 - exp(-sample_r.nmb_ot * cycle_length) # from CNMB to come off treatment
sample_pr.mb_ot <- 1 - exp(-sample_r.mb_ot * cycle_length) # from MB to come off treatment
# Create dataframe of lists of parameter values for state transition probabilities
l.probr <- list()
for (i in 1:n_samples){
sample.l.probr <- list(pr.ne_vr = sample_pr.ne_vr[i],
pr.ne_mb = sample_pr.ne_mb[i],
pr.ne_nmb = sample_pr.ne_nmb[i],
pr.ne_ne = 1 - sum(sample_pr.ne_vr[i],
sample_pr.ne_mb[i],
sample_pr.ne_nmb[i]),
pr.nmb_ot = sample_pr.nmb_ot[i],
pr.mb_ot = sample_pr.mb_ot[i])
l.probr[i] <- list(sample.l.probr)
}
df.probr_riv <- as.data.frame(do.call(rbind, l.probr))
```
```{r}
##Warfarin
#Health States
#Event rates per person-6months converted to person-year:
sample_r.ne_vr <- runif(n_samples, min = 0, max = 0.026 + 0.50) * 2 # rate of no event to VTE recurrence
sample_r.ne_mb <- runif(n_samples, min = 0, max = 0.090 + 0.05) * 2 # rate of no event to major bleed
sample_r.ne_nmb <- runif(n_samples, min = 0, max = 0.308 + 0.05) * 2 # rate of no event to CRNMB
sample_r.nmb_ot <- runif(n_samples, min = 0, max = 0.20) * 2 # rate of CNMB to come off treatment
sample_r.mb_ot <- runif(n_samples, min = 0, max = 0.20) * 2 # rate of MB to come off treatment
#Transition probabilities (per month cycle):
sample_pr.ne_vr <- 1 - exp(-sample_r.ne_vr * cycle_length) # from no event to VTE recurrence
sample_pr.ne_mb <- 1 - exp(-sample_r.ne_mb * cycle_length) # from no event to major bleed
sample_pr.ne_nmb <- 1 - exp(-sample_r.ne_nmb * cycle_length) # from no event to CRNMB
sample_pr.nmb_ot <- 1 - exp(-sample_r.nmb_ot * cycle_length) # from CNMB to come off treatment
sample_pr.mb_ot <- 1 - exp(-sample_r.mb_ot * cycle_length) # from MB to come off treatment
# Create dataframe of lists of parameter values for state transition probabilities
l.probr <- list()
for (i in 1:n_samples){
sample.l.probr <- list(pr.ne_vr = sample_pr.ne_vr[i],
pr.ne_mb = sample_pr.ne_mb[i],
pr.ne_nmb = sample_pr.ne_nmb[i],
pr.ne_ne = 1 - sum(sample_pr.ne_vr[i],
sample_pr.ne_mb[i],
sample_pr.ne_nmb[i]),
pr.nmb_ot = sample_pr.nmb_ot[i],
pr.mb_ot = sample_pr.mb_ot[i])
l.probr[i] <- list(sample.l.probr)
}
df.probr_war <- as.data.frame(do.call(rbind, l.probr))
```
Set MC draws for cost parameters: gamma distributions
```{r}
#' @param v.cost is a numeric vector containing cost values for each state
##Rivaroxaban Costs (per 6 months (year/2))
riv.c.ne <- rgamma(n_samples, shape=10, rate=1/14) / 12
#Cost for staying healthy in rivaroxaban arm
riv.c.vr <- rgamma(n_samples, shape=10, rate=1/21.5) / 12
#Cost for VTE recurrence in Rivaroxaban arm ($215 annually)
riv.c.mb <- rgamma(n_samples, shape=2, rate= 1/245) / 12
#Cost for major bleed in Rivaroxaban arm ($490 annually)
riv.c.nmb <- rgamma(n_samples, shape=2, rate= 1/152) / 12
#Cost for CRNMB in Rivaroxaban arm ($304 annually)
riv.c.ot <- rgamma(n_samples, shape=1, rate=1/107.5) / 12 #Cost of off treatment ($0 annually)
##Warfarin Costs (per 6 months (year/2))
war.c.ne <- rgamma(n_samples, shape=10, rate=1/16.5) / 12
#Cost for staying healthy in Warfarin arm ($165 annually)
war.c.vr <- rgamma(n_samples, shape=10, rate=1/26.7) / 12
#Cost for VTE recurrence in Warfarin arm ($267 annually)
war.c.mb <- rgamma(n_samples, shape=3, rate= 1/217.5) / 12
#Cost for major bleed in Warfarin arm ($435 annually)
war.c.nmb <- rgamma(n_samples, shape=3, rate= 1/166) / 12
#Cost for CRNMB in Warfarin arm ($332 annually)
war.c.ot <- rgamma(n_samples, shape=1, rate=1/133.5) / 12 #Cost of off treatment ($0 annually)
##Vectors for cost values
riv.v.cost <- cbind(riv.c.ne, riv.c.vr, riv.c.mb, riv.c.nmb, riv.c.ot)
war.v.cost <- cbind(war.c.ne, war.c.vr, war.c.mb, war.c.nmb, war.c.ot)
l.costr <- list()
l.costw <- list()
for (i in 1:n_samples){
l.costr[[i]] = unlist(riv.v.cost[i, ], use.names=F)
l.costw[[i]] = unlist(war.v.cost[i, ], use.names=F)
}
```
Set MC draws for utility parameters: logit distributions
```{r}
#' @param v.util is a numeric vector containing utility values for each state
##Rivaroxaban Values
riv.u.ne <- rlogis(n_samples, location=0.825, scale=0.02) #Utility of No Event
riv.u.vr <- rlogis(n_samples, location=0.76, scale=0.063774925) #Utility of recurrence
riv.u.mb <- rlogis(n_samples, location=0.61, scale=0.274593292) #Utility of major bleed
riv.u.nmb <- rlogis(n_samples, location=0.65, scale=0.115559069) #Utility of CRNMB
riv.u.ot <- rlogis(n_samples, location=0.68, scale=0.063774925) #Utility of off treatment
##Warfarin Values
war.u.ne <- rlogis(n_samples, location=0.825, scale=0.02) #Utility of No Event
war.u.vr <- rlogis(n_samples, location=0.76, scale=0.063774925) #Utility of recurrence
war.u.mb <- rlogis(n_samples, location=0.61, scale=0.274593292) #Utility of major bleed
war.u.nmb <- rlogis(n_samples, location=0.65, scale=0.115559069) #Utility of CRNMB
war.u.ot <- rlogis(n_samples, location=0.68, scale=0.063774925) #Utility of off treatment
##Vectors for utility values
riv.v.util <- cbind(riv.u.ne, riv.u.vr, riv.u.mb, riv.u.nmb, riv.u.ot)
war.v.util <- cbind(war.u.ne, war.u.vr, war.u.mb, war.u.nmb, war.u.ot)
l.utilr <- list()
l.utilw <- list()
for (i in 1:n_samples){
l.utilr[[i]] = unlist(riv.v.util[i, ], use.names=F)
l.utilw[[i]] = unlist(war.v.util[i, ], use.names=F)
}
```
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(rep(d.c, n_samples)),
d.e = as.numeric(rep(d.e, n_samples)),
v.n = rep(list(v.n), n_samples),
v.cost = l.costr, ### SENS
v.util = l.utilr, ### SENS
v.M_1 = rep(list(v.M_1), n_samples),
pr.ne_ne = df.probr_riv$pr.ne_ne,
pr.ne_vr = df.probr_riv$pr.ne_vr,
pr.ne_mb = df.probr_riv$pr.ne_mb,
pr.ne_nmb = df.probr_riv$pr.ne_nmb,
pr.mb_ot = df.probr_riv$pr.mb_ot,
pr.nmb_ot = df.probr_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(rep(d.c, n_samples)),
d.e = as.numeric(rep(d.e, n_samples)),
v.n = rep(list(v.n), n_samples),
v.cost = l.costw, ### SENS
v.util = l.utilw, ### SENS
v.M_1 = rep(list(v.M_1), n_samples),
pr.ne_ne = df.probr_war$pr.ne_ne,
pr.ne_vr = df.probr_war$pr.ne_vr,
pr.ne_mb = df.probr_war$pr.ne_mb,
pr.ne_nmb = df.probr_war$pr.ne_nmb,
pr.mb_ot = df.probr_war$pr.mb_ot,
pr.nmb_ot = df.probr_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_08_MC_CostsUtility.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")
```
Worst CERs (highest cost / lowest effectiveness)
```{r}
data.frame(
Strategy = c("Warfarin", "Rivaroxaban"),
max_cost = c(max(psa$cost$Warfarin), max(psa$cost$Rivaroxaban)),
min_effectiveness = c(min(psa$effectiveness$Warfarin), min(psa$effectiveness$Rivaroxaban)),
worst_CER = c(max(psa$cost$Warfarin)/min(psa$effectiveness$Warfarin),
max(psa$cost$Rivaroxaban)/min(psa$effectiveness$Rivaroxaban))
)
```