%macro Gatekeeper(dataset,outdata);
/*
Computes adjusted p-values for the Bonferroni-based parallel gatekeeping procedure.

Written by Alex Dmitrienko. 

Inputs:

DATASET = Data set defining the structure of the problem, including the family definition, 
raw p-values and weights.

OUTDATA = Output data set with adjusted p-values.
*/
proc iml;    
    use &dataset;
    read all var {family weight raw_p} into data;            
    data=t(data);
    family=data[1,];
    weight=data[2,];
    raw_p=data[3,];
    nhyps=ncol(data); nfams=family[ncol(data)]; nints=2**nhyps-1;
    h=j(nints,nhyps,0);
    do i=1 to nhyps;
        do j=0 to nints-1;
            k=floor(j/2**(nhyps-i));
            if k/2=floor(k/2) then h[j+1,i]=1;
        end;
    end;    
    v=j(nints,nhyps,0); modv=j(nints,nhyps,0); 
    hyp=j(nints,nhyps,0); adjp=j(nhyps,1,0);
    start bonf(p,w);
        bonfp=min(p[loc(w)]/w[loc(w)]);
        return(bonfp);
    finish bonf; 
    do i=1 to nints;        
        r=1; tempv=j(1,nhyps,0); tempmodv=j(1,nhyps,0);        
        do j=1 to nfams;            
            window=(family=j);            
            sumw=sum(window#weight#h[i,]);                     
            if (j<nfams) & sumw>0 then do; * Parallel testing;
                tempv=r#window#weight#h[i,];                
                if sum(h[i,]#(family>j))=0 
                then tempmodv=r#window#weight#h[i,]/sumw; 
                else tempmodv=tempv;
            end; 
            if (j=nfams) & sumw>0 then do; * Last family;
                tempv=r#window#weight#h[i,]/sumw; tempmodv=tempv;
            end;                                                     
            if sumw>0 then do; 
                r=r-sum(tempv); v[i,]=v[i,]+tempv; modv[i,]=modv[i,]+tempmodv;            
            end;            
        end;        
        hyp[i,]=h[i,]*bonf(raw_p,v[i,]);                    
    end;
    do i=1 to nhyps; adjp[i]=max(hyp[,i]); end;    
    create adjp from adjp[colname={adjp}];
    append from adjp;    
    quit;
data &outdata;
    merge &dataset adjp;
    format adjp 6.4;
    run;
%mend Gatekeeper;
* Example (study with three families of hypotheses);
data study;
    input hyp $ family weight raw_p;
    datalines;
    H11 1 0.5 0.010
    H12 1 0.5 0.005
    H21 2 0.5 0.030
    H22 2 0.5 0.015
    H31 3 0.5 0.035
    H32 3 0.5 0.010
    ;
%Gatekeeper(dataset=study,outdata=out);
proc print data=out;
    title "Bonferroni-based parallel gatekeeping procedure";
    run;