*The macro is designed to perform bootstrap resampling and analysis of variance (ANOVA) on repeated measures data (simulation study on power);

%macro initial(n1,te,nboot,nsimul);
%do f = &n1 %to &n1;
%do g = &te %to &te;
%do h = &nboot %to &nboot;
%do j = &nsimul %to &nsimul;

dm 'clear log';
dm 'odsresults;select all;clear;';

proc iml;

*Number of repeated measures;
mr=3;

*Distribution (Skewness = 1 and Kurtosis = 1.5);
b= 0.953076898;
c= 0.163194276;
d= 0.006597370;*/

*Unstructured Covariance Matrix, epsilon=0.80;
s1={1.2272416 0.0804696 0.4057249,
    0.0804696 1.1129622 0.6468167,
    0.4057249 0.6468167 1.0578321};

r=(nrow(s)*ncol(s1));
ss=t((shape(s1,1,r)));

*Alternative Hypothesis H1 effect_size te=1;
tee=&te;
beta=(tee*{1 1.27 1.54}); /*Means*/


*Output;
varnames='h1';create nmr from mr [colname=varnames];append from mr;close nmr;
varnames='h2';create ss_1 from ss [colname=varnames];append from ss;close ss_1;
varnames='h3';create bb from b [colname=varnames];append from b;close bb;
varnames='h4';create cc from c [colname=varnames];append from c;close cc;
varnames='h5';create dd from d [colname=varnames];append from d;close dd;
varnames='h6';create te from tee [colname=varnames];append from tee;close te;

n1=&n1;n=n1;
nx=n1;

do i=1 to ncol(nx);nx1=j(nx[i],1,i);
if i=1 then x0=nx1;
else x0=x0//nx1;
end;

*Variance-covariance matrices with the desired heterogeneity;  
v1=s1;

*Data generation;
start ingre(va,na)global(v1,n1);
va=v1;na=n1;
finish;
call ingre(va,na);

start genera(yj,vj,nj)global(b,c,d,mr);
	am=1*-c;
	nvar=nrow(vj);
	var=diag(vj);
	rho=inv(sqrt(var))*vj*inv(sqrt(var));
	rhovec=colvec(rho);
	corcnt=nrow(rhovec); 
	corpass=1;
	tcor=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);
	tcor=tcor//roots[1];
	corpass=corpass+1;
	end;
	tcor=tcor[2:corcnt];
	trho=shape(tcor,nvar,nvar);
	zj=j(nj,mr,.);
	do i=1 to nj;
      do j=1 to mr;
         zj[i,j]=rannor(31418);
         end;
      end;
	tx=zj*root(trho); 
	c0=am + b*tx + c*tx##2 + d*tx##3; 
	yj=c0*sqrt(var);
finish;
call genera(ya,va,na);
x1=repeat(beta[1,1:mr],n1,1);
ya1=ya+x1;
yy=(ya1);

*DATA MATRIX OUTPUT;
y=t(insert(t(yy),t(x0),1,0));
varnames='y0':'y3';
create Bootstraps from y [colname=varnames];
append from y;
close Bootstraps;

data Bootstrap (rename=(y0=group));set Bootstraps;run;

*RUN GLM PROCEDURE;
ods listing close; 
proc glm data=bootstrap plots=none;
 	model y1-y3 = / nouni;
	repeated time  3 / printe;
	ods output  ModelANOVA= probf_orig;
	ods output Epsilons = adjust;
	ods output  ModelANOVA= df;
run;
ods listing;

*BOOTSTRAP METHOD;

*Step 1. Center the data, subtracting the respective mean of each level of the repeated measure from each observation;
proc sql; 
 create table bootstrap_c as select *, 
 mean(y1) as meany1, y1-mean(y1) as c_y1,
 mean(y2) as meany2, y2-mean(y2) as c_y2,
 mean(y3) as meany3, y3-mean(y3) as c_y3
 from bootstrap 
 group by group; 
 quit;

