Predict out of sample using flexsurvreg in R

user_6043168 0
Viewed 2469

I have the following model in R 

library(flexsurv)

data(ovarian)

model = flexsurvreg(Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, dist='weibull')

model

predict(model,data = ovarian, type = 'response')

The model summary looks like this flexsurvreg model output

I am trying to predict the survival time using the predict function in R and get the following error

error while trying to predict

How can I predict expected lifetime using this flexsurvreg model?

I understand that the documentation mentions a totlos.fs function, but this data does not seem to have a trans variable that totlos.fs requires to provide an output.

If there is no other alternative to totlos.fs how can I create a trans variable in this data and handle it along with existing covariates?

Please advise.

survival-analysis predict r
2 Answers

Nik,

I know your question is an old one, but see below how I hacked a way to do it. It involves retrieving the shape and rate parameters from your fit of test data, then instead of predict, you use the qgompertz() from flexsurv. Please excuse the use of my own encapsulated example code, but you should be able to follow along.

# generate the training data "lung1" from data(lung) in survival package
# hacked way for truncating the lung data to 2 years of follow up
require(survival)
lung$yrs <- lung$time/365
lung1 <- lung[c("status", "yrs")]
lung1$status[ lung1$yrs >2] <- 1
lung1$yrs[ lung1$yrs >2]  <- 2

# from the training data build KM to obtain survival %s
s <- Surv(time=lung1$yrs, event=lung1$status)
km.lung <- survfit(s ~ 1, data=lung1)
plot(km.lung)

# generate dataframe to use later for plotting 
cut.length <- sum((km.lung$time <= 2)) # so I can create example test data
test.data <- data.frame(yrs = km.lung$time[1:cut.length] , surv=round(km.lung$surv[1:cut.length], 3))


##
##  doing the same as above with gompertz
##
require(flexsurv) #needed to run gompertz model
s <- Surv(time=lung1$yrs, event=lung1$status)
gomp <- flexsurvreg(s ~ 1, data=lung1, dist="gompertz") # run this to get shape and rate estimates for gompertz
gomp # notice the shape and rate values 

# create variables for these values
g.shape <- 0.5866
g.rate <- 0.5816


##
##  plot data and vizualize the gomperts
##
# vars for plotting
df1 <- test.data 
xvar <- "yrs"
yvar <- "surv"

extendedtime <- 3 # 
ylim1 <- c(0,1)
xlim1 <- c(0, extendedtime)

# plot the survival % for training data
plot(df1[,yvar]~df1[,xvar], type="S", ylab="", xlab="", lwd=3, xlim=xlim1, ylim=ylim1)
# Nik--here is where the magic happens... pay special attention to: qgompertz(seq(.01,.99,by=.01), shape=0.58656, rate = .5816) 
lines (qgompertz(seq(.01,.99,by=.01), shape=0.58656, rate = .5816) ,  seq(.99,.01,by=-.01) , col="red", lwd=2, lty=2  )

# generate a km curve from the testing data
s <- Surv(time=lung$yrs, event=lung$status)
km.lung <- survfit(s ~ 1, data=lung)
par(new=T)
# now draw remaining survival curve from the testing section
plot(km.lung$surv[(cut.length+1):length(km.lung$time)]~km.lung$time[(cut.length+1):length(km.lung$time)], type="S", col="blue", ylab="", xlab="", lwd=3, xlim=xlim1, ylim=ylim1)
user_12846170
Related