libname dat "M:\SOFIANE\Joane\Analyses mediation" ; 

data symptomes_dep26 ; set dat.symptomes_dep26 ; run ; 
data symptomes_dep26 ; set dat.symptomes_dep26 ; run ; 

proc freq data=symptomes_dep26 ; tables  infection_new ; run ; 
proc freq data=symptomes_dep26 ; tables  infection_new ; run ; 

*******************************************************************
*******************************************************************
*******************************************************************

************************************
************************************/

/* 1 : replicate table 2 begin  ***/
/* 1 : replicate table 2 begin  ***/
/* 1 : replicate table 2 begin  ***/


*** Independent variables ***

*FM_Sexe: sex
*FM_IncluAge: age at inclusion
*AQ_FOYVIE_Diplome_N: educational level
*AQ_FOYVIE_RevenuMont_N: household income
*tabac: smoking status
*BMI: body mass index
*sante_saprisT1: self-rated health
*diabete_snds: diabetes
*hypertension_snds: hypertension
*rank_PHQ9: PHQ-9 total score quartiles

*** Dependent variables ***

*at least: at least one persistent symptom at follow-up
*symp19bis1: fatigue
*symp2bis1: breathing difficulties
*symp20bis1: poor attention/concentration
*symp17bis1: anosmia

*** Selection variables ***

*infection_new: infected with SARS-CoV-2




*unadjested analyses ; 
%let var = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1  diabete_snds hypertension_snds rank_PHQ9 ; 

%macro unajusted ; 
%do i=1 %to 10 ; 
%let var1 = %scan(&var, &i) ; 
proc logistic data = symptomes_dep26 ; 
class  atleast  tabac BMI rank_PHQ9 /param=ref ref=first;
model atleast =    &var1
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;
%end ; 
%mend ; 

%unajusted ; 


*chez les infectés sans phQ9 = model 1*;
proc logistic data = symptomes_dep26  ;
class  atleast  tabac BMI /param=ref ref=first;
model atleast =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;


*chez les infectés + phQ9 = model 2*;
proc logistic data = symptomes_dep26  ;
class  atleast  tabac BMI rank_PHQ9 /param=ref ref=first;
model atleast =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;


*chez les infectés + phQ9 = model 3 (not in tbale 2) *;
proc logistic data = symptomes_dep26  ;
class  atleast  tabac BMI rank_PHQ9 /param=ref ref=first;
model atleast =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

/* supp test begin 1 */
/* diagnosis for multicol */
proc corr data=symptomes_dep26 nosimple;
    var FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9;
	  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;


/* Utilisation de PROC REG pour vérifier la multicollinéarité */
proc reg data=symptomes_dep26;
    model atleast = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N BMI sante_saprisT1 diabete_snds hypertension_snds / vif;
	 where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;




/* + examination of se values */

/* linearity test by plot */

proc logistic data=symptomes_dep26 outmodel=model_store noprint;
    class atleast tabac BMI rank_PHQ9 /param=ref ref=first;
    model atleast = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9 / expb;
    output out=work.predicted_values p=predicted  resdev=deviance_residuals ;
    where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1;
run;


proc sgplot data=work.predicted_values;
    scatter x=FM_IncluAge y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=FM_IncluAge y=predicted / smooth=0.5;
    xaxis label='FM_IncluAge';
    yaxis label='predicted';
    title "FM_IncluAge vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_Diplome_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_Diplome_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_Diplome_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_Diplome_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_RevenuMont_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_RevenuMont_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_RevenuMont_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_RevenuMont_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=sante_saprisT1 y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=sante_saprisT1 y=predicted / smooth=0.5;
    xaxis label='sante_saprisT1';
    yaxis label='predicted';
    title "sante_saprisT1 vs prdicted prob";
run;



/*** Normality of Residuals **/
proc univariate data=predicted_values;
    var deviance_residuals;
    histogram / normal;
run;

/* Homoscedasticity */
proc sgplot data=predicted_values;
    scatter x=predicted y=deviance_residuals;
    refline 0 / axis=y;
    xaxis label="Predicted Probabilities";
    yaxis label="Deviance Residuals";
    title "Deviance Residuals vs. Predicted Probabilities";
run;



/* supp test end 1 */

