*Repeated Measures ANOVA. Type I error. q=3 n=20 epsilon=0.70*/ %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*/ varcov={1.2316562 0.8707401 0.3559481, 0.8707401 1.1373441 0.1723334, 0.3559481 0.1723334 1.675346}; cpo=orpol(1:3,3); cp=cpo[1:3,2:3]; error=cp`*varcov*cp; e=eigval(error); /*Box epsilon*/ epsilon=(trace(error)**2)/((q-1)*(trace((error)**2))); simulacio=iter||epsilon; v =varcov; /*Cholesky decomposition*/ p1=root(v); /*Generation of normal data*/ z =j(n,q,.); do i=1 to n; do j=1 to q; z[i,j]=rand('NORMAL'); end; end; y0=z*p1; 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 uepsdef=HF uepsdef=HFL uepsdef =CM ; ods output GLM.Repeated.WithinSubject.ModelANOVA= table; ods output GLM.Repeated.WithinSubject.Epsilons = table1; ods output GLM.Repeated.MANOVA.Model.Error.Sphericity = table2; 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]; use table2; read all var {"Mauchly"} into new3; p7=new3[2,1]; P = p1||p2||p3||p4||p5||p6||p7; varnames='P1'||'P2'||'P3'||'P4'||'P5'||'P6'||'P7'; create new from P[colname=varnames]; append from P; close new; proc datasets; append base=firstdata; run; %end; data firstdata; set firstdata; proc iml; use firstdata; read all var('P1':'P7')into firstdata; data firstdata1(rename=(P1=F_Time P2=ProbF P3=ProbF_AdjGG P4=ProbF_AdjHF P5=Greenhouse_Geisser_Epsilon P6=Huynh_Feldt_Epsilon P7=Mauchly_epsilon)); set firstdata; run; proc means data=firstdata1; var F_Time ProbF ProbF_AdjGG ProbF_AdjHF Greenhouse_Geisser_Epsilon Huynh_Feldt_Epsilon Mauchly_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); /*Repeated Measures ANOVA. Power (alternative hypothesis) = {1.5 1.5 2} q=3 n=20 epsilon=0.70*/ %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.5 1.5 2}; proc iml; vd=1;k=vd*&q; r = I(&n); a1=j(1,&n,1); /*Unstructured covariance matrix UN*/ varcov={1.2316562 0.8707401 0.3559481, 0.8707401 1.1373441 0.1723334, 0.3559481 0.1723334 1.675346}; v =varcov; /*Cholesky decomposition*/ p1=root(v); /*Generation of normal data*/ z =j(&n,&q,.); do i=1 to &n; do j=1 to &q; z[i,j]=rand('NORMAL'); end; end; /*Alternative hypothesis*/; m1=t(a1)*&coe; y0= z*p1; 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 = sqrt(Omega2/(1-Omega2)); *Effect Size f; ncp = numdf*fvalue; fcrit = finv(1-alpha, numdf, dendf, 0); power = 1 - probf(fcrit, numdf, dendf, ncp); if fvalue <= fcrit then h1 = 0; else h1 = 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; proc freq data=power; title 'Results time effect'; tables h1 hF hGG hHF; run; proc means data=power; var power ncp Omega2 effect_size_f epsilon_GG epsilon_HF; 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);