%let path=/folders/myfolders/Memoire;
libname sor "&path./Output" access=readonly ;
libname don "&path./data/mlm" ;

/*Simulation dans le cas de la régression linéaire classique*/
%let N = 30;
%let NumSamples = 1000;
data SimG (replace=yes); 
call streaminit(123); 
do SampleID = 1 to &NumSamples;  /* Création de NumSamples échantillons */
   do iG = 1 to &N;               /* Echantillons de taille N */
      XiG = rand("Uniform", 0, 10); /* XiG ~ U(0,10) */
   	  EpG = rand("NORMal", 0, 1); 
      output;
   end;
end;
run;

/*Moyenne des XiG et EpG par échantillon
proc means data=SimG mean std ;
	by sampleId; 
	var xig epg; 
	output out=SimStat mean=MeanXiG MeanEpG
			Std=StdXiG StdEpG ;
run;

Moyenne des XiG et EpG sur l'ensemble des 1000 échantillons
proc means data=SimStat mean Std;
 var MeanXig meanepg stdxig stdepg;
run;*/
	
options mcompilenote=all; 
%macro mod (EP, Xi, bet, m);
/*Création des valeurs atypiques et des vecteurs di les caractérisant*/
data Sim&EP (replace=yes);
	set simG ;
	d1=0;
	d12=0;
	d30=0;
	if iG=1 then d1=1; /*L'observation 1 n'est pas atypique. Elle servira pour la comparaison*/
	if iG=12 then 
		do EpG= &ep;
			d12=1;
		end;
	if iG=30 then 
		do XiG=Ξ 
			EpG=&ep;
			d30=1;
		end;
	mu=&m; /*Intercept du mdèle*/
	beta=&bet; /*Pente*/
	yG=(beta*XiG)+mu+EpG; /* Création du vecteur réponse y */ 
run;

/*Ajustement du modèle linéaire puis stockage des résidus standardisés puis studentisés*/
proc reg data=sim&ep noprint;
	by SampleID; /*Un modèle par échantillon*/
	model yG = XiG;
	output out=RegStat rstudent=RES student=Rstand;
run;

/*Création d'une variable indicatrice si l'observation est identifiée comme aberrante*/ 
Data mem.regstat&ep;
	set regstat;
	if abs(res) lt 1.96 then abb="Non-aberrante" ; 
	else abb="Aberrante";						
	ti = rstand*rstand;							/*Calcul des RIS au carré*/
	if ti gt 1 then lrti= (27 * log(27/(28-ti))) - log (ti);	/*Calcul de la stat LRT*/
	else lrti=0; 
	if ti lt 3.77 then bti = "Non-aberrante";
		else bti = "Aberrante"  ;
	if lrti lt 1.60 then blrti= "Non-aberrante";
		else blrti = "Aberrante" ;
run;

/*VSOM pour les observations 1, 12 et 30 */
%ODSOff;
PROC MIXED data = Sim&ep;
	TITLE "VSOM1&12&30"; 
	by SampleID;
	CLASS iG ; 
	MODEL yG=XiG/ SOLUTION ; 
	RANDOM d1 d12 d30/ TYPE=VC SUBJECT=iG V G; 
	ods output CovParms=COV SolutionF=SF;
RUN;
%ODSOn;

/*Calcul de wi*/
Proc sql;
	create table mem.wtab&ep as
	Select *, 
	estimate/(select estimate 
				from cov as b 
				where CovParm="Residual" and a.sampleId=b.sampleId) as wi
	from cov as a;
quit;

*Renommer les coefficients associés aux effets du modèle;
data mem.sf&ep;
	set sf;
	if effect="XiG" then effect="Pente";
	else effect=effect;
run;

*Distribution de la variance des résidus;
proc univariate data=mem.wtab&ep;
	Title1 "VSOM simple, Epsilon=&ep";
	Title2 "Distribution de la variance des résidus" ;
	where Covparm="Residual" ;
	var Estimate; 
run;
title;

