使用ODS GRAPHICS和PROC LIFETEST繪製存活曲線圖 (Survival Curve)

太久沒寫有的沒有的SAS程式了，這次來寫寫關於Kaplan-Meier Estimate、Log-Rank Test及KM Plot的心得吧。
這篇主要的目的為兩個，如何使用PROC LIFETEST進行存活資料的的分析（跳過統計理論的部分）， 另一個為使用SAS 9.2新的ODS GRAPHICS，快速產生勘用的存活機率圖，並將圖檔自動放入RTF(Rich Text Format)檔中。

依照慣例，一定要有資料才有感覺。資料引用自SAS 9.2 說明檔中的Example 49.2，為骨髓移植病患的資料，依據將病患分為不同的三組，希望能描述不同組別病患其disease-free survival time（事件為死亡或疾病復發，以下我會以存活時間描述）。觀察並檢定三組存活機率(Kaplan-Meier Estimate)有無相同，常用的方法為Log-Rank Test，而依據不同研究型態Wilcoxon test也可能會是不錯的選擇。

data BMT; input Group T Status @@;label T='Disease Free Time'; datalines; 1 2081 0 1 1602 0 1 1496 0 1 1462 0 1 1433 0 1 1377 0 1 1330 0 1 996 0 1 226 0 1 1199 0 1 1111 0 1 530 0 1 1182 0 1 1167 0 1 418 1 1 383 1 1 276 1 1 104 1 1 609 1 1 172 1 1 487 1 1 662 1 1 194 1 1 230 1 1 526 1 1 122 1 1 129 1 1 74 1 1 122 1 1 86 1 1 466 1 1 192 1 1 109 1 1 55 1 1 1 1 1 107 1 1 110 1 1 332 1 2 2569 0 2 2506 0 2 2409 0 2 2218 0 2 1857 0 2 1829 0 2 1562 0 2 1470 0 2 1363 0 2 1030 0 2 860 0 2 1258 0 2 2246 0 2 1870 0 2 1799 0 2 1709 0 2 1674 0 2 1568 0 2 1527 0 2 1324 0 2 957 0 2 932 0 2 847 0 2 848 0 2 1850 0 2 1843 0 2 1535 0 2 1447 0 2 1384 0 2 414 1 2 2204 1 2 1063 1 2 481 1 2 105 1 2 641 1 2 390 1 2 288 1 2 421 1 2 79 1 2 748 1 2 486 1 2 48 1 2 272 1 2 1074 1 2 381 1 2 10 1 2 53 1 2 80 1 2 35 1 2 248 1 2 704 1 2 211 1 2 219 1 2 606 1 3 2640 0 3 2430 0 3 2252 0 3 2140 0 3 2133 0 3 1238 0 3 1631 0 3 2024 0 3 1345 0 3 1136 0 3 845 0 3 422 1 3 162 1 3 84 1 3 100 1 3 2 1 3 47 1 3 242 1 3 456 1 3 268 1 3 318 1 3 32 1 3 467 1 3 47 1 3 390 1 3 183 1 3 105 1 3 115 1 3 164 1 3 93 1 3 120 1 3 80 1 3 677 1 3 64 1 3 168 1 3 74 1 3 16 1 3 157 1 3 625 1 3 48 1 3 273 1 3 63 1 3 76 1 3 113 1 3 363 1 ; /*程式01*/ PROC LIFETEST DATA=BMT ; TIME T * Status(0); STRATA Group; RUN;

「程式 01」為PROC LIFETEST的基本語法，TIME指令表示如何定義存活時間，而TIME後面的指令以中文方式說明，可寫為【時間變相 ＊ 設限(設限值)】，在BMT資料中病患追蹤時間於變項 『T』；而病患的事件有無發生則於變項 『Status』，其中「0」表示病患至研究結束前都沒有死亡或復發，稱設限(censor)。STRATA用來定義病患的組別。
不過「程式 01」在樣本數大時，可能會發生output非常的冗長，因SAS會列出所有Product-Limit Estimates的結果。因此如在語法中多增加一個指令NOTABLE和PLOTS=S，會得到較精簡的結果以及存活曲線圖，於「程式 02」。