/* 1 : replicate table 2 end  ***/
/* 1 : replicate table 2 end  ***/
/* 1 : replicate table 2 end  ***/


/* 2 : replicate table 3 begin  ***/
/* 2 : replicate table 3 begin  ***/
/* 2 : replicate table 3 begin  ***/

*mediation fatigue*; 
proc logistic data = symptomes_dep26  ;
class  symp19bis1  tabac BMI rank_PHQ9/param=ref ref=first;
model symp19bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp19bis1  tabac BMI  /param=ref ref=first;
model symp19bis1 =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds 
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp19bis1  tabac BMI /param=ref ref=first;
model symp19bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds FM_Sexe*rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>.  & infection_new=1 ;
run;


/* linearity test by plot */

proc logistic data=symptomes_dep26 outmodel=model_store noprint;
    class symp19bis1 tabac BMI rank_PHQ9 /param=ref ref=first;
    model symp19bis1 = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9 / expb;
    output out=work.predicted_values p=predicted  resdev=deviance_residuals ;;
    where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1;
run;


proc sgplot data=work.predicted_values;
    scatter x=FM_IncluAge y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=FM_IncluAge y=predicted / smooth=0.5;
    xaxis label='FM_IncluAge';
    yaxis label='predicted';
    title "FM_IncluAge vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_Diplome_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_Diplome_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_Diplome_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_Diplome_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_RevenuMont_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_RevenuMont_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_RevenuMont_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_RevenuMont_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=sante_saprisT1 y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=sante_saprisT1 y=predicted / smooth=0.5;
    xaxis label='sante_saprisT1';
    yaxis label='predicted';
    title "sante_saprisT1 vs prdicted prob";
run;


/*** Normality of Residuals **/
proc univariate data=predicted_values;
    var deviance_residuals;
    histogram / normal;
run;

/* Homoscedasticity */
proc sgplot data=predicted_values;
    scatter x=predicted y=deviance_residuals;
    refline 0 / axis=y;
    xaxis label="Predicted Probabilities";
    yaxis label="Deviance Residuals";
    title "Deviance Residuals vs. Predicted Probabilities";
run;

*mediation attention*; 
proc logistic data = symptomes_dep26  ;
class  symp20bis1  tabac BMI rank_PHQ9/param=ref ref=first;
model symp20bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp20bis1  tabac BMI /param=ref ref=first;
model symp20bis1 =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp20bis1  tabac BMI /param=ref ref=first;
model symp20bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds FM_Sexe*rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>.  & infection_new=1 ;
run;


/* linearity test by plot */

proc logistic data=symptomes_dep26 outmodel=model_store noprint;
    class symp20bis1 tabac BMI rank_PHQ9 /param=ref ref=first;
    model symp20bis1 = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9 / expb;
    output out=work.predicted_values p=predicted  resdev=deviance_residuals ;;
    where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1;
run;


proc sgplot data=work.predicted_values;
    scatter x=FM_IncluAge y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=FM_IncluAge y=predicted / smooth=0.5;
    xaxis label='FM_IncluAge';
    yaxis label='predicted';
    title "FM_IncluAge vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_Diplome_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_Diplome_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_Diplome_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_Diplome_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_RevenuMont_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_RevenuMont_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_RevenuMont_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_RevenuMont_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=sante_saprisT1 y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=sante_saprisT1 y=predicted / smooth=0.5;
    xaxis label='sante_saprisT1';
    yaxis label='predicted';
    title "sante_saprisT1 vs prdicted prob";
run;

/*** Normality of Residuals **/
proc univariate data=predicted_values;
    var deviance_residuals;
    histogram / normal;
run;

/* Homoscedasticity */
proc sgplot data=predicted_values;
    scatter x=predicted y=deviance_residuals;
    refline 0 / axis=y;
    xaxis label="Predicted Probabilities";
    yaxis label="Deviance Residuals";
    title "Deviance Residuals vs. Predicted Probabilities";
run;

*mediation dyspnee*; 
proc logistic data = symptomes_dep26  ;
class  symp2bis1  tabac BMI rank_PHQ9/param=ref ref=first;
model symp2bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp2bis1  tabac BMI /param=ref ref=first;
model symp2bis1 =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp2bis1  tabac BMI /param=ref ref=first;
model symp2bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds FM_Sexe*rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>.  & infection_new=1 ;
run;


