load(file="Analysen_Section2/powercase/gini0.05power.RData")
load(file="Analysen_Section2/powercase/gini0.1power.RData")
load(file="Analysen_Section2/powercase/gini0.25power.RData")
load(file="Analysen_Section2/powercase/gini0.4power.RData")

load(file="Analysen_Section2/powercase/gini0.05power1.RData")
load(file="Analysen_Section2/powercase/gini0.1power1.RData")
load(file="Analysen_Section2/powercase/gini0.25power1.RData")
load(file="Analysen_Section2/powercase/gini0.4power1.RData")

load(file="Analysen_Section2/powercase/gini0.05power2.RData")
load(file="Analysen_Section2/powercase/gini0.1power2.RData")
load(file="Analysen_Section2/powercase/gini0.25power2.RData")
load(file="Analysen_Section2/powercase/gini0.4power2.RData")

load(file="Analysen_Section2/powercase/perm0.05power.RData")
load(file="Analysen_Section2/powercase/perm0.1power.RData")
load(file="Analysen_Section2/powercase/perm0.25power.RData")
load(file="Analysen_Section2/powercase/perm0.4power.RData")

load(file="Analysen_Section2/powercase/perm0.05power1.RData")
load(file="Analysen_Section2/powercase/perm0.1power1.RData")
load(file="Analysen_Section2/powercase/perm0.25power1.RData")
load(file="Analysen_Section2/powercase/perm0.4power1.RData")

load(file="Analysen_Section2/powercase/perm0.05power2.RData")
load(file="Analysen_Section2/powercase/perm0.1power2.RData")
load(file="Analysen_Section2/powercase/perm0.25power2.RData")
load(file="Analysen_Section2/powercase/perm0.4power2.RData")


load(file="Analysen_Section2/powercase/cforest0.05power.RData")
load(file="Analysen_Section2/powercase/cforest0.1power.RData")
load(file="Analysen_Section2/powercase/cforest0.25power.RData")
load(file="Analysen_Section2/powercase/cforest0.4power.RData")

load(file="Analysen_Section2/powercase/cforest0.05power1.RData")
load(file="Analysen_Section2/powercase/cforest0.1power1.RData")
load(file="Analysen_Section2/powercase/cforest0.25power1.RData")
load(file="Analysen_Section2/powercase/cforest0.4power1.RData")

load(file="Analysen_Section2/powercase/cforest0.05power2.RData")
load(file="Analysen_Section2/powercase/cforest0.1power2.RData")
load(file="Analysen_Section2/powercase/cforest0.25power2.RData")
load(file="Analysen_Section2/powercase/cforest0.4power2.RData")



####

load(file="Analysen_Section2/powercase_n10000/gini0.05powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.1powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.25powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.4powern10000.RData")

load(file="Analysen_Section2/powercase_n10000/gini0.05power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.1power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.25power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.4power1n10000.RData")

load(file="Analysen_Section2/powercase_n10000/gini0.05power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.1power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.25power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/gini0.4power2n10000.RData")

load(file="Analysen_Section2/powercase_n10000/perm0.05powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.1powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.25powern10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.4powern10000.RData")

load(file="Analysen_Section2/powercase_n10000/perm0.05power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.1power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.25power1n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.4power1n10000.RData")

load(file="Analysen_Section2/powercase_n10000/perm0.05power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.1power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.25power2n10000.RData")
load(file="Analysen_Section2/powercase_n10000/perm0.4power2n10000.RData")





pdf("Figure4.pdf",width=8,height=8)
par(mfrow=c(2,2),las=2)

# ROC n=500
giniinfo3<-matrix(0,100,4)
giniinfo1.5<-matrix(0,100,4)
gininoninfo<-matrix(0,100,4)
perminfo3<-matrix(0,100,4)
perminfo1.5<-matrix(0,100,4)
permnoninfo<-matrix(0,100,4)
cfinfo3<-matrix(0,100,4)
cfinfo1.5<-matrix(0,100,4)
cfnoninfo<-matrix(0,100,4)

