6
$\begingroup$

I have a Cox model looking at time to death, considering several different covariates. What I would like to do is estimate the survival probability at a given time (in this case, t=30 years) for each person in my data set (I am not trying to use the results of my model to find information about another, separate data set).

My understanding is that the estimated survival for a given person at time (t) from a Cox model is as follows:

S(t|x) = exp (-exp(Bx) * (Ho(t))

Where B refers to the coefficients, X are the covariate values, and Ho(t) is the cumulative hazard at time (t).

Now there is predict.coxph, however I'm not sure it really gives me what I want. For example

pred<-predict(survfit(coxph(....)),type='expected')
surv_prob<-exp(-pred)

This gives me the survival probabilities using the baseline hazard, but not the cumulative hazard at time 30. (The documentation says "Predictions of type "expected" incorporate the baseline hazard and are thus absolute instead of relative...")

Should my process be to find the cumulative hazard estimate at t=30 using basehaz()? Or is there another way by which I can find these estimates?

Thank you for any help you can provide. I am happy to provide clarifications if needed.

$\endgroup$
8
$\begingroup$

The Cox proportional hazards model can be described as follows:

$$h(t|X)=h_{0}(t)e^{\beta X}$$ where $h(t)$ is the hazard rate at time $t$, $h_{0}(t)$ is the baseline hazard rate at time $t$, $\beta$ is a vector of coefficients and $X$ is a vector of covariates.

As you will know, the Cox model is a semi-parametric model in that it is only partially defined parametrically. Essentially, the covariate part assumes a functional form whereas the baseline part has no parametric functional form (it's form is that of a step function).

Additionally, the survival curve of the Cox model is:

$$\begin{align} S(t|X)&=\text{exp}\bigg(-\int_{0}^{t}h_{0}(t)e^{\beta X}\,dt\bigg)\\ &=\text{exp}\big(-H_{0}(t)\big)^{\text{exp}(\beta X)}\\ &=S_{0}(t)^{\text{exp}(\beta X)} \end{align}$$

where $H_{0}(t)=\int_{0}^{t}h_{0}(t)\,dt$, $S_{0}(t)=\text{exp}\big(-H_{0}(t)\big)$, $S(t)$ is the survival function at time $t$, $S_{0}(t)$ is the baseline survival function at time $t$ and $H_{0}(t)$ is the baseline cumulative hazard function at time $t$.

The R function basehaz() provides the estimated cumulative hazard function, $H_{0}(t)$, defined above. For example, I can fit a Cox PH model with a single covariate sexn in R as follows:

f=formula(sv~factor(sexn))
cox.fit=coxph(f)

I can then extract (and plot) the underlying baseline cumulative hazard function as follows:

bh=basehaz(cox.fit)
plot(bh[,2],bh[,1],main="Cumulative hazard function",xlab="Time",ylab="H0(t)")

Cumulative hazard function plot

Now, because of the proportional nature of the Cox model, to obtain the survival curves of the two groups defined by their sexn value I can just raise the cumulative hazard function to the power of the estimated coefficient for sexn.

For example, for my variable $sexn=\{0,1\}$, the two survival curves would be:

$$S(t|X=0)=\text{exp}(-H_{0}(t))^{\text{exp}(\beta(0))}=\text{exp}(-H_{0}(t))$$

and

$$S(t|X=1)=\text{exp}(-H_{0}(t))^{\text{exp}(\beta(1))}=\text{exp}(-H_{0}(t))^{\text{exp}(\beta)}$$

If you want to see the relative survival, you can just plot the curves as follows:

plot(bh[,2],exp(-bh[,1])^(exp(cox.fit$coef)),xlim=c(40,85),ylim=c(0,1),
col="red",main="Survival curves for two groups",xlab="Time",ylab="S(t|X)")
par(new=TRUE)
plot(bh[,2],exp(-bh[,1]),xlim=c(40,85),ylim=c(0,1),
col="blue",main="Survival curves for two groups",xlab="Time",ylab="S(t|X)")
legend("topright",c("sexn=1","sexn=0"),lty=c(1,1),col=c(2,4))

Relative survival

Thus, you can see that the group with $sexn=1$ has a lower survival than the group with $sexn=0$. If you want to measure the relative survival of the two groups you can do so in many ways. You can say that for two individuals (differing in only $sexn$) that start at $\text{Time}=40$, the individual with $sexn=1$ has a lower probability of surviving to any time $t>40$ compared with the individual with $sexn=0$.

I believe what you are trying to achieve is to calculate the survival estimate:

$$S(t=30|X)$$

This can be achieved by fitting a Cox model to a given survival object and applying the estimated coefficients to each individual depending on their individual covariates. This will scale the baseline survival curve and give you the desired survival estimate for each of your individuals.

$\endgroup$
  • $\begingroup$ Thank you. So for example in the AML data that comes in the survival package, I can use fit <- coxph(Surv(time, status) ~ x, data = aml) then basehaz(fit) and, if I wanted to use t=30 days, the baseline hazard at this time is 0.75653467. The value of beta (coef) is 0.916. So if x=1, S(t=30|x)=exp(-exp(0.916)*0.75653467)=0.151. Am I understanding this correctly? $\endgroup$ – Scott Jackson Jul 2 '17 at 15:20
  • $\begingroup$ @ScottJackson Yes that looks correct. $\endgroup$ – StatsPlease Jul 3 '17 at 0:01
  • $\begingroup$ @ScottJackson The two points you are referring to are here. $\endgroup$ – StatsPlease Jul 3 '17 at 0:35
  • $\begingroup$ I'm thinking something is off. First, the KM estimate at t=30 is 0.442. That's quite a discrepancy. Also I tried another method using predict.coxph and I get an estimate that's much closer (0.461). This is kind of messy but code> fit <- coxph(Surv(time, status) ~ x, data = aml) > aml2<-aml > aml2$time<-30 > temp<-predict(fit,type='expected',newdata=aml2) > mean(exp(-temp)) [1] 0.4607259code $\endgroup$ – Scott Jackson Jul 3 '17 at 2:52
  • $\begingroup$ On the other hand, I didn't calculate S(t=30|x=0) initially which is exp(-exp(0)*0.75653467)=0.469. So the mean survival probability at t=30 years is 0.303 $\endgroup$ – Scott Jackson Jul 3 '17 at 3:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.