/**********************************************************************************************
 Program Name: wandel.sas
 Program Date: 28 February 2006
 Programmer: Peter Kraft

%wandel performs permutation adjustments for multiple testing, where each test is the 
likelihood ratio test from  a univariate, unconditional logistic regression. 

%wandel compares the sum of the k smallest observed p-values to the permutation distribution
of the k-smallest p vales. (All observations are assumed exchangeable under the null.) For k=1 
this is the standard permutation adjustment for the smallest observed p-value. For k>1 this is 
the Rank Tuncated Product Test [Dudbridge and Koelman Genet Epidemiol 25(4):360-6].

input   in:       input data set
        vars:     list of variables to be tested (e.g. counts of SNP minor alleles)
        nrep:     number of permutations
        k:        consider k smallest p-values
        caco:     Outcome variable name

output  OutTable data set with observed product of smallest k p-values and its permutation
        p-value (the proportino of permuted products that are smaller than or equal to the 
        observed product).
**********************************************************************************************/

/*****************************************************************/
/*                                 DEFINE MACROFUNCTION NUMARGS */
%macro numargs(arg);
   %if &arg= %then %do;
        0
      %end;
    %else %do;
       %let n=1;
       %do %until (%scan(&arg,%eval(&n),%str( ))=%str()); 
         %let n=%eval(&n+1);
       %end;
       %eval(&n-1)
    %end;
%mend numargs;

/*****************************************************************/
/*                                DEFINE MACROFUNCTION WANDEL    */
%macro wandel(in=,             /* Input data set                 */
              vars=,           /* List of variables to be tested */
              nrep=,           /* Number of permutations         */
              k=1,             /* Consider k smallest p-values   */
              caco=caco);      /* Outcome variable name          */
/*****************************************************************/
options nonotes;
ods select none;
%put;
%put ******************************************************************;
%put wandel macro;
%put Calculates rank truncated product test for multiple marginal tests;
%put Created and maintained by Peter Kraft (pkraft@hsph.harvard.edu);
%put Last modified 28 February 2006;
%put ******************************************************************;
%put;
/*                     Calcluate observed rank truncated product */
data ptemp; run;
%do i=1 %to %numargs(&vars);
  ods output GlobalTests=ptemptemp;
  proc logistic data=∈
     model &caco=%scan(&vars,&i);
  run;
  data ptemp; set ptemp ptemptemp(where=(test='Likelihood Ratio')); keep ProbChiSq; run;
%end;
data pobs; set ptemp; where ProbChiSq ne .; run;
proc sort data=pobs; by ProbChiSq; run;
data pobs; set pobs(obs=&k);
  retain LogProbChiSq;
  if _N_ eq 1 then LogProbChiSq=log(ProbChiSq);
     else LogProbChiSq=LogProbChiSq+log(ProbChiSq);
  if _N_ eq &k then output;
  keep LogProbChiSq;
run;
/*  Calcluate permutation distribution of rank truncated product */
data nophen; set ∈ drop &caco; run;
data phen; set ∈ keep &caco; run;
data pperm; run;
%do j=1 %to &nrep;
  data phen; set phen; dummysort=rannor(0); run;
  proc sort data=phen; by dummysort; run;
  data perm∈ merge nophen phen; run;
  data ptemp; run;
  %do i=1 %to %numargs(&vars);
    ods select none;
    ods output GlobalTests=ptemptemp;
    proc logistic data=perm∈
       model &caco=%scan(&vars,&i);
    run;
    data ptemp; set ptemp ptemptemp(where=(test='Likelihood Ratio')); keep ProbChiSq; run;
  %end;
  data ppermtemp; set ptemp; where ProbChiSq ne .; run;
  proc sort data=ppermtemp; by ProbChiSq; run;
  data ppermtemp; set ppermtemp(obs=&k);
    retain LogProbChiSq;
    if _N_ eq 1 then LogProbChiSq=log(ProbChiSq);
       else LogProbChiSq=LogProbChiSq+log(ProbChiSq);
    if _N_ eq &k then output;
    keep LogProbChiSq;
  run;
  data pperm; set pperm ppermtemp; run;
%end;
data pperm; set pperm(where=(LogProbChisq ne .)); run;
proc sort data=pperm; by LogProbChiSq; run;
data OutTable;
  set pobs pperm;
  retain target count (0 0);
  if _N_ eq 1 then target=LogProbChiSq;
  if _N_ gt 1 and LogProbChiSq le target then count=count+1;
  if _N_ gt 1 and LogProbChiSq gt target then do;
     ObsProdP=exp(target);
     PermPValue=count/&nrep;
	 output;
	 stop;
  end;
  keep ObsProdP PermPValue;
run;
proc datasets; delete nophen permuseme phen nophen ppermtemp ptemptemp pperm ptemp pobs; quit;
ods select all;
proc print data=OutTable noobs;
  title "Observed Rank Truncated Product and Permutation P-value (nreps=&nrep)";
run;
options notes;
%mend wandel;