Search code examples
rsurvival-analysis

predicted survival curve in R- Parametric method


I am trying to create predicted survival plots for survreg or flexsurvreg in R. But I am getting errors for the plot when I use more than one predictors in the survreg. I would like to try it using flexsurvreg or survreg, For lung data set, I used the following code to fit the model.

require(survival)
s <- with(lung,Surv(time,status))

sWei  <- survreg(s ~ as.factor(sex)+age+ph.ecog+wt.loss+ph.karno,dist='weibull',data=lung)

fitKM <- survfit(s ~ sex,data=lung)
plot(fitKM)

lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="blue")
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

I am getting errors when I use the above command for plotting. Please let me know where I am doing wrong when plotting for predicted survival curve.

> lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
Error in eval(expr, envir, enclos) : object 'age' not found

Solution

  • We need to give values in the list to every variable in your model, so that it can plot the curve.

    require(survival)
    s <- with(lung,Surv(time,status))
    
    sWei  <- survreg(s ~ as.factor(sex)+age+ph.ecog+wt.loss+ph.karno,dist='weibull',data=lung)
    
    fitKM <- survfit(s ~ sex,data=lung)
    plot(fitKM)
    
    lines(predict(sWei, newdata=list(sex=1, 
                                     age = 1, 
                                     ph.ecog = 1, 
                                     ph.karno = 90,
                                     wt.loss = 2),
                                     type="quantile",
                                     p=seq(.01,.99,by=.01)),
                                     seq(.99,.01,by=-.01),
                                     col="blue")
    lines(predict(sWei, newdata=list(sex=2, 
                                     age = 1, 
                                     ph.ecog = 1, 
                                     ph.karno = 90,
                                     wt.loss = 2),
                                     type="quantile",
                                     p=seq(.01,.99,by=.01)),
                                     seq(.99,.01,by=-.01),
                                     col="red")
    

    enter image description here