7
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ The name is R not r! $\endgroup$ Commented Jun 26, 2021 at 19:35

2 Answers 2

7
$\begingroup$

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:

  1. 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.

  2. 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.

$\endgroup$
1
  • 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 that cc and ss 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
7
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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