Solved – Comparison of Kaplan-Meier curves across ordered groups

kaplan-meierordinal-datasurvival

I am familiar with the log-rank test for comparing multiple Kaplan-Meier curves, but I am looking for a test that will compare across ordered groups (an ordinal variable). A significant result from the log rank test indicates that at least one of the groups has a survival curve different than the others. I want a test where a significant test indicates that there is a monotonic trend in the survival curves over the levels of the ordinal grouping variable.

As a more concrete example, consider survival curves for cancer of different stages, and I want to show not just that the curves differ, but that the survival worsens as stage increases.

As a bonus, pointers to implementations of the algorithm in R would be appreciated.

Edit to address @dardisco comment:

When I tried to find a suitable test, I did find a reference so something that seemed appropriate in the documentation for MedCalc software at http://www.medcalc.org/manual/kaplan-meier.php. Two quotes that sounded like what I wanted from that are

Options: Linear trend for factor levels: Allows testing for a linear trend across levels of the factor. It is appropriate if factor levels have a natural ordering (for example, factor codes represent doses applied to different groups). Kaplan-Meier assumes that the factor levels are equally spaced.

and

Logrank test for trend: If more than two survival curves are compared, and there is a natural ordering of the groups, then MedCalc can also perform the logrank test for trend. This tests the probability that there is a trend in survival scores across the groups.

But no specific algorithm for this test is stated.

Best Answer

Function comp in survMisc may be close to what you're after. From the docs:

Tests for trend are designed to detect ordered differences in survival curves.
That is, for at least one group: $$S_1(t) \geq S_2(t) \geq ... \geq S_K(t) \quad t \leq \tau$$ where $\tau$ is the largest $t$ where all groups have at least one subject at risk. The null hypothesis is that $$S_1(t) = S_2(t) = ... = S_K(t) \quad t \leq \tau$$ Scores used to construct the test are typically $s = 1,2,...,K$, but may be given as a vector representing a numeric characteristic of the group.
They are calculated by finding: $$Z_j(t_i) = \sum_{t_i \leq \tau} W(t_i)[e_{ji} - n_{ji} \frac{e_i}{n_i}], \quad j=1,2,...,K$$ The test statistic is: $$Z = \frac{ \sum_{j=1}^K s_jZ_j(\tau)}{\sqrt{\sum_{j=1}^K \sum_{g=1}^K s_js_g \sigma_{jg}}}$$ where $\sigma$ is the the appropriate element in the variance-covariance matrix (see COV).
If ordering is present, the statistic $Z$ will be greater than the upper $\alpha$-th percentile of a standard normal distribution.

For example:

library(survMisc)
data(larynx, package="KMsurv")
s4 <- survfit(Surv(time, delta) ~ stage, data=larynx)
comp(s4)

This will give tests for trend with various weights (as are used with the standard log-rank test):

$tests$trendTests
                                     Z       p
Log-rank                      3.718959 0.00010
Gehan-Breslow (mod~ Wilcoxon) 4.224765 0.00001
Tarone-Ware                   4.058010 0.00002
Peto-Peto                     4.129343 0.00002
Mod~ Peto-Peto (Andersen)     4.136319 0.00002
Trend F-H with p=1, q=1       2.396992 0.00827

The package has been changed slightly, and instead of tests or trendtests, use tft, which stands for... tests for trend. Furthermore, when using the comp() function, type it like this: comp(ten(<whatever your survival-fit object is>)).

Related Question