*The macro is designed to perform bootstrap resampling
and analysis of variance (ANOVA)on 2x3 split-plot design (simulation study on Type I error rate);

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

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

proc iml;

*Number of repeated measures and groups;
mr=3;gr=2;
heter=1;hetero=5;

*Distribution (Skewness = 0.4 and Kurtosis = 0.8);
b=0.933582234;
c=0.059238353;
d=0.020543667;

*Unstructured Covariance Matrix, epsilon=0.70;
s1={1.2316562 0.8707401 0.3559481,
    0.8707401 1.1373441 0.1723334,
    0.3559481 0.1723334  1.675346};

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

*Null Hypothesis H0 effect_size te=0;
tee=&te;
beta=(tee*{0 0 0, 0 0 0}); 

*Output;
varnames='h1';create nmr from mr [colname=varnames];append from mr;close nmr;
varnames='h2';create ngr from gr [colname=varnames];append from gr;close ngr;
varnames='h3';create ss_1 from ss [colname=varnames];append from ss;close ss_1;
varnames='h4';create heter_h from heter [colname=varnames];append from heter;close heter_h;
varnames='h5';create hetero_h from hetero [colname=varnames];append from hetero;close hetero_h;
varnames='h6';create bb from b [colname=varnames];append from b;close bb;
varnames='h7';create cc from c [colname=varnames];append from c;close cc;
varnames='h8';create dd from d [colname=varnames];append from d;close dd;
varnames='h9';create te from tee [colname=varnames];append from tee;close te;

n1=&n1;n2=&n2;n=n1+n2;
nx=n1||n2;

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=heter*(s1);v2=hetero*(s1);

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

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);call genera(yb,vb,nb);
x1=repeat(beta[1,1:mr],n1,1);
x2=repeat(beta[2,1:mr],n2,1);
ya1=ya+x1;yb1=yb+x2;
yy=(ya1//yb1);

*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;
	class group;
 	model y1-y3 = group/ nouni;
	repeated time 3/ printe;
	ods output  ModelANOVA= probf_orig;
	ods output Epsilons = adjust;
run;
ods listing;

*BOOTSTRAP METHOD;

*Step 1. Center the data, subtracting the respective mean of each level of the repeated measure from each observation.
This centering operation is performed for each group;

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 Time each group;
sasfile Bootstrap_c load;
proc surveyselect data=Bootstrap_c out=Bootstrap_b seed=31418 method=urs n=(&n1 &n2) outhits
rep=&nboot;
strata group;
ods output summary=n_rep;
run;
sasfile Bootstrap_c close;

proc sort data=Bootstrap_b;
by replicate;
run;

*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; 
	class group;
	model c_y1-c_y3 = group/ 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;t201=0;t301=0;
use probf_orig;read all variables {'ProbF' 'ProbFGG' 'ProbFHF' 'FValue'}into probf_orig;
use adjust;read all var {'nValue1'} into adjust;
q1=probf_orig[1,1]; *pvalue-Factor A;
q2=probf_orig[3,1]; *pvalue-Factor B;
q3=probf_orig[3,2]; *pvalueGG-Factor B;
q4=probf_orig[3,3]; *pvalueHF-Factor B;
q5=probf_orig[4,1]; *pvalue-interaction AB;
q6=probf_orig[4,2]; *pvalueGG-interaction AB;
q7=probf_orig[4,3]; *pvalueHF-interaction AB;
fa=probf_orig[1,4]; *Fvalue-Factor A;
fb=probf_orig[3,4]; *Fvalue-Factor B;
fi=probf_orig[4,4]; *Fvalue-Interaction AB;
gg=adjust[1,1];
hf=adjust[2,1];

*B-F values obtained by applying the Bootstrap method;
use Probf_boot;read all variables {'FValue'}into Probf_boot;
nb=nrow(probf_orig);
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];fbb=boot[3];fbi=boot[4];

*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;q8=t101/&nboot;
if fbb>=fb then do;t201=t201+1;end;q9=t201/&nboot;
if fbi>=fi then do;t301=t301+1;end;q10=t301/&nboot; 
end;

h=q1||q2||q3||q4||q5||q6||q7||q8||q9||q10||gg||hf;

varnames='h1':'h12';create BF from h[colname=varnames];append from h;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 ngr; read all variables ('h2')into ngr;
use ss_1;read all variables ('h3')into ss_1;
use heter_h;read all variables ('h4')into heter_h;
use hetero_h;read all variables ('h5')into hetero_h;
use bb;read all variables ('h6')into bb;
use cc;read all variables ('h7')into cc;
use dd;read all variables ('h8')into dd;
use te;read all variables ('h9') into te;

s1=shapecol(ss_1,nmr);

heter=heter_h;hetero=hetero_h;
n1=&n1;n2=&n2;n=n1+n2;
nx=n1||n2;

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=heter*(s1);v2=hetero*(s1);

*Null Hypothesis H0 effect_size te=0;
beta=(TE*{0 0 0,0 0 0}); 

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

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);
         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);call genera(yb,vb,nb);
