Estimating survival curves from Cox regression results

Question

I can understand that it is possible to estimate survival curves directly from the results of a Cox regression.

The way it can be done, mathematically, is furthermore very nicely explained in this video: https://www.coursera.org/lecture/simple-regression-analysis-public-health/estimating-survival-curves-from-cox-regression-results-iQRtN

I am currently trying to learn how to do this manually in R, but I have not been able to make it work so far. I was therefore hoping that some of you might be able to understand the mathemathics involved as well as being able to translate the proces into R code. I have tried different approaches but without any luck, so any help or input is greatly appreciated.

I am aware that there exist questions on here which share some similarities with this one, but they do not follow the approach in the video, and I do not believe that they answer the problem presented here, which is why I created this question.

As an example, I have created this Cox model

library(survival)
data("nwtco")
subcoh <- nwtco$in.subcohort
selccoh <- with(nwtco, rel==1|subcoh==1)
ccoh.data <- nwtco[selccoh,]
ccoh.data$subcohort <- subcoh[selccoh]
## central-lab histology
ccoh.data$histol <- factor(ccoh.data$histol,labels=c("FH","UH"))
## tumour stage
ccoh.data$stage <- factor(ccoh.data$stage,labels=c("I","II","III","IV"))
ccoh.data$age <- ccoh.data$age/12 # Age in years

fit.ccSP <- cch(Surv(edrel, rel) ~ stage + histol + age, data =ccoh.data,
                subcoh = ~subcohort, 
                id=~seqno, 
                cohort.size=4028, 
                method="LinYing")

from which I get the following results

> summary(fit.ccSP)
Case-cohort analysis,x$method, LinYing 
 with subcohort of 668 from cohort of 4028 

Call: cch(formula = Surv(edrel, rel) ~ stage + histol + age, data = ccoh.data, 
    subcoh = ~subcohort, id = ~seqno, cohort.size = 4028, method = "LinYing")

Coefficients:
          Coef    HR  (95%   CI)     p
stageII  0.693 1.999 1.453 2.751 0.000
stageIII 0.627 1.872 1.348 2.599 0.000
stageIV  1.300 3.668 2.529 5.320 0.000
histolUH 1.458 4.299 3.240 5.704 0.000
age      0.046 1.047 1.002 1.094 0.039

Based on these results, I am trying - as shown in the video - to estimate the survival curves for the different stages, i.e. "stageII", "stageIII" and "stageIV", but, as mentioned, I have not been able to make it work.

Answer 1

In general, survival as a function of time, $S(t)$, is conveniently expressed in terms of the cumulative hazard function $H(t)$:

$$ -\log S(t) = H(t) $$

where (trying to match the symbols in the linked video) $H(t) = \int_0^t \lambda(\tau) d\tau$ is the cumulative hazard function, integrating through time $t$ the instantaneous hazard function $\lambda(\tau)$. Applying this to your situation with a case-cohort Cox model (cch() function) requires establishing two additional facts.

First, the linked video only gives an equation for estimated survival as a function of time $t$ and covariate values $x$, $\hat S(t;x)$, from a Cox model when you know the baseline hazard function at some set of covariate reference values, $\lambda_0(t)$. It doesn't show how to estimate the baseline $\lambda_0(t)$ or $H_0(t)$, seeming to imply that one uses the survival curve restricted to a reference group for that.

This answer shows one standard way to estimate the baseline hazard from a Cox model, which includes information from all cases having events whether in a baseline reference group or not. The Breslow estimate extends the Nelson-Aalen estimate of cumulative hazard to include the covariate values in a Cox model:

$$\hat{H}_0(t)=\int_0^t\frac{dN(\tau)}{\sum_i \exp(\beta^T x_i)Y_i(\tau)},$$

where $Y_i(\tau)$ is the at-risk indicator (=1 for cases at risk at time $\tau$, 0 otherwise), $\beta$ is the vector of Cox regression coefficient estimates, $x_i$ is the vector of covariate values for case $i$, and $dN(\tau)$ is the number of events at time $\tau$. This is just a weighted cumulative sum of event numbers through time $t$, with the weight for events at each time the inverse of the summed relative risks of cases then at risk. You are welcome to calculate this baseline hazard "manually" if you wish.

Then for a particular new set of covariates $x_j$ you can calculate:

$$- \log\hat S(t;x_j) = \hat{H}_0(t) \exp (\beta^T x_j)$$

If you have a coxph model in R, this is all more conveniently done with the survfit.coxph() function, which can provide a different choice for estimating the baseline hazard (one that extends the Kaplan-Meier product-limit survival estimate instead of the Nelson-Aalen estimate) and 2 ways to handle tied event times, and plotting the result (with confidence limits, if desired).

Second, the case-control analysis provide by cch() doesn't seem to provide a function to extract the information you need to do this. Prentice proposed this form of survival analysis in 1986 to deal with situations in which events are rare or it's difficult/expensive to collect the required covariate values. You start with a defined cohort that is followed throughout the study for occurrence of the event, but you restrict covariate data at first to a sub-cohort. "Cases" are defined as those who experience events. If you evaluate survival only within the initial sub-cohort, that's just standard survival analysis on that subset of the total cohort.

The difference in the case-cohort analysis is that, when events occur in cohort members outside of the initial sub-cohort, you collect the covariate data on them and add them as "cases" to the survival analysis. This case-cohort analysis deliberately overweights the "cases" in the data collection, which must be taken into account in terms of bias and variance in the estimated Cox regression coefficients. The cch() function is designed to do that, but its output isn't a standard coxph object (although its inner workings use the coxph() function). So far as I can tell, it doesn't provide a pre-defined function to get the information needed to generate a baseline hazard or survival-curve estimates.

Prentice recommended an adaptation of the Breslow-Nelson-Aalen cumulative-hazard estimate for case-cohort analysis:

$$\hat{H}_0(t)= \frac{m}{n} \int_0^t\frac{d\bar N(\tau)}{\sum_{i\in C} \exp(\beta^T x_i)Y_i(\tau)},$$

where $C$ represents the initial sub cohort, $m$ is the initial sub-cohort size, $n$ is the total cohort size, and $\bar N$ includes all events within the entire cohort. That is, the weight for events at time $\tau$ is limited to cases within the initial sub cohort still at risk, whether or not an event occurs outside the initial sub-cohort. That overcomes a major source of bias, but the precision of the estimate of the baseline hazard is limited by the size of the initial sub-cohort. For getting this baseline hazard from a cch() model and estimating survival curves for specified covariate values, you might have to proceed "manually."