/*------------------------------------------------------------------*
   | MACRO NAME  : nesttest
   | SHORT DESC  : Likelihood ratio test for logistic or proportional
   |               hazards models
   *------------------------------------------------------------------*
   | CREATED BY  : Lennon, Ryan                  (06/26/2006 15:33)
   *------------------------------------------------------------------*
   | PURPOSE
   |
   | This program will produce a likelihood ratio test for the comparison
   | of two nested models.  The models may be logistic regression or Cox
   | proportional hazards models.
   *------------------------------------------------------------------*
   | MODIFIED BY : Lennon, Ryan                  (03/24/2009  8:46)
   |
   | 1. Changed input data set parameter from DS to DATA
   | 2. Allow user to name the output data set
   | 3. Macro only uses 2 footnotes, instead of 3
   *------------------------------------------------------------------*
   | OPERATING SYSTEM COMPATIBILITY
   |
   | UNIX SAS v8   :   YES
   | UNIX SAS v9   :   YES
   | MVS SAS v8    :
   | MVS SAS v9    :
   | PC SAS v8     :
   | PC SAS v9     :
   *------------------------------------------------------------------*
   | MACRO CALL
   |
   | %nesttest(
   |            DATA= ,
   |            event= ,
   |            time= ,
   |            e_val= ,
   |            x1= ,
   |            x2= ,
   |            fam= ,
   |            OUT=_test_,
   |            noprint=1,
   |            DS=
   |          );
   *------------------------------------------------------------------*
   | REQUIRED PARAMETERS
   |
   | Name      : DATA
   | Default   :
   | Type      : Dataset Name
   | Purpose   : Input data set name
   |
   | Name      : event
   | Default   :
   | Type      : Variable Name (Single)
   | Purpose   : The event indicator variable, must be numeric
   |
   | Name      : time
   | Default   :
   | Type      : Variable Name (Single)
   | Purpose   : The time to event variable (Only required if FAM=COX), must be numeric
   |
   | Name      : e_val
   | Default   :
   | Type      : Number (Single)
   | Purpose   : The value of EVENT which indicates an event has occurred.  Due to
   |             this, only binary logistic modelling is possible.
   |
   | Name      : x1
   | Default   :
   | Type      : Variable Name (List)
   | Purpose   : Explanatory variables to be found in both models. Must be numeric.
   |
   | Name      : x2
   | Default   :
   | Type      : Variable Name (List)
   | Purpose   : Explanatory variable found in only the large model. Must be numeric.
   |             So
   |             Model 1: Y = X1
   |             Model 2: Y = X1 X2
   |             The macro tests the significance of the partial effect (given X1)
   |             of the X2 variables.
   |
   | Name      : fam
   | Default   :
   | Type      : Text
   | Purpose   : BIN = Do logistic regression (binary)
   |             COX = Do Cox PH regression
   |
   *------------------------------------------------------------------*
   | OPTIONAL PARAMETERS
   |
   | Name      : OUT
   | Default   : _test_
   | Type      : Dataset Name
   | Purpose   : Name of output data set containing results
   |
   | Name      : noprint
   | Default   : 1
   | Type      : Number (Single)
   | Purpose   : 1=Do not print output from LOGISTIC/PHREG (default)
   |             0=Print LOGISTIC/PHREG output
   |
   | Name      : DS
   | Default   :
   | Type      : Dataset Name
   | Purpose   : The input data set. This parameter retained only for functionality
   |             with previous programs. If both DATA and DS are specified, the DS
   |             parameter is ignored.
   |
   *------------------------------------------------------------------*
   | RETURNED INFORMATION
   |
   | The returned output is a PROC PRINT with the Model1 and Model 2
   | parameter estimates and p-values given for all explanatory variables.
   | The last line has the -2 log likelihood values for both models, their
   | difference, the degrees of freedom, and resulting p-value.
   |
   *------------------------------------------------------------------*
   | ADDITIONAL NOTES
   |
   | 1. Data sets created: _tmp_ _est1 _est2 __est1 __est2 _temp_ _F &OUT.
   | 2. Data sets deleted: _tmp_ _est1 _est2 __est1 __est2 _temp_ _F.
   |
   |
   *------------------------------------------------------------------*
   | EXAMPLES
   |
   |
   | data a;
   |   do i=1 to 1000;
   |     var1=normal(0);
   |     var2=normal(0);
   |     var3=normal(0);
   |     P_e = (1/(1 + 1/(exp(-1 + var1 + 2*var2))));
   |     Y1 = ranbin(0,1,p_e);
   |     T = -log(ranuni(0))/exp(var2 + 2*var3);
   |     C = -log(ranuni(0));
   |     time=min(t,c);
   |     event=(t<=c);
   |     output;
   |   end;
   |  run;
   |
   | title1 'Nested Cox PH models';
   | %nesttest(data=a, event=event, time=time,  e_val=1,
   |           x1=var1 var3, x2=var2, fam=cox, noprint=0);
   |
   | title1 'Nested Logistic Regression Models';
   | %nesttest(data=a, event=Y1, e_val=1,
   |           x1=var1 var3, x2=var2, out=ntout, fam=bin, noprint=1);
   |
   |
   |
   *------------------------------------------------------------------*
   | REFERENCES
   |
   | Casella, G and RL Berger, Statistical Inference, Wadsworth Inc.,
   | Belmont 1990.
   |
   *------------------------------------------------------------------*
   | Copyright 2009 Mayo Clinic College of Medicine.
   |
   | This program is free software; you can redistribute it and/or
   | modify it under the terms of the GNU General Public License as
   | published by the Free Software Foundation; either version 2 of
   | the License, or (at your option) any later version.
   |
   | This program is distributed in the hope that it will be useful,
   | but WITHOUT ANY WARRANTY; without even the implied warranty of
   | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
   | General Public License for more details.
   *------------------------------------------------------------------*/
 
