Cluster analysis on a pivot table

The link of the pivot table is here

The increasing supremacy of JavaScript on both server side and client side seems a good news for those statistical workers who deal with data and model, and therefore always live in the darkness. They could eventually find a relatively easier way to show off their hard work on Web, the final destination of data. Here I show how to display the result of a cluster analysis on a web-based pivot table.

Back-end: cluster analysis

SAS has a FASTCLUS procedure, which implements a nearest centroid sorting algorithm and is similar to k-means. It has some time and space advantages over other more complicated clustering algorithms in SAS.

I still use the SASHELP.CLASS dataset and cluster the rows by weight and height. I specify 2 clusters and easily obtain the distances to the centroids by PROC FASTCLUS. The plot demonstrates thatweight=100 looks like the boundary to separate the two clusters. Next in DATA Step, I translate the SAS dataset to JSON format so that the browser can understand it.

************(1) Cluster the dataset*******;
proc fastclus data = sashelp.class maxclusters = 2 out = class;
   var height weight;
run;

proc sgplot data = class;
   scatter x = height y = weight /group = cluster;
   yaxis grid;
run;

************(2) Transform to JSON*********;
data toJSON;
   set class;
   length line $200.;
   array a[5] _numeric_; 
   array _a[5] $20.;
   do i = 1 to 5;
      _a[i] = cat('"',vname(a[i]),'":', a[i], ',');
   end;
   array b[2] name sex;
   array _b[2] $20.;
   do j = 1 to 2;
      _b[j] =cat('"',vname(b[j]),'":"', b[j], '",');
   end;
   line = cats('{', cats(of _:), '},');
   substr(line, length(line)-2, 1) = ' ';
   keep line;
run;

Front-end: pivot table

Pivot table is a nice way to present data, especially raw data. There are a few approaches to realize pivot table on web, such as Google's fusion table. Nicolas Kruchten developed a framework called PivotTable.js on github, which is very popular.

I embed the JSON data with the PivotTable.js to make the HTML file static, since the Blogger doesn't provide the function of HTTP server. The file content will be like:

Eventually we can view the cluster result on a pivot table. The audience can now interactively play with data.

The years to get green card for Indian and Chinese

The path toward a green card is especially difficult for Indian and Chinese who are working in the US and named as EB2 workers, since they have to wait years to submit an i485 form after being approved with their i140 form submission. Once the cutoff date set by USCIS each month catches the priority date that is a date of filling labor certification and decides the position in the long line, the applicants holding their priority date will be able to file the i485 form. However, the cutoff dates are unpredictable for the public, and “it is impossible to accurately estimate how long that may take”.

Situation
The recent cut-off dates data can be copied and pasted from Wiki. Then I adjusted the cut-off dates within the retrogression period. In the past decade, the waiting years range from about 0 year to 4.25 years. July 2007 looks like an amnesty for those who have priory date before that, otherwise people have to wait at least 1.85 years. Applicants from the adjacent 3 or 4 years usually wait in line for the door to be open. The door is recently shut off on June 2012 and hopefully will be open again one day in 2014.

data pd;
   input @1 _cmon $3. @5 _cyear $4. @13 _pdmon $3. @17 _pddate : @21 _pdyear $2.;
   format c_time pd_time date9.;
   c_time = input(cats('01', _cmon, _cyear), date9.);
   pd_time = input(cats(_pddate, _pdmon,_pdyear), date9.);
   dif_day = c_time - pd_time;
   dif_year = dif_day/365;
   drop _:;
