/* PaspalumRandomizationTest.sas */
/* paspalum fresh root weight example */
/* G1 inoculated data, G2 not inoculated */
/* NOTE: this program is resampling -- NOT all permutations */
proc iml;
G1={3.9, 4.3, 4.9, 5.2, 6.5, 7.6, 9.6, 10.0, 10.1, 12.3, 13.6, 19.7};
G2={6.2, 8.7, 11.0, 12.2, 12.3, 13.1, 13.6, 14.5, 15.4, 16.4, 16.7, 21.8};
obsdiff = mean(G1) - mean(G2);
* print obsdiff;          /* the observed difference between the sample means */         
* call randseed(12345);                             /* set random number seed */
alldata = G1 // G2;                          /* stack data in a single vector */
N1 = nrow(G1);  N = N1 + nrow(G2);                   /* find the sample sizes */
NRepl = 9999;                                 /* number of permutations less 1*/
nulldist = j(NRepl,1);                         /* initialize a vector of ones */
x = sample(alldata, N//NRepl, "WOR");                 /* create all resamples */
nulldist = x[,1:N1][,:] - x[,(N1+1):N][,:];   /* compute all mean differences */
/* Compute the randomization distribution of the difference between the 
   sample means, mark the location of the observed diference, and compute
   the P-value of the nondirectional hypothesis test. */
title "Histogram of Null Distribution -- Paspalum example";
title2 "The vertical line marks the observed value of the statistic. The P-value is nondirectional. 
The second row gives the t-test P-value as an approximation of the randomization P-value";
refline = "refline " + char(obsdiff) + " / axis=x lineattrs=(color=red);";
call Histogram(nulldist) other=refline;
pvalue = (1 + sum(abs(nulldist) >= abs(obsdiff))) / (NRepl+1);
print obsdiff pvalue;
title2;
/* We now use a t-test, which can be viewed as an approximation of the randomization 
   test,to find the P-value. */
sddiff = sqrt(((N1-1)*var(G1) + (N-N1-1)*var(G2))/(N-2));
sediff = sddiff*sqrt(N/(N1*(N-N1)));
tcalc = obsdiff/sediff;
abstcalc = abs(tcalc);
df=N-2;
tpvalue = 2*(1-cdf('T',abstcalc,df));
print "t-test pvalue as an approximation of randomization pvalue";
print obsdiff sediff df tcalc tpvalue;