data allel;options ls=130;
PROC IML WRKSPACE=400;
**************************************************;
*  Example for J = 3 groups & k = 5 variables  ***;
**************************************************;
*******      Set Number of Groups J  *************;
J=3;k=5;
*******      Set Sample Sizes n ******************;
n1= 23;n2= 54;n3= 66;nv=n1//n2//n3;nt=n1+n2+n3;
**************************************************;
***  Set Sample Variances  v for each of the k variables ***********;
v1={8 2 4 3 6};v2={4 4 4 3 2};v3={2 3 4 3 1};
**************************************************;
ALP={.05 .01 .001};**** the Alpha values to be tested *****;
**************************************************;
***  Set Sample Correlation Matrix for each group ***********;
r1={1.0 0.4 0.7 0.6 0.5,
    0.4 1.0 0.5 0.3 0.2,
    0.7 0.5 1.0 0.6 0.4,
    0.6 0.3 0.6 1.0 0.6,
    0.5 0.2 0.4 0.6 1.0};
r2={1.0 0.4 0.7 0.6 0.5,
    0.4 1.0 0.5 0.3 0.2,
    0.7 0.5 1.0 0.6 0.4,
    0.6 0.3 0.6 1.0 0.6,
    0.5 0.2 0.4 0.6 1.0};
r3={1.0 0.4 0.7 0.6 0.5,
    0.4 1.0 0.5 0.3 0.2,
    0.7 0.5 1.0 0.6 0.4,
    0.6 0.3 0.6 1.0 0.6,
    0.5 0.2 0.4 0.6 1.0};
**************************************************;
k=ncol(r1); ********* k = NUMBER OF VARIABLES   *************;
x1=j(n1,k,0);x2=j(n2,k,0);x3=j(n3,k,0);
fm=j(k,1,0);ones=j(k,1,1);
** Initialize other varibles    ******************;
rep=10000;**** Number of Replications ************;
typeI=j(1,3,0);FAMSIZE=j(1,3,0);WFAMSIZE=j(1,3,0);
wtypeI=j(1,3,0);adjalp=j(1,3,0);wadjalp=j(1,3,0);
rlab={Alpha, Emprical, Adjusted, C};
s1=v1##0.5;s2=v2##0.5;s3=v3##0.5;
**************************************************;
do rj=1 to rep;
do jj=1 to J;
if jj=1 then rx=r1;
if jj=1 then nn=n1;
if jj=1 then sd=s1;
if jj=2 then nn=n2;
if jj=2 then rx=r2;****Group 2 has a different correlation imposed **;
if jj=2 then sd=s2;****Group 2 has a different variance imposed **;
if jj=3 then nn=n3;
if jj=3 then rx=r3;****Group 3 has a different correlation imposed **;
if jj=3 then sd=s3;****Group 3 has a different variance imposed **;

*************************************************************;
*** PERFORM PRINCIPAL COMPONENTS ANALYSIS ON (rxx) **********;
*************************************************************;
lam=eigval(rx);***  LATENT ROOTS OF rxx  *******************;
lamrt=lam##0.5;
lsqrt=diag(lamrt);**  DIAGONAL MATRIX WITH THE SQUARE ROOT OF*;
********************                    THE EIGENVALUES  *******;
ev=eigvec(rx);** EIGENVECTORS OF rxx  *************************;
** CREATE FACTOR SCORE MATRIX (fr) BY MULTIPLYING EACH     ******;
**     EIGENVECTOR BY THE SQUARE ROOT OF ITS EIGENVALUE   ******;
           fr=ev*lsqrt;
