/*Repeated Measures ANOVA. Type I error. q=3 n=20 distribution = skewness 0.8 and kurtosis 0 epsilon = .90*/ 

%macro measurep(SIMULA=,options=NOPRINT);

%let iter  = 1;
%let clean = 1; 
%do %while (&iter<=&simula);
%put the value of iter is &iter;
%let iter = %eval(&iter + 1); 
 
proc iml;
	n=20;
	q=3;vd=1;k=vd*q;
	r = I(n);
	a1=j(1,n,1);

/*Unstructured covariance matrix UN epsilon = .90*/

varcov =     {5.0244188   0.75767 0.9551522,
              0.75767 1.2701189 0.5017649,
              0.9551522 0.5017649 3.3719781};
v =varcov; 
 
/*Fleishman's constants for non-normal distribution*/

b= 1.136856585;
c= 0.202117699;
d=-0.063802837;

*Vale and Maurelli (1983) procedure to obtain multivariate non-normality*;


a=-1*c;

nvar=nrow(v);
var=diag(v);
rho=inv(sqrt(var))*v*inv(sqrt(var));
rhovec=colvec(rho);
corcnt=nrow(rhovec);
corpass=1;
cor=0;
do until(corpass=corcnt+1);
	r=-1*rhovec[corpass,1];
coeff=(6*d**2)//(2*c**2)//(b**2+6*b*d+9*d**2)//t(r);
roots=polyroot(coeff);
nroots=nrow(roots);
cor=cor//roots[1];
corpass=corpass+1;
end;
cor=cor[2:corcnt];
rho=shape(cor,nvar,nvar);
	z=j(n,k,.);
	do i=1 to n;
  		do j=1 to k;
   		 z[i,j]=rannor(0);
  		end;
	end;
x1=z*root(rho);
y1=a + b*x1 + c*x1##2 + d*x1##3;
x1=y1*sqrt(var);
    y0= x1;
	w0=1:n;
	datos=t(w0)||y0;

 /*Data matrix and GLM analysis*/
	y=datos;
	varnames='y0':'y3';
	create mydata from y[colname=varnames];
	append from datos;
	data mydata;
	set mydata;
		subject=y0;
		    
run;
ods listing close; 
proc glm data=mydata plots=none;
 	model y1-y3 =   / nouni   ;
	repeated time  3 / printe  ;
	ods output  GLM.Repeated.WithinSubject.ModelANOVA= table;
	ods output  GLM.Repeated.WithinSubject.Epsilons = table1;
run;
ods listing;
proc iml;
use table;
read all var {"FValue" "ProbF" "ProbFGG" "ProbFHF"} into new1;
p1=new1[1,1];
p2=new1[1,2];
p3=new1[1,3];
p4=new1[1,4];
use table1;
read all var {"nValue1"} into new2;
p5=new2[1,1];
p6=new2[2,1]; 
P = p1||p2||p3||p4||p5||p6;
 
varnames='P1'||'P2'||'P3'||'P4'||'P5'||'P6';
create new3 from P[colname=varnames];
append from P;
close new3;

proc append base=firstdata;
run;

%end;

data firstdata;
set firstdata;
proc iml;
use firstdata;
read all var('P1':'P6')into firstdata;
data firstdata1(rename=(P1=F_Time P2=ProbF P3=ProbF_AdjGG P4=ProbF_AdjHF P5=Greenhouse_Geisser_Epsilon P6=Huynh_Feldt_Epsilon));
set firstdata; run;
proc means data=firstdata1;
var F_Time ProbF  ProbF_AdjGG ProbF_AdjHF Greenhouse_Geisser_Epsilon Huynh_Feldt_Epsilon;
title 'F_time and Adjusted F';
run;
data firstdata1; 
 
set firstdata1;
if ProbF <=0.05  then H_Time = 1; else H_Time = 0; 
if ProbF_AdjGG <=0.05 then H_ProbF_AdjGG = 1; else H_ProbF_AdjGG = 0; 
if ProbF_AdjHF <=0.05 then H_ProbF_AdjHF = 1;  else H_ProbF_AdjHF = 0;
proc freq data=firstdata1; 
title 'Results time effect';
tables  H_Time  H_ProbF_AdjGG  H_ProbF_AdjHF ;
run;

%if (&clean) %then %do; 
      proc datasets library=work; 
      delete firstdata firstdata1 mydata mydata1 new3 table table1;
      run; 
%end;
quit;
%mend measurep;
options mprint;
%measurep(simula=10000);