options nodate;
%let seed=2345;
%let nsims=500;
*%let nsims=2000;  * takes a while;
%let sigma=4;
*%let sigma=2;

data diets;
   *length treatment $19.;
   *array tarray{5} $19 ("Control" "Control+level-1-A" "Control+level-2-A" "Control+level-1-B" "Control+level-2-B");
   array tmean{5} (9,15,15,15,19);
   *array tmean{5} (12,15,15,15,16);
   do simnumber=1 to &nsims;
   do treatment=1 to 5; 
   do turkey=1 to 6; 
      *treatment=tarray{i};
      wtgain = tmean{treatment} + &sigma*rannor(&seed);
      output;
   end;
   end;
   end;
   keep wtgain treatment turkey simnumber; 
run;
proc mixed;
   title "Analysis of first simulated dataset";
   where simnumber=1;
   class treatment;
   model wtgain=treatment;
   lsmeans treatment/adj=tukey;
run;
ods listing close;  /* turns off output, we only want to see summary of all simulations */
proc mixed;
   by simnumber;
   class treatment;
   model wtgain=treatment;
   lsmeans treatment/adj=tukey;
   ods output diffs=mydiffs tests3=t3;  /* creates two new datasets with pairwise comps and overall ftests */
run;

ods listing ;  /* turns output back on */


proc print data=t3 (obs=10); 
   title "t3 dataset (contains) fratios and pvalues for overall test  of equality of treatment means";
run;

data mydiffs;
   length comparison $3;   /* need to specify length 3 or may have trailing blank spaces */
   set mydiffs;
   lsd_sig = (probt < 0.05);
   tukey_sig = (adjp < 0.05);
   comparison=trim(left(treatment))||"-"||left(_treatment);
   *if comparison in ('2-3','2-4','3-4') then equal=1; *else equal=0;
   if comparison in ('2-3','2-4','3-4');   *these comparisons involve EQUAL means;
   keep comparison simnumber estimate stderr tvalue probt adjustment adjp lsd_sig tukey_sig;
run;
proc print data=mydiffs (obs=10); 
   title "mydiffs dataset";
run;

proc sort data=mydiffs;  
   by simnumber;
run;
proc means data=mydiffs noprint;
   by simnumber ;
   var lsd_sig tukey_sig;  /* compute sum of false postives for each simulation */
   output out=sigresults sum=lsd_sig tukey_sig;
run;

data results;
   merge sigresults t3;
   by simnumber;
   fsig=(probf<.05);
   label  lsd_sig="false positive count using Fisher's LSD";
   label  tukey_sig="false positive count using Tukey's HSD";
   *keep simnumber equal lsd_sig tukey_sig _freq_ fsig probf comparison;
   keep simnumber lsd_sig tukey_sig _freq_ fsig probf comparison;
run;

proc print data=results(obs=10);
   title  "results dataset";
run;

proc sort data=results;
   *by fsig equal;
   by fsig ;
run;

*options ls=125;

proc freq data=results;
   *by fsig;
   *tables lsd_sig*comparison;
   *tables tukey_sig*comparison;
   tables lsd_sig*fsig/norow;
   tables tukey_sig*fsig/norow;
run;