

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

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

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

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

for (i in 1:100)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,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.05<-c(gini0.05,my.importance[1:50,4])
 #gini0.1<-c(gini0.1,my.importance[51:100,4])
 #gini0.25<-c(gini0.25,my.importance[101:150,4])
 #gini0.4<-c(gini0.4,my.importance[151:200,4])
 #perm0.05<-c(perm0.05,my.importance[1:50,3])
 #perm0.1<-c(perm0.1,my.importance[51:100,3])
 #perm0.25<-c(perm0.25,my.importance[101:150,3])
 #perm0.4<-c(perm0.4,my.importance[151:200,3])
 my.cf<-cforest(formula=y~.,data=my.data,control=cforest_control(teststat = "quad",testtype = "Univ",mincriterion = 0, minbucket=7, maxdepth=2,replace=FALSE,ntree=500,mtry=round(sqrt(200))))
 my.importance<-varimp(my.cf)
 cforest0.05<-c(cforest0.05,my.importance[1:50])
 cforest0.1<-c(cforest0.1,my.importance[51:100])
 cforest0.25<-c(cforest0.25,my.importance[101:150])
 cforest0.4<-c(cforest0.4,my.importance[151:200])
 rm(my.importance)
 rm(my.cf)
 
 #save(gini0.05,file="nullcase/gini0.05.RData")
 #save(gini0.1,file="nullcase/gini0.1.RData")
 #save(gini0.25,file="nullcase/gini0.25.RData")
 #save(gini0.4,file="nullcase/gini0.4.RData")

 #save(perm0.05,file="nullcase/perm0.05.RData")
 #save(perm0.1,file="nullcase/perm0.1.RData")
 #save(perm0.25,file="nullcase/perm0.25.RData")
 #save(perm0.4,file="nullcase/perm0.4.RData")
 
 save(cforest0.05,file="nullcase/cforest0.05.RData")
 save(cforest0.1,file="nullcase/cforest0.1.RData")
 save(cforest0.25,file="nullcase/cforest0.25.RData")
 save(cforest0.4,file="nullcase/cforest0.4.RData")
 }