cards;
May 2012   Aug 15  07   
Apr 2012   May 1   10   
Mar 2012   May 1   10   
Feb 2012   Jan 1   10   
Jan 2012   Jan 1   09   
Dec 2011   Mar 15  08
Nov 2011   Nov 1   07   
Oct 2011   Jul 15  07
Sep 2011   Apr 15  07
Aug 2011   Apr 15  07
Jul 2011   Mar 8   07   
Jun 2011   Oct 15  06
May 2011   Jul 1   06   
Apr 2011   May 8   06   
Mar 2011   May 8   06   
Feb 2011   May 8   06   
Jan 2011   May 8   06   
Dec 2010   May 8   06   
Nov 2010   May 8   06   
Oct 2010   May 8   06   
Sep 2010   May 8   06   
Aug 2010   Mar 1   06   
Jul 2010   Oct 1   05   
Jun 2010   Feb 1   05   
May 2010   Feb 1   05   
Apr 2010   Feb 1   05   
Mar 2010   Feb 1   05   
Feb 2010   Jan 22  05
Jan 2010   Jan 22  05
Dec 2009   Jan 22  05
Nov 2009   Jan 22  05
Oct 2009   Jan 22  05
Sep 2009   Jan 8   05   
Aug 2009   Oct 1   03   
Jul 2009   Jan 1   00   
Jun 2009   Jan 1   00   
May 2009   Feb 15  04
Apr 2009   Feb 15  04
Mar 2009   Feb 15  04
Feb 2009   Jan 1   04   
Jan 2009   Jul 1   03   
Dec 2008   Jun 1   03   
Nov 2008   Jun 1   03   
Oct 2008   Apr 1   03   
Sep 2008   Aug 1   06   
Aug 2008   Jun 1   06   
Jul 2008   Apr 1   04   
Jun 2008   Apr 1   04   
May 2008   Jan 1   04   
Jul 2007   Jul 1   07
Jun 2007   Apr 1   04   
May 2007   Jan 8   03   
Apr 2007   Jan 8   03      
Mar 2007   Jan 8   03      
Feb 2007   Jan 8   03      
Jan 2007   Jan 8   03      
Dec 2006   Jan 8   03      
Nov 2006   Jan 8   03   
;;;
run;

proc sql;
   create table pd1 as
   select a.pd_time 'priority date', a.c_time, (select min(c_time) from pd as b 
      where a.pd_time le b.c_time and a.pd_time le b.pd_time) as i140_submit_date format date9.,
      (calculated i140_submit_date - a.pd_time) / 365 as i140_waiting_year 
   from pd as a 
   where year(pd_time) gt 2002
   order by a.pd_time
;quit;

proc sgplot data = pd1;
   series x = pd_time y = i140_submit_date;
   xaxis grid;
run;

proc sgplot data = pd1;
   series x = pd_time y = i140_waiting_year;
   yaxis grid;
run;

Possible change
The cut-off dates much depend on the changes of law and policy which reflects the economic environment. Therefore the US unemployment rate is possibly helpful in predicting the fluctuation of the cut-off dates. First I imported such data, and then transformed the difference between the priority date and the cut-off date to 6-month average. Then the two curves match well. I can even fit the relationship with a simple linear regression. As a conclusion, one percent of unemployment rate decrease may shorten 1/3 year of the waiting time from i140 to i485.

filename _infile url "http://research.stlouisfed.org/fred2/data/UNRATE.txt" debug lrecl=100;
data unempl;
   infile _infile missover firstobs = 22;
   format date date9.;
   input @1 date yymmdd10. @13 unemployement_rate 4.1;
run;

proc sql;
   create table combine as
   select a.dif_year, b.*
   from pd as a, unempl as b
   where a.c_time = b.date
;quit;

proc expand data = combine out = combine1 method=none;
   convert dif_year = ma_dif_year / transform = (movave 6);
run;

proc sgplot data = combine1;
   series x = date y = ma_dif_year;
   series x = date y = unemployement_rate;
   xaxis grid; yaxis grid label = ' ';
run;

proc reg data = combine1;
   model ma_dif_year = unemployement_rate;
   ods select ParameterEstimates FitPlot;
run;

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***************************************;

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***************************************;

Effectiveness of two mail list groups: SAS-L and R-help

Software’s strength depends on the cohesion of the community backing it. Though a commercial package comes with technique support guarantee, the speed and efficiency of telephone wired customer service may not suit the fast-evolving programming need. Especially for a statistical package, such as SAS and R, which typically deals with many small extracting, loading, transformation and analysis tasks, quick short answer to a tricky question is desired. Community based mail list is a fast approach to get question posted and solved. With the help of Google’s Gmail, huge volume of emails generated by such mail lists can be collected and sorted conveniently. As for me, R-help and SAS-L are probably the most popular user groups for R and SAS, separately. And as a learner, I constantly gain kind help from those SAS or R gurus in the two communities. To compare the two user groups, I gathered threads from Gmail up-to-March 7th this year and parsed them to digestible data.

R’s ever-increasing popularity is reflected by enormous number of topics to be discussed. Averagely 37.4 questions are posted on R-help every day, compared with meagerly 8.8 a day on SAS-L. R-help Users tend to discuss a wide range of questions, most likely focusing on modeling and visualization on a specific package, while SAS-L users are mainly interested in data integration and management, such as Proc SQL, data step and macro. The usage deviation may demonstrate the distinctive fields for SAS and R in daily practice. Most SAS-L users are more familiar with the background involved for the questions others ask: a question usually gets 5.0 follow-ups. On the contrary, R-help users expect 2.8 follow-ups for the question they posted, and the chance of no-reply-at-all is also high. In SAS-L, information concentrates on several senior SAS experts who are experienced and ready to provide clues. R-help users are more diverse, possibly because that R has more than 3000 packages and is difficult for an R user to be an all-around player.

