I'm trying to understand how cox.zph
function in R programming language works and I find myself not knowing what km
transform mean. I get the rank
transform and obviously the identity
one, but km
is still not clear to me and codes for this function did not make it clearer.
-
$\begingroup$ The name is R not r! $\endgroup$– kjetil b halvorsen ♦Commented Jun 26, 2021 at 19:35
2 Answers
km stands for Kaplan-Meier estimator.
$$\hat{S}(t) = \prod_{i: t_i \le t}\left(1-\frac{d_i}{n_i} \right)$$
with $t_{i}$ a time when at least one event happened, $d_i$ the number of events (i.e., deaths) that happened at time $t_{i}$ and ${\displaystyle n_{i}}$ the individuals known to have survived (have not yet had an event or been censored) up to time $t_{i}$.
Here is quote from the paper Cox Proportional-Hazards Regression for Survival Data:
Tests and graphical diagnostics for proportional hazards may be based on the scaled Schoenfeld residuals; these can be obtained directly as residuals(model, "scaledsch"), where model is a coxph model object. The matrix returned by residuals has one column for each covariate in the model. More conveniently, the cox.zph function calculates tests of the proportional-hazards assumption for each covariate, by correlating the corresponding set of scaled Schoenfeld residuals with a suitable transformation of time [the default is based on the Kaplan-Meier estimate of the survival function, $K(t)$].
To know why the choice of km
as the default, Dr. Kevin E. Thorpe cited Dr. Therneau's reply in the R-news:
There are 2 reasons for making the KM the default:
Safety: The test for PH is essentially a least-squares fit of line to a plot of f(time) vs residual. If the plot contains an extreme oulier in x, then the test is basically worthless. This sometimes happens with transform= identity or transform =log. It doesn't with transform='KM'.
As a default value for naive users, I chose the safe course.
A secondary reason is efficiency. In DY Lin, JASA 1991 Dan-Yu argues that this is a "good" test statistic under various assumptions about censoring. (His measure has the same score statistics as the KM option).
But #1 is the big one.
Terry T.
-
1$\begingroup$ Thank you @SiongThyeGoh, but still I do not understand where are the transform times come from, here's (a naive) example
library(survival) library(asaur) library(landest) data("pharmacoSmoking") fit<-coxph(Surv(ttr, relapse) ~ ageGroup4, data=pharmacoSmoking) times<-cox.zph(fit, transform = "identity")$x cc<−cox.zph(fit)$x ss<-surv.km(pharmacoSmoking$ttr,pharmacoSmoking$relapse, times)
. I was guessing thatcc
andss
should have sth to do with each other but I don't see it. What I don't understand is the meaning of based in "default is based on the Kaplan-Meier estimate" $\endgroup$– Nuta.Commented Jun 2, 2019 at 10:59
Like the original poster, I also wondered what, exactly, is the transformation "based on the Kaplan-Meier estimate" doing? Tracking this down proved to be more difficult than you would expect, as pretty much every source I found used some variant of the "based on Kaplan-Meier" language without further explanation. I eventually resorted to figuring it out from the source code.
The transformation is computed in lines 77-91 of cox.zph.R. The relevant code is:
if (is.character(transform)) {
tname <- transform
ttimes <- switch(transform,
'identity'= times,
'rank' = rank(times),
'log' = log(times),
'km' = {
temp <- survfitKM(factor(rep(1L, nrow(y))),
y, se.fit=FALSE)
# A nuisance to do left continuous KM
indx <- findInterval(times, temp$time, left.open=TRUE)
1.0 - c(1, temp$surv)[indx+1]
},
stop("Unrecognized transform"))
}
What is happening in the 'km'
branch is that the code is performing a Kaplan-Meier estimate $\hat{S}(t)$ of the survival function. The last two lines in the branch are implementing the transform $t \rightarrow \hat{S}(t)$. In other words, since the survival function is monotonic, we can replace the time coordinate with the corresponding value of the (K-M estimate of the) survival function. This has the advantage of fixing the x-values between 1 and 0 (vs. 0 and infinity for the raw time coordinate), so that observations at very large $t$ values don't have undue influence on the fit (as described in the R-news post described in the other answer).
Finally, in Applied Survival Analysis Using R, p. 99, the author does a comparison of the different transforms and notes:
The rank transformation yields a similar p-value to what we found with the "km" transformation.
Looking at the code above, the penultimate line of the "km" branch is looking up the intervals of the results of the K-M estimate corresponding to each time in the input. Since the intervals in the K-M estimate will be defined by the times in the input, in the absence of ties this lookup is just the rank of the times. So, we can view the KM transform as remapping equally spaced values from 1 to N into equally spaced values from 1 to $S(t_\max)$. Either way the linear fit to the residuals is going to be the same, which explains why the results of the two tests are nearly identical.