/*---------------------------------------------------------------------------------
  Macro mfp8

   Version 6.5 vom 21.08.2003

   written by:   Carolina Meier-Hirmer (cmh@fdm.uni-freiburg.de)
                 Carina Ortseifen

   idea and scientific advices: Willi Sauerbrei

-----------------------------------------------------------------------------------------------*/
%macro mfp8(
            dsname=,                          /* dataset                                       */
            yname=,                           /* dependent variable                            */
            xname=,                           /* positive, continuous covariate(s)             */
            xbin=,                            /* binary covariate(s)                           */
            xinclude=,                        /* covariate(s) forced in all models             */
            model=N,                          /* kind of model                                 */
            m=2,                              /* global maximum degree of the polynomial       */
            pw=-2 -1 -0.5 0 0.5 1 2 3,        /* global power values of the polynomial         */
            dp=,                              /* degree of the polynomial for single variables */
            power=,                           /* powers of the polynomial for single variables */
            censvar=,                         /* censoring variable for survival data          */
            censval=0,                        /* censoring value(s) for survival data          */
            strata=,                          /* stratifying variable(s) for survival data     */
            ties=efron,                       /* method for treating ties for survival data    */
            alpha=0.05,                       /* significance level                            */
            mselect=ra2,                      /* method for choosing best model                */
            maxcycle=10,                      /* maximum number of cycles                      */
            macpath=,                         /* path for SAS-macros                           */
            dsout=fpout,                      /* result dataset                                */
            showres=r                         /* results in output window                      */
            );

/*------------------------------------------------------------------------------
  ---
  --- Parameter checking
  ---
  ------------------------------------------------------------------------------*/

options nodate;

%let order=1;
%let pref=_;
%let xname=%upcase(&xname);
%let yname=%upcase(&yname);
%let xbin=%upcase(&xbin);
%let model=%upcase(&model);
%let mselect=%upcase(&mselect);
%let showres=%upcase(&showres);

* check whether the dataset is inside the working library -> the macro stops;
%let bib1=%qscan(&dsname, 1, ".");
%let bib2=%qscan(&dsname, 2, ".");
%if %length(%trim(&bib2)) = 0 or %upcase(&bib1)=WORK %then
    %goto ENDLIB;

* model: N, L oder S;
%if %index("NSL",&model)=0 %then
   %goto ENDM;

* showres: N, D oder R;
%if %index("NRD",&showres)=0 %then
   %let showres=R;

* mselect: SEQ oder RA2;
%if %index("SEQRA2",&mselect)=0 %then
   %let mselect=RA2;

/*------------------------------------------------------------------------------
  ---
  --- Macro call and parameters
  ---
  ------------------------------------------------------------------------------*/

%if %index("N",&showres)=0 %then %do;

   data _null_;
      file print;
      put ;
      put 'MFP8 called with following settings';
      put ;
      put "Dataset name       : &dsname";
      put "Dependent variable : &yname";
      put "Covariates         : &xname &xbin";
      put "Covar. always incl.: &xinclude";
      put "Model (N, S or L)  : &model";
      put "Maximal degree     : &m";
      put "Power values       : &pw";
      put "Alpha              : &alpha";
      put "Model selection    : &mselect";

      %if %upcase(&model)=S %then
         %do;
         put "Censoring variable : &censvar";
         put "Censoring value    : &censval";
         put "Strata variable    : &strata";
         put "Ties               : &ties";
      %end;

      put "Path to SAS macros : &MacPath";
      put "Number of cycles   : &maxcycle";
      put ;
   run;
%end;

/*------------------------------------------------------------------------------
  ---
  --- Preparations
  ---
  ------------------------------------------------------------------------------*/


%include "&MacPath./boxtid.sas";
%include "&MacPath./xtop.sas";
%include "&MacPath./xvars.sas";
%include "&MacPath./fpmodels.sas";
%include "&MacPath./datasave.sas";
%include "&MacPath./exlabb.sas";
%include "&MacPath./exinc.sas";
%include "&MacPath./labs.sas";
%include "&MacPath./brename.sas";
%include "&MacPath./funcfm.sas";

/* all files in the SAS working directory are deletetd */
proc datasets lib=work nolist kill mt=data;
run;

%let xnameo=&xname;
%let xname=&xname &xbin;

/* determination of the number of explanatory variables */
%if %length(&xname)=0 %then %goto endnox;

data _null_;
   %let i=1;
   %let x=%qscan(&xname,&i,%str( ));
   %do %while (%length(%trim(&x))>0);
      %let mname&i=&x;
      %let i=%eval(&i+1);
      %let x=%qscan(&xname,&i,%str( ));
   %end;
run;
%let xmax=%eval(&i-1);
%if &xmax=1 %then %let maxcycle=1;

data _null_;
   %let i=1;
   %let x=%qscan(&xbin,&i,%str( ));
   %do %while (%length(%trim(&x))>0);
      %let i=%eval(&i+1);
      %let dp=&dp &x.:0;
      %let x=%qscan(&xbin,&i,%str( ));
   %end;
run;
%let xmaxb=%eval(&i-1);


/* determination of the number of included variables */
%let xinclude=%upcase(&xinclude);

data _null_;
   %let i=1;
   %let x=%qscan(&xinclude,&i,%str( ));
   %do %while (%length(%trim(&x))>0);
      %let i=%eval(&i+1);
      %let x=%qscan(&xinclude,&i,%str( ));
   %end;
run;
%let xincmax=%eval(&i-1);

/* determination of the number of strata variables */
data _null_;
   %let i=1;
   %let x=%qscan(&strata,&i,%str( ));
   %do %while (%length(%trim(&x))>0);
      %let i=%eval(&i+1);
      %let x=%qscan(&strata,&i,%str( ));
   %end;
