Is it possible and how to predict individual survival curve after Cox regression?

Question

Taking the veteran dataset of a two-treatment, randomized trial for lung cancer in the R package survival as an example, where

• time is the survival time in days
• status is the censoring status (0 for surviving and 1 for dead)
• trt is the treatment type (1 or 2)
• celltype: 1=squamous, 2=small cell, 3=adeno, 4=large
• age: age in years

(There are more covariates but for simplicity I only included 3)

Adding in the above covariates/factors, the Cox regression can be run as follows:

library(survival)
cox <- coxph(Surv(time, status) ~ trt + celltype + age, data = veteran)
summary(cox)
Call:
coxph(formula = Surv(time, status) ~ trt + celltype + age, data = veteran)
n= 137, number of events= 128
coef exp(coef) se(coef)     z Pr(>|z|)
trt               0.179011  1.196034 0.201404 0.889    0.374
celltypesmallcell 1.080310  2.945592 0.274647 3.933 8.37e-05 ***
celltypeadeno     1.170470  3.223506 0.294727 3.971 7.15e-05 ***
celltypelarge     0.292624  1.339939 0.285504 1.025    0.305
age               0.004097  1.004106 0.009581 0.428    0.669
---
Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
              exp(coef) exp(-coef) lower .95 upper .95

trt                   1.196     0.8361    0.8059     1.775
celltypesmallcell     2.946     0.3395    1.7195     5.046
celltypeadeno         3.224     0.3102    1.8091     5.744
celltypelarge         1.340     0.7463    0.7657     2.345
age                   1.004     0.9959    0.9854     1.023
Concordance= 0.619  (se = 0.028 )
Likelihood ratio test= 26.04  on 5 df,   p=9e-05
Wald test            = 25.01  on 5 df,   p=1e-04
Score (logrank) test = 26.51  on 5 df,   p=7e-05

Now, these coefficients allow me to compute the hazard ratio between any two groups of patients, each with a set of trt, celltype, and age. For examples, if

• group A has trt = 1, celltype = squamous, and age = 49, and
• group B has trt = 2, celltype = smallcell, and age = 50,

then we can say that at any time point, the proportion of subjects from group B who have died from group B should be $1.196\times2.946\times1.004=3.538$ times of that from group A. But this does not allow us to tell the proportion of survival at any time point for either group A or B - we only know the ratio.

In broader term, is it possible to tell at any given time point, for a group of patients with known trt, celltype, and age, the proportion of survival? If all terms tested are factors, we can just subset the target population from the dataset and generate the Kaplan-Meier curve, but then age is present in the model. Maybe a relevant question is: what is the baseline hazard function $h_0(t)$ in this case?

Answer 1

The baseline hazard rate $h_0(t)$ is hard to estimate (takes smoothing), but the baseline cumulative hazard $H_0(t)$ is easy. Given $H_0(t)$ you can estimate an individual survival curve by $$\log S(t;z) = H(t;z) = H_0(t)e^{\beta z}$$

Basically any survival analysis software will both let you get individual survival estimates and also compute the Breslow-Aalen estimator of $H_0$ $$\hat{H}_0(t)=\int_0^t\frac{dN(s)}{\sum_i e^{\beta z_i}Y_i(s)}$$ This is like the Nelson-Aalen estimator of the cumulative hazard for a group, except that the denominator is not the number at risk but a weighted number at risk using the individual hazard ratios as weights.

You do need to check whether your specific software defines $H_0(t)$ as the cumulative hazard when all covariates are zero or (more conveniently) as the cumulative hazard when all covariates are equal to the baseline averages. In R, you want predict.coxph.