% \VignetteIndexEntry{caret Manual -- Variable Selection}
% \VignetteDepends{caret}
% \VignettePackage{caret}
\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage[pdftex]{graphicx}
\usepackage{color}
\usepackage{xspace}
\usepackage{fancyvrb}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage[lined,algonl,algoruled, linesnumbered]{algorithm2e}
\usepackage[
         colorlinks=true,
         linkcolor=blue,
         citecolor=blue,
         urlcolor=blue]
         {hyperref}
\usepackage{lscape}
\usepackage{ctable}
\usepackage{Sweave}
\SweaveOpts{keep.source=TRUE}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define new colors for use
\definecolor{darkgreen}{rgb}{0,0.6,0}
\definecolor{darkred}{rgb}{0.6,0.0,0}
\definecolor{lightbrown}{rgb}{1,0.9,0.8}
\definecolor{brown}{rgb}{0.6,0.3,0.3}
\definecolor{darkblue}{rgb}{0,0,0.8}
\definecolor{darkmagenta}{rgb}{0.5,0,0.5}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\bld}[1]{\mbox{\boldmath $#1$}}
\newcommand{\shell}[1]{\mbox{$#1$}}
\renewcommand{\vec}[1]{\mbox{\bf {#1}}}

\newcommand{\ReallySmallSpacing}{\renewcommand{\baselinestretch}{.6}\Large\normalsize}
\newcommand{\SmallSpacing}{\renewcommand{\baselinestretch}{1.1}\Large\normalsize}

\newcommand{\halfs}{\frac{1}{2}}

\setlength{\oddsidemargin}{-.25 truein}
\setlength{\evensidemargin}{0truein}
\setlength{\topmargin}{-0.2truein}
\setlength{\textwidth}{7 truein}
\setlength{\textheight}{8.5 truein}
\setlength{\parindent}{0.20truein}
\setlength{\parskip}{0.10truein}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\pagestyle{fancy}
\lhead{}
\chead{Variable Selection Using The {\tt caret} Package}
\rhead{}
\lfoot{}
\cfoot{}
\rfoot{\thepage\ of \pageref{LastPage}}
\renewcommand{\headrulewidth}{1pt}
\renewcommand{\footrulewidth}{1pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Variable Selection Using The {\tt caret} Package}
\author{Max Kuhn \\ max.kuhn@pfizer.com}

\begin{document}

\maketitle

\thispagestyle{empty}

\section{Models with Built--In Feature Selection}

Many models that can be accessed using \texttt{caret}'s \texttt{train} function produce prediction equations that do not necessarily use all the predictors. These models are thought to have built--in feature selection and include \texttt{rpart}, \texttt{gbm}, \texttt{ada}, \texttt{glmboost}, \texttt{gamboost}, \texttt{blackboost}, \texttt{ctree}, \texttt{sparseLDA}, \texttt{sddaLDA}, \texttt{sddaQDA}	\texttt{glmnet}, \texttt{lasso}, \texttt{lars}, \texttt{spls}, \texttt{earth}, \texttt{fda},   \texttt{bagEarth}, \texttt{bagFDA} , \texttt{pam} and others. Many of the functions have an ancillary method called \texttt{predictors} that returns a vector indicating which predictors were used in the final model.

In many cases, using these models with built--in feature selection will be more efficient than algorithms where the search routine for the right predictors is external to the model (see Section 2). Built--in feature selection typically couples the predictor search algorithm with the parameter  estimation and are usually optimized with a single objective function (e.g. error rates or likelihood). 


\section{Feature Selection Using Search Algorithms}

\subsection{Searching the Feature Space}

Many feature selection routines used a ``wrapper'' approach to find appropriate variables such that an algorithm that searches the feature space repeatedly fits the model with different predictor sets. The best predictor set is determined by some measure of performance (i.e. $R^2$, classification accuracy, etc). Examples of search functions are genetic algorithms, simulated annealing and forward/backward/stepwise selection methods. In theory, each of these search routines could converge to an optimal set of predictors.

An example of one search routine is backwards selection (a.k.a. recursive feature elimination). 

\subsubsection{Backwards Selection}

First, the algorithm fits the model to all predictors. Each predictor is ranked using it's importance to the model. Let $S$ be a sequence of ordered numbers which are candidate values for the number of predictors to retain ($S_1 > S_2, \ldots$). At each iteration of feature selection, the $S_i$ top raked predictors are retained, the model is refit and performance is assessed. The value of $S_i$ with the best performance is determined and the top $S_i$ predictors are used to fit the final model. Algorithm \ref{A:rfe} has a more complete definition.

The algorithm has an optional step (line \ref{A:rerank}) where the predictor rankings are recomputed on the model on the reduced feature set. Svetnik {\it el al} (2004) showed that, for random forest models, there was a decrease in performance when the rankings were re--computed at every step. However, in other cases when the initial rankings are not good (e.g. linear models with highly collinear predictors), re--calculation can slightly improve performance.


\begin{algorithm}
   \caption{Recursive feature elimination}
   \label{A:rfe}
   \SetAlgoLined
   \RestyleAlgo{plain}
   \DontPrintSemicolon


    \vspace*{3pt} Tune/train the model on the training set using all predictors\vspace*{3pt}\; 

    \vspace*{3pt} Calculate model performance\vspace*{3pt}\; 

     \vspace*{3pt} Calculate variable importance or rankings\vspace*{3pt} \; 

     \For{Each subset size $S_i$, $i=1\ldots S$}{

         \vspace*{3pt} Keep the $S_i$ most important variables \vspace*{3pt}\;
         
         \vspace*{3pt} [Optional] Pre--process the data\;
			
         \vspace*{3pt} Tune/train the model on the training set using $S_i$ predictors \nllabel{A:resample}\vspace*{3pt}\;

         \vspace*{3pt}	Calculate model performance \vspace*{3pt}\;

         \vspace*{3pt} [Optional] Recalculate the rankings for each predictor\vspace*{3pt}\nllabel{A:rerank}\;
	  				
      } % end resample

    \vspace*{3pt} Calculate the performance profile over the $S_i$ \vspace*{3pt}\;
    
    \vspace*{3pt} Determine the appropriate number of predictors\vspace*{3pt}\;

    \vspace*{3pt} Determine the final ranks of each predictor\vspace*{3pt}\;

    \vspace*{3pt} Fit the final model based on the optimal $S_i$ \vspace*{3pt}\;    
    
\end{algorithm}

One potential issue over--fitting to the predictor set such that the wrapper procedure could focus on nuances of the training data that are not found in future samples (i.e. over--fitting to predictors and samples). 

For example, suppose a very large number of uninformative predictors were collected and one such predictor randomly correlated with the outcome. The RFE algorithm would give a good rank to this variable and the prediction error (on the same data set) would be lowered. It would take a different test/validation to find out that this predictor was uninformative. The was referred to as ``selection bias'' by Ambroise and McLachlan (2002).

In the current RFE algorithm, the training data is being used for at least three purposes: predictor selection, model fitting and performance evaluation. Unless the number of samples is large, especially in relation to the number of variables, one static training set may not be able to fulfill these needs.

\subsection{Resampling and External Validation}

Since feature selection is part of the model building process, resampling methods (e.g. cross--validation, the bootstrap) should factor in the variability caused by feature selection when calculating performance. For example, the RFE procedure in Algorithm \ref{A:rfe} can estimate the model performance on line \ref{A:resample}, which during the selection process.  Ambroise and McLachlan (2002) and Svetnik {\it el al} (2004) showed that improper use of resampling to measure performance will result in models that perform poorly on new samples.

To get performance estimates that incorporate the variation due to feature selection, it is suggested that the steps in Algorithm \ref{A:rfe} be encapsulated inside an outer layer of resampling (e.g. 10--fold cross--validation). Algorithm \ref{A:rfe2} shows a version of the algorithm that uses resampling.

While this will provide better estimates of performance, it is more computationally burdensome. For users with access to machines with multiple processors, the first \texttt{For} loop in Algorithm \ref{A:rfe2} (line \ref{A:outer}) can be easily parallelized. Another complication to using resampling is that multiple lists of the ``best'' predictors are generated at each iteration. At first this may seem like a disadvantage, but it does provide a more probabilistic assessment of predictor importance than a ranking based on a single fixed data set. At the end of the algorithm, a consensus ranking can be used to determine the best predictors to retain.

\begin{algorithm}
   \caption{Recursive feature elimination incorporating resampling}
   \label{A:rfe2}
   \SetAlgoLined
   \RestyleAlgo{plain}
   \DontPrintSemicolon

   \For{Each Resampling Iteration\nllabel{A:outer}}{
	
     \vspace*{3pt} Partition data into training and test/hold--back set via resampling\vspace*{3pt} \; 

    \vspace*{3pt} Tune/train the model on the training set using all predictors\nllabel{A:fit0}\vspace*{3pt}\; 

    \vspace*{3pt} Predict the held--back samples\nllabel{A:pred0}\vspace*{3pt}\; 

     \vspace*{3pt} Calculate variable importance or rankings\nllabel{A:rank0}\vspace*{3pt} \; 

     \For{Each subset size $S_i$, $i=1\ldots S$}{

         \vspace*{3pt} Keep the $S_i$ most important variables \vspace*{3pt}\;	
         
         \vspace*{3pt} [Optional] Pre--process the data\;
			
         \vspace*{3pt} Tune/train the model on the training set using $S_i$ predictors\nllabel{A:fit1}\vspace*{3pt}\;

         \vspace*{3pt}	Predict the held--back samples \nllabel{A:pred1}\vspace*{3pt}\;

         \vspace*{3pt} [Optional] Recalculate the rankings for each predictor\nllabel{A:rank1}\vspace*{3pt}\;
	  				
      } % end resample
    }
    \vspace*{3pt} Calculate the performance profile over the $S_i$ using the held--back samples \nllabel{A:summary}\vspace*{3pt}\;
    
    \vspace*{3pt} Determine the appropriate number of predictors \nllabel{A:selectSize}\vspace*{3pt}\;
    
   \vspace*{3pt}  Estimate the final list of predictors to keep in the final model\nllabel{A:selectVar}\vspace*{3pt}\;

    \vspace*{3pt} Fit the final model based on the optimal $S_i$ using the original training set\nllabel{A:fit2} \vspace*{3pt}\;    
    
\end{algorithm}

\clearpage
\section{Recursive Feature Elimination via \texttt{caret}}

In \texttt{caret}, Algorithm \ref{A:rfe} is implemented by the function \texttt{rfeIter}. The resampling-based Algorithm \ref{A:rfe2} is in the \texttt{rfe} function. Given the potential selection bias issues, this document focuses on \texttt{rfe}.  There are several arguments:
\begin{itemize}
\item \texttt{x}, a matrix or data frame of predictor variables
\item \texttt{y}, a vector (numeric or factor) of outcomes
\item \texttt{sizes}, a integer vector for the specific subset sizes that should be tested (which need not to include \texttt{ncol(x)}) 
\item \texttt{rfeControl}, a list of options that can be used to specify the model and the methods for prediction, ranking etc.
\end{itemize}
For a specific model, a set of functions must be specified in \texttt{rfeControl$\$$functions}. Section \ref{S:helpers} below has descriptions of these sub--functions. There are a number of pre--defined sets of functions for several models, including: linear regression (in the object \texttt{lmFuncs}), random forests (\texttt{rfFuncs}), naive Bayes (\texttt{nbFuncs}), bagged trees (\texttt{treebagFuncs}) and functions that can be used with \texttt{caret}'s \texttt{train} function (\texttt{caretFuncs}). The latter is useful if the model has tuning parameters that must be determined at each iteration.


\subsection{An Example}

<<startUp, echo=FALSE>>=
library(caret)
library(mlbench)
library(Hmisc)
library(randomForest)

@ 

To test the algorithm, the ``Friedman 1'' benchmark (Friedman, 1991) was used. There are five informative variables generated by the equation
\[
y = 10 \sin(\pi x_1x_2) + 20 (x_3 - 0.5)^2 + 10 x_4 + 5 x_5 + N(0, \sigma^2)
\]
In the simulation used here:
<<simSettings>>=
n <- 100
p <- 40
sigma <- 1
set.seed(1)
sim <- mlbench.friedman1(n, sd = sigma)
colnames(sim$x) <- c(paste("real", 1:5, sep = ""),
                     paste("bogus", 1:5, sep = ""))
bogus <- matrix(rnorm(n * p), nrow = n)
colnames(bogus) <- paste("bogus", 5+(1:ncol(bogus)), sep = "")
x <- cbind(sim$x, bogus)
y <- sim$y

@ 
Of the \Sexpr{ncol(x)} predictors, there are \Sexpr{ncol(x)-5} pure noise variables: 5 are uniform on [0, 1] and \Sexpr{p} are random univariate standard normals. The predictors are centered and scaled:
<<preProc>>=
normalization <- preProcess(x)
x <- predict(normalization, x)
x <- as.data.frame(x)
subsets <- c(1:5, 10, 15, 20, 25)
@ 
The simulation will fit models with subset sizes of \Sexpr{paste(sort(subsets, decreasing = TRUE), collapse = ", ")}. 


As previously mentioned, to fit linear models, the \texttt{lmFuncs} set of functions can be used. To do this, a control object is created with the \texttt{rfeControl} function. We also specify that 10--fold cross--validation should be used in line \ref{A:outer} of Algorithm  \ref{A:rfe2}. The number of folds can be changed via the \texttt{number} argument to \texttt{rfeControl} (defaults to 10). The \texttt{verbose} option prevents copious amounts of output from being produced and the \texttt{returnResamp} argument specifies that the 10 performance estimates should be saved only for the optimal subset size.

<<lm>>=
set.seed(10)

ctrl <- rfeControl(functions = lmFuncs,
                   method = "cv",
                   verbose = FALSE,
                   returnResamp = "final")

lmProfile <- rfe(x, y,
                 sizes = subsets,
                 rfeControl = ctrl)

lmProfile
@ 
The output shows that the  best subset size was estimated to be \Sexpr{lmProfile$bestSubset} predictors. This set includes informative variables but did not include them all. The \texttt{predictors} function can be used to get a text string of variable names that were picked in the final model. The \texttt{lmProfile} is a list of class \texttt{"rfe"} that contains an object \texttt{fit} that is the final linear model with the remaining terms. The model can be used to get predictions for future or test samples. 

<<lmObjects>>=
predictors(lmProfile)
lmProfile$fit
lmProfile$resample
@ 

There are also several plot methods to visualize the results. \texttt{plot(lmProfile)} produces the performance profile across different subset sizes, as shown in Figure \ref{F:lmProfile}. Also the resampling results are stored in the sub--object \texttt{lmProfile$\$$resample} and can be used with several lattice functions. Univariate lattice functions (\texttt{densityplot}, \texttt{histogram}) can be used to plot the resampling distribution while bivariate functions (\texttt{xyplot}, \texttt{stripplot}) can be used to plot the distributions for different subset sizes. In the latter case, the option \texttt{returnResamp = "all"} in \texttt{rfeControl} can be used to save all the resampling results. See Figure \ref{F:rf2} for two examples.

\begin{figure}
   \begin{center}		
<<lmPlot, echo = FALSE, results = hide>>=
pdf("lm.pdf", width = 5, height = 7)
   trellis.par.set(caretTheme())
   plot1 <- plot(lmProfile, type = c("g", "o"))
   plot2 <- plot(lmProfile, type = c("g", "o"), metric = "Rsquared")
   print(update(plot1, 
                ylab = "Resampled RMSE"), 
         split=c(1,1,1,2), 
         more=TRUE)
   print(update(plot2, 
                ylab = "Resampled R^2"), 
         split=c(1,2,1,2))
dev.off()
@
      \includegraphics[clip]{lm}
      \caption{Performance profiles for recursive feature elimination using linear models. These images were generated by \texttt{plot(lmProfile)} and  \texttt{plot(lmProfile, metric = "Rsquared")}.}
      \label{F:lmProfile} 
    \end{center}
\end{figure} 
 

\subsection{Helper Functions}\label{S:helpers}

To use feature elimination for an arbitrary model, a set of functions
must be passed to \texttt{rfe} for each of the steps in Algorithm
\ref{A:rfe2}. 

This section defines those functions and uses the
existing random forest functions as an illustrative
example. \texttt{caret} contains a list called \texttt{rfFuncs}, but
this document will use a more simple version that will be better for
illustrating the ideas at play. The functions used here are collected
in a list called \texttt{rfRFE}.

<<rfeAlt, results = hide, echo = FALSE>>=
rfRFE <-  list(summary = defaultSummary,
                 fit = function(x, y, first, last, ...)
                 {
                   library(randomForest)
                   randomForest(x, y, importance = first, ...)
                 },
                 pred = function(object, x)
                 {
                   predict(object, x)
                 },
                 rank = function(object, x, y)
                 {
                   vimp <- varImp(object)

                   if(is.factor(y))
                     {
                       if(all(levels(y) %in% colnames(vimp)))
                         {
                           avImp <- apply(vimp[, levels(y), drop = TRUE],
                                          1,
                                          mean)
                           vimp$Overall <- avImp
                         }

                     }
                   
                   vimp <- vimp[
                                order(
                                      vimp$Overall,
                                      decreasing = TRUE)
                                ,,
                                drop = FALSE]
                   
                   vimp$var <- rownames(vimp)                  
                   vimp
                 },
                 selectSize = pickSizeBest,
                 selectVar = pickVars)
@ 


\subsubsection{The \texttt{summary} Function}\label{S:summary}

The \texttt{summary} function takes the observed and predicted values
and computes one or more performance metrics (see line
\ref{A:summary}). The input is a data frame with columns \texttt{obs}
and \texttt{pred}. The output should be a named vector of numeric
variables. Note that the \texttt{metric} argument of the \texttt{rfe}
function should reference one of the names of the output of
\texttt{summary}. The example function is:
<<rfSummary>>=
rfRFE$summary
@ 

Two functions in \texttt{caret} that can be used as the summary
funciton are \texttt{defaultSummary} and \texttt{twoClassSummary} (for
classification probelms with two classes). 


\subsubsection{The \texttt{fit} Function}\label{S:fit}

This function builds the model based on the current data set (lines\ref{A:fit0},  \ref{A:fit1} and \ref{A:fit2}). The arguments for the function must be:
\begin{itemize}
\item \texttt{x}: the current training set of predictor data with the appropriate subset of variables
\item \texttt{y}: the current outcome data (either a numeric or factor vector)
\item \texttt{first}: a single logical value for whether the current predictor set has all possible variables (e.g. line \ref{A:fit0})
  \item \texttt{last}: similar to \texttt{first}, but \texttt{TRUE} when the last model is fit with the final subset size and predictors. (line \ref{A:fit2}) 
\item $\ldots$: optional arguments to pass to the fit function in the call to \texttt{rfe}
\end{itemize}
The function should return a model object that can be used to generate predictions. For random forest, the fit function is simple:
<<rfModel>>=
rfRFE$fit
@ 
For feature selection without re--ranking at each iteration, the random forest variable importances only need to be computed on the first iterations when all of the predictors are in the model. This can be accomplished using  \texttt{importance = first}.

\subsubsection{The \texttt{pred} Function}\label{S:pred}

This function returns a vector of predictions (numeric or factors) from the current model (lines \ref{A:pred0} and \ref{A:pred1}). The input arguments must be
\begin{itemize}
\item \texttt{object}: the model generated by the \texttt{fit} function
\item \texttt{x}: the current set of predictor set for the held--back samples
\end{itemize}
For random forests, the function is a simple wrapper for the predict function:
<<rfPredict>>=
rfRFE$pred
@ 
For classification, it is probably a good idea to ensure that the resulting factor variables of predictions has the same levels as the input data.

\subsubsection{The \texttt{rank} Function}

This function is used to return the predictors in the order of the most important to the least important (lines \ref{A:rank0} and \ref{A:rank1}). Inputs are:
\begin{itemize}
\item \texttt{object}: the model generated by the \texttt{fit} function
\item \texttt{x}: the current set of predictor set for the training samples
\item \texttt{y}: the current training outcomes
\end{itemize}
The function should return a data frame with a column called \texttt{var} that has the current variable names. The first row should be the most important predictor etc. Other columns can be included in the output and will be returned in the final \texttt{rfe} object.

For random forests, the function below uses \texttt{caret}'s \texttt{varImp} function to extract the random forest importances and orders them. For classification, \texttt{randomForest} will produce a column of importances for each class. In this case, the default ranking function orders the predictors by the averages importance across the classes.
<<rfRank>>=
rfRFE$rank
@ 

\subsubsection{The \texttt{selectSize} Function}

This function determines the optimal number of predictors based on the resampling output (line \ref{A:selectSize}). Inputs for the function are:
\begin{itemize}
\item \texttt{x}: a matrix with columns for the performance metrics and the number of variables, called \texttt{Variables}
\item \texttt{metric}: a character string of the performance measure to optimize (e.g. RMSE, Accuracy)
\item \texttt{maximize}: a single logical for whether the metric should be maximized
\end{itemize}
This function should return an integer corresponding to the optimal subset size.

\texttt{caret} comes with two examples functions for this purpose: \texttt{pickSizeBest} and \texttt{pickSizeTolerance}. The former simply selects the subset size that has the best value. The latter takes into account the whole profile and tries to pick a subset size that is small without sacrificing too much performance. For example, suppose we have computed the RMSE over a series of variables sizes:
<<tolerance>>=
  example <- data.frame(RMSE = c(
                          3.215, 2.819, 2.414, 2.144, 
                          2.014, 1.997, 2.025, 1.987, 
                          1.971, 2.055, 1.935, 1.999, 
                          2.047, 2.002, 1.895, 2.018),
                        Variables = 1:16)
example
@ 
These are depicted in Figure \ref{F:tolerance}. The solid circle identifies the subset size with the absolute smallest RMSE. However, there are many smaller subsets that produce approximately the same performance but with fewer predictors. In this case, we might be able to accept a slightly larger error for less predictors.

The \texttt{pickSizeTolerance} determines the absolute best value then the percent difference of the other points to this value. In the case of RMSE, this would be
\[
RMSE_{tol} = 100\times \frac{RMSE - RMSE_{opt}}{RMSE_{opt}}
\]
where $RMSE_{opt}$ is the absolute best error rate. These ``tolerance'' values are plotted in the bottom panel of Figure  \ref{F:tolerance}. The solid triangle is the smallest subset size that is within 10$\%$ of the optimal value.

This approach can produce good results for many of the tree based models, such as random forest, where there is a plateau of good performance for larger subset sizes. For trees, this is usually because unimportant variables are infrequently used in splits and do not significantly affect performance.


\begin{figure}
   \begin{center}
<<tolerancePlot, echo = FALSE, results = hide>>=

smallest <- pickSizeBest(example, metric = "RMSE", maximize = FALSE)

within10Pct <- pickSizeTolerance(example, metric = "RMSE", tol = 10, maximize = FALSE)

minRMSE <- min(example$RMSE)
example$Tolerance <- (example$RMSE - minRMSE)/minRMSE * 100   

pdf("tolerance.pdf", width = 4, height = 5)
par(mfrow = c(2, 1), mar = c(0, 4, 4, 2))

plot(example$Variables[-c(smallest, within10Pct)], 
     example$RMSE[-c(smallest, within10Pct)],
     ylim = extendrange(example$RMSE),
     ylab = "RMSE", xlab = "Variables")

points(example$Variables[smallest], 
       example$RMSE[smallest], pch = 16, cex= 1.3)

points(example$Variables[within10Pct], 
       example$RMSE[within10Pct], pch = 17, cex= 1.3)

 
  with(example, plot(Variables, Tolerance))
     abline(h = 10, lty = 2, col = "darkgrey")
dev.off()
@
      \includegraphics[clip, width = 4.5in]{tolerance}
      \caption{An example of where a smaller subset sizes is not necessarily the best choice. The solid circle in the top panel indicates the subset size with the absolute smallest RMSE. If the percent differences from the smallest RMSE  are calculated (lower panel), the user may want to accept a pre--specified drop in performance as long as the drop is within some limit of the optimal.}         
      \label{F:tolerance} 
    \end{center}
\end{figure}

\subsubsection{The \texttt{selectVar} Function}

After the optimal subset size is determined, this function will be used to calculate the best rankings for each variable across all the resampling iterations (line \ref{A:selectVar}). Inputs for the function are:
\begin{itemize}
\item \texttt{y}: a list of variables importance for each resampling iteration and each subset size (generated by the user--defined \texttt{rank} function). In the example, each each of the cross--validation groups the output of the \texttt{rank} function is saved for each of the \Sexpr{length(subsets)+1} subset sizes (including the original subset). If the rankings are not recomputed at each iteration, the values will be the same within each cross--validation iteration.
\item \texttt{size}: the integer returned by the  \texttt{selectSize} function
\end{itemize}
This function should return a character string of predictor names (of length \texttt{size}) in the order of most important to least important

For random forests, only the first importance calculation (line \ref{A:rank0}) is used since these are the rankings on the full set of predictors. These importances are averaged and the top predictors are returned.
<<rfSelectVar>>=
rfRFE$selectVar
@ 
Note that if the predictor rankings are recomputed at each iteration (line \ref{A:rank1}) the user will need to write their own selection function to use the other ranks.

\subsubsection{The Example}

For random forest, we fit the same series of model sizes as the linear model. The option to save all the resampling results across subset sizes was changed for this model and are used to show the lattice plot function capabilities in Figure \ref{F:rf2}.
<<rf>>=
ctrl$functions <- rfRFE
ctrl$returnResamp <- "all"
set.seed(10)
rfProfile <- rfe(x, y,
                 sizes = subsets,
                 rfeControl = ctrl)
print(rfProfile)
@ 

\begin{figure}
   \begin{center}		
<<rfPlot, echo = FALSE, results = hide>>=
pdf("rf.pdf", width = 5, height = 7)
   trellis.par.set(caretTheme())
   plot1 <- plot(rfProfile, type = c("g", "o"))
   plot2 <- plot(rfProfile, type = c("g", "o"), metric = "Rsquared")
   print(update(plot1, 
                ylab = "Resampled RMSE"), 
         split=c(1,1,1,2), 
         more=TRUE)
   print(update(plot2, 
                ylab = "Resampled R^2"), 
         split=c(1,2,1,2))
dev.off()
@
      \includegraphics[clip]{rf}
      \caption{Performance profiles for random forest.}
      \label{F:rf} 
    \end{center}
\end{figure}

\begin{figure}
   \begin{center}		
<<rfPlot, echo = FALSE, results = hide>>=
pdf("rf2.pdf", width = 6, height = 7.5)
   trellis.par.set(caretTheme())
   plot1 <- xyplot(rfProfile, type = c("g", "p", "smooth"), ylab = "RMSE CV Estimates")
   plot2 <- densityplot(rfProfile, subset = Variables < 5, adjust = 1.25, as.table = TRUE, xlab = "RMSE CV Estimates")
   print(plot1, split=c(1,1,1,2), more=TRUE)
   print(plot2, split=c(1,2,1,2))
dev.off()
@
      \includegraphics[clip]{rf2}
      \caption{Resampling RMSE estimates for random forests across different subset sizes. These plots were generated using \texttt{xyplot(rfProfile)} and \texttt{densityplot(rfProfile, subset = Variables < 5)}}
      \label{F:rf2} 
    \end{center}
\end{figure}

\clearpage

\section{Feature Selection Using Univariate Filters}

Another approach to feature selection is to pre--screen the predictors using simple univariate statistical methods then only use those that pass some criterion in the subsequent model steps. Similar to recursive selection, cross-validation of the subsequent models will be biased as the remaining predictors have already been evaluate on the data set. Proper performance estimates via resampling should include the feature selection step.

As an example, it has been suggested for classification models, that predictors can be filtered by conducting some sort of $k$--sample test (where $k$ is the number of classes) to see if the mean of the predictor is different between the classes. Wilcoxon tests, $t$-tests and ANOVA models are sometimes used. Predictors that have statistically significant differences between the classes are then used for modeling.

The caret function \texttt{sbf} (for selection by filter) can be used
to cross--validate such feature selection schemes. Similar to
\texttt{rfe}, functions can be passed into \texttt{sbf} for the
computational components: univariate filtering, model fitting, prediction and performance summaries (details are given below).

The function is applied to the entire training set and also to different resampled versions of the data set. From this, generalizable estimates of performance can be computed that properly take into account the feature selection step. Also, the results of the predictor filters can be tracked over resamples to understand the uncertainty in the filtering.

\subsection{Basic Syntax}

Similar to the \texttt{rfe} function, the syntax for \texttt{sbf} is:
\begin{Verbatim}
   sbf(predictors, outcome, sbfControl = sbfControl(), ...)
   ## or 
   sbf(formula, data, sbfControl = sbfControl(), ...)
\end{Verbatim}
In this case, the details are specificed using the \texttt{sbfControl}
function. Here, the argument \texttt{functions} dictates what the
different components should do. This argument should have elements
called \texttt{filter}, \texttt{fit}, \texttt{pred} and \texttt{summary}.

\subsubsection{The \texttt{score} Function}\label{S:filter}

This function takes as inputs the predictors and the outcome in
objects called \texttt{x} and \texttt{y}, respectively. The output
should be a named vector of scores where the names correspond to the
column names of \texttt{x}. 

There are two built--in functions called \texttt{anovaScores} and
\texttt{gamScores}. \texttt{anovaScores} treats the outcome as the
independent variable and the predictor as the outcome. In this way,
the null hypothesis is that the mean predictor values are equal across
the different classes. For regression, \texttt{gamScores} fits a smoothing spline in
the predictor to the outcome  using a generalized additive model and tests to see if there is any
functional relationship between the two. In each function the
p--value is used as the score.

\subsubsection{The \texttt{filter} Function}\label{S:filter}

This function takes as inputs the scores coming out of the
\texttt{score} function (in an argument called \texttt{score}). The
function also has the training set data as inputs (arguments are called \texttt{x} and \texttt{y}). 
The output
should be a named logical vector where the names correspond to the
column names of \texttt{x}. Columns with values of \texttt{TRUE} will
be used in the subsequent model.


\subsubsection{The \texttt{fit} Function}

The component is very similar to the function described in Section
\ref{S:fit}. For \texttt{sbf}, there are no \texttt{first} or
\texttt{last} arguments. The function should have arguments
\texttt{x}, \texttt{y} and \texttt{...}. The data within \texttt{x}
have been filtered using the \texttt{filter} function described
above. The output of the \texttt{fit} function should be a fitted model.


With some data sets, no predictors will survive the filter. In these cases, a model with predictors cannot be computed, but the lack of viable predictors should not be ignored in the final results. To account for this issue, \texttt{caret} contains a model function called \texttt{nullModel} that fits a simple model that is independent of any of the predictors. For problems where the outcome is numeric, the function predicts every sample using the simple mean of the training set outcomes. For classification, the model predicts all samples using the most prevalent class in the training data.

This function can be used in the \texttt{fit} component function to
``error-trap'' cases where no predictors are selected. For example,
there are several built--in functions for some models. The object
\texttt{rfSBF} is a set of functions that may be useful for fitting
random forest models with filtering. The \texttt{fit} function here
uses \texttt{nullModel} to check for cases with no predictors:
<<rfSBFfit>>=
rfSBF$fit
@ 

\subsubsection{The \texttt{summary} and \texttt{pred} Functions}

The \texttt{summary} function is used to calculate model performance on held--out samples. The \texttt{pred} function is used to predict new samples using the current
predictor set. The arguments and outputs for these two functions are identical to the previously discussed \texttt{summary} and \texttt{pred} functions in Sections \ref{S:summary} and \ref{S:pred}, respectively.

\subsection{The Example}

Returning to the example from (Friedman, 1991), we can fit another
random forest model with the predictors pre--filtered using the
generalized additive model approach described in Section \ref{S:filter}.

<<sbf>>=
set.seed(10)
rfWithFilter <- sbf(x, y,
                    sbfControl = sbfControl(
                      functions = rfSBF,
                      method = "cv",
                      verbose = FALSE))
print(rfWithFilter)
@ 

In this case, the training set indicated that
\Sexpr{length(predictors(rfWithFilter))} should be used in the random
forest model, but the resampling results indicate that there is some
variation in this number. Some of the informative predictors are
used, but a few others are erroneous retained.

Similar to \texttt{rfe}, there are methods for \texttt{predictors},
\texttt{densityplot}, \texttt{histogram} and \texttt{varImp}.

\section{Visualizing Resampling Results}

As shown in the main \texttt{caret} vignette, the \texttt{resamples}
class can be used to compare the resampling results of objects
resulting from either \texttt{train()}, \texttt{rfe()} or \texttt{sbf()}.

First, the resample results are collected together:
<<resamps1>>=
bootValues <- resamples(
                        list(lmRFE = lmProfile, 
                             rfRFE = rfProfile, 
                             rfFilter = rfWithFilter))
@ 
After this, there are several Lattice plot methods: \texttt{xyplot},
\texttt{dotplot}, \texttt{densityplot}, \texttt{bwplot},
\texttt{splom} and \texttt{parallel}. Figure \ref{F:resamps1} show a
parallel coordinate of the bootstrapped $R^2$ values. 

The \texttt{diff} function can be applied to these objects to produce
confidence intervals and hypothesis testing on the differences:
<<resamps2>>=
differences <- diff(bootValues)
summary(differences)
@ 
There are several plot methods for the differences: 
\texttt{dotplot}, \texttt{densityplot}, \texttt{bwplot} and
\texttt{levelplot}. 

\begin{figure}
   \begin{center}		
<<resamps3, echo = FALSE, results = hide>>=
pdf("resamps1.pdf", width = 5, height = 5)
   trellis.par.set(caretTheme())
   print(splom(bootValues, metric = "Rsquared"))
dev.off()
@
      \includegraphics[clip]{resamps1}
      \caption{Scatterplot matrix of the bootstrapped $R^2$ values
        from the random forests model with univariate filters and the RFE models using linear least squares and random forests. These plots were generated using \texttt{parallel(splom, metric = "Rsquared")}}
      \label{F:resamps1} 
    \end{center}
\end{figure}


\section{Session Information}

<<session, echo=FALSE, results=tex>>=
toLatex(sessionInfo())
@ 

\section{References}

\begin{description}    
  \item Ambroise, C. and McLachlan, J. H. (2002) ``Selection bias in gene extraction on the basis of microarray gene-expression data,'' {\it Proceedings of the National Academy of Science}, 99, 6562--6566
  
  \item Friedman, J. H. (1991) ``Multivariate adaptive regression splines (with discussion),'' {\it Annals of Statistics}, 19, 1--141
  \item  Svetnik, V., Liaw, A. , Tong, C amd Wang, T. (2004) ``Application of Breiman's random forest to modeling structure-activity relationships of pharmaceutical molecules,'' {\it Multiple Classier Systems, Fifth International Workshop}, 3077, 334--343
  \end{description}

\end{document}