run;
%let smax=%eval(&i-1);

/* determination whether there is a censoring variable */
%let censmax=0;
%if %index("S",&model)=0 %then %let censvar=;
%if %length(%trim(&censvar))>0 %then %let censmax=1;

%if %upcase(&model)=S and &censmax=0 %then %goto ends2;

/* allvar contains the names of all used variables */
%let allvar=&xinclude &xname &yname &strata &censvar;

/* only the required variables are kept */
data work.dsname;
   set &dsname;
   keep &allvar;
run;

/* determination of the number of observations */
%let nfirst=0;
data work.check;
   set sashelp.vtable;
   if libname="WORK" & memname="DSNAME" & memtype="DATA" then do;
      call symput("nfirst", put(nobs, 6.));
      stop;
   end;
run;

/* if there is no observation the program stops */
%if &nfirst=0 %then %goto ENDNOOBS;

/* if there is a survival time smaller than 0 the program stops */
%let nsecond=0;
%if %upcase(&model)=S %then %do;
   data work.test;
      set work.dsname;
      if &yname >= 0 then delete;
   run;
   data work.check;
   set sashelp.vtable;
      if libname="WORK"& memname="TEST" & memtype="DATA" then do;
         call symput("nsecond", put(nobs, 6.));
         stop;
      end;
   run;
%end;
%if &nsecond ne 0 %then %goto ENDS;
%let nsecond=0;

/* if there is a logit response with values out of {0,1} the program stops */
%if %upcase(&model)=L %then %do;
   data work.test;
      set work.dsname;
      if &yname = 0 then delete;
      if &yname = 1 then delete;
   run;
   data work.check;
   set sashelp.vtable;
      if libname="WORK" & memname="TEST" & memtype="DATA" then do;
         call symput("nsecond", put(nobs, 6.));
         stop;
      end;
   run;
%end;
%if &nsecond ne 0 %then %goto ENDL;

/* examination whether all specified variables are found in the data set */
%let nsecond=0;
data work.check;
   set sashelp.vtable;
   if libname="WORK" & memname="DSNAME" & memtype="DATA" then do;
      call symput("nsecond", put(nvar, 5.));
      stop;
   end;
run;

%if &nsecond ne %eval(&xmax+&xincmax+&smax+&censmax+1) %then %goto ENDX;

/* missing values are deleted */
data work.dsname;
   set work.dsname;
   array c &allvar;
   do i=1 to %eval(&xmax+&xincmax+&smax+&censmax+1);
      if c(i) =. then delete;
   end;
   drop i;
run;

%let nsecond=0;
data work.check;
   set sashelp.vtable;
   if libname="WORK" & memname="DSNAME" & memtype="DATA" then do;
      call symput("nsecond", put(nobs, 6.));
      stop;
   end;
run;

%if &nfirst ne &nsecond %then %do;
   data _null_;
      file print;
      put ;
      put "    There were %eval(&nfirst-&nsecond) observation(s) deleted because of missing values." //;
      put ;
   run;
%end;

%if &nsecond=0 %then %goto endnoobs2;

%let nfirst=0;
/* no explanatory variable with value zero or less is allowed */
%if %length(&xnameo)>0 %then %do;
   %let help=%eval(&xmax.-&xmaxb);
   data work.check;
      set work.dsname;
      keep &xnameo;
      array c &xnameo;
      do i=1 to &help;
         if c(i)<=0 & c(i) ne . then do;
            call symput('nfirst','1');
            stop;
         end;
      end;
   run;
   %if &nfirst=1 %then %goto ENDZ;
%end;

* separate Alpha for each covariable;
%if %length(%qscan(&alpha, 2,%str( ))) ne 0 %then %do;
   %do i=1 %to &xmax;
      %let x=%qscan(&alpha, &i, %str( ));
      %let alpha&i=&x;
   %end;
%end;
%else %do;
   %do i=1 %to &xmax;
      %let alpha&i=&alpha;
   %end;
%end;

data work.a;
    length varname $8;
    %let i=1;
    %do %while (&i <= &xmax);
       varname = "&&mname&i";
       alpha = &&alpha&i;
       output;
       %let i=%eval(&i+1);
    %end;
run;

/*------------------------------------------------------------------------------
  ---
  --- basic model (without transformation)
  ---
  ------------------------------------------------------------------------------*/

%if %index("N",&showres)=0 %then %do;
   title 'MFP8: Model without transformation';

   %if %upcase(&model) = N %then
      %do;
      proc reg data=work.dsname;
         model &yname=&xinclude &xname;
      run;
   %end;

   %else %if %upcase(&model)=S %then
      %do;
      options nolabel;
      proc phreg data=work.dsname;
         model &yname*&censvar(&censval)=&xinclude &xname / ties=&ties risklimits;
         %if &strata ne %then %do;
            strata &strata;
         %end;
      run;
      options label;
   %end;

   %else %do;
      proc logistic data=work.dsname descending;
         model &yname=&xinclude &xname;
      run;
   %end;
%end;

/*------------------------------------------------------------------------------
  ---
  --- Ordering of covariable list XNAME if ORDER=1
  ---
  ------------------------------------------------------------------------------*/

%if &order=1 and &xmax>1 %then %do;