/* linearity test by plot */

proc logistic data=symptomes_dep26 outmodel=model_store noprint;
    class symp2bis1 tabac BMI rank_PHQ9 /param=ref ref=first;
    model symp2bis1 = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9 / expb;
    output out=work.predicted_values p=predicted  resdev=deviance_residuals ;;
    where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1;
run;


proc sgplot data=work.predicted_values;
    scatter x=FM_IncluAge y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=FM_IncluAge y=predicted / smooth=0.5;
    xaxis label='FM_IncluAge';
    yaxis label='predicted';
    title "FM_IncluAge vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_Diplome_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_Diplome_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_Diplome_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_Diplome_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_RevenuMont_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_RevenuMont_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_RevenuMont_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_RevenuMont_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=sante_saprisT1 y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=sante_saprisT1 y=predicted / smooth=0.5;
    xaxis label='sante_saprisT1';
    yaxis label='predicted';
    title "sante_saprisT1 vs prdicted prob";
run;

/*** Normality of Residuals **/
proc univariate data=predicted_values;
    var deviance_residuals;
    histogram / normal;
run;

/* Homoscedasticity */
proc sgplot data=predicted_values;
    scatter x=predicted y=deviance_residuals;
    refline 0 / axis=y;
    xaxis label="Predicted Probabilities";
    yaxis label="Deviance Residuals";
    title "Deviance Residuals vs. Predicted Probabilities";
run;


*mediation anosmia *; 
proc logistic data = symptomes_dep26  ;
class  symp17bis1  tabac BMI rank_PHQ9/param=ref ref=first;
model symp17bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp17bis1  tabac BMI /param=ref ref=first;
model symp17bis1 =    FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1 ;
run;

proc logistic data = symptomes_dep26  ;
class  symp17bis1  tabac BMI /param=ref ref=first;
model symp17bis1 = rank_PHQ9   FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N
 tabac BMI sante_saprisT1  diabete_snds hypertension_snds FM_Sexe*rank_PHQ9
  /expb  ;
  where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. &
sante_saprisT1>. & tabac>. & infection_new>.  & infection_new=1 ;
run;


/* linearity test by plot */

proc logistic data=symptomes_dep26 outmodel=model_store noprint;
    class symp17bis1 tabac BMI rank_PHQ9 /param=ref ref=first;
    model symp17bis1 = FM_Sexe FM_IncluAge AQ_FOYVIE_Diplome_N AQ_FOYVIE_RevenuMont_N tabac BMI sante_saprisT1 diabete_snds hypertension_snds rank_PHQ9 FM_Sexe*rank_PHQ9 / expb;
    output out=work.predicted_values p=predicted  resdev=deviance_residuals ;;
    where FM_Sexe>. & FM_IncluAge>. & AQ_FOYVIE_Diplome_N>. & AQ_FOYVIE_RevenuMont_N>. & BMI>. & sante_saprisT1>. & tabac>. & infection_new>. & infection_new=1;
run;


proc sgplot data=work.predicted_values;
    scatter x=FM_IncluAge y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=FM_IncluAge y=predicted / smooth=0.5;
    xaxis label='FM_IncluAge';
    yaxis label='predicted';
    title "FM_IncluAge vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_Diplome_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_Diplome_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_Diplome_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_Diplome_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=AQ_FOYVIE_RevenuMont_N y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=AQ_FOYVIE_RevenuMont_N y=predicted / smooth=0.5;
    xaxis label='AQ_FOYVIE_RevenuMont_N';
    yaxis label='predicted';
    title "AQ_FOYVIE_RevenuMont_N vs prdicted prob";
run;
proc sgplot data=work.predicted_values;
    scatter x=sante_saprisT1 y=predicted / markerattrs=(symbol=CircleFilled) transparency=0.5;
	loess x=sante_saprisT1 y=predicted / smooth=0.5;
    xaxis label='sante_saprisT1';
    yaxis label='predicted';
    title "sante_saprisT1 vs prdicted prob";
run;

/*** Normality of Residuals **/
proc univariate data=predicted_values;
    var deviance_residuals;
    histogram / normal;
run;