%MACRO NESTTEST(data=, event=, time=, x1=, x2=, e_val=, fam=,
                noprint=1, out=_test_, ds=);
 
 
%LOCAL ERRORFLG I XV1 J XV2 XV NN NE RM NOBS F FPLUS1 IP2;
 
%LET ERRORFLG=0;
 
%IF (&DATA= ) & (&DS= ) %THEN %DO;
    %PUT "ERROR - No input data set identified in DATA parameter.
          Macro will stop executing.";
    %LET ERRORFLG=1;
  %END;
%ELSE %IF (&DATA^= ) & (&DS^= ) &
         (%UPCASE(&DATA)^=%UPCASE(&DS)) %THEN %PUT
 WARNING: Different input data sets were identified by the DATA and DS
input parameters. Data set &DATA (DATA) will be used.;
 
%IF (&DATA= ) %THEN %LET DATA=&DS;
 
%IF %upcase(&FAM)=COX and &TIME=  %THEN %DO;
  %PUT  ERROR - Variable <time> not defined for Cox model;
  %LET ERRORFLG=1;
 %END;
%IF &EVENT=  %THEN %DO;
  %PUT  ERROR - Variable <event> not defined;
  %LET ERRORFLG=1;
 %END;
%IF &E_VAL=  %THEN %DO;
  %PUT  ERROR - Variable <event value> not defined;
  %LET ERRORFLG=1;
 %END;
%IF &X1=  %THEN %DO;
  %PUT  ERROR - Variable(s) <X1> not defined;
  %LET ERRORFLG=1;
 %END;
%IF &X2=  %THEN %DO;
  %PUT  ERROR - Variable(s) <X2> not defined;
  %LET ERRORFLG=1;
 %END;
%IF &FAM=  %THEN %DO;
  %PUT  ERROR - Variable <family> not defined;
  %LET ERRORFLG=1;
 %END;
 
