Visually Comparing the Kaplan-Meier Curve to the Cox PH Model Curve

Question

I am conducting a survival analysis and have a few questions regarding it's interpretation with respect to the Cox Proportional Hazards Model:

Why does the inclusion of different covariates change the shape of the plotted survival function? When I compare the Kaplan-Meier survival curve with a plotted survival curve from the cox model, the shape often looks drastically different. Furthermore, the cox survival curve shape changes most when I add in covariates that have the largest hazard ratio of the generated model. Obviously these factors are related, but how?

I apologize if this is a simple question but I am new to survival analysis (and statistics in general!). If beneficial, i'd be happy to supply some R code that I am working with to increase comprehension.

Answer 1

A survival curve generated from the Cox model estimates 1 baseline hazard for all combinations of covariates in the model unless stratification is used. The instantaneous hazard for each group or covariate value combination is specified as: h(t)*hazard ratio, where h(t) is the instantaneous hazard at any specified time, t. The instantaneous hazard can be transformed to a cumulative hazard and a survival estimate. By virtue of this calculation, the KM curves, which are empirical, will not be the same as the Cox survival estimates. The Cox hazard estimates and survival estimates are multiples of the baseline hazard function, while the KM estimates are unrestricted and can take any form.

Some R code that may help:

library(survival)
library(ggplot2)

surv.obj <- with(jasa, Surv(futime, fustat ))
surg.fit <- survfit(surv.obj ~ surgery, data = jasa)
km.surg <- data.frame(time = surg.fit$time, surv = surg.fit$surv, 
                      surg = c(rep(names(surg.fit$strata[1]), surg.fit$strata[1]), rep(names(surg.fit$strata[2]), surg.fit$strata[2]) )
                      )
ggplot(km.surg, aes(time, surv, group = surg)) + geom_step()

surg.cox <- coxph(surv.obj ~ surgery, data = jasa)
base.haz <- rbind(basehaz(surg.cox), basehaz(surg.cox) )
base.haz$surg <- c(rep(0,88),rep(1,88))
base.haz$haz.hr <- exp(base.haz$surg * surg.cox$coefficients[1])
base.haz$surv <- exp(-base.haz$haz.hr*base.haz$hazard)

ggplot(base.haz, aes(time, surv, group = surg )) + geom_step()