/*程式02*/ PROC LIFETEST DATA=BMT NOTABLE PLOTS=S; TIME T * Status(0); STRATA Group; RUN;

「程式 02」的結果，在「Summary of the Number of Censored and Uncensored Values」可先檢查每個組的總人數(Total)、事件人數(Failed)、設限人數(Censored)。



我的習慣會先看存活曲線圖，由圖A中可以「看到」三組病人的存活表現應該是有差異的。在「Test of Equality over Strata」的Log-Rank和Wilcoxon test的P值均顯示，在此份資料中三組病患的存活機率有顯著的差異(僅為三組不全等)。

 



在SAS 9.2對ODS Graphics有大幅度的增修，以下內容將會介紹我常用的指令。
在「程式02」上面只要增加ODS GRAPHICS ON指令(程式03)，便可繪製較精緻的統計圖(圖B)。

/*程式03*/ ODS GRAPHICS ON; PROC LIFETEST DATA=BMT NOTABLE PLOTS=SURVIVAL; TIME T * Status(0); STRATA Group;

RUN;  

在PLOTS=SURVIVAL後面的刮號內的指令可以將存圖增加些有趣且實用的修改，例如NOCENSOR指令會將原本圖形上設限個案的符號(十字, + )移除。TEST指令則會將Log-Rank Test的P值呈現於圖形中。而ATRISK=<數列>則是將在列出各subjects at risk的人數列出，例如在圖C的第0時間點，三組at risk人數分別為38、54、45人，而到第2500時間點，三組at risk人數分別為0、2、1人。

/*程式04*/ PROC LIFETEST DATA=BMT NOTABLE PLOTS=SURVIVAL(NOCENSOR TEST ATRISK=0 TO 2500 BY 500); TIME T * Status(0); STRATA Group; RUN;

    

ODS GRAPHICS可以快速做出一張一張美美的統計圖，但如果有10張統計圖需要輸出，人工一張一張輸出還滿煩的（至少我很懶……）。遇此需求，ODS RTF或ODS PDF均可給分析者很大的幫助。

在「程式05」主要以ODS RTF將圖片一次輸出至Word檔案中。ODS RTF後面的指令FILE、STYLE和SELECT用途分別為，FILE宣告Word檔輸出的位置和檔案名稱（副檔名用RTF或DOC都可以）；STYLE設定輸出文件的風格，STATISTICAL風格還滿漂亮的～；SELECT為選擇ODS的輸出物件(output objects)，這裡只選擇存活圖的部分。在最後別忘了ODS RTF CLOSE。

/*程式05*/ ODS GRAPHICS ON; ODS RTF FILE="e:\ KMPLOTS.rtf" STYLE=STATISTICAL; ODS RTF SELECT SurvivalPlot; PROC LIFETEST DATA=BMT NOTABLE PLOTS=SURVIVAL; TIME T * Status(0); STRATA Group; RUN; ODS RTF SELECT SurvivalPlot; PROC LIFETEST DATA=BMT NOTABLE PLOTS=SURVIVAL( NOCENSOR ATRISK=0 TO 2500 BY 500 TEST ); TIME T * Status(0); STRATA Group; RUN; ODS RTF CLOSE; ODS GRAPHICS OFF;

這應該是我第二篇分享SAS 9.2關於ODS GRAPHICS的心得，正式場合所需的統計圖，傳統的繪圖程式，如PROC GPLOT、PROC GCHART等，還是相當重要的。
延伸閱讀：
1. Survival Analysis with SAS (http://www.ats.ucla.edu/stat/SAS/seminars/sas_survival/default.htm)
2. Tests of proportionality in SAS, Stata, R and SPLUS (http://www.ats.ucla.edu/stat/SAS/faq/test_proportionality.htm)

---2018/05/27 更新---
如果要修改此篇所繪製出的圖，可以參考我另一篇文章「修改PROC LIFETEST以ODS所繪製存活曲線圖」。