PROC SGPLOT DATA=mem.wtab&ep NOBORDER;
	where Covparm="Residual" ;
	histogram Estimate;
	YAXIS label="Pourcentage d'échantillons" ;
	XAXIS label="Variance des résidus" ;
RUN; 
%mend;

%mod(EP=5, Xi=20, bet=0.25, m=26.25);

*Changement des coefficients du modèles;
*%mod(EP=5, Xi=20, bet=0.5, m=26.25);
*%mod(EP=5, Xi=20, bet=1, m=26.25);
*%mod(EP=5, Xi=20, bet=0.25, m=40);
*%mod(EP=5, Xi=20, bet=1, m=40);

*Changement des valeurs aberrantes pour Xi;
*%mod(EP=5, Xi=10, bet=0.25, m=26.25);
*%mod(EP=5, Xi=1, bet=0.25, m=26.25);
*%mod(EP=5, Xi=0, bet=0.25, m=26.25);

/*Distribution des paramètres du modèle*/
%macro Eff(ep, EF);
ods graphics / noborder;
proc univariate data=mem.sf&ep ;
	Title "Distribution de &EF, Epsilon=&ep" ;
	where Effect="&EF";
	var Estimate;
run;
title;

*Histogramme des paramètres du modèle; 
PROC SGPLOT DATA=mem.sf&ep NOBORDER;
	where Effect="&EF";
	histogram Estimate;
	YAXIS label="Pourcentage d'échantillons" ;
	XAXIS label="&ef";
RUN;
%mend Eff;
%Eff(ep=5, EF=Intercept);
%Eff(ep=5, EF=Pente);

options mcompilenote=all;
%macro res (ep, id);
/*Résumé de la distribution des résidus stockés*/
proc univariate data=mem.regstat&ep cibasic MU0=-1.96 LOCCOUNT;
	Title1 "Regression simple, Epsilon=&ep";
	Title2 "Distribution des résidus de l'observation &Id" ;
	var RES Rstand ti lrti;
   	where iG=&id;
   	histogram ti lrti;
run;
TITLE;

*Histogramme des RES; 
ods graphics / noborder;
PROC SGPLOT DATA =mem.regstat&ep NOBORDER;
	where iG=&id;
	histogram RES /*Rstand*/ ;
	YAXIS label="Pourcentage d'échantillons" ;
	XAXIS label="RES" VALUES= (-2 TO 7);
RUN;

*Histogramme des ti^2; 
ods graphics / noborder;
PROC SGPLOT DATA =mem.regstat&ep NOBORDER;
	where iG=&id;
	histogram ti;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 50) ;
	XAXIS label="RIS au carré" VALUES= (0 TO 18);
RUN;

/*Proportion d'échantillon détectant chacune des observations aberrantes*/
/*Selon les RES*/
ods graphics / noborder;
PROC SGPLOT DATA = mem.regstat&ep NOBORDER;
	where iG=&ID;
	VBAR abb / STAT=percent datalabel ;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 1) ;
	XAXIS label="Valeur aberrante si RES supérieur à 1,96" ;	
RUN;

/*Selon les ti^2*/
PROC SGPLOT DATA = mem.regstat&ep NOBORDER;
	where iG=&id;
	VBAR bti / STAT=percent datalabel ;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 1);
	XAXIS label="Valeur aberrante si ti^2 supérieur à 3,77" ;	
RUN;

/*Selon les LRTi*/
ods graphics / noborder;
PROC SGPLOT DATA = mem.regstat&ep NOBORDER;
	where iG=&id;
	VBAR blrti / STAT=percent datalabel ;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 1);
	XAXIS label="Aberrante si LRTi supérieure à 1,60" ;	
RUN;

/*Analyse des distributions des paramètres de variance et des facteurs d'inflation associés */
ods graphics / noborder;
proc univariate data=mem.wtab&ep;
	Title1 "VSOM simple, Epsilon=&ep" ;
	Title2 "Distribution de la variance de d&Id" ;
	where Covparm="d&Id" ;
	var Estimate wi;
run;
title;

