1

Please consider the following:

My aim is to draw random survival times from a flexible parametric multi-state model fitted with flexsurvreg (more specifically msfit.flexsurvreg) and applying some hazard ratio (HR, in this example set to 0.2).

I found an example to generate random survival times using any hazard function here. This is also were I apply the HR.

Problem
With the actual data, I receive an error once the HR is below the value of 0.2: Error in uniroot((function(x) { : no sign change found in 1000 iterations.

This does not happen in the reproducible example below.

Questions
Is there another, better way than the one below to draw random survival times while applying a HR?

Can someone indicate why the "no sign change" error may occur and how this can be fixed?

Any help is greatly appreciated!

Minimal reproducible example

# Load package
library(flexsurv)
#> Loading required package: survival

# Load data
data("bosms4")

# Define hazard ratio
hr <- 0.2

# Fit model (weibull)
crwei <- flexsurvreg(formula = Surv(years, status) ~ trans + shape(trans), 
                     data = bosms3, dist = "weibull")

# Create transition matrix
Q <- rbind(c(NA,1,2),c(NA,NA,3), c(NA,NA,NA))

# Capture parameters
pars <- pars.fmsm(crwei, trans=Q, newdata=data.frame(trans=1:3))

# Code from https://eurekastatistics.com/generating-random-survival-times-from-any-hazard-function/ ----

inverse = function(fn, min_x, max_x){
        # Returns the inverse of a function for a given range.
        # E.g. inverse(sin, 0, pi/2)(sin(pi/4)) equals pi/4 because 0 <= pi/4 <= pi/2
        fn_inv = function(y){
                uniroot((function(x){fn(x) - y}), lower=min_x, upper=max_x)[1]$root
        }
        return(Vectorize(fn_inv))
}

integrate_from_0 = function(fn, t){
        int_fn = function(t) integrate(fn, 0, t)
        result = sapply(t, int_fn)
        value  = unlist(result["value",])
        msg    = unlist(result["message",])
        value[which(msg != "OK")] = NA
        return(value)
}

random_survival_times = function(hazard_fn, n, max_time=10000){
        # Given a hazard function, returns n random time-to-event observations.
        cumulative_density_fn = function(t) 1 - exp(-integrate_from_0(hazard_fn, t))
        inverse_cumulative_density_fn = inverse(cumulative_density_fn, 0, max_time)
        return(inverse_cumulative_density_fn(runif(n)))
}

# Run with data ----
random_survival_times(hazard_fn = function(t){crwei$dfns$h(t, pars[[1]][1],                                        pars[[1]][2]) * hr}, n = 100)
#> Error in integrate(fn, 0, t): non-finite function value

# Adapt random_survival time function replacing 0 with 0.1 ----

random_survival_times <- function(hazard_fn, n, max_time=10000){
        # Given a hazard function, returns n random time-to-event observations.
        cumulative_density_fn = function(t) 1 - exp(-integrate_from_0(hazard_fn, t))
        inverse_cumulative_density_fn = inverse(cumulative_density_fn, 0.1, max_time) 
        return(inverse_cumulative_density_fn(runif(n)))
}

# Run again with data ----
random_survival_times(hazard_fn = function(t){crwei$dfns$h(t, pars[[1]][1], pars[[1]][2]) * hr}, n = 100)
#> Error in uniroot((function(x) {: f() values at end points not of opposite sign

# Adapt inverse adding extendedInt = "yes" ----

inverse <- function(fn, min_x, max_x){
        # Returns the inverse of a function for a given range.
        # E.g. inverse(sin, 0, pi/2)(sin(pi/4)) equals pi/4 because 0 <= pi/4 <= pi/2
        fn_inv <- function(y){
                uniroot((function(x){fn(x) - y}), lower=min_x, upper=max_x,
                        extendInt = "yes" # extendInt added because of error on some distributions: "Error in uniroot((function(x) { : f() values at end points not of opposite sign. Solution found here: https://stackoverflow.com/questions/38961221/uniroot-solution-in-r
                )[1]$root
        }
        return(Vectorize(fn_inv))
}

# Run again with data ----
res <- random_survival_times(hazard_fn = function(t){crwei$dfns$h(t, pars[[1]][1],                                        pars[[1]][2]) * hr}, n = 100)

res[1:5]
#> [1]   1.074281  13.688663  30.896637 159.643827  15.805103

Created on 2022-10-18 with reprex v2.0.2

Frederick
  • 810
  • 8
  • 28
  • The cited authority uses `random_survival_times` without a prior library call to load a package with a function of that name. Also unclear is why you are using NAs in a transition matrix. – IRTFM Oct 19 '22 at 08:04
  • The function `random_survival_times()` is defined at the end of the cited blogpost. I agree that this is confusing in the blogpost, but in the reproducible example above this is well defined. Using `NA` in the transitions matrix is used to "fill" the matrix where no transitions can take place. See `?flexsurv::pars.fmsm()` for the `trans` argument and `?mstate::msprep()`. – Frederick Oct 19 '22 at 09:27
  • I’m puzzled by the use of NA instead of 0 for an impossible transition. One generally checks the rowSums all equaling unity to ensure validity. – IRTFM Oct 19 '22 at 15:02

1 Answers1

1

This method of sampling survival times basically works by sampling a random uniform(0,1) number p and finding x for which the survival probability is p. The uniroot step is used to solve S(x) = p by a numerical search. In your case, it is having difficulty finding a solution after 1000 steps.

I've seen this happen, and fixed by adding, e.g. uniroot(..., maxiter=10000) to tell it to try a bit harder to find the solution. That's always been enough in my tests, though those may be limited. If that doesn't work, I'd advise digging in by hand and examining the survival curve that you are trying to invert - it may be invalid due to some parameter value being extreme.

(This kind of thing is done in the function qgeneric in the flexsurv package, though it borrows a vectorised version of uniroot from the rstpm2 package which is faster.)

Chris Jackson
  • 871
  • 5
  • 7
  • Thank you for your reply. We tried scaling-up the maximum iterations like you suggested without any success unfortunately. – Frederick Oct 26 '22 at 14:37