Wednesday, 23 February 2011

FCMP: BKf_BH Controlling the FDR


 Only once did I see the process for incorporation of cutting edge statistical theory into a new sas procedure and actually know the people. A long time ago, when fish just started to explore the idea of dry land living, I studied under Yoav Benjamini who together with Yosi Hochberg suggested (and named) the FDR as an alternative to the FEW. As is with many revolutionary concepts is was not accepted immediately. But once a few articles were written Yoav was contacted by sas who asked a few questions. It was some time ago so I am a little hazy on the details; but if my memory serves me right, I saw the actual letter. Subsequently Wolfinger (I Think) attended the MCP conference in Tel-Aviv and met Yoav, Yosi and Daniel Yekutieli.

 Sometime later, I was working at sas/Austria (or was it already sas/Denmark?) when a new proc was introduced MultTest. I was excited as it was relevant to my research and also included measures suggested by people I knew – Yosi Hochberg and Yoav Benjamini. I later discovered that Yosi’s measure may also be used in the mean statement for proc Anova and GLM (GT2).

 One thing jarred with me at the time. The option in proc MultTest to use the step-down FDR controlling procedure was called FDR. To my opinion it should have been BH for the Benjamini-Hochberg procedure. Apart for giving them their due and being consisted with the naming of the options such as Tukey, Dunnet etc., I knew of at least one other procedure at the time that controlled the FDR and expected more. Moreover, the naming confused between the procedure and the measure.

 Nowadays the FDR measure is mainstream, especially after its relevance to BioInformatics was recognized. More powerful procedures to control the FDR were proposed and some are implemented in sas 9.2. But I still like the elegance of the BH procedure.

 It was just natural that chose to explore the FCMP functionality through the prism of the FDR (download).

To try it run:

%Let n=3;
data test;
 * Array parameters to subroutine calls must be temporary arrays;
 array a(&n.) _temporary_;
 array b(&n.);
 array c(&n.) _temporary_;
 array d(&n.);
 input b1 b2 b3;
 do i=1 to &n.; a[i]=b[i]; end;
 call BKs_BH(a,c);
 do i=1 to &n.; d[i]=c[i]; end;
 datalines;
 0.05 0.01 0.95
 5.00 0.10 0.01
 0.05 0.05 0.05
 0.03 0.02 0.01
 run;
proc print;run;

Output
Obs     b1      b2      b3       d1      d2      d3

1     0.05    0.01    0.95    0.075    0.03    0.95
2     5.00    0.10    0.01    1.000    0.15    0.03
3     0.05    0.05    0.05    0.050    0.05    0.05
4     0.03    0.02    0.01    0.030    0.03    0.03


using sas V9.2/Base

Monday, 17 January 2011

FCMP: BKf_GateKeeper Function


By Eli Kling using sas V9.2/Base

From time to time the developers in Cary make an addition that excites even sas–dinosaurs like myself. Not everyone sees the potential of, or even the need for proc FCMP. However, I think it signals a major change in the way sas is used and developed and opens the door for cutting edge statistical modelling also on platforms that have ‘only’ Base/Stat/Graph.

The obvious challenge I first tackled as a learning exercise was calculating the Benjamini-Hochberg FDR correction to a vector of p-values. It was the focus of my dissertation a thousand years ago. I will discuss it in a following entry.

Good practice calls to first check that the vector passed in is appropriate and abort the processing otherwise. The BKf_GateKeepr function below checks that all the elements of a vector contain values between 0 and 1 and returns 0 for OK and 1 for Failed. If the check failed it also puts a message in the log.

As I intend to expand the functionality of this function I included a ‘Check’ switch indicating what check to perform. In the version below only ‘PVALUE’ is programmed.

You will also notice that the sas editor has a problem to correctly colour the key words. I am sure that will be fixed in later versions.

proc delete data=sasuser.BusinessKenLibrary; run;
options cmplib=sasuser.BusinessKenLibrary;

proc fcmp outlib=sasuser.BusinessKenLibrary.dev;

/* RC=BKf_GateKeeper(Vector,Check)
 a utility to be used to check vector inputs
 Input Vector - the vector of values to check
       Check  - what to check
       pvalue: verify all values between 0 and 1
 output RC - 0 = passed check
             1 = failed check
*/
 function BKf_GateKeeper(vector(*),Check$);
 put "check" Check;
 N=dim(Vector);
 RC=0;
  if Upcase(Check)="PVALUE" then do;
    do i=1 to N;
       if not(0<=vector[i]<=1) then RC=1;
      end;
    if RC=1 then put "<< BK Values Gate Keeper: A vector of pvalues was expected but at least one value was not in the interval [0,1]. Note, it might have been a missing value. >>";
   end;
 return(RC);
 endsub;

run;


/********/
/* Test */
/********/
data test;
 * Array parameters to subroutine calls must be temporary arrays;
 array a(3) _temporary_;
 a[1]=5;
 a[2]=0.01;
 a[3]=0.96;
 RC=BKf_GateKeeper(a,"Pvalue");
 run;
proc print;run;
/*****************************************************************/
/* Code End                                                      */
/*****************************************************************/