boxcox14a.sas

 /*------------------------------------------------------------------*
  *    Name: boxcox.sas                                              *
  *   Title: Power transformations by Box-Cox method                 *
  *          with graphic display of maximum likelihood solution,    *
  *               t-values for model effects, and influence of       *
  *               observations on choice of power.                   *
  *     Doc: http://www.math.yorku.ca/SCS/sasmac/boxcox.html         *
  *------------------------------------------------------------------*
  *  Author:  Michael Friendly            <friendly@YorkU.ca>        *
  * Created:  23 Oct 1991 10:20:14                                   *
  * Revised:  21 Aug 2002 11:46:59                                   *
  * Version:  1.4                                                    *
  *  1.2  Added ADD= parameter (additive constant); Allow no MODEL=  *
  *  -handle special missing values (.A-.Z). Fixed V8 lc varname bug *
  *  1.3  Added folded-power transformations for variables bounded   *
  *   on [0, &FOLD], specified by FOLD=                              *
  *  1.4  Added plot of logL vs. lambda                              *
  *       Changed default to PPLOT=NONE. Fixed %do-%end buglet       *
  *                                                                  *
  *------------------------------------------------------------------*/

 /*=
=Description:
 
 The BOXCOX macro Computes the maximum likelihood power transformation 
 (or folded-power transformation) for a linear regression with zero or 
 more predictors.

=Usage:

 The BOXCOX macro is defined with keyword parameters.  The RESP=
 parameter is required.
 The arguments may be listed within parentheses in any order, separated
 by commas. For example: 
 
	%boxcox(data=duncan, resp=prestige, model=income educ, id=job);
 
==Parameters:

* RESP=       Name of the response variable for analysis

* MODEL=      The independent variables in the regression, i.e., the 
              terms on the right side of the = sign in the MODEL statement 
				  for PROC REG.  May be empty.

* ID=         Name of an ID variable for observations used as labels in
              the INFL plot.

* DATA=       The name of the data set holding the response and predictor
              variables.  [Default: DATA=_LAST_]

* OUT=        Output dataset with transformed resp. Contains all original
              variables, with the transformed response replacing the 
				  original variable. [Default: OUT=_DATA_]

* OUTPLOT=    Output dataset for plot of powers Contains _RMSE_, and 
              t-values for each effect in the model, with one observation
				  for each power value tried. [Default: OUTPLOT=_PLOT_]

* PPLOT=      Printer plots: One or more of RMSE, LOGL, EFFECT, and INFL
              [Default: PPLOT=NONE]

* GPLOT=      Graphic plots: One or more of RMSE, LOGL, EFFECT, and INFL 
              [Default: GPLOT=RMSE EFFECT INFL]

* ADD=        Additive constant for response

* FOLD=       Upper bound on response (FOLD>0)

* LOPOWER=    Low value for power [Default: LOPOWER=-2]

* HIPOWER=    High value for power [Default: HIPOWER=2]

* NPOWER=     Number of power values in interval [Default: NPOWER=21]

* CONF=       Confidence coefficient of CI on power [Default: CONF=0.95]

 
==Credits:

 This program incorporates portions of program code from the
 macro ADXTRANS in the ADX system for experimental design
 distributed as part of the SAS/QC product.

=References:

  Box, G.E.P. and Cox, D.R. (1964), "An Analysis of Transformations,"
  Journal of the Royal Statistics Society, B-26, 211 -252.

 =*/


