#setwd("Z:/Manuskripte/Papers/Preparation/Gini/Analysen_Final/Analysen_chi2/FirstLayer")

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

niter<-100

#########################
####  replace=TRUE ####
#########################

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

for (i in 1:niter)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,beta0=0,maf=rep(c(0.05,0.1,0.25,0.4),each=50))
 my.<-randomForest(x=my.data[,-1],y=my.data[,1],replace=TRUE,maxnodes=2)
 rfvar<-my.$forest$bestvar[1,]
 rf0.05[i]<-sum(is.element(rfvar,1:50))
 rf0.1[i]<-sum(is.element(rfvar,51:100))
 rf0.25[i]<-sum(is.element(rfvar,101:150))
 rf0.4[i]<-sum(is.element(rfvar,151:200))
 
 save(rf0.05,file="nullcase/rf0.05.RData")
 save(rf0.1,file="nullcase/rf0.1.RData")
 save(rf0.25,file="nullcase/rf0.25.RData")
 save(rf0.4,file="nullcase/rf0.4.RData")

 }


 
#######

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

niter<-100

for (i in 1:niter)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,beta0=0,maf=rep(c(0.05,0.1,0.25,0.4),each=50))
 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))))                    
 freqsel<-tabulate(unlist(lapply(my.cf@ensemble,function(x)x[[5]][[1]])),nbins=200)
 cforest0.05[i]<-sum(freqsel[1:50])
 cforest0.1[i]<-sum(freqsel[51:100])
 cforest0.25[i]<-sum(freqsel[101:150])
 cforest0.4[i]<-sum(freqsel[151:200])
 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")
 }


#### n=10000 ###
 
rf0.05n10000<-c()
rf0.1n10000<-c()
rf0.25n10000<-c()
rf0.4n10000<-c()

for (i in 1:niter)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,n.obs=10000,beta0=0,maf=rep(c(0.05,0.1,0.25,0.4),each=50))
 my.<-randomForest(x=my.data[,-1],y=my.data[,1],replace=TRUE,maxnodes=2)
 rfvar<-my.$forest$bestvar[1,]
 rf0.05n10000[i]<-sum(is.element(rfvar,1:50))
 rf0.1n10000[i]<-sum(is.element(rfvar,51:100))
 rf0.25n10000[i]<-sum(is.element(rfvar,101:150))
 rf0.4n10000[i]<-sum(is.element(rfvar,151:200))
 
 save(rf0.05n10000,file="nullcase_n10000/rf0.05n10000.RData")
 save(rf0.1n10000,file="nullcase_n10000/rf0.1n10000.RData")
 save(rf0.25n10000,file="nullcase_n10000/rf0.25n10000.RData")
 save(rf0.4n10000,file="nullcase_n10000/rf0.4n10000.RData")

 }


 
#######

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

for (i in 1:niter)
 {
 set.seed(i)
 print(i)
 my.data<-simSNP(n.snp=200,n.obs=10000,beta0=0,maf=rep(c(0.05,0.1,0.25,0.4),each=50))
 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))))                    
 freqsel<-tabulate(unlist(lapply(my.cf@ensemble,function(x)x[[5]][[1]])),nbins=200)
 cforest0.05n10000[i]<-sum(freqsel[1:50])
 cforest0.1n10000[i]<-sum(freqsel[51:100])
 cforest0.25n10000[i]<-sum(freqsel[101:150])
 cforest0.4n10000[i]<-sum(freqsel[151:200])
 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")
 }
