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