Unable to calculate the shape and scale parameters of weibull distribution from data simulated through either gems or simsurv in R

Question

I am trying to simulate survival data from a weibull distribution with shape = 1.3 and scale = 1.1. I am then fitting a weibull survival model to the data, to try and calculate back the values of 1.3 and 1.1. I am doing this in R.

To simulate the data, the packages I have tried to use to generate the survival times are gems and simsurv. I want to use these packages, as opposed to rweibull, as eventually I will be simulating for more complex multistate data, and the hazard of the survival function must also be dependent on baseline covariates, etc. For comparison, I have also simulated survival times using rweibull.

To analyse the data, I am using the flexsurvreg function from the package flexsurv.

When I simulate data using gems, both the shape and scale are calculated incorrectly when fitting a weibull model back to the data. When I simulate data using simsurv, the shape is calculated correctly, but the scale is calculated incorrectly. When I simulate data using rweibull, both the shape and scale are estimated correctly. This final result makes me think it is the data that is being generated incorrectly, rather than flexsurvreg function not working properly.

I can see no reason why this should be the case, does anybody have an explanation for this? Does anybody have any other suggestions for generating survival data in R that allows different distributions to be defined and hazards to be adjusted at baseline?

A reproducible example is below. Note the packages will need to be installed if you do not already have them, which can be done by removing the # on the first four lines of code.

#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 #
#####################################################

## 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 structure is as planned


## Define the transitions as weibull
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 <- data.frame([email protected])
colnames(gems.cohort) <- c("state1","state2")
gems.cohort$status <- 1

head(gems.cohort)


#####################################################
# Generate survival times using simsurv #
#####################################################

## Creaet an empty dataset of baseline variables (used to define number of observations also)
bl  <- data.frame(id = 1:npat)

## Generate the data using simsurv (note that lambda is actually defined as 1/lambda)
simsurv.data <- simsurv(lambdas = (1/1.1), gammas = 1.3, x = bl, maxt = 30)
head(simsurv.data)



#####################################################
# Generate survival times using rweibull #
#####################################################

## Generate directly from rweibull
rweibull.data <- data.frame("eventtime" = rweibull(npat, shape = 1.3, scale = 1.1), "status" = rep(1, npat))



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



simsurv.model <- flexsurvreg(Surv(eventtime, status) ~ 1, 
                           data = simsurv.data,dist = "weibull")



rweibull.model <- flexsurvreg(Surv(eventtime, status) ~ 1, 
                             data = rweibull.data,dist = "weibull")


#####################################################
# Report the shape and scale #
#####################################################
gems.model
simsurv.model
rweibull.model

# Seems odd, but need to exponentiate the $coefficients to get the values reported from the model output
# This is odd, because these are not multiplicative effects on the baseline hazard
exp(gems.model$coefficients)
exp(simsurv.model$coefficients)
exp(rweibull.model$coefficients)

> exp(gems.model$coefficients)
   shape    scale 
1.357903 1.141540 
> exp(simsurv.model$coefficients)
   shape    scale 
1.288051 1.073805 
> exp(rweibull.model$coefficients)
   shape    scale 
1.302091 1.100632 

Answer 1

The issue when using simsurv is that the scale is actually defined like this: scale_sm = 1/(scale_rw^shape_sm), where _sm refers to the shape/scale according to simsurv, and _rw refers to the shape/scale according to rweibull or flexsurvreg.

Using the following code to generate the data:

simsurv.data <- simsurv(lambdas = 1/(1.1^1.3), gammas = 1.3, x = bl, maxt = 30) 

The issue when using the gems package is answered in this question: Generating weibull survival times using the built in rweibull function vs manually defining the hazard (using the gems package in R)