-
Notifications
You must be signed in to change notification settings - Fork 42
/
Appendix D_online_supp.R
253 lines (207 loc) · 12.4 KB
/
Appendix D_online_supp.R
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
############################################################################################
################# Microsimulation modeling using R: a tutorial #### 2018 ###################
############################################################################################
# This code forms the basis for the microsimulation model of the article:
#
# Krijkamp EM, Alarid-Escudero F, Enns EA, Jalal HJ, Hunink MGM, Pechlivanoglou P.
# Microsimulation modeling for health decision sciences using R: A tutorial.
# Med Decis Making. 2018;38(3):400–22.
#
# Please cite the article when using this code
#
# See GitHub for more information or code updates
# https://github.com/DARTH-git/Microsimulation-tutorial
#
#
# To program this tutorial we made use of
# R: 3.3.0 GUI 1.68 Mavericks build (7202)
# RStudio: Version 1.0.136 2009-2016 RStudio, Inc.
############################################################################################
################# Code of Appendix D #######################################################
############################################################################################
rm(list = ls()) # remove any variables in R's memory
##################################### Model input #########################################
# Model input
n.i <- 100000 # number of simulated individuals
n.t <- 30 # time horizon, 30 cycles
v.n <- c("H","S1","S2","D") # the model states: Healthy (H), Sick (S1), Sicker (S2), Dead (D)
n.s <- length(v.n) # the number of health states
v.M_1 <- rep("H", n.i) # everyone begins in the healthy state
d.c <- d.e <- 0.03 # equal discounting of costs and QALYs by 3%
v.Trt <- c("No Treatment", "Treatment") # store the strategy names
# Transition probabilities (per cycle)
p.HD <- 0.005 # probability to die when healthy
p.HS1 <- 0.15 # probability to become sick when healthy
p.S1H <- 0.5 # probability to become healthy when sick
p.S1S2 <- 0.105 # probability to become sicker when sick
rr.S1 <- 3 # rate ratio of death when sick vs healthy
rr.S2 <- 10 # rate ratio of death when sicker vs healthy
r.HD <- -log(1 - p.HD) # rate of death when healthy
r.S1D <- rr.S1 * r.HD # rate of death when sick
r.S2D <- rr.S2 * r.HD # rate of death when sicker
p.S1D <- 1 - exp(- r.S1D) # probability to die when sick
p.S2D <- 1 - exp(- r.S2D) # probability to die when sicker
# Cost and utility inputs
c.H <- 2000 # cost of remaining one cycle healthy
c.S1 <- 4000 # cost of remaining one cycle sick
c.S2 <- 15000 # cost of remaining one cycle sicker
c.Trt <- 12000 # cost of treatment (per cycle)
u.H <- 1 # utility when healthy
u.S1 <- 0.75 # utility when sick
u.S2 <- 0.5 # utility when sicker
u.Trt <- 0.95 # utility when sick(er) and being treated
##################################### Functions ###########################################
# THE NEW samplev() FUNCTION
# efficient implementation of the rMultinom() function of the Hmisc package ####
samplev <- function (probs, m) {
d <- dim(probs)
n <- d[1]
k <- d[2]
lev <- dimnames(probs)[[2]]
if (!length(lev))
lev <- 1:k
ran <- matrix(lev[1], ncol = m, nrow = n)
U <- t(probs)
for(i in 2:k) {
U[i, ] <- U[i, ] + U[i - 1, ]
}
if (any((U[k, ] - 1) > 1e-05))
stop("error in multinom: probabilities do not sum to 1")
for (j in 1:m) {
un <- rep(runif(n), rep(k, n))
ran[, j] <- lev[1 + colSums(un > U)]
}
ran
}
# The MicroSim function for the simple microsimulation of the 'Sick-Sicker' model keeps track of what happens to each individual during each cycle.
MicroSim <- function(v.M_1, n.i, n.t, v.n, d.c, d.e, TR.out = TRUE, TS.out = TRUE,
Trt = FALSE, seed = 1) {
# Arguments:
# v.M_1: vector of initial states for individuals
# n.i: number of individuals
# n.t: total number of cycles to run the model
# v.n: vector of health state names
# d.c: discount rate for costs
# d.e: discount rate for health outcome (QALYs)
# TR.out: should the output include a Microsimulation trace? (default is TRUE)
# TS.out: should the output include a matrix of transitions between states? (default is TRUE)
# Trt: are the n.i individuals receiving treatment? (scalar with a Boolean value, default is FALSE)
# seed: starting seed number for random number generator (default is 1)
# Makes use of:
# Probs: function for the estimation of transition probabilities
# Costs: function for the estimation of cost state values
# Effs: function for the estimation of state specific health outcomes (QALYs)
v.dwc <- 1 / (1 + d.c) ^ (0:n.t) # calculate the cost discount weight based on the discount rate d.c
v.dwe <- 1 / (1 + d.e) ^ (0:n.t) # calculate the QALY discount weight based on the discount rate d.e
# Create the matrix capturing the state name/costs/health outcomes for all individuals at each time point
m.M <- m.C <- m.E <- matrix(nrow = n.i, ncol = n.t + 1,
dimnames = list(paste("ind", 1:n.i, sep = " "),
paste("cycle", 0:n.t, sep = " ")))
m.M[, 1] <- v.M_1 # indicate the initial health state
set.seed(seed) # set the seed for every individual for the random number generator
m.C[, 1] <- Costs(m.M[, 1], Trt) # estimate costs per individual for the initial health state
m.E[, 1] <- Effs (m.M[, 1], Trt) # estimate QALYs per individual for the initial health state
for (t in 1:n.t) {
m.p <- Probs(m.M[, t]) # calculate the transition probabilities at cycle t
m.M[, t + 1] <- samplev( prob = m.p, m = 1) # sample the next health state and store that state in matrix m.M
m.C[, t + 1] <- Costs(m.M[, t + 1], Trt) # estimate costs per individual during cycle t + 1 conditional on treatment
m.E[, t + 1] <- Effs( m.M[, t + 1], Trt) # estimate QALYs per individual during cycle t + 1 conditional on treatment
cat('\r', paste(round(t/n.t * 100), "% done", sep = " ")) # display the progress of the simulation
} # close the loop for the time points
tc <- m.C %*% v.dwc # total (discounted) cost per individual
te <- m.E %*% v.dwe # total (discounted) QALYs per individual
tc_hat <- mean(tc) # average (discounted) cost
te_hat <- mean(te) # average (discounted) QALYs
if (TS.out == TRUE) { # create a matrix of transitions across states
TS <- paste(m.M, cbind(m.M[, -1], NA), sep = "->") # transitions from one state to the other
TS <- matrix(TS, nrow = n.i)
rownames(TS) <- paste("Ind", 1:n.i, sep = " ") # name the rows
colnames(TS) <- paste("Cycle", 0:n.t, sep = " ") # name the columns
} else {
TS <- NULL
}
if (TR.out == TRUE) {
TR <- t(apply(m.M, 2, function(x) table(factor(x, levels = v.n, ordered = TRUE))))
TR <- TR / n.i # create a distribution trace
rownames(TR) <- paste("Cycle", 0:n.t, sep = " ") # name the rows
colnames(TR) <- v.n # name the columns
} else {
TR <- NULL
}
results <- list(m.M = m.M, m.C = m.C, m.E = m.E, tc = tc, te = te, tc_hat = tc_hat, te_hat = te_hat, TS = TS, TR = TR) # store the results from the simulation in a list
return(results) # return the results
} # end of the MicroSim function
#### Probability function
# The Probs function that updates the transition probabilities of every cycle is shown below.
Probs <- function(M_it) {
# M_it: health state occupied by individual i at cycle t (character variable)
m.p.it <- matrix(NA, n.s, n.i) # create vector of state transition probabilities
rownames(m.p.it) <- v.n # assign names to the vector
# update the v.p with the appropriate probabilities
m.p.it[,M_it == "H"] <- c(1 - p.HS1 - p.HD, p.HS1, 0, p.HD) # transition probabilities when healthy
m.p.it[,M_it == "S1"] <- c(p.S1H, 1- p.S1H - p.S1S2 - p.S1D, p.S1S2, p.S1D) # transition probabilities when sick
m.p.it[,M_it == "S2"] <- c(0, 0, 1 - p.S2D, p.S2D) # transition probabilities when sicker
m.p.it[,M_it == "D"] <- c(0, 0, 0, 1) # transition probabilities when dead
ifelse(colSums(m.p.it) == 1, return(t(m.p.it)), print("Probabilities do not sum to 1")) # return the transition probabilities or produce an error
}
### Costs function
# The Costs function estimates the costs at every cycle.
Costs <- function (M_it, Trt = FALSE) {
# M_it: health state occupied by individual i at cycle t (character variable)
# Trt: is the individual being treated? (default is FALSE)
c.it <- 0 # by default the cost for everyone is zero
c.it[M_it == "H"] <- c.H # update the cost if healthy
c.it[M_it == "S1"] <- c.S1 + c.Trt * Trt # update the cost if sick conditional on treatment
c.it[M_it == "S2"] <- c.S2 + c.Trt * Trt # update the cost if sicker conditional on treatment
c.it[M_it == "D"] <- 0 # update the cost if dead
return(c.it) # return the costs
}
### Health outcome function
# The Effs function to update the utilities at every cycle.
Effs <- function (M_it, Trt = FALSE, cl = 1) {
# M_it: health state occupied by individual i at cycle t (character variable)
# Trt: is the individual treated? (default is FALSE)
# cl: cycle length (default is 1)
u.it <- 0 # by default the utility for everyone is zero
u.it[M_it == "H"] <- u.H # update the utility if healthy
u.it[M_it == "S1"] <- Trt * u.Trt + (1 - Trt) * u.S1 # update the utility if sick conditional on treatment
u.it[M_it == "S2"] <- u.S2 # update the utility if sicker
u.it[M_it == "D"] <- 0 # update the utility if dead
QALYs <- u.it * cl # calculate the QALYs during cycle t
return(QALYs) # return the QALYs
}
##################################### Run the simulation ##################################
# START SIMULATION
p = Sys.time()
sim_no_trt <- MicroSim(v.M_1, n.i, n.t, v.n, d.c, d.e, Trt = FALSE) # run for no treatment
sim_trt <- MicroSim(v.M_1, n.i, n.t, v.n, d.c, d.e, Trt = TRUE) # run for treatment
comp.time = Sys.time() - p
comp.time
################################# Cost-effectiveness analysis #############################
# store the mean costs (and the MCSE) of each strategy in a new variable C (vector costs)
v.C <- c(sim_no_trt$tc_hat, sim_trt$tc_hat)
sd.C <- c(sd(sim_no_trt$tc), sd(sim_trt$tc)) / sqrt(n.i)
# store the mean QALYs (and the MCSE) of each strategy in a new variable E (vector effects)
v.E <- c(sim_no_trt$te_hat, sim_trt$te_hat)
sd.E <- c(sd(sim_no_trt$te), sd(sim_trt$te)) / sqrt(n.i)
delta.C <- v.C[2] - v.C[1] # calculate incremental costs
delta.E <- v.E[2] - v.E[1] # calculate incremental QALYs
sd.delta.E <- sd(sim_trt$te - sim_no_trt$te) / sqrt(n.i) # Monte Carlo Squared Error (MCSE) of incremental costs
sd.delta.C <- sd(sim_trt$tc - sim_no_trt$tc) / sqrt(n.i) # Monte Carlo Squared Error (MCSE) of incremental QALYs
ICER <- delta.C / delta.E # calculate the ICER
results <- c(delta.C, delta.E, ICER) # store the values in a new variable
# Create full incremental cost-effectiveness analysis table
table_micro <- data.frame(
c(round(v.C, 0), ""), # costs per arm
c(round(sd.C, 0), ""), # MCSE for costs
c(round(v.E, 3), ""), # health outcomes per arm
c(round(sd.E, 3), ""), # MCSE for health outcomes
c("", round(delta.C, 0), ""), # incremental costs
c("", round(sd.delta.C, 0),""), # MCSE for incremental costs
c("", round(delta.E, 3), ""), # incremental QALYs
c("", round(sd.delta.E, 3),""), # MCSE for health outcomes (QALYs) gained
c("", round(ICER, 0), "") # ICER
)
rownames(table_micro) <- c(v.Trt, "* are MCSE values") # name the rows
colnames(table_micro) <- c("Costs", "*", "QALYs", "*", "Incremental Costs", "*", "QALYs Gained", "*", "ICER") # name the columns
table_micro # print the table