/*
   PARAMETERS:
      data:   Name of anaylsis dataset.  Default is last dataset created.
      time:   Name of time variable. (Required)
      event:  Name of event variable. (Required)
              Values must be 0(censored) 1(uncensored).
      xvars:  Space-delimited list of independent variables. (Required)
      strata: Stratification variable.
      df:     Degrees of freedom for smoothing process.  Default is 4.
      pvars:  Variables to plot.  Default is all xvars.
      std:    Indicator to control plotting of standard error range.
              Values are yes(default) and no.
      points: Indicator to control plotting of data points.
              Values are yes(default) and no.
      rug:    Indicator to control plotting of rug of x values.
              Values are yes and no(default).
  */
%macro coxplot(data=_last_,
               time= ,
               event= ,
               xvars= ,
               strata= ,
               df=4,
               pvars=,
               std=yes,
               points=yes,
               rug=no);
   %let bad=0;
   %if &time= %then %do;
      %put ERROR: NO TIME VARIABLE GIVEN.;
      %let bad=1;
   %end;
   %if &event= %then %do;
      %put ERROR: NO EVENT VARIABLE GIVEN.;
      %let bad=1;
   %end;
   %if &xvars= %then %do;
      %put ERROR: NO XVARS GIVEN.;
      %let bad=1;
   %end;
   %if &df<3 | &df>7 %then %do;
      %put ERROR: DF MUST BE BETWEEN 3 AND 7.;
      %let bad=1;
   %end;
   %if %upcase(&std)^=YES & %upcase(&std)^=NO %then %do;
      %put ERROR: STD MUST BE YES OR NO.;
      %let bad=1;
   %end;
   %if %upcase(&points)^=YES & %upcase(&points)^=NO %then %do;
      %put ERROR: POINTS MUST BE YES OR NO.;
      %let bad=1;
   %end;
   %if %upcase(&rug)^=YES & %upcase(&rug)^=NO %then %do;
      %put ERROR: RUG MUST BE YES OR NO.;
      %let bad=1;
   %end;
   %let nxvar=0;
   %do %until(%scan(&xvars,&nxvar+1,' ')= );
      %let nxvar=%eval(&nxvar+1);
   %end;
   %if &pvars^= %then %let plotvars=&pvars;
   %else %let plotvars=&xvars;
   %let nplotvar=0;
   %do %until(%scan(&plotvars,&nplotvar+1,' ')= );
      %let nplotvar=%eval(&nplotvar+1);
   %end;
   %do i=1 %to &nplotvar;
      %let pv=%scan(&plotvars,&i,' ');
      %let found=0;
      %do j=1 %to &nxvar;
         %if &pv=%scan(&xvars,&j,' ') %then %let found=1;
      %end;
      %if &found=0 %then %do;
         %put ERROR: PLOT VARIABLE &PV NOT ON XVARS LIST.;
         %let bad=1;
      %end;
   %end;
   %if &bad=0 %then %do;
      proc phreg data=&data outest=__fit;
       model &time * &event(0) = &xvars/ties=breslow;
       output out=__temp1 resmart=rr/order=data;
       %if &strata^= %then %do;
          strata &strata;
       %end;
      data _null_;
       set __fit;
         %do i=1 %to &nplotvar;
            %let v=%scan(&plotvars,&i);
            call symput("beta&i",&v);
         %end;
      data __temp2;
       merge __temp1 &data;
       keep &event rr expected &plotvars;
         expected=&event-rr;
      %do i=1 %to &nplotvar;
         %let x=%scan(&plotvars,&i);
         %if %length(&x)=8 %then %let xx=%substr(&x,1,7);
         %else %let xx=&x;
         data __temp3;
          set __fit;
            keep beta;
            beta=&x;
         data __temp4;
          set __temp2;
            if _n_=1 then set __temp3;
            xbeta=&x*&&beta&i;
            off=log(expected)-xbeta;
         %daspline(&x,nk=&df)
         data __temp5;
          set __temp4;
            &&_&xx
         %global _print_;
         %let _print_=off;
         proc genmod data=__temp5;
         %let ndummies=%eval(&df-2);
          model &event=&x &xx.1-&xx.&ndummies/offset=off dist=poisson
                                            obstats;
          make 'OBSTATS' out=__temp6;
         data __temp7;
          merge __temp5 __temp6;
            lower=xbeta-2*std;
            upper=xbeta+2*std;
            res=xbeta+resdev;
         proc sort; by &x;
         %if %upcase(&rug)=YES %then %do;
            data __anno;
             set __temp7;
            keep x xsys y ysys function;
               x=&x;
               xsys='2';
               ysys='1';
               function='move';
               y=0;
               output;
               function='draw';
               y=5;
               output;
         %end;
         proc gplot data=__temp7
         %if %upcase(&rug)=YES %then %do;
            annotate=__anno
         %end;
         ;
         %if %upcase(&std)=YES %then %do;
            plot xbeta*&x=1
                 lower*&x=2
                 upper*&x=3
                 %if %upcase(&points)=YES %then %do;
                    res*&x=4
                 %end;
             /overlay vaxis=axis1 haxis=axis2;
         %end;
         %if %upcase(&std)=NO %then %do;
            plot xbeta*&x=1
                 %if %upcase(&points)=YES %then %do;
                    res*&x=4
                 %end;
             /overlay vaxis=axis1 haxis=axis2;
         %end;
          symbol1 i=join l=1;
          symbol2 i=join l=2;
          symbol3 i=join l=2;
          symbol4 i=none v=dot    h=.1 cm;
          axis1 label=(r=0 a=90 "smooth(&x)  df=&df");
          axis2 label=("&x");
          run; quit;
      %end;
   %end;
%mend coxplot;