-
Notifications
You must be signed in to change notification settings - Fork 0
/
01_CEA_1yr_Riv.Rmd
264 lines (208 loc) · 10.1 KB
/
01_CEA_1yr_Riv.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
---
title: "CEA anticoag - rivaroxaban"
author: "Emily O'Neill, Andrew Huang"
date: "2023-06-08"
output: html_document
---
```{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(ggplot2)
library(gridExtra)
# Load Markov model function for both rivaroxaban and warfarin
source(here::here("functions", "function_markov.R"))
source(here::here("functions", "function_plot-trace.R"))
```
Base case simulation
Model hyperparameters
```{r}
##Model Input
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%
#The model states:
#(1) No Event = NE;
#(2) VTE Recurrence = VR;
#(3) Major Bleed = MB;
#(4) Clinically relevant Non-Major Bleed = NMB;
#(5) Off Treatment = OT
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
```
Model parameters for Rivaroxaban
```{r}
###Rivaroxaban###
##Health States
#Event rates per person-6months converted to person-year:
r.ne_vr <- 0.000 * 2 # rate of no event to VTE recurrence
r.ne_mb <- 0.019 * 2 # rate of no event to major bleed
r.ne_nmb <- 0.125 * 2 # rate of no event to CRNMB
r.nmb_ot <- 0.100 * 2 # rate of CNMB to come off treatment
r.mb_ot <- 0.100 * 2 # rate of MB to come off treatment
#Transition probabilities (per month cycle):
pr.ne_vr <- 1 - exp(-r.ne_vr * cycle_length) # from no event to VTE recurrence
pr.ne_mb <- 1 - exp(-r.ne_mb * cycle_length) # from no event to major bleed
pr.ne_nmb <- 1 - exp(-r.ne_nmb * cycle_length) # from no event to CRNMB
pr.ne_ne <- 1 - sum(pr.ne_vr, pr.ne_mb, pr.ne_nmb) # from no event to no event
pr.nmb_ot <- 1 - exp(-r.nmb_ot * cycle_length) # from CNMB to come off treatment
pr.mb_ot <- 1 - exp(-r.mb_ot * cycle_length) # from MB to come off treatment
##Costs (per month cycle (year/12))
c.ne <- 142 / 12 #Cost for staying healthy in rivaroxaban arm
c.vr <- 215 / 12 #Cost for VTE recurrence in Rivaroxaban arm ($215 annually)
c.mb <- 490 / 12 #Cost for major bleed in Rivaroxaban arm ($490 annually)
c.nmb <- 304 / 12 #Cost for CRNMB in Rivaroxaban arm ($304 annually)
c.ot <- 0 / 12 #Cost of off treatment ($0 annually)
##Utility Values
u.ne <- 0.825 #Utility of No Event
u.vr <- 0.76 #Utility of recurrence
u.mb <- 0.61 #Utility of major bleed
u.nmb <- 0.65 #Utility of CRNMB
u.ot <- 0.68 #Utility of off treatment
##Vectors for cost and utility values
v.cost <- c(c.ne, c.vr, c.mb, c.nmb, c.ot)
v.util <- c(u.ne, u.vr, u.mb, u.nmb, u.ot)
```
```{r}
##Run the Simulation
sim_markov_riv <- Markov(n.t, d.c, d.e, v.n, v.cost, v.util, v.M_1,
pr.ne_ne, pr.ne_vr, pr.ne_mb, pr.ne_nmb, pr.mb_ot, pr.nmb_ot)
sim_markov_riv
```
Plot cohort trace
```{r}
# Plot cohort trace
plot_trace(n.t, v.n, sim_markov_riv$m.TRr) + ggtitle("Rivaroxaban base case")
```
Save base case results
```{r}
saveRDS(sim_markov_riv, here::here("output","results_1yr_riv_noMC.rds"))
```
Now try using MC simulation for varying parameter values
```{r}
### CEA with Monte Carlo Simulation ###
#Number of MC Samples
n_samples <- 1000
#Set the random number seed for reproducibility
set.seed(123)
##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 <- as.data.frame(do.call(rbind, l.probr))
#df.probr <- df.probr %>% rowwise() %>% mutate(sum_test = sum(pr.ne_ne, pr.ne_vr, pr.ne_mb, pr.ne_nmb))
##Costs (per month cycle (year/12))
c.ne <- 142 / 12 #Cost for staying healthy in rivaroxaban arm
c.vr <- 215 / 12 #Cost for VTE recurrence in Rivaroxaban arm ($73)
c.mb <- 490 / 12 #Cost for major bleed in Rivaroxaban arm ($348)
c.nmb <- 304 / 12 #Cost for CRNMB in Rivaroxaban arm ($162)
c.ot <- 0 / 12 #Cost of off treatment
##Utility Values
u.ne <- 0.825 #Utility of No Event
u.vr <- 0.76 #Utility of recurrence
u.mb <- 0.61 #Utility of major bleed
u.nmb <- 0.65 #Utility of CRNMB
u.ot <- 0.59 #Utility of off treatment
##Vectors for cost and utility values
v.cost <- c(c.ne, c.vr, c.mb, c.nmb, c.ot)
v.util <- c(u.ne, u.vr, u.mb, u.nmb, u.ot)
## Generate an input dataframe containing parameters for each simulation
input_params <- as.data.frame(cbind(nsim = seq(1, n_samples),
n.t = rep(n.t, n_samples),
d.c = rep(d.c, n_samples),
d.e = rep(d.e, n_samples),
v.n = rep(list(v.n), n_samples),
v.cost = rep(list(v.cost), n_samples),
v.util = rep(list(v.util), n_samples),
v.M_1 = rep(list(v.M_1), n_samples),
pr.ne_ne = df.probr$pr.ne_ne,
pr.ne_vr = df.probr$pr.ne_vr,
pr.ne_mb = df.probr$pr.ne_mb,
pr.ne_nmb = df.probr$pr.ne_nmb,
pr.mb_ot = df.probr$pr.mb_ot,
pr.nmb_ot = df.probr$pr.nmb_ot))
input_params$nsim = as.character(input_params$nsim)
input_params$n.t = as.numeric(input_params$n.t)
input_params$d.c = as.numeric(input_params$d.c)
input_params$d.e = as.numeric(input_params$d.e)
input_params$pr.ne_ne = as.numeric(input_params$pr.ne_ne)
input_params$pr.ne_vr = as.numeric(input_params$pr.ne_vr)
input_params$pr.ne_mb = as.numeric(input_params$pr.ne_mb)
input_params$pr.ne_nmb = as.numeric(input_params$pr.ne_nmb)
input_params$pr.mb_ot = as.numeric(input_params$pr.mb_ot)
input_params$pr.nmb_ot = as.numeric(input_params$pr.nmb_ot)
```
```{r}
# Plot simulated values
grid.arrange(
ggplot(input_params, aes(x=pr.ne_ne)) + geom_histogram() + xlab("Probability of no event to no event"),
ggplot(input_params, aes(x=pr.ne_vr)) + geom_histogram() + xlab("Probability of no event to vr"),
ggplot(input_params, aes(x=pr.ne_mb)) + geom_histogram() + xlab("Probability of no event to mb"),
ggplot(input_params, aes(x=pr.ne_nmb)) + geom_histogram() + xlab("Probability of no event to nmb"),
ncol=1
)
grid.arrange(
ggplot(input_params, aes(x=pr.mb_ot)) + geom_histogram() + xlab("Probability of mb to off-treatment"),
ggplot(input_params, aes(x=pr.nmb_ot)) + geom_histogram() + xlab("Probability of nmb to off-treatment"),
ncol=1
)
```
```{r}
#Initialize dataframe to store results
results_riv <- data.frame(nsim = as.character(seq(1:n_samples)))
for (i in 1:n_samples){
sim_markov_riv <- Markov(n.t = input_params$n.t[[i]],
d.c = input_params$d.c[[i]],
d.e = input_params$d.e[[i]],
v.n = input_params$v.n[[i]],
v.cost = input_params$v.cost[[i]],
v.util = input_params$v.util[[i]],
v.M_1 = input_params$v.M_1[[i]],
pr.ne_ne = input_params$pr.ne_ne[[i]],
pr.ne_vr = input_params$pr.ne_vr[[i]],
pr.ne_mb = input_params$pr.ne_mb[[i]],
pr.ne_nmb = input_params$pr.ne_nmb[[i]],
pr.mb_ot = input_params$pr.mb_ot[[i]],
pr.nmb_ot = input_params$pr.nmb_ot[[i]])
results_riv$trace[i] <- list(sim_markov_riv$m.TRr)
results_riv$tc_hat[i] <- sim_markov_riv$tc_hat
results_riv$te_hat[i] <- sim_markov_riv$te_hat
}
```
```{r}
results <- cbind(input_params, results_riv[,-1])
results
```
```{r}
summary(results$tc_hat)
summary(results$te_hat)
```
Save sim results
```{r}
saveRDS(results, here::here("output","results_1yr_riv.rds"))
```