Monday, June 24, 2013

Merging Data -- SAS, R, and Python

On analyticbridge, the question was posed about moving an inner join from Excel (which was taking many minutes via VLOOKUP()) to some other package.  The question asked what types of performance can be expected in other systems.  Of the list given, I have varying degrees of experience in SAS, R, and Python.

In the question, the user has about 80,000 records that match between 2 tables and need to be merged by a character key.

To start, let's create 2 tables in SAS.  Each will have 150,000 records and there will be 80,000 overlapping records.  We will randomize the tables (obviously the merge is trivial when the tables are sorted with the matching records on top) to be more "real world."

data A;
format address $32. a $8.;
a = "1";
do i=1 to 150000;
       rand = ranuni(123);
       if i <=80000 then do;
              address = put(i,z32.);
       end;
       else do;
              address = "AAA" || put(i,z29.);
       end;
       output;
end;
run;

data B;
format address $32. b $8.;
b = "1";
do i=1 to 150000;
       rand = ranuni(234);
       if i <=80000 then do;
              address = put(i,z32.);
       end;
       else do;
              address = "BBB" || put(i,z29.);
       end;
       output;
end;
run;

proc sort data=a;
by rand;
run;

proc sort data=b;
by rand;
run;

Now SAS pages tables to the hard drive.  Good and bad.  Python and R will be starting with the tables in memory, so we use SASFILE to load the tables into the main memory.  Also note that SAS is writing the result table to the HD.  We will do the same in Python and R.

594  sasfile a load;
NOTE: The file WORK.A.DATA has been loaded into memory by the SASFILE statement.
595  sasfile b load;
NOTE: The file WORK.B.DATA has been loaded into memory by the SASFILE statement.
596
597  proc sql noprint;
598  create table temp.tableA_B as
599  select a.address,
600         a.a,
601         b.b
602      from a inner join
603           b
604        on a.address = b.address;
NOTE: Table TEMP.TABLEA_B created, with 80000 rows and 3 columns.

605  quit;
NOTE: PROCEDURE SQL used (Total process time):
      real time           0.09 seconds
      cpu time            0.09 seconds

So SAS took 0.09 seconds.  Much faster than the many minutes in Excel.

In R this is trivial.  I wrote the tables to csv files (not shown).  So we will read them in, do the merge, and then save the result as an .rda file.

> system.time({
+   a = read.delim('c://temp//a.csv',header=T,sep=',')
+   b = read.delim('c://temp//b.csv',header=T,sep=',')
+ })
   user  system elapsed
   9.03    0.06    9.09

> system.time({
+   m = merge(a[c("address","a")],b[c("address","b")])
+ })
   user  system elapsed
   2.15    0.02    2.17

> system.time({
+   save(m, file="c://temp//tableA_B.rda")
+ })
   user  system elapsed
   0.21    0.00    0.20 

R took 2.17 seconds for the merge and 0.20 seconds to write.  A total of 2.37 seconds.

I am least familiar with the optimal way to do this in Python.  I have a question to my Python guru about the optimal way to do the merge.  For the time being, my attempt is here.  The basics is to read the files into a Dictionary with the address string as the key and a basic object as the value.  Then iterate over the keys in 1 table and see if they are in the second table.  If so, add the merged data to a new dictionary.

Outputting the merged data might be faster in a list instead of a dictionary.  The index hash is not being built in SAS or R.

C:\Temp>c:\Python27\python.exe merge.py
2.4200000763 Starting Merge
Took 0.2660000324 seconds to merge, 1.1210000515 total with pickle

Python took 0.27 seconds for the merge and a total of 1.12 seconds for the merge and write.

A final note, here is a screen grab of the resulting file sizes.  R wins hands down -- I assume there is some compression going on there.

This example is fairly trivial.  Hopefully, someone will find it useful while trying to learn one of these languages.

Sunday, June 2, 2013

Logistic Regression in PROC MODEL

PROC MODEL in SAS/ETS provides an open interface for fitting many different types of linear and nonlinear models.  PROC MODEL allows users to simulate from fitted models and allows users to register models with SAS Risk Dimensions (PROC RISK et al).  Risk Dimensions is SAS's market and credit risk analysis package used by many large banks and trading organisations.

Being a user of both ETS and occasional user of Risk Dimensions, I've often wanted to incorporate logistic models into my simulations.  It is possible to fit the logistic regressions through numerous PROCs in SAS/STAT and incorporate the fitted parameters into MODEL or RISK statements.  It would be nice to be able to fit a logistic regression in PROC MODEL to save the macro programming hoops necessary to fit in something else and incorporate in MODEL.

First there are two resources I used for this example.

  1. Wikipedia on Logistic Regression
  2. A well written paper by Scott A. Czepiel
Now, let's generate some binomial data
data test;
a = 5;
b = 10;

do i=1 to 10000;
       rand = rannor(123);
       xb = a + b*rand ;
       p = 1 / (1+ exp(-xb));
       if ranuni(123) < p then
              y = 1;
       else
              y = 0;
       output;
end;
run;
Next, let's fit the model in PROC LOGISTIC
proc logistic data=test descending;
model y = rand;
run;
The interesting part of the output:
                                      The LOGISTIC Procedure

                            Analysis of Maximum Likelihood Estimates

                                              Standard          Wald
               Parameter    DF    Estimate       Error    Chi-Square    Pr > ChiSq

               Intercept     1      4.9155      0.1517     1050.0422        <.0001
               rand          1      9.7196      0.2896     1126.4620        <.0001
The fitted values are within reason of the true values.

Now let's use the eq 9. on page 5 of the Czepiel paper.  Also of note is the documentation on the general likelihood function in PROC MODEL for use in ML estimation.
proc model data=test;
parameters a b;
reg = a + b*rand;
y = 1 / (1 + exp(-reg));

llik = -(y*reg - log(1 + exp(reg)));

errormodel y ~ general(llik);

fit y / fiml ;
quit;
First, we define a function, "reg," which is our regression equation.  Next, we define the equation for y.  Using Eq. 9, we define the likelihood function.  Note, PROC MODEL uses a minimization of the sum of the likelihood function, so we negate the value.

The key is to now define the ERRORMODEL for y as a general likelihood function.  The FIT statement with the "/ fiml" option tell SAS to use a ML estimation.  Because the general likelihood function is defined, this is used for the fitting.

Here is the relevant output from PROC MODEL:
                                       The MODEL Procedure

                         Nonlinear Liklhood Summary of Residual Errors

                       DF       DF                                                        Adj
    Equation        Model    Error         SSE         MSE    Root MSE    R-Square       R-Sq

    y                   2     9998       359.5      0.0360      0.1896      0.8335     0.8335


                              Nonlinear Liklhood Parameter Estimates

                                                 Approx                  Approx
                   Parameter       Estimate     Std Err    t Value     Pr > |t|

                   a               4.915868      0.1517      32.42       <.0001
                   b               9.719905      0.2895      33.57       <.0001


                      Number of Observations       Statistics for System

                      Used             10000    Log Likelihood        -1179
                      Missing              0