/*----------------------------------------------------------------------

The NLEST macro allows you to estimate one or more linear or nonlinear
combinations of parameters from any model for which you can save the model
parameters and their variance-covariance matrix. Most modeling procedures which
offer ESTIMATE, CONTRAST, or LSMEANS statements only provide for estimating or
testing linear combinations of model parameters. However, common estimation
problems often involve nonlinear combinations, particularly in generalized
models with nonidentity link functions such as logistic and Poisson models.

In addition to Base SASŪ, SAS/STATŪ is required to run this macro.

For full documentation of the NLEST macro, including descriptions of all 
available options and examples, see http://support.sas.com/kb/58775.

------------------------------------------------------------------------

DISCLAIMER:

THIS INFORMATION IS PROVIDED BY SAS INSTITUTE INC. AS A SERVICE
TO ITS USERS.  IT IS PROVIDED "AS IS".  THERE ARE NO WARRANTIES,
EXPRESSED OR IMPLIED, AS TO MERCHANTABILITY OR FITNESS FOR A
PARTICULAR PURPOSE REGARDING THE ACCURACY OF THE MATERIALS OR CODE
CONTAINED HEREIN.

----------------------------------------------------------------------*/

%macro NLEST(version, instore=, inest=, incovb=, where=, f=, fdata=, 
       label=, score=, outscore=NLEst, df=, alpha=0.05, listnames=yes,
       title=, print=yes, null=0, covdrop=);