/*---  Different procedures depending on model ---*/

   %if %upcase(&model)=N %then
      %do;
      proc reg data=work.dsname outest=work._fpout outseb noprint;
         model &yname=&xinclude &xname;
      run;

      proc transpose data=work._fpout(drop=_RMSE_ Intercept) out=work._fpout;
      run;

      data work._fpout;
         set work._fpout;
         if index("&xname.",upcase(trim(_NAME_)))>0;
         rename COL1=est COL2=std;
      run;
   %end;

   %if %upcase(&model)=S %then
      %do;
      %let fobs=%eval(&xincmax+2);
      proc phreg data=work.dsname noprint outest=work._fpout(keep=_name_ &xname.) covout;
         model &yname*&censvar(&censval)=&xinclude &xname / ties=&ties;
         %if &strata ne %then %do;
            strata &strata;
         %end;
      run;

      data work._fpstd;
         set work._fpout(firstobs=&fobs);
         array cov &xname;
         std=sqrt(cov[_n_]);
         keep _name_ std;
      run;

      proc transpose data=work._fpout(obs=1) out=work._fpout(rename=(&yname=est));
      run;

      data work._fpout;
         merge work._fpout work._fpstd;
      run;
   %end;

   %if %upcase(&model)=L %then
      %do;
      %let fobs=%eval(&xincmax+3);
      proc logistic data=work.dsname noprint descending outest=work._fpout(keep=_name_ &xname.) covout;
          model &yname=&xinclude &xname / include=&xincmax;
      run;

      data work._fpstd;
         set work._fpout(firstobs=&fobs);
         array cov &xname;
         std=sqrt(cov[_n_]);
         keep _name_ std;
      run;

      proc transpose data=work._fpout(obs=1) out=work._fpout(rename=(&yname=est));
      run;

      data work._fpout;
         merge work._fpout work._fpstd;
      run;
   %end;

   data work._fpout;
      set work._fpout;
      rename _name_=name;
      ts=abs(est/std);
   run;

   data work._fpout;
      merge work._fpout a(keep=alpha);
   run;

   proc sort data=work._fpout;
      by descending ts;
   run;

   proc sql noprint;
      SELECT name INTO :xname SEPARATED BY " "
      FROM work._fpout;
      SELECT alpha INTO :alpha1-:alpha&xmax
      FROM _fpout;
   quit;

   %let xname=%trim(&xname);

   data _null_;
      %let i=1;
      %let x=%qscan(&xname,&i,%str( ));
      %do %while (%length(%trim(&x))>0);
         %let mname&i=&x;
         %let i=%eval(&i+1);
         %let x=%qscan(&xname,&i,%str( ));
      %end;
   run;

   data work.a;
     length varname $8;
     %let i=1;
     %do %while (&i <= &xmax);
        varname = "&&mname&i";
        alpha = &&alpha&i;
        output;
        %let i=%eval(&i+1);
     %end;
   run;

   title1;

   %if %index("N",&showres)=0 %then %do;
      data _null_;
         file print;
         put ;
         put "MFP8: IMPORTANT: The new ordering of the covariates based on p-values is:";
         put ;
         put "&xname";
         put ;
         put "It is used in the following calculations and outputs.";
      run;
   %end;
%end;

/*------------------------------------------------------------------------------
  ---
  --- Renaming the dataset and the variables, labeling the variables
  ---
  ------------------------------------------------------------------------------*/

* Macrovariables MVXi for the xmax variables Xi, initialised with Xi;
%do i=1 %to &xmax;
   %let mvx&i=X&i;
%end;

* Macrovariables degXi for the chosen degree of the FP, initialised with 0;
%do i=1 %to &xmax;
   %let degx&i=0;
%end;

%if %index("NR",&showres)=0 and &xmax ne 1 %then %do;
   data _null_;
      file print;
      put ;
      put "MFP8: Macrovariable settings:";
      put ;
      %do i=1 %to &xmax;
         put "MV X&i              : &&mvx&i";
      %end;
      put ;
   run;
%end;

data work._fp;
   set work.dsname;
   rename &yname=y;
   %let tlabel=label;
   %do i=1 %to &xmax;
      %let x=%qscan(&xname,&i,%str( ));
      rename &x=X&i;
      %let tlabel=&tlabel X&i="&x.";
   %end;
   %do i=1 %to &xincmax;
      %let x=%qscan(&xinclude,&i,%str( ));
      rename &x=inc&i;
      %let tlabel=&tlabel inc&i="&x.";
   %end;
run;

%let newinc=;
%do i=1 %to &xincmax;
%let newinc=&newinc inc&i;
%end;

data work._fp;
     set work._fp;
     &tlabel.;
run;

/*------------------------------------------------------------------------------
  ---
  --- Calculation of macro variables PMAX (= number of power values) and
  --- P1-P(PMAX) (= power values)
  ---
  ------------------------------------------------------------------------------*/

%let pmax=1;
%let x=%qscan(&pw,&pmax,%str( ));

%do %while (%length(%trim(&x))>0);
   %let p&pmax=&x;
   %if &x=1 %then %let plin=&pmax;
   %let pmax=%eval(&pmax+1);
   %let x=%qscan(&pw,&pmax,%str( ));
%end;
%let pmax=%eval(&pmax-1);

* Macrovariables mdpi for the allowed degree of the FP;

%do i=1 %to &xmax;
   %let mdp&i=&m;
%end;

%if %length(%trim(&dp)) > 0 %then %do;
   %let i=1;
   %let y=%qscan(&dp,&i,%str( ));
   %do %while (%length(%trim(&y))>0);
      %do k=1 %to &xmax;
         %let x=%qscan(&xname,&k,%str( ));
         %if %index(%upcase(&y),%upcase(&x)) ne 0 %then %let mdp&k = %qscan(&y,2,":");
      %end;
      %let i=%eval(&i+1);
      %let y=%qscan(&dp,&i,%str( ));
   %end;
%end;

* Macrovariables mpwi for the allowed powers of the FP;
%do i=1 %to &xmax;
   %let mpw&i=;
