How to generate covariate-adjusted cox survival/hazard functions?

Question

I'm using the survminer package to try to generate survival and hazard function graphs for a longitudinal student-level dataset that has 5 subgroups of interest.

I've had success creating a model that shows the survival functions without adjusting for student-level covariates using ggsurvplot.

ggsurvplot(survfit(Surv(expectedgr, sped) ~ langstatus_new, data=mydata), pvalue=TRUE)

Output example

However, I cannot manage to get these curves adjusted for covariates. My aim is to create graphs like these. As you can see, these are covariate-adjusted survival curves according to some factor variable. Does anyone how such graphs can be obtained in R?

Answer 1

Although correct, I believe that the method described in the answer of Dion Groothof is not what is usually of interest. Usually, researchers are interested in visualizing the causal effect of a variable adjusted for confounders. Simply showing the predicted survival curve for one single covariate combination does not really do the trick here. I would recommend reading up on confounder-adjusted survival curves. See https://arxiv.org/abs/2203.10002 for example.

Those type of curves can be calculated in R using the adjustedCurves package: https://github.com/RobinDenz1/adjustedCurves

In your example, the following code could be used:

library(survival)
library(devtools)

# install adjustedCurves from github, load it
devtools::install_github("/RobinDenz1/adjustedCurves")
library(adjustedCurves)

# "event" needs to be binary
lung$status <- lung$status - 1

# "variable" needs to be a factor
lung$ph.ecog <- factor(lung$ph.ecog)

fit <- coxph(Surv(time, status) ~  ph.ecog + age + sex, data=lung,
             x=TRUE)

# calculate and plot curves
adj <- adjustedsurv(data=lung, variable="ph.ecog", ev_time="time",
                    event="status", method="direct",
                    outcome_model=fit, conf_int=TRUE)
plot(adj)

Producing the following output:

These survival curves are adjusted for the effect of age and sex. More information on how this adjustment works can be found in the documentation of the adjustedCurves package or the article I cited above.