** GENERATE A DATA MATRIX (x) OF RANDOM NORMAL VARIABLES  ******;
xx=rannor(j(nn,k,0));** 0 IS A CLOCK GENERATED PSEUDORANDOM SEED **;
******************************************************************;
*****************************************************************;
** TRANSFORM THE UNCORRELATED RANDOM NORMAL VARIABLES TO    *****;
** CORRELATED VARIBLES THROUGH EQUATION (1)    ******************;
**  fr MULTIPLIED BY THE TRANSPOSE OF X THEN RE-TRANSPOSE *******;
zh=fr*xx`;xx=t(zh);
************************************************************;
do ii=1 to k;
  xx(|,ii|)=xx(|,ii|)#sd(|,ii|);
end;
if jj=1 then x1=xx;
if jj=2 then x2=xx;
if jj=3 then x3=xx;
end;
************************************************;
xtot=x1//x2//x3;
GM=xtot(|:,|);** Calculate Grand mean ****;
***   Calculate group Means   ******************;
xm1=x1(|:,|);xm2=x2(|:,|);xm3=x3(|:,|);
XM=xm1//xm2//xm3;
************************************************;
***   Calculate Within-Group Sum of Squares  ***;
ss1=x1-repeat(xm1,n1,1);ss1=ss1##2;ss1=ss1(|+,|);
ss2=x2-repeat(xm2,n2,1);ss2=ss2##2;ss2=ss2(|+,|);
ss3=x3-repeat(xm3,n3,1);ss3=ss3##2;ss3=ss3(|+,|);
ssm=ss1//ss2//ss3;
************************************************;
**** Calculate F-ratios **********************;
dvv=nv`*((xm-repeat(GM,J,1))##2);
dfd=nt-J;*** Denominator degrees of freedom  ****;
FM= (dvv/(J-1))/(ssm(|+,|)/(dfd));
*****************************************************************;
Fmax=MAX(FM);*** Find the Maximum F to estimate the probability ***;
****************        of at least one Type I error            ***;
pmax=1-probf(Fmax,2,dfd);** Calculate p-values  *******;
if pmax < .05 then typeI(|1,1|)= typeI(|1,1|) + 1;*** Accumulate    *;
if pmax < .01 then typeI(|1,2|)= typeI(|1,2|) + 1;*** Rejections of *;
if pmax < .001 then typeI(|1,3|)= typeI(|1,3|) + 1;** the Null Hyp. *;
*****************************************************************;
****** Calculate Group Variance for calculate of Welch Test ***;
vv1=ss1/(n1-1);vv2=ss2/(n2-1);vv3=ss3/(n3-1);
vvm=vv1//vv2//vv3;wwm=ssm;
***** Calculate adjusted grand mean for Welch Test ********;
do ii=1 to k;
wwm(|,ii|)=nv/vvm(|,ii|);end;
wg=wwm#xm;WG=wg(|+,|);WG=wg/(wwm(|+,|));
*****************************************************;
***  Calculate Welch Test       *********************;
dvw=wwm#((xm-repeat(WG,3,1))##2);
dvw=dvw(|+,|);
delta= (((wwm / (repeat(wwm(|+,|),3,1)))-1)##2);
delta=(delta / (repeat((nv-1),1,k)));
delta=delta(|+,|);delta=(3/8)#delta;dfww=1/delta;
www=(dvw/2)/(((2/3)#delta)+1);
wpval=www;
********* Calculate p-values for Welch Test ***************;
do ii=1 to k;
wpval(|1,ii|)=1 - (probf(www(|1,ii|),2,dfww(|1,ii|)));
end;
******************************************************************;
WMAX=min(wpval);*** Find minimum p-value for Welch Test  **********;
if wmax < .05 then wtypeI(|1,1|)= wtypeI(|1,1|) + 1;*** Accumulate    *;
if wmax < .01 then wtypeI(|1,2|)= wtypeI(|1,2|) + 1;*** Rejections of *;
if wmax < .001 then wtypeI(|1,3|)= wtypeI(|1,3|) + 1;** the Null Hyp. *;
******************************************************************;

 end;
typeI=typeI/rep;
wtypeI=wtypeI/rep;
******************************;
do ii=1 to 3;
FAMSIZE(|1,ii|)=(log((1-typei(|1,ii|))))/(log((1-alp(|1,ii|))));
end;
do ii=1 to 3;
ADJALP(|1,ii|)=1-((1-alp(|1,ii|))##(1/FAMSIZE(|1,ii|)));
end;
results=ALP//typeI//ADJALP//FAMSIZE;
***************************************************            ;
print 'Simulation adjustments for Type I errors';
print 'ANOVA assuming homogeneous variances across groups';
print 'number of variables =' k;
print results(|rowname=rlab format=8.6|);
***************************************************            ;
do ii=1 to 3;
WFAMSIZE(|1,ii|)=(log((1-wtypei(|1,ii|))))/(log((1-alp(|1,ii|))));
end;
do ii=1 to 3;
WADJALP(|1,ii|)=1-((1-alp(|1,ii|))##(1/wFAMSIZE(|1,ii|)));
end;
wresults=ALP//wtypeI//wADJALP//WFAMSIZE;
***************************************************            ;
print 'Simulation adjustments for Type I errors';
print 'ANOVA NOT assuming homogeneous variances across groups';
print 'number of variables =' k;
print wresults(|rowname=rlab format=8.6|);
***************************************************            ;