%macro boxcox(
    resp=,                   /* name of response variable            */
    model=,                  /* independent variables in regression  */
    id=,                     /* ID variable for observations         */
    data=_last_,             /* input dataset                        */
    out=_data_,              /* output dataset with transformed resp */
    outplot=_plot_,          /* output dataset for plot of powers    */
    pplot=NONE,              /* printer plots: RMSE, LOGL, EFFECT, INFL*/
    gplot=RMSE EFFECT INFL,  /* graphic plots: RMSE, LOGL, EFFECT, INFL*/
	 add=0,                   /* additive constant for response       */
	 fold=0,                  /* upper bound on response (FOLD>0)     */
    lopower=-2,              /* low value for power                  */
    hipower=2,               /* high value for power                 */
    npower=21,               /* number of power values in interval   */
    conf=.95);               /* confidence coefficient of CI on power*/

	%*-- Reset required global options;
	%if &sysver >= 7 %then %do;
		%local o1 o2;
		%let o1 = %sysfunc(getoption(notes));
		%let o2 = %sysfunc(getoption(validvarname,keyword));
		options nonotes validvarname=V6;
		%end;
	%else %do;
	   options nonotes;
		%end;

   %if (%length(&resp) = 0) %then %do;
       %put ERROR: BOXCOX: RESP= must be specified;
       %goto exit;
       %end;
 
   %let nmod = %numwords(&model);
	%if %upcase(&data)=_LAST_ %then %let data = &syslast; 

 /*
 /  Get the number of non-missing observations; quick exit if none.
 /  Check for valid data (positive).
 /-------------------------------------------------------------------*/
   %let FLAG = 0;
   data _nomiss_; 
		set &data end=eof;
		drop nwarn;
		%if (&add>0) %then %do;
			&resp = &resp + &add;
			%end;
      if (.Z < &resp <= 0) then do;
			nwarn+1;
         if nwarn<11 then 
			put "WARNING: Non-positive value of dependent variable &resp"
             "(obs=" _n_ ") " &resp ;
	      call symput('FLAG','1');
         end;
		if eof then do;
         if nwarn then
				put 'WARNING:' nwarn " non-positive values of dependent variable &resp";
			end;
      if &resp > .Z;
   run;

   proc contents noprint data=_nomiss_ out=_NOBS_;
   data _null_; set _NOBS_;
      call symput('NOBS',left(put(nobs,best20.)));
   run;
   %if (&nobs <= 0) %then %do;
		%put ERROR: There are no non-missing values of the response &resp;
		%goto done;
		%end;
 
   %if (&FLAG) %then %do;
        %put BOXCOX: Cannot compute Box-Cox transformation for non-positive &resp. Use ADD=;
        %goto done;
        %end;

	%let bcpower=Box-Cox power;
	%if &fold>0 %then %let bcpower=Folded power [0,&fold];
   %put BOXCOX: Computing &bcpower transforms of &resp ...;
 
 /*
 /  Transpose the data set; transform all values for each power into
 /  variables NEW1, NEW2, ...; and then transpose back again.  Values
 /  are centered and scaled so as to be approximately in the same scale
 /  as originally and so that the transformed value of the geometric
 /  mean is the geometric mean.
 /-------------------------------------------------------------------*/
   data _tmp_;
      set _nomiss_; keep &resp;
   proc transpose data=_tmp_ out=_tmp_;
      var &resp;
   data _tmp_; set _tmp_;
      array col{&NOBS};
      array new{&NOBS};
      keep _name_ new1-new&NOBS;
 
      gmean = 0;
      do i = 1 to &NOBS; 
			%if &fold=0 %then %do;
				gmean = gmean + log(col{i}); 
				%end;
				%else %do;
