/*Repeated Measures ANOVA. Power (alternative hypothesis) = {1   1.48   1.96} q=3 n=20 distribution = skewness 0.4 and kurtosis 0.8  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);
%let n = 20;
%let q = 3;
%let coe ={1 1.48 1.96}; /*Means*/


proc iml;
	
	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=0.933582234;
c=0.059238353;
d=0.020543667;


/*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 &q;
   		 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;

/*Alternative hypothesis*/
m1=t(a1)*&coe;
	y0 = y0 + m1;
	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 =   /  ss3 effectsize alpha=0.1 ;
	repeated time  3  / printe printm summary ;
	ods output  GLM.Repeated.WithinSubject.ModelANOVA = table;
	ods output  GLM.ANOVA.y1.OverallEffectSize = table1;
    ods output  GLM.ANOVA.y2.OverallEffectSize = table2;
    ods output  GLM.ANOVA.y3.OverallEffectSize = table3;
	ods output  GLM.Repeated.WithinSubject.Epsilons = table4;
    ods output  GLM.Repeated.WithinSubject.ModelANOVA = table5;
run;
ods listing;
proc iml;
use table;
read all var {"DF" "FValue" "ProbF"} into new; 
p1=new[1.1]; 
p2=new[4.1];
p3=new[2.1];
p4=new[3.1];
use table1;
read all var {"Eta2"} into new1;
p5 = new1[4.1];
use table2;
read all var {"Eta2"} into new2;
p6 = new2[4.1];
use table3;
read all var {"Eta2"} into new3;  
p7 = new3[4.1];
use table1;
read all var {"Omega2"} into new4;
p8 =new4[5.1];
use table2;
read all var {"Omega2"} into new5;
p9 =new5[5.1];
use table3;
read all var {"Omega2"} into new6;
p10 =new6[5.1];
use table4;
read all var {"nValue1"} into new7;
p11=new7[1.1];
p12=new7[2.1];
use table5;
read all var {"MS"} into new8;
p13=new8[1.1];
p14=new8[2.1];
use table5;
read all var {"SS"} into new9;
p15=new9[1.1];
use table;
read all var {"ProbFGG"} into new10;
p16=new10[1.1];
use table;
read all var {"ProbFHF"} into new11;
p17=new11[1.1];

P = p1||p2||p3||p4||p5||p6||p7||p8||p9||p10||p11||p12||p13||p14||p15||p16||p17;

varnames='P1'||'P2'||'P3'||'P4'||'P5'||'P6'||'P7'||'P8'||'P9'||'P10'
             ||'P11'||'P12'||'P13'||'P14'||'P15'||'P16'||'P17';

create new12 from P[colname=varnames];
append from P;
proc append base=firstdata;
run;
%end;

/*Power estimation*/

data power;
set firstdata;
proc iml;
use power;
read all var('P1':'P17')into power;
options dkricond=nowarning;

data power;
set power;
     alpha = 0.05;
     corry1y2=P8;
     corry1y3=P9;
     corry2y3=P10;
     corr_mean=mean(corry1y2, corry1y3, corry2y3);
     epsilon_GG=P11; 
	 epsilon_HF=P12;
	 CMM= p13;
	 CME=p14;
	 SCM=p15;
     numdf=P1; 
     dendf=P2; 
     fvalue=P3; 
     probf=P4;
	 pF=probf;
	 pGG = p16;;
	 pHF = p17;;
     Eta2=P5 + P6 + P7;
	 Omega2 = (numdf*(CMM-CME))/((SCM)+(&n*&q-numdf)*CME);

	 *effect size f;
     ncp_F = numdf*fvalue;
     fcrit = finv(1-alpha, numdf, dendf, 0);
     power_F = 1 - probf(fcrit, numdf, dendf, ncp_F);
     if fvalue <= fcrit then h1F = 0; else h1F = 1;

	 if pF <= 0.05 then hF = 1; else hF = 0;
     if pGG <= 0.05 then hGG = 1; else hGG = 0;
     if pHF <= 0.05 then hHF = 1; else hHF = 0;

	 *Power Muller & Barton GG;
     numdf_GG = numdf*epsilon_GG;
	 dendf_GG = dendf*epsilon_GG;
     ncp_GG = numdf_GG*fvalue;
     fcrit_GG = finv(1-alpha, numdf_GG, dendf_GG, 0);
     power_GG_MullerBarton = 1 - probf(fcrit_GG, numdf_GG, dendf_GG, ncp_GG);
     if fvalue <= fcrit_GG then h1GG = 0; else h1GG = 1;
  
     *Power Muller & Barton HF;
     numdf_HF = numdf*epsilon_HF;
	 dendf_HF = dendf*epsilon_HF;
     ncp_HF = numdf_HF*fvalue;
     fcrit_HF = finv(1-alpha, numdf_HF, dendf_HF, 0);
     power_HF_MullerBarton = 1 - probf(fcrit_HF, numdf_HF, dendf_HF, ncp_HF);
     if fvalue <= fcrit_HF then h1HF = 0; else h1HF = 1;

proc freq data=power; 
title 'Results time effect';
tables hF hGG hHF;
run;
proc means data=power;
var power_F power_GG_MullerBarton power_HF_MullerBarton;
         
	 
output out=MeanEst;
ods results off;
run;

%if (&clean) %then %do; 
      proc datasets library=work; 
      delete firstdata meanest mydata new12 power table table1 table2 table3 table4 table5;
      run; 
      run; 
%end;
quit;
%mend measurep;
options mprint;
%measurep(simula=10000);