/* Homoscedasticity */
proc sgplot data=predicted_values;
    scatter x=predicted y=deviance_residuals;
    refline 0 / axis=y;
    xaxis label="Predicted Probabilities";
    yaxis label="Deviance Residuals";
    title "Deviance Residuals vs. Predicted Probabilities";
run;

/* 2 : replicate table 3 end  ***/
/* 2 : replicate table 3 end  ***/
/* 2 : replicate table 3 end  ***/


*merging to keep the sero positive in data*;



proc sort data=symptomes_dep26;
by proj_isp;run;

proc sort data=testcedric;
by proj_isp;
run; 

data symptomes_dep26;
set symptomes_dep26;
if interpretation_elisa_spike='pos' or sante4q02m_1=1 then infection_newbis=1;
if interpretation_elisa_spike='neg' and sante4q02m_1 ne 1 then infection_newbis=0;
run; 

*mediation analyses*; 

data symptomes_dep28 (keep=PROJ_ISP v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16);
set symptomes_dep26;
rename 
fm_sexe=v1
fm_incluage=v2
AQ_FOYVIE_Diplome_N=v3
AQ_FOYVIE_RevenuMont_N=v4
tabac=v5
BMI=v6
sante_saprisT1=v7
diabete_snds=v8
hypertension_snds=v9
atleast=v10
rank_PHQ9=v11
symp2bis1=v12
symp19bis1=v13
symp20bis1=v14
symp17bis1=v15
infection_newbis=v16

;
label 
fm_incluage='age'
fm_sexe='sexe'
AQ_FOYVIE_RevenuMont_N='income'
AQ_FOYVIE_Diplome_N='diploma'
tabac='smoking'
BMI='body mass index'
sante_saprisT1='self-rated health pandemic'
diabete_snds='diabetes'
atleast='at least one symptom'
hypertension_snds='hypertension'
rank_PHQ9='quartiles phq9'
symp2bis1='dyspnee'
symp19bis1='fatigue'
symp20bis1='attention'
symp17bis1='anosmia'
infection_new='infection covid'
;
run;




proc calis data=symptomes_dep28 covariance residual modification TOTEFF method=MLM ;
  lineqs
 
v10= pv10v11 v11 + pv10v1 v1 +  pv10v2 v2+ pv10v3 v3  + pv10v4 v4 + pv10v5 v5 + pv10v6 v6 +
pv10v7 v7 + pv10v8 v8 + pv10v9 v9
+ e10,


v11=pv11v1 v1 + e11; 

std
	e11=vare11,
	e10=vare10,


	v1=varv1,
	v2=varv2,
	v3=varv3,
	v4=varv4,
	v5=varv5,
	v6=varv6,
	v7=varv7,
	v8=varv8,
	v9=varv9;

cov

	v2 v3=cv2v3,
	v2 v4=cv2v4,
	v2 v5=cv2v5,
	v2 v6=cv2v6,
	v2 v7=cv2v7,
	v2 v8=cv2v8,
	v2 v9=cv2v9,
	v2 v1=cv2v1,

	v3 v4=cv3v4,
	v3 v5=cv3v5,
	v3 v6=cv3v6,
	v3 v7=cv3v7,
	v3 v8=cv3v8,
	v3 v9=cv1v9,
	v3 v1=cv1v1,
	
	v4 v5=cv4v5,
	v4 v6=cv4v6,
	v4 v7=cv4v7,
	v4 v8=cv4v8,
	v4 v9=cv4v9,
	v4 v1=cv4v1,
	
	v5 v6=cv5v6,
	v5 v7=cv5v7,
	v5 v8=cv5v8,
	v5 v9=cv5v9,
	v5 v1=cv5v1,
	
    v6 v7=cv6v7,
	v6 v8=cv6v8,
	v6 v9=cv6v9,
	v6 v1=cv6v1,
	
	v7 v8=cv7v8,
	v7 v9=cv7v9,
	v7 v1=cv7v1,
	
	v8 v9=cv8v9,
	v8 v1=cv8v1,

	v9 V1=cv9v1
	;

var  V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 ;

 where  V1>. & v2>. & v3>. & v4>. & v5>. & v6>. 
& v7>.
 & v8>. & v9>. &
v10>.  & v11>. & v16=1;

run;




*fatigue*;

proc calis data=symptomes_dep28 covariance residual modification TOTEFF method=MLM ;
  lineqs
 
