##############  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