SAS and R also have other active mail lists and forum for technical discussion. For those who love SAS-L and R-help mailing lists and wish to become a two-way statistical programmer, getting to know some aspects of the two user groups could benefit us in using them effectively.

****(0) EXPORT THE HEADS FROM MY GMAIL TO FLAT TEXT FILES BY DIFFERENT GROUP TAGS****;
********(1) PREPARE A MACRO TO EXTRACT KEY WORDS FROM TEXTS***************;
%macro extract(group);
   /*(1.1) INTEGRATE LINES FROM RAW TEXT*/
   data &group.ug;
      infile "H:\Regular_expression\&group._group.txt" truncover dsd  lrecl=400;
       input string $400.;
   run;
   /*(1.2) REMOVE GMAIL-TAG-CAUSING REDUNDANT LINES AND INDEX EACH LINES BY THREAD*/
   data &group.ug_c;
      set &group.ug;
      where string is not missing
              and string not in ('LinkedIn', 'R-Group', 'Inbox', 'SAS(noreply)', 'SAS-L');
      thread = ceil(_n_/3);
   run;
   /*(1.3) TRANSPOSE TOPIC|POST AUTHORS|TIME*/
   proc transpose data=&group.ug_c out=&group.ug_s(rename=(col1 = writer col2 = topic col3 = time));
      by thread;
      var string;
   run;
   /*(1.4) PARSE TARGET INFORMATION FOR TOPIC|POST AUTHORS|TIME*/
   data &group.ug1;
      set &group.ug_s(drop=thread _name_);
      /*(1.4.1) PARSE POST AUTHOR   */
       position = prxmatch("/\(\d+/",writer);
       if position ne  0 then  
            rep_num = input(compress(substr(writer, position+1, 2), ')'), 2.) ;
       else rep_num = 0; 
      /*(1.4.2) PARSE THREAD TOPIC*/
      idx = index(topic, '? -');
       topic_str= input(substr(topic, 1, idx), $200.);
      /*(1.4.3) PARSE TIME*/
      length time1 $6.;
      if sum(index(time, 'pm'), index(time, 'am')) > 0 then time1 = 'Mar 7';
      else time1 = strip(time);
      date = input(cats(substr(time1, 5, 2), substr(time1, 1, 3), '2011'), date9.);
      format date date9.;
      drop idx  writer topic position time:;
   run;
   /*(1.5)* SUMMARIZE BY EACH DATE*/
   proc sql;
      create table &group.ug2 as
      select mean(rep_num) as &group._mean_reply format=3.1,
            sum(ifn(rep_num=0, 1, 0)) as &group._zero_reply,
            max(rep_num) as &group._max_reply, count(topic_str) as &group._topicnum, 
            calculated &group._zero_reply / calculated &group._topicnum   
                  as &group._zero_reply_ratio format=3.1,
            date
      from &group.ug1
      group by date
   ;quit;
%mend extract;

**********(2) EXECUATE MACRO TO BUILD DATASETS FOR 2 GROUPS*****;
option mprint symbolgen;
%extract(group=sas);
%extract(group=R);

*********(3) CONCATENATE DATA AND REPORT****;
****(3.1) CONCATENATE 2 DATASETS FROM SAS AND R******; 
proc sql;
   create table reportdata as
   select mean(rep_num) as avg_rep_num, count(topic_str) as topicnum, date, 'SAS' as ug length=3
   from SASug1
   group by date
   union 
   select mean(rep_num) as avg_rep_num, count(topic_str) as topicnum, date, 'R' as ug length=3
   from Rug1
   group by date
   order by date
;quit;

****(3.2) REPORT MEAN AND STD FOR MEAN NUMBERS OF TOPIC AND REPLY****;
proc report data=reportdata nowd headline split='|';
   column ug topicnum,(mean std) avg_rep_num,(mean std);
   define ug/group 'USER|GROUP' width=6;
   define topicnum/'AVERAGE|TOPIC NUMBER';
   define avg_rep_num/'AVERAGE|REPLY NUMBER';
   define mean/format=5.1 'MEAN';
   define std/format=5.2 'STD';
run;

***********(4) VISUALIZE THE COMPARISON RESULT IN TIME SERIES*******;
****(4.1) PREPARE DATASET FOR PLOTTING******;
proc sql;
   create table combine as
   select a.date, a.*, b.*
   from rug2 as a left join sasug2  as b
   on a.date = b.date
