/*** Full Compositional Analysis  05.20.97 ***/

%Let sims=999; * input number of simulations;
options nodate ls=78 ps=100 nonumber formdlim=' ' stimer nosource;

filename rwilks 'rwilks.txt';
filename owilks 'owilks.txt';
filename use_avl 'radvsmcp.txt';  * raw data as u1 u2 ... uD a1 a2 ... aD;

data numb;
/*** Scan the data to see how many habs and obs there are ***/
 infile use_avl end=last missover firstobs=2;
 numb=_n_;
 do cols=1 to 50;
   input x @; habs=cols/2;
   if x ne . then output;
 end;
run;

data numb;
/*** Turning the # of habs and # of obs into macros ***/
 set numb end=last;
 if last then do;
   call symput('numb',left(put(numb,3.)));
   call symput('habs',left(put(habs,3.)));
 end;
run;

data raw;
/*** observed RAW data (with zero entries adjusted) ***/
 infile use_avl firstobs=2;
 %Let sdim=%eval(&habs-1);
 %Let ldim=%eval(&habs*&habs);
 input u1-u&habs a1-a&habs;
run;

data obspref;
/*** Constructing the observed antisymmetric matrix ***/
 array avail{&habs} a1-a&habs;
 array used{&habs} u1-u&habs;
 array compo{&habs,&habs};
 set raw;
 do i=1 to &habs;
   do j=i+1 to &habs;
     compo[i,j]=log(used[i]/used[j])-log(avail[i]/avail[j]);
     compo[j,i]=-compo[i,j];
   end;
 end;
 keep compo1-compo&ldim;
run;

proc means data=obspref mean stderr t noprint;
/*** Getting the mean of the observed log ratios (i.e. summing over the birds) ***/
 output out=obsmean mean=mean1-mean&ldim stderr=stderr1-stderr&ldim t=t1-t&ldim;
run;

data obsrank;
/*** Finding the number of positive rows in the observed matrix ***/
 array compo{&habs,&habs} mean1-mean&ldim;
 array prank{&habs} prank1-prank&habs;
 do i=1 to &habs;
   prank[i]=0;
 end;
 set obsmean;
 do i=1 to &habs;
   do j=1 to &habs;
     if compo[i,j] > 0 then do;
      prank[i]=prank[i]+1;
     end;
   end;
 end;
 drop mean1-mean&ldim stderr1-stderr&ldim i j _type_ _freq_;
run;

data obsvec;
/*** Reading off the OBSERVED compositions from the last column ***/
 array compo{&habs,&habs};
 array ds{&sdim} d1-d&sdim;
 set obspref;
 do i=1 to &habs-1;
   j=&habs;
   ds[i]=compo[i,j];
 end;
run;

proc printto new print=owilks;
/*** Routing the output to owilks.txt ***/
run;

proc glm data=obsvec;
/*** The actual REAL MANOVA ***/
 model d1-d&sdim = / ss3 nouni;
 manova H=intercept;
run;

proc printto print=print;
/*** Routing the output back to the output window ***/
run;

data obslam;
/*** Grepping the actual MANOVA output for info ***/
 infile owilks firstobs=20;
 input @ 'Overall INTERCEPT Eff' test$ @ "Wilks' Lambda"
  wilko fval df1 df2 tpval;
 lratioo=-&numb*log(wilko);
 output;
 keep wilko lratioo df1 tpval;
run;

/*** Initializing all counts to one ***/
%macro initial;
 %Global lcount;
 %Let lcount=1;
 %do j=1 %to &ldim;
    %Global tcount&j;
    %let tcount&j=1;
 %end;
%mend;
%initial;

/*** Macro loop starts here ***/
%macro loopy;
%do sim=1 %to &sims;

data ranmat;
/*** Randomizing the antisymmetric matrix ***/
 array compo{&habs,&habs};
 array random{&habs,&habs};
 merge obspref;
 retain seed1 %eval(1613217781+&sim);
   do i=1 to &habs;
     do j=i+1 to &habs;
       random[i,j]=compo[i,j]*(2*int(2*ranuni(seed1))-1);
       random[j,i]=-random[i,j];
     end;
   end;
 drop i j seed1;
run;

data ranvec;
/*** Reading off the RANDOM compositions from the last column ***/
 set ranmat;
 array random{&habs,&habs};
 array ns{&sdim} n1-n&sdim;
 do i=1 to &habs-1;
   j=&habs;
   ns[i]=random[i,j];
 end;
 keep n1-n&sdim;
run;

proc printto new print=rwilks;
/*** Routing the output to rwilks.txt ***/
run;