x1=repeat(beta[1,1:nmr],n1,1);
x2=repeat(beta[2,1:nmr],n2,1);
ya1=ya+x1;yb1=yb+x2;
yy=(ya1//yb1);

*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;
	class group;
 	model y1-y3 = group/ nouni;
	repeated time 3/ printe;
	ods output  ModelANOVA= probf_orig;
	ods output Epsilons = adjust;
run;
ods listing;

*BOOTSTRAP METHOD;

*Step 1. Center the data, subtracting the respective mean of each level of the repeated measure from each observation.
This centering operation is performed for each group;
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 Time each group;
sasfile Bootstrap_c load;
proc surveyselect data=Bootstrap_c out=Bootstrap_b seed=31418 method=urs n=(&n1 &n2) outhits
rep=&nboot;
strata group;
ods output summary=n_rep;
run;
sasfile Bootstrap_c close;

proc sort data=Bootstrap_b;
by replicate;
run;

*Step 3. Fit the General Linel 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; 
	class group;
	model c_y1-c_y3 = group/ 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;t201=0;t301=0;

use probf_orig;read all variables {'ProbF' 'ProbFGG' 'ProbFHF' 'FValue'}into probf_orig;
use adjust;read all var {'nValue1'} into adjust;
q1=probf_orig[1,1]; *pvalue-Factor A;
q2=probf_orig[3,1]; *pvalue-Factor B;
q3=probf_orig[3,2]; *pvalueGG-Factor B;
q4=probf_orig[3,3]; *pvalueHF-Factor B;
q5=probf_orig[4,1]; *pvalue-interaction AB;
q6=probf_orig[4,2]; *pvalueGG-interaction AB;
q7=probf_orig[4,3]; *pvalueHF-interaction AB;
fa=probf_orig[1,4]; *Fvalue-Factor A;
fb=probf_orig[3,4]; *Fvalue-Factor B;
fi=probf_orig[4,4]; *Fvalue-Interaction AB;
gg=adjust[1,1];
hf=adjust[2,1];

*B-F values obtained by applying the Bootstrap method;
use Probf_boot;read all variables {'FValue'}into Probf_boot;
nb=nrow(probf_orig);
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];fbb=boot[3];fbi=boot[4];

*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;q8=t101/&nboot;
if fbb>=fb then do;t201=t201+1;end;q9=t201/&nboot;
if fbi>=fi then do;t301=t301+1;end;q10=t301/&nboot; 
end;
h=q1||q2||q3||q4||q5||q6||q7||q8||q9||q10||gg||hf;

varnames='h1':'h12';create BF from h[colname=varnames];append from h;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=(
h1=ProbF1_Group 
h2=ProbF1_Time 
h3=ProbF1_Time_GG 
h4=ProbF1_Time_HF 
h5=ProbF1_Interaction 
h6=ProbF1_Interaction_GG 
h7=ProbF1_Interaction_HF 
h8=ProbF1_Group_BOOT 
h9=ProbF1_Time_BOOT 
h10=ProbF1_Interaction_BOOT 
h11=GG_Epsilon 
h12=HF_Epsilon));
set salida;run;

*Proportion of false rejections of H0;
data salida1; set salida1;
if ProbF1_Group<=0.05 then ProbF_Group=1;			if ProbF1_Group>0.05 then ProbF_Group=0;
if ProbF1_Time<=0.05 then ProbF_Time=1;			    if ProbF1_Time>0.05 then ProbF_Time=0;
if ProbF1_Time_GG<=0.05 then ProbF_Time_GG=1;		if ProbF1_Time_GG>0.05 then ProbF_Time_GG=0;
if ProbF1_Time_HF<=0.05 then ProbF_Time_HF=1;		if ProbF1_Time_HF>0.05 then ProbF_Time_HF=0;
if ProbF1_Interaction<=0.05 then ProbF_Interaction=1;			if ProbF1_Interaction>0.05 then ProbF_Interaction=0;
if ProbF1_Interaction_GG<=0.05 then ProbF_Interaction_GG=1;		if ProbF1_Interaction_GG>0.05 then ProbF_Interaction_GG= 0;
if ProbF1_Interaction_HF<=0.05 then ProbF_Interaction_HF=1;		if ProbF1_Interaction_HF>0.05 then ProbF_Interaction_HF= 0;
if ProbF1_Group_BOOT<=0.05 then ProbF_Group_BOOT=1;	        if ProbF1_Group_BOOT>0.05 then ProbF_Group_BOOT=0;
if ProbF1_Time_BOOT<=0.05 then ProbF_Time_BOOT=1;	        if ProbF1_Time_BOOT>0.05 then ProbF_Time_BOOT=0;
if ProbF1_Interaction_BOOT<=0.05 then ProbF_Interaction_BOOT=1;	if ProbF1_Interaction_BOOT>0.05 then ProbF_Interaction_BOOT=0;

proc freq data = salida1; 
table ProbF_Group ProbF_Time ProbF_Time_GG ProbF_Time_HF ProbF_Interaction ProbF_Interaction_GG ProbF_Interaction_HF ProbF_Group_BOOT ProbF_Time_BOOT ProbF_Interaction_BOOT;
ods output onewayfreqs = significance; 
run;


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