*				col{j} = col{j} / &fold;
				gmean = gmean + log(col{i} * (&fold-col{i})); 
				%end;
				end;
      gmean = exp(gmean / &NOBS);
 
      inc = (&hipower-&lopower)/(&npower-1);
      do i = 0 to &npower-1;
         lambda = &lopower + (i * inc);

         if (abs(lambda) < (inc/2)) then do;
            z1 = log(gmean);
				%if &fold=0 %then %do;
            do j = 1 to &NOBS;
               new{j} = gmean + (log(col{j})-z1)*gmean;
               end;
				%end;
				%else %do;
            do j = 1 to &NOBS;
               new{j} = (log(col{j} / (&fold-col{j})))*gmean / &fold;
               end;
				%end;
            end;
				
         else do;
            z1 = gmean**lambda;
				%if &fold=0 %then %do;
            z2 = lambda*(gmean**(lambda-1));
            do j = 1 to &NOBS;
               new{j} = gmean + ((col{j}**lambda)-z1)/z2;
               end;
				%end;
				%else %do;
				z2 = 0;
				do j = 1 to &NOBS;
					z2 = z2 + log( col{j}**(lambda-1) + (&fold-col{j})**(lambda-1));
					end;
				z2 = lambda * exp(z2/ &NOBS);
            do j = 1 to &NOBS;
               new{j} = gmean + ((col{j}**lambda - (&fold-col{j})**lambda)-z1)/z2;
               end;
				%end;
            end;

         _name_ = "_tf" || left(put(i+1,best20.));
         output;
         end;

   proc transpose data=_tmp_ out=_tmp_;
      var new1-new&NOBS;
   data _tmp_;
      merge _nomiss_ _tmp_; drop _NAME_;
   run;
 
   %put BOXCOX: Regressing on transforms for &resp ...;
 /*
 /  Perform the regression on all transformed variates at once.
 /-------------------------------------------------------------------*/
   proc reg data=_tmp_ outest=_reg_ noprint covout;
      model _tf1-_tf&npower = &model;
   run;
 
   %put BOXCOX: Computing optimal transformation and confidence interval;
 
 /*
 /  Divide estimates by their standard errors and get rid of the
 /  covariance matrix.  If no degrees of freedom left for error, exit.
 /-------------------------------------------------------------------*/
   %let nmodp1 = %eval(&nmod+1);
   %let FLAG = 0;
   data _reg_;
      drop _i_ _j_ t1-t&nmodp1;
      drop _model_ _name_ _type_ _tf1-_tf&npower _depvar_;
      array e{&nmodp1} intercep &model;
      array t{&nmodp1};
      do _i_ = 1 to &npower;
         pnt = (_i_-1)*(1+&nmodp1) + 1;
         set _reg_ point = pnt;
         adxlam = &lopower+((_i_-1)*((&hipower-&lopower)/
                  (&npower-1)));
         if (_rmse_ ^= .) then do;
            do _j_=1 to &nmodp1; t{_j_} = e{_j_}; end;
            do _j_=1 to &nmodp1;
               pnt = (_i_-1)*(2+&nmod) + 1 + _j_;
               set _reg_ point = pnt;
               t{_j_} = t{_j_} / sqrt(e{_j_});
               end;
            do _j_=1 to &nmodp1; e{_j_} = t{_j_}; end;
            _like_ = -&NOBS*log(_rmse_);
            output;
            end;
         else do;
            call symput('FLAG','1');
            end;
         end;
      stop;
   run;
   %if (&FLAG) %then %do;
      %put BOXCOX: No degrees of freedom left to estimate error.;
      %put %str(         )Transformation cannot be estimated.;
      %goto done;
      %end;
 
 /*
 /  Compute the optimal transform, the one with minimum RMSE, and an
 /  approximate confidence interval.  The approximate .95 confidence
 /  interval is based on the fact that
 /           2{ L(lambda-hat) - L(lambda) } <= cinv(.95,1) = 3.84
 /-------------------------------------------------------------------*/
   proc transpose data=_reg_ out=_reg_; *proc print;
   data _reg_;
      set _reg_; keep _name_ _label_ col1-col&npower;
      array col{&npower};
      array rmse{&npower};
      retain imin rmse1-rmse&npower;
      if (trim(_name_) = '_RMSE_') then do;
         rmse{1} = col{1};
         minrmse = rmse{1};
         imin   = 1;
         do i = 2 to &npower;
            rmse{i} = col{i};
            if (rmse{i} < minrmse) then do;
               minrmse = rmse{i};
               imin    = i;
               end;
            end;
         call symput('_imin_',left(put(imin,best20.)));
         lambda = &lopower+((imin-1)*(&hipower-&lopower)/(&npower-1));
         call symput('elambda',left(put(lambda,best20.)));
         end;
      else if (trim(_name_) = '_LIKE_') then do;
         call symput('maxlike',left(put(col{imin},best20.)));
         end;

   proc transpose data=_reg_ out=_reg_;
