A macro calls R in SAS for paneled 3d plotting

SAS and R could complement each other. SAS is a versatile ETL (extraction, transformation and loading) machine and its statistical procedures based on generalized linear model are impeccable. R would bring cutting-edge data mining and data visualization technologies at low cost (or no cost). Although the two packages dwell in distinctive ecosystems (for example: different OS/ETL/database/reporting layers) [Ref. 1], mixed programming by combining them together would make an analytics shop invincible.

Some SAS programmers like to use SAS/IML to call R’s functions [Ref. 2]. However, it seems that SAS/IML fails to work with the latest versions of R since 2.12 [Ref. 3]. Others tend to play tricks to call R into SAS’s data step to meet their daily needs [Ref. 4]. In this post, the macro below would call the ‘lattice’ package of R in SAS, on a PC platform, to draw paneled three dimension images, since currently SAS’s SG procedures don’t own such an option. The good thing is that there is no need to check the version of R installed before running it. And the modification of this macro can be extended to other applications to call R in SAS.

Reference:
1. ‘Keep an Eye on the emerging Open-Source Analytics Stack’. Revolution R Blog.
2. Zhengping Ma. ‘Data mining in SAS with open source software’. SAS Global 2011.
3. ‘SAS/IML incompatible with latest releases of R’. SAS-L. 11APR2011.
4. Liang Xie. ‘Regularized Discriminant Analysis’. www.sas-programming.com

/*******************READ ME*********************************************
* - A MACRO CALLS R IN SAS FOR PANELED 3D PLOTTING  -
*
* VERSION:     SAS 9.2(ts2m0), windows 64bit
* DATE:        25apr2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME*********************j********************/

