options ls=75 nodate;

%let seed=0;
*%let sigmaB=2;  /* block effect varcomp=4 */
%let sigmaB=2.5;  /* block effect varcomp=9 */
*%let sigmaB=3;  /* block effect varcomp=9 */
*%let sigmaB=4;  /* block effect varcomp=16 */
%let sigma=2.5;  /* error varcomp=6 */
%let numsims=100;

%macro datagen;

data one; 
   array trt t1-t4 (36,17,24,26);
   do simindex=1 to &numsims;
      do subject=1 to 9;
         block=normal(&seed)*&sigmaB;
         do ivsystem=1 to 4;
            Eij=normal(&seed)*σ
            truemean=trt{ivsystem};
            y=truemean+block+Eij;
            output;
         end;
      end;
   end;
run;

ods trace on;
proc mixed ;
   by simindex;
   class ivsystem block;
   model y=ivsystem block;
   *model y=ivsystem/ddfm=satterth; 
   *random block;
   lsmeans ivsystem/cl;
   ods output lsmeans=lsm;

proc print data=lsm;run;
proc append base=two data=lsm;
run;
%mend datagen;

%datagen;

data two;
   array trt t1-t4 (36,17,24,26);
   set two;
   truemean=trt{ivsystem};
   tstat=abs(estimate-truemean)/Stderr;
   pvalue=2*(1-probt(tstat,DF));
run;

proc sort;  by ivsystem pvalue;run;

proc print data=two;
   title "output from PROC MIXED macro";
   var simindex ivsystem truemean Estimate StdErr DF lower upper tstat pvalue;
run;

proc univariate;
   histogram pvalue;
run;