*Step 2. Generate B bootstrap samples with replacements by random sampling;
ods listing close;
sasfile Bootstrap_c load;
proc surveyselect data=Bootstrap_c out=Bootstrap_b seed=31418 method=urs n=(&n1) outhits
rep=&nboot;
strata group;
ods output summary=n_rep;
run;
sasfile Bootstrap_c close;
proc sort data=Bootstrap_b;
by replicate;
run;
ods listing;

*Step 3. Fit the General Lineal Model (GLM) to the bootstrap data and repeat steps 2 and 3 B times;
ods listing close; 
proc glm data=bootstrap_b plots=none;
	by replicate; 
 	model c_y1-c_y3 = / nouni;
	repeated time  3 / printe;
	ods output  ModelANOVA= probf_boot;
run;
ods listing;

*ESTIMATION OF THE SIGNIFICANCE LEVEL USING THE BOOTSTRAP METHOD;
proc iml;
t101=0;

use nmr; read all variables ('h1')into nmr;n=&n1;m=n*nmr;
use probf_orig;;read all variables {'ProbF' 'ProbFGG' 'ProbFHF' 'FValue'}into probf_orig;
use adjust;read all var {'nValue1'} into adjust;
use df;read all var {'df'} into df;
q1=probf_orig[1,1];
q2=probf_orig[1,2];
q3=probf_orig[1,3];
fa=probf_orig[1,4];
gg=adjust[1,1];
hf=adjust[2,1];
df1=df[1,1];
df2=df[2,1];

EtaSq_Orig=(df1*(Fa))/(df1*(Fa)+df2);
EtaSq_GG=((gg*df1)*(Fa))/((gg*df1)*(Fa)+(gg*df2));
EtaSq_HF=((hf*df1)*(Fa))/((hf*df1)*(Fa)+ (hf*df2));
effect_size_f_Orig=sqrt((df2*EtaSq_Orig)/(m*(1-EtaSq_Orig)));
effect_size_f_GG=sqrt((gg*df2*EtaSq_GG)/(m*(1-EtaSq_GG)));
effect_size_f_HF=sqrt((hf*df2*EtaSq_HF)/(m*(1-EtaSq_HF)));

ncp_orig=df1*(Fa);
fcrit_orig=finv(1-.05,df1,df2,0);
power_orig=1-probf(fcrit_orig,df1,df2,ncp_orig);
ncp_gg=(gg*df1)*(Fa);
fcrit_gg=finv(1-.05,gg*df1,gg*df2,0);
power_gg=1-probf(fcrit_gg,gg*df1,gg*df2,ncp_gg);
ncp_hf=(hf*df1)*(Fa);
fcrit_hf=finv(1-.05,hf*df1,hf*df2,0);
power_hf=1-probf(fcrit_hf,hf*df1,hf*df2,ncp_hf);

