打ち切りをヒゲにした Kaplan-Meier plot

LIFETEST Procedure のグラフ機能とか、GPLOT Procedure の symbol statement の value option の指定だけでは、打ち切りをヒゲにして出力できない。
ということで annotate を使って実装してみた。

proc datasets lib = work kill; run;
option linesize = 130 pagesize = 9999 mprint;
dm 'log; clear; output; clear';

%let execpath = " ";
%let Path = " ";
%macro setexecpath;
  %let execpath = %sysfunc(getoption(sysin));
  %if %length(&execpath) = 0 %then
    %let execpath = %sysget(sas_execfilepath);
  data _null_;
    do i = length("&execpath") to 1 by -1;
      if substr("&execpath", i, 1) = "\" then do;
        call symput("Path", substr("&execpath", 1, i));
        stop;
      end;
    end;
  run;
%mend setexecpath;
%setexecpath;
libname Out "&Path";

/*****************************************************************/
/* pseudo random number generation */
data data;
  call streaminit(945678);
  do stra = 1, 2;
    do j = 1 to 100;
      if stra = 1 then lambda = 0.2;
      if stra = 2 then lambda = 0.3;
      time = rand('Exponential') / lambda;
      /* random censoring */
      if (time > 10) & rand('Uniform') > 0.3 then censor = 1;
      else censor = 0;
      output;
  end; end;
  keep stra time censor lambda;
run;



/*****************************************************************/
/* Graph setting */
%let data = data;
%let time = time;
%let censor = censor;
%let cens_val = 1;             
%let stra = stra;
%let stra_val1 = 1;
%let stra_val2 = 2;

/*****************************************************************/
ods exclude all;
proc lifetest data = &data outsurv = outs;
  time &time * &censor(&cens_val);
  strata &stra;
run;
ods select all;

proc sort data = outs; by &stra; run;
data outs2;
  set outs; retain Svt 0;
  if Survival = . then Survival = Svt;
  else Svt = Survival;
  drop STRATUM SDF_LCL SDF_UCL;
run;
data outs2;
  set outs2(where=(_CENSOR_ IN(&cens_val.)))
      outs2(where=(_CENSOR_ ^= &cens_val | _CENSOR_ = .) drop=Svt);
run;
proc sort data = outs2; by &time; run;

data Graph;
  merge
    outs2(where=(&stra=&stra_val1) rename=(Survival=Sv1) drop=Svt)
    outs2(where=(&stra=&stra_val2) rename=(Survival=Sv2) drop=Svt)
    outs2(where=(&stra=&stra_val1&Svt1^=.) rename=(Svt=Svt1))
    outs2(where=(&stra=&stra_val2&Svt2^=.) rename=(Svt=Svt2));
  by &time; drop &stra;
  if Svt1 ^= . then Svtx = Svt1;
  else if Svt2 ^= . then Svtx = Svt2;
  else Svtx = .;
run;
proc sort data = Graph; by &time; run;


/*****************************************************************/
/* Print */
goptions reset = all;
goptions ftext = "Times New Roman" ftitle = "Times New Roman";
goptions vsize = 6 in hsize = 6 in htitle = 1.4 htext = 1.4;
options linesize = 98 pagesize = 200;

filename grafout "&Path.KMPlot2.emf";
goptions device = emf gsfname = grafout gsfmode = replace;

%let color1 = cx5ECD22;
%let color2 = cxFAA55C;

data anno;
  length color function $8;
  retain xsys ysys '2' hsys '4' when 'a' color 'black' size 2;
  set Graph;
  if Svtx ^= . then do;
    if Sv1 ^= . then color = "&color1";
    if Sv2 ^= . then color = "&color2";
    function='move'; xsys='2'; ysys='2'; x=&time; y=Svtx; output;
    function='draw'; xsys='2'; ysys='2'; x=&time; y=Svtx+0.015;  output;
  end;
run;

proc gplot data = Graph;
  plot (Sv1 Sv2) * &time / vaxis = axis1 haxis = axis2
    anno = anno overlay skipmiss noframe
    autovref cautovref = cxE9DECA;
  axis1 label = (a=90 'Propotion Surviving') minor = none;
  axis2 label = ('Time') order = (0 to 22 by 2) minor = none;
  symbol1 i = steplj c="&color1" w = 2;
  symbol2 i = steplj c="&color2" w = 2;
run; quit;

勾配が急だとヒゲが全く見えないので、他の記号にした方が方がいい。

annotate が実は便利だということに気づいて

こんなものも作っていたり (No. of at risk と一部の統計量をマクロで表示する)。