--- title: "binomialRF Feature Selection Vignette" author: "Samir Rachid Zaim" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{"binomialRF Feature Selection Vignette"} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} library('randomForest') library('data.table') library('stats') library('binomialRF') knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` The $\textit{binomialRF}$ is a $\textit{randomForest}$ feature selection wrapper (Zaim 2019) that treats the random forest as a binomial process where each tree represents an iid bernoulli random variable for the event of selecting $X_j$ as the main splitting variable at a given tree. The algorithm below describes the technical aspects of the algorithm. \begin{algorithm} \caption{binomialRF Feature Selection Algorithm } \label{alg1} \begin{algorithmic} \For{i=1:N} \State Grow $T_{i}$ \State $S_{ij} = \left\{ \begin{array}{ll} 1 & X_j \text{ is splitting variable at root node by } T_i \\ 0 & \text{otherwise} \\ \end{array} \right. $ \State $S \gets S_{ij}$ \EndFor \State $S_j = \sum_{i=1}^N S_{ij}= \sum_{i=1}^N I[X_j \in root(T_i)]$ \State $S_j \sim \text{binomial}(N,p_0), \hspace{6mm} \text{where } p_0=\bigg(1 - \prod_{i=1}^m\frac{L-i}{L-(i-1)} \bigg)\frac{1}{m}$. \State Test for Significance \State Adjust for Multiple Comparisons \end{algorithmic} \end{algorithm} ![binomialRF Algorithm.](algorithm.png) # Simulating Data Since $\textit{binomialRF}$ is a wrapper algorithm that internally calls and grows a randomForest object based on the inputted parameters. First we generate a simple simulated logistic data as follows: * $X_{10}\sim MNV(0, I_{10})$, * $p(x) = \frac{1}{1+e^{-X\beta}}$, and * $y \sim Binom(10,p)$. where $\beta$ is a vector of coefficients where the first 2 coefficients are set to 3, and the rest are 0. $$\beta = \begin{bmatrix} 3 & 3 & 0 & \cdots & 0 \end{bmatrix}^T$$ ## Simulated Data ```{r echo=T, warning=F, message=F} set.seed(324) ### Generate multivariate normal data in R10 X = matrix(rnorm(1000), ncol=10) ### let half of the coefficients be 0, the other be 10 trueBeta= c(rep(3,2), rep(0,8)) ### do logistic transform and generate the labels z = 1 + X %*% trueBeta pr = 1/(1+exp(-z)) y = rbinom(100,1,pr) ``` To generate data looking like this: ```{r, echo=FALSE, results='asis'} knitr::kable(head(cbind(round(X,2),y), 10)) ``` ## Generating the Stable Correlated Binomial Distribution Since the binomialRF requires a correlation adjustment to adjust for the tree-to-tree sampling correlation, we first generate the appropriately-parameterized stable correlated binomial distribution. Note, the correlbinom function call can take a while to execute for large number trials (i.e., trials > 1000). ```{r echo=T, warning=F, message=F} require(correlbinom) rho = 0.33 ntrees = 250 cbinom = correlbinom(rho, successprob = 1/ncol(X), trials = ntrees, precision = 1024, model = 'kuk') ``` ## binomialRF Function Call Then we can run the binomialRF function call as below: ```{r echo=T, warning=F, message=F} binom.rf <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = .6, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) print(binom.rf) ``` # Tuning Parameters ## Percent_features Note that since the binomial exact test is contingent on a test statistic measuring the likelihood of selecting a feature, if there is a dominant feature, then it will render all remaining 'important' features useless as it will always be selected as the splitting variable. So it is important to set the $percent_features$ parameter < 1. The results below show how setting the parameter to a fraction between .6 to 1 can allow other features to stand out as important. ```{r echo=F, warning=F, message=F} # set.seed(324) binom.rf <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = 1, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) cat('\n\nbinomialRF 100%\n\n') print(binom.rf) binom.rf <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = .8, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) cat('\n\nbinomialRF 80%\n\n') print(binom.rf) binom.rf <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = .6, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) cat('\n\nbinomialRF 60%\n\n') print(binom.rf) ``` ## ntrees We recommend growing at least 500 to 1,000 trees at a minimum so that the algorithm has a chance to stabilize, but also recommend choosing ntrees as a function of the number of features in your dataset. The ntrees tuning parameter must be set in conjunction with the percent_features as these two are inter-connectedm as well as the number of true features in the model. Since the correlbinom function call is slow to execute for ntrees > 1000, we recommend growing random forests with only 500-1000 trees. ```{r echo=F, warning=F, message=F} set.seed(324) binom.rf1000 <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = .5, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) rho = 0.33 ntrees = 500 cbinom = correlbinom(rho, successprob = 1/ncol(X), trials = ntrees, precision = 1024, model = 'kuk') binom.rf500 <- binomialRF::binomialRF(X,factor(y), fdr.threshold = .05, ntrees = ntrees,percent_features = .5, fdr.method = 'BY', user_cbinom_dist = cbinom, sampsize = round(nrow(X)*.33)) cat('\n\nbinomialRF 250 trees\n\n') print(binom.rf500) cat('\n\nbinomialRF 500 trees \n\n') print(binom.rf1000) ```