-
Notifications
You must be signed in to change notification settings - Fork 43
/
Appendix B_online_supp-samplev.R
273 lines (223 loc) · 13.8 KB
/
Appendix B_online_supp-samplev.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
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
############################################################################################
################# 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
rp.S1S2 <- 0.2 # increase of the mortality rate with every additional year being sick
# 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
ru.S1S2 <- 0.03 # decrease in utility of treated sick individuals with every additional year being sick/sicker
v.x <- runif(n.i, 0.95, 1.05) # vector capturing individuals' effect modifier at baseline
##################################### 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, X = NULL, 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
# X: vector or matrix of individual characteristics
# 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
# create the dur variable that stores the number of consecutive cycles the individual occupies either when sick or sicker
dur <- rep(0, n.i) # the individual start without history
m.C[, 1] <- Costs(m.M[, 1], Trt) # estimate costs per individual for the initial health state
m.E[, 1] <- Effs (m.M[, 1], dur, Trt, X = X) # estimate QALYs per individual for the initial health state
for (t in 1:n.t) { # t <- 3
m.p <- Probs(m.M[, t], dur) # 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], dur, Trt, X = X) # estimate QALYs per individual during cycle t + 1 conditional on treatment
dur <- ifelse(m.M[, t + 1] == "S1" | m.M[, t + 1] == "S2",
dur + 1,
0)
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, dur) {
# M_it: health state occupied by individual i at cycle t (character variable)
# dur: the duration of being sick (sick/sicker)
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 probabilities of death after first converting them to rates and applying the rate ratio
r.S1D <- - log(1 - p.S1D)
r.S2D <- - log(1 - p.S2D)
p.S1D <- 1 - exp(- r.S1D * (1 + dur[M_it == "S1"] * rp.S1S2)) # calculate p.S1D conditional on duration of being sick/sicker
p.S2D <- 1 - exp(- r.S2D * (1 + dur[M_it == "S2"] * rp.S1S2)) # calculate p.S2D conditional on duration of being sick/sicker
# 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"] <- rbind(p.S1H, 1- p.S1H - p.S1S2 - p.S1D, p.S1S2, p.S1D) # transition probabilities when sick
m.p.it[,M_it == "S2"] <- rbind(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, dur, Trt = FALSE, cl = 1, X = NULL) {
# M_it: health state occupied by individual i at cycle t (character variable)
# dur: the duration of being sick/sicker
# Trt: is the individual treated? (default is FALSE)
# cl: cycle length (default is 1)
# X: the vector or matrix of individual characteristics (optional)
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"] <- X * Trt * (u.Trt - dur * ru.S1S2) + (1 - Trt) * u.S1 # update the utility if sick conditional on treatment and duration of being sick/sicker
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, X = v.x, d.c, d.e, Trt = FALSE, seed = 100) # run for no treatment
sim_trt <- MicroSim(v.M_1, n.i, n.t, v.n, X = v.x, d.c, d.e, Trt = TRUE, seed = 100) # 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