#setwd("Z:/Manuskripte/Papers/Preparation/Gini/Analysen_Section2")

########################################################
# Bias in variable importance measure: null case study #
########################################################
library(randomForest)
library(party)
source("fun.r")

gini0.05n10000<-c()
gini0.1n10000<-c()
gini0.25n10000<-c()
gini0.4n10000<-c()

perm0.05n10000<-c()
perm0.1n10000<-c()
perm0.25n10000<-c()
perm0.4n10000<-c()

cforest0.05n10000<-c()
cforest0.1n10000<-c()
cforest0.25n10000<-c()
cforest0.4n10000<-c()

for (i in 1:100)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,n.obs=10000,maf=rep(c(0.05,0.1,0.25,0.4),each=50),beta0=0)
 my.rf<-randomForest(x=my.data[,-1],y=my.data[,1],maxnodes=2,importance=TRUE)
 my.importance<-importance(my.rf,scale=FALSE)
 gini0.05n10000<-c(gini0.05n10000,my.importance[1:50,4])
 gini0.1n10000<-c(gini0.1n10000,my.importance[51:100,4])
 gini0.25n10000<-c(gini0.25n10000,my.importance[101:150,4])
 gini0.4n10000<-c(gini0.4n10000,my.importance[151:200,4])
 perm0.05n10000<-c(perm0.05n10000,my.importance[1:50,3])
 perm0.1n10000<-c(perm0.1n10000,my.importance[51:100,3])
 perm0.25n10000<-c(perm0.25n10000,my.importance[101:150,3])
 perm0.4n10000<-c(perm0.4n10000,my.importance[151:200,3])
 my.cf<-cforest(formula=y~.,data=my.data,control=cforest_control(ntree=500,maxdepth=2,replace=FALSE,mtry=round(sqrt(200))))
 my.importance<-varimp(my.cf)
 cforest0.05n10000<-c(cforest0.05n10000,my.importance[1:50])
 cforest0.1n10000<-c(cforest0.1n10000,my.importance[51:100])
 cforest0.25n10000<-c(cforest0.25n10000,my.importance[101:150])
 cforest0.4n10000<-c(cforest0.4n10000,my.importance[151:200])
 rm(my.importance)
 rm(my.cf)


save(gini0.05n10000,file="nullcase_n10000/gini0.05n10000.RData")
save(gini0.1n10000,file="nullcase_n10000/gini0.1n10000.RData")
save(gini0.25n10000,file="nullcase_n10000/gini0.25n10000.RData")
save(gini0.4n10000,file="nullcase_n10000/gini0.4n10000.RData")

save(perm0.05n10000,file="nullcase_n10000/perm0.05n10000.RData")
save(perm0.1n10000,file="nullcase_n10000/perm0.1n10000.RData")
save(perm0.25n10000,file="nullcase_n10000/perm0.25n10000.RData")
save(perm0.4n10000,file="nullcase_n10000/perm0.4n10000.RData")

save(cforest0.05n10000,file="nullcase_n10000/cforest0.05n10000.RData")
save(cforest0.1n10000,file="nullcase_n10000/cforest0.1n10000.RData")
save(cforest0.25n10000,file="nullcase_n10000/cforest0.25n10000.RData")
save(cforest0.4n10000,file="nullcase_n10000/cforest0.4n10000.RData")
     }