%end;
%if %length(%trim(&power)) > 0 %then %do;
   %let i=1;
   %let y=%qscan(&power,&i,"]");
   %do %while (%length(%trim(&y))>0);
      %do k=1 %to &xmax;
         %let x=%qscan(&xname,&k,%str( ));
         %if %index(%upcase(&y),&x) ne 0 %then %do;
            %let l=1;
            %let help=%qscan(&y,2,"[");
            %let help2=%qscan(&help,&l,%str( ));
            %do %while (%length(%trim(&help2))>0);
               %do r=1 %to &pmax;
                  %if &&p&r=&help2 %then %let mpw&k=&&mpw&k &r;
               %end;
               %let l=%eval(&l+1);
               %let help2=%qscan(&help,&l,%str( ));
            %end;
         %end;
      %end;
      %let i=%eval(&i+1);
      %let y=%qscan(&power,&i,"]");
   %end;
%end;

/*------------------------------------------------------------------------------
  ---
  --- Calculation of Xi**P
  ---
  ------------------------------------------------------------------------------*/

data work._fp;
   set work._fp;
   %do i=1 %to &xmax;
       %xtop(X&i,&pref,&&mdp&i);
   %end;
run;

/*------------------------------------------------------------------------------
  ---
  --- Generation of formats
  ---    FPPW   for the power values
  ---    FPFORM for functional forms
  ---    FPSIGN for significance result
  ---
  ------------------------------------------------------------------------------*/

%let name=WORK.FORMATS.FPPW.FORMAT;
proc format;
   value fppw
   %do i=1 %to &pmax;
      &i=&&p&i
   %end;
   ;
run;

%let name=WORK.FORMATS.FPFFORM.FORMAT;
proc format;
   value fpfform -1="Omitted"
                  0="Linear"
                  1="First Degree"
                  2="Second Degree"
                  3="Third Degree"
   ;
run;

%let name=WORK.FORMATS.FPSIGN.FORMAT;
proc format;
   value fpsign 1="*"
                0="-"
   ;
run;

/*------------------------------------------------------------------------------
  ---
  --- Calculation of different models for all Xi combinations
  ---
  ------------------------------------------------------------------------------*/

