SAS Programming for Data Mining Applications

Multi-Threaded Principle Component Analysis

2012-01-31T21:45:00.001-05:00

SAS used to not support multithreading in PCA, then I figured out that its server version supports this functionality, see here. Today, I found this mutlithreading capability is finally available in PC SAS v9.22.

The figure above indicates that all 4 threads in my PC are utilized. FYI, My PC uses an Intel 2core 4threads CPU. This multi-threading capability directly help any work relying on SVD due to the direct relationshipbetween SVD and PC, see here.

Notice that in order to observe the effect of mutli-threading by comparing Real User Time and CPU Time, I/O should not be a bottleneck, that is why in the code, all outputs, either to screen or to data sets, are suppressed.



options fullstimer;
data _junky;
     length id x: 8;
  array x{800};
  do id=1 to 5E3;
     do j=1 to dim(x);
     x[j]=ranuni(0);
  end;
  drop j; output;
  end;
run;

proc princomp data=_junky noprint;
      var x:;
run;

Random Number Seeds: NOT only the first one matters!

2012-01-30T11:32:00.001-05:00

Today, Rick (blog @ here) wrote an article about random number seed in SAS to be used in random number functions in DATA Step. Rick noticed when multiple random number functions are called using different seeds, only the first one matters.

This is so true. In fact, SAS Manual also has a comprehensive writting on this issue, namely, how to control seeds for each iteration of the random number generation process and how to generatie multiple statistically independent streams of random numbers, see here. In fact, sasCommunity.org also has an article about this issue, see here.

To echo Rick's post, there is a way to control the seed so that NOT only the first one matters: use CALL routines which by theory will generate computationally independent random number sequence. But if the two seeds are too close, the generated sequences may not be statistically independent. Again, refer to the SAS manual for details.



data normal;
   seed1 = 11111;
   seed2 = 22222;
   seed3 = 383333;
   do i = 1 to 1000;
      call rannor(seed1, x1);
      call rannor(seed2, x2);   
      call rannor(seed3, x3);
   x4=rannor(seed2);
   x5=rannor(seed3);
      output;
   end;
run;

proc sgscatter data=normal;
   matrix x1-x5/ markerattrs = (size = 1);
run;

Using PROC CANCORR to solve large scale PLS problem

2011-11-16T15:12:00.001-05:00

Partial Least Square (PLS) is a powerful tool for discriminant analysis with large number of predictors [1].

