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

%Let sims=999; * input number of simulations;
options nodate ls=78 ps=100 nonumber formdlim=' ' 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;

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;

data ranmat;
/*** Randomizing the antisymmetric matrix many times ***/
 array compo{&habs,&habs};
 array random{&habs,&habs};
 merge obspref;
 retain seed1 1613217181;
 do sim=1 to &sims;
   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;
   output;
 end;
 drop i j seed1;
run;

proc sort data=ranmat out=sranmat;
 by sim;
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 sim;
run;

proc sort data=ranvec out=sranvec;
 by sim;
run;

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

proc glm data=sranvec;
/*** Doing a ton of MANOVAs ***/
 model n1-n&sdim = / ss3 nouni;
 manova H=intercept;
 by sim;
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=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;

proc means data=sranmat t noprint;
/*** Doing the t-tests many times ***/
 var random1-random&ldim;
 output out=ranmean t=rt1-rt&ldim;
 by sim;
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;

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;

data monlam;
/*** Calculating the randomized pvalues for lambda ***/
 plow=1/(&sims+1);
 call symput('plowm',left(plow));
 set obslam;
 do sim=1 to &sims;
   retain count 1;
   set ranlam end=last;
   if lratio >= lratioo then count = count+1;
    opval=count/(&sims+1);
    if last then output;
 end;
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 the randomized pvalues for the t-tests ***/
 array obsv{&ldim} t1-t&ldim;
 array ts{&ldim} rt1-rt&ldim;
 array pv{&ldim} pv1-pv&ldim;
 array count{&ldim};
 set obsrank;
 do sim = 1 to &sims;
   retain count 1;
   set ranmean end=last;
   do i=1 to &ldim;
     if abs(ts[i]) >= abs(obsv[i]) then count[i] = count[i]+1;
       pv[i]=count[i]/(&sims+1);
     if obsv[i]=. then pv[i]=.;
   end;
   if last then output;
 end;
 keep t1-t&ldim pv1-pv&ldim prank1-prank&habs;
run;

data nicetab;
/*** Generating an organized output table ***/
 array tvals{&ldim} t1-t&ldim;
 array pvals{&ldim} pv1-pv&ldim;
 array prank{&habs} prank1-prank&habs;
 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;