/*  Program to simulate exact randomization
    distribution of an F-statistic for a 
    randomized block design.  This code 
    assumes that each treatment is used
    exactly once in each block */


proc iml;
start rbsim;

/* Enter the subject responses in each block
   with one row for each block. */

   re = {1.4 2.8 1.8, 
         2.6 1.4 1.6,
         3.5 1.9 2.4,
         3.4 2.2 2.8,
         4.9 3.4 3.6,
         5.9 6.1 4.4,
         2.4 1.3 1.6,
         4.3 2.4 3.7,
         1.5 2.5 1.2,
         4.6 2.2 2.9};

/* Compute number of blocks */ 

	nb = nrow(re);

/* Compute number of units per block */

	ns = ncol(re);
    nt = nb*ns;

/* Construct the model matrix and 
	projection operator */

   x = j(nt,1)||(i(nb)@j(ns,1))||(j(nb,1)@i(ns));
   xb = j(nt,1)||(i(nb)@j(ns,1));
   px = x*ginv(t(x)*x)*t(x);
   pb = xb*ginv(t(xb)*xb)*t(xb); 

    
/* Randomly sample nsim=5000 of the possible
   random assignments of subjects to 
   treatments within blocks */

   nsim=10000;
   f = j(nsim,1,0);
 
   do ii=1 to nsim;
   y=j(nt,1,0);
   nc=0;
    do i = 1 to nb;
	  rv = j(ns,1,0);
      do j = 1 to ns;
        rv[j,1]=ranuni(0);
        end;
  	  rvr=rank(rv);
      do j = 1 to ns;
	     nc=nc+1;
		 y[nc,1] = re[i,rvr[j,1]];
         end;
      end;
      yhat = px*y;
	  r = y-yhat;
	  f[ii,1] = (t(y)*(px-pb)*y/(ns-1))/((t(r)*r)/(nt-nb-ns+1));
    end;
   
  create newf from f;
    append from f; 

/* Compute f-value for original data 
	assuming each column of the data
	array corresponds to a different
	treatment */

y = shape(re,nt,1);
  yhat = px*y;
  r = y-yhat;
  mse = (t(r)*r)/(nt-nb-ns+1);
  mst = t(y)*(px-pb)*y/(ns-1);
  fvalue = (t(y)*(px-pb)*y/(ns-1))/((t(r)*r)/(nt-nb-ns+1));
    
/* Compute p-value */

pvalue=nsim;
do i = 1 to nsim;
  if(f[i,1] < fvalue) then pvalue=pvalue-1;
  end;

pvalue=pvalue/nsim;

print fvalue pvalue;

finish;

run rbsim;


proc univariate data=newf; 
  var col1;   run;


goptions  device=WIN target=WINPRTC;

  axis1 label = (h=2 r=0 a=90 f=swiss 'Percent')
        length = 5 in
        value = (h=2.0 f=swiss);

 axis2 label = (h=2 f=swiss 'F-values')
        length = 5 in
        value = (h=2.0 f=swiss);


  proc gchart data=newf;
     vbar col1 / type=percent space=0 frame
                  raxis = axis1 gaxis=axis2;
     title h=3.0 f=swiss "Randomization F-values";
     run;


 data set1;
   input block diet y;
   datalines;
   1 1 1.4
   1 2 2.8
   1 3 1.8
   2 1 2.6
   2 2 1.4
   2 3 1.6
   3 1 3.5
   3 2 1.9
   3 3 2.4
   4 1 3.4
   4 2 2.2
   4 3 2.8
   5 1 4.9
   5 2 3.4
   5 3 3.6
   6 1 5.9
   6 2 6.1
   6 3 4.4
   7 1 2.4
   7 2 1.3
   7 3 1.6
   8 1 4.3
   8 2 2.4
   8 3 3.7
   9 1 1.5
   9 2 2.5
   9 3 1.2
   10 1 4.6
   10 2 2.2
   10 3 2.9
   run;

   proc glm data=set1;
     class block diet;
	 model y = block diet;
	 run;