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([email protected])
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([email protected])
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
