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>)).
The time at which the Kaplan-Meier survival curve crosses
the 50% line is the non-parametric estimate of the median
survival time. If the Kaplan-Meier curve does not cross the 50% line, then the non-parametric estimate is not defined. At this stage, I can see two simple options:
- use another quantile (e.g. 0.75) to compare the two groups;
- approximate the survival curve by means of a parametric fit and derive the median survival time using the model.
As a side note, the typical underlying assumption in survival analysis is that any subject will eventually experience the event of interest if he is observed for a sufficient long period of time. In other words, the survival curve will reach zero (and in particular cross the 50% line) if the observation period lasts a sufficient long time. In contrast, cure modelling deals with situations where the survival curve eventually reaches a plateau. In that case, it is not correct to use a parametric model to extrapolate the survival curve below the plateau.
Best Answer
One thing to keep in mind with the Kaplan-Meier survival curve is that it is basically descriptive and not inferential. It is just a function of the data, with an incredibly flexible model that lies behind it. This is a strength because this means there is virtually no assumptions that might be broken, but a weakness because it is hard to generalise it, and that it fits "noise" as well as "signal". If you want to make an inference, then you basically have to introduce something that is unknown that you wish to know.
Now one way to compare the median survival times is to make the following assumptions:
Now the "most conservative" way to use these assumptions is the principle of maximum entropy, so you get:
$$p(T_{i}|t_{i})= K exp(-\lambda T_{i})$$
Where $K$ and $\lambda$ are chosen such that the PDF is normalised, and the expected value is $t_{i}$. Now we have:
$$1=\int_{0}^{\infty}p(T_{i}|t_{i})dT_{i} =K \int_{0}^{\infty}exp(-\lambda T_{i})dT_{i} $$ $$=K \left[-\frac{exp(-\lambda T_{i})}{\lambda}\right]_{T_{i}=0}^{T_{i}=\infty}=\frac{K}{\lambda}\implies K=\lambda $$ and now we have $E(T_{i})=\frac{1}{\lambda}\implies \lambda=t_{i}^{-1}$
And so you have a set of probability distributions for each state.
$$p(T_{i}|t_{i})= \frac{1}{t_{i}} exp\left(-\frac{T_{i}}{t_{i}}\right)\;\;\;\;\;(i=1,\dots,N)$$
Which give a joint probability distribution of:
$$p(T_{1},T_{2},\dots,T_{N}|t_{1},t_{2},\dots,t_{N})= \prod_{i=1}^{N}\frac{1}{t_{i}} exp\left(-\frac{T_{i}}{t_{i}}\right)$$
Now it sounds like you want to test the hypothesis $H_{0}:T_{1}=T_{2}=\dots=T_{N}=\overline{t}$, where $\overline{t}=\frac{1}{N}\sum_{i=1}^{N}t_{i}$ is the mean median survivial time. The severe alternative hypothesis to test against is the "every state is a unique and beautiful snowflake" hypothesis $H_{A}:T_{1}=t_{1},\dots,T_{N}=t_{N}$ because this is the most likely alternative, and thus represents the information lost in moving to the simpler hypothesis (a "minimax" test). The measure of the evidence against the simpler hypothesis is given by the odds ratio:
$$O(H_{A}|H_{0})=\frac{p(T_{1}=t_{1},T_{2}=t_{2},\dots,T_{N}=t_{N}|t_{1},t_{2},\dots,t_{N})}{ p(T_{1}=\overline{t},T_{2}=\overline{t},\dots,T_{N}=\overline{t}|t_{1},t_{2},\dots,t_{N})}$$ $$=\frac{ \left[\prod_{i=1}^{N}\frac{1}{t_{i}}\right] exp\left(-\sum_{i=1}^{N}\frac{t_{i}}{t_{i}}\right) }{ \left[\prod_{i=1}^{N}\frac{1}{t_{i}}\right] exp\left(-\sum_{i=1}^{N}\frac{\overline{t}}{t_{i}}\right) } =exp\left(N\left[\frac{\overline{t}}{t_{harm}}-1\right]\right)$$
Where
$$t_{harm}=\left[\frac{1}{N}\sum_{i=1}^{N}t_{i}^{-1}\right]^{-1}\leq \overline{t}$$
is the harmonic mean. Note that the odds will always favour the perfect fit, but not by much if the median survival times are reasonably close. Further, this gives you a direct way to state the evidence of this particular hypothesis test:
assumptions 1-3 give maximum odds of $O(H_{A}|H_{0}):1$ against equal median survival times across all states
Combine this with a decision rule, loss function, utility function, etc. which says how advantageous it is to accept the simpler hypothesis, and you've got your conclusion!
There is no limit to the amount of hypothesis you can test for, and give similar odds for. Just change $H_{0}$ to specify a different set of possible "true values". You could do "significance testing" by choosing the hypothesis as:
$$H_{S,i}:T_{i}=t_{i},T_{j}=T=\overline{t}_{(i)}=\frac{1}{N-1}\sum_{j\neq i}t_{j}$$
So this hypothesis is verbally "state $i$ has different median survival rate, but all other states are the same". And then re-do the odds ratio calculation I did above. Although you should be careful about what the alternative hypothesis is. For any one of these below is "reasonable" in the sense that they might be questions you are interested in answering (and they will generally have different answers)
Now one thing which has been over-looked here is correlations between states - this structure assumes that knowing the median survival rate in one state tells you nothing about the median survival rate in another state. While this may seem "bad" it is not to difficult to improve on, and the above calculations are good initial results which are easy to calculate.
Adding connections between states will change the probability models, and you will effectively see some "pooling" of the median survival times. One way to incorporate correlations into the analysis is to separate the true survival times into two components, a "common part" or "trend" and an "individual part":
$$T_{i}=T+U_{i}$$
And then constrain the individual part $U_{i}$ to have average zero over all units and unknown variance $\sigma$ to be integrated out using a prior describing what knowledge you have of the individual variability, prior to observing the data (or jeffreys prior if you know nothing, and half cauchy if jeffreys causes problems).