決定木によるカプランマイヤー生存時間分析

決定木使ったシンプルな生存分析です。

分類や回帰でも使用されるrpart()関数ですが、目的変数がSurv()関数になっています。
Surv()のパラメータはイベント発生までの時間とイベントフラグです。

survival.R

library(survival)
library(rpart)
library(partykit)

head(stagec)
# pgtime pgstat age eet    g2 grade gleason     ploidy
# 1    6.1      0  64   2 10.26     2       4    diploid
# 2    9.4      0  62   1    NA     3       8  aneuploid
# 3    5.2      1  59   2  9.99     3       7    diploid
# 4    3.2      1  62   2  3.57     2       4    diploid
# 5    1.9      1  64   2 22.56     4       8 tetraploid
# 6    4.8      0  69   1  6.14     3       7    diploid

dat <- stagec
tfit <- rpart(Surv(pgtime,pgstat)~., data = dat)

# Decision Tree
p.fit <- as.party(tfit)
plot(p.fit)

各ノードの曲線の値が知りたい場合は、ノード番号でデータを抽出してカプランマイヤー推定の関数を実行すると詳細が出てきます

# 左から2番目のノードの詳細がみたい
> # survival rate in nodes
> d5 <- dat[tfit$where==5,]
> summary(survfit(Surv(d5$pgtime,d5$pgstat)~1))
Call: survfit(formula = Surv(d5$pgtime, d5$pgstat) ~ 1)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  2.9     19       1    0.947  0.0512        0.852            1
  4.9     17       1    0.892  0.0724        0.760            1
  7.6      8       1    0.780  0.1220        0.574            1
  8.0      5       1    0.624  0.1703        0.366            1

全ノードの曲線を同じグラフで見たい時

dat$node = as.factor(tfit$where)
plot(survfit(Surv(dat$pgtime, event = dat$pgstat) ~ dat$node), col = dat$node)
legend("topright", legend = levels(dat$node), col = as.numeric(levels(dat$node)), lty = 1)

参照

https://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf
https://stackoverflow.com/questions/30700627/using-a-survival-tree-from-the-rpart-package-in-r-to-predict-new-observations