SAS

SASでカプランマイヤープロットを作図する(5)RMSTの可視化

今年のSASユーザー総会でRMST(Restricted mean survival time) の計算方法が紹介されてからカプランマイヤープロットにRMSTを追加する方法が発表されていますが、
ちょっとこれは・・・と思ったので今年最後の記事として私もRMSTの可視化の記事を書きたいと思います。

RMSTとは

RMST(Restricted mean survival time)は要はカプランマイヤープロットの指定した時点(tau)までのAUC(曲線下面積)のことです。
算出方法の詳細はSASユーザー総会でとても良い発表があるのでこちらをみれば十分でしょう。

可視化方法

RMSTはAUCなのでbandplotで作図すればよいでしょう。下側下限値は0で一定となるため、時点=tau以下の場合は0を格納した変数を用意すればOKです。
時点=tauにおける生存率が必要なので、proc lifetestのtimelistオプションでtauを指定して生存率を取得します。

またリスク集合は0になった時点以降はデータが存在しないので、時点をオプションで指定してもリスク集合は欠損となってしまいます。
それを防ぐには作図用データとマスターデータをマージしてデータを生成し、欠損値を置換します。
これはわりと重要なのでカプランマイヤープロットを作成する人はおぼえておいた方が良いと思います。

統計量テーブルをSASで作成する方法は公式ではlayout griddedを使用する方法が紹介されていますが、それだと表示する値をすべてマクロ変数に格納しないとならず
表示するデータが多いと現実的ではありません。

annotationで作成する方法も個々の文字列の座標をユーザーが制御しなければならず非常に大変です。特に文字を左揃えにといった調整が大変で、実務でこれをやるのは
お勧めしません。

annotationをお勧めしない理由は以下の記事でも紹介しています。annotationを使わないといけないシーンは現在はかなり限られるはずです。

SASのクソな部分をぶった切る!(3)作図環境SASのクソな部分シリーズはなぜかアクセスされるようなので、第3回目を投稿します。 このブログは一応GTLに関する記事も投稿してい...

表形式で表示させたいのであればaxistableかscatterplotを使う方が良いと思います。

あまり知られていないかもしれませんが、layout overlayはネストすることができます。つまりlayout overlayの中にlayout overlayを配置できます。
レイアウトの配置位置や大きさはオプションで指定可能ですので、これを利用すればプロットエリアの任意の場所にテーブルを配置できます。

以下のコードはカプランマイヤープロットにリスク集合テーブルとRMST推定値テーブルを配置しています。
時点=tauの位置に参照線を追加しています。

*----環境設定---------; 
*tau; 
%let tau=1500; 
*----KM曲線とRMSTの計算----------; 
ods output ProductLimitEstimates=est RMST=rmst survivalplot=surv; 
proc lifetest data=sashelp.BMT 
     timelist=(&tau.) 
     RMST(BC tau=&tau.) 
     plots=survival(atrisk=0 to 3000 by 500); 

time T * Status(0); 
strata Group; 
run; 

*-----作図データの前処理---------------; 
*時点=tauにおける生存率を追加する; 
*mergeステートメントで変数を上書きしているので嫌だったら別の方法でもいいかも・・・; 

data km1; 
merge surv est(keep=stratum timelist survival rename=(timelist=time stratum=stratumnum)); 
by stratumnum time; 
run; 

*atrisk=0となった時点以降はatriskはすべて欠損となってしまうため、マスターデータをマージする; 
proc sort data=surv nodupkey out=_master(keep=stratumnum stratum); 
by stratumnum;run; 

data master; 
set _master; 
do time =0 to 3000 by 500; 
   tatrisk=time; 
   output; 
end; 
run; 

data km2; 
merge master km1 ; 
by stratumnum time; 

*number at riskの欠損値を0に置き換える; 
if atrisk=. then atrisk=0; 

*auc作成のためbandplot境界を作成する; 
*GTLで作成できるけど、他の前処理の工程が多いのでデータステップで一緒に実施; 
if time<=&tau. then do; 
   lower=0; 
   upper=survival; 
end; 
else do; 
   lower=.; 
   upper=.; 
end; 
run; 

*RMSTテーブル; 
data rmst2; 
length stat $100 ; 
set rmst; 
stat=strip(put(round(estimate, 0.1),5.1)) || " (" || strip(put(round(stderr, 0.1),5.1)) || ")"; 
*axistable用軸変数として連番を作成; 
no+1; 
keep no group stat; 
run; 

data plot; 
merge km2 rmst2; 
run; 

*層ごとに線の色を指定したい場合はattribute mapを作成すること; 
proc template ; 
define statgraph rmst; 
begingraph; 
   layout lattice/ rows=2 rowweights=(0.9 0.1); 
      layout overlay / xaxisopts=(linearopts=(tickvaluesequence=(start=0 end=3000 increment=500))); 
      *生存率曲線; 
         stepplot x=time y=survival / 
            group=stratum 
            lineattrs=(pattern=solid) 
            name="step"; 
       *censor; 
          vectorplot x=time y=eval(ifn(censored^=., censored+0.02, .)) 
             xorigin=time yorigin=censored / 
             arrowheads=false 
             group=stratum; 
       *AUC;
          bandplot x=time limitupper=lower limitlower=upper / 
             group=stratum 
             type=step 
             fillattrs=(transparency=0.8); 
       *tau; 
          referenceline x=&tau. / 
             lineattrs=(color=black pattern=shortdash thickness=1);

      *凡例; 
          discretelegend "step" / 
             location=outside 
             valign=bottom; 
       *統計量テーブル 画像サイズを変更することを考慮しパーセンテージで位置を指定する; 

       layout overlay / width=35pct height=20pct 
                        halign=right valign=top 
                        pad=1pct outerpad=(top=2pct right=2pct) 
                        walldisplay=none opaque=false 
                        border=true 
           yaxisopts=(display=none reverse=true) 
           xaxisopts=(display=none); 

           innermargin/align=left; 
              axistable y=NO value=group /
                 valueattrs=(size=10) 
                 display=(values) 
                 title="RMST (SE)" titleattrs=(weight=bold); 
           *padで文字列の位置を微調整する; 
              axistable y=NO value=stat / 
                 valueattrs=(size=10) 
                 display=(values) 
                 pad=(left=6pct); 
             endinnermargin; 
         endlayout; 
     endlayout; 
     layout overlay / 
        xaxisopts=(display=none) 
        walldisplay=none yaxisopts=(display=none); 
        *リスク集合テーブル; 
         axistable x=tatrisk value=atrisk / 
             class=stratum; 
      endlayout; 
   endlayout; 
endgraph; 
end; 
run; 

/* 画像出力 */ 
ods graphics /imagename="RMST" imagefmt=png noborder ; 
ods listing gpath="出力パス" ; 

proc sgrender data=plot template=rmst; 
run;
実行結果