<p>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 <em>adjusted</em> 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 <a href="https://arxiv.org/abs/2203.10002" rel="nofollow noreferrer">https://arxiv.org/abs/2203.10002</a> for example.</p>
<p>Those type of curves can be calculated in R using the <code>adjustedCurves</code> package: <a href="https://github.com/RobinDenz1/adjustedCurves" rel="nofollow noreferrer">https://github.com/RobinDenz1/adjustedCurves</a></p>
<p>In your example, the following code could be used:</p>
<pre><code>library(survival)
library(devtools)

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

# &quot;event&quot; needs to be binary
lung$status &lt;- lung$status - 1

# &quot;variable&quot; needs to be a factor
lung$ph.ecog &lt;- factor(lung$ph.ecog)

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

# calculate and plot curves
adj &lt;- adjustedsurv(data=lung, variable=&quot;ph.ecog&quot;, ev_time=&quot;time&quot;,
                    event=&quot;status&quot;, method=&quot;direct&quot;,
                    outcome_model=fit, conf_int=TRUE)
plot(adj)
</code></pre>
<p>Producing the following output:</p>
<p><a href="https://i.stack.imgur.com/UkZp7.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/UkZp7.png" alt="Output" /></a></p>
<p>These survival curves are adjusted for the effect of <code>age</code> and <code>sex</code>. More information on how this adjustment works can be found in the documentation of the <code>adjustedCurves</code> package or the article I cited above.</p>
How to generate covariate-adjusted cox survival/hazard functions?
<p>I'm using the <code>survminer</code> package to try to generate survival and hazard function graphs for a longitudinal student-level dataset that has 5 subgroups of interest.</p>
<p>I've had success creating a model that shows the survival functions <em>without</em> adjusting for student-level covariates using <code>ggsurvplot</code>.</p>
<pre><code>ggsurvplot(survfit(Surv(expectedgr, sped) ~ langstatus_new, data=mydata), pvalue=TRUE)
</code></pre>
<p><a href="https://i.stack.imgur.com/6r3ze.png" rel="nofollow noreferrer">Output example</a></p>
<p>However, I cannot manage to get these curves adjusted for covariates. My aim is to create <a href="https://i.stack.imgur.com/OMaga.jpg" rel="nofollow noreferrer">graphs like these</a>. As you can see, these are covariate-adjusted survival curves according to some factor variable. Does anyone how such graphs can be obtained in <code>R</code>?</p>