/* EnergyPermutationTest.sas */
/* energy usage example */
/* part 1: houses of either design type */
/* G1 extra insulation data, G2 standard insulation data */
/* NOTE: this program is resampling -- NOT all permutations */
proc iml;
/* extra */
G1={8.3, 11.7, 12.7, 13.0,13.4, 13.6, 13.7, 13.7, 13.8, 14.6, 15.3, 15.6, 16.0, 18.8}; 
/*standard */
G2={11.4, 13.9, 13.9, 14.0, 15.3, 18.0, 18.0, 18.1, 18.9, 19.0, 19.0, 21.7};
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 -- Energy example either design";
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 permutation 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 permutation pvalue";
print obsdiff sediff df tcalc tpvalue;