/**************************************************************************/
/* Read in HIV-1 SAS dataset, [y:X] */ 

options nonotes;
%let row = 317;    /* number of patients */
%let col = 283;    /* first column corresponds to RC outcome, remaining columns to PR and RT codons */
%let boots = 7500; /* number of bootsrap samples */
%let k = 5;        /* allowed number of false positives for gFWER AMTP */
%let q = 0.05;     /* allowed proportion of false positives for TPPFP AMTP */
%let nt= 282;      /* number of tested hypotheses */

libname resample "c:\hivexample";

/**************************************************************************/
/* lmt: t-statistics for univariate linear regression model */

%macro lmt;
  %do a = 2 %to &col;
    proc reg data=resample.hivdata noprint outest=outest&a tableout; 
        model VAR1=VAR&a;     
    data outest&a(rename=(VAR&a=t)); set outest&a; where _type_='T'; keep VAR&a;
    proc append base=tstats data=outest&a;
  %end;
%mend;
%lmt;
proc print data=tstats; title "t-statistics for data";
run;

/**************************************************************************/
/* boot and bootnull: Non-parametric bootstrap null shift-transformed test statistics null distribution */

%macro boot;
%do j=1 %to &boots;
proc surveyselect noprint data=resample.hivdata out=datanew
        seed=&j
        method=urs         
        rep=&row
        sampsize=1 
        stats;
%do a = 2 %to &col;
proc reg data=datanew noprint outest=outest&a tableout; 
        model VAR1=VAR&a;     
data outest&a(rename=(VAR&a=t)); set outest&a; where _type_='T'; keep VAR&a;  
proc append base=tstatsB data=outest&a;
%end;
%end;
%mend;
%boot;
quit;

/***********************/
%macro bootnull;
proc iml;
        use tstatsB;
        read all var {t} into x;
        close tstatsB;
        nt=&nt; 
        bt=&boots;
        tB=shape(x,bt,nt);
        b=tB[:,]; 
        bb=J(bt,nt,0);
        do i=1 to bt;
                bb[i,]=(tB[i,] - b);
        end;
        bb2=shape(bb,bt*nt,1); 
        create Qo from bb2; 
        append from bb2;  
        close Qo;
        quit;
%mend;
%bootnull;

/**************************************************************************/
/* ssmaxT: FWER-controlling single-step maxT procedure */

%macro ssmaxT;
proc iml;
        use Qo;
        read all into bb;
        close Qo; 
        nt=&nt; 
        bt=&boots;
        Qo=shape(bb,bt,nt);
        mx=J(1,bt,0);
        do i=1 to bt;
                mb=abs(Qo[i,]);
                mx[,i]=max(mb); 
        end;
        c=mx; 
        brank=rank(mx); 
        mx[brank]=c;
        use tstats;
        read all var {t} into t;
        close tstats;
        print t;
        pval=J(1,nt,0);
        do j=1 to nt;
                tmp=J(1,bt,0);
                do i=1 to bt;
                        if abs(t[j]) < mx[i] then tmp[i]=1;
                end;
                st=sum(tmp);
                pval[j]=st/bt;
        end;
        create fwer from pval; 
        append from pval;  
        close fwer;
        quit;
proc print data=fwer;
run;
%mend;
%ssmaxT;

/**************************************************************************/
/* gfwer: gFWER-controlling augmentation multing testing procedure */

%macro gfwer;
proc iml;
        use fwer;
        read all into pval;
        close fwer;
        k=&k; 
        nt=&nt;
        gp=J(1,nt,0);
        j=k+1;
        c=pval;
        brank=rank(pval); 
        pval[brank]=c;
        do i=1 to (nt-k);
                gp[j]=pval[i];
                j=j+1;
        end;
        gg=gp[brank];
        gpval=shape(gg,1,nt);
        create gfwer from gpval; 
        append from gpval;  
        close gfwer;
        quit;
proc print data=gfwer;
run;
%mend;
%gfwer;

/**************************************************************************/
/* tppfp: TPPFP-controlling augmentation multing testing procedure */

%macro tppfp;
proc iml;
        use fwer;
        read all into pval;
        close fwer;
        q=&q; 
        nt=&nt; 
        c=pval;
        brank=rank(pval); 
        pval[brank]=c;
        tp=J(1,nt,0);
        do i=1 to nt;
                m=ceil(i*(1-q));
                tp[i]=pval[m];
        end;
        tt=tp[brank];
        tpval=shape(tt,1,nt);
        create tppfp from tpval; 
        append from tpval;  
        close tppfp;
        quit;
proc print data=tppfp;
run;
%mend;
%tppfp;

/*************************************************************************/