
####################################################
# Bias in variable importance measure: power study #
####################################################

library(randomForest)
library(party)
source("fun.r")
#gini0.05power<-c()
#gini0.1power<-c()
#gini0.25power<-c()
#gini0.4power<-c()

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

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

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

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

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

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

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

cforest0.05power2<-c()
cforest0.1power2<-c()
cforest0.25power2<-c()
cforest0.4power2<-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),beta=rep(c(c(log(3),log(1.5)),rep(0,48)),4))
 #my.rf<-randomForest(x=my.data[,-1],y=my.data[,1],importance=TRUE)
 #my.importance<-importance(my.rf,scale=FALSE)

 #gini0.05power<-c(gini0.05power,my.importance[3,4])
 #gini0.1power<-c(gini0.1power,my.importance[53,4])
 #gini0.25power<-c(gini0.25power,my.importance[103,4])
 #gini0.4power<-c(gini0.4power,my.importance[153,4])

 #gini0.05power1<-c(gini0.05power1,my.importance[1,4])
 #gini0.1power1<-c(gini0.1power1,my.importance[51,4])
 #gini0.25power1<-c(gini0.25power1,my.importance[101,4])
 #gini0.4power1<-c(gini0.4power1,my.importance[151,4])
 
 #gini0.05power2<-c(gini0.05power2,my.importance[2,4])
 #gini0.1power2<-c(gini0.1power2,my.importance[52,4])
 #gini0.25power2<-c(gini0.25power2,my.importance[102,4])
 #gini0.4power2<-c(gini0.4power2,my.importance[152,4])

 #perm0.05power<-c(perm0.05power,my.importance[3,3])
 #perm0.1power<-c(perm0.1power,my.importance[53,3])
 #perm0.25power<-c(perm0.25power,my.importance[103,3])
 #perm0.4power<-c(perm0.4power,my.importance[153,3])
 
 #perm0.05power1<-c(perm0.05power1,my.importance[1,3])
 #perm0.1power1<-c(perm0.1power1,my.importance[51,3])
 #perm0.25power1<-c(perm0.25power1,my.importance[101,3])
 #perm0.4power1<-c(perm0.4power1,my.importance[151,3])

 #perm0.05power2<-c(perm0.05power2,my.importance[2,3])
 #perm0.1power2<-c(perm0.1power2,my.importance[52,3])
 #perm0.25power2<-c(perm0.25power2,my.importance[102,3])
 #perm0.4power2<-c(perm0.4power2,my.importance[152,3])

 my.cf<-cforest(formula=y~.,data=my.data,control=cforest_control(teststat = "quad",testtype = "Univ",mincriterion = 0, minbucket=7, replace=FALSE,ntree=500,mtry=round(sqrt(200))))
 my.importance<-varimp(my.cf)
 cforest0.05power<-c(cforest0.05power,my.importance[3])
 cforest0.1power<-c(cforest0.1power,my.importance[53])
 cforest0.25power<-c(cforest0.25power,my.importance[103])
 cforest0.4power<-c(cforest0.4power,my.importance[153])

 cforest0.05power1<-c(cforest0.05power1,my.importance[1])
 cforest0.1power1<-c(cforest0.1power1,my.importance[51])
 cforest0.25power1<-c(cforest0.25power1,my.importance[101])
 cforest0.4power1<-c(cforest0.4power1,my.importance[151])
 
 cforest0.05power2<-c(cforest0.05power2,my.importance[2])
 cforest0.1power2<-c(cforest0.1power2,my.importance[52])
 cforest0.25power2<-c(cforest0.25power2,my.importance[102])
 cforest0.4power2<-c(cforest0.4power2,my.importance[152])

 rm(my.importance)
 rm(my.cf)
 
 

save(cforest0.05power,file="powercase/cforest0.05power.RData")
save(cforest0.1power,file="powercase/cforest0.1power.RData")
save(cforest0.25power,file="powercase/cforest0.25power.RData")
save(cforest0.4power,file="powercase/cforest0.4power.RData")

save(cforest0.05power1,file="powercase/cforest0.05power1.RData")
save(cforest0.1power1,file="powercase/cforest0.1power1.RData")
save(cforest0.25power1,file="powercase/cforest0.25power1.RData")
save(cforest0.4power1,file="powercase/cforest0.4power1.RData")

save(cforest0.05power2,file="powercase/cforest0.05power2.RData")
save(cforest0.1power2,file="powercase/cforest0.1power2.RData")
save(cforest0.25power2,file="powercase/cforest0.25power2.RData")
save(cforest0.4power2,file="powercase/cforest0.4power2.RData")


 }



#save(gini0.05power,file="powercase/gini0.05power.RData")
#save(gini0.1power,file="powercase/gini0.1power.RData")
#save(gini0.25power,file="powercase/gini0.25power.RData")
#save(gini0.4power,file="powercase/gini0.4power.RData")

#save(gini0.05power1,file="powercase/gini0.05power1.RData")
#save(gini0.1power1,file="powercase/gini0.1power1.RData")
#save(gini0.25power1,file="powercase/gini0.25power1.RData")
#save(gini0.4power1,file="powercase/gini0.4power1.RData")

#save(gini0.05power2,file="powercase/gini0.05power2.RData")
#save(gini0.1power2,file="powercase/gini0.1power2.RData")
#save(gini0.25power2,file="powercase/gini0.25power2.RData")
#save(gini0.4power2,file="powercase/gini0.4power2.RData")

#save(perm0.05power,file="powercase/perm0.05power.RData")
#save(perm0.1power,file="powercase/perm0.1power.RData")
#save(perm0.25power,file="powercase/perm0.25power.RData")
#save(perm0.4power,file="powercase/perm0.4power.RData")

#save(perm0.05power1,file="powercase/perm0.05power1.RData")
#save(perm0.1power1,file="powercase/perm0.1power1.RData")
#save(perm0.25power1,file="powercase/perm0.25power1.RData")
#save(perm0.4power1,file="powercase/perm0.4power1.RData")

#save(perm0.05power2,file="powercase/perm0.05power2.RData")
#save(perm0.1power2,file="powercase/perm0.1power2.RData")
#save(perm0.25power2,file="powercase/perm0.25power2.RData")
#save(perm0.4power2,file="powercase/perm0.4power2.RData")

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



