Why is the survival probability decreasing the more variables I use?

Question

I'm using coxph in combination with survfit and I observed that, e.g., when I use the following syntax

coxph(Surv(time, event) ~ 1, data = df)
coxph(Surv(time, event) ~ position, data = df)
coxph(Surv(time, event) ~ position + country, data = df)

that the survival functions move more and more leftwards as more and more variables are used, thus, decreasing in time a certain survival probability has achieved. I wonder why this as the data, especially time and event, do not change when considering more variables?!!

The data comprises ~4500 devices which have components built in described by levels of position. Furthermore, they are used across the world with 20-30 countries.

Answer 1

If you're including factor variables in your coxph model and not specifically accounting for that when you produce the estimated survival curves, the curves are probably not what you want. As @Jarle notes, the survfit function produces an estimate of the survival curve for an observation with average covariate values. This is in general not equal to the average survival curve of the whole population.

This is especially not equal to the average survival curve if there are factor variables included. Instead surfit will use the reference value of the factor variables.

Here's some example code illustrating the difference. The first plot is an estimate of the average survival curve, and it comes from the model with no covariate. The second model includes the group covariate, and the naive survival curve is significantly shifted to the left. The second plot is not an estimate of the overall survival curve (despite what the title of the plot says!), but is instead an estimate of the survival curve for the reference group. This is obvious when we plot the survival curves for both groups at the same time in plot 3, we see that plot 2 and the red line for plot 3 match exactly. Also, the estimate of the average survival curve from plot 1 is between the red and blue lines in plot 3, as expected.

library(survival)
library(survminer)
n <- 1000
group <- sample(1:2,n,replace=TRUE)
survtime <- rexp(n, rate = 0.7 - (0.3*group))
group <- factor(group)
censortime <- runif(n,0,10)
event <- censortime > survtime
time <- pmin(survtime,censortime)
dat <- data.frame(time,event,group)
mod0 <- coxph(Surv(time,event) ~ 1,
             data = dat)
mod1 <- coxph(Surv(time,event) ~ group,
              data = dat)
newdat <- data.frame(group = factor(c(1,2)))
conditional_mod1 <- survfit(mod1,newdata = newdat)
ggsurvplot(survfit(mod0),
           data = dat)

ggsurvplot(survfit(mod1),
           data = dat)

ggsurvplot(conditional_mod1,
           data=newdat)

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

Note that there is no reason the curve should shift to the left specifically when adding the factor variable, it could shift to the right if the reference group has a longer survival time than average.

EDIT: The role of the newdata is explained in the documentation of survfit.coxph,

newdata
a data frame with the same variable names as those that appear in the coxph formula. It is also valid to use a vector, if the data frame would consist of a single row.

The curve(s) produced will be representative of a cohort whose covariates correspond to the values in newdata. Default is the mean of the covariates used in the coxph fit.

So, the relevant change is to use survfit(mod, newdata = newdat) instead of survift(mod). That makes the fitted curves use the values in newdat, instead of using the average covariate value and reference level of factor variables.

Because ggsurvplot requires a data argument I also used the newdat in that function, but that's secondary to the survfit function.

Answer 2

The admonition to "read the manual" particularly applies to survival-curve estimates from Cox proportional hazards (PH) models. If you are, say, using the R survminer package, make sure that you read and understand the default settings.*

You can think about the baseline survival curve** similarly to the intercept in a linear regression. An intercept is the outcome estimate when all covariates are at their reference levels. The value of an intercept depends on how covariates are coded: the choices of reference levels for categorical covariates and whether there's been centering of continuous covariates.

This gets further complicated with Cox models. First, instead of a single intercept you have a whole survival curve as a function of time. Second, the software can make unexpected choices for "reference levels" of covariates. I'm hesitant to say anything about the details of the baseline curves you got, as they depend on what software you used (perhaps unconsciously, depending on what other functions your software invoked).

The R survfit.coxph() package reports baseline survival at some "average" values even of categorical covariates! Even the manual page says: "The resulting curve(s) almost never make sense." Once you specify particular values of covariates you get appropriate covariate adjustment for corresponding conditional survival curves, but those baseline curves aren't useful on their own. Other software packages, like rms, can provide realistic baseline curves that are more representative of the data.

So don't worry about apparent differences in default baseline survival curves. Make sure that the Cox model is appropriate for your data and then examine predictions for realistic scenarios.

---

*I infer that's one of the tools you're using, based on other questions you've posted recently. That has some strengths but I've seen a few questions here that had more to do with not understanding the survminer defaults than with survival analysis per se. Such helper packages don't necessarily use the same default settings as the underlying survival package. You need to check.

**The baseline is inferred after fitting the regression; the survfit() default is the Breslow estimate