simulate function

[mclements]

Terry,

I would like to implement a simulate method for the survreg class. I have a draft implementation below. Would you consider including this in the survival package? (If not, then I could include this method in the microsimulation package).

Note that I have used an argument t0 for delayed entry. Is this a reasonable name?

Sincerely, Mark Clements.

#' Simulate event times from a survreg object
#' @param object survreg object
#' @param nsim number of simulations per row in newdata
#' @param seed random number seed
#' @param newdata data-frame for defining the covariates for the simulations. Required.
#' @param t0 delayed entry time. Defaults to NULL (which assumes that t0=0)
#' @param ... other arguments (not currently used)
#' @return vector of event times with nsim repeats per row in newdata
#' @importFrom stats simulate predict runif
#' @importFrom survival survreg.distributions
#' @rdname simulate
#' @export
#' @examples
#' library(survival)
#' fit <- survreg(Surv(time, status) ~ ph.ecog + age + sex + strata(sex),
#'                data = lung)
#' nd = transform(expand.grid(ph.ecog=0:1, sex=1:2), age=60)
#' simulate(fit, seed=1002, newdata=nd)
#' simulate(fit, seed=1002, newdata=nd, t0=500)
simulate.survreg = function(object, nsim=1, seed=NULL, newdata, t0=NULL, ...) {
    stopifnot(inherits(object, "survreg"),
              is.list(newdata))
    if (!is.null(seed)) set.seed(seed)
    lp = predict(object, newdata=newdata, type="lp")
    lp = lp[rep(1:length(lp), each=nsim)]
    n = length(lp)
    if (!is.null(strata <- attr(object$terms, "specials")$strata)) {
        scale = object$scale[eval(attr(object$terms,"variables")[[strata+1]], newdata)]
        scale = rep(scale, each=nsim)
    }
    else scale = rep(object$scale,n)
    if (is.character(object$dist)) 
        dd <- survival::survreg.distributions[[object$dist]]
    else dd <- object$dist
    if (is.null(dd$itrans)) {
        trans <- function(x) x
        itrans <- function(x) x
    }
    else {
        trans <- dd$trans
        itrans <- dd$itrans
    }
    if (!is.null(dd$dist)) 
        dd <- survival::survreg.distributions[[dd$dist]]
    if (!is.null(t0)) {
        stopifnot(length(t0) %in% c(1, length(newdata[[1]])))
        if (length(t0)==1) t0=rep(t0,n)
        else t0 = rep(t0,each=nsim)
        F0 = dd$density((trans(t0)-lp)/scale,object$parm)[,1]
    } else F0 = rep(0,n)
    itrans(lp+scale*dd$quantile(1-(runif(n, F0, 1)-F0), object$parm))
}

[therneau]

1. I think this is a sensible idea.
2. I would not include a t0 argument, since the survreg models do not allow for delayed entry. Someone can modify the returned values, if they desire.
3. You need to add code and a help file, separately. I am very much not a fan of roxygen, and you won't find it in my code base.
4. I think it makes the most sense to invoke rsurvreg for computation, since it exists and has been tested. (There is also dsurvreg, psurvreg, qsurvreg, in parallel to pnorm, qnorm, rnorm, etc. )
5. Use of [ on a terms object is not reliable; there are bugs in the [.terms function. I have tried to get them repaired, but so far without success. For example, add and offset() in the formula before the strata(). Your useage might be okay though in that respect. But, you need to allow for the person who has two strata() terms in their call. (Rare, but people do it.)