%LET I=1;
%LET XV1=%SCAN(&X1, &I);
%DO %WHILE (&XV1^= );
 %LET J=1;
 %LET XV2=%SCAN(&X2, &J);
 %DO %WHILE (&XV2^= );
   %IF &XV1=&XV2 %THEN %DO;
      %PUT ERROR - Variable &XV1 occurs in both <X1> and <X2>.;
      %LET ERRORFLG=1;
     %END;
   %LET J=%EVAL(&J+1);
   %LET XV2=%SCAN(&X2, &J);
  %END;
 %LET I=%EVAL(&I+1);
 %LET XV1=%SCAN(&X1, &I);
%END;
 
%IF &ERRORFLG=1 %THEN %DO;
  data _null_;
  error 'ERROR - detected in the input parameters for NESTTEST';
  %go to exit;
 %END;
 
 
data _temp_;
  set &DATA end=last;
  retain n nevents rm 0;
  %IF %upcase(&FAM)=BIN %THEN %DO;
    __time=.;
    if &event=. then do; rm+1; delete;  end;
   %END;
  %IF %upcase(&FAM)=COX %THEN %DO;
    __time=&time;
    if nmiss(&event,&time)>0 then do; rm+1; delete; end;
   %END;
  %LET I=1; %LET XV=%SCAN(&X1 &X2, &I);
  %DO %WHILE (&XV^= );
    if nmiss(&XV)>0 then do; rm+1; delete; end;
    %LET I=%EVAL(&I+1);
    %LET XV=%SCAN(&X1 &X2, &I);
   %END;
  n+1;
  if &event=&e_val then do;
     nevents+1; __ind=1;
    end;
  else __ind=0;
  keep &x1 &x2 __ind  __time;
  if last=1 then do;
    call symput('NN', trimn(left((put(n,6.)))));
    call symput('NE', trimn(left((put(nevents,6.)))));
    call symput('RM', trimn(left((put(rm,6.)))));
   end;
 run;
%PUT NOTE: &RM observations were removed due to missing values;
%PUT NOTE: in either the response or explanatory variables.;
 
%IF %upcase(&FAM)=BIN %THEN %DO;
 proc logistic data=_temp_ descending outest=_est1 covout
   %IF &NOPRINT=1 %THEN %DO; noprint %END;  ;
  model __ind = &X1 &X2;
  run;
 
 proc logistic data=_temp_ descending outest=_est2 covout
   %IF &NOPRINT=1 %THEN %DO; noprint %END;  ;
  model __ind = &X1;
  run;
%END;
 
%IF %upcase(&FAM)=COX %THEN %DO;
 proc phreg data=_temp_ outest=_est1 covout
   %IF &NOPRINT=1 %THEN %DO; noprint %END;  ;
   model __time*__ind(0) = &x1 &x2;
  run;
 
 proc phreg data=_temp_ outest=_est2 covout
   %IF &NOPRINT=1 %THEN %DO; noprint %END;  ;
   model __time*__ind(0) = &x1;
  run;
%END;
 
 
%DO J=1 %TO 2;
 proc transpose data=_est&J out=_tmp_;
  run;
 
 data _null_; set _tmp_ end=last;
   if last=1 then call symput('NOBS', left(put(_N_,3.)));
  run;
 %DO I=1 %TO &NOBS; %LOCAL VAR&NOBS; %END;
 data _null_; set _tmp_;
   %DO I=1 %TO &NOBS;
   if _N_=&I then call symput("VAR&I", _name_);
   %END;
  run;
 
 data __est&J;
   attrib _label_ length=$40;
   set _tmp_;
   retain _order 0; _order+1;
   keep _name_ _est&J _se&J _pv&J
     %IF &J=1 %THEN %DO;  _label_ _order %END;  ;
 
   %IF  %upcase(&FAM)=COX %THEN %DO;
     _esti_ = __time;
    %END;
 
   %IF  %upcase(&FAM)=BIN %THEN %DO;
     _esti_ = __ind;
    %END;
 
   if _NAME_='_LNLIKE_' then do;
     _est&J=-2*_esti_; _se&J=.; _pv&J=.;
     %IF &J=1 %THEN %DO; _label_='-2*LogLikelihood'; %END;
    end;
   else do;
     %DO I=1 %TO %EVAL(&NOBS-1);
       if _name_="&&VAR&I" then _se&J=sqrt(&&VAR&I);
      %END;
     _est&J=_esti_;
     wald = (_esti_*_esti_)/(_se&J*_se&J);
     _pv&J=1-probchi(wald, 1);
    end;
  run;
 