%do ic=1 %to &maxcycle;
   %if %index("N",&showres)=0 %then %do;
      data _null_;
         file print;
         put ">>> -------------------- CYCLE &ic --------------------";
      run;
   %end;

   %do ix=1 %to &xmax;

      %let xbase=%unquote(%xvars(&xmax,&ix,%nrstr(&mvx),));
      %fpmodels(&model,y,X&ix,&pref,&newinc &xbase,&&mdp&ix,&strata);

         *---;
         *--- Determination of standard errors of parameter estimates;
         *---;

      %let help=&newinc &xbase X&ix.;
      %if &&mdp&ix>0 %then %let help=&help X&ix.&pref.1 -X&ix&pref%eval(&pmax*&&mdp&ix);

      %if %upcase(&model)=N %then %do;

         data work._fpstd work._fpout;
            set work._fpout;
            if _type_="PARMS" then output _fpout;
            else output _fpstd;
         run;

         data work._fpout;
             set work._fpout;
             if _sse_<0.000000000001 then _sse_=0.000000000001;
         run;

         data work._fpstd;
            set work._fpstd;
            i=input(substr(_model_,6,3),3.);
            keep i intercept X:;
         run;

         proc transpose data=work._fpstd out=work._fpstd(where=(col1 ne .));
               by i;
         run;

         %if &ix>9 %then %let dez=3;
         %else %let dez=2;

         data work._fpstd;
            set work._fpstd;
            retain j;
            by i;
            if first.i then j=0;
            if substr(_NAME_,1,&dez) = "X&ix" and length(scan(_NAME_,1,"&pref")) = &dez
               then do;
                  j=j+1;
                  _NAME_=compress("SE"||put(j,3.));
               end;
            else _NAME_=compress(substr(compress("SE"||_NAME_)||"      ",1,8)," ");
            drop j _LABEL_;
         run;
         proc transpose data=_fpstd out=_fpstd(drop=i _NAME_);
            by i;
         run;
         data work._fpout;
            merge work._fpout work._fpstd;
         run;

      %end;

      %if %upcase(&model)=L %then %do;
  
         data work._fpstd work._fpout;
            set work._fpout;
            retain i 0;
            if _type_="PARMS" then i=i+1;
            if _type_="PARMS" then output _fpout;
            else output _fpstd;
         run;
   
         proc transpose data=work._fpstd(drop=_LINK_ _TYPE_ _STATUS_ _LNLIKE_) out=work._fpstd;
            by i;
         run;

         data work._fpstd;
            retain j 0;
            set work._fpstd;
            by i;
            if first.i then j=1;
            else j=j+1;
         run;

         data work._fpstd;
            set work._fpstd;
            array v interc &help;
            vari=v[j];
            if vari ne .;
            se=vari**0.5;
            drop interc &help vari;
         run;

         %if &ix>9 %then %let dez=3;
         %else %let dez=2;

         data work._fpstd;
            retain j 0;
            set work._fpstd(drop=j);
            by i;
            if first.i then j=0;
            if substr(_name_,1,&dez) = "X&ix" and length(scan(_name_,1,"&pref")) = &dez
               then do;
                  j=j+1;
                  _name_=compress("SE"||put(j,3.));
               end;
            else _NAME_=compress(substr(compress("SE"||_NAME_)||"      ",1,8)," ");
            drop j;
         run;

         proc transpose data=work._fpstd out=work._fpstd(drop=i _name_);
            by i;
         run;

         data work._fpout;
            merge work._fpout work._fpstd;
         run;
      %end;

      %if %upcase(&model)=S %then %do;
         data work._fpstd work._fpout;
            set work._fpout;
            retain i 0;
            if _type_="PARMS" then i=i+1;
            if _type_="PARMS" then output _fpout;
            else output work._fpstd;
         run;

         proc transpose data=work._fpstd(drop=_status_ _type_ _ties_ _lnlike_) out=work._fpstd;
            by i;
         run;

         data work._fpstd;
            retain j 0;
            set work._fpstd;
            by i;
            if first.i then j=1;
            else j=j+1;
         run;

         data work._fpstd;
            set work._fpstd;
            array v &help;
            vari=v[j];
            if vari ne .;
            se=vari**0.5;
            drop &help vari;
         run;

         %if &ix>9 %then %let dez=3;
         %else %let dez=2;

         data work._fpstd;
            retain j 0;
            set work._fpstd(drop=j);
            by i;
            if first.i then j=0;
            if substr(_name_,1,&dez) = "X&ix" and length(scan(_name_,1,"&pref")) = &dez
               then do;
                  j=j+1;
                  _name_=compress("SE"||put(j,3.));
               end;
            else _name_=compress(substr(compress("SE"||_name_)||"      ",1,8)," ");
            drop j;
         run;

         proc transpose data=work._fpstd out=work._fpstd(drop=_name_);
            by i;
         run;
         data work._fpout;
            merge work._fpout work._fpstd;
            by i;
         run;

      %end;

      *---;
      *--- Determination of deviance;
      *---;

      data work._fp1;
         set work._fpout;
         format deviance BEST12.;
         %if %upcase(&model)=N %then
            %do;
            n=_p_+_edf_;
            deviance=n*(1+log(2*arcos(-1)*_sse_/n));
         %end;
         %else
            %do;
            deviance=(-2)*_lnlike_;
         %end;

         %if %upcase(&model)=S %then
            %do;
            intercept=0;
         %end;

         keep intercept &help deviance se:;
      run;

      *---;
      *--- Differences in Deviance;
      *---;

      %let help2= p1;
      %if &&mdp&ix>1 %then %let help2=&help2.-p&&mdp&ix;
      %let help3= beta1;
      %if &&mdp&ix>1 %then %let help3=&help3.-beta&&mdp&ix;

      %if &&mdp&ix>0 %then %do;
         data work._fp2;
            set work._fp1;
            array vx X&ix._1-X&ix._%eval(&pmax*&&mdp&ix);
            array vb b1-b%eval(&pmax*&&mdp&ix);

            do i=1 to %eval(&pmax*&&mdp&ix);
               if vx[i]=. then
                  do;
                  vb[i]=0;
                  vx[i]=0;
               end;
               else
                  do;
                  vb[i]=vx[i];
                  vx[i]=1;
               end;
            end;
            m=sum(of X&ix._1-X&ix._%eval(&pmax*&&mdp&ix));

            array p p1-p&m;            *--- Variables p1 to pm;
            j=1;i=1;

            if m>0 then do;
               do until (i>%eval(&pmax*&&mdp&ix) or m<j);
                  if vx[i]=1 then do;
                     p[j]=mod(i-1,&pmax)+1;
                     j=j+1;
                  end;
                  i=i+1;
               end;
            end;

            if _n_=1 then
               do;
               m=-1;
               call symput('devwox1',deviance);
            end;
            if X&ix ne . then call symput('devlin',deviance);
         run;
     %end;
     %else %do;
         data work._fp2;
            set work._fp1;
            m=_n_-2;
            array p p1 p&m;
            j=1;i=1;
            if _n_=1 then
               do;
               m=-1;
               call symput('devwox1',deviance);
            end;
            if X&ix ne . then call symput('devlin',deviance);
         run;
      %end;

      %if %index("NR",&showres)=0 %then %do;
         %let adevlin=%substr(&devlin,1,7);
         %let adevwox1=%substr(&devwox1,1,7);
         data _null_;
            file print;
            put ;
            put "MFP8:";
            put "*** Deviance of linear model     = &adevlin ***";
            put "*** Deviance of model without x&ix = &adevwox1 ***";
            put ;
         run;
      %end;

      /* the models which are not wanted are deleted */
      %if %length(%trim(&&mpw&ix)) > 0 %then %do;
         %do k=1 %to &pmax;
            %if %index(&&mpw&ix,&k)=0 %then %do;
               data work._fp2;
                  set work._fp2;
                  if p1 = &k then delete;
                  if p2 = &k then delete;
               run;
            %end;
         %end;
      %end;

      data work._fp2;
         set work._fp2;
         format difflin difwox1 BEST12.;
         devwox1=input(symget('devwox1'),BEST12.);
         devlin=input(symget('devlin'),BEST12.);
         difflin=devlin-deviance;
         difwox1=devwox1-deviance;
         rename intercept=b0;
         drop i j;
      run;

      %if %index("NR",&showres)=0 %then %do;
         title 'MFP8: Deviances of all models';
         proc print data=work._fp2;
            var m &help2 deviance difwox1 difflin;
            format &help2 fppw.;
         run;
      %end;

   /*------------------------------------------------------------------------------
     ---
     --- Determination of Best Functions per Degree
     ---
     ------------------------------------------------------------------------------*/

      proc sort data=work._fp2;
         by m descending difflin;
      run;

      /* determination of the residual degree of freedom in case of normal distribution */
      %let dfres=%eval(&nsecond-1-&xincmax);
      %if %upcase(&model)=N %then %do;
         %do i=1 %to &xmax;
            %if &ix ne &i %then %do;
               %if &&degx&i = 0 %then %let dfres=%eval(&dfres-1);
               %if &&degx&i = 1 %then %let dfres=%eval(&dfres-2);
               %if &&degx&i = 2 %then %let dfres=%eval(&dfres-4);
            %end;

            data work._fp2;
               set work._fp2;
               dfres=&dfres;
               n=&nsecond;
            run;
         %end;
      %end;

      data work._fp3;
         set work._fp2;
         by m descending difflin;
         if first.m;
            diffseq=abs(lag(deviance)-deviance);
      run;

      /* within the linear model the F-test is used */
      %if %upcase(&mselect) = SEQ and %upcase(&model)=N %then %do;
         data work._fp3;
            set work._fp3;
            if diffseq ne . then do;
               z=(exp(diffseq/n)-1)*(dfres-max(m,1))/max(m,1);
               if z<0.00001 then z=0;
               pdiffdev=1-probf(z, max(m,1), dfres-max(m,1));
               if pdiffdev<&&alpha&ix then sign=1;
               else sign=0;
            end;
         run;
      %end;

      %if %upcase(&mselect) = SEQ and %upcase(&model) ne N %then %do;
         data work._fp3;
            set work._fp3;
            if diffseq ne . then do;
               if diffseq<0.00001 then z=0;
               else z=diffseq;
               pdiffdev=1-probchi(z,max(m,1));
               if pdiffdev<&&alpha&ix then sign=1;
               else sign=0;
            end;
         run;
      %end;

      %if %upcase(&mselect) = RA2 %then %do;

         proc sql noprint;
            SELECT MAX(m) INTO :maxm
            FROM work._fp3;
            SELECT deviance INTO :deviance
            FROM work._fp3
            WHERE m = &maxm;
         quit;

         data work._fp3;
            set work._fp3;
            diffra2=abs(deviance-&deviance);
            df=&maxm-m;
         run;

         %if &maxm=2 %then %do;
            data work._fp3;
               set work._fp3;
               df=df+1;
            run;
         %end;

         %if %upcase(&model)=N %then %do;
         data work._fp3;
            set work._fp3;
               z=(exp(diffra2/n)-1)*(dfres-max(df,1))/max(df,1);
               if z<0.00001 then z=0;
               pdiffdev=1-probf(z, max(df,1), dfres-max(df,1));
               if pdiffdev<&&alpha&ix then sign=1;
               else sign=0;
            run;
         %end;
         %if %upcase(&model) ne N %then %do;
         data work._fp3;
            set work._fp3;
            if diffra2<0.00001 then z=0;
            else z=diffra2;
            pdiffdev=1-probchi(z,max(df,1));
            if pdiffdev<&&alpha&ix then sign=1;
            else sign=0;
         run;
         %end;
      %end;

   /*------------------------------------------------------------------------------
     ---
     --- Determination of estimates BETA
     ---
     ------------------------------------------------------------------------------*/

     %if &&mdp&ix=0 %then %do;
         data work._fp4;
            set work._fp3;
            beta1=X&ix;
            fform=m;
         run;
     %end;
     %else %do;
         data work._fp4;
            set work._fp3;
            fform=m;
            if m>-1 then
            do;
               if m=0 then beta1=x&ix;
               else do;
                  array x X&ix.&pref.1-X&ix.&pref%eval(&pmax*&&mdp&ix);
                  array b b1-b%eval(&pmax*&&mdp&ix);
                  array beta &help3;
                  j=1;
                  do i=1 to %eval(&pmax*&&mdp&ix);
                     if x[i]=1 then
                     do;
                        beta[j]=b[i];
                        j=j+1;
                     end;
                  end;
               end;
            end;
         drop i j;
         run;
      %end;

   /*------------------------------------------------------------------------------
     ---
     --- Showing the Results in Output Window
     --- (at first creating the labels)
     ---
     ------------------------------------------------------------------------------*/

      %let lab=;
      %do j=1 %to &xmax;
         %let lab=&lab %labs(&j, &&degx&j, &&mvx&j, &ix, &xname, 2);
      %end;

      %if %index("N",&showres)=0 %then %do;
         title1 "MFP8: Variable -&&mname&ix- ";
         title2 "Best Functions for Different Degrees m ";

         %if %upcase(&mselect) = SEQ %then %do;
         proc print data=work._fp4 label noobs;
            var fform m p1-p&m deviance diffseq pdiffdev;
            format fform fpfform. p1-p&m fppw. deviance diffseq 10.3;
            label fform='Function' ;
         run;
         %end;
         %else %do;
         proc print data=work._fp4 label noobs;
            var fform m p1-p&m deviance diffra2 pdiffdev;
            format fform fpfform. p1-p&m fppw. deviance diffra2 10.3;
            label fform='Function' ;
         run;
         %end;
      %end;

      %if %index("NR",&showres)=0 %then %do;

         title1 "MFP8: Variable  -&&mname&ix- ";
         title2 "Estimates for Best Functions for Different Degrees m";

         proc print data=work._fp4 label noobs;
            var fform m p1-p&m
            %if &model ne S %then b0;
            &newinc &xbase &help3 se:;
            format fform fpfform. p1-p&m fppw. &newinc &xbase
            %if &model ne S %then b0;
            &help3 se: 13.6;
            label fform="Function" %exinc(&xinclude, &xincmax, 1) &lab %exlabb(&xname, &ix, &&mdp&ix, 1)
            %if &model ne S %then b0="Intercept" SEINTERC="Intercept(SE)";;
         run;
      %end;

   /*------------------------------------------------------------------------------
     ---
     --- Decision which Function is best
     ---
     --- MSELECT=ra2 (Stata 6, Def.): m=2 vs. m=-1, m=2 vs. m=0, m=2 vs. m=1
     --- MSELECT=seq (Stata 6, Seq.): m=2 vs. m=1 , m=1 vs. m=0, m=0 vs. m=-1
     ---
     ------------------------------------------------------------------------------*/

      %if %upcase(&mselect) = SEQ %then %do;
         proc sort data=work._fp4 out=work._fp5;
            by descending m;
         run;
      %end;
      %else %do;
         data work._fp5;
            set work._fp4;
         run;
      %end;

      %if %upcase(&mselect)=RA2 %then %do;
         data _null_;
            set work._fp5;
            if sign=0 then
               do;
               call symput("dm",m);
               stop;
            end;
         run;
      %end;

      %if %upcase(&mselect)=SEQ %then %do;
         %let dm=x;
         data _null_;
            set work._fp5;
            if sign=1 then
               do;
               call symput("dm",m);
               stop;
            end;
         run;
         %if &dm=x %then %let dm=-1;
      %end;

      data work._fp6;
         set work._fp5(where=(m=&dm));
      run;

      %if %index("NR",&showres)=0 %then %do;
         title1 "MFP8: Best Functions found";
         proc print data=work._fp6 label noobs;
            var fform m p1-p&m
            %if &model ne S %then b0;
            &newinc &xbase &help3 se:;
            format fform fpfform. p1-p&m fppw. &newinc &xbase
            %if &model ne S %then b0;
            &help3 se: 13.6;
            label fform="Function" %if &model ne S %then b0="Intercept" SEINTERC="Intercept(SE)"; &lab %exlabb(&xname, &ix, &&mdp&ix, 1) %exinc(&xinclude, &xincmax, 1);
         run;
      %end;

   /*------------------------------------------------------------------------------
     ---
     --- Generation of result data set
     ---
     ------------------------------------------------------------------------------*/

      data work._fp7;
         sim=&ic;
         length xneu $20 varname $3;
         varname="X&ix";
         set work._fp6;

         if m=-1 then status="out";
         else
            do;
            status="in";
            if m=0 then xneu=varname;
            else if m=1 then xneu=compress(varname||"&pref"||p1," ");
            else if m=2 then
               do;
               temp=p1+&pmax;
               if p1=p2 then
                  xneu=compress(varname||"&pref"||p1," ")||" "||compress(varname||"&pref"||temp," ");
               else
                  xneu=compress(varname||"&pref"||p1," ")||" "||compress(varname||"&pref"||p2," ");
               drop temp;
            end;
         end;
         call symput("degx&ix",m);
         call symput("MVX&ix",xneu);
      run;

      %if &ix=1 and &ic=1 %then
         %do;
         data work._fp8;
            set work._fp7(drop=X&ix.:);
         run;
      %end;
      %else
         %do;
         data work._fp8;
            set work._fp8 work._fp7(drop=X&ix.:);
            if sim=&ic-2 then delete;
         run;
      %end;
      %if %index("NR",&showres)=0 and &xmax ne 1 %then %do;
         data _null_;
            file print;
            put ;
            put "MFP8: Macrovariable settings, cycle: &ic, variable &ix";
            put ;
            %do i=1 %to &xmax;
               put "MV X&i              : &&mvx&i";
            %end;
            put ;
         run;
      %end;
   %end;


   /*------------------------------------------------------------------------------
     ---
     --- Decision of ongoing iteration (only for more than one covariable)
     ---
     ------------------------------------------------------------------------------*/

   title1;

   %if &ic ne 1 %then
      %do;
      proc sort data=work._fp8(keep=sim varname xneu where=(sim=&ic or sim=&ic-1)) out=work._fp9;
         by varname sim;
      run;

      proc transpose data=work._fp9 out=work._fp9;
         by varname;
         var xneu;
      run;

      data work._fp9;
         set work._fp9 end=eof;;
         retain diff 0;
         if col1 ne col2 then diff+1;
         if eof then
            do;
            if diff=0 then simdiff="E";
            else simdiff="D";
         end;
         call symput("simdiff",simdiff);
      run;

      %if &simdiff=E %then %do;
         %if %index("N",&showres)=0 %then %do;
            data _null_;
               file print;
               put ;
               put ">>> ---------------------------------------------------" ;
               put ;
               put ">>> MFP8 converged after &ic cycles." //;
               put ;
            run;
         %end;
      %goto ende;
      %end;
   %end;

