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

library(randomForest)
library(party)
source("fun.r")


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

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

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

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

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

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

niter<-100

for (i in 1:niter)
 {
 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),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.05powern10000<-c(gini0.05powern10000,my.importance[3,4])
 gini0.1powern10000<-c(gini0.1powern10000,my.importance[53,4])
 gini0.25powern10000<-c(gini0.25powern10000,my.importance[103,4])
 gini0.4powern10000<-c(gini0.4powern10000,my.importance[153,4])

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

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

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


 rm(my.rf)
 rm(my.importance)

save(gini0.05powern10000,file="powercase_n10000/gini0.05powern10000.RData")
save(gini0.1powern10000,file="powercase_n10000/gini0.1powern10000.RData")
save(gini0.25powern10000,file="powercase_n10000/gini0.25powern10000.RData")
save(gini0.4powern10000,file="powercase_n10000/gini0.4powern10000.RData")

save(gini0.05power1n10000,file="powercase_n10000/gini0.05power1n10000.RData")
save(gini0.1power1n10000,file="powercase_n10000/gini0.1power1n10000.RData")
save(gini0.25power1n10000,file="powercase_n10000/gini0.25power1n10000.RData")
save(gini0.4power1n10000,file="powercase_n10000/gini0.4power1n10000.RData")

save(gini0.05power2n10000,file="powercase_n10000/gini0.05power2n10000.RData")
save(gini0.1power2n10000,file="powercase_n10000/gini0.1power2n10000.RData")
save(gini0.25power2n10000,file="powercase_n10000/gini0.25power2n10000.RData")
save(gini0.4power2n10000,file="powercase_n10000/gini0.4power2n10000.RData")

save(perm0.05powern10000,file="powercase_n10000/perm0.05powern10000.RData")
save(perm0.1powern10000,file="powercase_n10000/perm0.1powern10000.RData")
save(perm0.25powern10000,file="powercase_n10000/perm0.25powern10000.RData")
save(perm0.4powern10000,file="powercase_n10000/perm0.4powern10000.RData")

save(perm0.05power1n10000,file="powercase_n10000/perm0.05power1n10000.RData")
save(perm0.1power1n10000,file="powercase_n10000/perm0.1power1n10000.RData")
save(perm0.25power1n10000,file="powercase_n10000/perm0.25power1n10000.RData")
save(perm0.4power1n10000,file="powercase_n10000/perm0.4power1n10000.RData")

save(perm0.05power2n10000,file="powercase_n10000/perm0.05power2n10000.RData")
save(perm0.1power2n10000,file="powercase_n10000/perm0.1power2n10000.RData")
save(perm0.25power2n10000,file="powercase_n10000/perm0.25power2n10000.RData")
save(perm0.4power2n10000,file="powercase_n10000/perm0.4power2n10000.RData")

}