*Histogramme des variances des effets aléatoires di; 
PROC SGPLOT DATA=mem.wtab&ep NOBORDER;
	where Covparm="d&Id" ;
	histogram Estimate;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 80) ;
	XAXIS label="Variance de l'effet aléatoire d&Id" VALUES= (0 TO 70);
RUN;

/*Histogramme des facteurs d'inflation wi*/ 
PROC SGPLOT DATA=mem.wtab&ep NOBORDER;
	where Covparm="d&Id";
	histogram wi;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 80) ;
	XAXIS label="Facteur d'inflation associé à d&ID" VALUES= (0 TO 128) ;
RUN;
%mend;
options mcompilenote=none;
%res (ep=5, Id=1);
%res (ep=5, Id=12);
%res (ep=5, Id=30);

*Changement des valeurs aberrantes pour epsilon;
*%mod(EP=1, Xi=20, bet=0.25, m=26.25);
%res (ep=1, Id=1);
%res (ep=1, Id=12);
%res (ep=1, Id=30);

*%mod(EP=2, Xi=20, bet=0.25, m=26.25);
%res (ep=2, Id=1);
%res (ep=2, Id=12);
%res (ep=2, Id=30);

*%mod(EP=3, Xi=20, bet=0.25, m=26.25);
%res (ep=3, Id=1);
%res (ep=3, Id=12);
%res (ep=3, Id=30);

*%mod(EP=4, Xi=20, bet=0.25, m=26.25);
%res (ep=4, Id=1);
%res (ep=4, Id=12);
%res (ep=4, Id=30);

/*Simulation dans le cas du modèle mixte à un cluster aberrant*/

/*Pour copier les bases entre machines*/
*proc copy inlib=sor outlib=don noclone datecopy ;
*run;

/*Ces macros sont créés pour désactiver l'afficahage des résultats des procédures qui ne supportent pas l'option no print.
La proc mixed par exemple*/
%macro ODSOff(); 
ods graphics off;
ods exclude all;
ods noresults;
%mend;
%macro ODSOn(); 
ods graphics on;
ods exclude none;
ods results;
%mend;

/*Simulation dans le cas d'un modèle mixte à un effet aléatoire*/
/*Création d'une base de données comportant 2 études se déroulant dans 3 pays différents 
	appliquant 2 traitements et recrutant 5 patients par traitement soit 10 patients par 
	étude par pays donc un échantillon de 60 individus.*/
data Simulation;
	y = 0; 
	do country = 1 to 10;
		do treatment = 1 to 2;
			do subject = 1 to 50;
				output;
			end;
		end;
	end;
run;

/*Rajouter une colonne identifiant pour chaque individu*/
Data Simulation;
	Set Simulation;
	Id=_n_;
	If country=1 and treatment=1
		then stuid="A";
	else if country=1 and treatment=2
		then stuid="B";
	If country=2 and treatment=1
		then stuid="C";
	else if country=2 and treatment=2
		then stuid="D";
	If country=3 and treatment=1
		then stuid="E";
	else if country=3 and treatment=2
		then stuid="F";
	If country=4 and treatment=1
		then stuid="G";
	else if country=4 and treatment=2
		then stuid="H";
	If country=5 and treatment=1
		then stuid="I";
	else if country=5 and treatment=2
		then stuid="J";
	If country=6 and treatment=1
		then stuid="K";
	else if country=6 and treatment=2
		then stuid="L";
	If country=7 and treatment=1
		then stuid="M";
	else if country=7 and treatment=2
		then stuid="N";
	If country=8 and treatment=1
		then stuid="O";
	else if country=8 and treatment=2
		then stuid="P";
	If country=9 and treatment=1
		then stuid="Q";
	else if country=9 and treatment=2
		then stuid="R";
	If country=10 and treatment=1
		then stuid="S";
	else if country=10 and treatment=2
		then stuid="T";
run;
	
/*Extraction des matrices X et Z correspondant au design précédemment construit*/
%ODSOff;
Proc glimmix data=Simulation outdesign(names)=Mat;
	class country treatment;
	model y= treatment;
	random country;
	ods select ColumnNames;
	run;
