1
$\begingroup$

In order to run Lasso and elastic net multiple regressions on my company's SAS server (which doesn't support R), I've been working on a coordinate descent macro for performing least squares regressions (as described in the 2010 paper "Regularization Paths for Generalized Linear Models via Coordinate Descent" by Jerome Friedman, Trevor Hastie, and Rob Tibshirani).

Ideally, I would like the coefficient estimates from my SAS algorithm to match the outcomes from the glmnet package in R written by Friedman, et al. which also implements coordinate descent for least squares regression.

I've decided to test the algorithm on the Fitness data from SAS documentation, with Oxygen as the response variable:

data fitness;
  input Age Weight Oxygen RunTime RestPulse RunPulse MaxPulse @@;
  datalines;
   44 89.47 44.609 11.37 62 178 182   40 75.07 45.313 10.07 62 185 185
   44 85.84 54.297  8.65 45 156 168   42 68.15 59.571  8.17 40 166 172
   38 89.02 49.874  9.22 55 178 180   47 77.45 44.811 11.63 58 176 176
   40 75.98 45.681 11.95 70 176 180   43 81.19 49.091 10.85 64 162 170
   44 81.42 39.442 13.08 63 174 176   38 81.87 60.055  8.63 48 170 186
   44 73.03 50.541 10.13 45 168 168   45 87.66 37.388 14.03 56 186 192
   45 66.45 44.754 11.12 51 176 176   47 79.15 47.273 10.60 47 162 164
   54 83.12 51.855 10.33 50 166 170   49 81.42 49.156  8.95 44 180 185
   51 69.63 40.836 10.95 57 168 172   51 77.91 46.672 10.00 48 162 168
   48 91.63 46.774 10.25 48 162 164   49 73.37 50.388 10.08 67 168 168
   57 73.37 39.407 12.63 58 174 176   54 79.38 46.080 11.17 62 156 165
   52 76.32 45.441  9.63 48 164 166   50 70.87 54.625  8.92 48 146 155
   51 67.25 45.118 11.08 48 172 172   54 91.63 39.203 12.88 44 168 172
   51 73.71 45.790 10.47 59 186 188   57 59.08 50.545  9.93 49 148 155
   49 76.32 48.673  9.40 56 186 188   48 61.24 47.920 11.50 52 170 176
   52 82.78 47.467 10.50 53 170 172
   ;
run;

Here's my first attempt at writing the code for a simple OLS model. (I realize running a data set inside a macro loop is bad form & slows down execution times - this is just a first pass at the problem.)