PLS extracts latent factors that maximize the covariance between independent variables and dependent variables. This process is equivalent to Generalized Eigenvalue Decomposition of the following formula [2]:
$$X'HXw =\phi X'Xw $$. For PLS $$H=Y'Y$$ Note that Canonical Correlation Analysis (CCA) follows the same generalized eigenvalue decomposition problem, specifically, for CCA, $$H=Y'(YY')^{-1}Y$$.

In SAS, PROC PLS implements 2 forms of PLS, namely the original NIPALS [2] and SIMPLS [3]. When there is only one dependent variable, the two algorithms generate the same output.PLS is a computationally very demanding algorithm. While powerful, when the dimension of the problem at hand becomes very large, PROC PLS will encounter issues such as insufficient memory and very long computing time. There is a rescue when only one dependent variable Y presents. In this case, CCA and PLS differ only up to a fixed scale parameter. Therefore we can use PROC CANCORR, which is very scalable and multithreaded, to solve the PLS problem. The obtained weights and loadings will not be the same but the difference is only up to a fixed scale parameter.

In the follow log, we demonstrate the behavior of PROC PLS and PROC CANCORR on a server with 4GB accessible memory when the number of independent variable is 5000 and sample size is 100K. PROC PLS reported insufficient memory and stopped computing in 45seconds after exhausting all accessible memory, while PROC CANCORR continued and finished computation in slightly more than 7 minutes. Both procedures used up 3.92GB memory available to SAS from the system. Also note that PROC CANCORR used more than 33minutes of CPU time, indicating its very good scale up capability in a multi-core environment.

Referece:
[1] Barker, M and Rayens, W (2003), "Partial Least Squares for Discrimination", Journal of Chemometrics, 17, 166-173

[2] Sun, L; Ji, S; Yu, S; and Ye, J (2009), "On the Equivalence Between Canonical Correlation Analysis and Orthonormalized Partial Least Squares", In Proceedings of the 21st International Joint Conference on Artificial Intelligence (IJCAI 2009).

[3] Wold, H. (1966), “Estimation of Principal Components and Related Models by Iterative Least Squares,” in P. R. Krishnaiah, ed., Multivariate Analysis, New York: Academic Press.

[4] de Jong, S. (1993), “SIMPLS: An Alternative Approach to Partial Least Squares Regression,” Chemometrics and Intelligent Laboratory Systems, 18, 251–263.


NOTE: PROCEDURE PRINTTO used (Total process time):
      real time           0.00 seconds
      user cpu time       0.00 seconds
      system cpu time     0.00 seconds
      Memory                            3920274k
      OS Memory                         3926280k
      Timestamp            11/16/2011  2:11:22 PM
      Page Faults                       0
      Page Reclaims                     9
      Page Swaps                        0
      Voluntary Context Switches        1
      Involuntary Context Switches      1
      Block Input Operations            0
      Block Output Operations           0
      

16         options fullstimer;
17         data x;
18              length id y x: 8;
19           array x{5000};
20              do id=1 to 1E5;
21              y=rannor(0);
22              do j=1 to dim(x);
23              x[j]=rannor(0);
24           end;
25           output;
26           drop j;
27           end;
28         run;

NOTE: The data set WORK.X has 100000 observations and 5002 variables.
NOTE: Compressing data set WORK.X increased size by 0.06 percent. 
      Compressed is 100066 pages; un-compressed would require 100009 pages.
NOTE: DATA statement used (Total process time):
      real time           1:07.58
      user cpu time       1:02.41
      system cpu time     5.16 seconds
      Memory                            3920274k
      OS Memory                         3926280k
      Timestamp            11/16/2011  2:12:29 PM
      Page Faults                       0
      Page Reclaims                     510
      Page Swaps                        0
      Voluntary Context Switches        63
      Involuntary Context Switches      107
      Block Input Operations            0
      Block Output Operations           0
      

29         
30         
31         proc pls data=x method=simpls noprint;
32              model y =x1-x5000;
33         run;

ERROR: The SAS System stopped processing this step because of insufficient memory.
NOTE: There were 100000 observations read from the data set WORK.X.
NOTE: PROCEDURE PLS used (Total process time):
      real time           45.89 seconds
      user cpu time       40.11 seconds
      system cpu time     5.77 seconds
      Memory                            3920284k
      OS Memory                         3926280k
      Timestamp            11/16/2011  2:13:15 PM
      Page Faults                       0
      Page Reclaims                     978075
      Page Swaps                        0
      Voluntary Context Switches        91
      Involuntary Context Switches      162
      Block Input Operations            0
      Block Output Operations           0
      
35         proc cancorr data=x noprint;
36              var y;
37           with x1-x5000;
38         run;

NOTE: PROCEDURE CANCORR used (Total process time):
      real time           7:02.85
      user cpu time       33:13.85
      system cpu time     4.40 seconds
      Memory                            3920284k
      OS Memory                         3926280k
      Timestamp            11/16/2011  2:20:18 PM
      Page Faults                       0
      Page Reclaims                     126339
      Page Swaps                        0
      Voluntary Context Switches        2096
      Involuntary Context Switches      83359
      Block Input Operations            0
      Block Output Operations           0
      

39    
40         proc printto; run;

Obtain Trace of the Projection Matrix in a Linear Regression

2011-10-06T14:21:00.006-04:00

Recently, I am working on coding in SAS for a set of regularized regressions and need to compute trace of the projection matrix:
$$ S=X(X'X + \lambda I)^{-1}X' $$.

Wikipedia has a well written introduction to Trace @ here.

To obtain the inverse of matrix (X'X + \lambda I) in SAS/STAT, there are multiple ways:
1. Build the SSCP matrix first, then inverse it using the following methods:
    1.1. Using SVD based method, shown @ here but not demonstrated in this post;
    1.2. Using PROC REG, shown @ here , and also demonstrated in the code below;
    1.3. Using SWEEP operator, shown @ here under Method 0 and Method 1;
2. Since (X'X + \lambda I) is the SSCP for a Ridge Regression with ridge parameter=\lambda, it is handy to directly use the ODS OUTPUT InvXPX= statement to obtain the inversed matrix, when X is appended by a diagonal matrix of \lambda*I_{p, p}. SAS code is demonstrated below;

Trace of the project matrix S is a key concept in modern regression analysis. For example, the effective degree of freedom of the model in a regularized linear regression is trace(S)/N, see [1] for details. For another example, in approximating leave-one-out-cross-validation using GCV, trace(S)/N is a key component (formula 7.52 of ~[1]).

Check the reference and the cited works therein for more information.

Ref:

1. Hastie et al., The Elements of Statistical Learning, 2nd Ed.


/* Use the SocEcon example data created in
   Example A.1: A TYPE=CORR Data Set Produced by PROC CORR
   on page 8153 of SAS/STAT 9.22 User Doc
*/
data SocEcon;
input Pop School Employ Services House;
datalines;
5700 12.8 2500 270 25000
1000 10.9 600 10 10000
3400 8.8 1000 10 9000
3800 13.6 1700 140 25000
4000 12.8 1600 140 25000
8200 8.3 2600 60 12000
1200 11.4 400 10 16000
9100 11.5 3300 60 14000
9900 12.5 3400 180 18000
9600 13.7 3600 390 25000
9600 9.6 3300 80 12000
9400 11.4 4000 100 13000
;
run;

%let depvar=HOUSE;
%let covars=pop school employ services;
%let lambda=1;
/* Below is the way to obtain trace(S), where S is the project matrix in a (regularized) linar regression. 
   For further information, check pp.68, pp.153 of Elements of Statistical Learning,2nd Ed.
*/

/* For details about TYPE=SSCP special SAS data, consult:
  Appendix A: Special SAS Data Sets, pp.8159 of SAS/STAT 9.22 User's Guide
*/
proc corr data=SocEcon sscp out=xtx(where=(_TYPE_='SSCP'))  noprint;
     var &covars;
run;


data xtx2;
     set xtx;
  array _n{*} _numeric_;
  array _i{*} i1-i5 (5*0);
  do j=1 to 5;
     if j=_n_ then _i[_n_]=λ
  else _i[j]=0;
  end;
  _n[_n_]=_n[_n_]+1;
  drop j _TYPE_  _NAME_;
run;

/* Obtain the inverse of (XTX+\lambda*I)
   Note that we explicitly specified Intercept term in the 
   covariate list and fit a model without implicit intercept 
   term in the model.
*/
proc reg data=xtx2  
         outest=S0(type=SSCP
                   drop=i1-i5 _MODEL_  _DEPVAR_  _RMSE_)
         singular=1E-17;
     model i1-i5 = Intercept &covars / noint   noprint;
run;quit;

data S0;
     set S0; 
  length _NAME_ $8;
  _NAME_=cats('X_', _n_);
run;

proc score data=SocEcon  score=S0  out=XS0(keep=X_:)  type=parms;
     var &covars;
run;     

data XS0X;
     merge XS0  X;
  array _x1{*} X_:;
  array _x0{*} intercept pop school employ services;
  do i=1 to dim(_X1);
     _x1[i]=_x1[i]*_x0[i];
  end;
  rowSum=sum(of _x1[*]);
  keep rowSum;
run;
proc means data=XS0X  noprint;
     var rowSum;
  output out=trace  sum(rowSum)=Trace;
run;

Verify the result using R:


> socecon<-read.csv('c:/socecon.csv', header=T)
> x<-as.matrix(cbind(1, socecon[,-5]))
> xtx<-t(x)%*%x
> phi<-xtx+diag(rep(1, 5))
> 
> # method 1. 
> S<-x%*%solve(phi)*x
> sum(S)
[1] 4.077865
> # method 2. 
> S<-(x%*%solve(phi))%*%t(x)
> sum(diag(S))
[1] 4.077865
>

Of course, except for method 1.2 shown above, we can also use Method 2 mentioned above, and obtain the same inverse matrix:



data SocEcon2 /view=SocEcon2;
     set x  end=eof;
  array x{5} Intercept &covars (5*0);
  Intercept=1;
  output;
  if eof then do;
     do j=1 to dim(x);
     x[j]=λ
     output;
     drop j;
     x[j]=0;
  end;
  end;
run;

ods select none;
ods output InvXPX=S1;
proc reg data=SocEcon2  singular=1E-17;
     model y = Intercept &covars /noint  i;
run; quit;
ods select all;

As a side point, it is also of interests to compare the computational performance of both illustrated methods. SVD-based approach and SWEEP operator based approach are omitted.

Using a SAS data set of 1001 covariates (including Intercept term) and 1E6 observations, total 7.6GB on a windows PC, the test was done on two dramatically different machines: one WindowsXP laptop with mediocre HDD and a 2-core T7300 CPU; the other one is a high-end Server running Linux64 with Disk Array and fast CPUs totaled 16 cores.

On the PC,:
-> Using method 1.2, the time used decomposition is listed below:
PROC CORR: real time: 25:47.72; CPU time: 15:26.15
DATA step on XTX: real time: 2.28 sec; CPU time: 0.07 sec
PROC REG : real time: 15.45 sec; CPU time: 27.31 sec;
Total real time: 26:05.45

-> Using method 2, the time used is:
PROC REG: real time: 48:28.44; CPU time: 1:16:46.41

On the server:
-> Using method 1.2, the time used decomposition is listed below:
PROC CORR: real time: 5:40.61; CPU time: 5:40.58
DATA step on XTX: real time: 0.05 sec; CPU time: 0.05 sec
PROC REG: real time 1.71 sec; CPU time: 5.27 sec
Total real time: 5:42.37

-> Using method 2, the time used is:
PROC REG: real time: 6:01.46; CPU time: 19.13.49

The performance is summarized below:

Benchmark Regression Procedures using OLS Regression

2011-08-18T15:18:00.007-04:00

Rick Wicklin discussed in his blog the performance in solving a linear system using SOLVE() function and INV() function from IML.

Since regression analysis is an integral part of SAS applications and there are many SAS procedures in SAS/STAT that are capable to conduct various regression analysis, it would be interesting to benchmark their relative performance using OLS regression, the fundamental regression analysis of all.

The analysis will compare REG, GLMSELECT, GENMOD, MIXED, GLIMMIX, GLM, ORTHOREG, HPMIXED and TRANSREG on 10 OLS regressions with 100 to 1000 variables, incremental at 100, and with the number of observations twice the number of variables to avoid possible numerical issues. HPMIXED uses sparse matrix techniques and will be put into great disadvantage in this comparison using large dense matrices. A macro wraps them together:



%macro wrap;
proc printto log='c:\testlog.txt';run;

%let t0=%sysfunc(datetime(), datetime.);
%let procnames=GLM REG GLMSELECT ORTHOREG MIXED GLIMMIX GENMOD ;
%let nproc=%sysfunc(countW(&procnames));
%put Test &nproc PROCEDURES;
%do i=1 %to 10;
    %let nobs=%sysevalf(&i*100);
    options nonotes;
    data _temp;
         array x{&nobs};
	 do i=1 to 2*&nobs;
	    do j=1 to &nobs;
	       x[j]=rannor(0);
            end;
	    y=rannor(0);
	    drop i j;
	    output;
	 end;		 
     run;
     options notes;
     sasfile _temp load;
     ods select none;

     %do j=1 %to &nproc;
         %let proc=%scan(&procnames, &j);
	 %put &proc;
	 proc &proc data=_temp;
	      model y = x1-x&nobs;
	 run;
    %end;
    %put TRANSREG;
    proc transreg data=_temp;
         model identity(y) = identity(x1-x&nobs);
    run;
    sasfile _temp close;
    ods select all;
%end;
proc printto; run;
%mend;
%wrap;

After running all iterations, the SAS log is parsed to obtain procedure names and corresponding real time and CPU time. The following SAS code does this job:



data proc_compare;
     infile "c:\testlog.txt";
	 input;
	 retain procedure ;
	 retain realtime cputime  realtime2 ; 
	 length procedure $12.;
	 length realtime  cputime $24.;
	 if _n_=1 then id=0;
	 x=_infile_;
	 if index(x, 'PROCEDURE')>0 then do;
	    procedure=scan(_infile_, 3);		
		if procedure="REG" then id+1;		
	 end;
	
	 if index(x, 'real time')>0 then do;
	    _t1=index(_infile_, 'real time');
	    _t2=index(_infile_, 'seconds');
	    if _t2=0 then _t2=length(_infile_);
            realtime=substr(_infile_, _t1+9, _t2-_t1-9);
	    if index(realtime, ':')>0 then do;
 	       realtime2=scan(realtime, 1, ':')*60;
	       sec=input(substr(realtime, index(realtime, ':')+1), best.);
	       realtime2=realtime2+sec;		 
	    end;
	    else realtime2=input(compress(realtime), best.);
	 end;
	 if index(x, 'cpu time')>0 then do;
	    _t1=index(_infile_, 'cpu time');
	    _t2=index(_infile_, 'seconds');
	    if _t2=0 then _t2=length(_infile_);
	    cputime=substr(_infile_, _t1+8, _t2-_t1-8);
	    if index(cputime, ':')>0 then do;
 	       cputime2=scan(cputime, 1, ':')*60;
	       sec=input(substr(cputime, index(cputime, ':')+1), best.);
	       cputime2=cputime2+sec;
	    end;
	    else cputime2=input(compress(cputime), best.);
	    keep id size  procedure cputime2 realtime2 ;
	    size=id*100;
	    if compress(procedure)^="PRINTTO" then output;
	end;
run;

We then visualize the results using the following code:



title "Benchmark Regression PROCs using OLS";
proc sgpanel data=proc_compare;
     panelby procedure /rows=2;
     series y=cputime2  x=size/ lineattrs=(thickness=2);
	 label cputime2="CPU Time (sec)"
	       size="Problem Size"
		   ;;	
	 colaxis grid;
	 rowaxis grid;
run;
title;

title "Closer Look on REG vs. GLM vs. GLMSELECT";
proc sgplot data=proc_compare  uniform=group;
     where procedure in ("GLMSELECT", "REG", "GLM");
     series x=size y=cputime2/group=procedure  curvelabel lineattrs=(thickness=2);
	 label cputime2="CPU Time (sec)"
	       size="# of Variables"
		   ;;
     yaxis grid ;
     xaxis grid ;
run;
title;

It is found that PROC GLM and GLMSELECT beat all other procedures with large margin while HPMIXED is the slowest followed by GLIMMIX. Surprisingly, REG is slower than both GLM and GLMSELECT even though it utilized multi-threading technique while GLMSELECT does not:

************ Partial LOG of the last iteration ********
NOTE: PROCEDURE REG used (Total process time):
real time 6.79 seconds
cpu time 9.36 seconds

NOTE: There were 2000 observations read from the data set WORK._TEMP.
NOTE: PROCEDURE GLMSELECT used (Total process time):
real time 3.06 seconds
cpu time 2.96 seconds
********************************************************

The performance gap between REG and GLM/GLMSELECT is getting larger when the number of variables increases to be more than 700.

Both REG and GLMSELECT are developed by the same group of developers in SAS, as far as I know.

********************* PS : ****************************
Rick and Charlie pointed out that real time is a more fair measure, which I agree.

The reading of real computing time has large variance from run to run because the testing enviornment is not very clean and there are many background window programs running. Below is part of the log file of another run with 2000 variables and 4000 records:

NOTE: PROCEDURE REG used (Total process time):
      real time           2.26 seconds
      cpu time            7.76 seconds


NOTE: PROCEDURE GLM used (Total process time):
      real time           3.57 seconds
      cpu time            4.58 seconds

NOTE: There were 2000 observations read from the data set WORK._TEMP.
NOTE: PROCEDURE GLMSELECT used (Total process time):
      real time           3.50 seconds
      cpu time            3.44 seconds

We see that REG has lower real time comparing to GLM/GLMSELECT, even though cpu time is about twice the average of GLM/GLMSELECT. In a case where BY-processing is used, GLMSELECT will use multi-threading as specified in PERFORMANCE statement, and the gap in real time between REG and GLMSELECT will be eliminated. In a collaborating environment, more CPU time also means competing for more resources. Below we show the real time of the same run as in above CPU Time figure.

Note that CPU Time difference and pattern is pretty consistent.Below is the mean CPU Time and its 90% C.I. of 100 runs using REG /GLM /GLMSELECT on different size of problems.

Rolling Window Regression of Time Series

2011-08-10T07:38:00.001-04:00

More often than not, we encounter a problem where an OLS over a rolling time window is required, see [1], [2], [3], [4], [5], [6], [7], for a few examples.

One solution is to resort to SAS MACRO, but it is extremely inefficient and can't handle large dataset in reality, [8]. This method is shown below as Method[4]. It couldn't finish the test in allowable time using the sample data below.

The other common solution is to use the BY-processing capability of PROC REG after re-shaping the data into appropriate format, see [9], [10]. This method is demonstrated as Method[3] below. While certainly much better than above one, it is still not the fastest and requires more memory.

The third solution comes to play by recognizing that in OLS, what you need is the SSCP and you can easily build up the SSCP over rolling time window by resorting to PROC EXPAND. This is demonstrated as Method[2] below. This approach will further improve the speed but still requires large amount of memory if the data is big and many rolling windows are generated.

Since what we need to do is to build the SSCP matrix and obtain the coefficient estimates based on the informaiton in SSCP, we can certainly code this in a DATA Step using ADJUST operator, which provides a solution that is both fast and low memory occupancy. See [11] for an introduction to ADJUST operator. To make this even faster, a modification of ADJUST operator, the SWEEP operator, can be used. For an introduction to SWEEP operator, see [11], [12]. In the code below, Method[0] implements the ADJUST operator, while Method[1] implements the SWEEP operator.

The experiment results are shown below:

                         Real Time    |        CPU Time          | Memory
=====================================================

Method 0 |    1.01 (seconds)   |    1.01 (seconds)        |    611K

Method 1 |   0.25 (seconds)   |    0.24 (seconds)        |    432K

Method 2 |    1.61 (seconds)   |    0.94 (seconds)        | 50381K

Method 3 | 80.54 (seconds)   |   79.61 (seconds)       | 2322K

Method 4 |         Failed           |           Failed               |    Failed
=====================================================

Reference:
[1]MYSAS.NET, http://www.mysas.net/forum/viewtopic.php?f=4&t=8070
[2]MYSAS.NET, http://www.mysas.net/forum/viewtopic.php?f=4&t=7898
[3]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0604D&L=sas-l&P=R32485
[4]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0704C&L=sas-l&P=R3305
[5]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0802C&L=sas-l&P=R9746
[6]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0801C&L=sas-l&P=R14671
[7]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0810A&L=sas-l&P=R19135
[8]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0802C&L=sas-l&P=R13489
[9]Michael D Boldin, "Programming Rolling Regressions in SAS", Proceedings of NESUG, 2007
[10]SAS-L, http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0604D&L=sas-l&D=0&P=56926
[11]J. H. Goodnight, "The Sweep Operator: Its Importance in Statistical Computing", SAS Tech Report R-106, 1978
[12]Kenneth Lange, "Numerical Analysis for Statisticians", Springer, 1998


proc datasets library=work kill; run;


options fullstimer;
data test;
     do seq=1 to 500000;
          x1=rannor(9347957);
          *x2=rannor(876769)+0.1*x1;
          epsilon=rannor(938647)*0.5;
          y = 1.5 + 0.5*x1 +epsilon;
          output;
     end;
run;

/* Method 0.*/
sasfile test load;
data res0;
        set test;
  array _x{3,3} _temporary_ ;
  array _a{3,3} _temporary_ ;
  array _tempval{5, 20} _temporary_ ;
  m=mod(_n_-1, 20)+1;
  _tempval[1, m]=x1; 
  _tempval[2, m]=y;
  _tempval[3, m]=x1**2;
  _tempval[4, m]=x1*y;
  _tempval[5, m]=y**2;

  link filler;
  if _n_>=20 then do;
        if _n_>20 then do; 
                   m2=mod(_n_-20, 20)+1;
       _x[1,2]+(-_tempval[1, m2]);
       _x[1,3]+(-_tempval[2, m2]);
       _x[2,2]+(-_tempval[3, m2]);
       _x[2,3]+(-_tempval[4, m2]);
       _x[3,3]+(-_tempval[5, m2]);
     end;
        do i=1 to dim(_a, 1);
         do j=1 to dim(_a, 2);
          _a[i, j]=_x[i, j];
      end;
     end;
            
              do k=1 to dim(_a, 1)-1;
         link adjust;
              end;
     Intercept=_a[1,3]; beta=_a[2,3];
     keep seq   intercept  beta;
     output;
  end;

  return;
filler:
   _x[1,1]=20; _x[1,2]+x1; _x[1,3]+y;
   _x[2,2]+_tempval[3,m];  _x[2,3]+_tempval[4,m]; _x[3,3]+_tempval[5,m];
   _x[2,1]=_x[1,2]; _x[3,1]=_x[1,3]; _x[3,2]=_x[2,3]; 
return;

adjust:
    B=_a[k, k];
 do j=1 to dim(_a, 2);
     _a[k, j]=_a[k, j]/B;
 end;
 do i=1 to dim(_a, 1);
     if i ^=k then do;
          B=_a[i, k];
    do j=1 to dim(_a, 2);
        _a[i, j]=_a[i, j]-B*_a[k, j];
    end;
  end;
 end;
return;

run;
sasfile test close;



/* Method 1.*/

sasfile test load;
data rest0;
        set test;
  array _x{4} _temporary_;
  array _a{2,20}  _temporary_;
  m=mod(_n_-1, 20)+1;
  _a[1, m]=x1; _a[2,m]=y;
  link filler;

  m2=mod(_n_-20, 20)+1;
  if _n_>=20 then do;
    if _n_>20 then do;
              link deduct;
    end;
    beta=(_x[2]-_x[1]*_x[4]/20)/(_x[3]-_x[1]**2/20);
    intercept=_x[4]/20 - beta*_x[1]/20;
    keep  seq   intercept  beta  ;
    output;
  end;
  return;       
filler:
     _x[1]+x1;
  _x[2]+x1*y;
  _x[3]+x1**2;
  _x[4]+y;
return;
deduct:
     _x[1]=_x[1]-_a[1,m2]; 
  _x[2]=_x[2]-_a[1,m2]*_a[2,m2];
  _x[3]=_x[3]-_a[1,m2]**2;
  _x[4]=_x[4]-_a[2,m2];
return;
run;
sasfile test close;



/* Method 2.*/

%macro wrap;
%let window=20;
%let diff=%eval(&window-0);
data testv/view=testv;
     set test;
       xy=x1*y;  
run;

proc expand data=testv  method=none  out=summary(keep=seq sumxy  sumx1  sumy  ussx1  may  max);
       convert  x1=sumx1/transformout=(movsum &diff);
       convert  xy=sumxy/transformout=(movsum &diff);
       convert  x1=ussx1/transformout=(movuss &diff);
       convert  y =sumy /transformout=(movsum &diff);
       convert  y =may / transformout=(movave &diff);
       convert  x1 =max / transformout=(movave &diff);  
run;

data result1;
     set summary(firstobs=&window);
       beta = (sumxy - sumx1*sumy/&window)/(ussx1 - sumx1/&window.*sumx1);  
       alpha= may - beta*max;
       keep seq  beta  alpha; 
run;
%mend;

%let t0=%sysfunc(datetime(), datetime24.);
*options nosource nonotes;
%wrap;
options source notes;
%let t1=%sysfunc(datetime(), datetime24.);
%put Start @ &t0;
%put End   @ &t1;

 

/* Method 3.*/
%let t0=%sysfunc(datetime(), datetime.);
 
data test2v/view=test2v;
       set test;
       array _x{2, 20} _temporary_ (20*0 20*0);
       k=mod(_n_-1, 20)+1;
       _x[1, k]=x1; _x[2, k]=y;
       if _n_>=20 then do;
          do j=1 to dim(_x, 2);
               x=_x[1, j]; y=_x[2, j];
               output;
               keep seq x y;
            end;
       end;
run;

ods select none;
proc  reg data=test2v  outest=res2(keep=seq x intercept);
         by seq;
         model y = x;
run;quit;
ods select all;

%let t1=%sysfunc(datetime(), datetime.);
%put Start @ &t0;
%put End   @ &t1;


/* Method 4. */
%macro wrap;
options nonotes;
ods select none;
%do i=20 %to 500000;
       %let fo=%eval(&i-19);
       proc reg data=test(firstobs=&fo  obs=&i)  outest=_xres(keep=x1 intercept);
           model y =x1;
       run;quit;
      %if %eval(&i=20) %then %do;
          data res3; set _xres; run;
      %end;
      %else %do;
        proc append base=res3  data=_xres; run;
      %end;
%end;

ods select all;
data res3;
       set res3;
       time=19+_n_;
run;
options notes;
%mend;

%let t0=%sysfunc(datetime(), datetime.);
%wrap;
%let t1=%sysfunc(datetime(), datetime.);
%put Start @ &t0;
%put End   @ &t1;

Numerical variables profiling in very large data set

2011-06-21T00:25:00.005-04:00

Profiling numerical variables is an integral part of data analytics, which generally consists of obtaining standard descriptive statistics such as quantiles, first central moments as well as missing ratio.

It is easily obtainable by using PROC MEANS (or PROC SUMMARY). But when we face very large data with many varibles, we will hit memory wall and very long processing time using default options in PROC MEANS. The most time consuming and memory intensive descriptive statistics is quantile calculation. The default method uses ordered statistics, which is the most accurate but also the most memory intensive and time consuming one.

There are two methods available to handle this scenario, both uses the one-pass method of Jain R. and Chlamtac I [1]. This method is able to obtain fairly accurate estimate of quantiles between P25 and P75, but for quantiles outside this range, the estimates will be more rough but given very large sample, maybe acceptable. Another problem is that this estimator is sensitive to the distribution of underlying data [2]. There are newer methods available that are insensitive to the distribution, such as [3].

In PROC MEANS, we can specify QMETHOD=P2 or QMETHOD=HIST to call this one pass method. Or we can use PROC STDIZE, specifying PCTLMTD=ONEPASS to invoke P2 method to calculate quantiles. Each has its pro and con.

PROC MEANS:
PRO: Multithreaded / rich set of descriptive statistics;
CON: Only pre-defined quantiles, output data set is not user friendly, in need of further manipulation;

PROC STDIZE:
PRO: quantiles of 0 to 100 are available, statistics output data is in a user friendly format;
CON: No multithreads, lack of higher central moments statistics, need to surpress data output explicitly;

Of couse, due to the way the output statistics are organized from PROC STDIZE OUTSTAT=, it can be parallelized very easily, and by examining the actual CPU time, PROC STDIZE is more efficient than PROC MEANS. Besides, higher central moments can be calculated from obtained first 2 central moments and this is done in the small output data set.

Below the difference between default method and the two new approaches is illustrated.


options fullstimer;
data test;
        length id 8;
  array x{100} ;
  do id=1 to 1e5;
      do j=1 to dim(x); x[j]=rannor(0)*1.2; end;
   output;
   drop j;
  end;
run;


proc means data=test noprint q1  qmethod=os;
        var x1-x100;
  output  out=_mean
             mean=mean1-mean100
                   std=std_x1-std_x100
                   q1=q1_x1-q1_x100
                   p95=p95_x1-p95_x100;
run;

proc means data=test noprint q1  qmethod=p2   qmarkers=211;
        var x1-x100;
  output  out=_mean
             mean=mean1-mean100
                   std=std_x1-std_x100
                   q1=q1_x1-q1_x100
                   p95=p95_x1-p95_x100;
run;

ods select none;
proc stdize data=test  
                   out=_null_  outstat=_stat  pctlmtd=onepass  
       nmarkers=211
                   pctlpts=1 5  10  25 50  75  90  95   99;
         var x1-x100;
run;
ods select all;

Take a look at the log:

Due to smaller sample size, the most gain is on the memory side. Default approach of PROC MEANS used 643MB memory, which specifying QMETHOD=P2, the memory usage reduced to only 7.3MB. The most memory efficient is PROC STDIZE with PCTLMTD=ONEPASS, only 0.68MB memory was used.

We also examine the difference on quantile estimates using the Ordered Statistics (OS) method and P2 method in PROC MEANS. There are observable differences but not that significant:

Reference:
[1]. Jain R. and Chlamtac I. (1985), "The Algorithm for Dynamic Calculation of Quantiles and Histograms Without Storing Observations," Communications of the ACM, 28(10), 1076–1085.
[2]. SAS/STAT(R) 9.2 User's Guide, Second Edition
[3]. Alsabti, Khaled; Ranka, Sanjay; and Singh, Vineet (1997), "A One-Pass Algorithm for Accurately Estimating Quantiles for Disk-Resident Data". L.C. Smith College of Engineering and Computer Science - Former Departments, Centers, Institutes and Projects. Paper 4. http://surface.syr.edu/lcsmith_other/4

%HPGLIMMIX macro on large scale HMM

2011-06-07T09:25:00.027-04:00

PROC GLIMMIX is good tool for generalized linear mixed model (GLMM), when the scale is small to medium. When facing a large scale GLMM, such as modeling all ZIPs nested in Counties nested in all 51 States in US, a 64-bit machine with extremely large memory is required and the computing may last for months! In a strictly nested hierarchical model, the variance covariance matrix is very sparse, and taking advantage of this property can accelerate computing by many folds.

The %HPGLIMMIX SAS macro is made for large scale Hierarchical Mixed Models. As an example, a sample data using Gamma Regression is shown below, with all ZIPs in AK, AL, AR, AZ with 2-level hierarchies: State and ZIP within State, total 4 blocks with max 693 columns per block. The reason not all ZIPs and all states are used is simply because PROC GLIMMI blows up on the machine.

Copmaring the estimates and std errors from both runs, they are the same, but drastically different running time of 71sec using %HPGLIMMIX v.s. 35min39sec using GLIMMIX.



1715
1716  options nomprint nomlogic;
1717  %hpglimmix(data=temp2,
1718             stmts=%str(
1719                    class zip  zip_state;
1720                    model y = x ;
1721                    random int zip/subject=zip_state;
1722                      ),
1723             error=gamma,
1724             link=LOG,
1725             options=NOTEST);
NOTEST

       The HPGLIMMIX Macro

Data Set           : WORK.TEMP2
Error Distribution : GAMMA
Link Function      : LOG
Response Variable  : Y


Job Starts at : 06JUN2011:15:51:19
    HPGLIMMIX Iteration History

Iteration    Convergence criterion
    1            0.0081058432  13 sec
    2            0.0004213646  13 sec
    3            2.7137935E-7  13 sec
    4            3.1854799E-9  12 sec

Output from final Proc HPMixed run:
Job Ends at : 06JUN2011:15:52:30
1726  options nomprint nomlogic;
1727
1728  proc glimmix data=temp2;
1729       class zip  zip_state;
1730       model y = x /s  dist=gamma;
1731       random int zip /subject=zip_state;
1732  run;



NOTE: Convergence criterion (PCONV=1.11022E-8) satisfied.
NOTE: PROCEDURE GLIMMIX used (Total process time):
      real time           35:38.93
      cpu time            34:30.90

GLIMMIX output:

              Covariance Parameter Estimates

                                                  Standard
             Cov Parm     Subject      Estimate     Error

             Intercept    zip_state    0.000152    0.000125
             zip          zip_state    0.000066    3.105E-6
             Residual                  0.000405    2.281E-6

             Solutions for Fixed Effects

                                Standard
     Effect       Estimate       Error       DF    t Value    Pr > |t|

    Intercept      6.5873    0.006180        3    1065.95      <.0001
    x              0.003436    0.000218    62634      15.79      <.0001

%HPGLIMMIX output:

                _cov data from HPGLIMMIX

               Obs    CovParm       Subject     Estimate

                1     Intercept    zip_state    0.000152
                2     zip          zip_state    0.000066
                3     Residual                  0.000405

               _soln data from HPGLIMMIX

               Obs    Effect       Estimate      StdErr

                1     Intercept      6.5873      0.006180
                2     x              0.003436    0.000218

Regularized Discriminant Analysis

2011-04-10T18:57:00.008-04:00

Demo SAS implementation of Regularized (Linear) Discriminate Analysis of J. Friedman (1989)[1]. Simpler introduction can be found at [2]. Regularized QDA follows similarly.

To save coding, I called R within SAS to finish the computation. For details to see how to call R within SAS, check here.

Reference:
1. Friedman, J. (1989). Regularized discriminant analysis, Journal of the American Statistical Association 84: 165-175.
2. Friedman, J; Hastie, T; Tibshirani, R (2008). The Elements of Statistical Learning, section 4.3.1

ElasticNet in SAS

2011-04-01T09:11:00.009-04:00

Try out Elastic Net [1] in normal linear regression, using Naive algorithm. Exploring possibilities for GLM Elastic Net in SAS.

1. Zou, H and Hastie, T (2005). Regularization and variable Selection via the Elastic Net, Journal Of The Royal Statistical Society Series B.

Bayesian Computation with SAS (2)

2011-02-01T10:47:00.038-05:00

In Chapter 2 of the book "Bayesian Computation with R" by Jim Albert, the philosoph of Bayesian statistics is introduced using an example where a parameter regarding proportion of heavy sleeping college students, i.e. those have more than 8 hours sleep a day, is studied. According to Bayes'Rule, the posterior is proportional to the product of prior probability and likelihood given prior, where the proportion factor is the likelihood of data over the whole parameter space. As we can see, the key point here is the specification of prior probability.

In order to illustrate this idea, Jim introduced three types of prior probability:
1. Prior probability over a discrete grid of choice, given a weight which serves as hyper-parameter;
2. Prior probability by a chosen probability. In most cases, the chosen probability is a so called conjugate prior becuase the posterior is the same family as the prior. In the book, a beta probability is chosen and it is conjugate in the sense that the posterior is also a beta distribution, only with different parameter;
3. A dense grid of prior probability choices which is called histogram prior;

Using a discrete grid of prior choices, we first generate a list of prior probabilities, in data set PRIOR, that we believe the true proportion of heavy sleepers should be, and assign a weight, i.e. the magnitude of credibility, for each choice in a data set P, which serves as hyper-parameter. We then consolidate the two data sets into one and plot the prior against its corresponding weight.


/* Chapter 2 of Bayesian Computation with SAS */
/* In this chapter, Jim talked about applying 
   Bayesian Rule to learn posterior probability
   using an example where the proportion of 
   students sleeping over 8 hours a day is studies.
   Bayes Rule:
     Posterior(p|data) \prop Prior(p)* likelihood(data|p)
*/
options fullstimer formchar='|----|+|---+=|-/\<>*'  error=3;
/***************************************************/
data prior;
     input prior @@;
cards;
2 4 8 8 4 2 1 1 1 1
;
run;

data p;
     do p=0.05 to 0.95 by 0.1;
     output;
  end;
run;

data p;
     retain sum_prior 0;
     do until (eof1);
        set prior end=eof1;
        sum_prior+prior;
     end;
     do until (eof2);
        merge p  prior  end=eof2;
        prior=prior/sum_prior;
        drop sum_prior;
        output;
     end;
run;

goptions reset=all;
symbol  interpol=Needle  value=none  line=1 color=black;
axis1   label=(angle=90  'Prior Probability') order=(0 to 0.3 by 0.1) minor=none;
axis2   label=('p')   order=(0 to 1 by 0.2)  minor=none;
proc gplot data=p;
      plot prior*p /vaxis=axis1  haxis=axis2;
run;quit;

Next, we calculate the posterior probability given prior and likelihood conditional on prior. In the code below, we observe 11 heavy sleepers out of a total sample of 27. The first DATA STEP calculates the likelihood conditional on each prior choice, the posterior is output in the data set P2, where the normalization factor is the sum of conditional likelihood, scaled by max likelihood for numerical stability. You can observe the shape of posterior against hyper-parameter using GPLOT. The INTERPOL=NEEDLE option in SYMBOL statement minic the histogram style line of PLOT function in R.


%let n0=16;
%let n1=11;

data p;      
     set p  end=eof1  nobs=ntotal;
     if p>0 & p<1 then do;
        log_like=log(p)*&n1 + log(1-p)*&n0;   
     end;
     else 
        log_like=-999*(p=0)*(&n1>0) + (p=1)*(&n0>0);   
     log_posterior=log(prior) + log_like; 
     if log_posterior>max_like then max_like=log_posterior;     
run;

proc means data=p  noprint;
     var log_posterior;
     output out=_max  max(log_posterior)=max;
run;

data p2;
     set _max;
     sum_post=0;
     do until (eof);
        set p  end=eof;
        sum_post+exp(log_posterior-max);         
     end;
     do until (eof2);
        set p end=eof2;
        posterior=exp(log_posterior-max)/sum_post;
        output;
        keep p prior  posterior;
     end;
run;

goptions reset=all;
symbol value=NONE interpol=Needle;
axis1  label=(angle=90  'Posterior Probability');
axis2  label=('p');
proc gplot data=p2;
      plot  posterior*p /vaxis=axis1  haxis=axis2;
run;quit;

After tried discrete grid of prior, beta conjugate prior is introduced. Jim mentioned that by matching the believed 50% and 90% percentiles, the parameter of prior beta distribution can be obtained by try-and-error approach. We note that the believed percentiles data indicates that the beta distribution skewed to 0, so that parameter alpha is less than parameter beta and we can use a binary search to nail down the values, which I came down to approximately 3.4 and 7.5. But I will use the {3.4, 7.4} vector as in the book. Since beta distribution is a conjugate prior for proportion, the computation is very easy. The curves for prior, likelihood and posterior are visualized.


/* beta prior */

%let a=3.4;
%let b=7.4;
%let s=11;
%let f=16;

data p;
      do i=1 to 500;
      p=i/500;
   prior=pdf('beta', p, &a, &b);
   like=pdf('beta', p, &s, &f);
   post=pdf('beta', p, &a+&s, &b+&f);
   drop i;
   output;
   end;
run;

goptions reset=all;
symbol1 interpol=j  color=red  width=1  value=none;
symbol2 interpol=j  color=blue width=2  value=none;
symbol3 interpol=j  color=gree width=3  value=none;
axis1   label=(angle=90  'Density')  order=(0 to 5 by 1) minor=none;
axis2   label=('p')      order=(0 to 1 by 0.2)  minor=none;;
legend label=none   position=(top right inside)  mode=share;
proc gplot data=p;
      plot post*p   prior*p  like*p  /overlay  
                                      vaxis=axis1  
                                      haxis=axis2  
                                      legend=legend;
run;quit;

Since we obtain closed form for posterior distribution, we can conduct inference based on exact calculation; alternatively, we can use simulated random sample.


/* generate 1000 random obs from posterior distribution */
data rs;
      call streaminit(123456);
      do i=1 to 1000;
         x=rand('beta', &a+&s, &b+&f);
   output;
   end;
run;

proc univariate data=rs  noprint;
     var x;
     histogram /midpoints=0 to 0.8 by 0.05  
                cbarline=blue  cfill=white  
                outhistogram=hist;
     title "Histogram of ps";
run;
title;

axis1 label=(angle=90 "Frequency");
proc gbarline data=hist;
      bar _midpt_/discrete  sumvar=_count_  space=0  axis=axis1 ;
   title "Histogram of ps";
run;
title;

/* hist prior */
data midpt;
     do pt=0.05 to 0.95 by 0.1;
        output;
     end;
run;

In the book, Jim also mentioned using so called Histogram prior. The Histogram prior is a densed-version of discrete prior, with the distance between discrete choices filled in with dense choices of the same value, separated only by a small distance. It provides a more refined description of prior believe than the discrete prior. We plot the prior and posterior as well as random sample from posterior to give a visual impression. Note that we use PROC STDIZE method=SUM to normalize prior variable.


data prior;
      input prior @@;
datalines;
2 4 8 8 4 2 1 1 1 1
;
run;
proc stdize data=prior  method=sum  out=prior; var prior; run;

data _null_;
     set midpt nobs=ntotal1;
     call symput('nmid', ntotal1);
     set prior nobs=ntotal2;
     call symput('nprior', ntotal2);
     stop;
run;

data histprior;
      array _midpt{&nmid} _temporary_;
      array _prior{&nprior} _temporary_;
      if _n_=1 then do;
         k=1;
         do until (eof1); 
            set midpt  end=eof1;
            _midpt[k]=pt; if k=2 then lo=_midpt[k]-_midpt[k-1];
            k+1;
         end;
         do k=1 to &nmid;
            _midpt[k]=round(_midpt[k]-lo/2, 0.00001);
            put _midpt[k]=;
         end; 
         k=1;
         do until (eof2);
            set prior  end=eof2;
            _prior[k]=prior;
            put k=  _prior[k]=;
            if k<&nprior then k=k+1;     
         end;
         drop k  lo;   
      end;
   
      do p=1/500 to 1 by 1/500;
         sk=0;
         do k=1 to &nmid;
            if p>=_midpt[k] then sk=min(&nmid, sk+1); 
         end;
         histprior=_prior[sk];
         output;
     end;
     stop;
run;


goptions reset=all;
symbol interpol=hiloj;
axis1  label=(angle=90  "Prior Density")  order=(0 to 0.3 by 0.05);
axis2  label=("p")                        order=(0 to 1   by 0.2);
proc gplot data=histprior;
      plot histprior*p /vaxis=axis1  haxis=axis2;
run;quit;

data histposterior;
     set histprior;
     like=pdf('BETA',p,  &s+1, &f+1);
     post=like*histprior;
run;


%put &s;
goptions reset=all;
symbol interpol=hiloj;
axis1 label=(angle=90  "Posterior Density")  order=(0 to 1 by 0.2);
proc gplot data=histposterior;
      plot   post*p/vaxis=axis1;
run;quit;

proc means data=histposterior  sum;
     var post;
run;

proc stdize data=histposterior  method=sum  out=data1;
      var post;
run;
proc means data=data1  sum;
      var post;
run;

In R, the SAMPLE function is able to do a URS sampling with given probability. SAS's PROC SURVEYSELECT seems can't do this with method=URS. But URS sampling with given probability is very easy to code in SAS DATA STEP with the help of FORMAT.


data fmt;
     merge histposterior(rename=(post=start))  
     histposterior(firstobs=2   in=_2    rename=(post=end));
     if _2;
run;

data fmt;
     set data1   end=eof;
     retain start;
     retain fmtname 'ursp'  type 'n';
     if _n_=1 then start=post;
     else do;
        end=start+post;  label=p;
        keep  fmtname  start  end  label;
        output;
        start=end;
     end;
     if eof then do;
        hlo='O';
        label=.;
        output;           
     end;
run;

proc sort data=fmt  nodupkey  out=fmt; 
     by start end;
run;

proc format cntlin=fmt; run;

data  samp;
      do i=1 to 1000;
         p=input(put(ranuni(78686), ursp.), best.);
         output;
      end;
run;

proc univariate data=samp noprint;
     title "Histogram of ps";
     histogram  p/midpoints=0.15 to 0.65 by 0.05  
                  outhist=samp_hist  
                  vaxislabel="Frequency"  vscale=count;      
run;
title;

Once we obtain posterior distribution of a variable, we can easily predict the value of some statistics given new samples. In the book, depending on the choice of weighting density, two types of predictions are discussed: 1. prior prediction; 2. posterior prediction. The calculation is straight-forward given formula.


/* prediction */
data prior;
     input prior @@;
cards;
2 4 8 8 4 2 1 1 1 1
;
run;

data p;
     do p=0.05 to 0.95 by 0.1;
     output;
  end;
run;

data p;
      retain sum_prior 0;
   do until (eof1);
      set prior end=eof1;
   sum_prior+prior;
   end;
   do until (eof2);
      merge p  prior  end=eof2;
   prior=prior/sum_prior;
   drop sum_prior;
   output;
   end;
run;

/* prediction using discrete prior */
%let m=20;
data pred;
        do ys=0 to &m;         
      pred=0;
            do i=1 to 10;
                set p   point=i;  
    pred = pred + prior * pdf('BINOMIAL', ys, p, &m);
    put p=  prior=  pred=  ;
   end;
   output;   
  end;
  stop;
run;

/* prediction using beta(a, b) prior */
%let m=20;
%let a=3.4;
%let b=7.4;
data pred;
     do ys=0 to &m;
       pred=comb(&m, ys) *exp(logbeta(&a + ys, &b+&m -ys) - logbeta(&a, &b));
       output;
     end;
run;

/* prediction by simulation */
%let n=1000;
data pred;
        do i=1 to &n;
      p=rand('BETA', &a, &b);
   ys=rand('BINOMIAL', p, &m);
   drop i;
   output;
    end;
run;

proc freq data=pred  noprint;
        table ys /out=y_freq;
run;
/*
proc gbarline data=y_freq;
        bar y /discrete sumvar=percent;
run;quit;
*/
data y_freq;
        set y_freq;
  percent=percent/100;
  label ys='y';
run;

goptions reset=all  border;
axis1  label=(angle=90  "Predictive Probability")  order=(0 to 0.14 by 0.02);
axis2  label=("y")   order=(0 to 16) offset=(15 points, 15 points) minor=none;
symbol  interpol=needle;
proc gplot data=y_freq;
        plot  percent * ys/vaxis=axis1  haxis=axis2  ;
run;quit;

/* summarize discrete outcome with given coverage probability */
/* input dsn is y_freq from PROC FREQ. The following can be used as gauge.
data y_freq;
        input y percent;
datalines;
0  0.013
1  0.032
2  0.065
3  0.103
4  0.102
5  0.115
6  0.114
7  0.115
8  0.095
9  0.083
10 0.058
11 0.036
12 0.029
13 0.014
14 0.015
15 0.006
16 0.005
;
run;
*/
%let covprob=0.9;
proc sort data=y_freq  out=y_freq2;
        by  percent;
run;

data y_freq2;
        set y_freq2;
  retain cumprob  0;
  cumprob+percent;
  if cumprob<(1-&covprob) then select=0;
  else select=1;
run;

proc means data=y_freq2  sum;
       where select=1;
    var percent;
run;

In Summary:

0. There different choices for specifying Prior believe that can be incorporated into Bayesian computation:

a.) A discrete grid of possible values, each value with some probability

b.) An appropriate distribution, where the parameters is determined from existing information

c.) A dense grid of possible values that span the whole range

1. We demonstrated calculation under these three scenarios and their impact on posterior distribution of statistics of interests;

2. We can conduct Sampling With Replacement with given sampling probability by using FORMAT and uniform random variable generator

3. We demonstrated that prediction in a Bayesian framework is straighforward, but depending on whether we use prior distribution or posterior distribution

Bayesian Computation with SAS (1).

2011-01-04T07:18:00.021-05:00

The book "Bayesian Computation with R" by Jim Albert is an easy to read entry level book on applied Bayesian Statistics. While the book was written for R users, it is not difficult to translate the languages between R and SAS and I believe it is a good way to show Bayesian capability of SAS. In the next several months, I am going to translate all R code in the book into SAS. Readers are encouraged to buy this book to understand what's behind the code. The posts here will only spend minimum effects to explain the statistics underlying while most resource will still be denoted to SAS coding.

This post will cover Chapter 1.

We first follow section 1.2.2, using PROC IMPORT to read in TAB deliminated file, which corresponds to read.table(file, sep='\t', header=T) in R. We also designated PDF as Output Destination.


/* Chapter 1 of Bayesian Computation with R */

%let folder=C:\Documents\SAS for Bayesian Compu with R;
%let datafolder=&folder\data;

libname data "&datafolder";
options missing='.' formchar = "|----|+|---+=|-/\<>*"  errors=2 fullstimer;

/*@@@@@@ Section 1.2.2 @@@@@*/
/*--> Read Tab Deliminated Data into SAS */
proc import datafile="&datafolder/studentdata.txt"  out=student dbms=dlm   replace;
      datarow=2;       /* redundent with GETNAMES= statement */
   guessingrows=2;  /* Only use first 2 rows of data to guess the type of data */
      delimiter='09'x;  /* this is redudent if we specify DBMS=TAB in PROC IMPORT option */
   getnames=YES;  /* This is the same as 'header=True' in R's table.read function */
run;

ods pdf file="&folder/Chapter1.pdf";
/*-->  To see Variable Names and print first row of this data 
   Use SPLIT="*" to specify the split character being "*", which controls 
   line breaks in column headings
*/
proc print data=student(obs=1)  split="*"; run;

SAS output looks like:

Section 1.2.3, summarize frequency tables and use Barplot to visualize the counts.

PROC FREQ directly corresponds to table() function in R, but provides much richer functionality. PROC GBARLINE is used to answer what barplot() does in R. The LABEL= statement in SAS is a very handy tool for annotation purpose for PRINT or GRAPH.

PROC MEANS somehow corresponds to summary() function to obtain summary statistics for quantitative variables. You can sepecify MAXDEC= to keep desired number of decimals in print out.
Histogram can be obtained in SAS via either PROC UNIVARIATE or newly introduced PROC SGPLOT. PROC UNIVARIATE provides more control though. As shown below, you can specify desired mid-points while you can't do so in SGPLOT. On the other hand, SGPLOT is very convenient if you want to overlay density curve, either a normal fit or kernel one.



/*@@@@ section 1.2.3 @@@*/
/*--> To summarize and graph a single Batch */
proc freq data=student;
     table Drink;
run;

/*--> Barplot on single variable */
title "Original idea, summarize and plot";
proc freq data=student  noprint;
      table Drink /out=t;
run;
goptions reset=all noborder;
proc gbarline data=t;
      bar  Drink /discrete space=4  sumvar=count;
run;quit;
title;
/*--> Actually, in SAS, it is better to directly use 
    PROC GBARLINE */
goptions reset=all noborder;
title "In SAS, PROC GBARLINE directly summarize";
proc gbarline data=student;
      bar Drink /discrete space=3 ;
run;quit;
title;

data student;
      set student;
   hours_of_sleep=WakeUp-ToSleep;
   /* Label can be useful for annotation purpose */
   label hours_of_sleep="Hours Of Sleep"; 
run;

proc means data=student 
           /* Request same statistics as in the book */
           min q1 median q3 max  nmiss             
           /* MAXDEC= specifies keeping 2 decimals */
           maxdec=2;                          
    var hours_of_sleep;
run;

/*--> Histogram using PROC UNIVARIATE */
proc univariate data=student  noprint;
      var hours_of_sleep;
   histogram /midpoints=2.5 to 12.5 by 1;
run;

/* new SGPLOT is a handy way to draw histogram and density */
proc sgplot data=student;
      histogram hours_of_sleep/scale=count;     
      density   hours_of_sleep/scale=count  type=normal;
run;

Compare SAS outputs:

Boxplot to visualize distribution of quantitative variables by some classification variable. Here, formula such as var1*var2 style corresponds to var1~var2 style in R.



/*@@@ section 1.2.4 @@@@*/
/*---> Compare Batches using Boxplot */
proc sort data=student out=studentsorted; by Gender; run;
proc boxplot data=studentsorted;
      plot  hours_of_sleep*Gender;
run;   

proc means data=student  
            min q1 median q3 max nmiss  
            nway  maxdec=2;
     class Gender;
  var  Haircut;
run;

SAS output looks like:

R has built-in jitter function, however, we have to make our own in SAS. This verision doesn't apply Fuzzy yet, only for demo purpose.



%macro jitter(var, newname, data=last, factor=1, amount=);
/* follow the JITTER function in R, no fuzzy applied yet . 
   if amount is given, then use this value for disturbance
   else if amount is 0, then use the range of &var for disturbance
   else if amount is NULL, then use the smallest difference among
        distinct values of &var for disturbance. if all values are 
        the same, then use the value obtained from range

   if range of &var is 0, then use the lower range for disturbance
   otherwise 
*/
%let blank=%str( );
%if &data=' ' %then %let data=last;
%local fid;

data _null_;
     fid=round(ranuni(0)*10000, 1);
  call symput('fid', compress(fid));
run;
proc means data=&data  noprint;
     var &var;
  output  out=_util&fid  
          range(&var)=range  
          min(&var)=min  
          max(&var)=max;
run;
data _util&fid;
     set _util&fid;
  if range^=0 then z=range; else z=min;
  if z=0 then z=1;
run;

%if %eval(&amount=&blank) %then %do;
    proc sort data=&data.(keep=&var  where=(&var^=.))   
              out=_xuti&fid  nodupkey;
       by &var;
 run;
 data _duti&fid;
      set _xuti&fid  nobs=ntotal  end=eof;      
   array _x{2} _temporary_;
   if ntotal=1 then do;
      amount=&factor/50*abs(&var);
   keep amount;
   output;
   stop;
   end;
   else do;
      if _n_=1 then do; 
         _x[1]=&var; _x[2]=constant('BIG'); 
      end;
   else do;
               _x[2]=min(_x[2], &var - _x[1]);
      _x[1]=&var;
      if eof then do;
         amount=&factor/5*abs(_x[2]);
      keep amount;
      output;
      end;
   end;
   end;
  run;
    
%end;
%else %if %eval(&amount=0) %then %do;
    data _duti&fid;
      set _util&fid;
   amount=&factor*z/50;
   keep amount;
   output;
 run;     
%end;
%else %do;
    data _duti&fid;
      amount=&amount;
   keep amount;
   output;
 run;
%end;

proc sql noprint;
     select name into :keepvars separated by ' '
  from   sashelp.vcolumn
  where  libname='WORK' 
    and  memname=%upcase("&data") 
    and  memtype="DATA"
  ;
quit;
data &data;     
  array _x{1} _temporary_;
     if _n_=1 then do;
     set _duti&fid;
  _x[1]=amount;
  end;
     set &data;
  &newname=&var + ranuni(0)*(2*_x[1])-_x[1];
  label &newname="jitter(&var)";
  keep &keepvars  &newname;
run;
proc datasets library=work nolist;
     delete _duti&fid _xuti&fid _util&fid;
run;quit;
%mend;

Now we can apply jitter function to the two variables studied in the book and check them out. For graphics, we use axis statement to define how we want the two axes to appear in print out, and use symbol statement to define the appearance of symbols for data. Here we ask for black circles, not connected (interpol=none).



%jitter(hours_of_sleep, xhos, data=student, factor=1, amount=); 
%jitter(ToSleep, xts, data=student, factor=1, amount=);

goptions border  reset=all;
axis1  label=("jitter(ToSleep)")    
       major=(height=2  number=5)   minor=none;
axis2  label=(angle=90  "jitter(Hours of Sleep")     
       major=(height=2  number=5)  minor=none;
symbol1 value=circle  interpol=none  color=black;
proc gplot data=student;
      plot xhos*xts/haxis=axis1  vaxis=axis2;
run;quit;

In order to overlay a fitted regression line, we have multiple ways to do it in SAS, and the Statistical Plot procedures are particularly handy to generate this type of figures. Compare the code and corresponding figures (in the same order of LINEPRINTER, ODS GRAPHICS, Manually Drawn, SGPLOT). For manually drawn figures, since we ask for connected dots in order to obtain a line overlay on circles, we have to sort the data by X-axis variable to make the line as straight as possible. The manually drawn figure is also the one closest to the R figure in the book.




/* old school LINEPRINTER option*/
proc reg data=student  lineprinter;
      title "Line Printer Style";
      model   Hours_of_Sleep=ToSleep /noprint;
   var     xhos  xts;   
   plot pred.*xts='X'  xhos*xts/overlay  symbol='o';
run;quit;
title;

/* Using ODS GRAPHICS */
ods select none; /* to suppress PROC REG output */
title "ODS Graphics";
ods graphics on;
proc reg data=student ;
      model   Hours_of_Sleep=ToSleep;
   var     xhos  xts;
   output  out=student  pred=p;  
      ods select FitPlot; 
run;quit;
ods graphics off;
title;
ods select all;

/* Following standard R approach, tedious in SAS */
proc sort data=student  out=studentsorted;
      by ToSleep;
run;
goptions border  reset=all;
axis1  label=("jitter(ToSleep)")      order=(-2 to 6 by 2)  major=(height=2  number=5)   minor=none;
axis2  label=(angle=90  "jitter(Hours of Sleep")    offset=(, 2 pct)  major=(height=2  number=5)  minor=none;
symbol1 value=circle  interpol=none  color=black;
symbol2 value=none  interpol=line  color=black;
proc gplot data=studentsorted;
      plot xhos*xts  p*xts/overlay  haxis=axis1  vaxis=axis2;   
run;quit;      

/* Using build-in REGRESSION capability of GPLOT by specifying
   INTERPOL=R<><0><> option in SYMBOL statement. 
      In the first <>, you identify the type of regression: 
         L=linear, Q=quadratic, C=cubic
      In the second <>, you specify if intercept is omitted;
      In the third <>, you specify type of CI:
         CLM=CI of mean predicted value;
         CLI=CI of individual points;
 */
goptions border reset=all;
axis1  label=("jitter(ToSleep)")      order=(-2 to 6 by 2)  major=(height=2  number=5)   minor=none;
axis2  label=(angle=90 "jitter(Hours of Sleep")    offset=(, 2 pct)  major=(height=2  number=5)  minor=none;
symbol1 value=circle  interpol=rlclm95  color=black  ci=blue  co=blue;
proc gplot data=student;
     plot xhos*xts/regeqn /* haxis=axis1  vaxis=axis2*/;
run;quit;

/* New SGPLOT has build in regression line capability: REG, SPLINE, LOESS */
proc sgplot data=student;
      scatter  x=xts  y=xhos;
   reg  x=ToSleep  y=Hours_of_Sleep;
run;

Now, to study the robustness of T-Statistics which relies on several assumptions, such as normality, homoskedasticity, etc. First, we need to define a "function" to calculate T-statistics like in the book.



/* section 1.3 Robustness of T-Statistics */
/* section 1.3.2 Write a function to compute T-Statistics */
%macro tstatistics(vars, dsn=last, outdsn=, class=);
%if &class='' %then %let classstmt=;
%else %let classstmt=class &class;

%let nvars=%sysfunc(count(&vars, ' '));

%let var1=%scan(&vars, 1, ' ');
%let var2=%scan(&vars, 2, ' ');
proc means data=&dsn noprint  nway;
     &classstmt;
  var &vars;
  output  out=tdata(where=(_STAT_ in ('STD', 'MEAN', 'N')));
run;
data &outdsn;
     set tdata; by &class _TYPE_;
  retain  m  n  mean1 mean2  std:;
  select(compress(_STAT_));
          when('N') do;
        m=&var1;
     n=&var2;
    end;
    when('MEAN') do;
              mean1=&var1;
     mean2=&var2;
    end;
          when('STD') do;
        std1=&var1;
     std2=&var2;
    end;
    otherwise;
  end;
     if last._TYPE_ then do;
        sp=sqrt(((m-1)*std1**2 + (n-1)*std2**2)/(m+n-2));
     t= (mean1-mean2)/(sp * sqrt(1/m + 1/n));
  keep &class _TYPE_ m  n std:  mean:  sp t;
  output;
  end;
run;
%mend;

Just like the "source(file)" command in R, in SAS, we have a counterpart called %include(or in abv. %inc). Suppose you have a SAS file called "tstatistics.sas" in some folder &folder, the following statement will read in the source code from this file. Since this is a SAS file, suffix ".sas" is not necessary.


filename BCR "&folder";
%inc BCR(tstatistics);

The following code pieces conduct the four simulate scenarios illustrated in the book and draw the last figure.


data test;
      input x y;
cards;
1  5
4  4
3  7
6  6
5  10
;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=);

