Adjusted Survival Curves in R From Cox Model at Specific Covariate Values

Question

I'd like to plot adjusted survival curves from a Cox model at specific covariate values. The survfit function in the survival package and ggsurvplot in survminer allow one to easily plot adjusted survival curves from a model, but seem to only do so at mean values of covariates. I would like to plot curves at values that I specify, but can't find a way to easily do this in R. SAScan do this easily by using the BASELINEcommand in PROC PHREGand I'm looking to be able to do something like this in R.

Answer 1

You can do this "by hand". Here is one possibility using the example data from coxph(). Essentially, you need to make a dataset that holds constant all model variables at values you want to use and then vary time from the minimum to maximum in the data (here 0:4). Then, you can predict the survival probabilities for these data and plot.

  library(ggplot2)
  library(survival)
  test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1)) 
# Fit a stratified model 
m1 <- coxph(Surv(time, status) ~ x + strata(sex), test1) 

tmp1 <- data.frame(
  time = 0:4, 
  status=0, 
  x = 1, 
  sex = 0)

tmp1$fit <- predict(m1, newdata=tmp1, type="survival")

ggplot(tmp1, aes(x=time, y=fit)) + 
  geom_line() + 
  geom_point() + 
  theme_classic()

Below, we pick a different value of x, generate predictions and put the data together with the prediction data above. This could be done in a single step if you wanted.

tmp2 <- data.frame(
  time = 0:4, 
  status=0, 
  x = 2, 
  sex = 0)

tmp2$fit <- predict(m1, newdata=tmp2, type="survival")

tmp <- rbind(tmp1, tmp2)

Now, we could plot both sets of predictions on the same plot.

ggplot(tmp, aes(x=time, y=fit, colour=as.factor(x))) + 
  geom_line() + 
  geom_point() + 
  theme_classic() + 
  labs(colour="X")

^{Created on 2022-04-29 by the reprex package (v2.0.1)}