The data
We simulate a small data set with a covariance structure to explore the actions of RFlocalfdr. This package implements an empirical Bayes method of determining a significance level for the Random Forest MDI (Mean Decrease in Impurity) variable importance measure. See Dunne et al. (2022) for details
We
simulate the data set
fit a Random Forest model
use RFlocalfdr to estimate the significant variables
estimate the significant variables using
- Boruta (Kursa and Rudnicki (2010))
- Recursive Fearure Elimination
- AIR (Nembrini et al. (2018))
- PIMP (Altmann et al. (2010))
as comparators.
setup
We plot a small data set of \(300\times 508\) to explore the data generation method. This dataset consist of bands, with blocks of \(\{1, 2, 4, 8, 16, 32, 64\}\) of identical variables (see figure @ref(fig:simulation2). The variables are \(\in \{0,1,2\}\), a common encoding for genomic data where the numbers represent the number of copies of the minor allele. Only band 1 is used to calculate the \(y\) vector and \(y\) is 1 if any of \(X[, c(1, 2, 4, 8, 16, 32, 64)]\) is non-zero, so \(y\) is unbalanced, containing more 1’s than 0’s.
set.seed(13)
num_samples <- 300
num_bands <- 4
band_rank <- 6
num_vars <- num_bands * (2 ** (band_rank+1) -1)
print(num_vars)
X <- matrix(NA, num_samples , num_vars)
set.seed(123)
temp<-matrix(0,508,3)
var_index <- 1
for(band in 1:num_bands) {
for (rank in 0:band_rank) {
for (i in 1:2**rank) {
temp[var_index,]<-c(band , rank, var_index)
var_index <- var_index +1
}
}
}
#png("./supp_figures/small_simulated_data_set.png")
plot(temp[,1],ylim=c(0,10),type="p")
lines(temp[,2],type="b",col="red")
legend(0,10,c("band","rank"),pch=1,col=c(1,2))
#dev.off()
abline(v=c( 1, 2 , 4 , 8 ,16 ,32, 64 ))
table(temp[temp[,1]==1,2])
## # 0 1 2 3 4 5 6
## # 1 2 4 8 16 32 64 
A small simulated data set. Each band contains blocks of size {1, 2, 4, 8, 16, 32, 64}, and each block consists of correlated (identical variables).
data
We now generate the daa for the simulation. We have 50 bands and 300 observations so \(X\) is \(300 \times 6350\) with 127 non-null variables. We fit a standard RF to this dataset and record the resulting MDI importance score.
set.seed(13)
num_samples <- 300
num_bands <- 50
band_rank <- 6
num_vars <- num_bands * (2 ** (band_rank+1) -1)
print(num_vars)
X <- matrix(NA, num_samples , num_vars)
set.seed(123)
var_index <- 1
for(band in 1:num_bands) {
for (rank in 0:band_rank) {
# cat("rank=",rank,"\n")
var <- sample(0:2, num_samples, replace=TRUE)
for (i in 1:2**rank) {
X[,var_index] <- var
var_index <- var_index +1
}
# print(X[1:2,1:140])
# system("sleep 3")
}
}
y <- as.numeric(X[,1] > 1 | X[,2] > 1 | X[,4] > 1 | X[,8] > 1 | X[,16] > 1 |
X[,32] > 1 | X[,64] > 1 )