Egalware/admodelling.org
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
develop
admodelling.org/STATS-PROC/Script_plot_density_mixture.R
T
marco 49b1bc823d Carica file su 'STATS-PROC' 
Script R
2019-10-17 10:08:47 +02:00

583 lines
26 KiB
R
Raw Permalink Blame History

############## 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
Reference in New Issue View Git Blame Copy Permalink