-
Notifications
You must be signed in to change notification settings - Fork 8
/
rpwexp_enroll.R
143 lines (135 loc) · 4.92 KB
/
rpwexp_enroll.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
# Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
# All rights reserved.
#
# This file is part of the simtrial program.
#
# simtrial is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#' Generate piecewise exponential enrollment
#'
#' With piecewise exponential enrollment rate generation any enrollment rate
#' distribution can be easily approximated.
#' `rpwexp_enroll()` is to support simulation of both the Lachin and Foulkes (1986)
#' sample size method for (fixed trial duration) as well as the
#' Kim and Tsiatis(1990) method (fixed enrollment rates and either
#' fixed enrollment duration or fixed minimum follow-up);
#' see [gsDesign::nSurv()].
#'
#' @param n Number of observations.
#' Default of `NULL` yields random enrollment size.
#' @param enroll_rate A data frame containing period duration (`duration`)
#' and enrollment rate (`rate`). for specified enrollment periods.
#' If necessary, last period will be extended to ensure enrollment
#' of specified `n`.
#'
#' @return A vector of random enrollment times.
#'
#' @importFrom data.table ":=" .N as.data.table last
#'
#' @export
#'
#' @examples
#' # Example 1
#' # Piecewise uniform (piecewise exponential inter-arrival times) for 10k patients enrollment
#' # Enrollment rates of 5 for time 0-100, 15 for 100-300, and 30 thereafter
#' x <- rpwexp_enroll(
#' n = 1e5,
#' enroll_rate = data.frame(
#' rate = c(5, 15, 30),
#' duration = c(100, 200, 100)
#' )
#' )
#' plot(x, 1:1e5,
#' main = "Piecewise uniform enrollment simulation",
#' xlab = "Time",
#' ylab = "Enrollment"
#' )
#'
#' # Example 2
#' # Exponential enrollment
#' x <- rpwexp_enroll(
#' n = 1e5,
#' enroll_rate = data.frame(rate = .03, duration = 1)
#' )
#' plot(x, 1:1e5,
#' main = "Simulated exponential inter-arrival times",
#' xlab = "Time",
#' ylab = "Enrollment"
#' )
rpwexp_enroll <- function(
n = NULL,
enroll_rate = data.frame(duration = c(1, 2), rate = c(2, 5))) {
# Take care of the simple case first if it is exponential enrollment
if (nrow(enroll_rate) == 1) {
# Stop with error message if only 1 enrollment period and the
# enrollment rate is less or equal with 0
if (enroll_rate$rate <= 0) {
stop("rpwexp_enroll: please specify > 0 enrollment rate, otherwise enrollment cannot finish.")
}
# Otherwise, return inter-arrival exponential times
else {
ans <- cumsum(stats::rexp(n = n, rate = enroll_rate$rate))
return(ans)
}
}
# Build `y` summarizes the start/end time, period order, etc.
y <- as.data.table(enroll_rate)
y[, period := seq_len(.N)]
y[, finish := cumsum(duration)]
y[, lambda := duration * rate]
y[, origin := c(0, finish[-length(finish)])]
y[, N := stats::rpois(n = .N, lambda = lambda)]
# Deal with extreme cases where none randomized in fixed intervals
if (sum(y$N) == 0) {
if (is.null(n)) {
ans <- NULL
return(ans)
}
if (last(enroll_rate$rate) <= 0) {
# Stop with error message if enrollment has not finished but enrollment rate for last period is less or equal with 0
stop("rpwexp_enroll: please specify > 0 enrollment rate for the last period; otherwise enrollment cannot finish.")
} else {
# Otherwise, return inter-arrival exponential times
ans <- cumsum(stats::rexp(n = n, rate = last(enroll_rate$rate))) + last(y$finish)
return(ans)
}
}
# Generate sorted uniform observations for Poisson count for each interval
z <- y[, .(enroll_time = sort(stats::runif(n = N, min = origin, max = finish))), by = "period"]
# If n not specified, return generated times
if (is.null(n)) {
ans <- z$enroll_time
return(ans)
}
# If n already achieved, return first n observations
if (nrow(z) >= n) {
ans <- z$enroll_time[1:n]
return(ans)
}
# After specified finite intervals, add required additional observations with
# exponential inter-arrival times
n_add <- n - nrow(z)
# Stop with error message if enrollment has not finished but
# enrollment rate for last period is less or equal with 0
if (last(enroll_rate$rate) <= 0) {
stop("rpwexp_enroll: please specify > 0 enrollment rate for the last period; otherwise enrollment cannot finish.")
}
# Otherwise, return inter-arrival exponential times
else {
ans <- c(
z$enroll_time,
cumsum(stats::rexp(n_add, rate = last(enroll_rate$rate))) + last(y$finish)
)
return(ans)
}
}