-
Notifications
You must be signed in to change notification settings - Fork 0
/
04c_CEA_1yr_Sensitivity_Utility.Rmd
491 lines (433 loc) · 19.2 KB
/
04c_CEA_1yr_Sensitivity_Utility.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
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
---
title: "Sensitivity - utilities"
author: "Emily O'Neill, Andrew Huang"
date: "2023-08-21"
output: html_document
---
Testing ICER sensitivity to Utility ranges
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
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
#' @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 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
#' @param v.util is a numeric vector containing utility values for each state
##Rivaroxaban Values
riv.u.ne <- c(0.75, 0.825, 1.00) #Utility of No Event
riv.u.vr <- c(0.57, 0.76, 0.95) #Utility of recurrence
riv.u.mb <- c(0.15, 0.61, 0.86) #Utility of major bleed
riv.u.nmb <- c(0.51, 0.65, 0.68) #Utility of CRNMB
riv.u.ot <- c(0.57, 0.68, 0.83) #Utility of off treatment
##Warfarin Values
war.u.ne <- c(0.75, 0.825, 1.00) #Utility of No Event
war.u.vr <- c(0.57, 0.76, 0.95) #Utility of recurrence
war.u.mb <- c(0.15, 0.61, 0.86) #Utility of major bleed
war.u.nmb <- c(0.51, 0.65, 0.68) #Utility of CRNMB
war.u.ot <- c(0.57, 0.68, 0.83) #Utility of off treatment
##Vectors for utility values
riv.df.util <- expand.grid(riv.u.ne, riv.u.vr, riv.u.mb, riv.u.nmb, riv.u.ot)
war.df.util <- expand.grid(war.u.ne, war.u.vr, war.u.mb, war.u.nmb, war.u.ot)
n_samples <- nrow(riv.df.util)
riv.v.util <- list()
war.v.util <- list()
for (i in 1:n_samples) {
riv.v.util[[i]] = unlist(riv.df.util[i, ], use.names=F)
war.v.util[[i]] = unlist(war.df.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 = rep(list(riv.v.cost), n_samples),
v.util = riv.v.util, ### SENS
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(rep(d.c, n_samples)),
d.e = as.numeric(rep(d.e, n_samples)),
v.n = rep(list(v.n), n_samples),
v.cost = rep(list(war.v.cost), n_samples),
v.util = war.v.util, ### SENS
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_04c_Utility.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")
```
Create df for tornado plot
```{r}
# Rivaroxaban
riv_df <- data.frame(
cbind(
riv.df.util,
sim_results_cost$riv,
sim_results_effectiveness$riv
)
)
riv_df <- riv_df %>%
mutate(CER = sim_results_cost.riv / sim_results_effectiveness.riv)
riv_df_base <- riv_df %>%
filter(Var1 == riv.u.ne[2],
Var2 == riv.u.vr[2],
Var3 == riv.u.mb[2],
Var4 == riv.u.nmb[2],
Var5 == riv.u.ot[2]) %>%
mutate(Parameter = "Base",
ParamTest = "Base")
riv_df_1 <- riv_df %>%
filter(Var1 != riv.u.ne[2],
Var2 == riv.u.vr[2],
Var3 == riv.u.mb[2],
Var4 == riv.u.nmb[2],
Var5 == riv.u.ot[2]) %>%
mutate(Parameter = c("riv.u.ne", "riv.u.ne"),
ParamTest = c("LowerBound", "UpperBound"))
riv_df_2 <- riv_df %>%
filter(Var1 == riv.u.ne[2],
Var2 != riv.u.vr[2],
Var3 == riv.u.mb[2],
Var4 == riv.u.nmb[2],
Var5 == riv.u.ot[2]) %>%
mutate(Parameter = c("riv.u.vr", "riv.u.vr"),
ParamTest = c("LowerBound", "UpperBound"))
riv_df_3 <- riv_df %>%
filter(Var1 == riv.u.ne[2],
Var2 == riv.u.vr[2],
Var3 != riv.u.mb[2],
Var4 == riv.u.nmb[2],
Var5 == riv.u.ot[2]) %>%
mutate(Parameter = c("riv.u.mb", "riv.u.mb"),
ParamTest = c("LowerBound", "UpperBound"))
riv_df_4 <- riv_df %>%
filter(Var1 == riv.u.ne[2],
Var2 == riv.u.vr[2],
Var3 == riv.u.mb[2],
Var4 != riv.u.nmb[2],
Var5 == riv.u.ot[2]) %>%
mutate(Parameter = c("riv.u.nmb", "riv.u.nmb"),
ParamTest = c("LowerBound", "UpperBound"))
riv_df_5 <- riv_df %>%
filter(Var1 == riv.u.ne[2],
Var2 == riv.u.vr[2],
Var3 == riv.u.mb[2],
Var4 == riv.u.nmb[2],
Var5 != riv.u.ot[2]) %>%
mutate(Parameter = c("riv.u.ot", "riv.u.ot"),
ParamTest = c("LowerBound", "UpperBound"))
riv_tornado <- rbind(riv_df_base,
riv_df_1,
riv_df_2,
riv_df_3,
riv_df_4,
riv_df_5)
# Save tornado data
saveRDS(riv_tornado, here::here("output","df_tornado_riv_utility.rds"))
```
```{r}
# Warfarin
war_df <- data.frame(
cbind(
war.df.util,
sim_results_cost$war,
sim_results_effectiveness$war
)
)
war_df <- war_df %>%
mutate(CER = sim_results_cost.war / sim_results_effectiveness.war)
war_df_base <- war_df %>%
filter(Var1 == war.u.ne[2],
Var2 == war.u.vr[2],
Var3 == war.u.mb[2],
Var4 == war.u.nmb[2],
Var5 == war.u.ot[2]) %>%
mutate(Parameter = "Base",
ParamTest = "Base")
war_df_1 <- war_df %>%
filter(Var1 != war.u.ne[2],
Var2 == war.u.vr[2],
Var3 == war.u.mb[2],
Var4 == war.u.nmb[2],
Var5 == war.u.ot[2]) %>%
mutate(Parameter = c("war.u.ne", "war.u.ne"),
ParamTest = c("LowerBound", "UpperBound"))
war_df_2 <- war_df %>%
filter(Var1 == war.u.ne[2],
Var2 != war.u.vr[2],
Var3 == war.u.mb[2],
Var4 == war.u.nmb[2],
Var5 == war.u.ot[2]) %>%
mutate(Parameter = c("war.u.vr", "war.u.vr"),
ParamTest = c("LowerBound", "UpperBound"))
war_df_3 <- war_df %>%
filter(Var1 == war.u.ne[2],
Var2 == war.u.vr[2],
Var3 != war.u.mb[2],
Var4 == war.u.nmb[2],
Var5 == war.u.ot[2]) %>%
mutate(Parameter = c("war.u.mb", "war.u.mb"),
ParamTest = c("LowerBound", "UpperBound"))
war_df_4 <- war_df %>%
filter(Var1 == war.u.ne[2],
Var2 == war.u.vr[2],
Var3 == war.u.mb[2],
Var4 != war.u.nmb[2],
Var5 == war.u.ot[2]) %>%
mutate(Parameter = c("war.u.nmb", "war.u.nmb"),
ParamTest = c("LowerBound", "UpperBound"))
war_df_5 <- war_df %>%
filter(Var1 == war.u.ne[2],
Var2 == war.u.vr[2],
Var3 == war.u.mb[2],
Var4 == war.u.nmb[2],
Var5 != war.u.ot[2]) %>%
mutate(Parameter = c("war.u.ot", "war.u.ot"),
ParamTest = c("LowerBound", "UpperBound"))
war_tornado <- rbind(war_df_base,
war_df_1,
war_df_2,
war_df_3,
war_df_4,
war_df_5)
# Save tornado data
saveRDS(war_tornado, here::here("output","df_tornado_war_utility.rds"))
```
Generate tornado plot
```{r }
# Generate Tornado Plot for Rivaroxaban
base.value = riv_tornado[1,"CER"]
df_tornado <- riv_tornado %>%
filter(! Parameter == "Base") %>%
select(Parameter, ParamTest, CER) %>%
pivot_wider(
names_from = ParamTest,
values_from = CER
) %>%
mutate(Diff = abs(LowerBound - UpperBound))
# Credit: kikoralston, StackOverflow
# get order of parameters according to size of intervals
order.parameters <- df_tornado %>% arrange(Diff) %>%
mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
select(Parameter) %>% unlist() %>% levels()
# width of columns in plot (value between 0 and 1)
width <- 0.95
# get data frame in shape for ggplot and geom_rect
df_tornado.2 <- df_tornado %>%
# gather columns Lower_Bound and Upper_Bound into a single column using gather
gather(key='type', value='output.value', LowerBound:UpperBound) %>%
# just reordering columns
select(Parameter, type, output.value, Diff) %>%
# create the columns for geom_rect
mutate(Parameter=factor(Parameter, levels=order.parameters),
ymin=pmin(output.value, base.value),
ymax=pmax(output.value, base.value),
xmin=as.numeric(Parameter)-width/2,
xmax=as.numeric(Parameter)+width/2)
# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
ggplot() +
geom_rect(data = df_tornado.2,
aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin, fill=type)) +
theme_bw() +
theme(legend.position = 'bottom',
legend.title = element_blank()) +
geom_hline(yintercept = base.value) +
scale_x_continuous(breaks = c(1:length(order.parameters)),
labels = order.parameters) +
coord_flip() +
ylab("Cost-Effectiveness Ratio") +
xlab("Parameter") +
ggtitle("Tornado Plot for Rivaroxaban")
```
```{r }
# Generate Tornado Plot for Warfarin
base.value = war_tornado[1,"CER"]
df_tornado <- war_tornado %>%
filter(! Parameter == "Base") %>%
select(Parameter, ParamTest, CER) %>%
pivot_wider(
names_from = ParamTest,
values_from = CER
) %>%
mutate(Diff = abs(LowerBound - UpperBound))
# Credit: kikoralston, StackOverflow
# get order of parameters according to size of intervals
order.parameters <- df_tornado %>% arrange(Diff) %>%
mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
select(Parameter) %>% unlist() %>% levels()
# width of columns in plot (value between 0 and 1)
width <- 0.95
# get data frame in shape for ggplot and geom_rect
df_tornado.2 <- df_tornado %>%
# gather columns Lower_Bound and Upper_Bound into a single column using gather
gather(key='type', value='output.value', LowerBound:UpperBound) %>%
# just reordering columns
select(Parameter, type, output.value, Diff) %>%
# create the columns for geom_rect
mutate(Parameter=factor(Parameter, levels=order.parameters),
ymin=pmin(output.value, base.value),
ymax=pmax(output.value, base.value),
xmin=as.numeric(Parameter)-width/2,
xmax=as.numeric(Parameter)+width/2)
# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
ggplot() +
geom_rect(data = df_tornado.2,
aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin, fill=type)) +
theme_bw() +
theme(legend.position = 'bottom',
legend.title = element_blank()) +
geom_hline(yintercept = base.value) +
scale_x_continuous(breaks = c(1:length(order.parameters)),
labels = order.parameters) +
coord_flip() +
ylab("Cost-Effectiveness Ratio") +
xlab("Parameter") +
ggtitle("Tornado Plot for Warfarin")
```