0

I am trying to calculate the restricted mean survival time grouped by the covariates but I get an error. Could you please have a look at the code example code below and let me know what the problem is?

library(survival)
library(flexsurv)

data <- survival::lung
formula <- Surv(data$time, data$status) ~data$sex
fsmodel_exponential<-flexsurvreg(formula,dist = "exp")

#Produces the expected results
rate_exponential<-fsmodel_exponential$res[1]
rmst_exp <- rmst_exp(t = 30, rate = rate_exponential, start = 0)
rmst_exp

#Gives error for the covariate
rate_exponential_sex<- fsmodel_exponential$res[2]
rmst_exp2 <- rmst_exp(t = 30, rate = rate_exponential_sex, start = 0)
rmst_exp2
gibsonmate
  • 25
  • 4

1 Answers1

1

Your fsmodel_exponential$res[2] is negative, which causes an error. In the exponential model covariates on these parameters represent an accelerated failure time (AFT) model. You do get those when you invert the factor levels of sex, which results in positive maximum likelihood estimates.

library(survival)
library(flexsurv)
data <- lung
data$sex <- relevel(factor(data$sex, labels = c("M", "F")), "F")
formula <- with(data, Surv(time, status) ~ sex)
fsmodel_exponential <- flexsurvreg(formula, dist = "exp")
rmst_exp(t = 30, rate = fsmodel_exponential$res[1], start = 0)
#> [1] 29.23162
rmst_exp(t = 30, rate = fsmodel_exponential$res[2], start = 0)
#> [1] 1.998406

plot(fsmodel_exponential, col = c("blue", "red"), 
    lwd.obs = 2, xlab = "months", ylab = "Recurrence-free survival")
legend("topright", levels(data$sex), col=c("red", "blue"), lty=1)

Edit: To get the RMST, simply run:

summary(fsmodel_exponential, type = "rmst", t= 30 )
#> sex=M 
#>   time     est      lcl      ucl
#> 1   30 28.7467 28.49528 28.94901
#> 
#> sex=F 
#>   time      est      lcl      ucl
#> 1   30 29.23162 28.98571 29.40346
user12728748
  • 8,106
  • 2
  • 9
  • 14
  • Thanks for your reply! However, I'd like to know if there is a way to calculate the RMST for both curves at the same time? E.g. 29.23 for the females in your example, but how much is it for men? – gibsonmate Feb 26 '20 at 17:17
  • @gibsonmate: see edit in my answer - `summary(fsmodel_exponential, type = "rmst", t= 30)` would give it for `t=30`, for example. See `?summary.flexsurvreg` for more details. – user12728748 Feb 26 '20 at 20:35