;quit;

****(4.2) BUILD PLOTTING TEMPLATE*******; 
proc template;
  define statgraph compplot;
     begingraph / designwidth=1000px designheight=800px;
       entrytitle "Comparison of user groups between SAS and R";
          layout lattice / columns=1 columndatarange=union rowweights=(0.3 0.3 0.4);
            columnaxes;
               columnaxis / offsetmin=0.02 griddisplay=on;
            endcolumnaxes;
         /*(4.2.1)DEFINITION OF TOP PANEL*/
            layout overlay / cycleattrs=true yaxisopts=(griddisplay=on label=" " 
                             display=(line) displaysecondary=all);                                       
               seriesplot x=date y=R_topicnum / lineattrs=(thickness=2px) 
                   legendlabel="R's topics" name="d25"; 
               seriesplot x=date y=SAS_topicnum / lineattrs=(pattern=solid thickness=2px) 
                     legendlabel="SAS's topics" name="d50"; 
               discretelegend "d25" "d50" / across=1 border=on valign=top halign=left 
                     location=inside opaque=true;                            
            endlayout;
         /*(4.2.2)DEFINITION OF MIDDLE PANEL*/
            layout overlay / cycleattrs=true yaxisopts=(griddisplay=on label=" " 
                        display=(line) displaysecondary=all);                                       
                    seriesplot x=date y=R_zero_reply_ratio/ lineattrs=(thickness=2px) 
                     legendlabel="R's zero reply ratio" name="d2"; 
               seriesplot x=date y=SAS_zero_reply_ratio/ lineattrs=(pattern=solid thickness=2px) 
                     legendlabel="SAS's zero reply ratio" name="d5"; 
               discretelegend "d2" "d5" / across=1 border=on valign=top halign=left 
                     location=inside opaque=true;                
            endlayout;
         /*(4.2.3)DEFINITION OF BOTTOM PANEL*/
           layout overlay / yaxisopts=(griddisplay=on display=(line) label=" " displaysecondary=all)  
                        cycleattrs=true xaxisopts=(griddisplay=on);                      
               seriesplot x=date y=R_mean_reply / lineattrs=(thickness=2px) 
                    legendlabel="R's mean reply" name="ser1"; 
               seriesplot x=date y=sas_mean_reply / lineattrs=(pattern=solid thickness=2px) 
                    legendlabel="SAS's mean reply" name="ser2"; 
               seriesplot x=date y=R_max_reply / lineattrs=(thickness=2px) 
                    legendlabel="R's max reply" name="ser3"; 
               seriesplot x=date y=sas_max_reply / lineattrs=(pattern=solid thickness=2px) 
                    legendlabel="SAS's max reply" name="ser4"; 
               discretelegend "ser1" "ser2" "ser3" "ser4"/ across=1 border=on valign=top halign=left 
                    location=inside opaque=true;
           endlayout;
        endlayout;
     endgraph;
  end;
run;

****(4.3) RENDER IMAGES BY USING TEMPLATE******; 
proc sgrender data=combine template=compplot; 
run;

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

Music social network on DNA microarray

The incoming 2011 KDD Cup data mining competition [1] by Yahoo! Lab posts an interesting challenge to predict the users' ratings for individual songs out of this company’s huge music database. Unlike previous KDD Cups projects filled by tons of variables that make dimension reduction a serious concern, Yahoo! Lab provides few variables: artist/genre/album. No demographic or geographic information is disclosed. It is interesting to forecast the behavior of a web user by limited web records. Digging valuable clues out for potential following direct marketing is also rewarding. Especially while the competition datasets contain up to 1 million users, 600 thousand songs, the project is a real world web-scale analytics question.

To predict song’s rating may need two-stage modeling. The relationship looks straight-forward: rating = (genre + album + artist)* user. Building 1millinon multiply by 600 thousand models to suit each of the 1 million users on each of the 600K songs would be a formidable job. The first stage has to group users and songs to reasonable levels by unsupervised modeling methods, such as clustering. At the second stage, supervised models are to be trained for the songs and users separately. The relationship among rate, genre and artist, in 3D cubes, would decide to which group those songs or users are most likely to fit in. After scoring, in the test dataset, each song would have a song group ID and each user has a user group ID. The association between song group and user group predicts the rating values.