%end;

%if &xmax ne 1 %then %do;
   %let ic=%eval(&maxcycle-1);
   data _null_;
         file print;
         put ;
         put ">>> MFP8 finished without convergence" //;
         put ;
   run;
%end;
%else %let ic=1;

%ende:
   %let help=0;
   %do i=1 %to %eval(&xmax-1);
     %if &&degx&i=-1 %then %let help=%eval(1+&help);
   %end;

   %let lab=;
   %do i=1 %to &xmax;
      %let lab= &lab %labs(&i, %eval(&&degx&i), &&mvx&i, 0, &xname, 0);
   %end;

   title "MFP8: Final multivariable fractional polynomial model";

   data work.incl;
      length varname $8;
      %let i=1;
      %do %while (&i <= &xincmax);
         varname = "%qscan(&xinclude, &i, %str( ))";
         fform = 0;
         status="forced in";
         alpha = 1;
         output;
         %let i=%eval(&i+1);
      %end;
   run;

   data work.aus1;
      set work._fp8;
      where sim=%eval(&ic);
      keep fform p1-p&m status deviance;
   run;

   data work.aus1;
       merge work.aus1 work.a;
       if fform=0 then p1=&plin;
   run;

   %if &xincmax >0 %then %do;
      data work.aus1;
         set work.incl work.aus1;
         if fform=0 then p1=&plin;
      run;
   %end;

   data &dsout.1;
      set work.aus1;
      format fform fpfform. p1-p&m fppw.;
      label varname="Variable" fform="Function" alpha="Alpha" status="Status";
   run;

   %if %index("N",&showres)=0 %then %do;
       proc print data=work.aus1 label noobs;
         var varname fform alpha status p1-p&m;
         format fform fpfform. p1-p&m fppw.;
         label varname="Variable" fform="Function" alpha="Alpha" status="Status";
      run;
   %end;

   data work.std;
      set work._fp6;
      drop x&xmax:;
   run;

   data work.std;
      set work.std;
      keep %if &xincmax > 0 %then inc:;
      %if &help < %eval(&xmax.-1) %then x:;
      se: beta: %if &model ne S %then b0;;
      rename %brename(&&mvx&xmax, 1) %if &model ne S %then b0=intercep;;
   run;

   data work.par;
      set work.std;
      drop se:;
      label &lab  %if &xincmax > 0 %then %exinc(&xinclude, &xincmax, 0);
      %if &model ne S %then intercept="Intercept";;
   run;

   data work.std;
      set work.std;
      keep se:;
   run;

   proc transpose data=work.par out=work.par;

   %if %length(&lab)>0 or &xincmax>0 or &model ne S%then %do;
      data work.par1;
         set work.par;
         rename _label_=variable col1=Coeff;
      run;
   %end;
   %else %do;
      data work.par1;
      run;
   %end;

   proc transpose data=work.std out=work.std;
   data work.std1;
      set work.std;
      rename col1=StdErr;
   run;

   data work.aus2;
      merge work.par1 work.std1;
      drop _name_;
   run;

   title;

   %if &&degx&xmax=-1 %then %let help=%eval(&help+1);

   %if &model = S %then %do;
      data work.aus2;
         set work.aus2;
         hr=exp(coeff);
         pwert=1-probchi((coeff/stderr)**2,1);
         ll=exp(coeff-1.96*stderr);
         ul=exp(coeff+1.96*stderr);
      run;
   %end;
   %if &model = L %then %do;
      data work.aus2;
         set work.aus2;
         pwert=1-probchi((coeff/stderr)**2,1);
         ll=coeff-1.96*stderr;
         ul=coeff+1.96*stderr;
      run;
   %end;
   %if &model = N %then %do;
       data work.aus2;
           set work.aus2;
           pwert=2*(1-probt(abs(coeff)/stderr,%eval(&nsecond-1)));
           ll=coeff-1.96*stderr;
           ul=coeff+1.96*stderr;
      run;
   %end;

   data work.aus1;
       set work.aus1;
       %if &model = S %then number = _n_;
       %else number=_n_+1;;
   run;

   data work.aus2;
       set work.aus2;
       where coeff ne .;
   run;

   %funcfm(&pw, &pmax, &model);

   %if &help < %eval(&xmax.+&xincmax) or &model ne S %then %do;
      data &dsout.2;
         set work.aus2;
         keep variable coeff stderr pwert %if &model = S %then hr; ll ul;
         label variable="Variable" Coeff="Coefficient" StdErr="Standard Error"
             %if &model ne N %then pwert="Pr > ChiSq"; %else pwert="Pr > |t|";
             %if &model = S %then hr="Hazard Ratio";
             ll="%nrbquote(95%) Confidence" ul="Limits";
      run;

      %if %index("N",&showres)=0 %then %do;
         proc print data=work.aus2 label noobs;
            var variable coeff stderr pwert %if &model = S %then hr; ll ul;
            format Coeff StdErr 13.5 pwert %if &model = S %then hr; ll ul 13.4;
            label Variable="Variable" Coeff="Coefficient" StdErr="Standard Error"
               %if &model ne N %then pwert="Pr > ChiSq"; %else pwert="Pr > |t|";
               %if &model = S %then hr="Hazard Ratio";
               ll="%nrbquote(95%) Confidence" ul="Limits";
         run;
      %end;
   %end;

   data _null_;
      set work.aus1;
      if deviance ne . then do;
         loglik=-0.5*deviance;
         call symput("loglik", put(loglik, 10.2));
         stop;
      end;
   run;

   %if %index("N",&showres)=0 %then %do;
      data _null_;
         file print;
         put "Log(Likelihood) of the resulting model:  &loglik";
         put;
      run;
   %end;

   %goto schluss;