/* section 1.3.3 Monte Carlo study on the Robustness of T-Statistics */
%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
      retain seed1 99787  seed2 76568;
   do iter=1 to &Niter;
      do j=1 to max(&n, &m);
            if j<=&m then call rannor(seed1, x); else x=.;
   if j<=&n then call rannor(seed2, y); else y=.;
   keep iter x y;
   output;
   end;
    end;
run;

%tstatistics(x y, dsn=test, outdsn=tdata2, class=iter);

data _null_;
      set tdata2  end=eof;
   retain n_reject 0;
   if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
   if eof then do;
      call symput('n_reject', n_reject/_n_);
   end;
run;
%put &n_reject;


/* section 1.3.4 Study the robustness of T-Statistics */

%let alpha=0.1;
%let m=10;
%let n=10;
%let Niter=10000;
data test;
      retain seed1 987687  seed2 76568;
   array x[4];
   array y[4];
      do iter=1 to &Niter;
         do j=1 to max(&n, &m);
      do k=1 to 4;
      select (k);
         when(1) do;
                      if j<=&m then call rannor(seed1, x[k]); else x[k]=.;      
                      if j<=&n then call rannor(seed2, y[k]); else y[k]=.;
      end;
      when(2) do;
                      if j<=&m then call rannor(seed1, x[k]); else x[k]=.;
                      if j<=&n then call rannor(seed2, y[k]); else y[k]=.;
                      y[k]=y[k]*10;
      end;
      when(3) do;
                      if j<=&m then x[k]=rand('T', 4); else x[k]=.;
                      if j<=&n then y[k]=rand('T', 4); else y[k]=.;  
      end;
      otherwise do;
                      if j<=&m then call rannor(seed1, x[k]); else x[k]=.;
                      if j<=&n then call ranexp(seed2, y[k]); else y[k]=.;
                      x[k]=10+x[k]*2;  y[k]=y[k]/0.1;
      end;
     end;           
   end;
            keep iter x1-x4 y1-y4;
            output;
         end;
      end;