Song listeners get connected by their preference, purposely or not purposely. In a broad picture, those relationships formed an artist-centered social network and listeners scattered around those axis with various distances. As a result, the music fans interact each other and dwell in many social neighborhoods. In this project, the scoring card to predict song's rating will be like a huge matrix between song groups and user groups with rating as values. Biologists tackling the interaction among thousands of genes by DNA microarray would be familiar with this scenario. Usually DNA microarray is used to discover high-throughput information of relationship of multiple genes. Similarly, for music social network, this DNA microarray like matrix would show the affinity between users and songs. To increase the number of cells in this matrix, or overall number of song groups or user groups, is a good idea to bring about higher precision. However, prediction accuracy and computer resource consumption of the models are the first and foremost considerations for this particular question.

Statistics persons call data mining as statistical learning, while computer persons refer it as machine learning. Eventually both roads lead to the same way. The datasets for KDD cup 2011 data mining competition will be released on March 15. Though its 263m rating data for track1 and 63m rating data for track looks quite a lot, statistics guys probably can gear up with their high-level weapons, such as SAS and R, to explore this usually computer-geek-dominated field.

Reference: 1. 'KDD CUP 2011 from Yahoo! Lab'. http://kddcup.yahoo.com/

********(0) DOWNLOAD PREVIEW DATASET TRACK2*****;
http://kddcup.yahoo.com/dataset.track2.sample.tar.bz2

******(1) INTEGRATE PREVIEW SONG DATA*********;
data trackdata;
 infile 'C:\trackData2.txt' dlm='|' missover dsd;
 informat TrackId AlbumId ArtistId GenreId1-GenreId15 $6.;
 input TrackId AlbumId ArtistId GenreId1-GenreId15 ;
 if AlbumId = 'None' then call missing(AlbumId);
 if ArtistId = 'None' then call missing(ArtistId);
run;

proc sql; 
 select count(unique(ArtistId))
 from trackdata
;quit;

********(2) NORMALIZE SONG DATA*********;
data three;
 set trackdata (where=(missing(ArtistId)=0));
    array genre[15] GenreId1-GenreId15;
 do i=1 to 15;
  if missing(genre[i])= 0 then do;
   songgenre=genre[i]; output;
  end;
 end;
 keep ArtistId songgenre;
run;

*******(3) FAST CLUSTERING THE SONGS ONLY BY GENRES******;
proc sort data=three out=four;
 by ArtistId;
run;

proc sql;
 create table five as 
 select ArtistId, songgenre, count(songgenre) as freq
 from four
 group by ArtistId, songgenre
 order by ArtistId, songgenre
;quit;

proc transpose data=five out=five_t(drop=_name_) prefix=genre;
 by ArtistId;
 id songgenre;
 var freq;
run;

data six;
 set five_t;
 array x[*] _numeric_;
 do i=1 to dim(x);
  if missing(x[i]) = 1 then x[i]=0;
 end;
 drop i;
run;

proc fastclus data=six maxc=10 maxiter=100 out=clus;
  var genre:;
run;

*****(4) Simulate data for Genomic Heat Map********;
data music (keep=songgroup usergroup rate);
      do i=1 to 20; 
              songgroup = cats('song', i);
               do j=1 to 20;
                usergroup=cats('user', j);
                    rate= ranuni(0)*100;
                  output;
            end;
      end;
run;

proc template;
      define statgraph myheat.Grid;
      begingraph;
      layout overlay / border=true xaxisopts=(label='Song groups' ) 
        yaxisopts=(label='User groups' ) pad=(top=5px bottom=0px right=15px);
        scatterplot x=songgroup y=usergroup / markercolorgradient=rate 
       markerattrs=(symbol=squarefilled size=30) colormodel=threecolorramp name='s2';
        continuouslegend 's2' / orient=vertical location=outside valign=center
         halign=right valuecounthint=10;
      endlayout;
    endgraph;
    end;
run;

proc template;
  define Style HeatMapStyle; 
    parent = styles.harvest;
    style GraphFonts from GraphFonts "Fonts used in graph styles" /
      'GraphFootnoteFont' = (", ",8pt)
      'GraphLabelFont' = (", ",7pt) 
      'GraphValueFont' = (", ",10pt)
      'GraphDataFont' = (", ",5pt);              
    style GraphColors from graphcolors /
      "gcdata1" = cxaf1515
      "gcdata2" = cxeabb14
      "gcdata3" = cxffffff
      "gramp3cend" = cxaa081b                                              
      "gramp3cneutral" = cx000000                                          
      "gramp3cstart" = cx1ba808;
  end;
run;

ods html style=HeatMapStyle image_dpi=300 file='heatmap.html' path='d:\myfun';

ods graphics on / reset imagename='GTLHandout_Heatmaps' imagefmt=gif;
proc sgrender data=music template=myheat.grid;
run;

ods html close;