for (i in 1:100)
 {
 giniinfo3[i,]<-c(gini0.05power1[i],gini0.1power1[i],gini0.25power1[i],gini0.4power1[i])
 giniinfo1.5[i,]<- c(gini0.05power2[i],gini0.1power2[i],gini0.25power2[i],gini0.4power2[i])
 gininoninfo[i,]<-c(gini0.05power[i],gini0.1power[i],gini0.25power[i],gini0.4power[i])
 perminfo3[i,]<-c(perm0.05power1[i],perm0.1power1[i],perm0.25power1[i],perm0.4power1[i])
 perminfo1.5[i,]<- c(perm0.05power2[i],perm0.1power2[i],perm0.25power2[i],perm0.4power2[i])
 permnoninfo[i,]<-c(perm0.05power[i],perm0.1power[i],perm0.25power[i],perm0.4power[i])
 cfinfo3[i,]<-c(cforest0.05power1[i],cforest0.1power1[i],cforest0.25power1[i],cforest0.4power1[i])
 cfinfo1.5[i,]<- c(cforest0.05power2[i],cforest0.1power2[i],cforest0.25power2[i],cforest0.4power2[i])
 cfnoninfo[i,]<-c(cforest0.05power[i],cforest0.1power[i],cforest0.25power[i],cforest0.4power[i]) 
 }
library(ROCR)



