proc iml;
start samp(b,t,u,rmat0,seed, rmat);
rmat=j(b,t,.);
do i = 1 to b;
call ranuni(seed,u);
u=rank(u);
rmat[i,]=(rmat0[i,])[,u];
end;
finish;
start check(dsin,test,simsize,seed,cmat, b,t,rmat0,tu,flag);
flag=0;
if test^='FR' & test^='PA' & test^='MCA' &
test^='MCC' & test^='MCU' then do;
print "ERROR: Invalid TEST Choice";
flag=1;
return;
end;
if simsize<1 | floor(simsize)^=simsize then do;
print "ERROR: SIMSIZE Must Be A Positive Integer";
flag=1;
return;
end;
if seed<0 | floor(seed)^=seed then do;
print "ERROR: SEED Must Be A Nonnegative Integer";
flag=1;
return;
end;
call execute ("use ", dsin, " ;");
read all var{block} into bvec;
read all var{treatmnt} into tvec;
read all var{y} into y;
bu=unique(bvec);
tu=unique(tvec);
b=ncol(bu);
t=ncol(tu);
n=nrow(y);
if test='MCU' then do;
if cmat=. then do;
print "ERROR: Must Supply CMAT For MCU Test";
flag=1;
return;
end;
if ncol(cmat)^=t then do;
print "ERROR: # of Cols For CMAT Must = # of Trts.";
flag=1;
return;
end;
end;
if n^=b#t then do;
print "ERROR: Not An RCB Design";
flag=1;
return;
end;
if max(y)-min(y)<1e-12 then do;
print "ERROR: Y Is A Constant Value";
flag=1;
return;
end;
mat0=j(b,t,.);
do i = 1 to n;
bind=(abs(bu-bvec[i,1]))[,>:<];
tind=(abs(tu-tvec[i,1]))[,>:<];
mat0[bind,tind]=y[i,1];
end;
do h = 1 to b;
rmat0=rmat0//ranktie(mat0[h,]);
end;
if min(rmat0)=max(rmat0) then do;
print "Between-Treatment Variability Is 0";
print "All Tests Are (Trivially) Non-Significant";
flag=1;
return;
end;
finish;
start getstat(test,b,t,rmat,rbar,bot,ceemat, w);
rjvec=rmat[+,];
if test='FR' then w=(12#ssq(rjvec-b#rbar))/bot;
if test='PA' then w=rjvec*(1:t)`;
if test='MCA' | test='MCC' | test='MCU' then
w=abs(ceemat*rjvec`);
finish;
start getmat(t,o);
eye=i(t);
do j = 2 to o;
add=eye[j:t,]-eye[j-1,]@j(t-j+1,1,1);
cmat=cmat//add;
end;
return(cmat);
finish;
/*--------------------------------------------------------*/
/* The NPARRCB Module */
/* */
/* This module generates simulation-based significance */
/* levels and confidence intervals for five different */
/* tests in the Nonparametric Randomized Complete Block */
/* (NRCB) setting. The five tests are */
/* */
/* 1) Friedman's test for general treatment differences */
/* 2) Page's test for ordered treatment differences */
/* 3) Multiple Comparisons -- all pairwise differences */
/* 4) Multiple Comparsions -- all treatments vs. control */
/* 5) Multiple Comparsions -- user-specified differences */
/* */
/* Input: */
/* */
/* 1) DSIN -- input data set that MUST contain the */
/* three numeric variables TREATMNT, BLOCK, */
/* and Y (treatment values, block values, */
/* and dependent variable values). */
/* */
/* 2) TEST -- a character string that specifies the */
/* test to be performed. The choices are */
/* FR (Friedman), PA (Page), MCA (all */
/* pairwise diffs), MCC (treatments vs. */
/* control), MCU (user-defined diffs). */
/* */
/* 3) SIMSIZE -- a positive integer specifying the */
/* desired simulation size. */
/* */
/* 4) SEED -- a nonnegative ingeger used in generating */
/* the simulated data. If SEED is set to */
/* 0, then the clocktime is used. SEED */
/* should be non-zero when you want to be */
/* able to reproduce your results. */
/* */
/* 5) CMAT -- set to missing unless TEST='MCU'. If */
/* TEST='MCU', CMAT is a qxt matrix */
/* defining the q treatment comparisons */
/* of interest. */
/*--------------------------------------------------------*/
start nparrcb(dsin,test,simsize,seed,cmat);
reset noname fw=6;
test=upcase(test);
run check(dsin,test,simsize,seed,cmat, b,t,rmat0,tu,flag);
if flag=1 then return;
eps=1e-12;
bot=.;
rbar=.;
k=1;
if test='FR' then do;
print "Friedman Test with Simulation Size = "
simsize;
pval=0;
rbar=(t+1)/2;
bot1=b#t#(t+1);
bot2=0;
do i = 1 to b;
ui=unique(rmat0[i,]);
do j = 1 to ncol(ui);
tie=tie//sum(rmat0[i,]=ui[1,j]);
end;
bot2=bot2+sum(tie##3);
free tie;
end;
bot2=(bot2-b#t)/(t-1);
bot=bot1-bot2;
end;
if test='PA' then do;
print "Page Test with Simulation Size = "
simsize;
pval=0;
end;
ch1={"P-Value" "99% LCL" "99% UCL"};
if test='MCA' | test='MCC' | test='MCU' then do;
if test='MCA' then do;
dum=t;
ceemat=getmat(t,dum);
descrip={"(All Pairwise)"};
end;
if test='MCC' then do;
ceemat=getmat(t,2);
descrip={"(All To Control)"};
end;
if test='MCU' then do;
ceemat=cmat;
descrip={"(User-Defined)"};
end;
print "Multiple Comparisons" descrip;
print "Simulation Size = " simsize;
k=nrow(ceemat);
pval=j(k,1,0);
end;
run getstat(test,b,t,rmat0,rbar,bot,ceemat, w0);
print "Test Statistic(s): W = " w0;
u=j(1,t,-99);
do i = 1 to simsize;
run samp(b,t,u,rmat0,seed, rmat);
run getstat(test,b,t,rmat,rbar,bot,ceemat, w);
pval=pval+(max(w)+eps>=w0);
end;
pval=pval/simsize;
z=probit(0.995);
se=sqrt((pval#(1-pval))/simsize);
inc=z#se;
zer=j(k,1,0);
one=j(k,1,1);
lower=((pval-inc)||zer)[,<>];
upper=((pval+inc)||one)[,><];
table=pval||lower||upper;
if test='FR' | test='PA' then print
table[colname=ch1 format=6.4];
if test='MCA' | test='MCC' | test='MCU' then do;
q=nrow(ceemat);
rh1={"Contrast"};
if q>1 then rh1=rh1||j(1,q-1," ");
pre=j(1,t,"T");
ch2=compress(concat(pre,char(tu,4)));
print ceemat[rowname=rh1 colname=ch2]
table[colname=ch1 format=6.4];
end;
finish;