/*************************************************** 
If you just want to use this program for analysis, then all you
will need to do is look at the top part of the program (the
USER ENTRY SECTION):

-You can put your own title after SOURCE='....'
-Enter the number of factors as K=...
-Enter the data in standard (Yates) order as Y0=(....)
-Pick off a decision rule P-level from Table 1 in the paper and enter as
  P0=...
-Enter the number of permutations as JMAX=...  (Recommend at least 5000)
-Enter seed for the random number generator as SEED=...
  (1<= SEED <= 2**31 - 1, or use 0 for internally-generated seed)

GOOD LUCK!
***************************************************/ 

options ps=500 ls=78 nonumber nodate;
dm 'zoom off; log; clear; out; clear; ';
/* PERMUTATION TEST FOR SIGNIFICANT EFFECTS IN */
/* NONREPLICATED FACTORIAL DESIGNS */
/* LAST REVISION: THURSDAY, SEPT 14, 2000 */
/* (Fixed HTML glitch that ate several lines) */

proc iml;
file print;


/* =================USER ENTRY SECTION============================ */
/* "source"=DESCRIPTION OF EXAMPLE */
source='Montgomery Problem 9.11(Totals); EER=.20';

/* "k"=NUMBER OF VARIABLES */
k=4;

/* "y0"=ORIGINAL DATA IN STANDARD ORDER */
y0={3.62, 2.90, 2.88, 3.22, 2.61, 2.51, 2.88, 2.56, 4.23, 3.71, 3.74, 3.01, 
3.73, 2.63, 2.99, 2.73};

/* "p0"=DECISION RULE p-VALUE */
p0=0.135;

/* "jmax"=NO. OF SIMULATIONS TO GENERATE PERMUTATION DIST */
jmax=50;

/* "seed"=SEED FOR RANDOM NUMBER GENERATOR */
seed=1578585532;
/* ==============END OF USER ENTRY SECTION======================== */


put @1 'The source of the data is: ' source '.';
put @1 'This is a 2^' k 1.0 ' factorial design.';
put @1 'The permutation test uses ' jmax 5.0 ' random permutations.';
put @1 'The decision rule p-value is ' p0 5.3 '.';
put @1;

/* "n"=NUMBER OF RUNS */
n=2**k;

/* "z"=DESIGN MATRIX */
pi=3.14159265359;
z0=j(n,1,1);
z1=-sin((pi*(do(0,n-1,1) + 0.5)`)*(2##(-do(0,k-1,1))));
z=z0;
do j=1 to k;
 z=z||(z#z1[,j]);
end;
z=z/abs(z);

/* "id"=PARAMETER I.D. */
id=0;
do j=1 to k;
 id=id//(id*10 + j);
end;
*print id y0;

/* "bhat"=PARAMETER ESTIMATES FOR ORIGINAL DATA */
bhat=solve(z`*z,z`*y0);

/* "absrank"=RANKS OF ABSOLUTE ESTIMATES */
/* 1=MEAN, 2=SMALLEST, ... , n=LARGEST */
absrank={1}//(1+rank(abs(bhat[2:n])));

/* "p"=VECTOR OF p-VALUES */
p=j(n,1,1);

/* MACRO: "test(y)" GENERATES p-VALUE FOR RESPONSE VECTOR Y,
   "bhatmax"=MAXIMUM ABSOLUTE BHAT,
   "count"=NUMBER OF TIMES BHATMAX EXCEEDED
   "yperm"=VECTOR TO HOLD PERMUTED DATA */
%macro test(y);
bhatmax=max(abs(remove((1/n)*z`*&y,1)));
count=0;
do j=1 to jmax;
 yperm=&y[rank(uniform(j(n,1,seed)))];
 if sqrt((n-1)/(i-1))*max(abs(remove((1/n)*z`*yperm,1))) >= bhatmax
   then count=count+1;
end;
p[loc(absrank=i)]=1-(1-count/jmax)##((i-1)/(n-1));
*p[loc(absrank=i)]=count/jmax;
%mend test;

/* SYSTEMATICALLY REMOVE THE LARGEST ABSOLUTE ESTIMATES
   AND CALCULATE p-VALUE FOR EACH NEW MODEL. */
y=y0;
do i=n to 2 by -1;
 %test(y);
 i1=i-2;
 put @1 i1 2.0 ' more parameter(s) to go.';
 y=y-bhat[loc(absrank=i)]*z[,loc(absrank=i)];
end;

/* "lnpe"=LARGEST NONSIGNIFICANT PARAMETER ESTIMATE */
sig=p < p0;
if sum(sig) = 0 then lnpe=n;
if sum(sig) > 0 then lnpe=min(absrank[loc(sig)])-1;

/* "decision"=1 IF SIGNIFICANT, "decision"=0 IF NOT SIGNIFICANT */
decision=absrank>lnpe;

/* PRINT UNSORTED RESULTS */
print ('UNSORTED RESULTS (DECISION=1 MEANS SIGNIFICANT)'//source);
print bhat id absrank p[format=6.4] decision;

/* PRINT SORTED RESULTS */
id1=id;
id[absrank]=id1;
bhat1=bhat;
bhat[absrank]=bhat1;
p1=p;
p[absrank]=p1;
decis1=decision;
decision[absrank]=decis1;
hnperc={0}//probit((do(1,n-1,1)-0.5)/(2*n) + 0.5)`;

print ('SORTED RESULTS (DECISION=1 MEANS SIGNIFICANT)'//source);
print bhat id p[format=6.4] decision hnperc;

/* MAKE HALF-NORMAL PLOT OF PARAMETER ESTIMATES */
options ps=53 ls=78 nonumber nodate;
call pgraf(abs(bhat[2:n])||hnperc,cshape('*',1,n-1,1),
           'ABSOLUTE PARAMETER ESTIMATES','HALF-NORMAL PERCENTILES',
           concat('HALF-NORMAL PLOT: ',source));
call pgraf(abs(bhat[2:n])||hnperc,char(id[2:n],k),
           'ABSOLUTE PARAMETER ESTIMATES','HALF-NORMAL PERCENTILES',
           concat('HALF-NORMAL PLOT: ',source));