%macro Gatekeeper(dataset,test,outdata);
/*
Inputs:

DATASET = Data set to be analyzed.

TEST    = B, MB, or S for Bonferroni, modified Bonferroni
          or Simes gatekeeping procedure, respectively.

OUTDATA = Output data set with adjusted p-values.
*/
proc iml;    
    use &dataset;
    read all var {family testtype weight relimp raw_p} into data;        
    data=t(data);
    nhyps=ncol(data); nfams=data[1,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; 
    start simes(p,w);
        comb=t(p[loc(w)])//t(w[loc(w)]); temp=comb;   
        comb[,rank(p[loc(w)])]=temp;    
        simesp=min(comb[1,]/cusum(comb[2,]));
        return(simesp);
    finish simes;  
    do i=1 to nints;        
        r=1; tempv=j(1,nhyps,0); tempmodv=j(1,nhyps,0);        
        do j=1 to nfams;            
            window=(data[1,]=j);            
            sumw=sum(window#data[3,]#h[i,]);         
            testtype=sum(window#data[2,]);             
            if (testtype=0 & j<nfams) & sumw>0 then do; * Parallel testing;
                tempv=r#window#data[3,]#h[i,]#((1-data[4,])+data[4,]/sumw);                
                if sum(h[i,]#(data[1,]>j))=0 
                then tempmodv=r#window#data[3,]#h[i,]/sumw; 
                else tempmodv=tempv;
            end; 
            if (testtype>0 | j=nfams) & sumw>0 then do; * Serial testing or last family;
                tempv=r#window#data[3,]#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;        
        if &test='B' then hyp[i,]=h[i,]*bonf(data[5,],v[i,]);
        if &test='MB' then hyp[i,]=h[i,]*bonf(data[5,],modv[i,]);
        if &test='S' then hyp[i,]=h[i,]*simes(data[5,],modv[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;
    run;
%mend Gatekeeper;
data raexample;
    input hyp $ family testtype weight relimp raw_p;
    datalines;
    H11 1 0 0.5 0 0.010
    H12 1 0 0.5 0 0.005
    H21 2 0 0.5 0 0.030
    H22 2 0 0.5 0 0.015 
    H31 3 0 0.5 0 0.035
    H32 3 0 0.5 0 0.010 
    ;
%Gatekeeper(dataset=raexample,test='B',outdata=out);
proc print data=out;
    run;