ggplot2でsurvival::survfitさせた生存曲線を描画する

こちらの survival 版: ggplot2でforecastインスタンスを描画する - StatsFragments

同じように survfit 用の fortify を定義すればよい。lung(肺ガンデータ)を使って、男女別のKaplan-Meier曲線を描いてみる。

library(survival)

library(ggplot2)
library(scales)

d.survfit <- survival::survfit(survival::Surv(time, status) ~ sex,
                               data = lung)

fortify.survfit <- function(survfit.data) {
  data.frame(time = survfit.data$time,
             n.risk = survfit.data$n.risk,
             n.event = survfit.data$n.event,
             n.censor = survfit.data$n.censor,
             surv = survfit.data$surv,
             std.err = survfit.data$std.err,
             upper = survfit.data$upper,
             lower = survfit.data$lower,
             strata = rep(names(survfit.data$strata), survfit.data$strata))
}

head(ggplot2::fortify(d.survfit))
#   time n.risk n.event n.censor      surv    std.err     upper     lower strata
# 1   11    138       3        0 0.9782609 0.01268978 1.0000000 0.9542301  sex=1
# 2   12    135       1        0 0.9710145 0.01470747 0.9994124 0.9434235  sex=1
# 3   13    134       2        0 0.9565217 0.01814885 0.9911586 0.9230952  sex=1
# 4   15    132       1        0 0.9492754 0.01967768 0.9866017 0.9133612  sex=1
# 5   26    131       1        0 0.9420290 0.02111708 0.9818365 0.9038355  sex=1
# 6   30    130       1        0 0.9347826 0.02248469 0.9768989 0.8944820  sex=1

y軸はパーセント表示したいので、scale_y_continuous で設定する。

ggplot(data = d.survfit) +
  geom_line(aes_string(x = 'time', y = 'surv', colour = 'strata')) +
  geom_ribbon(aes_string(x = 'time', ymin = 'lower', ymax = 'upper', fill = 'strata'), alpha = 0.5) + 
  scale_y_continuous(labels = scales::percent)

f:id:sinhrks:20141004193545p:plain

censorを描画したい場合は、

fortified <- ggplot2::fortify(d.survfit)
ggplot(data = fortified) +
  geom_line(aes_string(x = 'time', y = 'surv', colour = 'strata')) +
  geom_ribbon(aes_string(x = 'time', ymin = 'lower', ymax = 'upper', fill = 'strata'), alpha = 0.5) +
  geom_point(data = fortified[fortified$n.censor > 0, ],
             aes_string(x = 'time', y = 'surv'), shape = '+', size = 3) + 
  scale_y_continuous(labels = scales::percent)

f:id:sinhrks:20141004194611p:plain