*B-F values obtained by applying the Bootstrap method;
use Probf_boot;read all variables {'FValue'}into Probf_boot;
nb=2;nk=(nrow(probf_boot))/nb;
a=(shape((Probf_boot`),nk));
cm=a;
f=1;
m=0;
do j=1 to ncol(nb);
	l=m+nb[j];
end;
	do k=1 to nrow(cm);
		ct=cm[k,f:l];
		boot=ct;	
fba=boot[1];

*The proportion of B-F values higher than the observed F-value represents the Bootstrap p-value;
if fba>=fa then do;t101=t101+1;end;q4=t101/&nboot;
end;

f=q1||q2||q3||q4||gg||hf||effect_size_f_Orig||effect_size_f_GG||effect_size_f_HF||power_orig||power_gg||power_hf;


varnames='f1':'f12';create BF from f[colname=varnames];append from f;close BF;

data boot;set BF;
nsimula1=1;
run;

*Macro loop;

%macro loop(simula);
%do i = 2 %to &simula;

dm 'clear log';
dm 'odsresults;select all;clear;';

proc iml;
use nmr; read all variables ('h1')into nmr;
use ss_1;read all variables ('h2')into ss_1;
use bb;read all variables ('h3')into bb;
use cc;read all variables ('h4')into cc;
use dd;read all variables ('h5')into dd;
use te;read all variables ('h6') into te;

s1=shapecol(ss_1,nmr);

n1=&n1;n=n1;
nx=n1;

do i=1 to ncol(nx);nx1=j(nx[i],1,i);
if i=1 then x0=nx1;
else x0=x0//nx1;
end;

*Variance-covariance matrices with the desired heterogeneity;   
v1=s1;

*Alternative Hypothesis H1 effect_size te=1;
tee=&te;
beta=(tee*{1 1.27 1.54}); /*Means*/

*Data generation;
start ingre(va,na)global(v1,n1);
va=v1;
na=n1;
finish;
call ingre(va,na);

start genera(yj,vj,nj)global(bb,cc,dd,nmr);
	am=1*-cc;
	nvar=nrow(vj);
	var=diag(vj);
	rho=inv(sqrt(var))*vj*inv(sqrt(var));
	rhovec=colvec(rho);
	corcnt=nrow(rhovec); 
	corpass=1;
	tcor=0;
	do until(corpass=corcnt+1);
		r=-1*rhovec[corpass,1];
    coeff=(6*dd##2)//(2*cc##2)//(bb##2+6*bb*dd+9*dd##2)//t(r);
	roots=polyroot(coeff);
	nroots=nrow(roots);
	tcor=tcor//roots[1];
	corpass=corpass+1;
	end;
	tcor=tcor[2:corcnt];
	trho=shape(tcor,nvar,nvar);
	zj=j(nj,nmr,.);
	do i=1 to nj;
      do j=1 to nmr;
         zj[i,j]=rannor(31418*(&I-1));
         end;
      end;
	tx=zj*root(trho); 
	c0=am + bb*tx + cc*tx##2 + dd*tx##3; 
	yj=c0*sqrt(var);
finish;
call genera(ya,va,na);
x1=repeat(beta[1,1:nmr],n1,1);
ya1=ya+x1;
yy=(ya1);

*DATA MATRIX OUTPUT;
y=t(insert(t(yy),t(x0),1,0));
varnames='y0':'y3';
create Bootstraps from y [colname=varnames];
append from y;
close Bootstraps;

data Bootstrap (rename=(y0=group));set Bootstraps;run;

*RUN GLM PROCEDURE;
ods listing close; 
proc glm data=bootstrap plots=none;
	model y1-y3 = / nouni;
	repeated time 3 / printe;
	ods output  ModelANOVA= probf_orig;
	ods output Epsilons = adjust;
	ods output  ModelANOVA= df;
run;
ods listing;

*BOOTSTRAP METHOD;

*Step 1. Center the data, subtracting the respective mean of each level of the repeated measure from each observation;
proc sql; 
 create table bootstrap_c as select *, 
 mean(y1) as meany1, y1-mean(y1) as c_y1,
 mean(y2) as meany2, y2-mean(y2) as c_y2,
 mean(y3) as meany3, y3-mean(y3) as c_y3
 from bootstrap 
 group by group; 
 quit;

*Step 2. Generate B bootstrap samples with replacements by random sampling;
ods listing close;
sasfile Bootstrap_c load;
proc surveyselect data=Bootstrap_c out=Bootstrap_b seed=31418 method=urs n=(&n1) outhits
rep=&nboot;
strata group;
ods output summary=n_rep;
run;
sasfile Bootstrap_c close;

proc sort data=Bootstrap_b;
by replicate;
run;
ods listing;

*Step 3. Fit the General Lineal Model (GLM) to the bootstrap data and repeat steps 2 and 3 B times;
ods listing close; 
proc glm data=bootstrap_b plots=none;
	by replicate; 
 	model c_y1-c_y3 = / nouni;
	repeated time  3 / printe;
	ods output  ModelANOVA= probf_boot;
run;
ods listing;

*ESTIMATION OF THE SIGNIFICANCE LEVEL USING THE BOOTSTRAP METHOD;
proc iml;
t101=0;

use nmr; read all variables ('h1')into nmr;n=&n1;m=n*nmr;
use probf_orig;read all variables {'ProbF' 'ProbFGG' 'ProbFHF' 'FValue'}into probf_orig;
use adjust;read all var {'nValue1'} into adjust;
use df;read all var {'df'} into df;
q1=probf_orig[1,1];
q2=probf_orig[1,2];
q3=probf_orig[1,3];
fa=probf_orig[1,4];
gg=adjust[1,1];
hf=adjust[2,1];
df1=df[1,1];
df2=df[2,1];

EtaSq_Orig=(df1*(Fa))/(df1*(Fa)+df2);
EtaSq_GG=((gg*df1)*(Fa))/((gg*df1)*(Fa)+(gg*df2));
EtaSq_HF=((hf*df1)*(Fa))/((hf*df1)*(Fa)+ (hf*df2));
effect_size_f_Orig=sqrt((df2*EtaSq_Orig)/(m*(1-EtaSq_Orig)));
effect_size_f_GG=sqrt((gg*df2*EtaSq_GG)/(m*(1-EtaSq_GG)));
effect_size_f_HF=sqrt((hf*df2*EtaSq_HF)/(m*(1-EtaSq_HF)));

ncp_orig=df1*(Fa);
fcrit_orig=finv(1-.05,df1,df2,0);
power_orig=1-probf(fcrit_orig,df1,df2,ncp_orig);

ncp_gg=(gg*df1)*(Fa);
fcrit_gg=finv(1-.05,gg*df1,gg*df2,0);
power_gg=1-probf(fcrit_gg,gg*df1,gg*df2,ncp_gg);

ncp_hf=(hf*df1)*(Fa);
fcrit_hf=finv(1-.05,hf*df1,hf*df2,0);
power_hf=1-probf(fcrit_hf,hf*df1,hf*df2,ncp_hf);

*B-F values obtained by applying the Bootstrap method;
use Probf_boot;read all variables {'FValue'}into Probf_boot;
nb=2;nk=(nrow(probf_boot))/nb;
a=(shape((Probf_boot`),nk));
cm=a;
f=1;
m=0;
do j=1 to ncol(nb);
	l=m+nb[j];
