% \VignetteIndexEntry{caret Manual -- Model Building}
% \VignetteDepends{caret}
% \VignettePackage{caret}
\documentclass[12pt]{article}
\usepackage{colortbl}
\usepackage{amsmath}
\usepackage[pdftex]{graphicx}
\usepackage{color}
\usepackage{xspace}
\usepackage{fancyvrb}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage{booktabs} 
\usepackage{longtable} 
\usepackage[boxed, linesnumbered]{algorithm2e}
\usepackage[
         colorlinks=true,
         linkcolor=blue,
         citecolor=blue,
         urlcolor=blue]
         {hyperref}
\usepackage{lscape}

\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{\codeheading}[1]{\mbox{\color{darkblue}\texttt{#1}}}
\newcommand{\code}[1]{\mbox{\footnotesize\color{darkblue}\texttt{#1}}}
\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
\renewcommand{\pkg}[1]{{\textsf{#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}{0truein}
\setlength{\parskip}{0.10truein}


\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontshape=sl,formatcom=\color{darkblue}} 
\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}

\fvset{fontsize=\footnotesize}
 

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

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


\begin{document}

\maketitle
\renewcommand{\baselinestretch}{.6}

\tableofcontents

\thispagestyle{empty}
	
\vspace{.2in}

\renewcommand{\baselinestretch}{1}
 
	
	
<<loadLibs, results = hide, echo = FALSE>>=
library(caret)
library(kernlab)
library(gbm)
library(ipred)
library(grid)
library(Hmisc)
library(randomForest)
data(BloodBrain)
data(mdrr)

getInfo <- function(what = "Suggests")
{
  text <- packageDescription("caret")[what][[1]]
  text <- gsub("\n", ", ", text, fixed = TRUE)
  text <- gsub(">=", "$\\\\ge$", text, fixed = TRUE)
  eachPkg <- strsplit(text, ", ", fixed = TRUE)[[1]]
  eachPkg <- gsub(",", "", eachPkg, fixed = TRUE)
  out <- paste("\\\\pkg{", eachPkg[order(tolower(eachPkg))], "}", sep = "")
  paste(out, collapse = ", ")
}
@


The \pkg{caret} package (short for
{\bf{\color{blue}{c}}}lassification {\bf{\color{blue}{a}}}nd
{\bf{\color{blue}{re}}}gression {\bf{\color{blue}{t}}}raining)
contains functions to streamline the model training process for
complex regression and classification problems. The package utilizes a
number of R packages but tries not to load them all at package
start-up\footnote{By adding formal package dependencies, the package
  startup time can be greatly decreased}. The package ``suggests''
field includes: \Sexpr{getInfo("Suggests")}. \pkg{caret} loads
packages as needed and assumes that they are installed. Install
\pkg{caret} using  
\begin{Verbatim}
install.packages("caret", dependencies = c("Depends", "Suggests"))
\end{Verbatim}
to ensure that all the needed packages are installed.

\section{Model Training and Parameter Tuning}\label{S:train}

\pkg{caret} has several functions that attempt to streamline the model building and evaluation process. 

The \pkg{caret} function can be used to
\begin{itemize}
   \item evaluate, using resampling, the effect of model tuning parameters on performance
   \item choose the ``optimal'' model across these parameters 
   \item estimate model performance from a training set
\end{itemize}

More formally:

\begin{algorithm}[H]
   \label{A:tune}
   \SetAlgoLined
   \RestyleAlgo{plain}
   \DontPrintSemicolon
  Define  sets of model parameter values to evaluate \nllabel{A:grid}\;
  \For{each parameter set}{

    \For{each resampling iteration}{
      Hold--out specific samples \nllabel{A:resample} \;
      
      [Optional] Pre--process the data\;
      Fit the model on the remainder\;
      
      Predict the hold--out samples\;
    }
    Calculate the average performance across hold--out predictions \nllabel{A:perf}
  }
  Determine the optimal parameter set \nllabel{A:best}\;
  
  Fit the final model to all the training data using the optimal
  parameter set\;  
\end{algorithm}


First, a specific model must be chosen. Currently,
\Sexpr{length(unique(modelLookup()$model))} are available using
\pkg{caret}; see Tables \ref{T:methods1},  \ref{T:methods2} and  \ref{T:methods3} for details. 

In Tables \ref{T:methods1},  \ref{T:methods2} and  \ref{T:methods3}, there are lists of tuning parameters that can
potentially be optimized.  The first step in tuning the model (line
\ref{A:grid} in Algorithm \ref{A:tune}) is to choose a set of
parameters to evaluate. For example, if fitting a Partial Least Squares
(PLS) model, the number of PLS components to evaluate must be specified. 

Once the model and tuning parameter values have been defined, the type
of resampling should be also be specified. Currently, $k$--fold
cross--validation (once or repeated), 
leave--one--out cross--validation and bootstrap
(simple estimation or the 632 rule) 
resampling methods can be used by \pkg{caret}. After resampling,
the process produces a profile of performance measures is available to
guide the user as to which tuning parameter values should be
chosen. By default, the function automatically chooses the tuning
parameters associated with the best value, although different
algorithms can be used (see Section \ref{S:finalMod}). 


\subsection{An Example}

As an example, the multidrug resistance reversal (MDRR) agent data is
used to determine a predictive model for the ``ability of a compound
to reverse a leukemia cell's resistance to adriamycin''
(\href{http://pubs.acs.org/cgi-bin/abstract.cgi/jcisd8/2005/45/i03/abs/ci0500379.html}{Svetnik
  et al, 2003}). For each sample (i.e. compound), predictors are
calculated that reflect characteristics of the molecular
structure. These molecular descriptors are then used to predict assay
results that reflect resistance.  

The data are accessed using \code{data(mdrr)}. This creates a data
frame of predictors called \code{mdrrDescr} and a factor vector with
the observed class called \code{mdrrClass}. 

To start, we will:
 
\begin{itemize}
   \item use unsupervised filters to remove predictors with
     unattractive characteristics (e.g. spare distributions or high
     inter--predictor correlations) 
   \item split the entire data set into a training and test set
\end{itemize}

See the package vignette ``caret Manual -- Data and Functions'' for more details about these operations.

\begin{small}
<<preProc>>=
print(ncol(mdrrDescr))
nzv <- nearZeroVar(mdrrDescr)
filteredDescr <- mdrrDescr[, -nzv]

print(ncol(filteredDescr))

descrCor <- cor(filteredDescr)
highlyCorDescr <- findCorrelation(descrCor, cutoff = .75)
filteredDescr <- filteredDescr[,-highlyCorDescr]

print(ncol(filteredDescr))

set.seed(1)
inTrain <- sample(seq(along = mdrrClass), length(mdrrClass)/2)

trainDescr <- filteredDescr[inTrain,]
testDescr <- filteredDescr[-inTrain,]
trainMDRR <- mdrrClass[inTrain]
testMDRR <- mdrrClass[-inTrain]

print(length(trainMDRR))
print(length(testMDRR))
@
\end{small}

\subsection{Basic Parameter Tuning}\label{S:basic}

By default, simple bootstrap resampling is used for line
\ref{A:resample} in Algorithm \ref{A:tune}. Others are availible, such
as repeated $K$--fold cross--validation. The function
\code{trainControl} can be used to specifiy the type of resampling:

\begin{small}
<<setup>>=
fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 10,
                           ## repeated three times
                           repeats = 3,
                           ## Save all the resampling results
                           returnResamp = "all")
@
\end{small}
More information about \code{trainControl} is given in Section
\ref{S:control}. 

The first two arguments to \pkg{caret} are the predictor and
outcome data objects, respectively. The third argument,
\code{method}, specifies the type of model (see Tables \ref{T:methods1},  \ref{T:methods2} and  \ref{T:methods3}). We will fit a boosted tree model via the \pkg{gbm}
package. The basic syntax for fitting this model using repeated
cross--vlaidation is shown below:

\begin{small}
<<mdrrModel1>>=
gbmFit1 <- train(trainDescr, trainMDRR, 
                 method = "gbm", 
                 trControl = fitControl,
                 ## This last option is actually one
                 ## for gbm() that passes through
                 verbose = FALSE)
gbmFit1
@
\end{small}

For a gradient boosting machine (GBM) model, there are three main
tuning parameters:
\begin{itemize}
\item number of iterations,  {\it i.e. trees},  (called \code{n.trees} in the
  \code{gbm} function)
\item complexity of the tree, called \code{interaction.depth}
\item learning rate: how quickly the algorithm adapts, called
  \code{shrinkage} 
\end{itemize}
The default values tested for this model are shown in the first two
columns (\code{shrinkage} is not shown beause the grid set of
candidate models all use a value of 0.1 for this tuning parameter). 
The column labeled ``\code{Accuracy}'' is the overall agreement rate
averaged over cross--validation iterations. The agreement standard
deviation is also calculated from the cross-validation results. The
column ``\code{Kappa}'' is Cohen's (unweighted) Kappa statistic
averaged across the resampling results. \pkg{caret} works with
specific models (see Tables \ref{T:methods1},  \ref{T:methods2} and  \ref{T:methods3}). For these models, \pkg{caret} can automatically
create a grid of tuning parameters. By default, if $p$ is the number
of tuning parameters, the grid size is $3^p$. For example, regularized
discriminant analysis (RDA) models have two  parameters
(\texttt{gamma} and \texttt{lambda}), both of which lie on $[0,
1]$. The default training grid would produce nine combinations in this
two--dimensional space. 

<<getSeqMods, echo = FALSE, results = hide>>=
seqModList <- paste(paste("\\\\pkg{", unique(subset(modelLookup(), seq)$model), "}", sep = ""), collapse = ", ")
@

\subsection{Notes}

\begin{itemize}
  
\item There is a formula interface (e.g. \code{train(y$\sim$., data =
    someData}) that can be used. One of the issues with a large
  number of predictors is that the objects related to the formula
  which are saved can get very large. In these cases, it is best to
  stick with the non--formula interface described above. 
  
\item The function determines the type of problem
  (classification or regression) from the type of the response
  given in the \code{y} argument.  
  
\item The \code{$\ldots$} option can be used to pass
  parameters to the fitting function. For example, in random
  forest models, you can specify the number of trees to be
  used in the call to \pkg{caret}. 
  
\item For regression models (i.e. a numeric outcome), a similar table would
be produced showing the average root mean squared error and average
$R^2$ value statistic across tuning parameters, otherwise known as
$Q^2$ (see the note below related to this calculation).
For regression models, the classical $R^2$ statistic
  cannot be compared between models that contain an intercept
  and models that do not. Also, some models do not have an
  intercept only null model.	 
  
\item [] To approximate this statistic across different
  types of models,  the square of the correlation between
  the observed and predicted outcomes is used. This means
  that the $R^2$ values produced by \pkg{caret} will not
  match the results of \code{lm} and other functions.  
  
\item [] Also, the correlation estimate does not take into
  account the degrees of freedom in a model and thus does
  not penalize models with more parameters. For some
  models (e.g random forests or on--linear support vector
  machines) there is no clear sense of the degrees of
  freedom, so this information cannot be used in $R^2$ if
  we would like to compare different models. 
  
  
\item The nearest shrunken centroid model of
  \href{http://projecteuclid.org/handle/euclid.ss/1056397488}{Tibshirani
    et al (2003)} is specified using \code{method =
    "pam"}. For this model, there must be at least two samples
  in each class. \pkg{caret} will ignore classes where
  there are less than two samples per class from every model
  fit during bootstrapping or cross--validation (this model
  only). 
  
\item For recursive partitioning models, an initial model is
  fit to all of the training data to obtain the possible
  values of the maximum depth of any node
  (\code{maxdepth}). The tuning grid is created based on
  these values. If \code{tuneLength} is larger than the
  number of possible \code{maxdepth} values determined by
  the initial model, the grid will be truncated to the
  \code{maxdepth} list. 
  
\item [] The same is also true for nearest shrunken centroid
  models, where an initial model is fit to find the range of
  possible threshold values, and MARS models (see the details
  below). 
  
\item For multivariate adaptive regression splines (MARS), the
  \texttt{earth} package is used with a model type of \code{mars}
  or \code{earth} is requested. The tuning parameters used by
  \pkg{caret} are \code{degree} and \code{nprune}. The
  parameter \code{nk} is not automatically specified and, if not
  specified, the default in the \code{earth} function is used.  

\item [] For example, suppose a training set with 40 predictors is
  used with \code{degree = 1} and \code{nprune = 20}. An initial
  model with \code{nk = 41} is fit and is pruned down to 20
  terms. This number includes the intercept and may include
  ``singleton'' terms instead of pairs.  

\item [] Alternate model training schemes can be used by passing
  \code{nk} and/or \code{pmethod} to the \code{earth}
  function. Also, using \code{method = 'gcvEearth'} will use the
  basic GCV pruning procedure and only tune the degree.

\item [] Also, there may be cases where the message such as
  ``specified 'nprune' 29 is greater than the number of available
  model terms 24, forcing 'nprune' to 24'' show up after the model
  fit. This can occur since the \code{earth} function may not
  actually use the number of terms in the initial model as specified
  by \code{nk}. This may be because the \code{earth} function
  removes terms with linear dependencies and  the forward pass
  counts as if terms were added in pairs (although singleton terms
  may be used). By default, the \pkg{caret} function fits and
  initial MARS model is used to determine the number of possible
  terms in the training set to create the tuning grid. Resampled
  data sets may produce slightly different models that do not have
  as many possible values of  \code{nprune}. 
  
\item For the \code{glmboost} and \code{gamboost} functions
  from the \texttt{mboost} package, an additional tuning parameter,
  \code{prune}, is used by train. If \code{prune = "yes"}, the
  number of trees is reduced based on the AIC statistic. If
  \code{"no"}, the number of trees is kept at the value specified
  by the \code{mstop} parameter. See the \texttt{mboost} package
  vignette for more details about AIC pruning. 
  
\item The partitioning model of Molinaro {\it et al.} (2010) has a
  tuning parameter that is the number of partitions in the
  model. The R function \code{partDSA} has the argument
  \code{cut.off.growth} which is described as ``the maximum number
  of terminal partitions to be considered when building the model.''
  Since this is the maximum, the user might ask for a model with $X$
  partitions but the model can only predict $Y < X$. In these cases,
  the model predictions will be based on the largest model available
  ($Y$). 

    \item For generalized additive models, a formula is generated from
      the data. First, predictors with degenerate distributions are
      excluded (via the \code{nearZeroVar} function). Then, the
      number of distinct values for each predictor is calculated. If
      this value is greater than 10, the predictor is entered into the
      formula via a smoothed term (otherwise a linear term is
      used). For models in the \texttt{gam} package, the smooth terms
      have the same amount of smoothing applied to them (i.e. equal
      \code{df} or \code{span} across all the smoothed predictors).
      
    \item For some models (\Sexpr{seqModList}), the \pkg{caret}
      function will fit a model that can be used to derive predictions
      for some sub-models. For example, for MARS (via the
      \code{earth} function), for a fixed degree, a model with a
      maximum number of terms will be fit and the predictions of all of
      the requested models with the same degree and smaller number of
      terms will be computed using \code{update.earth} instead of
      fitting a new model. When the \code{verboseIter} option of the
      \code{trainControl} function is
      used, a line is printed for the ``top--level'' model (instead of
      each model in the tuning grid). 

    \item There are \code{print} and \code{plot} methods for the
      \code{train} class. The
      \code{plot} method visualizes the profile of average resampled
      performance values across the different tuning parameters using
      scatter plots or level plots. See Figures \ref{f:plots1} and
      \ref{f:plots2} for examples. Functions that visualize the individual
      resampling results for \pkg{caret} objects are discussed in
      Section \ref{S:withinMod}.
      
    \item Using the first set of tuning parameters that are
      optimal (in the sense of accuracy or mean squared error),
      \pkg{caret} automatically fits a model with these
      parameters to the entire training data set. That model
      object is accessible in the \code{finalModel} object
      within \pkg{caret}. For example, \verb+gbmFit$finalModel+
      is the same object that would have been produced using a
      direct call to the \code{gbm} function with the final tuning
      parameters. 
      
  
\end{itemize}

There is additional functionality in \pkg{caret} that is described in the next section.

<<getModelData, echo = FALSE, results = hide>>=

models <- modelLookup()
modelKey <- read.delim(system.file("modelKey.txt", package = "caret"))
modelInfo <- merge(models, modelKey, all.y = TRUE)
modelInfo$parameter <- latexTranslate(modelInfo$parameter)

cranRef <- function(x) paste("{\\tt \\href{http://cran.r-project.org/web/packages/", x, "/index.html}{", x, "}}", sep = "")

makeTable <- function(x)
  {
    params <- paste("\\code{", as.character(x$parameter), "}", sep = "", collapse = ", ")
    params <- ifelse(params == "\\code{parameter}", "None", params)

    data.frame(method = paste("\\texttt{", as.character(x$model)[1], "}", sep = ""),
               Package = cranRef(as.character(x$Package)[1]),
               Parameters = params)
               

  }

forTable <- ddply(modelInfo, .(Order, ProblemType, Class, model), makeTable)
forTable$model <- NULL

dual <- subset(forTable, ProblemType == "Dual Use")[,-1]
reg <- subset(forTable, ProblemType == "Regression Only")[,-1]
clas <-  subset(forTable, ProblemType == "Classification Only")[,-1]


@ 

\pagestyle{plain}

{\footnotesize
  
<<dualTable, echo = FALSE, results = tex>>=
dualCats <- unique(as.character(dual$Class))
dualReps <- table(as.character(dual$Class))
dualReps <- dualReps[dualCats]

dual <- as.matrix(dual[, -(1:2)])
rownames(dual) <- NULL

latex(dual, 
      file = "",
      title = "Type",
      rowname = " ",
      rgroup = dualCats,
      n.rgroup = dualReps,
      label = "T:methods1",
      caption = "Dual use models (regression and classification) available in \\code{train}.",
      booktabs = TRUE, 
      landscape = TRUE, 
      longtable = TRUE)
@ 

<<regTable, echo = FALSE, results = tex>>=
regCats <- unique(as.character(reg$Class))
regReps <- table(as.character(reg$Class))
regReps <- regReps[regCats]

reg <- as.matrix(reg[, -(1:2)])
rownames(reg) <- NULL

latex(reg, 
      file = "",
      title = "Type",
      rowname = " ",
      rgroup = regCats,
      n.rgroup = regReps,
      label = "T:methods2",
      caption = "Regression only models available in \\code{train}.",
      booktabs = TRUE, 
      landscape = TRUE, 
      longtable = TRUE)
@ 


<<clasTable, echo = FALSE, results = tex>>=
clasCats <- unique(as.character(clas$Class))
clasReps <- table(as.character(clas$Class))
clasReps <- clasReps[clasCats]

clas <- as.matrix(clas[, -(1:2)])
rownames(clas) <- NULL

latex(clas, 
      file = "",
      title = "Type",
      rowname = " ",
      rgroup = clasCats,
      n.rgroup = clasReps,
      label = "T:methods3",
      caption = "Classification only models available in \\code{train}.",
      booktabs = TRUE, 
      landscape = TRUE, 
      longtable = TRUE)
@ 



}
\pagestyle{fancy}

\section{Customizing the Tuning Process}

There are a few ways to customize the process of selecting
tuning/complexity parameters and building the final model.

\subsection{Pre--Processing Options}\label{S:pp}

As previously mentioned, \pkg{caret} can pre--process the data in
various ways prior to model fitting. The \pkg{caret} function
\code{preProcess} is automatically used. This function can be used
for centering and scaling, imputation (see details below),
applying the spatial sign transformation and feature extraction via
principal component analysis or independent component
analysis. Options to the \code{preProcess} function can be passed
via the \code{trainControl} function.  

These processing steps would be applied during any predictions
generated using \code{predict.train}, \code{extractPrediction} or
\code{extractProbs} (see Section \ref{S:probs}  later in this document). The
pre--processing would {\bf not} be applied to predictions that
directly use the \verb+object$finalModel+ object.  

For imputation, there are two methods currently implemented:
\begin{itemize}
\item $k$--nearest neighbors takes a sample with missing values and
  finds the $k$ closest samples in the training set. The average of
  the $k$ training set values for that predictor are used as a
  substitute for the original data. When calculating the distances to
  the training set samples, the predictors used in the calculation are
  the ones with no missing values for that sample and no missing
  values in the training set.
\item another approach is to fit a bagged tree model for each
  predictor using the training set samples. This is usually a fairly
  accurate model and can handle missing values. When a predictor for a
  sample requires imputation, the values for the other predictors are
  fed through the bagged tree and the prediction is used as the new
  value. This model can have significant computational cost.
\end{itemize}

If there are missing values in the training set, PCA and ICA models
only use complete samples.


\subsection{Alternate Tuning Grids}\label{S:altTune}

The tuning parameter grid can be specified by the user. The argument
\code{tuneGrid} can take a data frame with columns for each tuning
parameter (see Tables \ref{T:methods1},  \ref{T:methods2} and  \ref{T:methods3} for specific details). The column
names should be the same as the fitting function's arguments with a
period preceding the name. For the previously mentioned RDA example, the names would be
\code{.gamma} and \code{.lambda}. \pkg{caret} will tune the
model over each combination of values in the rows. 
 

We can fix the learning rate and evaluate more than three values of
\code{n.trees}: 

\begin{small}
<<mdrrGrid>>=
gbmGrid <-  expand.grid(.interaction.depth = c(1, 3), 
                        .n.trees = c(10, 50, 100, 150, 200, 250, 300), 
                        .shrinkage = 0.1)
@ 

<<mdrrGridFit>>=
set.seed(3)
gbmFit2 <- train(trainDescr, trainMDRR, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 ## Now specify the exact models 
                 ## to evaludate:
                 tuneGrid = gbmGrid)
gbmFit2
@

\end{small}

\subsection{The \code{trainControl} Function}\label{S:control}

The function \code{trainControl} generates parameters that further
control how models are created, with possible values: 
\begin{itemize}
\item \code{method}: The resampling method: \code{boot}, \code{boot632},
  \code{cv}, \code{LOOCV}, \code{LGOCV}, \code{repeatedcv}  and \code{oob}. The
  last value, out--of--bag estimates, can only be used by random
  forest, bagged trees, bagged earth, bagged flexible discriminant
  analysis, or conditional tree forest models. GBM models are not
  included (the \texttt{gbm} package maintainer has indicated that
  it would not be a good idea to choose tuning parameter values
  based on the model OOB error estimates with boosted trees). Also,
  for leave--one--out cross--validation, no uncertainty estimates
  are given for the resampled performance measures. 
\item \code{number and repeats}: \code{number} controls with the 
  number of folds in $K$--fold cross--validation or number of
  resampling iterations for bootstrapping and leave--group--out
  cross--validation. \code{repeats} applied only to repeated  
  $K$--fold cross--validation. Suppose that \code{method = "repeatedcv"}, 
  \code{number = 10} and \code{repeats = 3}, then three separate
  10--fold cross--validations are used as the resampling scheme.
\item \code{verboseIter}: A logical for printing a training log.
\item \code{returnData}: A logical for saving the data into a slot
  called \code{trainingData}. 
\item \code{p}: For leave-group out cross-validation: the training
  percentage 
\item \code{ classProbs}: a logical value determining whether
  class probabilities should be computed for held--out samples
  during resample. Examples of using this argument are given in
  Section \ref{S:altPerf}. 
\item \code{index}: a list with elements for each resampling
  iteration. Each list element is the sample rows used for training
  at that iteration. When these values are not specified,
  \pkg{caret} will generate them. 
\item \code{summaryFunction}: a function to compute alternate
  performance summaries. See Section \ref{S:altPerf} for more
  details. 
\item \code{selectionFunction}: a function to choose the optimal
  tuning parameters. See Section \ref{S:finalMod} for more details
  and examples. 
\item \code{PCAthresh}, \code{ICAcomp} and \code{k}: these are
  all options to pass to the \code{preProcess} function (when used).
\item \code{returnResamp}: a character string containing one of
  the following values: \code{"all"}, \code{"final"} or
  \code{"none"}. This specifies how much of the resampled
  performance measures to save.  
\end{itemize}



\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .4\textwidth]{gbm1}
      \hspace*{.2 in}
      \includegraphics[clip, width = .4\textwidth]{gbm2}    

      \caption{ Examples of output from \code{plot.tain}. {\bf left} a plot produced using
        \code{plot(gbmFit3)} showing the relationship between the
        number of boosting iterations, the interaction depth and the
        resampled classification accuracy {\bf right} the same
        plot, but the Kappa statistic is plotted using
        \code{plot(gbmFit3, metric = "Kappa")}} 
      \label{f:plots1} 
    \end{center}
\end{figure} 

\begin{figure}
  \begin{center}	
    \includegraphics[clip, width = .6\textwidth]{gbm3}      
    \caption{For the boosted tree example in Section \ref{S:altTune}, using
      \code{plot(gbmFit metric = "Kappa", plotType = "level")} shows
      the relationship (using a \code{levelplot}) between the number of
      boosting iterations, the interaction depth and the resampled
      estimate of the Kappa statistic.  }
    \label{f:plots2}         
  \end{center}
\end{figure}   


\subsection{Alternate Performance Metrics}\label{S:altPerf}

The user can change the metric used to determine the best settings. By
default, RMSE and $R^2$ are computed for regression while accuracy and
Kappa are computed for classification. Also by default, the parameter
values are chosen using RMSE and accuracy, respectively  for
regression and classification. The \code{metric} argument of the
\pkg{caret} function allows the user to control which the
optimality criterion is used. For example, in problems where there are
a low percentage of samples in one class, using \code{metric =
  "Kappa"} can improve quality of the final model. 

If none of these parameters are satisfactory, the user can also
compute custom performance metrics. The \code{trainControl} function
has a argument called \code{summaryFunction} that specifies a
function for computing performance. The function should have these
arguments: 
\begin{itemize}
\item \code{data} is a reference for a data frame or matrix with
  columns called \code{obs} and \code{pred} for the observed and
  predicted outcome values (either numeric data for regression or
  character values for classification). Currently, class probabilities
  are not passed to the function. The values in data are the held--out
  predictions (and their associated reference values) for a single
  combination of tuning parameters. If the \code{classProbs}
  argument of the \code{trainControl} object is set to
  \code{TRUE}, additional columns in \code{data} will be present
  that contains the class probabilities. The names of these columns
  are the same as the class levels.
\item \code{lev} is a character string that has the outcome factor
  levels taken from the training data. For regression, a value of
  \code{NULL} is passed into the function. 
\item \code{model} is a character string for the model being used
  (i.e. the value passed to the \code{method} value of
  \pkg{caret}). 
\end{itemize}
The output to the function should be a vector of numeric summary
metrics with non--null names.  By default, \pkg{caret} evaluate classification models in terms of
the predicted classes. Optionally, class probabilities can also be
used to measure performance. To obtain predicted class probabilities
within the resampling process, the argument \code{classProbs} in
\code{trainControl} must be set to \code{TRUE}. This merges
columns of probabilities into the predictions generated from each
resample (there is a column per class and the column names are the
class names).

As shown in the last section, custom functions can be used to
calculate performance scores that are averaged over the
resamples. Another built--in function, \code{twoClassSummary}, will
compute the sensitivity, specificity and area under the ROC curve.

\clearpage
<<summaryFunc>>=
twoClassSummary
@ 

To rebuild the boosted tree model using this criterion, we
can see the relationship between the tuning parameters and
the area under the ROC curve using the following code: 
\begin{small}
<<setup>>=
fitControl <- trainControl(method = "repeatedcv",
                           number = 10,
                           repeats = 3,
                           returnResamp = "all",
                           ## Estimate class probabilities
                           classProbs = TRUE,
                           ## Evaluate performance using 
                           ## the following function
                           summaryFunction = twoClassSummary)
@
\end{small}
<<reTune>>=
set.seed(3)
gbmFit3 <- train(trainDescr, trainMDRR, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 tuneGrid = gbmGrid,
                 ## Specify which metric to optimize
                 metric = "ROC")
gbmFit3
@ 

In this case, the average area under the ROC curve associated with the
optimal tuning parameters was
\Sexpr{round(caret:::getTrainPerf(gbmFit3)[1,"TrainROC"], 3)} across
the \Sexpr{length(gbmFit3$control$index)} resamples.


\subsection{Choosing the Final Model}\label{S:finalMod}

Another method for customizing the tuning process is to modify the
algorithm that is used to select the ``best'' parameter values, given
the performance numbers. By default, the \pkg{caret} function
chooses the model with the largest performance value (or smallest, for
mean squared error in regression models). Other schemes for selecting
model can be used.  Breiman et al (1984) suggested the ``one standard
error rule'' for simple tree--based models. In this case, the model
with the best performance value is identified and, using resampling,
we can estimate the standard error of performance. The final model
used was the simplest model within one standard error of the
(empirically) best model. With simple trees this makes sense, since
these models will start to over-fit as they become more and more
specific to the training data. 

\pkg{caret} allows the user to specify alternate rules for
selecting the final model. The argument \code{selectionFunction}
can be used to supply a function to algorithmically determine the
final model. There are three existing functions in the package:
\code{best} is chooses the largest/smallest value, \code{oneSE}
attempts to capture the spirit of Breiman et al (1984) and
\code{tolerance} selects the least complex model within some percent
tolerance of the best value. See \code{?best} for more details. 

User--defined functions can be used, as long as they have the
following arguments: 
\begin{itemize}
  \item \code{x} is a data frame containing the tune parameters and
    their associated performance metrics. Each row corresponds to a
    different tuning parameter combination 
    \item \code{metric} a character string indicating which
      performance metric should be optimized (this is passed in
      directly from the \code{metric} argument of \pkg{caret}.  
      \item \code{maximize} is a single logical value indicating
        whether larger values of the performance metric are better
        (this is also directly passed from the call to
        \pkg{caret}). 
  \end{itemize}
  The function should output a single integer indicating which row in
  \code{x} is chosen. 

<<bestGBM, echo = FALSE, results = hide>>=
printSelected <- function(x)
  {
    tmp <- x$bestTune
    names(tmp) <- gsub(".", " ", names(tmp), fixed = TRUE)
    tmp <- paste(names(tmp), "=", tmp)
    paste(tmp, collapse = ", ")
  }
getTrainPerf <- function(x)
  {
    bst <- x$bestTune
    names(bst) <- substring(names(bst), 2)
    merge(bst, x$results)
  }
@   
  
As an example, if we chose the previous boosted tree model on the basis of
overall accuracy (Figure \ref{f:plots1}), we would choose:
\Sexpr{printSelected(gbmFit3)}. However, the scale in this plots is
fairly tight, with accuracy values ranging from
\Sexpr{round(min(gbmFit3$results$ROC), 3)} to
\Sexpr{round(max(gbmFit3$results$ROC), 3)}. A less complex model
(e.g. fewer, more shallow trees) might also yield acceptable
accuracy. 

The tolerance function could be used to find a less complex
model based on $(x-x_{best})/x_{best}\times 100$, which is the percent
difference. For example, to select parameter values based on a $2\%$ loss of performance: 
<<tolerance>>=
whichTwoPct <- tolerance(gbmFit3$results, "ROC", 2, TRUE)  
cat("best model within 2 pct of best:\n")
gbmFit3$results[whichTwoPct,]
@ 
This indicates that we can get a less complex model with and accuracy
of \Sexpr{round(gbmFit3$results[whichTwoPct,"ROC"], 3)} (compared
to the ``pick the best'' value of 
 \Sexpr{round(getTrainPerf(gbmFit3)$ROC, 3)}). 

The main issue with these functions is related to ordering the models
from simplest to complex. In some cases, this is easy (e.g. simple
trees, partial least squares), but in cases such as this model, the ordering of
models is subjective. For example, is a boosted tree model using 100
iterations and a tree depth of 2 more complex than one with 50
iterations and a depth of 8? The package makes some choices regarding
the orderings. In the case of boosted trees, the package assumes that
increasing the number of iterations adds complexity at a faster rate
than increasing the tree depth, so models are ordered on the number of
iterations then ordered with depth. See \code{?best} for more
examples for specific models. 

<<makePlots, results = hide, echo = FALSE>>=
pdf(paste(getwd(), "/gbm1.pdf", sep = ""), width = 5, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3))
dev.off()
pdf(paste(getwd(), "/gbm2.pdf", sep = ""), width = 5, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3, metric = "Spec"))
dev.off()
pdf(paste(getwd(), "/gbm3.pdf", sep = ""), width = 8, height = 4)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3, meric = "ROC", plotType = "level"))
dev.off()

@

\subsection{Parallel Processing}\label{S:parallel}

If a model is tuned using resampling, the number of model fits can
become large as the number of tuning combinations increases (see the
two loops in Algorithm \ref{A:tune}). To reduce the training time,
parallel processing can be used.  
For example, to train the gradient boosting machine model in Section
\ref{S:basic}, each of the \Sexpr{nrow(gbmFit1$results)} candidate
models was fit to \Sexpr{length(gbmFit1$control$index)} separate
resamples. Since each resample is independent of the other, these
\Sexpr{nrow(gbmFit1$results)*length(gbmFit1$control$index)}  models
could be computed in parallel. 

R has several packages that facilitates parallel processing when
multiple processors are available (see Schmidberger et al., 2009). 
\pkg{caret} can be used to build multiple models simultaneously. 
As of \pkg{caret} version 4.99, a new parallel processing framework
is used to increase the computational efficiency. The \pkg{foreach}
package allows parallel computations using several different
technologies. Although the execution times using \pkg{foreach} is
similar to the framework used prior to version 4.99, there are a few advantages:
\begin{itemize}
  \item the call to \code{train} (or \code{rfe} or \code{sbf}) does
    not change. Parallel ``backends'' are registered with
    \pkg{foreach} prior to the call to \code{train}
  \item compared to the kludgy techniques in \pkg{caret} prior to
    version 4.99, \pkg{foreach} does a much better job of managing memory.
  \item \pkg{foreach} code can be added in several places in the
    package and nest parallelism can be used.
\end{itemize}


For example, to use the \pkg{multicore} package to parallelize the
computations, invoking these commands prior to \code{train} would
split the computations into two workers:
\begin{small}
\begin{Schunk}
\begin{Sinput}
> library(doMC)
> registerDoMC(2)
\end{Sinput}
\end{Schunk}
\end{small}

One common metric used to assess the efficacy of parallelization is
$speedup = T_{seq}/T_{par}$, where $T_{seq}$ and $T_{par}$ denote the
execution times to train the model serially and in parallel,
respectively. Excluding systems with sophisticated shared memory
capabilities, the maximum possible speedup attained by parallelization
with $P$ processors is equal to $P$. Factors affecting the speedup
include the overhead of starting the parallel workers, data transfer,
the percentage of the algorithm's computations that can be done in
parallel, etc.   


\begin{figure}[t]
  \begin{center}		
    \includegraphics[clip, width = .4\textwidth]{min}    
    \includegraphics[clip, width = .4\textwidth]{speedup}
    \caption{Training time profiles using parallel processing via
      \pkg{caret} for a benchmarking data set run on a 16 core
      machine. The left panel shows the elapsed time to train
      a model using single or multiple
      processors. The panel on the right shows the ``speedup,''
      defined to be the time for serial execution divided by the
      parallel execution time. The ``old'' line corresponds to version
      4.98 and below while the ``new'' curve is using \pkg{foreach}.
    }
    \label{f:parallel}         
  \end{center}
\end{figure}

Figure \ref{f:parallel} shows the results of a benchmarking study run
on a 16 core machine. In the left panel, the actual training time for
a radial basis function SVM model for a data set with 5,000 samples
and 400 predictors. The model was tuned over 10 values of the cost
parameter using 50 bootstrap samples. The ``new'' curve corresponds to
the \pkg{foreach} infrastructure. 

One downside to parallel processing in this manner is that the dataset
is held in memory for every node used to train the model. For example,
if parallelism is used to compute the results from 50 bootstrap
samples using $P$ processors, $P$ data sets are held in memory. For
large datasets, this can become a problem if the additional processors
are on the same machines where they are competing for the same
physical memory. The old codebase starts to slow down around 10
workers due to exhausting the physical memory on the machine. The new
codebase does a better job at managing memory with no additional slowdown.


\section{Extracting Predictions and Class Probabilities}\label{S:probs}

As previously mentioned, objects produced by the \pkg{caret}
function contain the ``optimized'' model in the \code{finalModel}
sub--object. Predictions can be made from these objects as usual. In
some cases, such as \code{pls} or \code{gbm} objects, additional
parameters from the optimized fit may need to be specified. In these
cases, the \pkg{caret} objects uses the results of the parameter
optimization to predict new samples. 

For example, we can load the Boston Housing data:

\begin{small}
<<bhExample>>=
library(mlbench)
data(BostonHousing)
# we could use the formula interface too
bhDesignMatrix <-  model.matrix(medv ~. - 1, BostonHousing)
@
\end{small}

\noindent split the data into random training/test groups:

\begin{small}
<<bhSplit>>=
set.seed(4)
inTrain <- createDataPartition(BostonHousing$medv, p = .8, list = FALSE, times = 1)
trainBH <- bhDesignMatrix[inTrain,]
testBH <- bhDesignMatrix[-inTrain,]

trainMedv <- BostonHousing$medv[inTrain]
testMedv <- BostonHousing$medv[-inTrain]
@
\end{small}

\noindent fit partial least squares and multivariate adaptive regression spline models:

\begin{small}
<<bhModels>>=
set.seed(5)
plsFit <- train(trainBH, trainMedv, 
                "pls", 
                preProcess = c("center", "scale"),
                tuneLength = 10, 
                trControl = trainControl(verboseIter = FALSE,
                                         returnResamp = "all"))
set.seed(5)
marsFit <- train(trainBH, trainMedv, 
                 "earth", 
                 tuneLength = 10, 
                 trControl = trainControl(verboseIter = FALSE,
                                          returnResamp = "all"))

@
\end{small} 
To obtain predictions  for the MARS model, \code{predict.earth} can be
used. 

\begin{small}
<<marPrediction1>>=
marsPred1 <- predict(marsFit$finalModel, newdata = testBH)
head(marsPred1)
@ 
\end{small}
Alternatively, \code{predict.train} can be used to get a vector of predictions for the optimal model only:
\begin{small}
<<marsPrediction2>>=
marsPred2 <- predict(marsFit, newdata = testBH)
head(marsPred2)
@
\end{small}
Note that the \code{plsFit} object used pre--processing. In this
case, we cannot directly call \code{predict.mvr} and expect to get
the same answers as \code{predict.train}. The latter function knows
that centering and scaling is required and execute these calculations
on the new samples, whereas \code{predict.mvr} does not. For the
\code{pls} function, there is an argument called \code{scale} that
can be used instead of the pre--processing options in the
\pkg{caret} function.

For multiple models, the objects can be grouped using a list and predicted simultaneously:
\begin{small}
<<bhPrediction1>>=
bhModels <- list(pls = plsFit, mars = marsFit)

bhPred1 <- predict(bhModels, newdata = testBH)
str(bhPred1)
@
\end{small}
In some cases,observed outcomes and their associated predictions may
be needed for a set of models. In this case,
\code{extractPrediction} can be used. This function takes a list of
models and test and/or unknown samples as inputs and returns a data
frame of predictions: 
\begin{small}
<<bhPrediction1>>=
allPred <- extractPrediction(bhModels,
                             testX = testBH,
                             testY = testMedv)
testPred <- subset(allPred, dataType == "Test")
head(testPred)

ddply(testPred, .(model), defaultSummary)

@
\end{small}
The output of \code{extractPrediction} is a data frame with columns:
   \begin{itemize}   
      \item \code{obs}, the observed data
      \item \code{pred}, the predicted values from each model
      \item \code{model}, a character string (``\code{rpart}'', ``\code{pls}'' etc.)
      \item \code{dataType}, a character string for the type of data:
      \begin{itemize}
         \item ``\code{Training}'' data are the predictions on the training data from
            the optimal model,
         \item ``\code{Test}'' denote the predictions on the test set (if one is specified),
         \item ``\code{Unknown}'' data are the predictions on the unknown samples (if specified). 
         Only the predictions are produced for these data. Also, if the quick prediction of the unknowns
         is the primary goal, the argument \code{unkOnly} can be used to only process the unknowns.
      \end{itemize}
   \end{itemize}      

Some classification models can produce probabilities for each
class. The functions \code{predict.train} and \code{predict.list}
can be used with the \code{type = "probs"} argument to produce data
frames of class probabilities (with one column per class). Also, the
function \code{extractProbs} can be used to get these probabilities
from one or more models. The results are very similar to what is
produced by \code{extractPrediction} but with columns for each
class. The column \code{pred} is still the predicted class from the
model.  


\section{Evaluating Test Sets}

A function, \code{postResample}, can be used obtain the same
performance measures as generated by \pkg{caret} for regression or
classification. 

\subsection{Confusion Matrices}\label{S:confusion}

\pkg{caret} also contains several functions that can be used to
describe the performance of classification models. The functions
\code{sensitivity}, \code{specificity}, \code{posPredValue} and
\code{negPredValue} can be used to characterize performance where
there are two classes. By default, the first level of the outcome
factor is used to define the ``positive'' result (i.e. the event of
interest), although this can be changed.  

The function \code{confusionMatrix} can also be used to summarize
the results of a classification model: 

\begin{small}
<<mdrConfusion>>=
mdrrPredictions <- extractPrediction(list(gbmFit3), testX = testDescr, testY = testMDRR)
mdrrPredictions <- mdrrPredictions[mdrrPredictions$dataType == "Test",]
sensitivity(mdrrPredictions$pred, mdrrPredictions$obs)
confusionMatrix(mdrrPredictions$pred, mdrrPredictions$obs)
@
\end{small}

The ``no--information rate'' is the largest proportion of the observed
classes (there were more actives than inactives in this test set). A
hypothesis test is also computed to evaluate whether the overall
accuracy rate is greater than the rate of the largest class. Also, the
prevalence of the ``positive event'' is computed from the data (unless
passed in as an argument), the detection rate (the rate of true events
also predicted to be events) and the detection prevalence (the
prevalence of predicted events). 

Suppose a $2 \times 2$ table with notation

\begin{tabular}{r|c|c|}
                    & \multicolumn{2}{c}{{\bf Reference}}          \\
\cline{2-3}
         {\bf Predicted}  & Event     & No Event \\
\cline{2-3}
         Event      & A         & B        \\
\cline{2-3}
         No Event   & C         & D      \\
\cline{2-3}  
\end{tabular}

The formulas used here are:
\begin{align}
Sensitivity &= \frac{A}{A+C}\notag \\
Specificity &= \frac{D}{B+D}\notag \\
Prevalence &= \frac{A+C}{A+B+C+D}\notag \\
PPV &= \frac{sensitivity \times prevalence}{((sensitivity \times prevalence) + ((1-specificity) \times(1-prevalence))}\notag \\
NPV &= \frac{specificity \times (1-prevalence)}{((1-sensitivity) \times prevalence) + ((specificity) \times(1-prevalence))}\notag \\
Detection\: Rate &=  \frac{A}{A+B+C+D}\notag \\
Detection\: Prevalence &=  \frac{A + B}{A+B+C+D}\notag 
\end{align}

When there are three or more classes, \code{confusionMatrix} will
show the confusion matrix and a set of ``one--versus--all''
results. For example, in a three class problem, the sensitivity of the
first class is calculated against all the samples in the second and
third classes (and so on). 

Also, a resampled estimate of the training set can also be obtained
using \code{confusionMatrix.train}. For each resampling iteration,
a confusion matrix is created from the hold--out samples and these
values can be aggregated to diagnose issues with the model fit.

For example, in the two--class SVM model used in Section
\ref{S:basic}, we could use:

<<resampCM>>=
confusionMatrix(gbmFit3)
@ 

These values are the percentages that hold--out samples landed in the
confusion matrix during resampling. There are several methods for
normalizing these values. See \code{?confusionMatrix.train} for details.


\subsection*{Plotting Predictions and Probabilities}

Two functions, \code{plotObsVsPred} and \code{plotClassProbs}, are
interfaces to lattice to plot model results. For regression,
\code{plotObsVsPred} plots the observed versus predicted values by
model type and data (e.g. test). See Figure \ref{f:bhPredPlot} for example. For classification data,
\code{plotObsVsPred} plots the accuracy rates for models/data in a
dotplot.  


<<mdrrPlots, echo = FALSE, results = hide>>=
mdrrProbs <- extractProb(list(gbmFit3), testX = testDescr, testY = testMDRR)
mdrrProbs <- mdrrProbs[mdrrProbs$dataType == "Test",]

pdf(paste(getwd(), "/gbmProbs.pdf", sep = ""), width = 9, height = 7)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plotClassProbs(mdrrProbs))
dev.off()
@

\begin{figure}
   \begin{center}      
      \includegraphics[clip, width = .5\textwidth]{gbmProbs}   
      \caption{The predicted class probabilities from a gradient boosting machine fit for the MDRR test set. This plot was created using \code{plotClassProbs(mdrrProbs)}.}
      \label{f:mdrrProbs}     
   \end{center}
\end{figure}  

\begin{figure}
   \begin{center}      
<<bhPredPlot, echo = FALSE, results = hide, fig = FALSE>>=
pdf("bhPredPlot.pdf", width = 8, height = 5)
trellis.par.set(caretTheme(), warn = FALSE)
print(plotObsVsPred(testPred))
dev.off()
@         
      \includegraphics[clip, width = .75\textwidth]{bhPredPlot}   
      \caption{The results of using \code{ plotObsVsPred } to show
        plots of the observed median home price against the
        predictions from two models. The plot shows the training and
        test sets in the same Lattice plot} 
      \label{f:bhPredPlot}         
   \end{center}
\end{figure}


To plot class probabilities, \code{plotClassProbs} will display the
results by model, data and true class (for example, Figure
\ref{f:mdrrProbs}).      

\clearpage

\section{Exploring and Comparing Resampling Distributions}

\subsection{Within--Model}\label{S:withinMod}

There are several Lattice functions than can be used to explore
relationships between tuning parameters and the resampling results for
a specific model:
\begin{itemize}
\item \code{xyplot} and \code{stripplot} can be used to plot resampling statistics
  against (numeric) tuning parameters.
  \item \code{histogram} and \code{densityplot} can also be used
      to look at distributions of the tuning parameters across tuning parameters.
\end{itemize}

For example, the following statements produces the images in Figure \ref{f:marsDesns}.   

\begin{small}
\begin{Schunk}
\begin{Sinput}
> xyplot(marsFit, type= c("g", "p", "smooth"), degree = 2)
> densityplot(marsFit, as.table = TRUE, subset = nprune < 10)
\end{Sinput}
\end{Schunk}
\end{small}

<<makeResampPlots, results = hide, echo = FALSE>>=
pdf(paste(getwd(), "/marsXY.pdf", sep = ""), width = 7, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(xyplot(marsFit, type= c("g", "p", "smooth"), degree = 2))
dev.off()
pdf(paste(getwd(), "/marsDens.pdf", sep = ""), width = 7, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(densityplot(marsFit, as.table = TRUE, subset = nprune < 10))
dev.off()
@


\begin{figure}[hb]
   \begin{center}      
      \includegraphics[clip, width = .7\textwidth]{marsXY}   
      
      \vspace{.3 in}      

       \includegraphics[clip, width = .7\textwidth]{marsDens}   
      \caption{Scatter plots and density plots of the resampled RMSE by the number of
        retained terms for the MARS model fit to the Boston Housing data}
      \label{f:marsDesns}       
   \end{center}
\end{figure}  
\clearpage

\subsection{Between--Models}\label{S:betweenMod}


\pkg{caret} also includes functions to characterize the differences
between models (generated using \pkg{caret}, \code{sbf} or
\code{rfe}) via their resampling distributions. These functions are
based on the work of Hothorn et al. (2005) and Eugster et al (2008).

Using the blood-brain barrier data (see \code{?BloodBrain)}, three regression models were
created: an \code{rpart} tree, a conditional inference tree using
\code{ctree}, M5 rules using \code{M5Rules} and a MARS model using
\code{earth}. We ensure that the models use the same resampling data
sets. In this case, 100 leave--group--out cross--validation was
employed. 

<<bbbModelFits, eval = FALSE>>=
data(BloodBrain)
set.seed(1)
tmp <- createDataPartition(logBBB, p = 0.8, times = 100)
bbbControl <- trainControl(method = "LGOCV", index = tmp, timingSamps = 50)

rpartFit <- train(bbbDescr, logBBB,
                   "rpart", 
                   tuneLength = 16,
                   trControl = bbbControl)

ctreeFit <- train(bbbDescr, logBBB,
                 "ctree2",
                  tuneLength = 10,
                  trControl = bbbControl)

earthFit <- train(bbbDescr, logBBB,
                   "earth",
                   tuneLength = 20,
                   trControl = bbbControl)
m5Fit <- train(bbbDescr, logBBB,
                "M5Rules",
                trControl = bbbControl)  
@ 

<<>>=
print(getwd())
@ 

Given these models, can we make statistical statements about their
performance differences? To do this, we first collect the resampling
results using \code{resamples}. 
<<resamps, eval = FALSE>>=
resamps <- resamples(list(CART = rpartFit,
                          CondInfTree = ctreeFit,
                          MARS = earthFit,
                          M5 = m5Fit))
@ 
<<resampLoad, echo = FALSE, results = hide>>=
load("resamps.RData")
@ 
<<resampFuncs>>=
resamps
summary(resamps)
@ 

There are several Lattice plot methods that can be used to visualize
the resampling distributions: density plots, box--whisker plots,
scatterplot matrices and scatterplots of summary statistics. In the
latter case, the plot consists of a scatterplot between the two
models.  (See Figure
\ref{f:resampleScatter}). In Figure \ref{F:resampleDens}, density plots
of the data are shown. In this figure,  the $R^2$ distributions
indicate that M5 rules and MARS appear to be
similar to one another but different from the two tree--based
models. However, this pattern is inconsistent with the root
mean squared error distributions.

<<resamplePlots>>=
bwplot(resamps, metric = "RMSE")

densityplot(resamps, metric = "RMSE")

xyplot(resamps,
       models = c("CART", "MARS"),
       metric = "RMSE")

splom(resamps, metric = "RMSE")
@ 

<<resamplePlots2, results = hide, echo = FALSE>>=
pdf("resampleScatter.pdf", width = 6, height = 6)
trellis.par.set(caretTheme())
print(xyplot(resamps, models = c("CART", "MARS")))
dev.off()
pdf("resampleDens.pdf", width = 7, height = 4.5)
print(
      densityplot(resamps, 
                  scales = list(x = list(relation = "free")), 
                  adjust = 1.2, 
                  plot.points = FALSE, 
                  auto.key = list(columns = 2)))

dev.off() 
@ 


Since models are fit on the same versions of the training data, it
makes sense to make inferences on the differences between models. In
this way we reduce the within--resample correlation that may exist. We
can compute the differences, then use a simple $t$--test to evaluate
the null hypothesis that there is no difference between models. 

<<diffs>>=
difValues <- diff(resamps)

difValues

summary(difValues)
@ 
Note that these results are consistent with the patterns shown in
Figure \ref{F:resampleDens}; there are more differences in the $R^2$
distributions than in the error distributions.

Several Lattices methods also exist to plot the differences (density
and box--whisker plots) or the inferential results (level and dot
plots). Figures \ref{F:diffLevel} and \ref{F:diffDot} show examples of
level and dot plots.

<<diffPlots>>=
dotplot(difValues)

densityplot(difValues,
            metric = "RMSE",
            auto.key = TRUE,
            pch = "|")

bwplot(difValues,
       metric = "RMSE")

levelplot(difValues, what = "differences")
@ 

<<diffPlots2, results = hide, echo = FALSE>>=
pdf("diffLevel.pdf", width = 7.5, height = 6)
trellis.par.set(caretTheme())
print(levelplot(difValues, what = "differences"))
dev.off()

pdf("diffDot.pdf", width = 6, height = 6)
plotTheme <- caretTheme()
plotTheme$plot.symbol$pch <- 16
plotTheme$plot.line$col <- "black"
trellis.par.set(plotTheme)
print(dotplot(difValues))
dev.off()

@ 


\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .6\textwidth]{resampleScatter}    

      \caption{ Examples of output from \code{xyplot(resamps, models = c("CART", "MARS"))}. }
      \label{f:resampleScatter} 
    \end{center}
\end{figure} 

\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .9\textwidth]{resampleDens}    

      \caption{ Examples of output from
        \code{densityplot(resamps)}. Looking at $R^2$, M5 rules and MARS appear to be
        similar to one another but different from the two tree--based
        models. However, this pattern is inconsistent with the root
        mean squared error distributions.}
      \label{F:resampleDens} 
    \end{center}
\end{figure} 


\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .8\textwidth]{diffLevel}    

      \caption{ Examples of output from \code{levelplot(difValues,
          what = "differences")}. The pair--wise differences in RMSE
        are shown}
      \label{F:diffLevel} 
    \end{center}
\end{figure} 


\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .8\textwidth]{diffDot}    

      \caption{ Examples of output from
        \code{dotplot(difValues)}. The differences in RMSE and their
        associated confidence intervals are shown.}
      \label{F:diffDot} 
    \end{center}
\end{figure} 

\section{Custom Methods for \codeheading{train}}

Although there are currently more than
\Sexpr{length(unique(modelLookup()$model)) -(length(unique(modelLookup()$model)) %% 10)}
  methods available to \code{train}, there may be the need to create
  custom model functions (e.g. testing a new model etc). One
  application of custom models would be to create diverse ensembles of
  models. For example, a set of different classification models may be
  fit to the same data and a ''pick--the--winner'' approach can be
  taken (or the average of the class probabilities could be used, see
  Kuncheva (2004) or Seni and Elder (2010)). \code{train} already has
  a framework for resampling and tuning models and
  \code{predict.train} can be used to encapsulate the ensemble of
  models into one call for prediction.

\subsection*{Why Not Re--Write \codeheading{train} Altogether?}

One could make a strong argument that back--fitting the package to work with custom models is a bit kludgy. This is probably true, but there are several munches of the current code base that, if implemented in a more formal manner, would make the new code base overly complicated. 

For example, the current code base exploits models that can produce predictions for many tuning parameters using a single model object. For example PLS can fit a K component model and make predictions from models with $1\ldots K$ components. This potentially saves us $K-1$ model fits. PLS is fairly fast, but pother models with this same feature (gbm, cubist, etc) are computationally taxing and this saves a lot of time.  

Also, the package currently does many types of resampling and some, such as the bootstrap 632 estimator, are not as simple as others. In this case, an additional model must be fit for all the data to get the apparent error rate. Formalizing this (or starting over) with a new code base would be unnecessarily complex.

\subsection{How To Write Custom Methods}

The user will need to deviate from the standard call in two ways:
\begin{itemize}
  \item use \code{method = "custom"} in the call to \code{train}, and
  \item add the required functions for the model using \code{trainControl}
\end{itemize}

The \code{custom} argument of \code{trainControl} requires a list of
named functions with the following elements: \code{parameters},
\code{model}, \code{prediction}, \code{prob} and \code{sort}.

\subsubsection{The \codeheading{parameters} Argument}

This element is used to specify or generate the models tuning
parameters. This can be done either as a function to generate them or
a data frame of the actual parameters.

{\it Inputs}:
\begin{itemize}
  \item \code{data}: a data frame of the training set data. The
    outcome will be in a column labeled \code{.outcome}. If the
    formula method for \code{train} was invoked, the data passed into
    this function will have been processed (i.e. dummy variables have
    been created etc).
  \item \code{len} an optional parameter passed in form the
    \code{tuneLength} argument to \code{train}
\end{itemize}

{\it Outputs}: a data frame where
\begin{itemize}
  \item all columns start with a dot 
  \item there is at least one row
\end{itemize}

Instead of a function, the final data frame can be passed in

\subsubsection{The \codeheading{model} Function}

This element fits the model and any other functions
(e.g. pre--processing of the data)


{\it Inputs}:
\begin{itemize}
\item \code{data}: a data frame of the training set data. The outcome
  will be in a column labeled \code{.outcome}. If case weights were
  specified in the \code{train} call, these are in the column
  \code{.modelWeights}. If the formula method for \code{train} was
  invoked, the data passed into this function will have been processed
  (i.e. dummy variables have been created etc).
  \item \code{weights} case weights
  \item \code{parameter} a single row data frame with the current
    tuning parameter
  \item \code{levels}: either \code{NULL} or a character vector or
    factor labels
  \item \code{last} a logical vector for the final model fit with the
    selected tuning parameters and the full training set
  \item \code{...} arguments passed form \code{train} to this function
\end{itemize}

{\it Outputs}: a list with at least one element:
\begin{itemize}
  \item \code{fit}: the object corresponding to the trained model
\end{itemize}
Anything else can be attached to this object. If custom
pre--processing is required, this can be estimated in the \code{model}
function and attached to the output list. Subsequent calls to the
\code{prediction} and \code{probability} functions will have the
entire list available, so the processing can be applied to the new
data. 


\subsubsection{The \codeheading{prediction} Function}

This should be a function that produces either a number vector (for
regression) or a factor (or character) vector for classification.

{\it Inputs}:
\begin{itemize}
\item \code{object}: a list with two elements resulting from the
  model function
\item \code{newdata}: a matrix or data frame of predictors to be
  processed through the model (and possibly pre-processing routine)
\end{itemize}

The output should be either a numeric, character or factor vector. For
classification, factors are converted to character elsewhere to ensure
the proper levels are in the output.


\subsubsection{The \codeheading{probability} Argument}

For classification models, this function should generate a data frame of class probabilities. For regression, a value of \code{NULL} can be used. 

{\it Inputs}:
\begin{itemize}
\item \code{object}: a list with two elements resulting from the
  model function
\item \code{newdata}: a matrix or data frame of predictors to be
  processed through the model (and possibly pre-processing routine)
\end{itemize}

The output should be a data frame with these characteristics:
\begin{itemize}
\item as many columns are factor levels
\item column names are the same as the factor levels and in the same order
\end{itemize}

\subsubsection{The \codeheading{sort} Function}

There are cases where multiple tuning parameters yield the same level
of performance. In these situations, \code{train} will choose the
parameters associated with the most simplistic model. This function
should take the grid of tuning parameters and order them from least
complex to most complex.

The input is a data frame of tuning parameters (without the
preceding dot in the name). 

The output is the same data frame sorted appropriately. 

\subsection{An Example}

As an example, suppose we want to test out \code{rpart} models where
we tune over the complexity parameter and the minimum number of
samples in a node to do further splitting (a.k.a \code{minsplit}).

We'll use the Blood--brain barrier data in \pkg{caret} to illustrate.

First, we would need to create a training grid with the candidate
values of \code{cp} and \code{minsplit}. When using the nominal
\code{rpart} method in \code{train}, an initial \code{rpart} model is
created and the unique values of the complexity parameter are obtained
from the sub--object \code{cptable}. We will test two values of
\code{minsplit}: 10 and 30. First, we get the unique
$C_p$ values for \code{minsplit = 10}

<<getCp1>>=
## rpart requires a formula method
tmpData <- bbbDescr
tmpData$logBBB <- logBBB

cpValues10 <- rpart(logBBB ~ ., data = tmpData,
              control = rpart.control(minsplit = 10))$cptable[,"CP"]
cpValues30 <- rpart(logBBB ~ ., data = tmpData,
              control = rpart.control(minsplit = 30))$cptable[,"CP"]
head(cpValues10)
@ 

From these, we will create the tuning grid of \Sexpr{nrow(cpValues10)*2}
canidiate models (3 values of \code{minsplit} for each of the
\Sexpr{nrow(cpValues10)} possible $C_p$ values:
<<cpGrid>>=
rpartGrid <- data.frame(.cp = c(cpValues10, cpValues30),
                        .minsplit = 
                        c(rep(10, length(cpValues10)), 
                          rep(30, length(cpValues30))))

@ 

We can now write a model function:
<<modelFoo>>=
modelFunc <- function(data, parameter, levels, last, ...)
  {
    library(rpart)
    ctrl <- rpart.control(cp = parameter$.cp,
                          minsplit = parameter$.minsplit)
    
    list(fit = rpart(.outcome ~ ., data = data, control = ctrl))                             
  }
@ 
It is a good idea to load the \pkg{rpart} package and anything else needed within the function.

The prediciton function is simple:
<<predFoo>>=
predFunc <- function(object, newdata)
  {
    library(rpart)
    predict(object$fit, newdata)                            
  }
@ 

Sorting by complexity is somewhat subjective. Both parameters govern
how deep the tree can be. We will sort by \code{cp} then \code{minsplit}:
<<sortFoo>>=
sortFunc <- function(x)  x[order(x$cp, x$minsplit),]                           
@ 

Now we can create a control object for \code{train}:
<<trControl, eval = FALSE>>=
ctrl <- trainControl(custom = list(
                       parameters = rpartGrid,
                       model = modelFunc,
                       prediction = predFunc,
                       probability = NULL,
                       sort = sortFunc),
                     method = "cv")

set.seed(581)
customRpart <- train(bbbDescr, logBBB, "custom", trControl = ctrl)
@ 

The \code{predict}, \code{print}, \code{plot} and \code{resamples}
methods work with custom models. In the case of \code{plot.train}, the
axis and key labels will be the parameter names. However, \code{update}
can be used to make the labels more descriptive:

<<plotResults, eval = FALSE>>=
rpartPlot <- plot(customRpart, xTrans = log10)
rpartPlot <- update(rpartPlot, xlab = "Complexity Parameter")
@ 




\section{Session Information}

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


\section{References}

\begin{description}

  
\item Breiman, Friedman, Olshen, and Stone. (1984) {\it Classification
    and Regression Trees}. Wadsworth. 

\item Eugster et al. (2008),  ``Exploratory and inferential analysis
  of benchmark experiments, '' {\it Ludwigs-Maximilians-Universitat
    Munchen, Department of Statistics, Tech. Rep} vol. 30 
  
  
\item Hothorn et al. (2005),  ``The design and analysis of
  benchmark experiments, '' {\it Journal of Computational and
    Graphical Statistics}, 14,  675--699

  \item Kuncheva (2004), {\it Combining Pattern Classifiers: Methods
      and Algorithms}. Wiley.
    
  \item Molinaro et al. (2010), ``partDSA:
    deletion/substitution/addition algorithm for partitioning the
    covariate space in prediction,'' {\it Bioinformatics}, 26, 1357--1363

\item Rand (1971),  ``Objective criteria for the evaluation of
  clustering methods,'' {\it Journal of the American Statistical
    Association} 66, 846--850. 
    
\item Schmidberger et al. (2009), `` State-of-the-art in Parallel
  Computing with R,'' {\it Journal of Statistical Software}, 31
  
\item Seni and Elder (2010) {\it Ensemble Methods in Data Mining:
    Improving Accuracy Through Combining Predictions}.
  Morgan and Claypool Publishers
  
\item Svetnik, V., Wang, T., Tong, C., Liaw, A., Sheridan, R. P. and
  Song, Q. (2005), ``Boosting: An ensemble learning tool for compound
  classification and QSAR modeling,'' {\it Journal of Chemical
    Information and Modeling}, 45, 786 --799. 
  
\item Tibshirani, R., Hastie, T., Narasimhan, B., Chu, G. (2003),
  ``Class prediction by nearest shrunken centroids, with applications
  to DNA microarrays,'' {\it  Statistical Science}, 18, 104--117. 

\end{description}

\end{document}
