3
$\begingroup$

I am hoping to plot the probability of a patient surviving to 1 year without progression (status==0) vs. the value of a continuous variable (radiation.dose). This would be based on a univariate cox proportional hazards model. I've Googled this and can't find anyone who has done this before. What I imagine is similar axes to this logistic regression plot: Logistic probability vs. variable

Here is some sample code for you to work with:

library("survival")
require("survival")

days <- rpois(100, 365)
status <- rbinom(100,1,0.34)
radiation.dose <- sapply(status,function(x){ifelse(x==0,rnorm(1,80,20),rnorm(1,60,20))})

df<-data.frame(days,status,radiation.dose)
mod<-Surv(df$days,df$status)
cox.mod<-coxph(mod~radiation.dose,data=df)

Thank you for your help!

Improve this question
$\endgroup$

2 Answers 2

Reset to default
0
$\begingroup$

I'm just thinking out loud here, I'm posting this as an answer instead of a comment for better formatting.

According to the help file on predict.coxph, "The survival probability for a subject is equal to exp(-expected)." So we can predict the survival probability

cox.pred<- exp(-predict(cox.mod, type = "expected"))
plot(radiation.dose, cox.pred)

Predicted Probability of Survival

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for thinking through this with me! I think this looks like the right function. I'm far from a statistics expert, but I would have thought that since the predicted probabilities are based on an underlying mathematical model there wouldn't be so much variance in the probabilities. I guess I was expecting what looks like a curve sort of like that logistic regression plot I linked to. Does that make sense? What are your thoughts? $\endgroup$
    – JJM
    Commented Dec 18, 2019 at 19:02
  • $\begingroup$ The way the data set is created means both survivors and non-survivors have a 1 sd overlap of treatment between 60 and 80. Status doesn't depend on days, so why are you doing cox at all? My understanding: There are 100 patients. Each patient receives a treatment dose between 0 and 140. Patients then makes a single followup at ~365 days and are either in remission (status = 1) or relapsed (status = 0). Higher doses are more likely to cause remission. I am assuming relapse patients that present too early are likely to be misdiagnosed as in remission and that is how days should be incorporated? $\endgroup$
    – Ron Jensen
    Commented Dec 18, 2019 at 19:59
0
$\begingroup$

The answer given by @Ron Jensen is incorrect. By simply calling the predict function without specifying the newdata argument, you are predicting the survival probability given the radiation.dose and the observed value of time. You can easily see this because there is no clear relationship in your plot, which has to be the case since you used a cox-regression model.

If you want to plot the survival probability at t as a function of a continuous variable based on a cox-model, you can instead use the contsurvplot package (https://cran.r-project.org/package=contsurvplot) I created. Using your example, you can do it like this:

library(survival)
library(contsurvplot)

set.seed(42)

days <- rpois(100, 365)
status <- rbinom(100, 1, 0.34)
radiation.dose <- sapply(status,function(x){ifelse(x==0, rnorm(1,80,20), rnorm(1,60,20))})

df <- data.frame(days,status,radiation.dose)
mod <- cox.mod<-coxph(Surv(days, status) ~ radiation.dose, data=df, x=TRUE)

plot_surv_at_t(time="days",
               status="status",
               variable="radiation.dose",
               data=df,
               model=cox.mod,
               t=365)

contsurvplot

By specifying t = 365, you're telling the plot_surv_at_t function to use the survival probability at 1 year.

Since there are no other independent variables in this cox model, this is equivalent to using:

vals <- seq(min(df$radiation.dose), max(df$radiation.dose), 1)

newdata <- data.frame(radiation.dose=vals,
                      days=365,
                      status=0)

p <- predict(cox.mod, type="surv", newdata=newdata)

plot(vals, p)

base_r

If there is no specific reason to use one value of time, I would however recommend to use survival area plots or survival contour plots instead, as described in detail in this preprint: https://arxiv.org/pdf/2208.04644.pdf

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.