run;

%tstatistics(x1 y1, dsn=test, outdsn=tdata1, class=iter);

%tstatistics(x2 y2, dsn=test, outdsn=tdata2, class=iter);

%tstatistics(x3 y3, dsn=test, outdsn=tdata3, class=iter);

%tstatistics(x4 y4, dsn=test, outdsn=tdata4, class=iter);

data _null_;
      n_reject=0;
      do until (eof1);
         set tdata1  end=eof1  nobs=n1;
         retain n_reject 0;
         if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
         if eof1 then do;
             *call symput('n_reject', n_reject/_n_);
       p_reject=n_reject/n1;
       put "NOTE>> Prob(Rejection)=" p_reject  8.4;
         end;
  end;
  n_reject=0;
  do until (eof2);
         set tdata2  end=eof2  nobs=n2;
         retain n_reject 0;
         if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
         if eof2 then do;
             *call symput('n_reject', n_reject/_n_);
       p_reject=n_reject/n2;
       put "NOTE>> Prob(Rejection)=" p_reject  8.4;
         end;
  end;

  do until (eof3);
         set tdata3  end=eof3 nobs=n3;
         retain n_reject 0;
         if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
         if eof3 then do;
             *call symput('n_reject', n_reject/_n_);
       p_reject=n_reject/n3;
       put "NOTE>> Prob(Rejection)=" p_reject  8.4;
         end;
  end;

  do until (eof4);
         set tdata4  end=eof4  nobs=n4;
         retain n_reject 0;
         if abs(t)>quantile('T', 1-&alpha/2, n+m-2) then n_reject+1;
         if eof4 then do;
             *call symput('n_reject', n_reject/_n_);
       p_reject=n_reject/n4;
       put "NOTE>> Prob(Rejection)=" p_reject  8.4;
         end;
  end;
