proc iml;
/*----------------------------------------------------*/
/* The SIMSIZE Module */
/* */
/* This module computes the simulation size necessary */
/* to estimate a binomial proportion p within a given */
/* bound of error 99% of the time. This module is */
/* useful when p is the significance level of a */
/* simulation-based test. */
/* */
/* Input: */
/* 1) SIGLEV, a number in [0,1] that you think is */
/* close to the true binomial proportion (or */
/* significance level). If you have no idea, */
/* then you should use 0.5. */
/* 2) BOUND, your desired bound of error -- the */
/* amount added & subtracted to obtain lower */
/* and upper 99% confidence limits. */
/* */
/* Returns: */
/* 1) NSS, a positive integer that is the */
/* necessary simulation size. */
/*----------------------------------------------------*/
start simsize(siglev,bound, nss);
reset noname;
z=probit(0.995);
nss=ceil(siglev#(1-siglev)#(z/bound)##2);
print nss;
reset fw=2;
print "is the simulation size necessary to obtain a value within";
reset fw=5;
print bound "of" siglev "99% of the time";
finish;
start samp(n,seed, rvec);
rvec=j(1,n,.);
u=-99;
do j = 1 to n;
call ranuni(seed,u);
rvec[1,j]=u;
end;
rvec=rank(rvec);
finish;
start order(mat0,n,k, mat1);
d=mat0[,1:k];
nvec=d[+,];
add=j(n,1,.);
do i = 1 to n;
ind=d[i,<:>];
add[i,1]=10#nvec[1,ind]+ind;
end;
mat1=mat0||add;
mat2=mat1;
mat2[rank(add),]=mat1;
mat1=mat2[,1:k+1];
finish;
start getstat(test,n,k,d,r,nvec,cont,stderr, w);
if test='KW' then w=sum(((r`*d)##2)/nvec);
if test='JO' then do;
dr=hdir(d,r);
w=0;
do j = 1 to k-1;
drj=dr[,j];
do h = j+1 to k;
drh=dr[,h];
mannwitt=0;
do i = 1 to n;
mannwitt=mannwitt+
sum((drj[i,1]>1e-2 & drj[i,1]1e-2 & drj[i,1]=drh));
end;
w=w+mannwitt;
end;
end;
end;
if test='MC' then do;
mean=cont*((d`*r)/nvec`);
w=abs(mean/stderr);
end;
finish;
/*----------------------------------------------------*/
/* The ONEWAY Module */
/* */
/* This module obtains simulation-based significance */
/* levels with 99% confidence limits for either the */
/* Kruskal-Wallis test or the Jonckheere test in the */
/* one-way layout. */
/* */
/* Input: */
/* 1) TEST, a character string that denotes which */
/* test you want to perform. The choices are */
/* 'KW' (for Kruskal-Wallis) and 'JO' (for */
/* Jonckheere). */
/* 2) NN, a positive integer specifying desired */
/* simulation size. */
/* 3) MAT, the Nx2 data matrix -- the 1st column */
/* specifies numeric class variable levels, the */
/* 2nd specifies the dependent variable values. */
/* 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. */
/*----------------------------------------------------*/
start oneway(test,nn,mat,seed);
reset noname fw=6;
test=upcase(test);
pval=0; eps=1e-12;
n=nrow(mat);
matu=design(mat[,1])||ranktie(mat[,2]);
k=ncol(matu)-1;
run order(matu,n,k, mat0);
d=mat0[,1:k]; r=mat0[,k+1];
nvec=d[+,];
run getstat(test,n,k,d,r,nvec,.,., w0);
do i = 1 to nn;
run samp(n,seed, rvec);
run getstat(test,n,k,d,rvec`,nvec,.,., w);
if w+eps>=w0 then pval=pval+1;
end;
pval=pval/nn;
z=probit(0.995);
se=sqrt((pval#(1-pval))/nn);
inc=z#se;
l=max(pval-inc,0); u=min(pval+inc,1);
table=nn||pval||l||u;
chead={"Sim. Size" "P-Value" "99% LCL" "99% UCL"};
rhead=test;
print table[rowname=rhead colname=chead];
finish;
/*----------------------------------------------------*/
/* The MULTCOMP Module */
/* */
/* This module obtains simulation-based significance */
/* levels for multiple comparisions with an overall */
/* type I error rate in the K-sample one-way layout. */
/* */
/* Input: */
/* 1) NN, a positive integer specifying desired */
/* simulation size. */
/* 2) MAT, the Nx2 data matrix -- the 1st column */
/* specifies numeric class variable levels, the */
/* 2nd specifies the dependent variable values. */
/* 3) CONT, a QxK matrix defining the Q linear */
/* combinations of interest */
/* 4) SEED, a nonnegative integer 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. */
/*----------------------------------------------------*/
start multcomp(nn,mat,cont,seed);
reset noname fw=6;
q=nrow(cont);
pval=j(q,1,0); eps=1e-12;
n=nrow(mat);
matu=design(mat[,1])||ranktie(mat[,2]);
k=ncol(matu)-1;
run order(matu,n,k, mat0);
d=mat0[,1:k]; r=mat0[,k+1];
nvec=d[+,];
stderr=sqrt(vecdiag(cont*diag(1/nvec)*cont`));
run getstat('MC',n,k,d,r,nvec,cont,stderr, w0);
do i = 1 to nn;
run samp(n,seed, rvec);
run getstat('MC',n,k,d,rvec`,nvec,cont,stderr, w);
crit=max(w);
pval=pval+((crit+eps)>=w0);
end;
pval=pval/nn;
ch1={"Contrast"}; ch2={"P-Value"};
print cont[colname=ch1] pval[colname=ch2];
finish;