SASでカプランマイヤープロットを作図する（３）Number at riskテーブルの挿入

前回に引き続きカプランマイヤープロットの作図です。今回はNumber at riskのテーブルの挿入です。
Number at riskのテーブルデータは手動でやるとめんどくさいですが、Lifetestプロシジャから出力されるテーブルをそのまま使うと簡単に作図できます。

テーブルをプロットエリア内に設置する

LIFETESTプロシジャで特に指定しない場合の体裁と同じように、Number at riskのテーブルをX軸の上に配置します。

データセットは前回同様odsでsurvival plotを参照して作図データを出力します。plots=survival(atrisk)を指定するとデータセットにatriskとそれに対応する時間が格納された変数tatriskが新規に作成されます。この時maxlenオプションで系列名の文字列以上の数値を指定します。また以下のように時間のリストを追加すると指定した時間のnumber at riskのみが出力されます。

plots=survival(atrisk(maxlen=24)=0 to 2750 by 250)

これを手動でやるのはめんどくさいのでここは素直にプロシジャに頼ったほうが良いでしょう。

テーブルはinnermerginブロック内にaxistableを設置することで作成します。この時X軸の変数は時間timeではなくnumber at risk用に作成されたtatriskを指定するのがポイントです。また複数の系列がある場合は、classオプションにクラス変数を指定します。

ods output survivalplot=plotdata; 
proc LIFETEST data=Bmt2 plots=survival(atrisk(maxlen=24)=0 to 2750 by 250) notable; 
time T*Status(0); 
strata Group ; 
run; 

proc template; 
define statgraph surv_risk_inner; 
begingraph; 
discreteattrmap name="group" / discretelegendentrypolicy=attrmap; 
   value "ALL" / lineattrs=GraphData1(pattern=solid); 
   value "AML-Low Risk" / lineattrs=GraphData2(pattern=solid); 
   value "AML-High Risk" / lineattrs=GraphData3(pattern=solid); 
enddiscreteattrmap; 

discreteattrvar var=stratum attrmap="group" attrvar=_grp; 

   layout overlay / 
      xaxisopts=(label="Time" 
         linearopts=(viewmin=0 viewmax=2750 
            tickvaluesequence=(start=0 end=2750 increment=250)
         ) 
      ) 
      yaxisopts=(label="survival probability" 
         offsetmin=0.05 offsetmax=0.05 
         linearopts=(viewmin=0 viewmax=1 
            tickvaluesequence=(start=0 end=1 increment=0.25)
         )
      );

 /*生存曲線*/ 
   stepplot x=time y=survival / 
      name="surv" 
      group=_grp; 

/*打ち切り*/ 
   vectorplot x=time y=eval(ifn(censored^=., censored+0.02, .)) 
        xorigin=time yorigin=censored / 
      group=_grp 
      arrowheads=false; 

/*凡例*/ 
    entry "| : Censored " / 
      autoalign=(topright topleft) 
      textattrs=(size=10) 
      border=true 
      pad=5; 

   discretelegend "surv" / 
       location=outside ; 

/*リスク集合*/ 
innermargin / align=bottom; 
   axistable x=tatrisk value=atrisk / 
      class=_grp valueattrs=(color=black); 
endinnermargin; 
endlayout; 
endgraph; 
end; 
run; 

proc sgrender data=plotdata template=surv_risk_inner; 
run;

テーブルをプロットエリア外に設置する

Number at riskテーブルをプロットエリア外に設置するにはlatticeレイアウトを利用してグラフエリアを分割します。今回はグラフエリアを上下に分割し、下の区画にNumber at riskテーブルを設置します。

proc template; 
define statgraph surv_risk_outer; 
begingraph; 
discreteattrmap name="group" / discretelegendentrypolicy=attrmap; 
   value "ALL" / lineattrs=GraphData1(pattern=solid); 
   value "AML-Low Risk" / lineattrs=GraphData2(pattern=solid); 
   value "AML-High Risk" / lineattrs=GraphData3(pattern=solid); 
enddiscreteattrmap; 
discreteattrvar var=stratum attrmap="group" attrvar=_grp; 

   layout lattice / rows=2 rowweights=(0.85 0.15); 

/* グラフエリア */ 
   layout overlay / 
      xaxisopts=(label="Time" 
         linearopts=(viewmin=0 viewmax=2750 
            tickvaluesequence=(start=0 end=2750 increment=250)
         ) 
      ) 
      yaxisopts=(label="survival probability" 
         linearopts=(viewmin=0 viewmax=1 
            tickvaluesequence=(start=0 end=1 increment=0.25)
         )
      ); 

      stepplot x=time y=survival / 
         name="surv" 
         group=_grp; 

      vectorplot x=time y=eval(ifn(censored^=., censored+0.02, .)) 
           xorigin=time yorigin=censored / 
         group=_grp 
         arrowheads=false; 

      entry "| : Censored " / 
         autoalign=(topright topleft) 
         textattrs=(size=10) 
         border=true pad=5; 
 
      discretelegend "surv" / 
         location=outside ; 
   endlayout; 

/*Number at risk */ 
   layout overlay / walldisplay=none pad=(top=10px) 
      xaxisopts=(display=none) 
      yaxisopts=(display=none) ; 

      axistable x=tatrisk value=atrisk / 
         class=_grp 
         valueattrs=(color=black) 
         title="Number at risk" 
         titleattrs=(weight=bold); 
   endlayout; 
endlayout; 
endgraph; 
end; 
run; 

proc sgrender data=plotdata template=surv_risk_outer; 
run;

at riskテーブルの欠損値を置換する

上記の作図結果をよく見るとNumber at riskのテーブルの層＝ALL行で、値が表示されていない箇所があります。これは最初にat risk=0になった時点以降の at riskはすべて欠損値となるためです。そもそもこの時点のレコードが存在しないため、proc lifetestやグラフの設定で修正することはできません。本来は欠損値ではなく０となりますが、これを正しく表示させるには作図データセットを加工する必要があります。

具体的にはat riskを表示する時点と層の組み合わせを格納したデータセットと作図データセットをマージし、at riskが欠損値の場合に0で置換する方法を採用します。

/*層変数は連番にフォーマットを適用した文字列として出力*/ 
data master; 
length stratum $20; 
do stratumnum=1 to 3; 
do time=0 to 2750 by 250; 
   stratum=put(stratumnum,bmtfmt.); 
   output; 
end; 
end; 
keep stratum time; proc sort; by stratum time; 
run; 

proc sort data=plotdata; 
by stratum time;
run; 

data plotdata2; 
merge master plotdata; 
by stratum time; 
if atrisk=. then do; 
   atrisk=0; 
   tatrisk=time; 
end;

run;

加工した後のデータセットで作図すると以下のようになります。number at riskテーブルの値が表示されていなかった箇所に0が表示されています。