-
Notifications
You must be signed in to change notification settings - Fork 28
/
simulate_flexsurvreg.R
93 lines (93 loc) · 3.78 KB
/
simulate_flexsurvreg.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
##' Simulate censored time-to-event data from a fitted flexsurvreg model
##'
##' @param object Object returned by \code{\link{flexsurvreg}}.
##'
##' @param nsim Number of simulations per row in \code{newdata}.
##'
##' @param seed Random number seed. This is returned with the result of this
##' function, as described in \code{\link{simulate}} for the \code{lm} method.
##'
##' @param newdata Data frame defining alternative sets of covariate values to simulate with.
##' If omitted, this defaults to the data originally used to fit the model.
##'
##' @template start
##'
##' @param censtime A right-censoring time, or vector of times matching the rows
##' of \code{newdata}. If \code{NULL} (the default) then uncensored times to events
##' are simulated.
##'
##' @param tidy If \code{TRUE} then a "tidy" or "long"-format data frame is
##' returned, with rows defined by combinations of covariates and simulation
##' replicates. The simulation replicate is indicated in the column named \code{i}.
##'
##' If \code{FALSE}, then a data frame is returned with one row per set of
##' covariate values, and different columns for different simulation
##' replicates. This is the traditional format for `simulate` methods in base
##' R.
##'
##' In either case, the simulated time and indicator for whether the time is
##' an event time (rather than a time of right-censoring) are returned in
##' different columns.
##'
##' @param ... Other arguments (not currently used).
##'
##' @return A data frame, with format determined by whether \code{tidy} was specified.
##'
##' @examples
##' fit <- flexsurvreg(formula = Surv(futime, fustat) ~ rx, data = ovarian, dist="weibull")
##' fit2 <- flexsurvspline(formula = Surv(futime, fustat) ~ rx, data = ovarian, k=3)
##' nd = data.frame(rx=1:2)
##' simulate(fit, seed=1002, newdata=nd)
##' simulate(fit, seed=1002, newdata=nd, start=500)
##' simulate(fit2, nsim=3, seed=1002, newdata=nd)
##' simulate(fit2, nsim=3, seed=1002, newdata=nd, start=c(500,1000))
##'
##' @export
simulate.flexsurvreg <- function(object, nsim=1, seed=NULL,
newdata=NULL, start=NULL, censtime=NULL, tidy=FALSE,...) {
if (is.null(newdata))
newdata <- model.frame(object)
else if (!is.data.frame(newdata))
stop("`newdata` should be a data frame")
if (!is.null(seed)) {
set.seed(seed)
attr(seed, "kind") <- as.list(RNGkind())
seed_attr <- seed
} else seed_attr <- .Random.seed
nd <- nrow(newdata)
if (!is.numeric(nsim) || (nsim < 1)) stop("`nsim` should be a number >= 1")
n <- nd*nsim
if (!is.null(start)) {
if(!(length(start) %in% c(1, nd)))
stop(sprintf("`start` of length %s, should be of length 1 or %s = nrow(newdata)"),
length(start), nd)
start <- rep(start, length.out=n)
} else start <- 0
U <- runif(n, 0, 1)
newdata <- newdata[rep(1:nd, each=nsim), , drop=FALSE]
time <- summary(object, newdata=newdata, start=start,
type="quantile", quantiles=U,
cross=FALSE, ci=FALSE, se=FALSE, tidy=TRUE)$est
if (is.null(censtime))
censtime <- rep(Inf, nd)
if (length(censtime) ==1)
censtime <- rep(censtime, nd)
if (length(censtime)!=nd)
stop(sprintf("`censtime` of length %s, should be of length %s = nrow(newdata)",
length(censtime), nd))
censtime <- rep(censtime, each = nsim)
event <- as.numeric(time <= censtime)
time <- ifelse(event, time, censtime)
if (tidy) {
res <- cbind(newdata, i=rep(1:nsim, nd), time=time, event=event)
rownames(res) <- NULL
} else {
time <- matrix(time, nrow=nd, ncol=nsim, byrow=TRUE)
event <- matrix(event, nrow=nd, ncol=nsim, byrow=TRUE)
colnames(time) <- paste0("time_",1:nsim)
colnames(event) <- paste0("event_",1:nsim)
res <- as.data.frame(cbind(time, event))
}
attr(res, "seed") <- seed_attr
res
}