RMSTのグラフを描いてみた話

とても反抗的で手の付けられない後輩が，SASユーザー総会でRMSTの発表をしてました．
RMSTオプションは良いなぁと思うのですが，ods gprahicsででてくるKM plotは特にいつも通りのものしかでません

Pythonのlifelinesライブラリとかだと割と簡単で

from matplotlib import pyplot as plt

from lifelines.utils import restricted_mean_survival_time

from lifelines.plotting import rmst_plot

bmt=impsas(folder,"bmt")

ix=bmt['Group']=='ALL'

T, E = bmt['T'], bmt['Status']

time_limit = 1500

kmf_exp = KaplanMeierFitter().fit(T[ix], E[ix], label='exp')

rmst_plot(kmf_exp, t=time_limit)

それっぽいものがでて，視覚的にもわかりやすいのですが，これをSASでやるのは面倒です

すでにSASoneさんの方で作図例があるのですが
http://sasonediver.blog.fc2.com/blog-entry-253.html 

annotation使わずに，描けないかなぁと思いました

Number at Risk表つくってる場合，通常，境界時間τは，そんな中途半端な時点にならずにat Riskだしてるポイントと被ることが多いので
それなら，信頼区間バンドつけるときの下限を0にしちゃえばいいんじゃねってことで

%let tau=1500;

ods output Survivalplot=SurvivalPlotData;

ods output RMST=RMST(rename=(Stratum=StratumNum));

proc lifetest data=sashelp.bmt  rmst(tau=&tau) plots=survival(atrisk=0 to 2500 by 500 );

    time T*STATUS(0);

    strata group;

run; 

data SurvivalPlotData_1;

  set SurvivalPlotData;

  if 0 then set RMST;

  if _N_=1 then do;

    declare hash h1(dataset:"RMST");

    h1.definekey("StratumNum");

    h1.definedata("Estimate");

    h1.definedone();

  end;

  retain _Survival;

  if ^missing(Survival) then _Survival=Survival;

  markerchar=" ";

  if ^missing(Censored) then hige=Censored+0.03;

  if time<=&tau then limit_time=time;

  h1.find();

  if tAtRisk eq &tau then do;

    tau=&tau;

    rmst=cats("RMST=",round(Estimate,0.1));

    ry=_Survival -0.05;

  end;

run;

proc sgplot data=SurvivalPlotData_1 noautolegend;

styleattrs datacontrastcolors=(black black black) datacolors =(blue yellow red);

  step x=time y=survival /CURVELABEL group=stratum lineattrs=(thickness=3);

  scatter x=time y=censored /NOERRORCAPS yerrorupper=hige errorbarattrs=(pattern=1 thickness=3)  markerchar=markerchar GROUP=stratum;

  band  x=limit_time upper=_survival lower=0/ group=stratum type=step transparency=0.7;

  text x=tau y=ry text=rmst/textattrs=(size=15);

  xaxistable atrisk / x=tatrisk class=stratum location=outside colorgroup=stratum valueattrs=(size=10);

  refline tau/axis=x;

  yaxis label="Survival Probability" min=0 offsetmax=0.1;

  xaxis label="Day" min=0 offsetmin=0.01 ;

run;

おー，なんかそれっぽい