/
utilmod.R
112 lines (107 loc) · 5.24 KB
/
utilmod.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
#' Create utility model
#'
#' Create a model for health state utility given values of utility by
#' health state, treatment, and time sampled from a probability distribution.
#' @param n The number of random observations of the parameters to draw.
#' @param struct A \code{\link{model_structure}} object.
#' @param patients A data table returned from \code{\link{create_patients}}.
#' @param ae_probs An "ae_probs" object as returned by \code{\link{ae_probs}}.
#' @param params_utility Parameter estimates for health state utilities and
#' adverse event disutilities in the same format as \code{\link{params_utility}}.
#' @param ae_duration Duration of time over with disutility from adverse events
#' should accrue. If \code{"month"}, then disutility accrues over the first
#' month of treatment; if \code{progression} then disutility accrues until
#' disease progression (i.e., over the entire duration of time spent in
#' stable disease).
#' @return An object of class "StateVals" from the
#' \href{https://hesim-dev.github.io/hesim/}{hesim} package.
#' @examples
#' # Treatment sequences
#' txseq1 <- txseq(first = "erlotinib",
#' second = c("osimertinib", "PBDC"),
#' second_plus = c("PBDC + bevacizumab", "PBDC + bevacizumab"))
#' txseq2 <- txseq(first = "gefitinib",
#' second = c("osimertinib", "PBDC"),
#' second_plus = c("PBDC + bevacizumab", "PBDC + bevacizumab"))
#' txseqs <- txseq_list(seq1 = txseq1, seq2 = txseq2)
#'
#' # Patient population
#' pats <- create_patients(n = 2)
#'
#' ## Model structure
#' struct <- model_structure(txseqs, dist = "weibull")
#'
#' ## Utility model
#' n_samples <- 2
#' ae_probs <- ae_probs(n = n_samples, struct = struct)
#' utilmod <- create_utilmod(n = n_samples, struct = struct, patients = pats,
#' ae_probs = ae_probs)
#' @export
create_utilmod <- function(n = 100, struct, patients,
ae_probs,
params_utility = iviNSCLC::params_utility,
ae_duration = c("month", "progression")){
ae_duration <- match.arg(ae_duration)
# hesim data
strategies <- data.table(strategy_id = 1:length(struct$txseqs))
hesim_dat <- hesim::hesim_data(strategies = strategies,
patients = patients)
# Utility by health state
states <- create_states(struct)[get("state_name") != "D"]
state_utility_tbl <- merge(states, params_utility$state_utility, by = "state_name")
setorderv(state_utility_tbl, "state_id")
state_utility_dist <- matrix(stats::rnorm(n * nrow(state_utility_tbl),
state_utility_tbl$mean,
state_utility_tbl$se),
nrow = n, byrow = TRUE)
# Disutilities from adverse events
disutility_ae_dist <- matrix(stats::rnorm(n * nrow(params_utility$ae_disutility),
params_utility$ae_disutility$mean,
params_utility$ae_disutility$se),
nrow = n, byrow = TRUE)
colnames(disutility_ae_dist) <- params_utility$ae_disutility$ae_abb
## Weight by adverse events
indices <- match(params_utility$ae_disutility$ae_abb, names(ae_probs))
if (any(is.na(indices))){
stop(paste0("The adverse event abbreviations in 'params_utility$ae_disutility' do not match ",
"the names of the adverse events in 'ae_probs'."))
}
expected_disutility <- vector(mode = "list", length = length(ae_probs))
names(expected_disutility) <- names(ae_probs)
for (i in 1:length(ae_probs)){
name_i <- names(ae_probs)[i]
expected_disutility[[i]] <- ae_probs[[i]] * disutility_ae_dist[, name_i]
}
expected_disutility <- Reduce('+', expected_disutility)
# Create utility model
tbl1 <- tbl2 <- vector(mode = "list", length = nrow(strategies))
for (i in 1:nrow(strategies)){
if (ae_duration == "month"){
tbl1[[i]] <- data.table(strategy_id = i,
state_id = rep(states$state_id, each = n),
sample = rep(1:n, times = nrow(states)),
value = c(state_utility_dist) - expected_disutility[, i],
time_start = 0)
tbl2[[i]] <- data.table(strategy_id = i,
state_id = rep(states$state_id, each = n),
sample = rep(1:n, times = nrow(states)),
value = c(state_utility_dist),
time_start = 1/12)
} else {
tbl1[[i]] <- data.table(strategy_id = i,
state_id = rep(1, each = n), # State S1
sample = 1:n,
value = state_utility_dist[, 1] - expected_disutility[, i])
tbl2[[i]] <- data.table(strategy_id = i,
state_id = rep(states$state_id[-1], each = n),
sample = rep(1:n, times = nrow(states[-1])),
value = c(state_utility_dist[, -1]))
}
}
tbl1 <- rbindlist(tbl1)
tbl2 <- rbindlist(tbl2)
tbl <- rbind(tbl1, tbl2)
tbl <- hesim::stateval_tbl(tbl, dist = "custom", hesim_data = hesim_dat)
utilmod <- hesim::create_StateVals(tbl, n = n)
return(utilmod)
}