%ODSOn; 

/*Simulation des donnée*/
options mcompilenote=all;
%MACRO DAT(be, alp);	
Proc iml;
	FixedVAR = "_X1":"_X3"; /*Extraction des matrices X et Y*/
	RandomVar = "_Z1":"_Z3";
	use Mat;
	read all var FixedVar into X;
	read all var RandomVar into Z;
	close Mat;
	beta = {&be, &alp, 0}; 			/*Coefficients du modèle : intercept et pente*/
	eta = X*beta;
	varb = 0.1; 						/*Variance intra-groupe*/
	sigma2 = 1; 						/*variance des résidus*/
	G = diag( repeat (varb,3)); 		/*Matrice G de la variance de l'effet aléatoire*/
	call randseed (1); 								/*Seed*/
	zero = repeat (0, nrow(G));
	NumSamples = 100; 								/*Nombres d'échantillons*/
	gamma = RandNormal(NumSamples, zero, G);		/*Remplissage de la marice G*/
	eps = J(nrow(X), NumSamples); 					/*Matrice des residus epsilon*/
	call randgen (eps, "Normal", 0, sqrt(sigma2));	/*Remplissage de la marice des résidus*/
	Y = eta + Z*gamma` + eps;						/*Création de la variale réponse y*/
	yNames = "y1":("y"+strip(char(NumSamples)));
	create Sim from Y[C=yNames];					/*Création du dataset comportant y*/
	append from Y;
	close;
quit;
%mend; 
options mcompilenote=none; 

%dat(be=26.25 , alp=0.25);
*%dat(be=40, alp=1);

Data Sim (replace=yes drop=y);
	Merge Simulation Sim; 
Run;

Proc transpose data=sim 
				out=Lsim (replace=yes rename=(col1=y));
	by Id country stuid treatment subject; 
run;

data lsim;
  set lsim;
  nums=input(substr(_name_, 2), 5.);
  drop _name_;
run; 

Proc sort data=LSim;
	by nums;
run;

%macro grp(pl);
/*Création du groupe abberrant*/
Data don.lSimm&pl (replace=yes); 
	Set lSim;
	If country=1 and treatment=1
		then do; 
			y = y + &pl;
			dA=1;
		end;
	else do y=y;
			dA=0;
		end;
	If country=3 and treatment=1
		then dE=1;
	else dE=0;
Run;

/*VSOM pour les groupes 1 et 6*/ 
%odsoff;
PROC MIXED data = don.lSimm&pl;
	TITLE "VSOM mixte pl=&pl"; 
	by nums;
	CLASS Treatment Country ; 
	MODEL y=Treatment/ SOLUTION ; 
	RANDOM Country dA dE/ TYPE=VC SUBJECT=Id; 
	ods output CovParms=don.COV&pl SolutionF=don.SF&pl;
run;

PROC MIXED data = don.lsimm&pl;
	by nums;
	TITLE "Modèle mixte de base"; 
	CLASS Treatment Country stuid; 
	MODEL y=Treatment/ SOLUTION outp=don.PredBase&pl influence; 
	RANDOM Country/ TYPE=VC SUBJECT=Id; 
	ods output SolutionF=don.SFBaseE&pl CovParms=DON.COVBaseE&pl 
				influence=don.inf&pl; 
RUN;
%odson;

/*Calcul des wi*/
Proc sql;
	create table don.wtab&pl as
	Select *, 
	estimate/(select estimate 
				from don.cov&pl as b 
				where CovParm="Residual" and a.nums=b.nums) as wi
	from don.cov&pl as a;
quit;

*Distribution des coefficients du modèle mixte;
ods graphics / noborder;
proc univariate data=don.SF&pl cibasic ;
	TITLE1 "VSOM mixte pl=&pl"; 
	Title2 "Distribution de l'intercept";
	where Effect="Intercept";
	var Estimate ;
   	histogram Estimate /NOFRAME 
   				VAXISLABEL="Pourcentage"/*d'échantillons*/ ;
   	output out=StatInt mean=MeanInt
   						std=StdInt;
run;
	
proc univariate data=don.SF&pl cibasic ;
	TITLE1 "VSOM mixte pl=&pl"; 
	Title2 "Distribution de la pente";
	where Effect="treatment" and treatment=1;
	var Estimate ;
   	histogram Estimate /NOFRAME 
   				VAXISLABEL="Pourcentage"/*d'échantillons*/ ;
   	output out=Stattreat mean=Meantreat
   						std=Stdtreat;
run;
%grp(pl=10);
%grp(pl=1);

options mcompilenote= all; 
%macro disres (id, pl);
*Résidus externes studentisés;
%odsoff;
PROC MIXED data = don.lsimm&pl;
	by nums;
	WHERE STUID ne "&id" ;
	TITLE "Modèle mixte de base sans &id"; 
	CLASS Treatment Country stuid; 
	MODEL y=Treatment/ SOLUTION /*outp=sor.PredS&id.&pl*/; 
	RANDOM Country/ TYPE=VC SUBJECT=Id; 
	ods output SolutionF=sor.SFBs&id.&pl CovParms=sor.COVBs&id.&pl; 
RUN;
%odson; 

*Calcul de la valeur moyenne prédite pour yj; 
Proc sql;
	create table lsim&ID (drop=intercept pente) as 
	select *, 
			(select estimate from don.sfbs&id.&pl as b where effect="Intercept" and a.nums=b.nums) as intercept, 
			(select estimate from don.sfbs&id.&pl as c where effect="treatment" and treatment=1 and a.nums=c.nums) as pente,
			calculated intercept + calculated pente as pred,
			(select stderr from don.sfbs&id.&pl as d  where effect="Intercept" and a.nums=d.nums) + 
				(select stderr from don.sfbs&id.&pl as e where effect="treatment" and treatment=1 and a.nums=e.nums) as stderrpred
	from don.lsimm&pl as a
	where stuid="&id";
quit; 
	
Proc sql ;
	create table don.Tpreds&id as 
	select *, 
			avg(y) as ymean, 
			avg(pred) as ypred,
			(avg(stderrpred))**2 as varui,
			calculated ymean - calculated ypred as ei,
			(select estimate from don.covbs&id.&pl as b where covparm="Residual" and a.nums=b.nums ) as vi,
			(select estimate from don.covbs&id.&pl as c where covparm="country" and a.nums=c.nums) as taui,
			sqrt(calculated varui + calculated vi + calculated taui) as den,
			calculated ei/calculated den as RES
	from tpres&id as a
	where stuid="&id";
quit;

Data don.Tpreds&id;
	set don.Tpreds&id;
	if abs(res) lt 1.96 then abb="Non-aberrante" ; 
	else abb="Aberrante";	
RUN;

Proc sql;
create table infdi&id as 
	select avg(student) as ti,
			avg(rstudent) as res,
			calculated ti**2 as ti2,
			(select stuid from don.Predbase&pl as b 
							where a.observed=b.y and a.nums=b.nums) as stuid
	from don.inf&pl as a
	group by nums, stuid ; 
quit;

data don.infdi&id.&pl;
	set infdia;
	if ti2 gt 1 
		then lrti = ((1000-2-1) * log((1000-2-1)/(1000-2-ti2))) - log (ti2) ;
		else lrti=0;
	if lrti lt 2.01 then blrti= "Non-aberrante";
		else blrti = "Aberrante" ;
run; 


ods graphics / noborder;
proc univariate data=don.infdi&id.&pl cibasic;
	Title1 "Résidus MIXTE, study=&id";
	where STUID="&id";
	var RES ti2 lrti;
   	histogram RES ti2 lrti/NOFRAME 
   				VAXISLABEL="Pourcentage d'échantillons" ;
run;
title;

*Histogramme des RIS; 
ods graphics / noborder;
PROC SGPLOT DATA =don.infdi&id.&pl NOBORDER;
	where STUID = "&id";
	histogram ti2;
	YAXIS label="Pourcentage d'échantillons" /*VALUES= (0 TO 30)*/;
	XAXIS label="RIS au carré";
RUN;

*Histogramme des LRTi; 
ods graphics / noborder;
PROC SGPLOT DATA =don.infdi&id.&pl NOBORDER;
	where STUID = "&id";
	histogram lrti;
	YAXIS label="Pourcentage d'échantillons" /*VALUES= (0 TO 100)*/ ;
	XAXIS label="Statistique du rapport de vraisemblance LRTj";
RUN;

*Distribution des RES;
ods graphics / noborder;
proc univariate data=don.tpreds&id cibasic;
	Title1 "RES MIXTE, Etude &id";
	where STUID = "&id";
	var RES;
   	histogram RES/NOFRAME 
   				VAXISLABEL="Pourcentage"/*d'échantillons*/ ;
run;
title;

*Histogramme des RES; 
ods graphics / noborder;
PROC SGPLOT DATA =don.tpreds&id NOBORDER;
	where STUID = "&id";
	histogram RES;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 11) ;
	XAXIS label="RES";
RUN;

*Distribution des wi;
ods graphics / noborder;
proc univariate data=DON.wtab&pl cibasic;
	Title1 "Wi MIXTE, Etude &id";
	where CovParm = "d&id";
	var wi;
   	histogram wi/NOFRAME VAXISLABEL="Pourcentage d'échantillons";
run;
title;

*Histogramme des wi; 
ods graphics / noborder;
PROC SGPLOT DATA =DON.wtab&pl NOBORDER;
	where CovParm = "d&id";
	histogram wi;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 60) ;
	XAXIS label="Facteur d'inflation associé à la variance des effets aléatoires dj";
RUN;

/*Proportion d'échantillon détectant chacune des observations aberrantes */
/*selon les RES*/
ods graphics / noborder;
PROC SGPLOT DATA = don.tpreds&id NOBORDER;
	where STUID = "&id";
	VBAR abb / STAT=percent datalabel ;
	YAXIS label="Pourcentage d'échantillons" VALUES= (0 TO 1) ;
	XAXIS label="Valeur aberrante si RES supérieur à 1,96" ;	
RUN;

/*selon les LRTi*/
ods graphics / noborder;
PROC SGPLOT DATA = nou.infdia2 NOBORDER;
	where STUID = "&id";
	VBAR blrti / STAT=percent datalabel ;
	YAXIS label="Pourcentage d'échantillons" /*VALUES= (0 TO 1)*/ ;
	XAXIS label="Valeur aberrante si LRTj > 2,01" ;	
RUN;
%mend;
options mcompilenote= none;  
%disres (id=A, pl=10);
%disres (id=E, pl=10);

%disres (id=A, pl=2);
%disres (id=E, pl=2);

%disres (id=A, pl=1);
%disres (id=E, pl=1);

*Distribution des paramètres de variance du VSOM;
options mcompilenote=all;
%macro uni(s, pl);
ods graphics / noborder;
proc univariate data=sor.Cov&pl cibasic;
	Title1 "VSOM mixte pl=&pl";
	Title2 "Distribution de la variance de &s";
	where Covparm="&s";
	var Estimate;
   	histogram Estimate /NOFRAME 
   				VAXISLABEL="Pourcentage"/*d'échantillons*/ ;
run;
%mend;
options mcompilenote=none;

%uni(s=dA, pl=10);
%uni(s=dE, pl=10);
%uni(s=Residual, pl=10);
%uni(s=country, pl=10);
ods rtf close; 

%abb(pl=1);
ods rtf file="Z:\M2\Q4\Mémoire\sortie1.rtf";
%dis(pl=1);
%uni(s=dA, pl=1);
%uni(s=dE, pl=1);
%uni(s=Residual, pl=1);
%uni(s=country, pl=1);
ods rtf close;
 
%abb(pl=2);
ods rtf file="Z:\M2\Q4\Mémoire\sortie2.rtf";
%dis(pl=2);
%uni(s=dA, pl=2);
%uni(s=dE, pl=2);
%uni(s=Residual, pl=2);
%uni(s=country, pl=2);
ods rtf close;