/*   RANOMV IS A TEST FOR THE EQUALITY
     OF K VARIANCES BASED ON K INDEPENDENT
     SAMPLES OF SIZE N   */
/*******************************
        INPUT DATA
********************************/
data data1 ;
        pops = 7;
        samp = 10;
        ns =1000;
        do tool=1 to pops;
                do j=1 to samp;
                diam = tool*rannor(-1) + 100;
                do shuf = 1 to ns;
                        shufno = uniform(-1);
                        output ;
        end; end; end;
        drop j ;



/*******************************
       DEFINE MACRO
********************************/

%macro ranomv(
   k=,     /* the number of populations being compared */
   n=,     /* the sample size                          */
   alpha=, /* level of significance                    */
   ds=,    /* the data set containing the observations */
   varname=,         /* the variable name for the observations   */
   classvar=,    /* the variable name for the populations    */
   ns =,   /* the number of shufles */
   tops=);       /*  1 = randanomv-R     */




/***************************
Print Data Set
***************************/



data basedat;
        set data1;
        if  shuf > 1 then delete;
        title 'basedat: original data set';
        proc print;  var &classvar &varname;
run;

/********************************
   DETERMINE RANOMV CRITICAL VALUES
********************************/
title 'initial data set';
       /* proc print data = &ds;*/
run;



data shufdat1;
      set &ds;
      proc sort; by shuf shufno;
      title 'shuffled data';
      /* proc print;*/
run;





data shufdat2;
        set shufdat1;
        drop &classvar;
        title 'shuffled data without tool';
        /*proc print;*/
run;




data origdat0;
        set data1;
        proc sort; by shuf;
        title 'original data with diam';
        /*proc print;*/



data origdat;
        set data1; drop &varname;
        proc sort; by shuf;
        title 'original data without diam';
        /* proc print; */



data shufdat3;
       merge shufdat2 origdat;
       title 'merged data';
       /* proc print;*/
run;


proc means data = shufdat3 noprint;
        by shuf &classvar; var &varname ; output out = stats1 var = varx;
        title 'variance of each shuffled sample';
        /*proc print data = stats1;*/
run;



proc means data=stats1 noprint;
        by shuf; var varx; output out = stats2 sum = sumvarx;
        /*proc print data = stats2;*/
run;



data vardat;
        merge stats1 stats2; by shuf;
        varrat = varx/sumvarx;
        title 'vardat: variances of shuffled samples along with sums and ratios';
        /* proc print data = vardat; */
run;




data admaxwk;
        set vardat;
        proc means noprint;
        by shuf; var varrat; output out = admax max = maxrat;
        /* proc sort data = admax ; by maxrat;*/
        title 'admax: admax distribution';
        /* proc print data = admax;*/
run;



data adminwk;
        set vardat;
        proc means noprint;
        by shuf; var varrat; output out = admin min = minrat;
        title 'admin: admin distribution';
        /* proc sort data = admin; by minrat; */
        /*proc print data = admin;*/
run;




proc rank data=admin   out=rkadmin;
        var minrat;
        ranks rminrat;
        /* proc print data = rkadmin; */
run;



data critlow;
        set rkadmin;
        ranklow = &ns-floor((&ns+1)*(1-&alpha/2)-1)-1;
        if rminrat > ranklow then delete; if rminrat < ranklow then delete;
        lowcrit = minrat;
        title 'crit low';
        /*proc print data = critlow;*/
run;



proc rank data=admax   out=rkadmax;
        var maxrat ;
        ranks rmaxrat;
        title 'rkadmax:  ranked admax';
        /*proc print data = rkadmax; */
run;




data crithi;
        set rkadmax;
        rankhi = &ns-floor((&ns+1)*(&alpha/2)-1);
        if rmaxrat > rankhi then delete; if rmaxrat < rankhi then delete;
        hicrit = maxrat;
        title 'crithi: crit high';
        /*proc print;*/
run;




data critdum1;
        merge crithi critlow;
        codex =1;
        title 'critdum1: critcial values';
        /*proc print;*/
run;
data c2;set critdum1;codex=2;
data c3;set critdum1;codex=3;
data c4;set critdum1;codex=4;
data c5;set critdum1;codex=5;
data c6;set critdum1;codex=6;
data c7;set critdum1;codex=7;
data c8;set critdum1;codex=8;
data c9;set critdum1;codex=9;
data c10;set critdum1;codex=10;
data c11;set critdum1;codex=11;
data c12;set critdum1;codex=12;
data c13;set critdum1;codex=13;
data c14;set critdum1;codex=14;
data c15;set critdum1;codex=15;
data c16;set critdum1;codex=16;
data c17;set critdum1;codex=17;
data c18;set critdum1;codex=18;
data c19;set critdum1;codex=19;
data c20;set critdum1;codex=20;


