% \VignetteIndexEntry{caret Manual -- Data and Functions}
% \VignetteDepends{caret}
% \VignettePackage{caret}
\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage[pdftex]{graphicx}
\usepackage{color}
\usepackage{xspace}
\usepackage{fancyvrb}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage{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{\ReallySmallSpacing}{\renewcommand{\baselinestretch}{.6}\Large\normalsize}
\newcommand{\SmallSpacing}{\renewcommand{\baselinestretch}{1.1}\Large\normalsize}

\providecommand{\SetAlgoLined}{\SetLine}

\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{The {\tt caret} Package -- Data Sets and Miscellaneous Functions}
\rhead{}	
\lfoot{}
\cfoot{}
\rfoot{\thepage\ of \pageref{LastPage}}
\renewcommand{\headrulewidth}{1pt}
\renewcommand{\footrulewidth}{1pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Data Sets and Miscellaneous Functions in the {\tt caret} Package}
\author{Max Kuhn \\ max.kuhn@pfizer.com}


\begin{document}

\maketitle

\thispagestyle{empty}
	
\section{Introduction}

<<loadLibs, results = hide, echo = FALSE>>=
library(caret)
library(kernlab)
library(gbm)
library(ipred)
library(grid)
library(pls)
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]]
  
  out <- paste("\\\\texttt{", eachPkg[order(tolower(eachPkg))], "}", sep = "")
  paste(out, collapse = ", ")
}

@

