Solved – How to plot adjusted Kaplan-Meier Curves

cox-modeldata visualizationkaplan-meiersurvival

I am trying to plot adjusted Kaplan-Meier curves. I know publications like to see something graphical. But using R, I don't know how to go about adjusting for something like age, gender, income when graphing a survival curve.

Otherwise my curves will always be just crude and unadjusted, which I'm guessing people will not like. Any ideas?

I have 4 groups to compare.

Best Answer

The only way to provide differential survival with true KM curves is to generate new curves for different groups. You could then display a curve for all persons of group 3, for example. The number of units in each group will decrease as the number of strata increase. However, this method is empiric and does not truly adjust the sample to some chosen set of values.

I am most familiar with methods for obtaining adjusted curves derived from Cox or parametric survival models. Generally speaking, the role of an adjusted curve is to graphically display the expected mortality (or mortality transformed to survival) of the sample if a single or combination of values is set to some fixed value or set of values, respectively.

For example, one might find from a Cox model the hazard ratio for blood pressure is 1.1. Thus, for each 1 unit increase in blood pressure, the average hazard at a given time point multiplies by 1.1. Now, if we wanted to display the mortality curve for all units under analysis (sample) adjusted to a blood pressure of 1 standard deviation above the mean, we could display an adjusted curve.

Here is a self-contained example using group for your reference. Note that the final, adjusted curve is for the mean of group which, for most applications, would have no real meaning. Also note that transformation from mortality to survival are required for this method.

library("survival")
require("survival")

days <- rpois(100, 3)
status <- rbinom(100,1,0.34)
group <- sample(c(1,2,3,4), 100, replace=TRUE)
df <- data.frame(days, status)

#overall survival
surv <- survfit(Surv(days, status)~1)
summary(surv)
plot(surv)

#survival by group
kmsurv <- survfit(Surv(days,status) ~ strata(group), df)
plot(kmsurv)

#survival adjusted to group effect
cox <- coxph(Surv(days,status)~group, df)
summary(cox)
plot(survfit(cox))

Essential information on the R code can be found here.

Lastly, it is my opinion that adjusted survival analysis is generally statistical sleight of hand as the adjustment process can 1) be used for unrealistic patterns of covariates, 2) fool the reader into believing that non-significant effects can result in some displayed survival/mortality pattern and 3) be confused with empiric curves, leading readers to believe you have more events or information for each subgroup/pattern than you actually possess. I would carefully consider why adjusted curves are desirable over adjusted hazard ratios before spending too much time.

Related Solutions

Solved – Comparison of Kaplan-Meier curves across ordered groups

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

Solved – Estimating median survival times from Kaplan-Meier plot inspection

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

Related Solutions

Solved – Comparison of Kaplan-Meier curves across ordered groups

Solved – Estimating median survival times from Kaplan-Meier plot inspection

Related Question