proc glm data=ranvec;
/*** Doing a random MANOVA ***/
 model n1-n&sdim = / ss3 nouni;
 manova H=intercept;
run;

proc printto print=print;
/*** Routing the output back to the output window ***/
run;

data ranlam;
/*** Grepping the MANOVA output for Wilks lambda ***/
 infile rwilks firstobs=20;
 input @ 'Overall INTERCEPT Eff' test$ @ "Wilks' Lambda" wilky;
 lratio=-&numb*log(wilky);
 output;
run;

proc means data=ranmat t noprint;
/*** Doing a random t-test ***/
 var random1-random&ldim;
 output out=ranmean t=rt1-rt&ldim;
run;

data monlam;
/*** Comparing the two L.R.S.s (lambdas) ***/
 merge obslam ranlam;
 if lratio >= lratioo then countl = 1;
 if lratio < lratioo then countl = 0;
 output;
 call symput('lcount', left(put(countl,3.))+&lcount);
run;

data mont;
/*** Comparing all of the t statistics ***/
 array obsv{&ldim} t1-t&ldim;
 array ts{&ldim} rt1-rt&ldim;
 array countt{&ldim} countt1-countt&ldim;
 merge obsrank ranmean;
   do i=1 to &ldim;
     if abs(ts[i]) >= abs(obsv[i]) then countt[i] = 1;
     if abs(ts[i]) < abs(obsv[i]) then countt[i] = 0;
     if obsv[i]= . then countt[i]=.;
   end;
   output;
   %do j=1 %to &ldim;
      i=&j;
      call symput('tcount'||left(&j),left(put(countt[i],3.))+&&tcount&j);
   %end;
keep t1-t&ldim prank1-prank&habs countt1-countt&ldim;
run;

/*** Macro loop finishes here ***/
%end;
%mend loopy;
/*** Invoking Macro loop ***/
%loopy;

data monlam;
/*** Calculating randomized pvalue for lambda ***/
 set monlam;
 plow=1/(&sims+1);
 call symput('plowm',left(plow));
 opval=&lcount/(&sims+1);
run;

proc print noobs data=monlam;
/*** Examining Wilks lambda, the LRS and pvalues ***/
 var wilko lratioo df1 tpval opval;
 title 'Observed Wilks with theoretical (tpval) and randomized (opval) pvalues';
 title2 'Note that the smallest randomized p value you can obtain is' &plowm;
run;

data mont;
/*** Calculating randomized pvalue for the t-tests ***/
 array pvals{&ldim} pv1-pv&ldim;
 set mont;
 %macro tpval;
 %do j=1 %to &ldim;
    i=&j;
    pvals[i]=&&tcount&j/(&sims+1);
 %end;
 %mend;
 %tpval;
 keep t1-t&ldim prank1-prank&habs countt1-countt&ldim pv1-pv&ldim;
run;

data nicetab;
/*** Generating an organized output table ***/
 array tvals{&ldim} t1-t&ldim;
 array prank{&habs} prank1-prank&habs;
 array pvals{&ldim} pv1-pv&ldim;
 array tvec{&habs} ts1-ts&habs;
 array pvec{&habs} pvs1-pvs&habs;
 array ind{&habs} $ inds1-inds&habs;
 plow=1/(&sims+1);
 call symput('plowm',left(plow));
 set mont;
 do i=1 to &habs;
   do j=(i-1)*&habs+1 to &habs*i;
     k=j-(i-1)*&habs;
     tvec[k] = tvals[j];
     pvec[k] = pvals[j];
     if tvec[k]>=0 and pvec[k] < 0.05 then ind[k]="+++";
     if tvec[k]>=0 and pvec[k] >= 0.05 then ind[k]="+";
     if tvec[k]<0 and pvec[k] < 0.05 then ind[k]="---";
     if tvec[k]<0 and pvec[k] >= 0.05 then ind[k]="-";
     if tvec[k]=. or pvec[k]=. then ind[k]=".";
   end;
   crank=prank[i]; output;
 end;
 keep ts1-ts&habs pvs1-pvs&habs inds1-inds&habs crank;
run;

proc print noobs data=nicetab;
/*** Observed t-stats ***/
 title1 'Observed t stats';
 var ts1-ts&habs;
run;

proc print noobs;
/*** Printing randomized p values ***/
 title1 'Randomized p vals';
 title2 'Note that the smallest p value you can obtain is' &plowm;
 var pvs1-pvs&habs;
run;

proc print noobs;
/*** Examining preferences (with symbols) and the simplified ranks ***/
 var inds1-inds&habs crank;
 title1 'Simplified ranks based on observed t stats and randomized p vals';
 title2 '(Large rank indicates more preferred)';
run;


quit;