%endlib:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    The specified dataset is located in the working directory.";
      put "    This directory will be emptied before the macro starts.";
      put "    Please save the dataset within another library!" //;
      put ;
   run;
   %goto schluss;

%endz:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    At least one explanatory variable which is used has value zero or less. ";
      put "    Please transform this variable!" //;
      put ;
   run;
   %goto schluss;

%endm:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    Invalid parameter for Macro variable MODEL. ";
      put "    Valid values are N, S and L!" //;
      put ;
   run;
   %goto schluss;

%ends:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    The survival time must be greater than or equal to zero!" //;
      put ;
   run;
   %goto schluss;

%ends2:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    A censoring variable must be specified!" //;
      put ;
   run;
   %goto schluss;

%endl:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    The response must take values in {0, 1}!" //;
      put ;
   run;
   %goto schluss;

%endnox:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    At least one response variable has to be specified!" //;
      put ;
   run;
%goto schluss;

%endnoobs:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    The dataset was not found.";
      put "    Probably the library is not properly assigned. ";
      put ;
   run;
%goto schluss;

%endnoobs2:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    After deleting observations with missing values no observations remained.";
      put ;
   run;
%goto schluss;

%endx:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 cancelled!" //;
      put "    Not all variables you have specified are found in the data!" //;
      put ;
   run;
   %goto schluss;

%schluss:
   data _null_;
      file print;
      put ;
      put ">>> MFP8 has finished!";
      put ;
   run;
%mend mfp8;