option ls=76 nocenter;

%let seed      = 71853048;
%let var_cage  = 5.03;
%let var_cc    = -6.24;
%let var_resid = 31.56; 
%let nsim      = 1000;
%let condm1    = 56; %let condm2    = 53;     /* Input the condition means */
%let condm3    = 59; %let condm4    = 60;
                      /* A: To generate the cc and residual covariance matrix cov
                            and its Cholesky decomposition half(cov)              */
proc iml;
   cov = &var_cc * j(3,3) + &var_resid * i(3);
   covh = half(cov)`;
   print cov covh;
run;
	             /* END A  */			
data noise;
   retain i j;
   do sim = 1 to ≁
      i = 1;                                 /* i = condition index                */
      j = 2;                                 /* j = cage index                     */
      do cage = 1 to 6;                      /* cage is the index k in the model   */
         u_cage = sqrt(&var_cage) * rannor(&seed);  /* For the cage component      */ 
         do condition = i,j;                 /* B: For the cc and residual component */
            e1 = rannor(&seed);
            e2 = rannor(&seed);
            e3 = rannor(&seed);
            diet = "normal  ";
            noise = u_cage + 5.0318983*e1;
            output;
            diet = "restrict";
            noise = u_cage - 1.240089*e1 + 4.8766977*e2;
            output;
            diet = "suppleme";
            noise = u_cage - 1.240089*e1 - 1.594895*e2 + 4.6085237*e3;
            output;	
         end;                                       /* END B  */ 
         j = j + 1;                          /* Here is to arrange the incomplete block index */
         if (j = 5) then do;                 /* see output4.25.sas input data                 */
            i = i + 1;
            j = i + 1;
         end;
      end;
   end;
run;
               /* End of the simulation. It contains &nsim (=1000) simulated experiments      */       
               /* How does the output look like? See output12_6a.sas */

DATA non_null;SET noise;                                    /* Create non_null data */
y=noise;
IF condition=1 THEN y=y+&condm1;
IF condition=2 THEN y=y+&condm2;
IF condition=3 THEN y=y+&condm3;
IF condition=4 THEN y=y+&condm4;
  
ods exclude all;      /* These two ods statements are for the next proc mixed  */ 
ods noresults;
                                     
proc mixed data=non_null nobound;   /* Note: DATA=non_null                */ 
   by sim;
   class cage condition diet;
   model y=condition|diet;
   random cage cage*condition;
                                     /*parms &var_cage &var_cc &var_resid;                         */      
                                     /* These parms will help the iteration process in computation */
                                     /* It is unnecessary. See output12_6b.sas                     */ 
   ods output covparms=cpsim tests3=t3sim;
run;

ods exclude none;
ods results;

proc sort data=cpsim;
   by covparm;
run;

proc means data=cpsim;   /* output the mean variance components */
   by covparm;
   var estimate;
   TITLE 'Output the variance components'; /* added statement */
run;

data t3sim;
   set t3sim;
   if probf < 0.05 then sig = 1;
   else sig = 0;
proc sort data=t3sim;
   by effect;
proc means data=t3sim mean;
   by effect;
   var sig;
   TITLE 'Output the number of rejection under the null'; /*added statement */ 
run;

PROC PRINT DATA=non_null (obs=10);VAR cage condition diet noise y;