run;

/* draw the overlay density curves of empirical T-statistics and T-distribution */
proc sort data=tdata2;
      by t;
run;
data tdata2;
      set tdata2;
   dt=pdf('T', t, 18);
run;

ods select none;
ods output Controls=Controls;
ods output UnivariateStatistics=UniStat;
proc kde data=tdata2  ;
      univar  t /out=density  method=snr  unistats;
run;
ods select all;

data density;
      set density;
   dt=pdf('T', value, 18);
run;

data _null_;
      set UniStat;
   if Descr='Bandwidth' then call symput ('bw', compress(round(t, 0.0001)));
run;

goptions reset=all border;

title "Comparing Empirical and Theoretical Densities";
legend label=none   position=(top right inside)  mode=share;
axis1  label=("N=&Niter Bandwidth=&bw") 
       order=(-4 to 8 by 2) offset=(1,1)  
       major=(height=2) minor=none;
axis2  label=(angle=90  "Density")     
       order=(0 to 0.4 by 0.1)  offset=(0, 0.1) 
       major=(height=2)  minor=none;
symbol1  interpol=join  color=black  value=none  width=3;
symbol2  interpol=join  color=grey   value=none  width=1;
proc gplot data=density;
      plot density*value  dt*value/overlay  legend=legend
                                   haxis=axis1  vaxis=axis2  ;
      label dt='t(18)' density='Exact';   
run;quit;

ods pdf close;

Key Points Summary:

0. We can use FORMCHAR option to eliminate weird characters

1. We can use PROC IMPORT to import deliminated data and it provided flexible options to control the output;

2. We can use PROC FREQ to summarize categorical data

3. We can use PROC MEANS and PROC UNIVARIATE to generate descriptive statistics for numerical data

4. We can study the dispersion of numerical variable by categorical variable using BOXPLOT

5. We provided a JITTER function for better graphical presentation of heavily overlapped data

6. We showed that there are multiple ways to visualize scatter plots together with fitted regression lines:

6.1 LINEPRINTER in PROC REG
6.2 Using PROC REG to obtain regression line and overlay it on the scatter plot
6.3 Using ODS Graphics in PROC REG
6.4 Using INTERPOL=R<><><> option in SYMBOL statement
6.5 Using PROC SGPLOT

7. We demonstrated how to conduct robust T-test in SAS following the code in the book.

Summarize Numerical Data in a Rolling Window

2010-10-29T03:51:00.006-04:00

Obtain summary statistics over a rolling window for a given data, usually on a time dimension, is not quit easy in SAS, especially the rolling window may contain different number of records and the maximum number is unknown without pass the data once. For example, given a transaction data over several days, a business analyst wants to summarize the data for each 24 hour period. This is actually a recent question on SAS-L.

There are several approaches. Typically people use an array that is large enough to handle a reasonable guess of maximum number of records within the given interval. Or use a hash table so to manage the data cells dynamically. The coder basically needs to build a stack which is first in first out. Both of these methods are not easy to code and error prone for a beginner.

This rolling window problem can be efficiently solved by using MultiLabel Format, a feature tend to be ignored. Here are two example. We want to summarize for a rolling window of size 5. Pay attention to the second example.


data fmt;
      retain fmtname 'rollwindow'  type 'n'  hlo 'M';
      do start=1 to 10;
         end=start+5; 
        label=cats('time', start);
        output;
      end; 
      hlo='O'; label='Out-Of-Bound';
      output;
run;

data dsn;
      do time=1 to 20;
         x=rannor(0);
         y=ranuni(0);
         output;
      end;
run;

proc format cntlin=fmt; run;

proc means data=dsn noprint;
      class time /preloadfmt  mlf; 
      format time rollwindow.;
      var x y;
      output  out=summary_roll mean(x y)=  std(x y)= /autoname;
run;
    

data dsn2;
      do time=1 to 20;
         k=ranpoi(10, 10);
         do j=1 to k;
            time=time+j/(k+1);
            x=rannor(0); y=ranuni(0);
            output;
         end;
      end;
      drop k j;
run;
proc means data=dsn2  noprint;
      class time/preloadfmt mlf  exclusive;
      format time rollwindow.;
      var x y;
      output  out=summary_roll2  mean()=  std()=/autoname;
run;


/**************************************************** 
          non-rolling but shrinking time window, 
          similar for growing time window 
*****************************************************/
data fmt2;
      retain fmtname 'winx'  type 'n'  hlo 'M';
      do start=1 to 10;
         end=18; 
        label=cats('time', start);
        output;
      end; 
      hlo='O'; label='Out-Of-Bound';
      output;
run;      

proc format cntlin=fmt2  cntlout=fmt_all;
run;

proc means data=dsn2  noprint;
      class time/preloadfmt mlf  exclusive;
      format time winx.;
      var x y;
      output  out=summary_roll3  mean()=  std()=/autoname;
run;