For the example here I'm fitting a model for a single value of lambda and alpha in an elastic net model. I'm achieving the closest match to glmnet output when I standardize the six predictor variables using proc standard. Initial values for the coefficients are fit via proc reg. Output coefficient values are then converted back to the original unstandardized scale (scroll to bottom of the code below).

            /* Calculate mean and stnd dev values for standardizing fitness variables. */
            proc means data=fitness mean std;
               var Oxygen Age Weight RunTime RestPulse RunPulse MaxPulse;
               output out=fitness_mean_std;
            run;

            data fitness_mean_std (drop=_TYPE_ _FREQ_);
            set fitness_mean_std;
               if _STAT_ in ('MEAN','STD');
            run;

            %let t=7;
            data _null_;
            set fitness_mean_std;
               if _STAT_='MEAN' then do;
                 array mean[1:&t] Oxygen Age Weight RunTime RunPulse RestPulse MaxPulse;
                 do m = 1 to &t;
                    call symputx(cats('mean',m),mean[m],'g');
                 end;
               end;
               else if _STAT_='STD' then do;
                 array std[1:&t] Oxygen Age Weight RunTime RunPulse RestPulse MaxPulse;
                 do s = 1 to &t;
                     call symputx(cats('std',s),mean[s],'g');
                 end;
               end;
            run;

            /* Create input dataset for coordinate descent macro. */
            proc standard data=fitness mean=0 std=1 out=fitness_stnd;
               var Age Weight RunTime RunPulse RestPulse MaxPulse;
            run;

            proc reg data=fitness_stnd outest=params_stnd;
               model Oxygen = Age Weight RunTime RunPulse RestPulse MaxPulse;
            run;
            quit;

            %let t=6;
            data _null_;
            set params_stnd;
               array x[0:&t] Intercept Age Weight RunTime RunPulse RestPulse MaxPulse;
               do _n_ = 0 to &t;
                  call symputx(cats('p',_n_),x[_n_],'g');
               end;
            run;
            %put &p0 &p1 &p2 &p3 &p4 &p5 &p6;

            %macro assignvar(k);
            data fitness_array (drop=Oxygen Age Weight RunTime RunPulse RestPulse MaxPulse);
            set fitness_stnd;
               y=Oxygen;
               array a[6] Age Weight RunTime RunPulse RestPulse MaxPulse;
               array x[6];
               %do i=1 %to 6;
                 x[&i]=a[&i];
               %end;
            run;
            %mend;
            %assignvar(6)    

            /* Coordinate descent macro. */
            %macro test(dataset=, numvars=, numiter=, lambda=, alpha=);
               %do i=1 %to &numiter;
                 %do j=1 %to &numvars;
                    data &dataset (keep=y x1-x&numvars);
                    set &dataset end=end_data;
                       array x[&numvars] x1-x&numvars;
                       %let gamma = %sysevalf(&lambda*&alpha);

                       /* Calculate partial residuals for fitting coefficients.*/
                       yhat_&j = &p0 - &&p&j*x[&j];
                       %do k=1 %to &numvars;
                            yhat_&j = yhat_&j + &&p&k*x[&k];
                       %end;
                       if _n_=1 then z_&j = x&j*(y - yhat_&j);                                           else z_&j = x&j*(y - yhat_&j) + z_&j; end;
                       if end_data then do;
                          z_avg_&j = z_&j/_n_;
                          if (z_avg_&j > 0 and &gamma < abs(z_avg_&j)) then do;
                             p&j = (z_avg_&j - &gamma)/(1 + &lambda - &gamma);
                             call symputx("p&j", p&j, 'g');
                          end;  
                          else if (z_avg_&j < 0 and &gamma < abs(z_avg_&j)) then do;
                             p&j = (z_avg_&j + &gamma)/(1 + &lambda - &gamma);
                             call symputx("p&j", p&j, 'g');
                          end;
                          else if &gamma >= abs(z_avg_&j) then do;
                             p&j = 0;
                             call symputx("p&j", p&j, 'g');
                          end;
                       end;
                       retain z_&j;
                    run;
                 %end;
               %end;
               %put _user_;
            %mend;
            %test(dataset=fitness_array, numvars=6, numiter=50, lambda=.1, alpha=.5)            

            /* Return regression coefficients in original scale. */
                %let p0_unstand = %sysevalf(&p0-(&p1*&mean2/&std2)-(&p2*&mean3/&std3)-(&p3*&mean4/&std4)-(&p4*&mean5/&std5)-(&p5*&mean6/&std6)-(&p6*&mean7/&std7));
                %let p1_unstand = %sysevalf(&p1/&std2);
                %let p2_unstand = %sysevalf(&p2/&std3);
                %let p3_unstand = %sysevalf(&p3/&std4);
                %let p4_unstand = %sysevalf(&p4/&std5);
                %let p5_unstand = %sysevalf(&p5/&std6);
                %let p6_unstand = %sysevalf(&p6/&std7);

                %put 
                p0_unstand = &p0_unstand 
                p1_unstand = &p1_unstand 
                p2_unstand = &p2_unstand 
                p3_unstand = &p3_unstand 
                p4_unstand = &p4_unstand 
                p5_unstand = &p5_unstand 
                p6_unstand = &p6_unstand;
                /*
                p0_unstand = 107.068671308076 
                p1_unstand = -0.24178321947146 
                p2_unstand = -0.05008520720235
                p3_unstand = -2.47736090772018 
                p4_unstand = -0.16124847253703 
                p5_unstand = -0.03822018686055
                p6_unstand = 0.06524084784678
                */

The corresponding commands in R:

> elasticnet_fit = glmnet(x, y, family="gaussian", lambda=.1, alpha=.5); 
> coef(elasticnet_fit);   
7 x 1 sparse Matrix of class "dgCMatrix"
                      s0
(Intercept) 105.10824556
x1           -0.22996264
x2           -0.05775625
x3           -2.64766834
x4           -0.01998125
x5           -0.26222028
x6            0.18003526

My question: The coefficients output from the SAS macro doesn't match the output from glmnet, although the values are close. Is there a flaw in my code or should I not be too concerned? Thanks!

$\endgroup$
  • $\begingroup$ You should not be concerned. It is common to see a small difference even between 2 SAS procedures or using different but equivalent syntax within the same SAS procedure. They may be using a different weight, or different optimization procedure. $\endgroup$ – Peter Oct 8 '14 at 12:02
  • 1
    $\begingroup$ How about checking if the KKT conditions are fulfilled? stat.cmu.edu/~ryantibs/papers/lassounique.pdf describes them for the lasso--but basically the $X^T(Y-X\beta)$ needs to lie in the subgradient of the penalty: either exactly equal to $\lambda \mbox{sign} \beta_j + 2\lambda \alpha \beta_j$ or in an interval when $\beta_j = 0$. $\endgroup$ – Andrew M Oct 10 '14 at 2:00
  • $\begingroup$ @AndrewM - Thanks for the link - should I worry about fulfilling the KKT conditions if p < n? $\endgroup$ – RobertF Oct 10 '14 at 19:43
  • 1
    $\begingroup$ Yes, the KKT conditions tell you that you've actually found the unique optimum, so they hold regardless of $p<n$ or not. They are the analog to testing if your putative solution is at critical point of the function, when you've got a non-differentiable, but convex function. $\endgroup$ – Andrew M Oct 10 '14 at 23:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.