0
$\begingroup$

I am trying to generate survival times from a weibull distribution using the gems package, and calculate the shape and scale parameters of that weibull distribution using flexsurvreg

The gems package gives you the option to use a built in weibull function based on rweibull, or define the hazard manually. When I use the built in rweibull function option to generate data, flexsurvreg successfully calculates the shape and scale parameters of the weibull distribution. When I use the manual hazard function option, flexsurvreg cannot calculate the shape and scale of the generating distribution.

I have specified the hazard function of the weibull distribution to be (shape/scale)*(t/scale)^(shape - 1). Is this the wrong hazard function associated with the weibull distribution? Or is the numerical approximation procedure used to calculate the survival distribution (which only gets used when the hazard function is defined manually), not accurate? Or something else?

This question is very similar to another I posted, but this question is specifically about the definition of the hazard function for a weibull survival distribution, so I wanted to separate it out: Unable to calculate the shape and scale parameters of weibull distribution from data simulated through either gems or simsurv in R

Many thanks for any insight on this. Full reproducible example below.

#install.packages("survival")
#install.packages("gems")
#install.packages("flexsurv")
#install.packages("simsurv")

library(survival)
library(gems)
library(flexsurv)
library(simsurv)

### Set seed
set.seed(101)

### Set number of patients
npat <- 10000

#######################################################################
# Generate survival times using gems, with a manual hazard function #
#######################################################################

### Going to set up a simple model that goes from state 1, to state 2

## Generate an empty hazard matrix
hf <- generateHazardMatrix(2)

## The states named "impossible" are the ones which I need to change the name of
hf

## Define the transitions as weibull for now
## I am leaving the transition from 1 -> 3 as "impossible' to that the strcture is as planned


## Define the transitions as weibull using the manual hazard definition
hf[[1, 2]] <- function(t, shape, scale) {
  (shape/scale)*(t/scale)^(shape - 1)}

## Generate an empty parameter matrix
par <- generateParameterMatrix(hf)

## Use the vector of scales in each transition hazard
par[[1, 2]] <- list(shape = 1.3, scale = 1.1)

## Generate the cohort
cohort <- simulateCohort(transitionFunctions = hf, parameters = par,
                         cohortSize = npat, to = 30)

## Turn event times into a dataframe and make the colnames not have any spaces in them, and 
## add a status variable
gems.cohort.manual <- data.frame(cohort@time.to.state)
colnames(gems.cohort.manual) <- c("state1","state2")
gems.cohort.manual$status <- 1

head(gems.cohort.manual)


############################################################################
# Generate survival times using gems, with a pre-specified built in hazard function #
############################################################################

### Going to set up a simple model that goes from state 1, to state 2

## Generate an empty hazard matrix
hf <- generateHazardMatrix(2)

## The states named "impossible" are the ones which I need to change the name of
hf

## Define the transitions as weibull for now
## I am leaving the transition from 1 -> 3 as "impossible' to that the strcture is as planned


## Define the transitions as weibull using the prespecified built in function
hf[[1, 2]] <- "Weibull"


## Generate an empty parameter matrix
par <- generateParameterMatrix(hf)

## Use the vector of scales in each transition hazard
par[[1, 2]] <- list(shape = 1.3, scale = 1.1)

## Generate the cohort
cohort <- simulateCohort(transitionFunctions = hf, parameters = par,
                         cohortSize = npat, to = 30)

## Turn event times into a dataframe and make the colnames not have any spaces in them, and 
## add a status variable
gems.cohort.prespec <- data.frame(cohort@time.to.state)
colnames(gems.cohort.prespec) <- c("state1","state2")
gems.cohort.prespec$status <- 1

head(gems.cohort.prespec)



#####################################################
# Fit a parametric weibull model to each dataset #
#####################################################
gems.model.manual <- flexsurvreg(Surv(state2, status) ~ 1, 
                          data = gems.cohort.manual,dist = "weibull")


gems.model.prespec <- flexsurvreg(Surv(state2, status) ~ 1, 
                          data = gems.cohort.prespec,dist = "weibull")


#####################################################
# Report the models, in particular the shape and scale #
#####################################################
gems.model.manual

Call:
flexsurvreg(formula = Surv(state2, status) ~ 1, data = gems.cohort.manual, 
    dist = "weibull")

Estimates: 
       est      L95%     U95%     se     
shape  1.35790  1.33748  1.37864  0.01050
scale  1.14154  1.12429  1.15905  0.00887

N = 10000,  Events: 10000,  Censored: 0
Total time at risk: 10437.52
Log-likelihood = -9756.79, df = 2
AIC = 19517.58


gems.model.prespec

Call:
flexsurvreg(formula = Surv(state2, status) ~ 1, data = gems.cohort.prespec, 
    dist = "weibull")

Estimates: 
       est      L95%     U95%     se     
shape  1.32685  1.30668  1.34733  0.01037
scale  1.09413  1.07725  1.11128  0.00868

N = 10000,  Events: 10000,  Censored: 0
Total time at risk: 10068.98
Log-likelihood = -9499.818, df = 2
AIC = 19003.64
$\endgroup$
1
$\begingroup$

I contacted the authors of the package who responded very quickly. The problem was that when generating the data, the default number of sampler steps in the numerical approximation procedure was not high enough for this scenario (default is 1000). I upped this to 10000, and the estimates of shape and scale were accurate.

How to generate the data correctly:

## Generate the cohort
cohort <- simulateCohort(transitionFunctions = hf, parameters = par,
                         cohortSize = npat, to = 30, sampler.steps = 10000)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.