The \texttt{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")}. \texttt{caret} loads packages as needed and assumes that they are installed. Install \texttt{caret} using 
\begin{Verbatim}
install.packages("caret", dependencies = c("Depends", "Suggests"))
\end{Verbatim}
to ensure that all the needed packages are installed.


\section{Example Data}

There are a few data sets included in \texttt{caret}. The first four
are computational chemistry problems where the object is to relate the
molecular structure of compounds (via molecular descriptors) to some
property of interest (\href{http://dx.doi.org/10.1016/S1359-6446(99)01451-8}{Blake (2000)}, Clark and Pickett (2000). Similar data sets can be found in the
\texttt{QSARdata} R pacakge. 

\subsection{Blood--Brain Barrier Data}


\href{http://www.springerlink.com/content/72j377175n536768/?p=f546488cc8fa4ec7a3d4911eb20adb3c&pi=0} {Mente and Lombardo (2005)} developed models to predict the log of the ratio of the concentration of a compound in the brain and the concentration in blood. For each compound, they computed three sets of molecular descriptors: MOE 2D, rule--of--five and Charge Polar Surface Area (CPSA). In all, \Sexpr{dim(bbbDescr)[2]} descriptors were calculated.  Included in this package are \Sexpr{dim(bbbDescr)[1]} non--proprietary literature compounds. The vector \texttt{logBBB} contains the log concentration ratio and the data fame \texttt{bbbDescr} contains the descriptor values.

\subsection{COX-2 Activity Data}

From  \href{http://pubs.acs.org/cgi-bin/abstract.cgi/jmcmar/2004/47/i22/abs/jm0497141.html}{Sutherland, O'Brien, and Weaver (2003)}: ``A set of 467 cyclooxygenase-2 (COX-2) inhibitors has been assembled from the published work of a single research group, with in vitro activities against human recombinant enzyme expressed as IC50 values ranging from 1 nM to >100 uM (53 compounds have indeterminate IC50
values).'' A set of 255 descriptors (MOE2D and QikProp) were generated. To classify the data, we used a cutoff of $2^{2.5}$ to determine activity.

Using \texttt{data(cox2)} exposes three R objects: \texttt{cox2Descr} is a data frame with the descriptor data, \texttt{cox2IC50} is a numeric vector of IC50 assay values and \texttt{cox2Class} is a factor vector with the activity results.


\subsection{DHFR Inhibition}
  
\href{http://www.springerlink.com/content/q5m5xp1q356p2071/}{Sutherland and Weaver (2004)} discuss QSAR models for dihydrofolate reductase (DHFR) inhibition. This data set contains values for 325 compounds. For each compound, 228 molecular descriptors have been calculated. Additionally, each samples is designated as ``active'' or``inactive''.

The data frame \texttt{dhfr} contains a column called \texttt{Y} with the outcome classification. The remainder of the columns are molecular descriptor values.


\subsection{Multidrug Resistance Reversal (MDRR) Agent Data}

\href{http://pubs.acs.org/cgi-bin/abstract.cgi/jcisd8/2005/45/i03/abs/ci0500379.html}{Svetnik et al (2005)} describe these data: ``Bakken and Jurs studied a set of compounds originally discussed by Klopman et al., who were interested in multidrug resistance reversal (MDRR) agents. The original response variable is a ratio measuring the ability of a compound to reverse a leukemia cell's resistance to adriamycin. However, the problem was treated as a classification problem, and compounds with the ratio >4.2 were considered active, and those with the ratio <= 2.0 were considered inactive. Compounds with the ratio between these two cutoffs were called moderate and removed from the data for two class classification, leaving a set of 528 compounds (298 actives and 230 inactives). (Various other arrangements of these data were examined by Bakken and Jurs, but we will focus on this particular one.) We did not have access to the original descriptors, but we generated a set of 342 descriptors of three different types that should be similar to the original descriptors, using the DRAGON software.'' The data and R code are in the Supplemental Data file for the article.

Using \texttt{data(mdrr)} exposes two R objects: \texttt{mdrrDescr} is a data frame with the descriptor data and \texttt{mdrrClass} is a factor vector with the activity results.

\subsection{Tecator NIR Data}

These data can be found in the datasets section of \href{http://lib.stat.cmu.edu/datasets/tecator} {StatLib}. The data consist of 100 near infrared absorbance spectra used to predict the moisture, fat and protein values of chopped meat.

From \href{http://lib.stat.cmu.edu/datasets/tecator} {StatLib}:

``These data are recorded on a Tecator Infratec Food and Feed Analyzer 
working in the wavelength range 850 -- 1050 nm by the Near Infrared 
Transmission (NIT) principle. Each sample contains finely chopped pure 
meat with different moisture, fat and protein contents.

If results from these data are used in a publication we want you to 
mention the instrument and company name (Tecator) in the publication. 
In addition, please send a preprint of your article to:  Karin Thente, Tecator AB, Box 70, S-263 21 Hoganas, Sweden.''

One reference for these data is Borggaard and Thodberg (1992).

Using \texttt{data(tecator)} loads a 215 $\times$ 100 matrix of absorbance spectra and a 215 $\times$ 3 matrix of outcomes.


\subsection{Fatty Acid Composition Data}

\href{http://dx.doi.org/10.1016/j.chemolab.2004.04.011}{Brodnjak-Voncina et al. (2005)} describe a set of data where seven fatty acid compositions were used to classify commercial oils as either pumpkin (labeled \texttt{A}), sunflower (\texttt{B}), peanut (\texttt{C}), olive (\texttt{D}), soybean (\texttt{E}), rapeseed (\texttt{F}) and corn (\texttt{G}). There were 96 data points contained in their Table 1 with known results. The breakdown of the classes is given in below:

\begin{small}
<<oilData>>=
data(oil)
dim(fattyAcids)
table(oilType)
@
\end{small}

As a note, the paper states on page 32 that there are 37 unknown samples while the table on pages 33 and 34 shows that there are 34 unknowns. 

\subsection{German Credit Data}

Data from Dr. Hans Hofmann of the University of Hamburg and stored at
the \href{http://archive.ics.uci.edu/ml/datasets/Statlog+%28German+Credit+Data%29}{UC Irvine Machine Learning Repository}.

These data have two classes for the credit worthiness: good or bad. There are predictors related to attributes, such as: checking account status, duration, credit history, purpose of the loan, amount of the loan, savings accounts or bonds, employment duration, Installment rate in percentage of disposable income, personal information, other debtors/guarantors, residence duration, property, age, other installment plans, housing, number of existing credits, job information, Number of people being liable to provide maintenance for, telephone, and foreign worker status.

Many of these predictors are discrete and have been expanded into several 0/1 indicator variables

\subsection{Kelly Blue Book resale data for 2005 model year GM cars}

\href{http://www.amstat.org/publications/jse/v16n3/datasets.kuiper.html}{Kuiper (2008)} collected data on Kelly Blue Book resale data for 804 GM cars (2005 model year). 

\texttt{cars} is data frame of the suggested retail price (column Price) and various characteristics of each car (columns Mileage, Cylinder, Doors, Cruise, Sound, Leather, Buick, Cadillac, Chevy, Pontiac, Saab, Saturn, convertible, coupe, hatchback, sedan and wagon)

\subsection{Cell Body Segmentation Data}

\href{http://www.biomedcentral.com/1471-2105/8/340}{Hill, LaPan, Li and Haney (2007)} develop models to predict which cells in a high content screen were well segmented. The data consists of 119 imaging measurements on 2019. The original analysis used 1009 for training and 1010 as a test set (see the column called Case).

The outcome class is contained in a factor variable called Class with levels ``PS'' for poorly segmented and ``WS'' for well segmented.



\section{Creating Dummy Variables}

The function \texttt{dummyVars} can be used to generate a complete (less than full rank parameterized) set of dummy variables from one or more factors. The function takes a formula and a data set and outputs an object that can be used to create the dummy variables using the predict method.

<<etitanic0, results = hide, echo = FALSE>>=
library(earth)
data(etitanic)
@ 
 

For example, the \texttt{etitanic} data set in the \texttt{earth} includes two factors: \texttt{pclass} with levels \Sexpr{paste(levels(etitanic$pclass),  collapse = ", ")} and \texttt{sex} with levels \Sexpr{paste(levels(etitanic$sex),  collapse = ", ")}. The base R function \texttt{model.matrix} would generate the following variables:
<<etitanic1>>=
library(earth)
data(etitanic)
head(model.matrix(survived ~ ., data = etitanic))
@ 
Using \texttt{dummyVars}:'
<<etitanic1>>=
dummies <- dummyVars(survived ~ ., data = etitanic)
head(predict(dummies, newdata = etitanic))
@ 
Note there is no intercept and each factor has a dummy variable for each level, so this parameterization may not be useful for some model functions, such as \texttt{lm}.


\section{Data Pre--Processing Functions}

\texttt{caret} includes several functions to pre--process the predictor data.  Assumes that all of the data are numeric (i.e. factors have been converted to dummy variables via \texttt{model.matrix}, \texttt{dummaryVars}  or other means).

\subsection{Zero-- and Near Zero--Variance Predictors}

In some situations, the data generating mechanism can create predictors that only have a single unique value (i.e. a ``zero--variance predictor"). For many models (excluding tree-based models), this may cause the model to crash or the fit to be unstable. 

Similarly, predictors might have only a handful of unique values that occur with very low frequencies. For example, in the drug resistance data, the \verb+nR11+ descriptor (number of 11-membered rings) data have a few unique numeric values that are highly unbalanced:
 
<<nzvDescr>>=
data.frame(table(mdrrDescr$nR11))
@

The concern here that these predictors may become zero--variance predictors when the data are split into cross--validation/bootstrap sub--samples or that a few samples may have an undue influence on the model. These ``near--zero--variance" predictors may need to be identified and eliminated prior to modeling.


To identify these types of predictors, the following two metrics can be calculated:
\begin{enumerate}
   \item the frequency of the most prevalent value over the second most frequent value (called the ``frequency ratio''), which would be near one for well--behaved predictors and very large for highly--unbalanced data
   \item the ``percent of unique values'' is the number of unique values divided by the total number of samples (times 100) that approaches zero as the granularity of the data increases
\end{enumerate}
If the frequency ratio is less than a pre--specified threshold {\em and} the unique value percentage is less than a threshold, we might consider a predictor to be near zero--variance.

We would not want to falsely identify data that have low granularity but are evenly distributed, such as data from a discrete uniform distribution. Using both criteria should not falsely detect such predictors.

Looking at the MDRR data, the \texttt{nearZeroVar} function can be used to identify near zero--variance variables (the \texttt{saveMetrics} argument can be used to show the details and usually defaults to \texttt{FALSE}):

\begin{small}
<<nzvEx>>=
nzv <- nearZeroVar(mdrrDescr, saveMetrics= TRUE)
nzv[nzv$nzv,][1:10,]
@
\end{small}

<<nzvReduce, results = hide, echo = FALSE>>=
nzv <- nearZeroVar(mdrrDescr)
filteredDescr <- mdrrDescr[, -nzv]
@

By default, \texttt{nearZeroVar(data)} will return the positions of the variables that are flagged to be problematic.

\subsection{Identifying Correlated Predictors}

While there are some models that thrive on correlated predictors (such as \texttt{pls}), other models may benefit from reducing the level of correlation between the predictors. 

Given a correlation matrix, the \texttt{findCorrelation} function uses the following algorithm to flag predictors for removal:

\begin{algorithm}[H]
	\SetAlgoLined
	\RestyleAlgo{boxed}
	\Repeat{no correlations are above a threshold}{
	Find the pair of predictors with the largest absolute correlation;
	
	For both predictors, compute the average correlation between each predictor and all of the other variables\;
	
	Flag the variable with the largest mean absolute correlation for removal\;
	
	Remove this column from the correlation matrix\;
	}
\end{algorithm}

<<tmpCorr, echo = FALSE, results = hide>>=
descrCor <- abs(cor(filteredDescr))
@

For the previous MDRR data, there are \Sexpr{sum(descrCor[upper.tri(descrCor)] > .999)} descriptors that are almost perfectly correlated (|correlation| > 0.999), such as the total information index of atomic composition (\texttt{IAC}) and the total information content index (neighborhood symmetry of 0-order) (\texttt{TIC0}) (correlation = 1). The code chunk below shows the effect of removing descriptors with absolute correlations above 0.75.

\begin{small}
<<corr>>=
descrCor <- cor(filteredDescr)
summary(descrCor[upper.tri(descrCor)])
highlyCorDescr <- findCorrelation(descrCor, cutoff = .75)
filteredDescr <- filteredDescr[,-highlyCorDescr]
descrCor2 <- cor(filteredDescr)
summary(descrCor2[upper.tri(descrCor2)])
@
\end{small}

\subsection{Linear Dependencies}

The function \texttt{findLinearCombos} uses the QR decomposition of a matrix to enumerate sets of linear combinations (if they exist). For example, consider the following matrix that is could have been produced by a less--than--full--rank parameterizations of a two--way experimental layout:

\begin{small}
<<linearCombData>>=
ltfrDesign <- matrix(0, nrow=6, ncol=6)
ltfrDesign[,1] <- c(1, 1, 1, 1, 1, 1)
ltfrDesign[,2] <- c(1, 1, 1, 0, 0, 0)
ltfrDesign[,3] <- c(0, 0, 0, 1, 1, 1)
ltfrDesign[,4] <- c(1, 0, 0, 1, 0, 0)
ltfrDesign[,5] <- c(0, 1, 0, 0, 1, 0)
ltfrDesign[,6] <- c(0, 0, 1, 0, 0, 1)
ltfrDesign
@
\end{small}

Note that columns two and three add up to the first column. Similarly, columns four, five and six add up the first column. \texttt{findLinearCombos} will return a list that enumerates these dependencies. For each linear combination, it will incrementally remove columns from the matrix and test to see if the dependencies have been resolved. \texttt{findLinearCombos} will also return a vector of column positions can be removed to eliminate the linear dependencies:

\newpage
\begin{small}
<<linearComb>>=
comboInfo <- findLinearCombos(ltfrDesign)
comboInfo
ltfrDesign[, -comboInfo$remove]
findLinearCombos(ltfrDesign[, -comboInfo$remove])
@
\end{small}

These types of dependencies can arise when large numbers of binary chemical fingerprints are used to describe the structure of a molecule.


\subsection{Centering and Scaling}

The \texttt{preProcess} class can be used for many operations on predictors, including centering and scaling. The function \texttt{preProcess} estimates the required parameters for each operation and \texttt{predict.preProcess} is used to apply them to specific data sets.


In the example below, the half of the MDRR data are used to estimate the location and scale of the predictors. The function \texttt{preProcess} doesn't actually pre--process the data. \texttt{predict.preProcess} is used to pre--process this and other data sets.
\begin{small}
<<centerScale>>=
set.seed(96)
inTrain <- sample(seq(along = mdrrClass), length(mdrrClass)/2)

training <- filteredDescr[inTrain,]
test <- filteredDescr[-inTrain,]
trainMDRR <- mdrrClass[inTrain]
testMDRR <- mdrrClass[-inTrain]

preProcValues <- preProcess(training, method = c("center", "scale"))

trainDescr <- predict(preProcValues, training)
testDescr <- predict(preProcValues, test)
@
\end{small}

The \texttt{preProcess} option \texttt{"ranges"} scales the data to the interval [0, 1].

\subsection{Imputation}

\texttt{preProcess} can be used to impute data sets based only on
information in the training set. One method of doing this is with
$K$--nearest neighbors. For an arbitrary sample, the $k$ closest
neighbors are found in the training set and the value for the
predictor is imputed using these values (e.g. using the mean). Using
this approach will automatically trigger \texttt{preProcess} to center
and scale the data, regardless of what is in the \texttt{method}
argument. Alternatively, bagged trees can also be used to impute. For
each predictor in the data, a bagged tree is created using all of the
other predictors in the training set. When a new sample has a missing
predictor value, the bagged model is used to predict the value. While,
in theory, this is a more powerful method of imputing, the
computational costs are much higher than the nearest neighbor technique.

\subsection{Transforming Predictors}

In some cases, there is a need to use principal component analysis (PCA) to transform the data to a smaller sub--space where the new variable are uncorrelated with one another. The \texttt{preProcess} class can apply this transformation by including \texttt{"pca"} in the \texttt{method} argument. Doing this will also force scaling of the predictors. Note that when PCA is requested, \texttt{predict.preProcess} changes the column names to \texttt{PC1},  \texttt{PC2} and so on.

Similarly, independent component analysis (ICA) can also be used to
find new variables that are linear combinations of the original set
such that the components are independent (as opposed to uncorrelated in
PCA). The new variables will be labeled as  \texttt{IC1},  \texttt{IC2} and so on.

The ``spatial sign" transformation (\href{http://pubs.acs.org/cgi-bin/abstract.cgi/jcisd8/2006/46/i03/abs/ci050498u.html}{Serneels {\it et al}, 2006}) projects the data for a predictor to the unit circle in $p$ dimensions, where $p$ is the number of predictors. Essentially, a vector of data is divided by its norm. Figures \ref{f:descSplom} and \ref{f:descSplomTrans} show two centered and scaled descriptors from the MDRR data before and after the spatial sign transformation. The predictors should be centered and scaled before applying this transformation.

Another option, \texttt{"BoxCox"} will estimate a Box--Cox transformation on the {\em predictors} if the data are greater than zero.

This transformation can also be applied using the \texttt{preProcess} class.

\begin{figure}[p]
   \begin{center}		
<<splomBefore, echo = FALSE, fig = TRUE>>=
plotSubset <- data.frame(scale(mdrrDescr[, c("nC", "X4v")]))
trellis.par.set(caretTheme(), warn = FALSE)
print(xyplot(nC ~ X4v, plotSubset, groups = mdrrClass, auto.key = list(columns = 2)))
@
      \caption{A plot of two descriptors from the MDRR training set (centered and scaled)}
      \label{f:descSplom}         
   \end{center}
\end{figure}

\begin{figure}[p]
   \begin{center}		
<<splomAfter, echo = FALSE, fig = TRUE>>=
trellis.par.set(caretTheme(), warn = FALSE)
print(xyplot(nC ~ X4v, spatialSign(plotSubset), groups = mdrrClass, auto.key = list(columns = 2)))
@
      \caption{A plot of the two descriptors from the MDRR training set after the spatial sign transformation}
      \label{f:descSplomTrans}         
   \end{center}
\end{figure}


\subsection{Data Splitting}

The function \texttt{createDataPartition} can be used to create balanced splits of the data. If the \texttt{y} argument to this function is a factor, the random sampling occurs within each class and should preserve the overall class distribution of the data. For example, to create a single 80$\%$$\backslash$20$\%$ split of the MDRR data:
\begin{small}
<<mbrrSplit1>>=
set.seed(3456)
set1index <- createDataPartition(mdrrClass, p = .8, list = FALSE, times = 1)
set1 <- mdrrClass[set1index]
round(table(set1)/length(set1), 2)
round(table(mdrrClass)/length(mdrrClass), 2)
@
The \texttt{list = FALSE} returns the data in an integer vector. This function also has an argument, \texttt{times}, that can create multiple splits at once; the data indices are returned in a list of integer vectors. Similarly, \texttt{createResample} can be used to make simple bootstrap samples and \texttt{createFolds} can be used to generate balanced cross--validation groupings from a set of data.
\end{small}

Also, the function \texttt{maxDissim} can be used to create sub--samples using a maximum dissimilarity approach (\href{http://www.liebertonline.com/doi/abs/10.1089/106652799318382}{Willett, 1999}). Suppose there is a data set $A$ with $m$ samples and a larger data set $B$ with $n$ samples. We may want to create a sub--sample from $B$ that is diverse when compared to $A$. To do this, for each sample in $B$, the function calculates the $m$ dissimilarities between each point in $A$. The most dissimilar point in $B$ is added to $A$ and the process continues. There are many methods in R to calculate dissimilarity. \texttt{caret} uses the \texttt{proxy} package. See the manual for that package for a list of available measures. Also, there are many ways to calculate which sample is ``most dissimilar". The argument \texttt{obj} can be used to specify any function that returns a scalar measure. \texttt{caret} includes two functions, \texttt{minDiss} and \texttt{sumDiss}, that can be used to maximize the minimum and total dissimilarities, respectfully.

As an example, Figure \ref{f:diss} shows a scatter plot of two chemical descriptors for the Cox2 data. Using an initial random sample of 5 compounds, we can select 20 more compounds from the data so that the new compounds are most dissimilar from the initial 5 that were specified. The panels in Figure \ref{f:diss} show the results using several combinations of distance metrics and scoring functions. For these data, the distance measure has less of an impact than the scoring method for determining which compounds are most dissimilar.

<<diss, echo = FALSE, results = hide>>=
pdf(paste(getwd(), "/diss.pdf", sep = ""), width = 7, height = 9.5)

data(cox2)

testing <- scale(cox2Descr[, c(102, 30)])
set.seed(1)
startSet <- sample(1:dim(testing)[1], 5)
samplePool <- testing[-startSet,]
start <- testing[startSet,]


par(mfrow=c(3,2))

plot(testing[,1], testing[,2], xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
title("Original Data")

subsetSize <- 20

########################################################

newSamp <- maxDissim(start, samplePool, n = subsetSize)

plot(samplePool[,1], samplePool[,2], type = "n", xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
points(samplePool[-newSamp,1], samplePool[-newSamp,2], col = "grey")

text(samplePool[newSamp,1], samplePool[newSamp,2], paste(1:subsetSize), col = "blue", cex = 1.2)
title("Euclidean Dist with MaxMin")

text(start[,1], start[,2], rep("S", dim(start)[1]), col = "red", cex = 1.2)


########################################################

newSamp <- maxDissim(start, samplePool, n = subsetSize, obj = sumDiss)

plot(samplePool[,1], samplePool[,2], type = "n", xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
points(samplePool[-newSamp,1], samplePool[-newSamp,2], col = "grey")

text(samplePool[newSamp,1], samplePool[newSamp,2], paste(1:subsetSize), col = "blue", cex = 1.2)
title("Euclidean Dist with MaxSum")

text(start[,1], start[,2], rep("S", dim(start)[1]), col = "red", cex = 1.2)

########################################################

newSamp <- maxDissim(start, samplePool, n = subsetSize, obj = sumDiss, randomFrac = 0.05)

plot(samplePool[,1], samplePool[,2], type = "n", xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
points(samplePool[-newSamp,1], samplePool[-newSamp,2], col = "grey")

text(samplePool[newSamp,1], samplePool[newSamp,2], paste(1:subsetSize), col = "blue", cex = 1.2)
title("Euclidean Dist with MaxSum and 5% Sampling")

text(start[,1], start[,2], rep("S", dim(start)[1]), col = "red", cex = 1.2)



########################################################


newSamp <- maxDissim(start, samplePool, n = subsetSize, method = "Manhattan")

plot(samplePool[,1], samplePool[,2], type = "n", xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
points(samplePool[-newSamp,1], samplePool[-newSamp,2], col = "grey")

text(samplePool[newSamp,1], samplePool[newSamp,2], paste(1:subsetSize), col = "blue", cex = 1.2)
title("Manhattan Dist with MaxMin")

text(start[,1], start[,2], rep("S", dim(start)[1]), col = "red", cex = 1.2)


########################################################


newSamp <- maxDissim(start, samplePool, n = subsetSize, method = "Manhattan", obj = sumDiss)

plot(samplePool[,1], samplePool[,2], type = "n", xlab = "Total positive partial charge", ylab = "Predicted Skin Permeability")
points(samplePool[-newSamp,1], samplePool[-newSamp,2], col = "grey")

text(samplePool[newSamp,1], samplePool[newSamp,2], paste(1:subsetSize), col = "blue", cex = 1.2)
title("Manhattan Dist with MaxSum")

text(start[,1], start[,2], rep("S", dim(start)[1]), col = "red", cex = 1.2)


dev.off()
@

\begin{figure}[p]
   \begin{center}	
      \includegraphics[clip, width = .65\textheight]{diss}
      \caption{Examples of selecting 20 dissimilar compounds from an initial set of 5 (denoted by ``S'') for different combinations of arguments to the function \texttt{maxDissim}.}
      \label{f:diss}         
   \end{center}
\end{figure}  

\section{Visualizing Data}

The \texttt{featurePlot} function is a wrapper for different \texttt{lattice} plots to visualize the data.

For example, Figures \ref{f:classEllipse}, \ref{f:densityPlot}, \ref{f:boxPlot} and \ref{f:scatPlot} shows the default plot for continuous outcomes generated using the \texttt{featurePlot} function. 

<<featurePlot, echo = FALSE, results = hide>>=
exDescr <- c("Wap", "G.N..N.", "SEigp")
mdrrPreProc1 <- preProcess(mdrrDescr[,exDescr])

mdrrProc <- predict(mdrrPreProc1, mdrrDescr[,exDescr])

pdf(paste(getwd(), "/fp1.pdf", sep = ""), width = 7, height = 7)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(featurePlot(mdrrProc, mdrrClass, "ellipse"))
dev.off()

pdf(paste(getwd(), "/fp2.pdf", sep = ""), width = 7, height = 3.5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(featurePlot(mdrrProc, mdrrClass, 
      "density", scales = list(x = list(relation="free"), 
      y = list(relation="free")), adjust = 1.5, pch = "|", layout = c(3,1), 
      auto.key = list(columns = 3)))
dev.off()

pdf(paste(getwd(), "/fp3.pdf", sep = ""), width = 7, height = 3.5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(featurePlot(mdrrProc, mdrrClass, "boxplot", 
      scales = list(y = list(relation="free")), layout = c(3,1)))
dev.off()

pdf(paste(getwd(), "/fp4.pdf", sep = ""), width = 7, height = 3.5)
   data(BloodBrain)
   trellis.par.set(caretTheme(), warn = FALSE)
   exDescr <- c("psa_npsa", "clogp", "vsa_acc")   
   print(featurePlot(bbbDescr[,exDescr], logBBB, layout = c(3,1)))
dev.off()
@

\begin{figure}[hb]
   \begin{center}	
      \vspace*{1 in}   	
      \includegraphics[clip, width = .6\textwidth]{fp1}
      \caption{An example plot of three example predictors from the MDRR data using the \texttt{featurePlot} function with density ellipses for each class generated with the arguments
      {\tt plot = "ellipse"} and a factor outcome.}
      \label{f:classEllipse}         
   \end{center}
\end{figure}  

\begin{figure}[p]
   \begin{center}		
      \includegraphics[clip, width = .6\textwidth]{fp2}
      \caption{An example of using the {\tt plot = "density"} argument in the \texttt{featurePlot} function with three example predictors from the MDRR data and a factor outcome.}
      \label{f:densityPlot}
      \vspace*{.25 in}             
      \includegraphics[clip, width = .6\textwidth]{fp3}
      \caption{An example of using the {\tt plot = "boxplot"} argument in the \texttt{featurePlot} function with three example predictors from the MDRR data and a factor outcome.}
      \label{f:boxPlot}      
      \vspace*{.25 in}                   
      \includegraphics[clip, width = .6\textwidth]{fp4}
      \caption{A set of scatter--plots for three example predictors from the blood--brain barrier data and a numeric outcome (on the y--axis).}            
      \label{f:scatPlot}         
   \end{center}
\end{figure}  

\clearpage 

\section{Other Functions}

\subsection{Processing Affy Arrays}

For Affymetrix gene chip data, RMA processing (\href{http://nar.oxfordjournals.org/cgi/content/abstract/31/4/e15}{Irizarry, 2003}) is a popular method of processing gene expression data. However, for predictive modeling, it has a drawback in that the processing is batch oriented; if an additional sample is collected, the RMA processing must be repeated using all the samples. This is mainly because of two steps in the processing: the quantile normalization process and the calculation of expression. Quantile normalization normalizes the data such that the between--chip distributions have equivalent quantiles and is very effective in improving the quality of the data. It is possible to let the samples in the training set define the reference distribution for normalization. For the expression calculation, a robust method such as a trimmed mean can be used to summarize the probe level data into a single summary metric per sample. For example,
\begin{Verbatim}
# first, let affy/expresso know that the method exists
normalize.AffyBatch.methods <- c(normalize.AffyBatch.methods, "normalize2Reference")
   
RawData <- ReadAffy(celfile.path=FilePath)
   
Batch1Step1 <- bg.correct(RawData, "rma")
Batch1Step2 <- normalize.AffyBatch.quantiles(Batch1Step1)
referencePM <- pm(Batch1Step2)[,1]
Batch1Step3 <- computeExprSet(Batch1Step2, "pmonly", "trimMean")  
   
Batch2Step1 <- bg.correct(RawData2, "rma")
Batch2Step2 <- normalize.AffyBatch.normalize2Reference(Batch2Step1, ref = referencePM)
Batch2Step3 <- computeExprSet(Batch2Step2, "pmonly", "trimMean")     
\end{Verbatim}


\subsection{Yet Another $k$--Nearest Neighbor Function}

\texttt{knn3} is a function for $k$--nearest neighbor classification. This particular implementation is a modification of the \texttt{knn} C code and returns the vote information for all of the classes (\texttt{knn} only returns the probability for the winning class). There is a formula interface via
\begin{Verbatim}
knn3(formula, data)
\end{Verbatim}
or by passing the training data directly
\begin{Verbatim}
# x is a matrix, y is a factor vector
knn3(x, y)
\end{Verbatim}
There are also \texttt{print} and \texttt{predict} methods.

Similarly, \texttt{caret} contains a $k$--nearest neighbor regression function, \texttt{knnreg}. It returns the average outcome for the neighbor.


\subsection{Partial Least Squares Discriminant Analysis}

The \texttt{plsda} function is a wrapper for the \texttt{plsr} function in the \texttt{pls} package that does not require a formula interface and can take factor outcomes as arguments. The classes are broken down into dummy variables (one for each class). These 0/1 dummy variables are modeled by partial least squares. 

From this model, there are two approaches to computing the class predictions and probabilities:
\begin{itemize}
\item the softmax technique can be used on a per-sample basis to normalize the scores so that they are more ``probability like'' (i.e. they sum to one and are between zero and one). For a vector of model predictions for each class $x$, the softmax class probabilities are computed as
\[
prob_k= \frac{exp(x_k)}{\sum_{k=1}^C exp(x_k)}
\]
The predicted class is simply the class with the largest model prediction, or equivalently, the largest class probability. This is the default behavior for \texttt{plsda}.
\item Bayes rule can be applied to the model predictions to form posterior probabilities. Here, the model predictions for the training set are used along with the training set outcomes to create conditional distributions for each class. When new samples are predicted, the raw model predictions are run through these conditional distributions to produce a posterior probability for each class (along with the prior). Bayes rule can be used by specifying \texttt{probModel = "Bayes"}. An additional parameter, \texttt{prior}, can be used to set prior probabilities  for the classes.
\end{itemize}
The advantage to using Bayes rule is that the full training set is used to directly compute the class probabilities (unlike the softmax function which only uses the current sample's scores). This creates more realistic probability estimates (See Figure \ref{f:plsProbs}). The disadvantage is that a separate Bayesian model must be created for each value of \texttt{ncomp}, which is more time consuming.  
 
For the MDRR data set, we can fit two PLS models using each technique and predict the class probabilities for the test set. Figure \ref{f:plsProbs}  shows the results.

Similar to \texttt{plsda}, \texttt{caret} also contains a function \texttt{splsda} that allows for classification using sparse PLS. A dummy 
matrix is created for each class and used with the \texttt{spls} function in the \texttt{spls} package. The same approach to estimating class probabilities is used for \texttt{plsda} and \texttt{splsda}.

<<plsProb, echo = FALSE, results = hide>>=
bTime <- system.time(Bayes <- plsda(trainDescr, trainMDRR, ncomp = 5, probMethod = "Bayes"))[3]
sTime <- system.time(softmax <- plsda(trainDescr, trainMDRR, ncomp = 5))[3]

bProbs <- as.data.frame(predict(Bayes, as.matrix(testDescr), type = "prob"))
colnames(bProbs) <- gsub(".ncomp5", "", colnames(bProbs), fixed = TRUE)
bProbs$method <- "Bayes Rule"
bProbs$obs <- testMDRR
sProbs <- as.data.frame(predict(softmax, as.matrix(testDescr), type = "prob"))
colnames(bProbs) <- gsub(".5 comps", "", colnames(bProbs), fixed = TRUE)
sProbs$method <- "Softmax"
sProbs$obs <- testMDRR
names(sProbs) <- names(bProbs)
both <- rbind(bProbs, sProbs)

bPerf <- confusionMatrix(predict(Bayes, as.matrix(testDescr)),
                         testMDRR)
sPerf <- confusionMatrix(predict(softmax, as.matrix(testDescr)),
                         testMDRR)
                                 
pdf(paste(getwd(), "/plsProbs.pdf", sep = ""), width = 7, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
print(
      histogram(
                ~Active|method*obs, data = both, 
                xlab = "Active Probability",
                as.table = TRUE))
dev.off()
@ 

\begin{figure}
   \begin{center}		
      \includegraphics[clip, width = .9\textwidth]{plsProbs}
\caption{A comparison of two methods for computing class probabilities for PLS classification models on the MDRR data. Using a 5 component PLS model, the test set sample class probabilities were predicted using either the softmax function or Bayes rule. The histograms show the distributions of the active class probabilities for both methods and the true class of the samples. While Bayes rule takes more time to compute (\Sexpr{round(bTime, 2)} seconds compared to \Sexpr{round(sTime, 3)} seconds for softmax), the probabilities using Bayes rule appear more appropriate since they encompass the full scale. The softmax probabilities are close to the middle of the distribution and rarely show high confidence in the predicted class. Both techniques show approximately the same accuracy (\Sexpr{round(bPerf$overall[1], 2)} for Bayes and \Sexpr{round(sPerf$overall[1], 2)} for softmax) and Kappa values (\Sexpr{round(bPerf$overall[2], 2)} for Bayes and \Sexpr{round(sPerf$overall[2], 2)} for softmax).}
      \label{f:plsProbs}       
   \end{center}
\end{figure}  

\subsection{Bagged MARS and FDA}

Multivariate adaptive regression splines (MARS) models, like classification/regression trees, are unstable predictors (Breiman, 1996). This means that small perturbations in the training data might lead to significantly different models. Bagged trees and random forests are effective ways of improving tree models by exploiting these instabilities. \texttt{caret} contains a function, \texttt{bagEarth}, that fits MARS models via the \texttt{earth} function. There are formula and non--formula interfaces. 

Also, flexible discriminant analysis is a generalization of linear discriminant analysis that can use non--linear features as inputs. One way of doing this is the use MARS--type features to classify samples. The function \texttt{bagFDA} fits FDA models of a set of bootstrap samples and aggregates the predictions to reduce noise.  

This function is deprecated in favor of the \texttt{bag} function.

\subsection{Neural Networks with a Principal Component Step}

Neural networks can be affected by severe amounts of multicollinearity
in the predictors. The function \texttt{pcaNNet} is a wrapper around
the \texttt{preProcess} and \texttt{nnet} functions that will run
principal component analysis on the predictors before using them as
inputs into a neural network. The function will keep enough components
that will capture some pre-defined threshold on the cumulative
proportion of variance (see the \texttt{thresh}
argument). For new samples, the same transformation is applied to the new predictor values (based on the loadings from the training set). The function is available for both regression and classification. 

This function is deprecated in favor of the \texttt{train} function using \texttt{method = "nnet"} and \texttt{preProc = "pca"}.

\subsection{Class Distance Calculations}

\texttt{caret} contains functions to generate new predictors variables based on distances to class centroids (similar to how linear discriminant analysis works). For each level of a factor variable, the class centroid and covariance matrix is calculated. For new samples, the Mahalanobis distance to each of the class centroids is computed and can be used as an additional predictor. This can be helpful for non--linear models when the true decision boundary is actually linear.

In cases where there are more predictors within a class than samples, the \texttt{classDist} function has arguments called \texttt{pca} and \texttt{keep} arguments that allow for principal components analysis within each class to be used to avoid issues with singular covariance matrices. 

\texttt{predict.classDist} is then used to generate the class distances. By default, the distances are logged, but this can be changed via the \texttt{trans} argument to \texttt{predict.classDist}.

As an example, we can used the Fisher/Anderson iris data. 
<<classDist>>=
inTrain <- seq(1, 150, 6)
centroids <- classDist(iris[inTrain, 1:4], iris$Species[inTrain])
distances <- predict(centroids, iris[-inTrain, 1:4])
head(distances)
@ 
Figure \ref{F:classDist} shows a scatterplot matrix of the class distances for the held--out samples. 

For data sets with numerical outcomes, the data can be split into groups based on quantiles of the data (more fine granularity can be achieved with the \texttt{groups} argument to the \texttt{classDist} function).

\begin{figure}
   \begin{center}		
<<classDistPlot, echo = FALSE, results = hide>>=
pdf("classDist.pdf", width = 5, height = 7)
   trellis.par.set(caretTheme())
   print(
         splom(distances, 
               groups = iris$Species[-inTrain],
               auto.key = list(columns = 3)))
dev.off()
@
      \includegraphics[clip]{classDist}
      \caption{Class distances for the help--out samples for the Fisher/Anderson iris computed using the \texttt{classDist} function. The distances have been logged.}
      \label{F:classDist} 
    \end{center}
\end{figure} 

For example, in the Fisher/Anderson iris data there are four predictors and a class variable describing three species of iris. These data are typically used to demonstrate how the four predictors can be used to predict the three species. 


\subsection{Independent Component Regression}

The \texttt{icr} function can be used to fit a model analogous to
principal component regression (PCR), but using independent component
analysis (ICA). The predictor data are centered and projected to the
ICA components. These components are then regressed against the
outcome. The user needed to specify the number of components to keep.

The model uses the \texttt{preProcess} function to compute the latent
variables using the  {\tt \href{http://cran.r-project.org/web/packages/fastICA/index.html}{fastICA}} package.

Like PCR, there is no guarantee that there will be a correlation
between the new latent variable and the outcomes.

\section{Session Information}

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

\section{References}

\begin{description}
   \item Blake, J. F. (2000), ``Chemoinformatics: predicting the physicochemical properties of 'drug--like' molecules,'' {\it Current Opinion in Biotechnology}, 11, 104--107.

   \item Borggaard, C., Thodberg, H. H. (1992), ``Optimal minimal neural interpretation of spectra,''  {\it Analytical Chemistry}, 64, 545--551.

   \item Breiman, L. (1996), ``Bagging predictors,'' {\it Machine Learning}, 24, 123--140.

	\item  Brodnjak-Voncina et al. (2005). ``Multivariate data analysis in classification of vegetable oils characterized
	by the content of fatty acids,'' {\it Chemometrics and Intelligent Laboratory Systems}, 75, 31--45.

   \item Clark, D. E., Pickett, S. D. (2000), ``Computational methods for the prediction of 'drug--likeness','' {\it Drug Discov Today}, 5, 49--58.

     
     \item Hill, LaPan, Li and Haney (2007). ``Impact of image
       segmentation on high-content screening data quality for SK-BR-3
       cells'', {\it BMC Bioinformatics}, Vol. 8, pg. 340
       
	\item Irizarry, RA, Bolstad BM, Collin, F, Cope, LM, Hobbs, B, and, Speed, TP (2003), ``Summaries of Affymetrix Genechip probe level data,''  {\it Nucleic Acids Research}, 31, e15

        \item Kuiper, S. (2008). ``Introduction to Multiple Regression:
          How Much Is Your Car Worth?'', {\it Journal of Statistics Education}, Vol. 16,  
          
   \item Mente, S. R., Lombardo, F. (2005), ``A recursive-partitioning model for blood--brain barrier permeation,'' {\it Journal of Computer--Aided Molecular Design}, 19, 465--481.


   \item Serneels, S., De Nolf, E. and Van Espen, P. J. (2006), ``Spatial sign preprocessing: a simple way to impart moderate robustness to multivariate estimators,'' {\it Journal of Chemical Information and Modeling}, 46, 1402--1409.

   \item Sutherland, J.J., O'Brien L.A. and Weaver D.F. (2004), ``A comparison of methods for modeling quantitative structure-activity relationships,'' {\it Journal of Medicinal Chemistry}, 47, 5541--5554.
     
     \item Sutherland, J.J. and Weaver,
       D.F. (2004). ``Three-dimensional quantitative
       structure-activity and structure-selectivity relationships of
       dihydrofolate reductase inhibitors'',{\it  Journal of Computer-Aided Molecular Design}, Vol. 18, pg. 309--331.

   \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 Willett, P. (1999), ``Dissimilarity-Based Algorithms for Selecting Structurally Diverse Sets of Compounds,'' {\it Journal of Computational Biology}, 6, 447--457.

\end{description}

\end{document}


      
