Closed form expression for survfit function used in R (Cox proportional hazard survival model)

Question

I have a dataset of patient survival times that I am trying to predict based on n covariates X1...Xn (patient co-morbidities, other diagnoses, lab results) using a Cox proportional hazard model. I'm using the survival and survminer libraries in R to implement the model.

I have

test.cox = coxph(Surv(Elapsed_Time,outcome) ~ Hypertension + Diabetes + Age + Gender, data = StageMat)

where outcome = 1,0 for death/survival, I get for coefficients (making some numbers up) 10 hypertension, 8 diabetes, 0.1 Age, 0.8 GenderMale, i.e. b1 = 10, b2 = 8, b3 = 0.1, b4 = 0.8. This corresponds to a hazard model

h(t) = h0(t)exp[b1*x1 + b2*x2 + b3*x3 + b4*x4]

where h0(t) is an (unknown) baseline hazard function assuming mean values for each covariate. If I then use these estimated values to fit new data, I implement

test.fit = survfit(test.cox, newdata = eval_data[1,])

My question is, what is the closed form expression that R uses to calculate time-dependent probabilities in survfit given that h0(t) is an unknown? There is an approximation based on the mean of the hazard function (or the hazard function of the mean) that I can't seem to get at.

The general form of a survival function is

 S(T) = Exp[- Integral_0_T h(t)]

where Integral_0_T h(t) is the definite integral of h(t) from 0 to T

I would like to write down a closed form exponential function for S(T) given b1...b4 from my model, but because h0(t) is unknown, I'm not sure what this expression would be.

For a proportional hazard model, S(T) can be rewritten as

S(T) = S0(T)^Exp[Sum_i bi*Xi]

but as for h0(t), the baseline S0(T) is unknown. I've seen some papers where S0(t) is replaced by the mean value of S(T) in the dataset to get a closed-form expression, but I don't think this is generally correct because the mean of a function is not necessarily equal to the function of the mean.

Any clarification of how to deal with baseline estimates to get a closed-form expression for the survival (and hazard) functions would help.

Answer 1

The (Aalen-Breslow) estimator of the baseline cumulative hazard is

$$\hat\Lambda_0(t)=\sum_{s\leq t} \frac{dN(s)}{\sum_{i}Y_i(s)\exp(Z_i(s)\beta)}$$ where $dN(s)$ is the number of deaths at time $s$, $Y_i(s)$ is 1 if individual $s$ is under observation at time $s$, and $Z_i(s)$ is the covariate vector for individual $i$ at time $s$. The sum is over times that an event occurs, because $dN$ is zero at other times

Heuristically (according to Breslow), it comes from Aalen's cumulative hazard estimator $$\hat\Lambda(t)=\sum_{s\leq t} \frac{dN(s)}{\sum_{i}Y_i(s)}$$ but in the Cox model each individual contributes $e^{Z\beta}$ rather than 1 to the denominator.

In fact, what's used in the computation is the increments of $\hat\Lambda$ $$d\hat\Lambda_0(t)=\frac{dN(s)}{\sum_{i}Y_i(s)\exp(Z_i(s)\beta)}$$

The survival for individuals is estimated by first computing the cumulative hazard for individuals $$\Lambda(t;z)=\sum_{s\leq t} \exp(Z_i(s)\beta)\, d\hat\Lambda_0(t)$$ and then using $S(t)=e^{-\Lambda(t)}$.

Finally, if you want the predicted survival curve for a cohort, you compute $S(t;z)$ for each individual and then average them. As you note, this is not the same as the survival curve for an average individual, though in practice the difference is typically not all that big.