1. 先設定好變數兩組，當然，也可以直接帶入 database 而略過這步。
>survival=rnorm(20,100,60)
>status=c(1,0,1,1,1,0,0,0,0,0,1,1,1,0,0,1,1,0,0,1)
把資料合在一起，看看長度是否正確。
>data=cbind(survival, status)
2.記得安裝package："survival"
3.
先來試試這個 formula（公式正不正確），此步可省略。
> Surv(survival,status)
重來了，要開始用 survfit 這個 function
> survfit(Surv(survival,status)~1)
Call: survfit(formula = Surv(survival, status) ~ 1)
records n.max n.start events median 0.95LCL 0.95UCL
20.0 20.0 20.0 10.0 99.6 90.4 NA
按Enter後會出現結果。
survfit這個function裡，有幾個要素。
formula後方要加上 "~1" ，表示分組的意思。
後面還有一項 "data=...." 這個下次找機會搞懂。
4.定義 myfit 為 survfit(.....)
再利用 summary 來秀 survival table
> myfit=survfit(Surv(survival,status)~1)
> summary(myfit)
Call: survfit(formula = Surv(survival, status) ~ 1)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
46.5 20 1 0.950 0.0487 0.859 1.000
61.5 18 1 0.897 0.0689 0.772 1.000
72.6 16 1 0.841 0.0844 0.691 1.000
79.9 14 1 0.781 0.0974 0.612 0.997
83.3 13 1 0.721 0.1069 0.539 0.964
90.4 10 1 0.649 0.1180 0.454 0.927
91.1 9 1 0.577 0.1250 0.377 0.882
91.7 8 1 0.505 0.1285 0.306 0.831
99.6 7 1 0.433 0.1288 0.241 0.775
111.6 4 1 0.324 0.1345 0.144 0.731
5.
不小心按錯的東西，但似乎會跑出來存活時間的一般統計資料，如：minimal, maximal, median, 1st Q, 3rd Q, mean....
> summary(myfit$time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
46.49 76.19 90.04 88.38 103.50 124.10
6.真正要鍵入的是下面這行，會出現event的survival time
> summary(myfit)$time
[1] 46.48571 61.54051 72.60968 79.86437 83.26426 90.41380 91.09167 91.65930 99.62062 111.61128
7.重點來囉！
又到 plot 出場的時候了。
按照plot所需，依序輸入：儲存survfit的變數，圖表抬頭，X軸抬頭，Y軸抬頭。
> plot(myfit,main="KM curve", xlab="time", ylab="proportion of survival")
就會在另一個視窗跳出 KM survival curve 了。
Reference:
http://stat.ethz.ch/education/semesters/ss2011/seminar/diez.pdf
先附上網路上找到的東西，改天再繼續研究。
張貼留言