**************** An example **************;
data TradeDate;
input TradeDate yymmdd10.;
format TradeDate yymmdd10.;
cards;
2007-01-04
2007-01-05
2007-01-08
2007-01-09
2007-01-10
2007-01-11
2007-01-12
2007-01-15
2007-01-16
2007-01-17
2007-01-18
2007-01-19
2007-01-22
2007-01-23
2007-01-24
2007-01-25
2007-01-26
2007-01-29
2007-01-30
2007-01-31
 ;
run;
data raw;
input id $ Date_S yymmdd10. +1 Date_e yymmdd10. Buy;
format Date_S Date_E yymmdd10.;
cards;
A001 2007-01-09 2007-01-24 24.5
A001 2007-01-12 2007-01-16 56.6
 ;
run;

/*------------- Desired Output ------------*
id      Date_S       Date_E        Buy    Hold_Days
A001 2007-01-09 2007-01-24 24.5    12
A001 2007-01-12 2007-01-30 56.6     3
-------------------------------------------*/
data fmt;
     set raw;
     retain fmtname 'tdate' type 'n'  hlo 'M';
     start=Date_S; end=Date_e;
     label=cats(ID, _n_);
run;
proc format cntlin=fmt  out=fmt_ref;
run;

proc means data=tradeDate noprint  nway;
     class TradeDate/mlf exclusive preloadfmt ;
     format TradeDate tdate.;
     var TradeDate;
     output  out=_test  n()=_freq_;
run;

When the data is getting large, there will be some computing difficulty and insufficient resource. To solve this problem, simply sort and divide the original data into smaller pieces with an overlap equals to the size of rolling window. When combine the summarized pieces, you should discard the overlap part from the rest pieces.

The multilabel approach has an overhead that is smaller overall to the other methods and can be easily changed to accommandate other sizes of rolling windows.

Low Rank Radial Smoothing using GLIMMIX and its Scoring

2010-09-30T19:13:00.004-04:00

Low Rank Radial Smoothing using GLIMMIX [1], a semiparametric approach to smooth curves [2]. Specifying TYPE=RSMOOTH option in RANDOM statement, we can implement this spline smooth approach. The best thing is that for future scoring, data preparation is extremely easy by using the OUTDESIGN= & NOFIT options in v9.2 PROC GLIMMIX, then use PROC SCORE twice on this design matrix to score the fixed effects design matrix X and the random effects design matrix Z, respective, add up together is the score from this radial smoothing method.

[Coming soon]


proc glimmix data=train_data  absconv=0.005;
     model y = &covars /s;
     random &z /s type=rsmooth  knotmethod=equal(20);
run;

proc glimmix data=test  nofit  outdesign=test2;
     model y=&covars /s;
     random &z /s type=rsmooth knotmethod=equal(20);
run;


proc score data=test2  score=beta_fix  type=parms  out=score_fix;
     var  &covars;
run;

proc score data=test2 score=beta_random type=parms  out=score_random;
     var _z:;
run;

Reference:

1. SAS Institute, Statistical Analysis with the GLIMMIX procedure Course Notes, SAS Press, SAS Institute
2. D Rupper, M.P. Wand, R.J. Carroll, Semiparametric Regression, Cambridge University Press, Cambridge, 2003

Generalized Discriminant Analysis

2010-09-30T08:41:00.006-04:00

Using SAS to implement Flexible Discriminant Analysis [Chapter 12, ESL].

Since a discriminant analysis is equivalent to a 2-step process, i.e. regress first then conduct discriminant analysis, it is easy to implement the so called generalized discriminant analysis shown in Ch.12.4--12.6 of Elements of Statistical Learning. The basic idea here is to use some regression analysis procedures, such as using PROC REG and its RIDGE= option in MODEL statement for the ridge regression, and then use the prediction from L2 regularized regression in next step's discriminant analysis. Using PROC GLMSELECT, we can replace L2 regularization with a L1 regularization.

A piece of prototype code looks like this:


PROC GLMMOD data=&yourdata  OUTDESIGN=&design;
     CLASS &dep_var;
     model X = &dep_var /noint;
RUN;
data &yourdata;
     merge &yourdata  &design;
     rename Col1-Col&k = Y1 -Y&k;
run;
%let deps= Y1-Y&k; /* for the case of 5-class problem */
PROC REG  DATA=&yourdata RIDGE=&minridge to &maxridge by 0.1  OUTEST=beta;
MODEL &deps = &covars ;
OUTPUT  OUT=predicted  PRED=&dep._HAT;
RUN;

PROC DISCRIM DATA=predicted &options;
CLASS  &dep;
VAR    &dep._HAT;
RUN;

VARIMAX rotation of PLS loadings

2010-08-09T16:09:00.006-04:00

Partial Least Square is one of several supervised dimension reduction techniques and attracts attention in recent years. In the one hand, PLS is able to generate a series of scores that maximize linear correlation between dependent variables and independent variables, on the other hand, the loading of PLS can be regarded as similar counterpart from factor analysis, hence we can rotate the loadings from PLS therefore eliminate some of the non-significant variable in terms of prediction.


%macro PLSRotate(Loading, TransMat, PatternOut, PatternShort, 
                 method=VARIMAX, threshold=0.25);
/* VARIMAX rotation of PLS loadings. Only variables having 
   large loadings after rotation will enter the final model. 

   Loading dataset contains XLoadings output from PROC PLS 
   and should have variable called NumberOfFactors
   TransMat is the generated Transformation matrix;
   PatternOut is the output Pattern after rotation;
   PatternShort is the output Pattern with selected variables
*/

%local covars;
proc sql noprint;
     select name into :covars separated by ' '
  from   sashelp.vcolumn
  where  libname="WORK" & memname=upcase("&Loading") 
        &   upcase(name) NE "NUMBEROFFACTORS"
  &   type="num"
  ;
quit;
%put &covars;

data &Loading.(type=factor);
         set &Loading;
         _TYPE_='PATTERN';
         _NAME_=compress('factor'||_n_);
run;
ods select none;
ods output OrthRotFactPat=&PatternOut;
ods output OrthTrans=&TransMat; 
proc factor  data=&Loading   method=pattern  rotate=&method  simple; 
         var &covars;
run;
ods select all;

data &PatternShort;
     set &PatternOut;
  array _f{*} factor:;
  _cntfac=0;
  do _j=1 to dim(_f);  
        _f[_j]=_f[_j]*(abs(_f[_j])>&threshold); _cntfac+(_f[_j]>0); 
     end;
  if _cntfac>0 then output;
  drop _cntfac _j;
run;
%mend;

Here I try to replicate the case study in [1] which elaborated how to do and properties of VARIMAX rotation to PLS loadings. The PROC PLS output, after various tweaks on convergence criteria and singularity conditions, is still a little different from the result reported in [1] for factors other than the leading one, therefore, I will directly use the U=PS matrix in pp.215.



data loading;
input factor1-factor3;
cards;
-0.9280  -0.0481  0.2750
0.0563  -0.8833  0.5306
-0.9296  -0.0450  0.2720
-0.7534  0.1705  -0.5945
0.5917  -0.0251  -0.6450
0.9082  0.3345    0.1118
-0.8086  0.4551  -0.3800
;
run;


proc transpose data=loading  out=loading2;
run;

data loading2(type=factor);
     retain _TYPE_ "PATTERN";
  set loading2;
run;


ods select none;
ods output OrthRotFactPat=OrthRotationOut;
ods output OrthTrans=OrthTrans; 
proc factor  data=Loading2   method=pattern  rotate=varimax  simple; 
         var col1-col7;
run;
ods select all;

Reference:
[1] Huiwen Wang; Qiang Liu , Yongping Tu, "Interpretation of PLS Regression Models with VARIMAX Rotation", Computational Statistics and Data Analysis, Vol.48 (2005) pp207 – 219

Table Look Up in SAS, practical problems

2010-08-04T12:06:00.008-04:00

One guy asked in a SAS forum about a typical table look up problem:
He has a data with two IDs:
id1 id2
a b
a e
b c
b e
c e
d e

and he wants to generate a new data set with the following structure according to above information :
id a b c d e
a 0 1 0 0 1
b 1 0 1 0 1
c 0 1 0 0 1
d 0 0 0 0 1
e 1 1 1 1 0

The real data is potentially big.
***************************;
At first look, this is a typical table look up problem SAS programmers facing almost everyday, that is duplicate keyed lookup table. It is a simple one because there is no inherent relationship among records.


data original;
   input id1 $ id2 $;
datalines;
a b
a e
b c
b e
c e
d e
;
run;

proc datasets library=work nolist;
     modify original;
  index create id1 id2;
quit;

proc sql;
     create table all_cases as
  select a.*, monotonic() as seq
  from (
  select distinct id1 as id
  from original
  union
  select distinct id2 as id
  from original
  ) as a
  order by a.id
  ;
quit;

proc sql noprint;
     select id into :idnames separated by ' '
  from   all_cases
  ;
quit;

data new;
  if _n_=1 then do;
     declare hash _h(dataset:'all_cases');
     _h.defineKey('id');
     _h.defineData('seq');
     _h.defineDone();
     end; 
     set all_cases;

  array _a{*} &idnames; 

  id1=id;  
  set original key=id1;      
  _mx_=%sysrc(_sok);
  
  do while (_iorc_=%sysrc(_sok));   
     rc=_h.find(key:id2); if rc=0 then _a[seq]=1;
  id1=id;
     set original key=id1;  
  
  end;
  _ERROR_=0;
  
  id2=id;  
  set original key=id2;      
  do while (_iorc_=%sysrc(_sok)); 
     rc=_h.find(key:id1); if rc=0 then _a[seq]=1;
  id2=id;
     set original key=id2;  
  end;
  _ERROR_=0;
  do j=1 to dim(_a); _a[j]=max(0, _a[j]); end;
  keep id &idnames;
run;

On the other hand, this problem can be solved in a more SASsy way like this:


data original;
   input id1 $ id2 $;
datalines;
a b
a e
b c
b e
c e
d e
;
run;

proc sql;
     create table newx as
     select a.id1, a.id2, (sum(a.id1=c.id1 & a.id2=c.id2)>0) as count
     from   
       (select a.id as id1, b.id as id2 
        from all_cases as a, all_cases as b) as a
left join   original as c
       on   a.id1=c.id1 or a.id2=c.id1
    group by a.id1, a.id2
    ;
quit;

proc transpose data=newx  out=_freq_t name=id2;
     by id1;
     var count;
     id id2;
run;

data _freq_t;
     set _freq_t;
     array _n{*} _numeric_;
     do i=1 to dim(_n);
        _n[i]=(_n[i]>0);
     end;
     drop i;
run;

proc transpose data=_freq_t(drop=id2) out=_freq_t2  name=id1;
     id id1;
run;

proc sql noprint;
     select id1, count(distinct id1) into :covars separated by ' ', :count
     from   _freq_t;  
quit;

data new2;
     set _freq_t;
     array _x{*} &covars;
     array _x2{&count} _temporary_;

     do j=1 to &count; _x2[j]=_x[j]; end;
     set _freq_t2;
     do j=1 to &count; _x[j]=(_x[j]+_x2[j]>0); end;
     drop j id2;
run;

An Economic Approach for a Class of Dimensionality Reduction Techniques

2010-07-30T00:19:00.009-04:00

Just back from KDD2010. In the conference, there are several papers that interested me.

On the computation side, Liang Sun et al.'s paper [1], "A Scalable Two-Stage Approach for a Class of Dimensionality Reduction Techniques" caught my eyes. Liang proves that a class of dimension reduction techniques, such as CCA, OPLS, LDA, etc, that relies on general eigenvalue decomposition, can be computed in a much cheaper way by decomposing the original computation into a least square problem and a much smaller scale eigenvalue decomposition problem. The equivalence of their two stage approach and direct eigenvalue decomposition is rigourously proved.

This technique is of particular interest to ppl like me that only have limited computing resources and I believe it would be good to implement their algorithm in SAS. For example, a Canonical Discriminant Analysis with above idea is demonstrated below. Note also that by specifing RIDGE= option in PROC REG, the regularized version can be implemented as well, besides, PROC REG is multi-threaded in SAS. Of course, the computing advantage is only appreciatable when the number of features is very large.

The canonical analysis result from reduced version PROC CANDISC is the same as the full version.

In fact, this exercise is the answer for Exercise 4.3 of The Elements of Statistical Learning [2]

[1]. Liang Sun, Betul Ceran, Jieping Ye, "A Scalable Two-Stage Approach for a Class of Dimensionality Reduction Techniques", KDD2010, Washington DC.

