/*
%Component macro

Performs inferences for individual endpoints after a significant global test
using the O'Brien OLS test, GLS test, modified GLS test and rank-sum test.

For more information about this macro, see Chapter 2 in

Dmitrienko, A., Molenberghs, G., Chuang-Stein, C., Offen, W. (2005). Analysis of 
Clinical Trials Using SAS: A Practical Guide. SAS Press: Cary, NC.

Written by Alex Dmitrienko. 

Inputs:
 
DATASET = Data set to be analyzed
 
GROUP   = Name of the group variable in the data set
 
NGROUPS = Number of groups in the data set
 
VARLIST = List of variable names corresponding to
          multiple endpoints
 
TEST    = OLS, GLS, MGLS or RS, for OLS test, GLS test, 
          modified GLS test or rank-sum test, respectively

Example:

%Component(dataset=oa,group=group,ngroups=2,varlist=pain physical,test="OLS");

*/
* Global test;
%macro globtest(dataset,group,ngroups,varlist,test);
    ods listing close;
    %if &test='RS' %then %do;
    proc rank data=&dataset out=ranks;
        var &varlist;    
    data comp;
        set ranks;
        array endp{*} &varlist;
        comp=0;
        do i=1 to dim(endp);
            comp=comp+endp{i};
        end;
    proc mixed data=comp;
        class &group;
        model comp=&group;
        ods output tests3=pval;
    data pval;
        set pval;
        format fvalue 5.2 ndf 3.0 ddf 3.0 adjp 6.4;
        ndf=numdf;
        ddf=dendf;
        adjp=probf;
        label fvalue='F-value' adjp='Global p-value';
        keep ndf ddf fvalue adjp;
    %end;
    %else %do;
    %let m=1;
    %let word=%scan(&varlist,&m);
    %do %while (&word ne);
        %let m=%eval(&m+1);
        %let word=%scan(&varlist,&m);
    %end;
    %let m=%eval(&m-1);
    data stand;
    %do i=1 %to &m;
        %let var=%scan(&varlist,&i);
        proc glm data=&dataset;
            class &group;
            model &var=&group;
            ods output FitStatistics=est(keep=rootmse depmean);
        data stand;
            set stand est;
    %end;
    data stand;
        set stand;
        if rootmse^=.;
    data _null_;
        set &dataset nobs=m;
        call symput('n',trim(put(m,5.0)));
    proc corr data=&dataset outp=corr(where=(_type_="CORR")) noprint;
        var &varlist;
    proc iml;
        use &dataset var{&varlist};
        read all into data;
        use stand var{depmean};
        read all into mean;
        meanmat=j(nrow(data),ncol(data),1)*diag(mean);
        use stand var{rootmse};
        read all into pooledsd;
        sdmat=inv(diag(pooledsd));
        use corr var{&varlist};
        read all into r;
        stand=(data-meanmat)*sdmat;
        rinv=inv(r);
        if &test='OLS' then comp=stand*j(&m,1);
        if &test='GLS' then comp=stand*rinv*j(&m,1);
        if &test='MGLS' then comp=stand*sqrt(diag(rinv))*j(&m,1);
        create comp from comp;
        append from comp;
    data comp;
        merge comp &dataset;
        comp=col1;
        drop col1;
    proc mixed data=comp;
        class &group;
        model comp=&group;
        ods output tests3=pval;
    data pval;
        set pval;
        format fvalue 5.2 ndf 3.0 ddf 3.0 adjp 6.4;
        ndf=&ngroups-1;
        ddf=&n-&m*&ngroups;
        adjp=1-probf(fvalue,ndf,ddf);
        label fvalue='F-value' adjp='Global p-value';
        keep ndf ddf fvalue adjp;
    %end;
    ods listing;
    proc print data=pval noobs label;
        run;
%mend globtest;

%macro component(dataset,group,ngroups,varlist,test);
    data hyp;
        length list $50;
        array h{*} h1-h4;
        do i=15 to 1 by -1;
            string=put(i,binary4.);
            hypnum=16-i;
            do j=1 to 4;
                h{j}=substr(string,j,1);
            end;
            list=" ";
            do j=1 to 4;
                if h{j}=1 then list=
                trim(list) || " " || scan("&varlist",j);
            end;
        output;
        end;
        keep hypnum h1-h4 list;    
    data pvalcomb;
    * Compute a global test p-value for each intersection hypothesis;
    %macro pval;
        %do j=1 %to 15;
        data _null_;
            set hyp;
            if hypnum=&j then call symput("list",trim(list));
            run;				
            %globtest(&dataset,&group,&ngroups,&list,&test);
        data pvalcomb;
            set pvalcomb pval;
            keep adjp;
            %end;
    %mend pval;
    %pval;
    data decrule;
        merge hyp pvalcomb(where=(adjp^=.));
        array h{*} h1-h4;
        array hyp{*} hyp1-hyp4;
        do i=1 to 4;
            hyp{i}=adjp*h{i};
        end;
        keep hyp1-hyp4;
    proc means data=decrule noprint;
        var hyp1-hyp4;
        output out=indadjp(keep=adjp1-adjp4)
        max(hyp1-hyp4)=adjp1-adjp4;
    data indadjp;
        set indadjp;
        format adjp1-adjp4 6.4;
    proc print data=indadjp noobs;
        title "Adjusted p-values for individual hypotheses";
        run;
%mend component;
* Example: Analysis of components in a rheumatoid arthritis trial with four endpoints;
data ra;
    input group $ sjc tjc pta pha @@;
    datalines;
    Placebo   -5  -9  -14  -21 Placebo   -3  -7   -5  -25
    Placebo   -7  -4  -28  -15 Placebo   -3   0  -17   -6
    Placebo   -4  -1   -5    5 Placebo    0   5   -8  -11
    Placebo   -3   1   15    0 Placebo    2   6   15   27
    Placebo   -1  -4  -11   -8 Placebo    0   1    8   12
    Placebo    2  -2    6   -9 Placebo    8   2   11   33
    Therapy   -7  -1  -21   -9 Therapy   -6 -11  -36  -12
    Therapy   -3  -7  -14  -21 Therapy   -4   2   10  -10
    Therapy  -11  -4  -28  -45 Therapy   -4  -1  -11  -23
    Therapy   -3  -1   -7  -15 Therapy   -5  -9  -36  -15
    Therapy   -4  -9  -35  -32 Therapy  -11 -10  -47  -31
    Therapy    3  -1    6   17 Therapy   -1  -9   -5  -27
    ;
%component(dataset=ra,group=group,ngroups=2,varlist=sjc tjc pta pha,test="OLS");