##############################
## 1. Plot: n=500, all MAFs ##
##############################
predgini3<-prediction(predictions=c(as.numeric(giniinfo3),as.numeric(gininoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfgini3<-performance(predgini3,"tpr","fpr")
aucgini3<-performance(predgini3,measure="auc")@y.values[[1]]

predgini1.5<-prediction(predictions=c(as.numeric(giniinfo1.5),as.numeric(gininoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfgini1.5<-performance(predgini1.5,"tpr","fpr")
aucgini1.5<-performance(predgini1.5,measure="auc")@y.values[[1]]

predperm3<-prediction(predictions=c(as.numeric(perminfo3),as.numeric(permnoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfperm3<-performance(predperm3,"tpr","fpr")
aucperm3<-performance(predperm3,measure="auc")@y.values[[1]]

predperm1.5<-prediction(predictions=c(as.numeric(perminfo1.5),as.numeric(permnoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfperm1.5<-performance(predperm1.5,"tpr","fpr")
aucperm1.5<-performance(predperm1.5,measure="auc")@y.values[[1]]

predcf3<-prediction(predictions=c(as.numeric(cfinfo3),as.numeric(cfnoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfcf3<-performance(predcf3,"tpr","fpr")
auccf3<-performance(predcf3,measure="auc")@y.values[[1]]

predcf1.5<-prediction(predictions=c(as.numeric(cfinfo1.5),as.numeric(cfnoninfo)),labels=c(rep(1,100*4),rep(0,400)))
perfcf1.5<-performance(predcf1.5,"tpr","fpr")
auccf1.5<-performance(predcf1.5,measure="auc")@y.values[[1]]


plot(perfperm1.5@x.values[[1]],perfperm1.5@y.values[[1]],col=2,type="l",ylab="sensitivity",xlab="1-specificity",main="all MAFs\n n=500" ,lty=2)
points(perfgini1.5@x.values[[1]],perfgini1.5@y.values[[1]],col=1,type="l" ,lty=2)
points(perfcf1.5@x.values[[1]],perfcf1.5@y.values[[1]],col=3,type="l" ,lty=2)
points(perfperm3@x.values[[1]],perfperm3@y.values[[1]],col=2,type="l",lty=1)
points(perfgini3@x.values[[1]],perfgini3@y.values[[1]],col=1,type="l",lty=1)
points(perfcf3@x.values[[1]],perfcf3@y.values[[1]],col=3,type="l",lty=1)
legend(x=0.40,y=0.47, legend=c( "gini - OR=3","perm rf - OR=3","perm cf - OR=3","gini - OR=1.5","perm rf - OR=1.5","perm cf - OR=1.5"),lty=c(1,1,1,2,2,2),col=c(1,2,3,1,2,3))

##################################
## 2. Plot: n=500, MAF=0.05,0.4 ##
##################################
predgini3MAF0.05MAF0.4<-prediction(predictions=c(gini0.05power1,gini0.4power1,gini0.05power,gini0.4power),labels=c(rep(1,200),rep(0,200)))
perfgini3MAF0.05MAF0.4<-performance(predgini3MAF0.05MAF0.4,"tpr","fpr")
aucgini3MAF0.05MAF0.4<-performance(predgini3MAF0.05MAF0.4,measure="auc")@y.values[[1]]

predperm3MAF0.05MAF0.4<-prediction(predictions=c(perm0.05power1,perm0.4power1,perm0.05power,perm0.4power),labels=c(rep(1,200),rep(0,200)))
perfperm3MAF0.05MAF0.4<-performance(predperm3MAF0.05MAF0.4,"tpr","fpr")
aucperm3MAF0.05MAF0.4<-performance(predperm3MAF0.05MAF0.4,measure="auc")@y.values[[1]]

predcf3MAF0.05MAF0.4<-prediction(predictions=c(cforest0.05power1,cforest0.4power1,cforest0.05power,cforest0.4power),labels=c(rep(1,200),rep(0,200)))
perfcf3MAF0.05MAF0.4<-performance(predcf3MAF0.05MAF0.4,"tpr","fpr")
auccf3MAF0.05MAF0.4<-performance(predcf3MAF0.05MAF0.4,measure="auc")@y.values[[1]]

predgini1.5MAF0.05MAF0.4<-prediction(predictions=c(gini0.05power2,gini0.4power2,gini0.05power,gini0.4power),labels=c(rep(1,200),rep(0,200)))
perfgini1.5MAF0.05MAF0.4<-performance(predgini1.5MAF0.05MAF0.4,"tpr","fpr")
aucgini1.5MAF0.05MAF0.4<-performance(predgini1.5MAF0.05MAF0.4,measure="auc")@y.values[[1]]

predperm1.5MAF0.05MAF0.4<-prediction(predictions=c(perm0.05power2,perm0.4power2,perm0.05power,perm0.4power),labels=c(rep(1,200),rep(0,200)))
perfperm1.5MAF0.05MAF0.4<-performance(predperm1.5MAF0.05MAF0.4,"tpr","fpr")
aucperm1.5MAF0.05MAF0.4<-performance(predperm1.5MAF0.05MAF0.4,measure="auc")@y.values[[1]]

predcf1.5MAF0.05MAF0.4<-prediction(predictions=c(cforest0.05power2,cforest0.4power2,cforest0.05power,cforest0.4power),labels=c(rep(1,200),rep(0,200)))
perfcf1.5MAF0.05MAF0.4<-performance(predcf1.5MAF0.05MAF0.4,"tpr","fpr")
auccf1.5MAF0.05MAF0.4<-performance(predcf1.5MAF0.05MAF0.4,measure="auc")@y.values[[1]]

plot(perfperm1.5MAF0.05MAF0.4@x.values[[1]],perfperm1.5MAF0.05MAF0.4@y.values[[1]],col=2,type="l",ylab="sensitivity",xlab="1-specificity",main="MAF=0.05 and MAF=0.4\n n=500" ,lty=2)
points(perfgini1.5MAF0.05MAF0.4@x.values[[1]],perfgini1.5MAF0.05MAF0.4@y.values[[1]],col=1,type="l" ,lty=2)
points(perfcf1.5MAF0.05MAF0.4@x.values[[1]],perfcf1.5MAF0.05MAF0.4@y.values[[1]],col=3,type="l" ,lty=2)
points(perfperm3MAF0.05MAF0.4@x.values[[1]],perfperm3MAF0.05MAF0.4@y.values[[1]],col=2,type="l",lty=1)
points(perfgini3MAF0.05MAF0.4@x.values[[1]],perfgini3MAF0.05MAF0.4@y.values[[1]],col=1,type="l",lty=1)
points(perfcf3MAF0.05MAF0.4@x.values[[1]],perfcf3MAF0.05MAF0.4@y.values[[1]],col=3,type="l",lty=1)
legend(x=0.40,y=0.47, legend=c( "gini - OR=3","perm rf - OR=3","perm cf - OR=3","gini - OR=1.5","perm rf - OR=1.5","perm cf - OR=1.5"),lty=c(1,1,1,2,2,2),col=c(1,2,3,1,2,3))


################################
## 3. Plot: n=10000, all MAFs ##
################################

# ROC n=10000
giniinfo3n10000<-matrix(0,100,4)
giniinfo1.5n10000<-matrix(0,100,4)
gininoninfon10000<-matrix(0,100,4)
perminfo3n10000<-matrix(0,100,4)
perminfo1.5n10000<-matrix(0,100,4)
permnoninfon10000<-matrix(0,100,4)

for (i in 1:100)
 {
 giniinfo3n10000[i,]<-c(gini0.05power1n10000[i],gini0.1power1n10000[i],gini0.25power1n10000[i],gini0.4power1n10000[i])
 giniinfo1.5n10000[i,]<- c(gini0.05power2n10000[i],gini0.1power2n10000[i],gini0.25power2n10000[i],gini0.4power2n10000[i])
 gininoninfon10000[i,]<-c(gini0.05powern10000[i],gini0.1powern10000[i],gini0.25powern10000[i],gini0.4powern10000[i])
 perminfo3n10000[i,]<-c(perm0.05power1n10000[i],perm0.1power1n10000[i],perm0.25power1n10000[i],perm0.4power1n10000[i])
 perminfo1.5n10000[i,]<- c(perm0.05power2n10000[i],perm0.1power2n10000[i],perm0.25power2n10000[i],perm0.4power2n10000[i])
 permnoninfon10000[i,]<-c(perm0.05powern10000[i],perm0.1powern10000[i],perm0.25powern10000[i],perm0.4powern10000[i])
 }
library(ROCR)


predgini3n10000<-prediction(predictions=c(as.numeric(giniinfo3n10000),as.numeric(gininoninfon10000)),labels=c(rep(1,100*4),rep(0,400)))
perfgini3n10000<-performance(predgini3n10000,"tpr","fpr")
aucgini3n10000<-performance(predgini3n10000,measure="auc")@y.values[[1]]

predgini1.5n10000<-prediction(predictions=c(as.numeric(giniinfo1.5n10000),as.numeric(gininoninfon10000)),labels=c(rep(1,100*4),rep(0,400)))
perfgini1.5n10000<-performance(predgini1.5n10000,"tpr","fpr")
aucgini1.5n10000<-performance(predgini1.5n10000,measure="auc")@y.values[[1]]

predperm3n10000<-prediction(predictions=c(as.numeric(perminfo3n10000),as.numeric(permnoninfon10000)),labels=c(rep(1,100*4),rep(0,400)))
perfperm3n10000<-performance(predperm3n10000,"tpr","fpr")
aucperm3n10000<-performance(predperm3n10000,measure="auc")@y.values[[1]]

predperm1.5n10000<-prediction(predictions=c(as.numeric(perminfo1.5n10000),as.numeric(permnoninfon10000)),labels=c(rep(1,100*4),rep(0,400)))
perfperm1.5n10000<-performance(predperm1.5n10000,"tpr","fpr")
aucperm1.5n10000<-performance(predperm1.5n10000,measure="auc")@y.values[[1]]


plot(perfperm1.5n10000@x.values[[1]],perfperm1.5n10000@y.values[[1]],col=2,type="l",ylab="sensitivity",xlab="1-specificity",main="all MAFs\n n=10000" ,lty=2)
points(perfgini1.5n10000@x.values[[1]],perfgini1.5n10000@y.values[[1]],col=1,type="l" ,lty=2)
points(perfperm3n10000@x.values[[1]],perfperm3n10000@y.values[[1]],col=2,type="l",lty=1)
points(perfgini3n10000@x.values[[1]],perfgini3n10000@y.values[[1]],col=1,type="l",lty=1)
legend(x=0.40,y=0.47, legend=c( "gini - OR=3","perm rf - OR=3","gini - OR=1.5","perm rf - OR=1.5"),lty=c(1,1,2,2),col=c(1,2,1,2))

####################################
## 4. Plot: n=10000, MAF=0.05,0.4 ##
####################################

giniinfo1.5n10000

predgini3MAF0.05MAF0.4n10000<-prediction(predictions=c(gini0.05power1n10000,gini0.4power1n10000,gini0.05powern10000,gini0.4powern10000),labels=c(rep(1,200),rep(0,200)))
perfgini3MAF0.05MAF0.4n10000<-performance(predgini3MAF0.05MAF0.4n10000,"tpr","fpr")
aucgini3MAF0.05MAF0.4n10000<-performance(predgini3MAF0.05MAF0.4n10000,measure="auc")@y.values[[1]]

predperm3MAF0.05MAF0.4n10000<-prediction(predictions=c(perm0.05power1n10000,perm0.4power1n10000,perm0.05powern10000,perm0.4powern10000),labels=c(rep(1,200),rep(0,200)))
perfperm3MAF0.05MAF0.4n10000<-performance(predperm3MAF0.05MAF0.4n10000,"tpr","fpr")
aucperm3MAF0.05MAF0.4n10000<-performance(predperm3MAF0.05MAF0.4n10000,measure="auc")@y.values[[1]]

predgini1.5MAF0.05MAF0.4n10000<-prediction(predictions=c(gini0.05power2n10000,gini0.4power2n10000,gini0.05powern10000,gini0.4powern10000),labels=c(rep(1,200),rep(0,200)))
perfgini1.5MAF0.05MAF0.4n10000<-performance(predgini1.5MAF0.05MAF0.4n10000,"tpr","fpr")
aucgini1.5MAF0.05MAF0.4n10000<-performance(predgini1.5MAF0.05MAF0.4n10000,measure="auc")@y.values[[1]]

predperm1.5MAF0.05MAF0.4n10000<-prediction(predictions=c(perm0.05power2n10000,perm0.4power2n10000,perm0.05powern10000,perm0.4powern10000),labels=c(rep(1,200),rep(0,200)))
perfperm1.5MAF0.05MAF0.4n10000<-performance(predperm1.5MAF0.05MAF0.4n10000,"tpr","fpr")
aucperm1.5MAF0.05MAF0.4n10000<-performance(predperm1.5MAF0.05MAF0.4n10000,measure="auc")@y.values[[1]]

plot(perfperm1.5MAF0.05MAF0.4n10000@x.values[[1]],perfperm1.5MAF0.05MAF0.4n10000@y.values[[1]],col=2,type="l",ylab="sensitivity",xlab="1-specificity",main="MAF=0.05 and MAF=0.4\n n=10000" ,lty=2)
points(perfgini1.5MAF0.05MAF0.4n10000@x.values[[1]],perfgini1.5MAF0.05MAF0.4n10000@y.values[[1]],col=1,type="l" ,lty=2)
points(perfperm3MAF0.05MAF0.4n10000@x.values[[1]],perfperm3MAF0.05MAF0.4n10000@y.values[[1]],col=2,type="l",lty=1)
points(perfgini3MAF0.05MAF0.4n10000@x.values[[1]],perfgini3MAF0.05MAF0.4n10000@y.values[[1]],col=1,type="l",lty=1)
legend(x=0.40,y=0.47, legend=c( "gini - OR=3","perm rf - OR=3","gini - OR=1.5","perm rf - OR=1.5"),lty=c(1,1,2,2),col=c(1,2,1,2))


dev.off()