[2]. Trevor Hastie, Robert Tibshirani, Jerome Friedman, "The Elements of Statistical Learning", 2nd Edition.




   proc format; 
      value specname 
         1='Setosa    ' 
         2='Versicolor' 
         3='Virginica '; 
   run; 
 
   data iris; 
      title 'Fisher (1936) Iris Data'; 
      input SepalLength SepalWidth PetalLength PetalWidth 
            Species @@; 
      format Species specname.; 
      label SepalLength='Sepal Length in mm.' 
            SepalWidth ='Sepal Width in mm.' 
            PetalLength='Petal Length in mm.' 
            PetalWidth ='Petal Width in mm.'; 
      symbol = put(Species, specname10.); 
      datalines; 
   50 33 14 02 1 64 28 56 22 3 65 28 46 15 2 67 31 56 24 3 
   63 28 51 15 3 46 34 14 03 1 69 31 51 23 3 62 22 45 15 2 
   59 32 48 18 2 46 36 10 02 1 61 30 46 14 2 60 27 51 16 2 
   65 30 52 20 3 56 25 39 11 2 65 30 55 18 3 58 27 51 19 3 
   68 32 59 23 3 51 33 17 05 1 57 28 45 13 2 62 34 54 23 3 
   77 38 67 22 3 63 33 47 16 2 67 33 57 25 3 76 30 66 21 3 
   49 25 45 17 3 55 35 13 02 1 67 30 52 23 3 70 32 47 14 2 
   64 32 45 15 2 61 28 40 13 2 48 31 16 02 1 59 30 51 18 3 
   55 24 38 11 2 63 25 50 19 3 64 32 53 23 3 52 34 14 02 1 
   49 36 14 01 1 54 30 45 15 2 79 38 64 20 3 44 32 13 02 1 
   67 33 57 21 3 50 35 16 06 1 58 26 40 12 2 44 30 13 02 1 
   77 28 67 20 3 63 27 49 18 3 47 32 16 02 1 55 26 44 12 2 
   50 23 33 10 2 72 32 60 18 3 48 30 14 03 1 51 38 16 02 1 
   61 30 49 18 3 48 34 19 02 1 50 30 16 02 1 50 32 12 02 1 
   61 26 56 14 3 64 28 56 21 3 43 30 11 01 1 58 40 12 02 1 
   51 38 19 04 1 67 31 44 14 2 62 28 48 18 3 49 30 14 02 1 
   51 35 14 02 1 56 30 45 15 2 58 27 41 10 2 50 34 16 04 1 
   46 32 14 02 1 60 29 45 15 2 57 26 35 10 2 57 44 15 04 1 
   50 36 14 02 1 77 30 61 23 3 63 34 56 24 3 58 27 51 19 3 
   57 29 42 13 2 72 30 58 16 3 54 34 15 04 1 52 41 15 01 1 
   71 30 59 21 3 64 31 55 18 3 60 30 48 18 3 63 29 56 18 3 
   49 24 33 10 2 56 27 42 13 2 57 30 42 12 2 55 42 14 02 1 
   49 31 15 02 1 77 26 69 23 3 60 22 50 15 3 54 39 17 04 1 
   66 29 46 13 2 52 27 39 14 2 60 34 45 16 2 50 34 15 02 1 
   44 29 14 02 1 50 20 35 10 2 55 24 37 10 2 58 27 39 12 2 
   47 32 13 02 1 46 31 15 02 1 69 32 57 23 3 62 29 43 13 2 
   74 28 61 19 3 59 30 42 15 2 51 34 15 02 1 50 35 13 03 1 
   56 28 49 20 3 60 22 40 10 2 73 29 63 18 3 67 25 58 18 3 
   49 31 15 01 1 67 31 47 15 2 63 23 44 13 2 54 37 15 02 1 
   56 30 41 13 2 63 25 49 15 2 61 28 47 12 2 64 29 43 13 2 
   51 25 30 11 2 57 28 41 13 2 65 30 58 22 3 69 31 54 21 3 
   54 39 13 04 1 51 35 14 03 1 72 36 61 25 3 65 32 51 20 3 
   61 29 47 14 2 56 29 36 13 2 69 31 49 15 2 64 27 53 19 3 
   68 30 55 21 3 55 25 40 13 2 48 34 16 02 1 48 30 14 01 1 
   45 23 13 03 1 57 25 50 20 3 57 38 17 03 1 51 38 15 03 1 
   55 23 40 13 2 66 30 44 14 2 68 28 48 14 2 54 34 17 02 1 
   51 37 15 04 1 52 35 15 02 1 58 28 51 24 3 67 30 50 17 2 
   63 33 60 25 3 53 37 15 02 1 
   ; 
   proc candisc data=iris out=outcan distance anova; 
      class Species; 
      var SepalLength SepalWidth PetalLength PetalWidth; 
   run;
 
  ods select none;
  proc glmmod data=iris  outdesign=H(keep=COL:);
           class  Species;
     model SepalLength=Species/noint;
  run;  

  data H;
          merge H   iris;
  run;

/**************************
for efficiency consideration, a view can also be used:
data H/view=H;
     set iris;
     array _S{*} Col1-Col3 (3*0);     
     do j=1 to dim(_S); _S[j]=0; end;
     _S[Species]=1;
     drop j;
run;
****************************/
  proc reg data=H  outest=beta;
          model Col1-Col3 = SepalLength SepalWidth PetalLength PetalWidth;
    output   out=P  p=yhat1-yhat3;
  run;quit;
  ods select all;


  proc candisc  data=P;
          class Species;
    var   yhat1-yhat3;
  run;

Implement Randomized SVD in SAS

2010-07-12T23:06:00.002-04:00

In the 2010 SASware Ballot®, a dedicated PROC for Randomized SVD was among the options. While an official SAS PROC will not be available in the immediate future as well as in older SAS releases, it is fairly simple to implement this algorithm using existing SAS/STAT procedures.

Randomized SVD will be useful for large scale, high dimension data mining problems, for instance Text Mining. In SAS/Base and SAS/STAT, lack of sparse matrix operation capability puts any serious Text Mining task at the edge of infeasibility, such as using LSI or NMF algorithms. Randomized SVD provides an economic alternate solution by sacrificing a little accuracy which is bounded under the three sampling schema proposed by the authors [1], while the code below demos sampling schema 1.



/* Randomized SVD with sampling schema 1. */
%let dim=2048;
%let nobs=1e4;
%let s=256;
data matrix;
     array _x{*} x1-x&dim;
     do id=1 to &nobs;
     do _j=1 to dim(_x); _x[_j]=sin(mod(id, _j))+rannor(id); end;
  output;
  drop _j;
  end;
run;

%let datetime_start = %sysfunc(TIME()) ;
%let time=%sysfunc(datetime(), datetime.); %put &time;
data seed;
     array _x{*} x1-x&dim;
     do _j=1 to dim(_x); _x[_j]=0; end;
  output;
  stop;
run;

proc fastclus data=matrix  seed=seed      out=norm(keep=ID DISTANCE)
              maxiter=0    maxclusters=1  noprint  replace=none;
     var x1-x&dim;
run;
data normv/ view=normv;
     set norm(keep=DISTANCE);
     DISTANCE2=DISTANCE**2;
     drop DISTANCE;
run;
proc means data=normv noprint;
     var DISTANCE2;
     output  out=matrixnorm  sum(DISTANCE2)=Frobenius_sqr;
run;
data prob;
     set matrixnorm ;
  retain Frobenius_sqr;
  do until (eof);
     set norm  end=eof;
  _rate_=DISTANCE**2/Frobenius_sqr;
  keep ID _rate_;
  output;
  end;
run;

data matrixv/view=matrixv;
     merge matrix  prob(keep=_rate_);
run;

ods select none;
proc surveyselect data=matrixv  out=matrixsamp(drop=SamplingWeight  ExpectedHits  NumberHits _rate_)  
                   sampsize=&s  method=pps_wr   outhits  ;
     size _rate_;
run;
ods select all;

proc transpose data=matrixsamp  out=matrixsamp;
     var x1-x&dim;
run;

proc princomp data=matrixsamp  outstat=testv(where=(_type_ in ("USCORE")))  
              noint  cov  noprint;
  var col1-col&s;
run;
data testV_t/view=testV_t;
     retain _TYPE_ 'PARMS';
  set testv(drop=_TYPE_);
run;     

proc score data=matrixsamp   score=testV_t  type=parms  
           out=SW(keep=ID Prin:);
    var   col1-col&s;
run;

data seed;
     array _s{*}  prin1-prin&s;
  do _j=1 to dim(_s); _s[_j]=0; end;
  drop _j; output; stop;
run;

proc fastclus data=SW   seed=seed  maxiter=0  maxc=1  replace=NONE   out=SW2(drop=CLUSTER)  noprint;
     var prin1-prin&s;
run;

data HHT;
     set SW2;
  array _x{*}  prin1-prin&s;
  do _j=1 to dim(_x); _x[_j]=(_x[_j]/distance)**2; end;
  drop _j  distance;
run;

proc transpose data=HHT  out=HHT2(drop=_LABEL_);      
run;
data HHT2; 
     _TYPE_='PARMS';
     set HHT2; 
  rename COL1-COL2048=x1-X2048;
run;

proc score data=matrix   score=HHT2  type=parms  out=P(drop=x:);
     var x:;
run;     

%let time=%sysfunc(datetime(), datetime.); %put &time;
%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss) ;
options notes source;

Reference:
[1], P. Drineas and M. W. Mahoney, "Randomized Algorithms for Matrices and Massive Data Sets", Proc. of the 32nd Annual Conference on Very Large Data Bases (VLDB), p. 1269, 2006.

"Entrywise" Norm calculation using PROC FASTCLUS

2010-06-25T21:48:00.012-04:00

In some data mining applications, matrix norm has to be calculated, for instance [1]. You can find a detailed explanation of Matrix Norm on Wiki @ Here

Instead of user written routine in DATA STEP, we can obtain "Entrywise" norm via PROC FASTCLUS efficiently and accurately.


data matrix;
     input X1-X5;
datalines;
1 2 4 5 6
7 8 9 0 1
2 3 4 5 6
3 4 5 6 7
7 8 9 0 2
2 4 6 8 0
;
run;

data seed;
     input X1-X5;
datalines;
0 0 0 0 0
;
run;

options nosource;
proc export data=matrix  outfile='c:\matrix.csv'  dbms=csv replace; run;
options source;

proc fastclus data=matrix  seed=seed      out=norm(keep=DISTANCE)
              maxiter=0    maxclusters=1  noprint  ;
     var x1-x5;
run;

/* 
In output file NORM, variable DISTANCE is the square root of Frobenius norm. If LEAST=P option is specified, then p-norm is calculated. In PROC FASTCLUS, you can specify p in the range of  [1, \inf].

Now what you got is vector norm for each row, taking the sum of squares of DISTANCE, you obtain the Frobenius norm of the data matrix, which can be easily obtained through PROC MEANS on a data view: 
*/
data normv/ view=normv;
     set norm(keep=DISTANCE);
     DISTANCE2=DISTANCE**2;
     drop DISTANCE;
run;
proc means data=normv noprint;
     var DISTANCE2;
     output  out=matrixnorm  sum(DISTANCE2)=Frobenius_sqr;
run;

You can use the following R code to verify the results;


mat <- read.csv('c:/matrix.csv', header=T)
#verify vector norm
vnorm <- apply(mat, 1, function(x){sqrt(sum(x^2))});
#verify norm of the matrix
x<-as.matrix(mat)
sqrt(sum(diag(t(x)%*%x)))

PS:
1. Of course, above process is designed for implementing the randomized SVD in [1]. If only the matrix Frobenius norm is of interests, you can also use the following code snippet:


data matrixv/view=matrixv;
     set matrix;
     array _x{*}  x1-x5;
     array _y{*}  y1-y5;
     do j=1 to dim(_x);  _y[j]=_x[j]**2; end;
     keep y1-y5;
run;

proc means data=matrixv  noprint;
     var y1-y5;
     output  out=_var(drop=_TYPE_  _FREQ_)   sum()=/autoname;
run;

data _null_;
     set _var;  
     norm=sqrt(sum(of _numeric_));
     put norm=;
run;
/* --LOG WRITES:
norm=28.635642127
NOTE: There were 1 observations read from the data set WORK._VAR.
*/

2. Using its built-in computing engine for Eucleadian Distance, PROC FASTCLUS is also a powerful tool to search for the data point in main table that is CLOEST to the a record in lookup table. This technique is shown Here and [2].

Reference:
[1], P. Drineas and M. W. Mahoney, "Randomized Algorithms for Matrices and Massive Data Sets", Proc. of the 32nd Annual Conference on Very Large Data Bases (VLDB), p. 1269, 2006.

[2], Dorfman, Paul M.; Vyverman, Koen; Dorfman, Victor P., "Black Belt Hashigana", Proc. of the 2010 SAS Global Forum, Seattle, WA, 2010

Boost to tackle nonlinearity

2010-06-01T17:20:00.004-04:00


data nonlinear;
     do x=1 to 627;
     p=(sin(x/100)+1)*0.45;
     do j=1 to 100;
     x1=x+(j-1)/100;
     if ranuni(8655645)<=p then y=1; else y=0;
     output; drop p j;
  end;
  end;
run;

proc rank data=nonlinear  out=nonlinearrank groups=20;
     var x1;
  ranks rank1;
run;

proc means data=nonlinearrank noprint;
     class rank1;
  var y x1;
  output out=_mean(where=(_type_=1))  mean(y)=y  mean(x1)=x1;
run;

%inc "C:\Documents and Settings\lxie\Desktop\SAS Prog and Docs\Boost macro2 ver3.0.sas";
%inc "C:\Documents and Settings\lxie\Desktop\SAS Prog and Docs\predict macro.sas";

%boost2(nonlinear, 1, outputds,  outwts, 100, 3);

%macro pred(niter);
%do i=1 %to &niter;
  %predict(nonlinear, 1, outputds, outlogds, out_pred, sumpred&i, 3, &i);
%end;
data sumpred_all;
     merge %do i=1 %to &niter;
            sumpred&i.(rename=(sum_pred=sum_pred&i))
     %end;;
run;
%mend;

%pred(100);

proc datasets library=work nolist;
     delete sumpred1-sumpred100;
quit;

options nosource;
proc export data=sumpred_all  outfile="c:\sumpred.csv"  dbms=csv replace; run;
options source;

In R, use this code piece to recover the animation:


library(caTools);
test<-read.csv('c:/sumpred.csv', header=T)
minmax<-quantile(as.matrix(test), c(0,1))
sz=314; p=100;
y<-sort(sample(1:62700, size=sz))
image=array(0, c(sz, p, 100))
for (i in (1:100)){
  pic=matrix(0, ncol=p, nrow=sz)
  trace=round((test[y,i]-minmax[1])/(minmax[2]-minmax[1])*p);
  for(j in (1:sz)){
     image[j, trace[j], i]=1
  }
}
write.gif(image, "c:/boost.gif", col=gray(1:2/2))
im = read.gif("c:/boost.gif")
for(i in 1:100){
   image(im$image[,,i], col=(im$col), 
         main=paste('Iter', i), y=1:100, x=1:314, 
         xlab="Index", ylab="Percentage of 1")
}

test2<-test[,-1]-test[-100]
test2<-cbind(test[,1], test2)

par(mfrow=c(2,1))
for (i in (1:100)){
    plot(test2[,i], type='l', ylim=c(-0.5, 1), 
         main=paste('Iteration', i, 'Delta'),
         ylab='Delta')
    plot(test[,i], type='l', ylim=minmax, 
         main=paste('Iteration', i, 'Score'),
         ylab='Score')
    Sys.sleep(0.1)
}

K-Nearest Neighbor in SAS

2010-05-05T18:18:00.007-04:00

K-Nearest-Neighbor, aka KNN, is a widely used data mining tool and is often called memory-based/case-based/instance-based method as no model is fit. A good introduction to KNN can be find at [1], or @ Wiki.

Typically, KNN algorithm relies on a sophisticated data structure called kd-Tree [2] to quickly find the cloeset points for a given point in high dimensional space. While you can find good pseudo-code for kd-Tree implementation and KNN online everywhere, for example [3], it is not trivial to implement your own in SAS [I mean an efficient one]. For most analysts, the first idea is to turn to your $200K-Initial-Fee-AND-$60K-Per-Year Enerprise Miner(R) for this method, however, it turns out that PROC DISCRIM in SAS/STAT is able to do the same thing! Note that annual license fee for SAS/STAT is a tiny fraction of EM and most analysts have ready access to SAS/STAT.