%macro NLEstimate(version, instore=, inest=, incovb=, where=, print=yes,
                  f=, fdata=, label=, score=, outscore=NLEst, null=0,
                  df=, alpha=0.05, listnames=yes, title=, covdrop=);

   %let time = %sysfunc(datetime());
   %local notesopt; 
   %let notesopt = %sysfunc(getoption(notes));
   %let _version=2.1;
   options notes;
     %if &version ne %then %put NOTE: NLEST macro Version &_version..;
   %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;;
   %let status=ok;
   %let newchk=1;
   %let PInit =;
   %let nrecdiv = 300;
   %let Expnt=;
   %let Expnt1=;
   %if %index(%upcase(&version),DEBUG) %then %do;
     options notes mprint
       %if %index(%upcase(&version),DEBUG2) %then mlogic symbolgen;
     ;  
     ods select all;
     %put _user_;
   %end;
   %else %do;
     %if %index(%upcase(&version),NOCHK) %then %let newchk=0;
     options nonotes nomprint nomlogic nosymbolgen;
     ods exclude all;
   %end;

   /* 
   /  Check for newer version. 
   /------------------------------------------------------------------*/
   %if &newchk %then %do;
      %let _notfound=0; %let _newver=0;
      filename _ver url 'http://ftp.sas.com/techsup/download/stat/versions.dat' 
               termstr=crlf;
      data _null_;
        infile _ver end=_eof;
        input name:$15. ver;
        if upcase(name)="&sysmacroname" then do;
          call symput("_newver",ver); stop;
        end;
        if _eof then call symput("_notfound",1);
        run;
      options notes;
      %if &syserr ne 0 or &_notfound=1 or &_newver=0 %then
        %put NOTE: Unable to check for newer version of the NLEST macro.;
      %else %if %sysevalf(&_newver > &_version) %then %do;
        %put NOTE: A newer version of the NLEST macro is available at;
        %put NOTE- this location: http://support.sas.com/ ;
      %end;
      %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;;
   %end;


   /*
   /  Input checking.
   /------------------------------------------------------------------*/
   %if (&instore= and &inest= and &incovb=) or
       (&instore ne and (&inest ne or &incovb ne)) %then %do;
     %put ERROR: Either INSTORE= or both INEST= and INCOVB= must be specified.;
     %let status=input_err;
     %goto exit;
   %end;
   %if (&inest ne and &incovb=) or (&inest= and &incovb ne) %then %do;
     %put ERROR: Both INEST= and INCOVB= must be specified.;
     %let status=input_err;
     %goto exit;
   %end;
   %if %index(&f,%str(%"))>0 or %index(&label,%str(%"))>0 %then %do;
     %put ERROR: Do not use double quotes (%str(%")) in F= or LABEL=.;
     %put ERROR- Use single quotes (%str(%')).;
     %let status=input_err;
     %goto exit;
   %end;
   %if %index(%quote(&title),%str(%")) or %index(%quote(&title),%str(%')) %then %do;
     %put ERROR: Do not use quotes (%str(%") or %str(%')) in TITLE=.;
     %let status=input_err;
     %goto exit;
   %end;
   %if %index(%upcase(&version),SHOWNAMES)=0 and &fdata= and %quote(&f)= %then %do;
     %put ERROR: Either or both of FDATA= and F= must be specified.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &inest ne %then %do;
     %if %sysfunc(exist(&inest)) ne 1 %then %do;
       %put ERROR: INEST= data set not found.;
       %let status=input_err;
       %goto exit;
     %end;
     %else %do;
       %let dsid=%sysfunc(open(&inest));
       %if &dsid %then %do;
         %if %sysfunc(varnum(&dsid,estimate))=0 %then %do;
           %put ERROR: Required variable, ESTIMATE, not found.;
           %let status=input_err;
         %end;
         %let rc=%sysfunc(close(&dsid));
         %if &status=input_err %then %goto exit;
       %end;
       %else %do;
         %put ERROR: Could not open INEST= data set.;
         %goto exit;
       %end;
     %end;
   %end;
   %if &incovb ne and %sysfunc(exist(&incovb)) ne 1 %then %do;
     %put ERROR: INCOVB= data set not found.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &score ne and %sysfunc(exist(&score)) ne 1 %then %do;
     %put ERROR: SCORE= data set not found.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &score ne and &fdata ne %then %do;
     %put ERROR: When SCORE= is specified FDATA= is not allowed. Specify F=.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &score ne and %quote(&f)= %then %do;
     %put ERROR: When SCORE= is specified F= must also be specified.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &instore ne and %sysfunc(exist(&instore,ITEMSTOR)) ne 1 %then %do;
     %put ERROR: INSTORE= item store not found.;
     %let status=input_err;
     %goto exit;
   %end;
   %if &fdata ne %then %do;
     %if %sysfunc(exist(&fdata)) ne 1 %then %do;
       %put ERROR: FDATA= data set not found.;
       %let status=input_err;
       %goto exit;
     %end;
     %else %do;
       %let dsid=%sysfunc(open(&fdata));
       %if &dsid %then %do;
         %if %sysfunc(varnum(&dsid,label))=0 %then %do;
           %put ERROR: Required variable, LABEL, not found.;
           %let status=input_err;
         %end;
         %if %sysfunc(varnum(&dsid,f))=0 %then %do;
           %put ERROR: Required variable, F, not found.;
           %let status=input_err;
         %end;
         %let rc=%sysfunc(close(&dsid));
         %if &status=input_err %then %goto exit;
       %end;
       %else %do;
         %put ERROR: Could not open FDATA= data set.;
         %goto exit;
       %end;
     %end;
   %end;
   %let listnames=%upcase(%substr(&listnames,1,1));
   %if &listnames ne Y and &listnames ne N %then %do;
     %put ERROR: The LISTNAMES= value must be either YES or NO.;
     %goto exit;
   %end;
   %if %quote(&df) ne %then %do;
     %if (%sysfunc(verify(%sysfunc(trim(%sysfunc(left(&df)))),'0123456789.'))>0) + 
         (%sysfunc(count(&df,.))>1) %then %do;
       %put ERROR: The DF= value must be a non-zero, positive value.;
       %goto exit;
     %end;
     %if  %sysevalf(&df=0 or &df=.) %then %do;
       %put ERROR: The DF= value must be a non-zero, positive value.;
       %goto exit;
     %end;
   %end;
   %if %sysevalf(&alpha<=0 or &alpha>=1) %then %do;
     %put ERROR: The ALPHA= value must be between 0 and 1.;
     %goto exit;
   %end;
   %if %quote(&null)= %then %let null=0;
   %else %do;
     %if (%sysfunc(verify(%sysfunc(trim(%sysfunc(left(&null)))),'0123456789.-'))>0) + 
         (%sysfunc(findc(%sysfunc(trim(%sysfunc(left(&null)))),'-'))>1) + 
         (%sysfunc(count(&null,-))>1) + (%sysfunc(count(&null,.))>1) %then %do;
           %put ERROR: The NULL= value must be a numeric value.;
           %goto exit;
     %end;
     %if %quote(&null)=. or %quote(&null)=%quote(-) %then %do;
         %put ERROR: The NULL= value must be a numeric value.;
         %goto exit;
     %end;
   %end;

   /*
   /  Put parameter estimates and covariance matrics into data sets.
   /------------------------------------------------------------------*/
   %if &instore ne %then %do;
     proc plm restore=&instore;
        show parameters covariance;
        ods output ParameterEstimates=_Parm
                   Cov               =_Cov;
     run;
     %if &syserr %then %do;
       %let status=plm_err;
       %goto exit;
     %end;
   %end;
   %else %do;
     data _Parm; 
       set &inest; 
       run;
     data _Cov; 
       set &incovb; 
       run;
   %end;
   
   data _Parm; 
      set _Parm; 
      %if %quote(&where) ne %then where &where;;
      run;
   data _Parm; 
      set _Parm nobs=nobs; 
      call symput("dim",cats(nobs));
      Row=cats("p",_N_);
      run;
   data _Cov; 
      set _Cov; 
      %if %quote(&where) ne %then where &where;;
      array v (*) _numeric_; 
      do i=1 to dim(v); 
        if v(i)=. then v(i)=0; 
      end; 
      %if &instore ne %then %do; keep col:; %end;
      %else %do; drop i &covdrop; %end;
      run; 

   /*
   /  Display parameter names for use in f= or fdata=.
   /------------------------------------------------------------------*/
   %if &score= %then %do;
     %if &listnames=Y %then %do;
       data _names; set _Parm;
          Name=cats("B_",Row);
       run;
       ods select all;
       proc print data=_names;
          id Name;
          var Estimate;
       run;
       ods exclude all;
     %end;
     %if %index(%upcase(&version),SHOWNAMES) %then %goto exit;
   %end;

   /*
   /  Create input data and list of variables in it.
   /------------------------------------------------------------------*/
   proc transpose data=_Parm(rename=(Row=_NAME_)) out=_tParm(drop=_NAME_);
      format Estimate best16.;
      var Estimate;
   run;
   proc transpose data=_tParm out=_ttParm;
   run;

   /* 
   /  Build list of Parameter names and values for NLMIXED.
   /------------------------------------------------------------------*/
   data _null_; set _ttParm;
      call symput("PInit",cats(symget("PInit"))||' '||
      'B_'||cats(_NAME_)||' '||cats(COL1));
   run;

   /* 
   /  Check compatibility and prepare covariance matrix.
   /------------------------------------------------------------------*/
   data _null_;
     set _cov nobs=nobs;
     array cols (*) _numeric_;
     if _n_=1 then 
     if nobs ne &dim or nobs ne dim(cols) then do;
       _cols=dim(cols);
       put "ERROR: Covariance matrix improper or incompatible with parameter
 vector.";
       put "ERROR- The number of parameters, &dim, does not equal the number
 of covariance matrix";
       put "ERROR- rows, " nobs ", and columns, " _cols ".";
       %if &inest ne and &incovb ne %then %do;
         put "ERROR- Specify any numeric variables not containing covariance
 matrix columns in COVDROP=.";
       %end;
       call symput("status","bad_cov");
     end;
   run;
   %if &status=bad_cov %then %goto exit;
   proc transpose data=_cov out=_tcov(keep=Col:); 
      var _numeric_;
   run;

   /* 
   /  Invert covariance matrix.
   /------------------------------------------------------------------*/
   data _addI;
     array idmx{&dim} I1-I&dim;
     set _tcov;
     do _j=1 to &dim;
        if _j=_N_ then idmx[_j]=1;
        else idmx[_j]=0;
     end;
     drop _j;
   run;
   proc reg data=_addI outest=_CovInv(keep=Col:) noprint;
     model I1-I&dim = Col: / noint;
   run;quit;
   
   /* 
   /  Gather covariance values.
   /------------------------------------------------------------------*/
   data _likepieces;
     set _CovInv;
     array c(&dim) Col1-Col&dim;
     do j=1 to _n_;
       covlow=c(j);
       lhs=cats("p",_n_);
       rhs=cats("p",j);
       if j = _n_ then factor = 1;
       else factor = 2;
       keep covlow lhs rhs factor;
       output;
      end;
   run;

   %if %index(%upcase(&version),DEBUG) %then %do;
      proc print data = _likepieces;
      run;
   %end;

   data _likepieces;
      set _likepieces end = last;
   
      if _N_ = 1 then call symputx('midx',1);
      else do;
         if mod(_N_-1, &nrecdiv) = 0 then do;
            call symputx('midx', symgetN('midx')+1);
            end;
         end;
      
      if(_N_ = 1 and last = 1) then do;
         call symput('Expnt'||cats(symgetN('midx')), 
           cats(symget('Expnt'||cats(symgetN('midx'))))||
           cats(factor)||'*'||'('||trim(translate(lhs,"_","*"))||'-B_'||
           trim(translate(lhs,"_","*"))||')'||'*'||cats(covlow)||'*'||
           '('||trim(translate(rhs,"_","*"))||'-B_'||
           trim(translate(rhs,"_","*"))||')');
         end;
      else do;
         if (mod(_N_-1, &nrecdiv) = 0) then do;
            call symput('Expnt'||cats(symgetN('midx')),
              cats(factor)||'*'||'('||trim(translate(lhs,"_","*"))||'-B_'||
              trim(translate(lhs,"_","*"))||')'||'*'||cats(covlow)||'*'||
              '('||trim(translate(rhs,"_","*"))||'-B_'||
              trim(translate(rhs,"_","*"))||')'||'+');         
            end;
         
         else if (mod(_N_, &nrecdiv) = 0) then do;
            call symput('Expnt'||cats(symgetN('midx')), 
              cats(symget('Expnt'||cats(symgetN('midx'))))||
              cats(factor)||'*'||'('||trim(translate(lhs,"_","*"))||'-B_'||
              trim(translate(lhs,"_","*"))||')'||'*'||cats(covlow)||'*'||
              '('||trim(translate(rhs,"_","*"))||'-B_'||
              trim(translate(rhs,"_","*"))||')');
            end;
         else if last then do;
            call symput('Expnt'||cats(symgetN('midx')), 
              cats(symget('Expnt'||cats(symgetN('midx'))))||
              cats(factor)||'*'||'('||trim(translate(lhs,"_","*"))||'-B_'||
              trim(translate(lhs,"_","*"))||')'||'*'||cats(covlow)||'*'||
              '('||trim(translate(rhs,"_","*"))||'-B_'||
              trim(translate(rhs,"_","*"))||')');
            end;   
         else do;
            call symput('Expnt'||cats(symgetN('midx')), 
              cats(symget('Expnt'||cats(symgetN('midx'))))||
              cats(factor)||'*'||'('||trim(translate(lhs,"_","*"))||'-B_'||
              trim(translate(lhs,"_","*"))||')'||'*'||cats(covlow)||'*'||
              '('||trim(translate(rhs,"_","*"))||'-B_'||
              trim(translate(rhs,"_","*"))||')'||'+');
            end;
         end;
      run;

   /* 
   /  Use NLMIXED to compute estimate(s).
   /------------------------------------------------------------------*/
   %macro sum;
     %do i = 1 %to &midx-1;
        &&Expnt&i +
     %end;
     &&Expnt&midx;
   %mend;                  

   %if &fdata ne %then %do;
     data _null_; 
        set &fdata end=last nobs=nobs;
        if index(label,'"')>0 or index(f,'"')>0 then 
           call symput("status","quotes");
        call symput(cats("l",_n_),trim(label));
        call symput(cats("f",_n_),trim(f));
        if last then call symput("nf",nobs);
        run;
     %if &status=quotes %then %do;
   %put ERROR: Do not use double quotes (%str(%")) in the F or LABEL variable.;
   %put ERROR- Use single quotes (%str(%')).;
        %goto exit;
     %end;
   %end;
   
   %if &score ne %then %do;
     data _tParm;
       set _tParm &score(in=inscore);
       _inscore=0;
       if inscore then _inscore=1;
   %end;

   %if (^%length(&label)) %then %let label = &f;
   %let dfopt=;
   %if ( %length(&df)   ) %then %let dfopt = df = &df;

   proc nlmixed data = _tParm;
    %if &score= %then 
      ods output AdditionalEstimates=Est;
    ;
      parms &PInit;
      _e =  %sum;
      _l = -0.5*_e;
      %if &score ne %then if _inscore=0 then;
      _dummy = 1;  
      model _dummy ~ general(_l);
    %if %quote(&f) ne %then %do;
      %if &score ne %then %do;
          predict &f out=&outscore(where=(_inscore=1)) &dfopt
            %if &alpha ne %then alpha=&alpha;
          ;
      %end;
      %else %do;
        estimate "&label" &f &dfopt
          %if &alpha ne %then alpha=&alpha;
        ;
      %end;
    %end;
    
    %if &fdata ne %then %do i=1 %to &nf;
      estimate "&&l&i" &&f&i &dfopt
        %if &alpha ne %then alpha=&alpha;
      ;    
    %end;

    run;
   %if &syserr=3000 %then %do;
      %put ERROR: PROC NLMIXED in SAS/STAT is required.;
      %goto exit;
   %end;

   %if &score ne %then %do;
      data &outscore;
        set &outscore nobs=nobs;
        call symput("nobs",cats(nobs));
        %if &null ne 0 %then %do;
          Null=&null; label null="Null Value";
        %end;
        %if &df = %then %do;
           Lower = Pred - probit(1-&alpha/2)*StdErrPred;
           Upper = Pred + probit(1-&alpha/2)*StdErrPred;
           ChiSq = (tValue-&null/StdErrPred)**2; Pr = 1-probchi(ChiSq,1);
           format Pr pvalue6.4;
           label ChiSq = "Wald Chi-Square" Pr = "Pr > ChiSq";
           drop df tValue Probt;
        %end;
        %else %do;
           tValue = tValue-&null/StdErrPred; 
           Probt = 2*(1-probt(abs(tValue),&df));        
        %end;
        drop _inscore;
        run;
      %if %sysfunc(exist(&outscore)) %then %do;
        options notes;
        %put NOTE: The data set %upcase(&outscore) has &nobs observations.;
        %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;;
      %end;
   %end;
   %else %do;
      ods select all;
      %if %quote(&title) ne %then %do; 
        title "%quote(&title)";
      %end;
      %else %do;
        title "Nonlinear Function Estimate";
      %end;
      %if &df = %then %do;
         data Est;
           set Est;
           %if &null ne 0 %then %do;
             Null=&null; label null="Null Value";
           %end;
           Lower = Estimate - probit(1-&alpha/2)*StandardError;
           Upper = Estimate + probit(1-&alpha/2)*StandardError;
           ChiSq = (tValue-&null/StandardError)**2; Pr = 1-probchi(ChiSq,1);
           format Pr pvalue6.4;
           label ChiSq = "Wald Chi-Square" Pr = "Pr > ChiSq";
           drop df tValue Probt;
           run;
      %end;
      %else %do;
         data Est;
           set Est;
           %if &null ne 0 %then %do;
             Null=&null; label null="Null Value";
           %end;
           tValue = tValue-&null/StandardError; 
           Probt = 2*(1-probt(abs(tValue),&df));
           run;
      %end;
      %if &syserr=0 and %sysfunc(exist(est)) %then %do;
        options notes;
        %put NOTE: The data set EST has been created.;
        %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;;
        %if %upcase(&print)=YES %then %do;
           proc print data=Est label;
             id label;
             var Estimate StandardError 
                 %if &null ne 0 %then Null; 
                 %if &df= %then ChiSq Pr;
                 %else df tValue Probt;
                 Alpha Lower Upper;
             run;
        %end;
      %end;
      ods exclude all;
   %end;

%exit:
   /* 
   /  Delete temporary data sets; turn output and notes back on.
   /------------------------------------------------------------------*/
   %if %index(%upcase(&version),DEBUG)=0 %then %do;
     %if &status=ok %then %do;
       proc datasets nolist; 
         delete _Parm _names _Cov  
         %if %index(%upcase(&version),SHOWNAMES)=0 %then
                _tParm _ttParm _tcov _addI _CovInv _likepieces;
         ;
         run; quit;
     %end;
     %else %if &status=bad_cov %then %do;
       proc datasets nolist; 
         delete _Parm _names _Cov _tParm _ttParm;
         run; quit;
     %end;
     %else %if &status=quotes %then %do;
       proc datasets nolist; 
         delete _Parm _names _Cov _tParm _ttParm _tcov _addI 
                _CovInv _likepieces;
         run; quit;
     %end;
   %end;
   %else %do;
     %put status = &status;
     %put midx   = &midx;
     %put Expnt1 = &Expnt1;   
     %put PInit  = &PInit;
     %put dim    = &dim;
   %end;   

   %if %index(%upcase(&version),DEBUG) %then %do;
    options nomprint nomlogic nosymbolgen;
    %put _user_;
   %end;
   ods select all;
   options &notesopt;
   title;
   %let time = %sysfunc(round(%sysevalf(%sysfunc(datetime()) - &time), 0.01));
   %put NOTE: The NLEST macro used &time seconds.;
%mend NLEstimate;

%NLEstimate(&version, instore=&instore, inest=&inest, incovb=&incovb, 
            where=&where, null=&null, f=&f, fdata=&fdata, label=&label,
            score=&score, outscore=&outscore, df=&df, alpha=&alpha, 
            listnames=&listnames, title=&title, print=&print, covdrop=&covdrop);
%mend NLEST;