2010年3月15日

使用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的結果。因此如在語法中多增加一個指令NOTABLEPLOTS=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)。

P01

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

  P00

P02

在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;   P03-1

在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;

     P04-1

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所繪製存活曲線圖」。

8 則留言:

  1. style除了STATISTICAL外還有什麼其他的選擇嗎 = =

    回覆刪除
  2. 你好..我見你的sas 很強呢...
    我對sas有點問題...希望你可指導一下..
    data a;
    input yr y;
    cards;
    2005 10
    2006 20
    2007 30
    2008 40
    run;

    data b;
    set a;
    run;

    proc print data=a;
    run;

    Obs r yr y
    1 . 2005 10
    2 10 2006 20
    3 20 2007 30
    4 30 2008 40

    我想問..如何在output中...
    可新增一行r是第2開數開始等於y的第1個數..
    如此類推...謝謝..

    回覆刪除
  3. 您的問題應該不難
    用lag就好了。
    r=lag(y);

    敢問您的問題是來自於課業嗎?

    回覆刪除
  4. 請問在ODS GRAPHICS ON之下, 是否可以更改survival curve 線條的顏色, 粗細以及種類? 謝謝

    回覆刪除
  5. 我知道的應該是不能修改,除非去修改版模。

    回覆刪除
  6. 您好~最近看到paper可以將median,95%跟censor summary擺在survival plot右上方的圖,不知道是如何用ods做到呢?找遍help都查不到阿~

    回覆刪除
  7. 請問怎麼畫出校正過變項的SURVIVAL CURVE?

    回覆刪除
    回覆
    1. 你可以參考
      http://support.sas.com/documentation/cdl/en/statug/68162/HTML/default/viewer.htm#statug_phreg_examples08.htm

      刪除