/*  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 1.9, 
         2.6 1.4 1.6 1.5,
         3.5 1.9 2.4 5.6,
         3.4 2.2 2.8 2.4,
         4.9 3.4 3.6 3.5,
         5.9 6.1 4.4 4.7,
         2.4 1.3 1.6 1.5,
         4.3 2.4 3.7 3.1,
         1.5 2.5 1.2 1.4,
         4.6 2.2 2.9 3.1};

/* 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=10000 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;

/* p-value from central F-distribution */

pvaluef = 1-probf(fvalue,ns-1,(ns-1)*(nb-1));

print pvaluef;
finish;

run rbsim;


proc univariate data=newf; 
  var col1;   run;


  /*  Trim some extreme values before making a histogram */

data newf; set newf;
  if(col1<0 or col1>9) then delete;
  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;