%IF  &J=1 %THEN %DO;
 data __est&J; set __est&J;  ____one=1;
   if _N_=1 then set _temp_(keep=&X1 &X2) point=____one;
   if _name_^='_LNLIKE_' then do;
      %DO I=1 %TO %EVAL(&NOBS-1);
        if _name_="&&VAR&I" then do;
        if _label_='' then call label(&&VAR&I, _label_);
       end;
      %END;
     end;
   drop &X1 &X2;
  run;
 %END;
 
proc datasets nolist; delete _tmp_; run;
proc sort data=__est&J; by _name_; run;
%END;
 
data &OUT;
  merge __est1(in=in1) __est2(in=in2);
  by _name_;
  if in1;
  x2=.; x2p=.; df=.;
  label _name_='Variable' _label_='Label'
        _est1='Model 1 Estimate' _est2='Model 2 Estimate'
        _se1='Model 1 S.E.'  _se2='Model 2 S.E.'
        _pv1='Model 1 P-value' _pv2='Model 2 P-value'
        x2='X^2 Statistic' x2p='X^2 P-value' df='DF';
 run;
 
proc sort data=&OUT; by _order;
data &OUT; set &OUT;
  drop  df_ _order;
  retain df_ 0;
  if _est2=. then df_+1;
  if _name_='_LNLIKE_' then do;
    x2=_est2-_est1; df=df_; x2p=1-probchi(x2, df);
   end;
 run;
 
 
*** FOOTNOTE STUFF ***;
proc sql ;
  create table work._f as select * from dictionary.titles where type='F';
  reset noprint;
  select nobs into :F from dictionary.tables
    where libname="WORK" & memname="_F";
 quit;
 
  %LET FPLUS1 = %EVAL(&F + 1);
  %LET FNOTE1= ;** Initialize at least one footnote **;
  %DO I=1 %TO &F; %LOCAL FNOTE&I; %END;
  data _null_; set _f;
    %IF (&F>=1) %THEN %DO I=1 %TO &F;
    if number=&I then call symput("FNOTE&I", trim(left(text))); %END;
   run;
 
 
proc print data=&OUT l noobs;
  var  _name_ _est1 _pv1 _est2 _pv2 x2 df x2p _label_;
  format _est1 _est2 9.3 _pv1 _pv2 x2p 6.4
         x2 8.3 df 2.0;
  footnote1 "DATA=&DATA, EVENT=&event, TIME=&TIME, E_VAL=&E_VAL, FAM=&FAM";
  footnote2 "Number of Observations: &NN, Number of Events: &NE";
  %IF (&F>=1) %THEN %DO I=1 %TO &F;
    %LET IP2=%EVAL(&I+2);
    footnote&IP2 "&&FNOTE&I";
   %END;
 run; footnote1; footnote2; footnote3;
 
  %IF (&F>=1) %THEN %DO I=1 %TO &F;
    footnote&I "&&FNOTE&I";
   %END;
 
 
 
proc datasets nolist; delete _est1 _est2 __est1 __est2 _temp_ _F;
 run;  quit;
 
%EXIT:
 
%MEND NESTTEST;