data critvals;
        set critdum1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20;
        if codex > &k then delete;
        title 'critvals: critical values';
        proc print data=critvals;  var  rankhi hicrit ranklow lowcrit ;


/********************************
     DETERMINE DECISION LINES
*********************************/







data basedat2;
        set origdat0;
proc means data = basedat2 noprint;
        by shuf &classvar;
        var &varname ;
        output out = stats4a var = var0 std = std0;
        title ' variances and standard deviations of original data shuf replicated ';
        /*proc print data = stats4a; */
run;



proc means data=stats4a noprint;
        by shuf;
        var var0;
        output out = stats5a sum = sumvar0;
        /* proc print data = stats5a; */
run;

data vardat0;
        merge stats4a stats5a; by shuf;
        varrat0 = var0/sumvar0;
        title 'variances of original data shuf repl samples along with sums and ratios';
       /* proc print data = vardat0;*/
run;




data vardat0i;
        set vardat0;
        if shuf > 1 then delete;
        avgvar = sumvar0/&k;
        title 'vardat01i: variance data for initial sample';
        /* proc print data = vardat0i;*/



data vardat1;
        merge vardat vardat0; by shuf &classvar;
        title 'variances of original data + shuf data repl samples along with sums and ratios';
        /* proc print data = vardat1;*/
run;
data adminwk0;
        set vardat1;
        proc means noprint;
        by shuf; var varrat0; output out = admin0 min = minrat0;
data admaxwk0;
        set vardat1;
        proc means noprint;
        by shuf; var varrat0; output out= admax0 max = maxrat0;




data dldat1;
        merge critvals vardat0i;
        UDLVAR = sumvar0*hicrit;
        CLVAR = avgvar;
        LDLVAR = sumvar0*lowcrit;
        UDL = sqrt(UDLVAR);
        CL = sqrt(CLVAR);
        LDL = sqrt(LDLVAR);
        nameit = &classvar;
        title 'dldat1: data to determine decision lines';
        /*proc print;*/
        title "RANDANOMV decision table for equality of &k variances";
        proc print data = dldat1; var &classvar UDL CL LDL std0 ;



run;



/********************************
     OUTPUT ANOMV DECISION CHART
*********************************/
 proc gplot data=dldat1 ;

        plot std0*&classvar=4
                 ldl*&classvar=1
                 cl*&classvar=2
                 udl*&classvar=3
                 /overlay
                 haxis=axis2
                /* annotate=bars */
                legend;
        symbol1 c=BLUE,i=join, l=14, v=none;
        symbol2 c=BLUE, i=join, l=1, v=none;
        symbol3 c=BLUE, i=join, l=2 v=none;
        symbol4 c=BLACK, i=none, v=star;
        axis2 order=(1 to &k by 1) offset=(2) label=(h=1.5);
        title1 "RANDANOMV Decision Chart for &varname";
        title2 "Alpha = &alpha and &ns Permutation Shuffles";
        title3 "Standard Deviation Plotted";
run;

/*********************
Calculate p-values
**********************/

data pvaldat;
        merge admax admin admin0 admax0;
        if maxrat > maxrat0 then sighi = 1; else sighi =0;
        if minrat < minrat0 then siglo = 1; else siglo = 0;
        sig = min(sighi+siglo,1);
        title 'pvalue data';
        /*proc print data=pvaldat; */


proc means data = pvaldat noprint;
        var sig; output out = pvaldat1 mean = emppval;
        /*title1 "p-value for test for equality of &k variances";
        titile2 "&n observations per group and &ns shuffles";*/
        /*proc print data = pvaldat1;  var emppval; */
proc means data = pvaldat noprint;
        var sighi; output out = pvaldat2 mean = pvalhi;
        /*title1 "p-value for test for equality of &k variances";
        titile2 "&n observations per group and &ns shuffles";*/
        /*proc print data = pvaldat2;  var pvalhi; */
proc means data = pvaldat noprint;
        var siglo; output out = pvaldat3 mean = pvallow;
        /*title1 "p-value for test for equality of &k variances";
        titile2 "&n observations per group and &ns shuffles";*/
        /*proc print data = pvaldat3;  var pvallow;*/
data pvalout;
        merge pvaldat2 pvaldat3;
        title 'Empirical p-values should should be compared to alpha/2';
        proc print data=pvalout; var pvallow pvalhi;
run;
%mend ranomv;
%ranomv(k=7,n=10,alpha=0.10,ds=data1,varname=diam,classvar=tool,ns=1000,tops=1);