Simply ask PROC DISCRIM to use nonparametric method by using option "METHOD=NPAR K=". Note that do not use "R=" option at the same time, which corresponds to radius-based of nearest-neighbor method. Also pay attention to how PROC DISCRIM treat categorical data automatically. Sometimes, you may want to change categorical data into metric coordinates in advance.

Because KNN is a memory-based method, when you score the test data or new data in production, you have to use the raw table in the scoring process. Test different values of K using Cross-Validation to select the best one [you need a macro loop then.]

PROC DISCRIM uses memory approximately proportional to the second order of number of variables and time usage, excluding I/O time, is roughly proportional to log(N)*(N*P) where N is the number of observations and P is the number of variables used as it uses the tree search algorithm in [4].

Therefore, time consumption is still a concern on large data set. As an experiment, I conducted KNN on data set with records ranging from 5000 to 50000 while scoring the same amount of records at the same time, using 20 features and let k=11, we observe:

obs        time (sec)     memory (KB)
5000        1.03            1068
10000      3.90           1949
15000      8.65            2830
20000    15.48            3709
25000 34.58            4587
30000    70.26            5472
35000 109.83            6350
40000 161.71            7229
45000 208.80            8108
50000 263.58            8987

Empirical time usage is roughly quadratic in the multiplier of per 5000 observations, which means to work on a data set with 300K observations and score the same number of records, SAS needs to take 3.75 hours! Large K value will only increase time used by a very small fraction, though.
Sample code using Iris data from SAS Online Doc for v9.1.3:


ods select none;
proc surveyselect data=iris  out=iris2  
                  samprate=0.5  method=srs  outall;
run;
ods select all;

%let k=5;
proc discrim data=iris2(where=(selected=1))   
             test=iris2(where=(selected=0))
             testout=iris2testout
             method=NPAR k=&k 
             listerr crosslisterr; 
      class Species; 
      var SepalLength SepalWidth PetalLength PetalWidth; 
      title2 'Using KNN on Iris Data'; 
run; 

proc freq data=iris2testout;
     table Species*_INTO_;
run;

Reference:
[1]. Trevor Hastie, Robert Tibshirani, Jerome Friedman, "Elements of Statistical Learning", Ed.2, Springer, 2008;
[2]. J. L. Bentley. "Multidimensional binary search trees used for associative searching", Communications of the ACM, 18(9):509-517, 1975;
[3]. http://simsearch.yury.name/tutorial.html ;
[4]. Friedman, J.H., Bentley, J.L., and Finkel, R.A., "An Algorithm for Finding Best Matches in Logarithmic Expected Time," ACM Transactions on Mathematical Software, 3, 209 - 226. 1977;

Next Project: Regularized Logistic Regression

2010-05-05T11:50:00.006-04:00

L1 Regularized Logistic Regression effectively handles large number of predictors and serves variable selection simultaneously. [1] indicates that L1 RLR can be implemented via IRLS-LARS algorithm. You can tweak PROC GLMSELECT in v9.2 for this.

L2 Reguarlized Logistic Regression can be used to approximate SVM solutions [2], and can be implemented via TR-IRLS as suggested by [3], which is a ridge LR.

Reference:
[1]Su-In Lee, Honglak Lee, Pieter Abbeel and Andrew Y. Ng, "Efficient L1 Regularized Logistic Regression", Proceedings of Annual Conference of American Association for Artificial Intelligence, 2006

[2] Jian Zhang, Rong Jin, Yiming Yang & Alex G. Hauptmann, "Modified Logistic Regression: An Approximation to SVM and Its Applications in Large-Scale Text Categorization", Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003), Washington DC, 2003.

[3]Paul Komarek, "Logistic Regression for Data Mining and High-Dimensional Classification", Ph.D Dissertation, Robotics Institute, Carnegie Mellon University, 2004

Conduct R analysis within SAS

2010-04-30T16:27:00.002-04:00

R is attractive to statistical analysts for its ease of use and ready access of packages implementing modern methodologies. If you have IML, you can submit R commands within SAS/IML enviornment, see Rick's post @ here. Unfortunately, not all analysts have licensed IML. To work around this limitation, I proposed the following technique to submit R statements in a SAS/Base enviornment.

To exchange data between R and SAS, I pushed SAS dataset into a CSV file that can be read by R. Instead, with some more coding, we can leverage the 'sas7bdat' R package to read SAS data directly into R, see Charlie Huang's blog @ here for a demonstration.


%macro RScript(Rscript);
data _null_;
     file "&Rscript";
     infile cards;
     input;
     put _infile_;
%mend;

 
%macro CallR(Rscript, Rlog);
systask command "C:\Progra~1\R\R-2.8.0\bin\R.exe CMD BATCH --vanilla --quiet
                    &Rscript  &Rlog "
        taskname=rjob1  wait  status=rjobstatus1;
%mend;

/****************************/
data a;
     length i 4;
     array _x{100} x1-x100;
     do i=1 to 300;
        do j=1 to dim(_x); _x[j]=rannor(98765); end;
        output;    drop j;
     end;
run;

proc export data=a  outfile="c:\a.csv"  dbms=csv; run;

%RScript(c:\rscript.r)
cards4;
dcsv <- read.csv('c:/a.csv', header=T);
dsvd<-svd(dcsv[,2:101]);
dsvd$u[1:5, 1:5];
dsvd$d[1:8];
;;;;
run;

%CallR(c:/rscript.r, c:/rlog1.txt);

data _null_;
     infile "c:\rlog1.txt";
     input;
     put _infile_;
run;

Improve the Boost macro from Prof. Rayens, W and Dr. Johnson, K

2010-04-21T16:52:00.032-04:00

In Chapter 2 of the book "Pharmaceutical Statistics Using SAS: A Practical Guide" (SAS Press), Prof. Rayens, W and Dr. Johnson K. presented their SAS implementation of boost algorithms, including AdaBoost, RealBoost, GentleBoost and LogitBoost. The original SAS macro can be found at Here.

Their macro uses PROC IML and is preferred by my colleagues since it is not I/O bounded as opposed to my DATA STEP implementation. One shortcoming of their macro is inefficiency because I think it is for DEMONSTRATION purpose rather than for serious applications. The way they wrote the program shows too many redundent calculations as well as inefficient matrix operation. In return, their macro is not able to handle a data set with more than 5000 observations, and on my PC [Intel E6750 2.66GHz, 7.4Gflops/core, 4G memory], it took 1m41s to go over 10 iterations on a data set with 2000 observations and 10 numerical variables and consumed 158MB memory. Note that the resource consumptions increase quadratically in the number of observations. In another experiment with 4000 observations, the macro consumed 752MB memory and took 6m53s to go over 10 iterations on 10 variables. My colleagues asked me if I could make some improvement so that this SAS implementation is usable in industrial applicaitons where data sets with >10K observations and hundreds of features are more than common.

In their implementation, an upper triangle matrix and a lower triangle matrix are used to obtain cumulative weighted sum from either direction of the sorted data, which increases memory consumption and calculation time in quadratic rates. But in PROC IML, there is a built-in function called CUSUM that can be used to obtain cumulative sum lightning fast. In order to use CUSUM function to replace the cumbersome matrix operation, we also need to pay attention to the dimension change. Since there is no huge triangle matrices involves, the consumption on RAM also significantly reduced, which in turn means we can process bigger data [see experiment below]. I did a test on the speed and memory consumption using AdaBoost on 2000 observation and 10 variables, with 10 iterations, before and after the improvement we see:

2000 Obs:

4000 Obs:

Consider the difference showing above, with 2000 observations, new macro used less than 0.3s and merely 1.1MB memory while the original one used 102.17s and 158.3MB memory. With 4000 observations, original macro used 752MB and 6:53s, as comparison, new macro used 0.5s and 1.9MB memory. We see that the memory and time consumption is no longer O(n^2) but rather O(n), where n is the number of observations. Both macro produced almost identical results in two experiments.

With 500K observations, the original macro was not able to proceed but the improved one finished successfully in about 70s and took up only 194MB memory:

With the much reduced processing time, on top of the boost algorithms, analysts are now able to apply more computationally intensive algorithms, such as Bagging.

I also tried to improve the REGSPLIT_IML subroutine. Yesterday, time usage was reduced 69% [1m21s vs 25s in 2K records] comparing to original macro but still increases quadratically with number of observation, memory consumption is much reduced and increase linearly with size of data. With careful study of the calculation and formula involved, I further decomposed the matrix operation into more efficient mathematical calculations, and now the time usage is further reduced to only 0.9% [329.84s vs 2.59s in 4K records] of original macro and increase only linearly with the number of records. Therefore we are practically able to use GentleBoost and LogitBoost [preferred over AdaBoost/RealBoost in many cases] for predictive modeling projects.

First Improvement

Second Improvement:

Using the most recent version of the macro, on a PC with E6320 1.86GHz (5.5Gflps/core) and 4GB memory, we observe the following performance benchmark (another 15% improvement over what's shown above):

Replace SPLIT_IML subroutine with the following code.


/*************************************************************************
This is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 or 3 of the License
(at your option).

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
**************************************************************************/
    start split_iml(x,y,w,g_info,out_type, y_pred);
        n = nrow(x);
        p = ncol(x);
        gini_min = 2;
        gini_var = 0;
        gini_cut = 0;
        y_pred = repeat(0,n,1);
        wsum = sum(w);
        ywsum = sum(y#w);  
        ywsum1 = wsum - ywsum; 
        do j=1 to p;
            x_curr = x[,j]||y||w;
            b=x_curr;  
            x_curr[rank(x[,j]),]=b;  free b;
            x_sort = x_curr[,1]; 
            y_sort = x_curr[,2]; 
            w_sort = x_curr[,3];
            yw_sort=(y_sort#w_sort);
            yw_sort1=(w_sort - yw_sort);        
            yw_cusum=cusum(yw_sort[1:(n-1)]);    

            lpwt = cusum(w_sort[1:(n-1)]);
            lpwt = lpwt#(lpwt >= 2*CONSTANT('SMALL')) + 
                   (lpwt < 2*CONSTANT('SMALL'))*2*CONSTANT('SMALL');
    
            p1_L = yw_cusum # (1/lpwt);
            gini = yw_cusum # (1-p1_L);
                       
            rpwt = wsum - lpwt; 
            rpwt = rpwt#(rpwt >= 2*CONSTANT('SMALL')) + 
                   (rpwt < 2*CONSTANT('SMALL'))*2*CONSTANT('SMALL');
    
            yw_cusum = ywsum - yw_cusum;
            p1_R = yw_cusum # (1/rpwt);
            
            gini = gini + yw_cusum # (1-p1_R);

            free lpwt  rpwt  yw_cusum  yw_sort1;

            g_min=gini[><];  g_loc=gini[>:<];

            if g_min < gini_min then do;
                gini_min=g_min;
                gini_var = j;
                gini_cut = (x_sort[g_loc] + x_sort[g_loc+1]) / 2;
                p1_RH = p1_R[g_loc];
                p0_RH = 1-p1_R[g_loc];
                p1_LH = p1_L[g_loc];
                p0_LH = 1-p1_L[g_loc];

                c_R = 0;
                if p1_RH > 0.5 then c_R = 1;
                c_L = 0;
                if p1_LH > 0.5 then c_L = 1;
            end;
        end;
        g_info = gini_var||gini_min||gini_cut||p0_LH||p1_LH||c_L||p0_RH||p1_RH||c_R;
        if out_type = 1 then 
           y_pred = (x[, gini_var] <=gini_cut)*c_L + 
                    (x[, gini_var] > gini_cut) *c_R
        ;
       
        if out_type=2 then
           y_pred[, 1] =( x[, gini_var]<=gini_cut) * ( (c_L=0)*(1-p0_LH) + (c_L=1)*p1_LH) +
                        ( x[, gini_var] >  gini_cut) * ( (c_R=0)*(1-p0_RH) + (c_R=1)*p1_RH)
        ;
 
 
    finish split_iml;

Replace REGSPLIT_IML subroutine with the following code.


/*************************************************************************
This is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 or 3 of the License
(at your option).

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
**************************************************************************/
start regsplit_iml(x,y,w,j_info,y_pred);
        n = nrow(x);
        p = ncol(x);
        min_css = 10000000000000;
        y_pred = repeat(0,n,1);
  wy2sum = sum( w#y#y );
        wsum = sum(w);
        ywsum = sum(y#w);
  ywsum1 = wsum - ywsum;
        do j=1 to p;
            x_curr = x[,j]||y||w;
            b=x_curr;
            x_curr[rank(x[,j]),]=b;   free b;
            x_sort = x_curr[,1];
            y_sort = x_curr[,2];
            w_sort = x_curr[,3];

   yw_sort=(y_sort#w_sort);
   yw_sort1=((1-y_sort)#w_sort);
   w_sort = (w_sort);

   yw_cusum = cusum(yw_sort[1:(n-1)]);

   lpwt = cusum(w_sort[1:(n-1)]);
   lpwt = lpwt# (lpwt>constant('SMALL')) + 
                         constant('SMALL')#(lpwt<=constant('SMALL'));
   p1_L = (yw_cusum # (1/lpwt));

   rpwt = wsum - lpwt;
   rpwt = rpwt#(rpwt>constant('MACEPS')) + 
                         constant('MACEPS')#(lpwt<=constant('MACEPS'));
   p1_R = ((ywsum - yw_cusum) # (1/rpwt)); 

   css=(1:n-1)*0;
   lpwt = cusum(w_sort); rpwt = cusum(yw_sort);

   css = wy2sum + p1_L##2#lpwt[1:(n-1)] + p1_R##2#(wsum - lpwt[1:(n-1)]) -
           2*(p1_L#rpwt[1:(n-1)] + p1_R#(ywsum - rpwt[1:(n-1)]));

   free  lpwt  rpwt  yw_cusum  yw_sort1;
   css_min=css[><];  css_loc=css[>:<];

            if css_min < min_css then do;
                min_css = css_min;
                cut_val = (x_sort[css_loc] + x_sort[css_loc+1]) / 2;
                reg_var = j;
                ypred_L = (sum(yw_sort[1:css_loc]))/sum(w_sort[1:css_loc]);
                ypred_R = (sum(yw_sort[css_loc+1:n]))/
                        sum(w_sort[css_loc+1:n]);
                y_pred = ypred_L*(x[,j] < cut_val) + ypred_R*(x[,j] >= cut_val);
                j_info = reg_var||min_css||cut_val||ypred_L||ypred_R;
            end;
        end;
    finish regsplit_iml;

Reference:
Dmitrienko, Alex, Christy Chuang-Stein, and Ralph D’Agostino. Pharmaceutical Statistics Using SAS®: A Practical Guide. Cary, NC: SAS Institute Inc. 2007