*--  Added code for plot of logL vs. lambda   8/02/02 MF;
   data &outplot;
      set _reg_ end=eof;
      drop _name_  adxlam _lolam_ _hilam_ _hirmse_ tval; ** _label_ removed;
		drop _cutoff_;
      retain _lolam_ 10 _hilam_ -10 _hirmse_ 0;
      _lambda_ = adxlam;
      label _lambda_ = 'Box-Cox Power (lambda)'
            _like_   = 'Log Likelihood'
            conf     = "&conf Confidence Interval";
	   _cutoff_ = &maxlike - (cinv(&conf,1)/2);
      if (_like_ < _cutoff_ )  /* was 1.92 */
         then conf = " ";
         else do;
              conf = "*";
              _lolam_ = min(_lolam_,_lambda_);
              _hilam_ = max(_hilam_,_lambda_);
              _hirmse_= max(_hirmse_,_rmse_);
              end;
      if eof then do;
         call symput('lline',left(put(_cutoff_,best20.)));
         call symput('lolam',left(put(_lolam_,best20.)));
         call symput('hilam',left(put(_hilam_,best20.)));
         call symput('hirms',left(put(_hirmse_,best20.)));
         tval = tinv(1-(1-&conf)/2, &nobs-&nmodp1 );
         call symput('tval', left(put(tval,best10.)));
         end;
   run;
   options notes;
   proc print data=&outplot label;
      var _lambda_ _like_ _rmse_ conf;
   run;
   %put BOXCOX: Estimated Power Transformation, Lambda = &elambda.;
 
 /*
 /  Plot likelihood and t-values for effects as functions of the
 /  power.
 /-------------------------------------------------------------------*/
   data &out; set _tmp_; drop _tf1-_tf&npower;
      &resp = _tf&_imin_;
      label &resp = "Transformed &resp (lambda=&elambda)";
 
   %let pplot = %upcase(&pplot);
   %if &pplot ^= NONE %then %do;
   %let sym = 1 2 3 4 5 6 7 8 9;
   %let sym = &sym A B C D E F G H I J K L M N O P Q R S T U V W X Y Z;
   %let sym = &sym a b c d e f g h i j k l m n o p q r s t u v w x y z;
   proc plot data=&outplot;
      %if %index(&pplot,RMSE)>0 %then %do;
         title 'RMSE for Box-Cox Power Transform';
         plot _rmse_ * _lambda_ = 'L'
            / href=&lolam &hilam vref=&hirms;
         label _rmse_ = "Root Mean Squared Error for &resp";
      run;
		%end;
    %if %index(&pplot,LOGL)>0 %then %do;
	 plot _like_ * _lambda_ = 'L' /
            / href=&lolam &hilam ;
      title "Box-Cox Power Transform for &resp";
      label _like_ = "Log Likelihood";
      %end;
      %if (&nmod>0) and %index(&pplot,EFFECT)>0 %then %do;
         plot
         %do i = 1 %to &nmod;
            %scan(&model,&i) * _lambda_ = "%scan(&sym,&i)"
            %end;
            / overlay href=&lolam &hilam ;
         title 't-values for Model Effects';
      run;
      %end;
   %end;
 
   %let gplot = %upcase(&gplot);
   %if &gplot ^= NONE %then %do;
   %let sym = dot star square circle triangle $ + hash x;
   %if (&nmod>0) and %index(&gplot,EFFECT)>0 %then %do;
   data _labels_;
      set &outplot end=eof;
      length function text $8;
      if eof then do;
         xsys='2'; ysys='2'; size=1.1;
         function='LABEL';  position='6';
         x = _lambda_;
         %do i = 1 %to &nmod;
            y = %scan(&model,&i) ;
            text = " %scan(&model,&i) "; output;
            %end;
         end;
   %end;
 