****************(1) MODULE-BUILDING STEP******************;
%macro scatter3dpanel(data = , x = , y = , z = , factor = , 
                      width = , height = , outfile = );
   /***********************************************************
   *  MACRO:      scatter3dpanel()
   *  PARAMETERS: data   = dataset for plotting
   *              x      = x-axis variable
   *              y      = y-axis variable
   *              z      = z-axis variable
   *              factor = partition factor variable
   *              width  = width of output graph
   *              height = height of output graph
   *              outfile= location of output image
   ***********************************************************/
  proc export data = &data outfile = "d:\tmp.csv" replace;
  run;
  
  proc sql;
    create table _tmp0 (string char(80));
    insert into _tmp0 
    set string = 'tmp=read.csv("d:/tmp.csv", header=T)'
    set string = 'attach(tmp)'
    set string = 'library(lattice)'
    set string = 'windows()'
    set string = 'cloud(sas_zvar~sas_xvar+sas_yvar|as.factor(sas_factor), pretty=T)'
    set string = 'dev.print(device=png, width=sas_width, height=sas_height, file="sas_file")';
  quit;
  
  data _tmp1;
    set _tmp0;
    string = tranwrd(string, "sas_xvar", propcase("&x"));
    string = tranwrd(string, "sas_yvar", propcase("&y"));
    string = tranwrd(string, "sas_zvar", propcase("&z"));
    string = tranwrd(string, "sas_factor", propcase("&factor"));
    string = tranwrd(string, "sas_width", "&width");
    string = tranwrd(string, "sas_height", "&height");
    string = tranwrd(string, "sas_file", translate("&outfile", "/", "\"));
  run;

  data _null_;
    set _tmp1;
    file "d:\callRinSAS.r";
    put string;
  run;

  options noxsync noxwait;
  x ' "d:\Program Files\R\R-2.12.1\bin\R.exe" CMD BATCH --vanilla --slave "d:\callRinSAS.r" ';
%mend;

****************(2) TESTING STEP******************;
%scatter3dpanel(data = sashelp.cars, x = length, y = wheelbase, z = horsepower, 
            factor = type, width = 1200, height = 600, outfile = d:\test1.png );

****************END OF ALL CODING***************************************;

Some analysis on university ranking by US News

The yearly US News best college ranking is an important tool in comparing schools for students and their eager parents. The latest data is publicly available (paying 20 bucks would get full access) [Ref.1]. And the methodology is easy to find and explain [Ref.2]: a score would be weighted by peer assessment, retention, faculty resources, student selectivity, graduation rate, etc; therefore the final ranking would be based on the scores of a number of colleges.

It is interesting to explore and dissect the ranking process by US News. Still the dirty job of data extraction, transformation and loading occupied 90% of the working time. Data crunching was performed with logistic regression (for private/public), and selective linear regression (for score), by the nice tools from SAS/STAT. Factor analysis and partial least square regression were used to minimize the multicollinearity that is widespread in this data.

The analysis leads to two conclusions. First, the ranking is relatively qualitative instead of quantitative. The ranking heavily depends on the reputation opinion form surveying institutions’ administrators and high schools’ counselors. Other variables just modify the result. Second, the ranking favors private universities. Being a private university would add 3 points to the overall score. The best public university, UC Berkeley, is ranked as 22nd. I didn’t find any reason why it is inferior to some private universities ahead. At the data level, the public universities and private ones are distinguishable. And apparently they target different customer groups. To be fair, the US News may divide the university ranking into two subsystems: public universities and private universities, which could be more helpful in understanding the universities' standing in their sectors.

References:
1. http://colleges.usnews.rankingsandreviews.com/best-colleges/rankings/national-universities/data
2. http://collegethrive.com/college-rankings-us-news-world-report-method

****************(1)CLUSTERING STEP******************;
ods listing close;
ods output variables = _varlist;
proc contents data = uscr11;
run;

proc sort data = _varlist;
   by num;
run;

proc sql;
   select variable into: num_vars separated by ' '
   from _varlist
   where lowcase(type) = "num" and num not in (4, 5)
;quit;

proc varclus data = uscr11 summary outtree=tree;   
   var &num_vars;
run;

ods html style = harvest;
ods graphics on;
goptions htext = 4pct ftext = "Albany AMT";
axis1 order = (0.5 to 1 by 0.1);
axis2 label = none;

proc tree horizontal haxis=axis1 vaxis=axis2;
   height _propor_;
   id _label_;
run;

proc sgscatter data = uscr11;
   matrix &num_vars /ellipse=(alpha=0.25) markerattrs=(size=1);
run;

****************(2)IMPUTATION STEP******************;
proc mi data = uscr11 nimpute = 1 round = .01 
        seed = 20110425 out = _tmp0;
   monotone regpmm(donaterate = score ugrepidx gradrate retention);
   var score ugrepidx gradrate retention donaterate;
run;

proc mi data = _tmp0 nimpute = 1 round = .01 
        seed = 20110425 out = imputed;
   monotone reg(top10fresh = score ugrepidx sat25p sat75p acceptrate);
   var score ugrepidx sat25p sat75p acceptrate top10fresh;
run;

****************(3)FACTOR ANALYSIS STEP******************;
proc factor data = imputed nfactors = 3 rotate=promax
            reorder out = factorized plots=(scree);
   var &num_vars;
run;

data _tmp1;
   set factorized;
   if type = "private" then do; shape = "club"; color = "blue"; end;
   else do; shape = "diamond"; color = "red"; end;
   keep shape color factor:;
run;

proc g3d data = _tmp1;
   scatter factor2*factor3 = factor1 / color = color shape = shape;
run;

****************(4)LOGISTIC REGRESSION STEP******************;
proc logistic data = imputed plots = (roc);
   model type = &num_vars / 
         selection = stepwise slentry = 0.3 slstay = 0.3;
run;

proc pls data = imputed plot = (corrloadplot variableimportanceplot);
   model score = &num_vars;
run;

proc sql;
   select variable into: vars separated by ' '
   from _varlist
   where num in (3, 6, 7, 8, 12, 13, 14, 15, 16)
;quit; 

****************(5)VARIABLE SELECTION STEP******************;
proc glmselect data = imputed plot = (coefficientpanel aseplot);
   partition fraction(validate = 0.5);
   class type;
   model score = &vars / 
         selection = stepwise(choose = validate select = sl);
run;
ods graphics off;
ods html close;
ods listing;

****************END OF ALL CODING***************************************;

A subroutine in SAS to simulate asset pricing paths

For matrix computation in SAS, SAS/IML is the choice. This module has its own syntax, functions and even plotting subsystem. Some statisticians used it to realize the algorithms beyond the reach of SAS’s procedures, for example, boosting [Ref. 1]. However, comparing with other popular matrix-based languages, such as R and Matlab, SAS/IML has no edge. SAS’s most valuable products are still its robust data step and statistical procedures. ‘Porting’ source codes from other languages into SAS has to rely on data step.

Asset prices can be estimated by Monte Carlo simulation. To generate a series of price-evolving paths with several most common parameters, Dr. Brandimarte codes a naive Matlab function to apply the standard Wiener process [Ref. 2]. With random seeds for normal distribution, multiple pricing mechanisms can be demonstrated and compared. In SAS 9.2, the workflow of simulation and visualization would be modularized as a data step subroutine. Later, such a pricing subroutine could be easily invoked under given circumstances.

It is said that SAS 9.3 is going to be released 2011Q3. Hope this time, the data step function compiler, Proc Fcmp, could be more dynamic and with more methods.

References:
1. Dmitrienko, Alex, Christy Chuang-Stein, and Ralph D’Agostino. Pharmaceutical Statistics Using SAS: A Practical Guide. SAS Publishing. 2007
2. Paolo Brandimarte. Numerical methods in finance. John Wiley & Sons. 2002.

/*******************READ ME*********************************************
* - A SUBROUTINE IN SAS TO SIMULATE ASSET PRICING PATHS -
*
* VERSION:     SAS 9.2(ts2m0), windows 64bit
* DATE:        17apr2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME*****************************************/

****************(1) MODULE-BUILDING STEP******************;
******(1.1) COMPILE FUNCTION-ACCOMPANYING MACRO*************;
option mstored sasmstore = sasuser;
%macro AssetPath_macro() / store source;
   data _tmp1;
      set _tmp;
      day  = _n_ - 1;
   run;

   proc transpose data = _tmp1 out = _tmp2 ;
      by day;
      var path:;
   run;

   data _tmp2;
      set _tmp2;
      label _name_ = 'Simulated paths';
   run;

   ods html style = money;
   proc sgplot data = _tmp2;
      series x = day y = col1 / group = _name_;
      yaxis grid label = 'Asset price';
      xaxis grid label = 'Change by days';
   run;
   ods html close;
%mend;

******(1.2) COMPILE SUBROUTINE FOR ASSET PRICING*************;
proc fcmp outlib = sasuser.astpth.funcs;
   subroutine AssetPath(S0,  mu, sigma, T,  NSteps,  NRepl);
   /*************************************************************
   * SUBROUTINE:    AssetPath()
   * PARAMETERS:    S0    = the initial price
   *                mu    = the drift
   *                sigma = the volatility
   *                T     = the horizontal time
   *                NSteps= the number of time steps
   *                NRepl = the number of replications
   *************************************************************/
      array SPaths[1, 1] / nosymbols;
      array Path[1, 1] / nosymbols;
      call dynamic_array(SPaths, NRepl, NSteps + 1);
      call dynamic_array(Path, NSteps + 1, NRepl);
      call zeromatrix(SPaths);
      do _row = 1 to NRepl;
         SPaths[_row, 1] = S0;
      end;
      dt = T / NSteps;
      nudt = (mu - 0.5*sigma**2) * dt;
      sidt = sigma * sqrt(dt);
      do _row = 1 to NRepl;
         do _col = 1 to Nsteps;
            SPaths[_row, _col + 1] = SPaths[_row, _col] * 
               exp(nudt + sidt*rannor(0));
         end;
      end;   
      call transpose(SPaths, Path);
      rc1 = write_array('_tmp', Path);
      rc2 = run_macro('AssetPath_macro');
    endsub;
quit;
****************END OF STEP (1)******************;

****************(2) TESTING STEP******************;
option cmplib = (sasuser.astpth) mstored sasmstore = sasuser;
data _null_;
   S0 = 50; mu = 0.1; sigma = 0.3; T = 1; NSteps = 365; NRepl = 3;
   call AssetPath(S0,  mu, sigma, T,  NSteps,  NRepl);
run;
****************END OF STEP (2)******************;

****************END OF ALL CODING***************************************;

Predict unemployment rate for Election 2012 by SAS

Since recently President Obama announced that he is seeking reelection, the unemployment rate on November 2012 would decide the result. The Wall Street Journal averaged 54 economists’ predication and concluded that the number is going to be 7.7%. Apparently, those economists rely on the historical data to forecast the future, together with more or less their subjective judgment. However, the newly released March data is surprisingly good: 8.8%, which means that this predication number has to be adjusted downwardly to be below 7.7%. Then what is the real-time prediction of the unemployment rate for this ‘big’ time?

SAS has one of the finest time-series packages in the world: SAS/ETS which includes a few predictive procedures such as the ARIMA procedure and the FORECAST procedure[Ref. 2]. And the economic data is updated by Federal Reserve and well accessible on their website. To predict unemployment rate like a real professional is possible with a notebook computer and SAS. Of course SAS’s procedures have tons of methods and parameters to tune. To simply this problem, in the SAS macro below, I chose a conservative method and an aggressive one, to give a rough estimation about the unemployment range. Just like what the WSJ said, the trend matters. The predication will be more approaching to the real number as time goes forward. Right now, my prediction for the unemployment rate on November 2012 is from 7.1% to 7.4%.

References:
1."Jobless Rate at 2012 Presidential Vote Forecast at 7.7%, Highest Since Carter-Ford, but the Trend May Matter Most". The Wall Street Journal, 13Mar2011.
2.SAS/ETS 9.2 User Guide. SAS Publishing, 2008.

/*******************READ ME*********************************************
* -- PREDICT UNEMPLOYMENT FOR ELECTION 2012 LIKE A PRO --
*
* VERSION:     SAS 9.2(ts2m0), windows 64bit
* DATE:        09apr2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME*****************************************/

****************(1) MODULE-BUILDING STEP******************;
%macro unemrate(predtime = );
   /***********************************************************
   *  MACRO:      unemrate()
   *  GOAL:       use time series based on latest FED data to
   *              predict unemployment rate in US and plot
   *  PARAMETERS: predtime = the time when unemployement rate 
   *                         is to be predicted
   *
   ***********************************************************/
   filename _infile url 
      "http://research.stlouisfed.org/fred2/data/UNRATE.txt" 
      debug lrecl=100;

   data raw;
      infile _infile missover firstobs = 22;
      format date date9.;
      input @1 date yymmdd10. @13 value 4.1;
   run;

   data _null_;
      set raw end = eof;
      if eof then do;
         interval = intck('month', date, input("&predtime", monyy7.));
         call symput('interval', interval);
         call symput('eodate', date);
         call symput('insert', 'Lastest data:' || put(value, 4.1) || 
                   '% on ' || put(date, monyy7.));
      end;
   run;

   %if %eval(&interval) le 0 %then %do; 
      %put ERROR: Predicted time must be greater than latest time FED posts data; 
      %goto finish; 
   %end;

   ods select none;
   proc forecast data = raw out = _predbyfc outfull
      method = stepar lead = &interval interval = month;
      id date;
      var value;
   run;

   proc arima data = raw;
      identify var = value;
      estimate p = 1 q = 12;
      forecast lead = &interval interval = month 
              id = date out = _predbyar;
   quit;
   ods select all; 

   proc sql;
      create table predicted0 as
      select a.date, a.value label = 'Real unemployment rate', 
            a.forecast as predbyar label = 'ARIMA model', 
            b.value as predbyfc label = 'STEPAR model'
      from _predbyar as a,
           _predbyfc (where = (lowcase(_type_) = 'forecast')) as b
      where a.date = b.date
   ;quit;

   data predicted1;
      set predicted0 end = eof;
      if date lt &eodate then call missing(predbyar, predbyfc);
      else if date eq &eodate then do;
         predbyar = value; 
         predbyfc = value;
      end;
      if eof then do;
         call symput('arlast', put(predbyar, 4.2));
         call symput('fclast', put(predbyfc, 4.2));
      end;
   run;
    
   ods html style = harvest;
   proc sgplot data = predicted1;
      where date ge '01jan2006'd;
      series x = date y = value;
      series x = date y = predbyar;
      series x = date y = predbyfc;
      refline &arlast / axis = y labelloc = inside
         label = "&arlast" transparency = 1;
      refline &fclast / axis = y labelloc = inside
         label = "&fclast" transparency = 1;
      xaxis grid label = ' ';
      yaxis grid label = 'Unemployment percentage %'
         values = (4 to 11 by 0.2);
      inset "Prediction ends on &predtime" / position = topright border;
      inset "&insert" / position = bottomright;
   run;
   ods html close;

   proc datasets nolist;
      delete _:;
   quit;

   %finish: ; 
%mend;

****************(2) TESTING STEP******************;
%unemrate(predtime = NOV2012);

****************END OF ALL CODING***************************************;

Using SAS to find the best k for k-Nearest Neighbor classification

Optimize many-to-one mapping by user-defined functions

In many occasions, fast access into a lookup table to find desired value is necessary. In computer science, linked list, associative array, and hash table are widely used to construct the relationship between values and keys. Hash function, like value <-- index = Function(key), is essential to build such a hash table. Improving the hash function’s performance is pretty challenging and rewarding [Ref. 1]. In SAS, macro may be utilized to substitute function. However, macro would be failed in front of some cases, such as f(x1) + g(x2) or f(g(x)). Function or functional programming is still a better choice. With the user-defined function complier, Proc Fcmp, building reliable functions in SAS to map many keys to a single value looks promising.

However, SAS’s hash object is not applicable for coding interactive or reusable hash functions. The concept of hash object is introduced since SAS version 9, and it is ‘the first truly runtime dynamic, memory-resident DATA step structure’ [Ref. 2]. Evidence shows that it is robust, and most importantly, faster than other hard-disk based lookup solutions. And this technology is vigorously growing: SAS’s hash object has more methods, and can even accept duplicate keys [Ref. 3]. The biggest problem for SAS’s hash object is that it is transient and not savable. Hash object has to be declared within a data step to load data. If there is any other query, it has to be loaded again. If with frequently invocation, the loading factor will be formidable. As the result, SAS’s hash object is particularly useful for batch processing of lookup tasks, say, finding many key-value pairs simultaneously. That is probably why SAS programmers like to use hash object in merging datasets.

Thus, mapping many keys to one value by user-defined functions has to rely on other methods. Senior SAS programmers may prefer SAS arrays with help of the POINT option. However, SAS’s array is still transient. SQL and Proc Format are the options available. SQL is interactive and does not rely on any particular data type. Format is a unique data type in SAS created by Proc Format, which is permanent and exchangeable with a common SAS dataset. To test the efficiency of many-to-one mapping by user-defined functions, first a dataset of a million records with 2 keys and 1 value was simulated, and three functions were defined and complied by Proc Fcmp. Then those functions were repeatedly called for 10000 times. For SQL-based function, invoking each time is equivalent to querying the raw dataset once. On my notebook, the total time cost of running SQL-based function is more than 3 minutes, which is intolerable. Adding indexes to both keys of the lookup table will largely decrease the time expense to 1/6 of the original time. In addition, loading the lookup table into the memory, by the SASFILE statement, would improve the efficiency of the function a little further. Since Proc Format only allows one key, the two keys need to be concatenated as a composite key before finalizing the lookup format. Amazingly, calling of Format-based function 10000 times only spends less than 2 seconds. As the result, the Format-based function is the champion of this test. Understandably, the memory consumption is proportional to the functions’ efficiency: faster speed means more memory requirement. Much more than other methods, the Format-based function used 200 MB memories. I guess that the ‘format’ data type may be loaded into memory and stayed there during the processing time, which slashed the loading factor. In conclusion, Proc Format produced ‘format’ is not only an alternative solution to merge large datasets, but also is a reliable foundation to define key-value pairs by many-to-one mapping functions.

References:
1. Hash table on Wikimedia. http://en.wikipedia.org/wiki/Hash_table
2. Elena Muriel. Hashing Performance Time with Hash Tables. SAS Global 2007.
3. Robert Ray and Jason Secosky. Better Hashing in SAS 9.2. SAS Global 2008.

/*******************READ ME****************************************
* --- OPTIMIZE MANY-TO-ONE MAPPING BY USER-DEFINED FUNCTIONS ----
*
* HARDWARE:    a notebook with amd64 cpu 2.0g, 3g ram 
* SAS:         SAS 9.2(TS2M0), pc 64bit
* 
* TEST PASSED: 31mar2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME************************************/

************(1) GENERATE TEST DATASET*********************;
******(1.1) SIMULATE A DATASET OF 1M RECORDS WITH 2 KEYS*********;
data raw;
  do key1 =1 to 1000;
    do key2 = 1 to 1000;
      value = ranuni(20100322) + rannor(20100322);
      output;
    end;
  end;
run;

******(1.2) DERIVE INDEXED DATASET*********;
proc sql;
  create table indexed as select * from raw;
  create index key1 on indexed(key1);
  create index key2 on indexed(key2);
quit;

******(1.3) TRANSFORM DATASET INTO FORMAT*********;
****(1.3.1) CAST VALUE FROM NUMBER TO CHARACTER*****;
proc format;
  picture key low-high = '0000'(fill = '0');
run;
****(1.3.2) COMBINE 2 KEYS TO 1 COMPOSITE KEY******;
data keycombined;
  set raw;
  length keystr $8;
  keystr = cats(put(key1, key.), put(key2, key.));
  drop key1 key2;
run;
****(1.3.3) DERIVE FORMAT DATASET*****;
data fmtds;
  set keycombined;
  rename keystr = start 
          value = label;
  fmtname = '$myfmt';
  type = 'C';
run;
****(1.3.4) INCORPORATE FORMAT DATASET TO BUILD FORMAT *****;
proc format cntlin = fmtds;
run;

**************END OF STEP (1)*****************;

**************(2) CREATE USER-DEFINED FUNCTIONS***********;
*******(2.1) BUILD THE FIRST FUNCTION************;
****(2.1.1) THE EMBEDDED MACRO************; 
option mstored sasmstore = sasuser;
%macro myfunc1_macro / store source;
  %let key1 = %sysfunc(dequote(&key1));
  %let key2 = %sysfunc(dequote(&key2));
   proc sql noprint;
      select value into: value
      from raw
      where key1 = &key1
          and key2 = &key2
  ;quit;
%mend;
****(2.1.2) CREATE THE FUNCTION AND OUTPUT***********;
proc fcmp outlib = sasuser.keyvalue.funcs;
   function myfunc1(key1, key2);
      rc = run_macro('myfunc1_macro', key1, key2, value);
      if rc eq 0 then return(value);
      else return(.);
   endsub;
run;

*******(2.2) BUILD THE SECOND FUNCTION************;
****(2.2.1) THE EMBEDDED MACRO************; 
option mstored sasmstore = sasuser;
%macro myfunc2_macro / store source;
  %let key1 = %sysfunc(dequote(&key1));
  %let key2 = %sysfunc(dequote(&key2));
   proc sql noprint;
      select value into: value
      from indexed
      where key1 = &key1
          and key2 = &key2
  ;quit;
%mend;
****(2.2.2) CREATE THE FUNCTION AND OUTPUT***********;
proc fcmp outlib = sasuser.keyvalue.funcs;
   function myfunc2(key1, key2);
      rc = run_macro('myfunc2_macro', key1, key2, value);
      if rc eq 0 then return(value);
      else return(.);
   endsub;
run;

*******(2.3) BUILD THE THIRD FUNCTION************;
proc fcmp outlib = sasuser.keyvalue.funcs;
  function myfunc3(key1, key2);
    key = put(cats(put(key1, key.), put(key2, key.)), $8.);
    value = put(key, $myfmt.);
    return(value);
   endsub;
run;

**************END OF STEP (2)*****************;

**************(3) TEST USER-DEFINED FUNCTIONS***********;
****(3.1) CREATE A TEST MACRO TO RUN FUNCTIONS 10000 TIMES******;
%macro test(num);
  option cmplib = (sasuser.keyvalue) mstored sasmstore = sasuser;
  data test;
     do x = 101 to 200;
        do y = 301 to 400;  
           z = myfunc&num(x, y);
           output;
        end;
     end;
  run;
%mend;

****(3.2) TEST THE FIRST FUNCTION********;
%test(1);

****(3.3) TEST THE SECOND FUNCTION********;
%test(2);

****(3.4) TEST THE SECOND FUNCTION WITH IN-MEMORY LOOKUP TABLE********;
sasfile indexed open;
%test(2);
sasfile indexed close;

****(3.5) TEST THE THIRD FUNCTION********;
%test(3);

**************END OF STEP (3)*****************;

*************************END OF ALL CODING***************************;

Array 2.0: matrix-friendly array in Proc Fcmp

Array is probably the only number-indexed data type in SAS. Interestingly, SAS and R, the major statistical softwares, use 1 instead of 0 to specify the first element. SAS programmers adopt array mostly for multiple-variable batch-processing. For example, longitudinal summation can be achieved by specifying a one-dimensional array and then adding all array elements together. On the contrary, general-purpose programming languages, like Java, Python, C++,  widely use array or list to store and index data. The programmers who newly switch into SAS often feel confused about the limitation of SAS’ array, which is transient and need output to the parental data step to save result. Since data step does not memorize the array’s original structure, multiple dimensional arrays are almost useless in creating more complex data structures. In SAS, those multiple dimensional arrays are occasionally implemented to reshape data [Ref. 1]. If a matrix, such as a big one with random numbers of 5000 rows by 5000 columns [Example 1], is needed, Other programming languages would use int[][] to declare a matrix. In SAS, 2D array, such as a [5000, 5000] array, may be intuitively chosen to perform this task. However, the result is not desired, and the right way is still to apply a one-dimensional array. 

Array in Proc Fcmp is an upgraded array data type other than data step array. This new version of array allows not only creating matrices but also indexing matrix. In addition, the communication between data step and Proc Fcmp is fast and convenient: this procedure supplies the READ_ARRAY() and WRITE_ARRAY() functions, which transform a dataset to an array and vice versa. For example, to test an interviewee’s knowledge on SAS’s data step array, a typical interview question is to ask her/him to write a 9X9 multiplication table [Example 2]. The expected solution is to place the OUTPUT statement between the inner Do-loop and outer Do-loop. Placing OUTPUT into the inner layer of Do-loops builds an 81*9 monster, while ignoring OUPTUT would only generate the last row of multiplication table. This task is much simpler in Proc Fcmp: just produce a matrix and write it to a dataset. Proc Fcmp is a nice alternative to SAS/IML, a specialized matrices language module in SAS. For instance, a typical SAS/IML kind of job, such as filling missing cells in a matrix (or a dataset) with its diagonal elements, can be fulfilled by arrays in Proc Fcmp [Example 3]. Proc Fcmp is shipped in SAS/BASE, which means that no extra out-of-pocket money is needed for another license. Another concern is the efficiency of SAS/IML module. It is reported that frequently calling of SAS/IML's functions would decrease system speed dramatically[Ref. 3].

The matrix feature of Proc Fcmp’s array benefits other SAS programming areas. For example, Proc Transpose and data step array usually reshape data structure from either long to wide or wide to long. However, in many cases that positions between observation and variable in a dataset have to be exchanged, the two methods may require several try-and-error steps. The TRANSPOSE() subroutine in Proc Fcmp solves such a problem easily [Example 4]. Currently there are 13 functions or subroutines available in Proc Fcmp for matrices operation[Ref. 2]. Some may complain they are still not enough for their particular need. Don’t forget: Proc Fcmp makes user-defined functions and subroutines! For example, to simulate set() function in Python, a deldup_array() function, based on encapsulating Proc SQL in a macro by RUN_MACRO(), can delete duplicate elements in array[Example 5]. Therefore, users of Proc Fcmp’s array can always construct and accumulate their tools to suit their purpose.

References: 1. UCLA Statistics Course. http://www.ats.ucla.edu/stat/sas/library/multidimensional_arrays.htm
2. SAS 9.2 Online Help. http://support.sas.com/documentation/cdl/en/proc/61895/HTML/default/viewer.htm#a003193719.htm
3. Wang, Songfeng; Zhang, Jiajia. Developing User-Defined Functions in SAS®: A Summary and Comparison. SAS Global 2011.

******(1) EXAMPLE 1: CREATE A 5000X5000 RANDOM MATRIX****;
****(1.1) WRONG 2D ARRAY SOLUTION******;
data matrix1;
   array v[5000, 5000];
   do i = 1 to 5000;
     do j = 1 to 5000;
       v[i, j] = ranuni(0);
     end;
     output;
   end;
run;  

****(1.2) CORRECT 1D ARRAY SOLUTION*****;
data matrix2;
   array x[5000];
   do i = 1 to 5000;
     do j = 1 to 5000;
       x[j] = ranuni(0);
     end;
     output;
   end;
run;

******(2) EXAMPLE 2: CREATE A MULTIPLICATION TABLE******;
****(2.1) DATA STEP ARRAY SOLUTION ********;
data mt1;
   array a[9] a1-a9;   
   do row = 1 to 9;           
      do col= 1 to 9;    
         if row ge col then a[col]=row*col;  
      end;
     output;        
   end;
   drop row col;
run;

****(2.2) PROC FCMP ARRAY SOLUTION*****;
proc fcmp;
   array a[9, 9] / nosymbols;
   do row =1 to 9;
      do col = 1 to 9;
         if row ge col  then a[row,  col] = row*col;
      end;
   end;
   rc1 = write_array('mt2', a);
quit;

****(3) EXAMPLE 3: FILL MISSING CELL WITH DIAGONAL ELEMENT******;
****(3.0) INPUT RAW DATA****;
data have;
   input x1-x4;
   datalines;
   1 . . .
   2 1 . . 
   3 4 1 .
   7 6 5 1
;
run;

****(3.1) PROC FCMP TRANSPOSITION******;
proc fcmp;
   array a[4, 4] / nosymbols;
   rc1 = read_array('have', a);
   do i = 1 to 4;
      do j = 1 to 4;
         if missing(a[i, j]) = 1 then a[i, j] = a[j, i];
      end;
   end;
   rc2 = write_array('want', a);
quit;

*****(4) EXAMPLE 4: RESHAPE SQUARE-SHAPE DATA********;
****(4.0) INPUT RAW DATA********;
data have1;
   input x1-x5;
   cards;
   1 . 0 1 1
   0 1 . 0 0
   .  . 1 1 1
   0 0 0 1 .
   . 0 0 1 1
;
run;

****(4.1) PROC TRANSPOSE SOLUTION******;
data trps1;
   set have1;
   obs = _n_;
run;

proc transpose data = trps1 out = trps2;
   by obs;
   var x1-x5;
run;

proc sort data = trps2 out = trps3;
   by _name_;
run;

proc transpose data = trps3 out = want1_1;
   by _name_;
   var col1;
run;

****(4.2) PROC FCMP SOLUTION********;
proc fcmp;
   array a[5, 5] / nosymbols;
   rc1 = read_array('have1', a);
   array b[5, 5] ;
   call transpose(a, b);
   rc2 = write_array('want1_2', b);
quit;

******(5) EXAMPLE 5: FCMP DEDUPLICATION FUNCTION FOR FCMP ARRAY*******;
****(5.1) ENCAPSULATE DEDUPLICATIONA UTILITY OF PROC SQL IN A MACRO ******;
%macro deldup_array;
   %let dsname = %sysfunc(dequote(&dsname));
   %let arrayname = %sysfunc(dequote(&arrayname));
   /*(5.1.1) GENERATE UNIQUE DATASET AND ITS OBSERVATION NUMBER*/
   proc sql noprint;
      select count(unique(&arrayname.1)) into: obs_num
      from &dsname;
      create table _temp as
      select distinct *
      from &dsname;
   quit;
   /*(5.1.2) USE TEMP DATASET TO REPLACE RAW DATASET*/
   data &dsname;
      set _temp;
   run;
   /*(5.1.3) DELETE TEMP DATASET*/ 
   proc datasets;
      delete _temp;
   run;
%mend deldup_array; 

****(5.2) ENCAPSULATE MACRO ABOVE IN A FUNCTION*****;
proc fcmp outlib=work.func.practice;
   function deldup_array(dsname $, arrayname $);
      rc = run_macro('deldup_array', dsname, arrayname, obs_num);
      if rc eq 0 then return(obs_num);
      else return(.);
   endsub;
run;

****(5.3) USE THIS FUNCION TO DELETE DUPLICATES IN ARRAY*****;
option cmplib = (work.func) mlogic mprint symbolgen; 
proc fcmp;
   array a[1000] /nosymbols;
   do j = 1 to 1000;
      a[j] = ceil((ranuni(12345)*100) + rannor(12345));
   end;
   
   dsname = 'numbers'; 
   rc1 = write_array(dsname, a);
   
   n = deldup_array(dsname,  %sysfunc(quote(a)));
   
   call dynamic_array(a, n);
   rc2 = read_array(dsname, a);
   put a = ;
quit;

*****************END OF ALL CODING****************************;