v13= pv13v11 v11 + pv13v1 v1 +  pv13v2 v2+ pv13v3 v3  + pv13v4 v4 + pv13v5 v5 + pv13v6 v6 +
pv13v7 v7 + pv13v8 v8 + pv13v9 v9
+ e13,


v11=pv11v1 v1 + e11; 

std
	e11=vare11,
	e13=vare13,


	v1=varv1,
	v2=varv2,
	v3=varv3,
	v4=varv4,
	v5=varv5,
	v6=varv6,
	v7=varv7,
	v8=varv8,
	v9=varv9;

cov

	v2 v3=cv2v3,
	v2 v4=cv2v4,
	v2 v5=cv2v5,
	v2 v6=cv2v6,
	v2 v7=cv2v7,
	v2 v8=cv2v8,
	v2 v9=cv2v9,
	v2 v1=cv2v1,

	v3 v4=cv3v4,
	v3 v5=cv3v5,
	v3 v6=cv3v6,
	v3 v7=cv3v7,
	v3 v8=cv3v8,
	v3 v9=cv1v9,
	v3 v1=cv1v1,
	
	v4 v5=cv4v5,
	v4 v6=cv4v6,
	v4 v7=cv4v7,
	v4 v8=cv4v8,
	v4 v9=cv4v9,
	v4 v1=cv4v1,
	
	v5 v6=cv5v6,
	v5 v7=cv5v7,
	v5 v8=cv5v8,
	v5 v9=cv5v9,
	v5 v1=cv5v1,
	
    v6 v7=cv6v7,
	v6 v8=cv6v8,
	v6 v9=cv6v9,
	v6 v1=cv6v1,
	
	v7 v8=cv7v8,
	v7 v9=cv7v9,
	v7 v1=cv7v1,
	
	v8 v9=cv8v9,
	v8 v1=cv8v1,

	v9 V1=cv9v1
	;

var  V1 V2 V3 V4 V5 V6 V7 V8 V9 V13 V11 ;

 where v16=1;

run;


symptom19=fatigue
symptom20=attention
symptom17=anosmia
symptome 2=respi

symp2bis1
;

*dyspnee*;
proc calis data=symptomes_dep28 covariance residual modification TOTEFF method=MLM ;
  lineqs
 
v12= pv12v11 v11 + pv12v1 v1 +  pv12v2 v2+ pv12v3 v3  + pv12v4 v4 + pv12v5 v5 + pv12v6 v6 +
pv12v7 v7 + pv12v8 v8 + pv12v9 v9
+ e13,



v11=pv11v1 v1 + e11; 
std
	e11=vare11,
	e13=vare13,


	v1=varv1,
	v2=varv2,
	v3=varv3,
	v4=varv4,
	v5=varv5,
	v6=varv6,
	v7=varv7,
	v8=varv8,
	v9=varv9;

cov

	v2 v3=cv2v3,
	v2 v4=cv2v4,
	v2 v5=cv2v5,
	v2 v6=cv2v6,
	v2 v7=cv2v7,
	v2 v8=cv2v8,
	v2 v9=cv2v9,
	v2 v1=cv2v1,

	v3 v4=cv3v4,
	v3 v5=cv3v5,
	v3 v6=cv3v6,
	v3 v7=cv3v7,
	v3 v8=cv3v8,
	v3 v9=cv1v9,
	v3 v1=cv1v1,
	
	v4 v5=cv4v5,
	v4 v6=cv4v6,
	v4 v7=cv4v7,
	v4 v8=cv4v8,
	v4 v9=cv4v9,
	v4 v1=cv4v1,
	
	v5 v6=cv5v6,
	v5 v7=cv5v7,
	v5 v8=cv5v8,
	v5 v9=cv5v9,
	v5 v1=cv5v1,
	
    v6 v7=cv6v7,
	v6 v8=cv6v8,
	v6 v9=cv6v9,
	v6 v1=cv6v1,
	
	v7 v8=cv7v8,
	v7 v9=cv7v9,
	v7 v1=cv7v1,
	
	v8 v9=cv8v9,
	v8 v1=cv8v1,

	v9 V1=cv9v1
	;

var  V1 V2 V3 V4 V5 V6 V7 V8 V9 V12 V11 ;

 where v16=1;

run;

