CODEWITHSUNDEEP
rplotsurvival-analysiscox-regression
stats_noob
stats_noob

Reputation: 5897

Predicting and Plotting Survival Curve with the CoxPH

I am trying to predict and plot the (estimated) survival curve for a new observation in R. Using the "survival" library and the "lung" data set, I first fit a cox proportional hazards model to the data. Then, I tried to predict and plot the survival curve for a hypothetical new observation (I entered the details for this hypothetical new observation in the "list" command). However, this is not working.

I have attached my code below:

#load library
library(survival)
data(lung)

#create survival object
s <- with(lung,Surv(time,status))

#create model
modelA  <- coxph(s ~ as.factor(sex)+age+ph.ecog+wt.loss+ph.karno,data=lung, model=TRUE)
summary(modelA)

#plot
plot(survfit(modelA), ylab="Probability of Survival",
     xlab="Time", col=c("red", "black", "black"))

#predict for a hypothetical NEW observation (here is where the error is)
lines(predict(modelA, newdata=list(sex=1, 
                                  age = 56, 
                                  ph.ecog = 1, 
                                  ph.karno = 50,
                                  wt.loss = 11),
               type="quantile",
               p=seq(.01,.99,by=.01)),
       seq(.99,.01,by=-.01),
       col="blue")
## Error in match.arg(type) : 
##  'arg' should be one of “lp”, “risk”, “expected”, “terms”, “survival”

Does anyone know what I am doing wrong? Thanks

Upvotes: 4

Views: 5198

Answers (2)

Allan Cameron
Allan Cameron

Reputation: 173793

This is what the survfit function is for. In your example, you plot the survfit for the model, but you can feed a newdata argument into this function and it will produce the estimated survival for these data.

If we reproduce your example:

library(survival)

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

modelA  <- coxph(s ~ as.factor(sex) + age + ph.ecog + wt.loss + ph.karno,
                 data = lung, model = TRUE)

plot(survfit(modelA), ylab = "Probability of Survival",
     xlab = "Time", col = c("red", "black", "black"))

Then we can create a survival curve given your specified covariates like this:

est <- survfit(modelA, newdata = data.frame(sex      = 1,
                                            age      = 56,
                                            ph.ecog  = 1, 
                                            ph.karno = 50,
                                            wt.loss  = 11))

Now est is an S3 object with members that include time and survival, so we can plot a blue line tracking the estimated survival of individuals with the given covariates like this:

lines(est$time, est$surv, col = 'blue', type = 's')

enter image description here

Or plot it on its own with a 95% confidence interval:

plot(est, ylab = "Probability of Survival",
     xlab = "Time", col = c("red", "black", "black"))

enter image description here

Created on 2022-05-26 by the reprex package (v2.0.1)

Upvotes: 4

Vasily A
Vasily A

Reputation: 8626

See the description of the predict() function (you can open it in R help by running ?predict.coxph, or here for example):

type - the type of predicted value. Choices are the linear predictor ("lp"), the risk score exp(lp) ("risk"), the expected number of events given the covariates and follow-up time ("expected"), and the terms of the linear predictor ("terms"). The survival probability for a subject is equal to exp(-expected).

You can see that your type="quantile" does not match expected input. If you call predict() without the type argument, in your case it will default to using lp (linear predictor).
When you call predict() function for your object modelA, it determines that it is of coxph class, so the predict.coxph() function is applied. The arguments like type="quantile" and p=seq(.01,.99,by=.01) are not acceptable for predict.coxph() (p is ignored, type raises error). They are used in another function, predict.survreg() - for it to be called, your modelA object must be of survreg class, i.e. it should be created using survreg() call instead of coxph() call.

Upvotes: 2

Related Questions