end;
	do k=1 to nrow(cm);
		ct=cm[k,f:l];
		boot=ct;	
fba=boot[1];

*The proportion of B-F values higher than the observed F-value represents the Bootstrap p-value;
if fba>=fa then do;t101=t101+1;end;q4=t101/&nboot;
end;

f=q1||q2||q3||q4||gg||hf||effect_size_f_Orig||effect_size_f_GG||effect_size_f_HF||power_orig||power_gg||power_hf;

varnames='f1':'f12';create BF from f[colname=varnames];append from f;close BF;

PROC IML;
data BF;set BF; 
nsimula = &i; 
run; 
data boot;set boot BF;run;
%end; 
%mend; 
%loop(&nsimul);
proc iml;
data salida;set boot;run;

data salida1(rename=(f1=ProbF1 f4=ProbF1_BOOT));
set salida;run;

*Proportion of true acceptances of H1;
data salida1; set salida1;
if ProbF1 <=0.05 then ProbF = 1;		if ProbF1>0.05 then ProbF = 0;
if ProbF1_BOOT <=0.05 then ProbF_BOOT=1;if ProbF1_BOOT>0.05 then ProbF_BOOT =0;

proc freq data = salida1; 
table ProbF ProbF_BOOT;
ods output onewayfreqs = significance; 
run;

%end;
%end;
%end;
%end;
%mend;
%initial(20,1,599,5000);*n1=20,effect_size=1,nboot=599,nsimul=5000;