*v14 attention*;
proc calis data=symptomes_dep28 covariance residual modification TOTEFF method=MLM ;
  lineqs
 
v14= pv14v11 v11 + pv14v1 v1 +  pv14v2 v2+ pv14v3 v3  + pv14v4 v4 + pv14v5 v5 + pv14v6 v6 +
pv14v7 v7 + pv14v8 v8 + pv14v9 v9
+ e13,


v11=pv11v1 v1 + e11; 

std
	e11=vare11,
	e13=vare13,


	v1=varv1,
	v2=varv2,
	v3=varv3,
	v4=varv4,
	v5=varv5,
	v6=varv6,
	v7=varv7,
	v8=varv8,
	v9=varv9;

cov

	v2 v3=cv2v3,
	v2 v4=cv2v4,
	v2 v5=cv2v5,
	v2 v6=cv2v6,
	v2 v7=cv2v7,
	v2 v8=cv2v8,
	v2 v9=cv2v9,
	v2 v1=cv2v1,

	v3 v4=cv3v4,
	v3 v5=cv3v5,
	v3 v6=cv3v6,
	v3 v7=cv3v7,
	v3 v8=cv3v8,
	v3 v9=cv1v9,
	v3 v1=cv1v1,
	
	v4 v5=cv4v5,
	v4 v6=cv4v6,
	v4 v7=cv4v7,
	v4 v8=cv4v8,
	v4 v9=cv4v9,
	v4 v1=cv4v1,
	
	v5 v6=cv5v6,
	v5 v7=cv5v7,
	v5 v8=cv5v8,
	v5 v9=cv5v9,
	v5 v1=cv5v1,
	
    v6 v7=cv6v7,
	v6 v8=cv6v8,
	v6 v9=cv6v9,
	v6 v1=cv6v1,
	
	v7 v8=cv7v8,
	v7 v9=cv7v9,
	v7 v1=cv7v1,
	
	v8 v9=cv8v9,
	v8 v1=cv8v1,

	v9 V1=cv9v1
	;

var  V1 V2 V3 V4 V5 V6 V7 V8 V9 V14 V11 ;
where v16=1;
 

run;


*v15 anosmie*; 

proc calis data=symptomes_dep28 covariance residual modification TOTEFF method=MLM ;
  lineqs
 
v15= pv15v11 v11 + pv15v1 v1 +  pv15v2 v2+ pv15v3 v3  + pv15v4 v4 + pv15v5 v5 + pv15v6 v6 +
pv15v7 v7 + pv15v8 v8 + pv15v9 v9
+ e13,



v11=pv11v1 v1 + e11; 
std
	e11=vare11,
	e13=vare13,


	v1=varv1,
	v2=varv2,
	v3=varv3,
	v4=varv4,
	v5=varv5,
	v6=varv6,
	v7=varv7,
	v8=varv8,
	v9=varv9;

cov

	v2 v3=cv2v3,
	v2 v4=cv2v4,
	v2 v5=cv2v5,
	v2 v6=cv2v6,
	v2 v7=cv2v7,
	v2 v8=cv2v8,
	v2 v9=cv2v9,
	v2 v1=cv2v1,

	v3 v4=cv3v4,
	v3 v5=cv3v5,
	v3 v6=cv3v6,
	v3 v7=cv3v7,
	v3 v8=cv3v8,
	v3 v9=cv1v9,
	v3 v1=cv1v1,
	
	v4 v5=cv4v5,
	v4 v6=cv4v6,
	v4 v7=cv4v7,
	v4 v8=cv4v8,
	v4 v9=cv4v9,
	v4 v1=cv4v1,
	
	v5 v6=cv5v6,
	v5 v7=cv5v7,
	v5 v8=cv5v8,
	v5 v9=cv5v9,
	v5 v1=cv5v1,
	
    v6 v7=cv6v7,
	v6 v8=cv6v8,
	v6 v9=cv6v9,
	v6 v1=cv6v1,
	
	v7 v8=cv7v8,
	v7 v9=cv7v9,
	v7 v1=cv7v1,
	
	v8 v9=cv8v9,
	v8 v1=cv8v1,

	v9 V1=cv9v1
	;

var  V1 V2 V3 V4 V5 V6 V7 V8 V9 V15 V11 ;

 
where v16=1;
run;

*descriptive stats*;