proc gplot data=&outplot;
	symbol1 i=spline v=dot c=black;   *-- make sure symbol1 is defined;
   %do i = 1 %to &nmod;
      symbol&i i=spline v=%scan(&sym,&i) c=black ;
      %end;
      axis1 label=(a=90);
      axis2 label=('Box-Cox Power (' f=cgreek 'l' ')' ) 
            offset=(3);
      axis3 label=(a=90 't-value');
      axis4 label=('Box-Cox Power (' f=cgreek 'l' ')' ) 
            offset=(2,9);
   %if %index(&gplot,RMSE)>0 %then %do;
      plot _rmse_ * _lambda_ = 1 /
            href=&lolam &hilam lhref=20 chref=red
            vref=&hirms        lvref=33 cvref=red
            vaxis=axis1 haxis=axis2 hminor=1 vminor=1
				des="BoxCox plot of RMSE * Lambda for &resp";
      title h=1.5 "Box-Cox Power Transform for &resp";
      label _rmse_ = "Root Mean Squared Error";
      run;
	%end;

    %if %index(&gplot,LOGL)>0 %then %do;
	 plot _like_ * _lambda_ = 1 /
			href=&lolam &hilam lhref=20 chref=red 
			vref=&lline lvref=33 cvref=red
		   vaxis=axis1 haxis=axis2 hminor=1 vminor=1
				des="BoxCox plot of Log Likelihood * Lambda for &resp";
      title h=1.5 "Box-Cox Power Transform for &resp";
      label _like_ = "Log Likelihood";
		run;
	%end;

	   %if %index(&gplot,EFFECT)>0 or %index(&gplot,INFL)>0 %then %do;
			%gskip;
			%end;
      %end;
   %if (&nmod>0) and %index(&gplot,EFFECT)>0 %then %do;
proc gplot data=&outplot;
      plot
         %do i = 1 %to &nmod;
            %scan(&model,&i) * _lambda_ = &i
            %end;
         / overlay anno=_labels_
            href=&lolam &hilam lhref=20 chref=red
            vref=&tval -&tval  lvref=33 cvref=red
           vaxis=axis3 haxis=axis4 hminor=1
			  des="BoxCox Effect plot for &resp";
       title h=1.5 "t-values for Model Effects on &resp";
       run;
	   %if %index(&gplot,INFL)>0 %then %do;
			%gskip;
			%end;
   %end;  /* if &gplot ^= NONE */
   %if &pplot ^= NONE or &gplot ^= NONE %then %do;
   title ;
   %end;
 
%done:
   proc datasets nolist;
      delete _NOBS_  _tmp_ _nomiss_ _reg_;
      %if (&nmod>0) and %index(&gplot,EFFECT)>0 %then %do;
      delete _labels_;
      %end;
   quit;

   %if %index(&gplot,INFL)>0 or %index(&pplot,INFL)>0
      %then %do;
         %boxinfl(data=&data, resp=&resp, model=&model, id=&id, fold=&fold,
                  gplot=&gplot, pplot=&pplot);
      %end;
%exit:
	%*-- Restore global options;
	%if &sysver >= 7 %then %do;
		options &o1 &o2;
		%end;
	%else %do;
	   options notes;
		%end;
%mend;
 
 /*------------------------------------------------------------------*/
 /* NAME:    BOXINFL                                                 */
 /* PURPOSE: Computes Atkinson's score-test for power transformation */
 /*          and produces a constructed-variable influence plot for  */
 /*          the impact of observations on the choise of power.      */
 /*------------------------------------------------------------------*/
