/*****************************************************************/
/*    %Samuelsen macro program for SAS system, release 9.1.3     */
/*****************************************************************/

%MACRO Samuelsen(
  data=,             /* Specify the data set to be used.  If omitted,
                             the last-created data set is used. */
  time=,             /* REQUIRED. The failure-time variable */
  censor=,          /* REQUIRED. The censoring variable */
  c_values=,        /* List of values that correspond to right censoring */ 
  match=,          /* The variable counting the number of matched controls
                              for the case */
  x=                 /* Exposure variable */
) ;

options nonotes ;  ods listing close ; ods noresults ;

data _null_ ;
  t_start=datetime() ;
  call symput('t_start',t_start) ;
run ;

* data processing for analysis ;
data surv_time ;   set &data. ;
  if &censor. ne . then 
  do ;
    select ;
      when(&censor. in (&c_values.)) do ;
        censor=2 ;
        &match.=0 ;
        end;
      otherwise censor=1;
  end ;  end ;
run ;

proc summary data=surv_time(keep=&time. censor
                            where=(&time. ne . & censor ne .)) ;
  output out=N_sample(keep=_FREQ_) ;
run ;

data _null_; 
  set n_sample;
  call symput("n", trim(left(_FREQ_))) ;
run ;

ods output Attributes=d_name(where=(Label1 in ("Data Set Name" "データセット名"))
                             keep=Label1 cValue1) ;
proc contents data=&data. ;
run ;
data _null_ ;  set d_name(keep=cValue1) ;
  call symput("data_name", trim(left(cValue1))) ;
run ;
options notes ;
%PUT ;
%PUT NOTE: データセット &data_name. 内の &N. オブザーベーションを対象に解析しました. ;
options nonotes ;

ODS OUTPUT CrossTabFreqs=surv_weight0(keep=censor &time. frequency 
                                      where=(censor^=. and frequency^=0 and &time.^=.) ) ;
proc freq data=surv_time(keep= &time. censor) ;
  tables &time.*censor /nopercent nocol norow ;
run;

proc sql ;
  create view v_surv_time as
  select * from surv_time
  order by &time., censor, &match. ;
quit ;

proc sql ;
  create view v_surv_weight as
  select A.*, B.Frequency
  from v_surv_time as A, surv_weight0 as B
  where A.&time.=B.&time. and A.censor=B.censor
  order by A.&time., A.censor ;
quit ;

data surv_weight ;  
  set v_surv_weight ;  
  by &time. censor &match. ;
  D=2-censor ;
  retain cum0 0 cum 0.0 product_p 1.0 ;
  if first.censor then cum0=1 ;  else cum0=0 ;
  cum=frequency*cum0+cum ;
  rs=&N.-cum+frequency ; * # of at risk ;
  if rs-frequency ne 0 then nonselect0=1-(&match.*D/(rs-frequency)) ;
    else if rs-frequency=0 then nonselect0=1 ; 
  nonselect=lag(nonselect0) ;
  if _n_=1 then nonselect=1 ;
  product_p=product_p*nonselect ;
  p=(D=0)*(1-product_p)+(D=1)*1 ;
  select;
    when (p^=0) w=1/p;
    when (p=0) w=0;
  end ;
  keep id &time. censor  w &match. ;
run ;

proc sql ;
  create view sort_surv_weight as
  select &time., censor,  w, &match., id
  from surv_weight
  order by id ;
quit ; 

proc sql ;
  create view sort_surv_time as
  select * from surv_time
  where &x.^=.
  order by id ;
quit ; 

data surv ;  merge sort_surv_time(in=a) sort_surv_weight ;
  by id ;
  if a;
run;

* analysis ;
options notes ;
%PUT NOTE: PHREGプロシジャを用いてSamuelsen推定量を求めました. ;
title3 "Samuelsen macro was used for this analysis" ;
title5 "-------------------------- Samuelsen estimator --------------------------" ;
ods exclude CensoredSummary ;
ods listing ;  ods results ;
ods output ParameterEstimates=result(keep= variable estimate stderr ProbChiSq
                                           HazardRatio HRLowerCL HRUpperCL);
proc phreg covs data=surv( keep= &time. censor &x. w) ;
  model &time.*censor(2)=&x./ RL ties=efron ;
  weight w /norm ;
run;


/******************** Clean up temporary data sets ******************/
%clean:
options nonotes ;
proc datasets nolist;
  delete SORT_SURV_TIME SORT_SURV_WEIGHT SURV SURV_WEIGHT
         SURV_WEIGHT0 V_SURV_TIME V_SURV_WEIGHT D_NAME
    /memtype=all ;
quit;
title3 ;

/******************** Counting total process time *******************/
data _null_ ;
  t_end=datetime() ;
  r_sec=round(t_end-&t_start, 0.001) ;
  call symput('r_sec', trim(left(r_sec))) ;
run ;
options notes ;
%PUT NOTE: SAMUELSEN MACRO used ;
%PUT NOTE: 総演算時間: &r_sec. 秒 ;
%PUT ;
%mend ;

%samuelsen;