diff --git a/STATS-PROC/Script_plot_density_mixture.R b/STATS-PROC/Script_plot_density_mixture.R
new file mode 100644
index 0000000..5ba62da
--- /dev/null
+++ b/STATS-PROC/Script_plot_density_mixture.R
@@ -0,0 +1,582 @@
+
+
+##############  16 July 2019 per Samuele
+
+#install.packages("mclust")
+#install.packages("flexmix")
+
+library(mclust)
+library(flexmix)
+
+
+
+#### Read data####
+rm(list=ls())
+dati<-read.csv("DB_INPUT.csv",header=T,sep=";",dec=",")
+dati<-read.csv("DB_INPUT_2cluster.csv",header=T,sep=";",dec=",")
+
+## N.B. OUPUT LAYOUT IS OPIMIZED FOR OUTCOME VARIABLE WITHIN RANGE [-100, 100], 
+### IF YOU HAVE VARIABLE WITH A DIFFERENT RANGE CONSIDER TO RESCALE IT
+
+
+ 
+
+
+### Mixture model### 
+#prima function da cui dipende numero clusters
+mod_variable<-densityMclust(dati[,4],modelnames=c("E","V"))
+mod_variable$G
+
+if (mod_variable$G==1){print("Data distribution does not present any subgroup: the mixture model estimates just one component")}
+if (mod_variable$G>=4){print("Data distribution presents many subgroups: the mixture model cannot be meaningfully adjusted for covariates")}
+
+if (mod_variable$G==3) {
+
+  #### Plot original density + points####
+
+par(bty = 'l') #"l", "7", "c", "u", or "]", "n"
+par(lwd=2.5)
+par(cex.axis=1.2)
+plot(mod_variable, what = "density",
+               data =dati[,4] ,xlab="My_variable",
+             breaks = 20,xlim=c(0,max(dati[,4]+3)))  # qui da settare i margini del plot
+
+points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+
+
+#### scatter mixture classification no covariate 
+colori<-mod_variable$classification
+colori[which(colori==2)]<-4
+colori[which(colori==1)]<-2
+plot(dati[,4],col=colori,ylab="My_variable Classification",pty=2,pch=19,cex=1.2, main="Scatter plot classification original unadjusted_variable")
+
+
+
+  #### ADJUSTMENT FOR COVARIATE APOE####
+
+mod_adj_apoe<-flexmix(dati[,4]~as.factor(dati$APOE_category_pos_neg),data=dati,k=mod_variable$G)
+# summary(mod_adj_apoe) 
+# parameters(mod_adj_apoe) 
+# table(mod_adj_apoe@cluster) 
+r_mod_adj_apoe<-refit(mod_adj_apoe) 
+#summary(r_mod_adj_apoe)
+
+
+
+### scatter mixture classification with ApoE covariate 
+colori_apoe<-clusters(mod_adj_apoe)
+colori_apoe[which(colori_apoe==2)]<-4
+colori_apoe[which(colori_apoe==1)]<-2
+plot(dati[,4],col=colori_apoe,ylab="Variable_classification_adjusted_for_ApoE",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted for APOE")
+
+
+## for obtaining density adjusted for APOE  
+pp_apoe<-parameters(mod_adj_apoe)
+
+    c1<-which(pp_apoe[1,]==min(pp_apoe[1,]))
+    c3<-which(pp_apoe[1,]==max(pp_apoe[1,]))
+    c2<-which((pp_apoe[1,]!=max(pp_apoe[1,]))&(pp_apoe[1,]!=min(pp_apoe[1,])))
+    m1<-mean(pp_apoe[1,c1]+pp_apoe[2,c1]*(as.numeric(dati$APOE_category_pos_neg)-1))
+    m2<-mean(pp_apoe[1,c2]+pp_apoe[2,c2]*(as.numeric(dati$APOE_category_pos_neg)-1))
+    m3<-mean(pp_apoe[1,c3]+pp_apoe[2,c3]*(as.numeric(dati$APOE_category_pos_neg)-1))
+    sd1<-pp_apoe[3,c1]
+    sd2<-pp_apoe[3,c2]
+    sd3<-pp_apoe[3,c3]
+    n1<-table(mod_adj_apoe@cluster)[c1]
+    n2<-table(mod_adj_apoe@cluster)[c2]
+    n3<-table(mod_adj_apoe@cluster)[c3] 
+  
+    y1_apoe<-rnorm(n1,m1,sd1)
+    y2_apoe<-rnorm(n2,m2,sd2)
+    y3_apoe<-rnorm(n3,m3,sd3)
+  
+    yy_apoe<-c(y1_apoe,y2_apoe,y3_apoe)
+  
+  
+
+
+  ### Plot Density ApoE adjusted with colored points non adjusted
+  mod_adj_apoe_dens<-densityMclust(yy_apoe,modelnames=c("E","V"))
+  plot(mod_adj_apoe_dens, what = "density",data=yy_apoe,xlab="Variable_apoe_adjusted",breaks=15, main="Density adjusted for APOE")
+
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+
+
+  
+  
+  
+  
+  
+  
+  #### ADJUSTMENT FOR  COVARIATE AGE ####
+  
+  mod_adj_age<-flexmix(dati[,4]~dati$AGE_years,data=dati,k=mod_variable$G)
+  r_mod_adj_age<-refit(mod_adj_age) 
+  
+  
+  
+  ###### scatter mixture classification with AGE covariate 
+  colori_age<-clusters(mod_adj_age)
+  colori_age[which(colori_age==2)]<-4
+  colori_age[which(colori_age==1)]<-2
+  plot(dati[,4],col=colori_age,ylab="Variable_classification_adjusted_for_Age",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Age")
+  
+  
+
+  ## for obtaining density adjusted for AGE
+  pp_age<-parameters(mod_adj_age)
+  
+  c1<-which(pp_age[1,]==min(pp_age[1,]))
+  c3<-which(pp_age[1,]==max(pp_age[1,]))
+  c2<-which((pp_age[1,]!=max(pp_age[1,]))&(pp_age[1,]!=min(pp_age[1,])))
+  m1<-mean(pp_age[1,c1]+pp_age[2,c1]*dati$AGE_years)
+  m2<-mean(pp_age[1,c2]+pp_age[2,c2]*dati$AGE_years)
+  m3<-mean(pp_age[1,c3]+pp_age[2,c3]*dati$AGE_years)
+  sd1<-pp_age[3,c1]
+  sd2<-pp_age[3,c2]
+  sd3<-pp_age[3,c3]
+  n1<-table(mod_adj_age@cluster)[c1]
+  n2<-table(mod_adj_age@cluster)[c2]
+  n3<-table(mod_adj_age@cluster)[c3] 
+  
+  y1_age<-rnorm(n1,m1,sd1)
+  y2_age<-rnorm(n2,m2,sd2)
+  y3_age<-rnorm(n3,m3,sd3)
+  
+  yy_age<-c(y1_age,y2_age,y3_age)
+  
+  
+  
+  
+  
+  ### Plot Density AGE adjusted with colored points non adjusted
+  mod_adj_age_dens<-densityMclust(yy_age,modelnames=c("E","V"))
+  plot(mod_adj_age_dens,data=yy_age, what = "density",xlab="Variable_age_adjusted",breaks=20,main="Density adjusted for Age")
+ 
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+ 
+  
+  
+  
+  
+  
+  #### ADJUSTMENT FOR COVARIATE SEX ####
+  
+  mod_adj_sex<-flexmix(dati[,4]~as.factor(dati$GENDER_0f_1m),data=dati,k=mod_variable$G)
+  r_mod_adj_sex<-refit(mod_adj_sex) 
+  
+  ### scatter mixture classification with Sex covariate 
+
+  colori_sex<-clusters(mod_adj_sex)
+  colori_sex[which(colori_sex==2)]<-4
+  colori_sex[which(colori_sex==1)]<-2
+  plot(dati[,4],col=colori_sex,ylab="Variable_classification_adjusted_for_Sex",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Sex")
+  
+ 
+  
+  ## for obtaining density adjusted for SEX
+  pp_sex<-parameters(mod_adj_sex)
+  
+  c1<-which(pp_sex[1,]==min(pp_sex[1,]))
+  c3<-which(pp_sex[1,]==max(pp_sex[1,]))
+  c2<-which((pp_sex[1,]!=max(pp_sex[1,]))&(pp_sex[1,]!=min(pp_sex[1,])))
+  m1<-mean(pp_sex[1,c1]+pp_sex[2,c1]*dati$GENDER_0f_1m)
+  m2<-mean(pp_sex[1,c2]+pp_sex[2,c2]*dati$GENDER_0f_1m)
+  m3<-mean(pp_sex[1,c3]+pp_sex[2,c3]*dati$GENDER_0f_1m)
+  sd1<-pp_age[3,c1]
+  sd2<-pp_age[3,c2]
+  sd3<-pp_age[3,c3]
+  n1<-table(mod_adj_sex@cluster)[c1]
+  n2<-table(mod_adj_sex@cluster)[c2]
+  n3<-table(mod_adj_sex@cluster)[c3] 
+  
+  y1_sex<-rnorm(n1,m1,sd1)
+  y2_sex<-rnorm(n2,m2,sd2)
+  y3_sex<-rnorm(n3,m3,sd3)
+  
+  yy_sex<-c(y1_sex,y2_sex,y3_sex)
+  
+ 
+  
+ 
+  ### Plot Density SEX adjusted with colored points non adjusted
+  
+  mod_adj_sex_dens<-densityMclust(yy_sex,modelnames=c("E","V"))
+  plot(mod_adj_sex_dens,data=yy_sex, what = "density",xlab="Variable_sex_adjusted",breaks=20,main="Density adjusted for Sex")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  
+  
+  ####  Plots density together ####
+  par(mfrow=c(2,2))
+  
+  # original plot density
+  par(bty = 'l') #"l", "7", "c", "u", or "]", "n"
+  par(lwd=2.5)
+  par(cex.axis=1.2)
+  plot(mod_variable, what = "density",
+       data =dati[,4] ,xlab="My_variable",
+       breaks = 20,xlim=c(0,max(dati[,4]+3)))  # qui da settare i margini del plot
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density ApoE adjusted with colored points non adjusted
+  mod_adj_apoe_dens<-densityMclust(yy_apoe,modelnames=c("E","V"))
+  plot(mod_adj_apoe_dens, what = "density",data=yy_apoe,xlab="Variable_apoe_adjusted",breaks=15, main="Density adjusted for APOE")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density AGE adjusted with colored points non adjusted
+  mod_adj_age_dens<-densityMclust(yy_age,modelnames=c("E","V"))
+  plot(mod_adj_age_dens,data=yy_age, what = "density",xlab="Variable_age_adjusted",breaks=20,main="Density adjusted for Age")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density SEX adjusted with colored points non adjusted
+  mod_adj_sex_dens<-densityMclust(yy_sex,modelnames=c("E","V"))
+  plot(mod_adj_sex_dens,data=yy_sex, what = "density",xlab="Variable_sex_adjusted",breaks=20,main="Density adjusted for Sex")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  #### Plots scatter together ####
+  
+  par(mfrow=c(2,2))
+  # unadjusted
+  colori<-mod_variable$classification
+  colori[which(colori==2)]<-4
+  colori[which(colori==1)]<-2
+  plot(dati[,4],col=colori,ylab="My_variable Classification",pty=2,pch=19,cex=1.2, main="Scatter plot classification original unadjusted_variable")
+  
+  #apoe adjusted
+  colori_apoe<-clusters(mod_adj_apoe)
+  colori_apoe[which(colori_apoe==2)]<-4
+  colori_apoe[which(colori_apoe==1)]<-2
+  plot(dati[,4],col=colori_apoe,ylab="Variable_classification_adjusted_for_ApoE",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted for APOE")
+  
+    # age adjusted
+  colori_age<-clusters(mod_adj_age)
+  colori_age[which(colori_age==2)]<-4
+  colori_age[which(colori_age==1)]<-2
+  plot(dati[,4],col=colori_age,ylab="Variable_classification_adjusted_for_Age",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Age")
+  
+  # sex adjusted
+  colori_sex<-clusters(mod_adj_sex)
+  colori_sex[which(colori_sex==2)]<-4
+  colori_sex[which(colori_sex==1)]<-2
+  plot(dati[,4],col=colori_sex,ylab="Variable_classification_adjusted_for_Sex",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Sex")
+  
+  
+  #### Print output della main function in un txt ####
+  out_flexmix_apoe1 <- capture.output(parameters(mod_adj_apoe) )
+  out_flexmix_apoe2 <- capture.output(table(mod_adj_apoe@cluster) )
+  out_flexmix_apoe3 <- capture.output(summary(r_mod_adj_apoe) )
+  cat("Flexmix output Adjustment for APOE_first_out ", out_flexmix_apoe1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for APOE_second_out ", out_flexmix_apoe2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for APOE_third_out ", out_flexmix_apoe3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+  out_flexmix_age1 <- capture.output(parameters(mod_adj_age) )
+  out_flexmix_age2 <- capture.output(table(mod_adj_age@cluster) )
+  out_flexmix_age3 <- capture.output(summary(r_mod_adj_age) )
+  cat("Flexmix output Adjustment for AGE_first_out ", out_flexmix_age1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for AGE_second_out ", out_flexmix_age2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for AGE_third_out ", out_flexmix_age3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+
+  out_flexmix_sex1 <- capture.output(parameters(mod_adj_sex) )
+  out_flexmix_sex2 <- capture.output(table(mod_adj_sex@cluster) )
+  out_flexmix_sex3 <- capture.output(summary(r_mod_adj_sex) )
+  cat("Flexmix output Adjustment for SEX_first_out ", out_flexmix_sex1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for SEX_second_out ", out_flexmix_sex2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for SEX_third_out ", out_flexmix_sex3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+  
+    } #END IF 3 cluster di mixture
+
+
+
+if (mod_variable$G==2) {
+  
+  #### Plot original density + points####
+  par(bty = 'l') #"l", "7", "c", "u", or "]", "n"
+  par(lwd=2.5)
+  par(cex.axis=1.2)
+  plot(mod_variable, what = "density",
+       data =dati[,4] ,xlab="My_variable",
+       breaks = 20,xlim=c(0,max(dati[,4]+3)))  # qui da settare i margini del plot
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  
+  #### scatter mixture classification no covariate 
+  colori<-mod_variable$classification
+  colori[which(colori==2)]<-4
+  colori[which(colori==1)]<-2
+  plot(dati[,4],col=colori,ylab="My_variable Classification",pty=2,pch=19,cex=1.2, main="Scatter plot classification original unadjusted_variable")
+  
+  
+  #### ADJUSTMENT FOR COVARIATE APOE####
+  
+  mod_adj_apoe<-flexmix(dati[,4]~as.factor(dati$APOE_category_pos_neg),data=dati,k=mod_variable$G)
+  r_mod_adj_apoe<-refit(mod_adj_apoe) 
+  
+  
+  ### scatter mixture classification with ApoE covariate 
+  colori_apoe<-clusters(mod_adj_apoe)
+  colori_apoe[which(colori_apoe==2)]<-4
+  colori_apoe[which(colori_apoe==1)]<-2
+  plot(dati[,4],col=colori_apoe,ylab="Variable_classification_adjusted_for_ApoE",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted for APOE")
+  
+  
+  ## for obtaining density adjusted for APOE  
+  pp_apoe<-parameters(mod_adj_apoe)
+  
+  c1<-which(pp_apoe[1,]==min(pp_apoe[1,]))
+  c2<-which(pp_apoe[1,]==max(pp_apoe[1,]))
+  #c2<-which((pp_apoe[1,]!=max(pp_apoe[1,]))&(pp_apoe[1,]!=min(pp_apoe[1,])))
+  m1<-mean(pp_apoe[1,c1]+pp_apoe[2,c1]*(as.numeric(dati$APOE_category_pos_neg)-1))
+  m2<-mean(pp_apoe[1,c2]+pp_apoe[2,c2]*(as.numeric(dati$APOE_category_pos_neg)-1))
+  #m3<-mean(pp_apoe[1,c3]+pp_apoe[2,c3]*(as.numeric(dati$APOE_category_pos_neg)-1))
+  sd1<-pp_apoe[3,c1]
+  sd2<-pp_apoe[3,c2]
+  #sd3<-pp_apoe[3,c3]
+  n1<-table(mod_adj_apoe@cluster)[c1]
+  n2<-table(mod_adj_apoe@cluster)[c2]
+  #n3<-table(mod_adj_apoe@cluster)[c3] 
+  
+  y1_apoe<-rnorm(n1,m1,sd1)
+  y2_apoe<-rnorm(n2,m2,sd2)
+  #y3_apoe<-rnorm(n3,m3,sd3)
+  
+  yy_apoe<-c(y1_apoe,y2_apoe)
+  
+  
+  
+  
+  ### Plot Density ApoE adjusted with colored points non adjusted
+  mod_adj_apoe_dens<-densityMclust(yy_apoe,modelnames=c("E","V"))
+  plot(mod_adj_apoe_dens, what = "density",data=yy_apoe,xlab="Variable_apoe_adjusted",breaks=15, main="Density adjusted for APOE")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  
+  
+  
+  #### ADJUSTMENT FOR  COVARIATE AGE ####
+  
+  mod_adj_age<-flexmix(dati[,4]~dati$AGE_years,data=dati,k=mod_variable$G)
+  r_mod_adj_age<-refit(mod_adj_age) 
+  
+  
+  ###### scatter mixture classification with AGE covariate 
+  colori_age<-clusters(mod_adj_age)
+  colori_age[which(colori_age==2)]<-4
+  colori_age[which(colori_age==1)]<-2
+  plot(dati[,4],col=colori_age,ylab="Variable_classification_adjusted_for_Age",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Age")
+  
+  
+  
+  ## for obtaining density adjusted for AGE
+  pp_age<-parameters(mod_adj_age)
+  
+  c1<-which(pp_age[1,]==min(pp_age[1,]))
+  c2<-which(pp_age[1,]==max(pp_age[1,]))
+  #c2<-which((pp_age[1,]!=max(pp_age[1,]))&(pp_age[1,]!=min(pp_age[1,])))
+  m1<-mean(pp_age[1,c1]+pp_age[2,c1]*dati$AGE_years)
+  m2<-mean(pp_age[1,c2]+pp_age[2,c2]*dati$AGE_years)
+  #m3<-mean(pp_age[1,c3]+pp_age[2,c3]*dati$AGE_years)
+  sd1<-pp_age[3,c1]
+  sd2<-pp_age[3,c2]
+  #sd3<-pp_age[3,c3]
+  n1<-table(mod_adj_age@cluster)[c1]
+  n2<-table(mod_adj_age@cluster)[c2]
+  #n3<-table(mod_adj_age@cluster)[c3] 
+  
+  y1_age<-rnorm(n1,m1,sd1)
+  y2_age<-rnorm(n2,m2,sd2)
+ # y3_age<-rnorm(n3,m3,sd3)
+  
+  yy_age<-c(y1_age,y2_age)
+  
+  
+  
+  
+  
+  ### Plot Density AGE adjusted with colored points non adjusted
+  mod_adj_age_dens<-densityMclust(yy_age,modelnames=c("E","V"))
+  plot(mod_adj_age_dens,data=yy_age, what = "density",xlab="Variable_age_adjusted",breaks=20,main="Density adjusted for Age")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  
+   
+  #### ADJUSTMENT FOR COVARIATE SEX ####
+  
+  mod_adj_sex<-flexmix(dati[,4]~as.factor(dati$GENDER_0f_1m),data=dati,k=mod_variable$G)
+  r_mod_adj_sex<-refit(mod_adj_sex) 
+  
+  ### scatter mixture classification with Sex covariate 
+  
+  colori_sex<-clusters(mod_adj_sex)
+  colori_sex[which(colori_sex==2)]<-4
+  colori_sex[which(colori_sex==1)]<-2
+  plot(dati[,4],col=colori_sex,ylab="Variable_classification_adjusted_for_Sex",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Sex")
+  
+  
+  
+  ## for obtaining density adjusted for SEX
+  pp_sex<-parameters(mod_adj_sex)
+  
+  c1<-which(pp_sex[1,]==min(pp_sex[1,]))
+  c2<-which(pp_sex[1,]==max(pp_sex[1,]))
+  #c2<-which((pp_sex[1,]!=max(pp_sex[1,]))&(pp_sex[1,]!=min(pp_sex[1,])))
+  m1<-mean(pp_sex[1,c1]+pp_sex[2,c1]*dati$GENDER_0f_1m)
+  m2<-mean(pp_sex[1,c2]+pp_sex[2,c2]*dati$GENDER_0f_1m)
+  #m3<-mean(pp_sex[1,c3]+pp_sex[2,c3]*dati$GENDER_0f_1m)
+  sd1<-pp_age[3,c1]
+  sd2<-pp_age[3,c2]
+  #sd3<-pp_age[3,c3]
+  n1<-table(mod_adj_sex@cluster)[c1]
+  n2<-table(mod_adj_sex@cluster)[c2]
+  #n3<-table(mod_adj_sex@cluster)[c3] 
+  
+  y1_sex<-rnorm(n1,m1,sd1)
+  y2_sex<-rnorm(n2,m2,sd2)
+  #y3_sex<-rnorm(n3,m3,sd3)
+  
+  yy_sex<-c(y1_sex,y2_sex)
+  
+  
+  
+  ### Plot Density SEX adjusted with colored points non adjusted
+  
+  mod_adj_sex_dens<-densityMclust(yy_sex,modelnames=c("E","V"))
+  plot(mod_adj_sex_dens,data=yy_sex, what = "density",xlab="Variable_sex_adjusted",breaks=20,main="Density adjusted for Sex")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ####  Plots density together ####
+  par(mfrow=c(2,2))
+  
+  # original plot density
+  par(bty = 'l') #"l", "7", "c", "u", or "]", "n"
+  par(lwd=2.5)
+  par(cex.axis=1.2)
+  plot(mod_variable, what = "density",
+       data =dati[,4] ,xlab="My_variable",
+       breaks = 20,xlim=c(0,max(dati[,4]+3)))  # qui da settare i margini del plot
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density ApoE adjusted with colored points non adjusted
+  mod_adj_apoe_dens<-densityMclust(yy_apoe,modelnames=c("E","V"))
+  plot(mod_adj_apoe_dens, what = "density",data=yy_apoe,xlab="Variable_apoe_adjusted",breaks=15, main="Density adjusted for APOE")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density AGE adjusted with colored points non adjusted
+  mod_adj_age_dens<-densityMclust(yy_age,modelnames=c("E","V"))
+  plot(mod_adj_age_dens,data=yy_age, what = "density",xlab="Variable_age_adjusted",breaks=20,main="Density adjusted for Age")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  ### Plot Density SEX adjusted with colored points non adjusted
+  mod_adj_sex_dens<-densityMclust(yy_sex,modelnames=c("E","V"))
+  plot(mod_adj_sex_dens,data=yy_sex, what = "density",xlab="Variable_sex_adjusted",breaks=20,main="Density adjusted for Sex")
+  
+  points(mod_variable$data[which(mod_variable$classification==1)],rep(0,sum(mod_variable$classification==1)),col=2,pch=21,lwd=2)
+  points(mod_variable$data[which(mod_variable$classification==2)],rep(0.002,sum(mod_variable$classification==2)),col=4,pch=23,lwd=2)
+  #points(mod_variable$data[which(mod_variable$classification==3)],rep(0.004,sum(mod_variable$classification==3)),col=3,pch=22,lwd=2)
+  
+  
+  #### Plots scatter together ####
+  
+  par(mfrow=c(2,2))
+  # unadjusted
+  colori<-mod_variable$classification
+  colori[which(colori==2)]<-4
+  colori[which(colori==1)]<-2
+  plot(dati[,4],col=colori,ylab="My_variable Classification",pty=2,pch=19,cex=1.2, main="Scatter plot classification original unadjusted_variable")
+  
+  #apoe adjusted
+  colori_apoe<-clusters(mod_adj_apoe)
+  colori_apoe[which(colori_apoe==2)]<-4
+  colori_apoe[which(colori_apoe==1)]<-2
+  plot(dati[,4],col=colori_apoe,ylab="Variable_classification_adjusted_for_ApoE",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted for APOE")
+  
+  # age adjusted
+  colori_age<-clusters(mod_adj_age)
+  colori_age[which(colori_age==2)]<-4
+  colori_age[which(colori_age==1)]<-2
+  plot(dati[,4],col=colori_age,ylab="Variable_classification_adjusted_for_Age",pty=2, pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Age")
+  
+  # sex adjusted
+  colori_sex<-clusters(mod_adj_sex)
+  colori_sex[which(colori_sex==2)]<-4
+  colori_sex[which(colori_sex==1)]<-2
+  plot(dati[,4],col=colori_sex,ylab="Variable_classification_adjusted_for_Sex",pty=2,pch=19,cex=1.2, main="Scatter plot classification_adjusted_for_Sex")
+  
+  #### Print output della main function in un txt ####
+  out_flexmix_apoe1 <- capture.output(parameters(mod_adj_apoe) )
+  out_flexmix_apoe2 <- capture.output(table(mod_adj_apoe@cluster) )
+  out_flexmix_apoe3 <- capture.output(summary(r_mod_adj_apoe) )
+  cat("Flexmix output Adjustment for APOE_first_out ", out_flexmix_apoe1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for APOE_second_out ", out_flexmix_apoe2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for APOE_third_out ", out_flexmix_apoe3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+  out_flexmix_age1 <- capture.output(parameters(mod_adj_age) )
+  out_flexmix_age2 <- capture.output(table(mod_adj_age@cluster) )
+  out_flexmix_age3 <- capture.output(summary(r_mod_adj_age) )
+  cat("Flexmix output Adjustment for AGE_first_out ", out_flexmix_age1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for AGE_second_out ", out_flexmix_age2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for AGE_third_out ", out_flexmix_age3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+  
+  out_flexmix_sex1 <- capture.output(parameters(mod_adj_sex) )
+  out_flexmix_sex2 <- capture.output(table(mod_adj_sex@cluster) )
+  out_flexmix_sex3 <- capture.output(summary(r_mod_adj_sex) )
+  cat("Flexmix output Adjustment for SEX_first_out ", out_flexmix_sex1, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for SEX_second_out ", out_flexmix_sex2, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  cat("Flexmix output Adjustment for SEX_third_out ", out_flexmix_sex3, file="output_flexmix_output_for_APOE_AGE_SEX.txt", sep="\n", append=TRUE)
+  
+  
+
+    }#END IF 2 cluster di mixture  and END SCRIPT
+