%macro boxinfl(
    resp=,                   /* name of response variable            */
    model=,                  /* independent variables in regression  */
    data=_last_,             /* input dataset                        */
    id=,                     /* ID variable for observations         */
	 fold=0,
    gplot=INFL,              /* Graphic plot?                        */
    pplot=INFL);             /* Printer plot?                        */
 
   data _cvar_;
      set &data;
		%if &fold=0 %then %do;
	      logy = log(&resp);       *-- geometric mean =mean log(y);
			%end;
		%else %do;
	      logy = log(&resp * (&fold-&resp));
			%end;
      %if %length(&id)=0 %then %do;
        %let id =_id_;
        _id_ = _n_;
      %end;
   proc means noprint;
      var logy;
      output out=_gm_ mean=gmean;
   data _cvar_;
      set _cvar_;
      drop gmean;
      if _n_=1 then do;
         set _gm_;
         gmean = exp(gmean);
         end;
		%if &fold=0 %then %do;
	      g = &resp * ( log ( &resp / gmean ) - 1);
			%end;
		%else %do;
	      g = &resp * log(&resp/gmean) - (&fold-&resp) * log(&fold-&resp/gmean) ;
			%end;
      label g='Constructed variable';
 
 /* produce values for partial-regression plot of residuals from
 /  &resp vs. residuals from constructed variable.  1-slope =
 /  power for transformation based on score test
 / -------------------------------------------------------------*/
   proc reg data=_cvar_ outest=_parm_ ;
     id &id;
    m0: model &resp=&model g;       * y = Xb + lambda g;
        output out=m0 rstudent=_resy_ cookd=_infl_;
   proc reg data=_cvar_ noprint;
     id &id;
    m1: model &resp=&model;
        output out=m1 r=_resx_;     * y - Xb;
    m2: model g =&model;
        output out=m2 r=_resg_;     * g - Xb;
   data _part_;
      keep _resx_ _resg_ &id _resy_ _infl_;
      merge m0 m1 m2;
   proc means noprint data=_part_;
      var _resg_;
      output out=_gm_ min=rg_min max=rg_max;
   data _slope_;
      set _parm_(keep=_model_ g);
      if (_model_='M0');
      xsys='1'; ysys='1';
      x=8; y=90;
      function = 'LABEL';
      position='3';  size=1.5;
      text = 'Slope:  '|| put(g,best5.);  output;
      position='F';
      text = 'Power:  '|| put(round(1-g,.25),best5.);  output;
      call symput('slope',put(g,best5.));
      call symput('power',put(round(1-g,.25),best5.));
      run;
   %if %index(&pplot,INFL)>0 %then %do;
      proc plot data=_part_;
         plot _resx_  * _resg_ = &id / vref=0;
         label _resx_ ="Partial &resp"
               _resg_ ='Partial Constructed Variable';
         title 'Partial Regression Influence plot for Box-Cox power';
         title2 "Slope: &slope   Power: &power";
      run;quit;
   %end;
 
%if %index(&gplot,INFL)>0 %then %do;
data _anno_;
   set _part_ ;
   length text $16;
   if _n_=1 then set _gm_;
   if abs(_resg_/(rg_max-rg_min)) > .5
    | abs(_resy_) > 3
    | _infl_> 4/(n-1);
   xsys='2'; ysys='2';
   x =  _resg_; y=_resx_ ;
   function='LABEL';
   if _resg_ > 0 then position='1';
                 else position='3';
   text=&id;
   data _anno_;
      set _anno_ _slope_;
   proc gplot data=_part_;
      plot _resx_  * _resg_ /
         vaxis=axis1 vminor=1 hminor=1
         vref=0 lvref=34 anno=_anno_
			name='boxinfl'
			des="BoxCox influence plot for &resp";

      axis1 label=(a=90);
      symbol i=rl h=1.3 v=circle ci=red;
      label _resx_ ="Partial &resp"
            _resg_ ='Partial Constructed Variable';
      title h=1.4 'Partial Regression Influence plot for Box-Cox power';
   run;quit;
	title;
%end;
   proc datasets nolist;
      delete _cvar_ _parm_ m0 m1 m2 _slope_;
   quit;
%mend;
 /*-------------------------------------------------------------------*/
 /* Copyright (c) 1987 by SAS Institute Inc. Cary NC 27511  USA       */
 /*                                                                   */
 /* NAME:    NUM(ber of )WORDS                                        */
 /* PURPOSE: Returns the number of words in a given list, with an     */
 /*          optional specification of word delimiters.               */
 /*-------------------------------------------------------------------*/
%macro numwords(lst,wordchar);
   %let i = 1;
   %if (%length(&wordchar)) %then %do;
      %let v = %scan(&lst,&i,&wordchar);
      %do %while (%length(&v) > 0);
         %let i = %eval(&i + 1);
         %let v = %scan(&lst,&i,&wordchar);
         %end;
      %end;
   %else %do;
      %let v = %scan(&lst,&i);
      %do %while (%length(&v) > 0);
         %let i = %eval(&i + 1);
         %let v = %scan(&lst,&i);
         %end;
      %end;
   %eval(&i - 1)
%mend;