|
0
|
1 def __template4Rnw():
|
|
|
2
|
|
|
3 template4Rnw = r'''%% Classification Modeling Script
|
|
|
4 %% Max Kuhn (max.kuhn@pfizer.com, mxkuhn@gmail.com)
|
|
|
5 %% Version: 1.00
|
|
|
6 %% Created on: 2010/10/02
|
|
|
7 %%
|
|
|
8 %% This is an Sweave template for building and describing
|
|
|
9 %% classification models. It mixes R and LaTeX code. The document can
|
|
|
10 %% be processing using R's Sweave function to produce a tex file.
|
|
|
11 %%
|
|
|
12 %% The inputs are:
|
|
|
13 %% - the initial data set in a data frame called 'rawData'
|
|
|
14 %% - a factor column in the data set called 'class'. this should be the
|
|
|
15 %% outcome variable
|
|
|
16 %% - all other columns in rawData should be predictor variables
|
|
|
17 %% - the type of model should be in a variable called 'modName'.
|
|
|
18 %%
|
|
|
19 %% The script attempts to make some intelligent choices based on the
|
|
|
20 %% model being used. For example, if modName is "pls", the script will
|
|
|
21 %% automatically center and scale the predictor data. There are
|
|
|
22 %% situations where these choices can (and should be) changed.
|
|
|
23 %%
|
|
|
24 %% There are other options that may make sense to change. For example,
|
|
|
25 %% the user may want to adjust the type of resampling. To find these
|
|
|
26 %% parts of the script, search on the string 'OPTION'. These parts of
|
|
|
27 %% the code will document the options.
|
|
|
28
|
|
|
29 \documentclass[14pt]{report}
|
|
|
30 \usepackage{amsmath}
|
|
|
31 \usepackage[pdftex]{graphicx}
|
|
|
32 \usepackage{color}
|
|
|
33 \usepackage{ctable}
|
|
|
34 \usepackage{xspace}
|
|
|
35 \usepackage{fancyvrb}
|
|
|
36 \usepackage{fancyhdr}
|
|
|
37 \usepackage{lastpage}
|
|
|
38 \usepackage{longtable}
|
|
|
39 \usepackage{algorithm2e}
|
|
|
40 \usepackage[
|
|
|
41 colorlinks=true,
|
|
|
42 linkcolor=blue,
|
|
|
43 citecolor=blue,
|
|
|
44 urlcolor=blue]
|
|
|
45 {hyperref}
|
|
|
46 \usepackage{lscape}
|
|
|
47 \usepackage{Sweave}
|
|
|
48 \SweaveOpts{keep.source = TRUE}
|
|
|
49
|
|
|
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
51
|
|
|
52 % define new colors for use
|
|
|
53 \definecolor{darkgreen}{rgb}{0,0.6,0}
|
|
|
54 \definecolor{darkred}{rgb}{0.6,0.0,0}
|
|
|
55 \definecolor{lightbrown}{rgb}{1,0.9,0.8}
|
|
|
56 \definecolor{brown}{rgb}{0.6,0.3,0.3}
|
|
|
57 \definecolor{darkblue}{rgb}{0,0,0.8}
|
|
|
58 \definecolor{darkmagenta}{rgb}{0.5,0,0.5}
|
|
|
59
|
|
|
60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
61
|
|
|
62 \newcommand{\bld}[1]{\mbox{\boldmath $$#1$$}}
|
|
|
63 \newcommand{\shell}[1]{\mbox{$$#1$$}}
|
|
|
64 \renewcommand{\vec}[1]{\mbox{\bf {#1}}}
|
|
|
65
|
|
|
66 \newcommand{\ReallySmallSpacing}{\renewcommand{\baselinestretch}{.6}\Large\normalsize}
|
|
|
67 \newcommand{\SmallSpacing}{\renewcommand{\baselinestretch}{1.1}\Large\normalsize}
|
|
|
68
|
|
|
69 \newcommand{\halfs}{\frac{1}{2}}
|
|
|
70
|
|
|
71 \setlength{\oddsidemargin}{-.25 truein}
|
|
|
72 \setlength{\evensidemargin}{0truein}
|
|
|
73 \setlength{\topmargin}{-0.2truein}
|
|
|
74 \setlength{\textwidth}{7 truein}
|
|
|
75 \setlength{\textheight}{8.5 truein}
|
|
|
76 \setlength{\parindent}{0.20truein}
|
|
|
77 \setlength{\parskip}{0.10truein}
|
|
|
78
|
|
|
79 \setcounter{LTchunksize}{50}
|
|
|
80
|
|
|
81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
82 \pagestyle{fancy}
|
|
|
83 \lhead{}
|
|
|
84 %% OPTION Report header name
|
|
|
85 \chead{Classification Model Script}
|
|
|
86 \rhead{}
|
|
|
87 \lfoot{}
|
|
|
88 \cfoot{}
|
|
|
89 \rfoot{\thepage\ of \pageref{LastPage}}
|
|
|
90 \renewcommand{\headrulewidth}{1pt}
|
|
|
91 \renewcommand{\footrulewidth}{1pt}
|
|
|
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
93
|
|
|
94 %% OPTION Report title and modeler name
|
|
|
95 \title{Classification Model Script using $METHOD}
|
|
|
96 \author{"Lynn Group with M. Kuhn, SCIS, JNU, New Delhi"}
|
|
|
97
|
|
|
98 \begin{document}
|
|
|
99
|
|
|
100 \maketitle
|
|
|
101
|
|
|
102 \thispagestyle{empty}
|
|
|
103 <<dummy, eval=TRUE, echo=FALSE, results=hide>>=
|
|
|
104 # sets values for variables used later in the program to prevent the \Sexpr error on parsing with Sweave
|
|
|
105 numSamples=''
|
|
|
106 classDistString=''
|
|
|
107 missingText=''
|
|
|
108 numPredictors=''
|
|
|
109 numPCAcomp=''
|
|
|
110 pcaText=''
|
|
|
111 nzvText=''
|
|
|
112 corrText=''
|
|
|
113 ppText=''
|
|
|
114 varText=''
|
|
|
115 splitText="Dummy Text"
|
|
|
116 nirText="Dummy Text"
|
|
|
117 # pctTrain is a variable that is initialised in Data splitting, and reused later in testPred
|
|
|
118 pctTrain=0.8
|
|
|
119 Smpling=''
|
|
|
120 nzvText1=''
|
|
|
121 classDistString1=''
|
|
|
122 dwnsmpl=''
|
|
|
123 upsmpl=''
|
|
|
124
|
|
|
125 @
|
|
|
126 <<startup, eval= TRUE, results = hide, echo = FALSE>>=
|
|
|
127 library(Hmisc)
|
|
|
128 library(caret)
|
|
|
129 library(pROC)
|
|
|
130 versionTest <- compareVersion(packageDescription("caret")$$Version,
|
|
|
131 "4.65")
|
|
|
132 if(versionTest < 0) stop("caret version 4.65 or later is required")
|
|
|
133
|
|
|
134 library(RColorBrewer)
|
|
|
135
|
|
|
136
|
|
|
137 listString <- function (x, period = FALSE, verbose = FALSE)
|
|
|
138 {
|
|
|
139 if (verbose) cat("\n entering listString\n")
|
|
|
140 flush.console()
|
|
|
141 if (!is.character(x))
|
|
|
142 x <- as.character(x)
|
|
|
143 numElements <- length(x)
|
|
|
144 out <- if (length(x) > 0) {
|
|
|
145 switch(min(numElements, 3), x, paste(x, collapse = " and "),
|
|
|
146 {
|
|
|
147 x <- paste(x, c(rep(",", numElements - 2), " and", ""), sep = "")
|
|
|
148 paste(x, collapse = " ")
|
|
|
149 })
|
|
|
150 }
|
|
|
151 else ""
|
|
|
152 if (period) out <- paste(out, ".", sep = "")
|
|
|
153 if (verbose) cat(" leaving listString\n\n")
|
|
|
154 flush.console()
|
|
|
155 out
|
|
|
156 }
|
|
|
157
|
|
|
158 resampleStats <- function(x, digits = 3)
|
|
|
159 {
|
|
|
160 bestPerf <- x$$bestTune
|
|
|
161 colnames(bestPerf) <- gsub("^\\.", "", colnames(bestPerf))
|
|
|
162 out <- merge(x$$results, bestPerf)
|
|
|
163 out <- out[, colnames(out) %in% x$$perfNames]
|
|
|
164 names(out) <- gsub("ROC", "area under the ROC curve", names(out), fixed = TRUE)
|
|
|
165 names(out) <- gsub("Sens", "sensitivity", names(out), fixed = TRUE)
|
|
|
166 names(out) <- gsub("Spec", "specificity", names(out), fixed = TRUE)
|
|
|
167 names(out) <- gsub("Accuracy", "overall accuracy", names(out), fixed = TRUE)
|
|
|
168 names(out) <- gsub("Kappa", "Kappa statistics", names(out), fixed = TRUE)
|
|
|
169
|
|
|
170 out <- format(out, digits = digits)
|
|
|
171 listString(paste(names(out), "was", out))
|
|
|
172 }
|
|
|
173
|
|
|
174 twoClassNoProbs <- function (data, lev = NULL, model = NULL)
|
|
|
175 {
|
|
|
176 out <- c(sensitivity(data[, "pred"], data[, "obs"], lev[1]),
|
|
|
177 specificity(data[, "pred"], data[, "obs"], lev[2]),
|
|
|
178 confusionMatrix(data[, "pred"], data[, "obs"])$$overall["Kappa"])
|
|
|
179
|
|
|
180 names(out) <- c("Sens", "Spec", "Kappa")
|
|
|
181 out
|
|
|
182 }
|
|
|
183
|
|
|
184
|
|
|
185
|
|
|
186 ##OPTION: model name: see ?train for more values/models
|
|
|
187 modName <- "$METHOD"
|
|
|
188
|
|
|
189
|
|
|
190 load("$RDATA")
|
|
|
191 rawData <- dataX
|
|
|
192 rawData$$outcome <- dataY
|
|
|
193
|
|
|
194 @
|
|
|
195
|
|
|
196
|
|
|
197 \section*{Data Sets}\label{S:data}
|
|
|
198
|
|
|
199 %% OPTION: provide some background on the problem, the experimental
|
|
|
200 %% data, how the compounds were selected etc
|
|
|
201
|
|
|
202 <<getDataInfo, eval = $GETDATAINFOEVAL, echo = $GETDATAINFOECHO, results = $GETDATAINFORESULT>>=
|
|
|
203 if(!any(names(rawData) == "outcome")) stop("a variable called outcome should be in the data set")
|
|
|
204 if(!is.factor(rawData$$outcome)) stop("the outcome should be a factor vector")
|
|
|
205
|
|
|
206 ## OPTION: when there are only two classes, the first level of the
|
|
|
207 ## factor is used as the "positive" or "event" for calculating
|
|
|
208 ## sensitivity and specificity. Adjust the outcome factor accordingly.
|
|
|
209 numClasses <- length(levels(rawData$$outcome))
|
|
|
210 numSamples <- nrow(rawData)
|
|
|
211 numPredictors <- ncol(rawData) - 1
|
|
|
212 predictorNames <- names(rawData)[names(rawData) != "outcome"]
|
|
|
213
|
|
|
214 isNum <- apply(rawData[,predictorNames, drop = FALSE], 2, is.numeric)
|
|
|
215 if(any(!isNum)) stop("all predictors in rawData should be numeric")
|
|
|
216
|
|
|
217 classTextCheck <- all.equal(levels(rawData$$outcome), make.names(levels(rawData$$outcome)))
|
|
|
218 if(!classTextCheck) warning("the class levels are not valid R variable names; this may cause errors")
|
|
|
219
|
|
|
220 ## Get the class distribution
|
|
|
221 classDist <- table(rawData$$outcome)
|
|
|
222 classDistString <- paste("``",
|
|
|
223 names(classDist),
|
|
|
224 "'' ($$n$$=",
|
|
|
225 classDist,
|
|
|
226 ")",
|
|
|
227 sep = "")
|
|
|
228 classDistString <- listString(classDistString)
|
|
|
229 @
|
|
|
230
|
|
|
231 <<missingFilter, eval = $MISSINGFILTEREVAL, echo = $MISSINGFILTERECHO, results = $MISSINGFILTERRESULT>>=
|
|
|
232 colRate <- apply(rawData[, predictorNames, drop = FALSE],
|
|
|
233 2, function(x) mean(is.na(x)))
|
|
|
234
|
|
|
235 ##OPTION thresholds can be changed
|
|
|
236 colExclude <- colRate > $MISSINGFILTERTHRESHC
|
|
|
237
|
|
|
238 missingText <- ""
|
|
|
239
|
|
|
240 if(any(colExclude))
|
|
|
241 {
|
|
|
242 missingText <- paste(missingText,
|
|
|
243 ifelse(sum(colExclude) > 1,
|
|
|
244 " There were ",
|
|
|
245 " There was "),
|
|
|
246 sum(colExclude),
|
|
|
247 ifelse(sum(colExclude) > 1,
|
|
|
248 " predictors ",
|
|
|
249 " predictor "),
|
|
|
250 "with an excessive number of ",
|
|
|
251 "missing data. ",
|
|
|
252 ifelse(sum(colExclude) > 1,
|
|
|
253 " These were excluded. ",
|
|
|
254 " This was excluded. "))
|
|
|
255 predictorNames <- predictorNames[!colExclude]
|
|
|
256 rawData <- rawData[, names(rawData) %in% c("outcome", predictorNames), drop = FALSE]
|
|
|
257 }
|
|
|
258
|
|
|
259
|
|
|
260 rowRate <- apply(rawData[, predictorNames, drop = FALSE],
|
|
|
261 1, function(x) mean(is.na(x)))
|
|
|
262
|
|
|
263 rowExclude <- rowRate > $MISSINGFILTERTHRESHR
|
|
|
264
|
|
|
265
|
|
|
266 if(any(rowExclude)) {
|
|
|
267 missingText <- paste(missingText,
|
|
|
268 ifelse(sum(rowExclude) > 1,
|
|
|
269 " There were ",
|
|
|
270 " There was "),
|
|
|
271 sum(colExclude),
|
|
|
272 ifelse(sum(rowExclude) > 1,
|
|
|
273 " samples ",
|
|
|
274 " sample "),
|
|
|
275 "with an excessive number of ",
|
|
|
276 "missing data. ",
|
|
|
277 ifelse(sum(rowExclude) > 1,
|
|
|
278 " These were excluded. ",
|
|
|
279 " This was excluded. "),
|
|
|
280 "After filtering, ",
|
|
|
281 sum(!rowExclude),
|
|
|
282 " samples remained.")
|
|
|
283 rawData <- rawData[!rowExclude, ]
|
|
|
284 hasMissing <- apply(rawData[, predictorNames, drop = FALSE],
|
|
|
285 1, function(x) mean(is.na(x)))
|
|
|
286 } else {
|
|
|
287 hasMissing <- apply(rawData[, predictorNames, drop = FALSE],
|
|
|
288 1, function(x) any(is.na(x)))
|
|
|
289 missingText <- paste(missingText,
|
|
|
290 ifelse(missingText == "",
|
|
|
291 "There ",
|
|
|
292 "Subsequently, there "),
|
|
|
293 ifelse(sum(hasMissing) == 1,
|
|
|
294 "was ",
|
|
|
295 "were "),
|
|
|
296 ifelse(sum(hasMissing) > 0,
|
|
|
297 sum(hasMissing),
|
|
|
298 "no"),
|
|
|
299 ifelse(sum(hasMissing) == 1,
|
|
|
300 "sample ",
|
|
|
301 "samples "),
|
|
|
302 "with missing values.")
|
|
|
303
|
|
|
304 rawData <- rawData[complete.cases(rawData),]
|
|
|
305
|
|
|
306 }
|
|
|
307
|
|
|
308 rawData1 <- rawData[,1:length(rawData)-1]
|
|
|
309 rawData2 <- rawData[,length(rawData)]
|
|
|
310
|
|
|
311 set.seed(222)
|
|
|
312 nzv1 <- nearZeroVar(rawData1)
|
|
|
313 if(length(nzv1) > 0)
|
|
|
314 {
|
|
|
315 nzvVars1 <- names(rawData1)[nzv1]
|
|
|
316 rawData <- rawData1[, -nzv1]
|
|
|
317 rawData$outcome <- rawData2
|
|
|
318 nzvText1 <- paste("There were ",
|
|
|
319 length(nzv1),
|
|
|
320 " predictors that were removed from original data due to",
|
|
|
321 " severely unbalanced distributions that",
|
|
|
322 " could negatively affect the model fit",
|
|
|
323 ifelse(length(nzv1) > 10,
|
|
|
324 ".",
|
|
|
325 paste(": ",
|
|
|
326 listString(nzvVars1),
|
|
|
327 ".",
|
|
|
328 sep = "")),
|
|
|
329 sep = "")
|
|
|
330
|
|
|
331 } else {
|
|
|
332 rawData <- rawData1
|
|
|
333 rawData$outcome <- rawData2
|
|
|
334 nzvText1 <- ""
|
|
|
335
|
|
|
336 }
|
|
|
337
|
|
|
338 remove("rawData1")
|
|
|
339 remove("rawData2")
|
|
|
340
|
|
|
341 @
|
|
|
342
|
|
|
343 The initial data set consisted of \Sexpr{numSamples} samples and
|
|
|
344 \Sexpr{numPredictors} predictor variables. The breakdown of the
|
|
|
345 outcome data classes were: \Sexpr{classDistString}.
|
|
|
346
|
|
|
347 \Sexpr{missingText}
|
|
|
348
|
|
|
349 \Sexpr{nzvText1}
|
|
|
350
|
|
|
351 <<pca, eval= $PCAEVAL, echo = $PCAECHO, results = $PCARESULT>>=
|
|
|
352
|
|
|
353 predictorNames <- names(rawData)[names(rawData) != "outcome"]
|
|
|
354 numPredictors <- length(predictorNames)
|
|
|
355 predictors <- rawData[, predictorNames, drop = FALSE]
|
|
|
356 ## PCA will fail with predictors having less than 2 unique values
|
|
|
357 isZeroVar <- apply(predictors, 2,
|
|
|
358 function(x) length(unique(x)) < 2)
|
|
|
359 if(any(isZeroVar)) predictors <- predictors[, !isZeroVar, drop = FALSE]
|
|
|
360 ## For whatever, only the formula interface to prcomp
|
|
|
361 ## handles missing values
|
|
|
362 pcaForm <- as.formula(
|
|
|
363 paste("~",
|
|
|
364 paste(names(predictors), collapse = "+")))
|
|
|
365 pca <- prcomp(pcaForm,
|
|
|
366 data = predictors,
|
|
|
367 center = TRUE,
|
|
|
368 scale. = TRUE,
|
|
|
369 na.action = na.omit)
|
|
|
370 ## OPTION: the number of components plotted/discussed can be set
|
|
|
371 numPCAcomp <- $PCACOMP
|
|
|
372 pctVar <- pca$$sdev^2/sum(pca$$sdev^2)*100
|
|
|
373 pcaText <- paste(round(pctVar[1:numPCAcomp], 1),
|
|
|
374 "\\\\%",
|
|
|
375 sep = "")
|
|
|
376 pcaText <- listString(pcaText)
|
|
|
377 @
|
|
|
378
|
|
|
379 To get an initial assessment of the separability of the classes,
|
|
|
380 principal component analysis (PCA) was used to distill the
|
|
|
381 \Sexpr{numPredictors} predictors down into \Sexpr{numPCAcomp}
|
|
|
382 surrogate variables (i.e. the principal components) in a manner that
|
|
|
383 attempts to maximize the amount of information preserved from the
|
|
|
384 original predictor set. Figure \ref{F:inititalPCA} contains plots of
|
|
|
385 the first \Sexpr{numPCAcomp} components, which accounted for
|
|
|
386 \Sexpr{pcaText} percent of the variability in the original predictors
|
|
|
387 (respectively).
|
|
|
388
|
|
|
389
|
|
|
390 %% OPTION: remark on how well (or poorly) the data separated
|
|
|
391
|
|
|
392 \setkeys{Gin}{width = 0.8\textwidth}
|
|
|
393 \begin{figure}[p]
|
|
|
394 \begin{center}
|
|
|
395
|
|
|
396 <<pcaPlot, eval = $PCAPLOTEVAL, echo = $PCAPLOTECHO, results = $PCAPLOTRESULT, fig = $PCAPLOTFIG, width = 8, height = 8>>=
|
|
|
397 trellis.par.set(caretTheme(), warn = TRUE)
|
|
|
398 if(numPCAcomp == 2)
|
|
|
399 {
|
|
|
400 axisRange <- extendrange(pca$$x[, 1:2])
|
|
|
401 print(
|
|
|
402 xyplot(PC1 ~ PC2,
|
|
|
403 data = as.data.frame(pca$$x),
|
|
|
404 type = c("p", "g"),
|
|
|
405 groups = rawData$$outcome,
|
|
|
406 auto.key = list(columns = 2),
|
|
|
407 xlim = axisRange,
|
|
|
408 ylim = axisRange))
|
|
|
409 } else {
|
|
|
410 axisRange <- extendrange(pca$$x[, 1:numPCAcomp])
|
|
|
411 print(
|
|
|
412 splom(~as.data.frame(pca$$x)[, 1:numPCAcomp],
|
|
|
413 type = c("p", "g"),
|
|
|
414 groups = rawData$$outcome,
|
|
|
415 auto.key = list(columns = 2),
|
|
|
416 as.table = TRUE,
|
|
|
417 prepanel.limits = function(x) axisRange
|
|
|
418 ))
|
|
|
419
|
|
|
420 }
|
|
|
421
|
|
|
422 @
|
|
|
423
|
|
|
424 \caption[PCA Plot]{A plot of the first \Sexpr{numPCAcomp}
|
|
|
425 principal components for the original data set.}
|
|
|
426 \label{F:inititalPCA}
|
|
|
427 \end{center}
|
|
|
428 \end{figure}
|
|
|
429
|
|
|
430
|
|
|
431
|
|
|
432 <<initialDataSplit, eval = $INITIALDATASPLITEVAL, echo = $INITIALDATASPLITECHO, results = $INITIALDATASPLITRESULT>>=
|
|
|
433
|
|
|
434 ## OPTION: in small samples sizes, you may not want to set aside a
|
|
|
435 ## training set and focus on the resampling results.
|
|
|
436
|
|
|
437 set.seed(1234)
|
|
|
438 dataX <- rawData[,1:length(rawData)-1]
|
|
|
439 dataY <- rawData[,length(rawData)]
|
|
|
440
|
|
|
441 Smpling <- "$SAAMPLING"
|
|
|
442
|
|
|
443 if(Smpling=="downsampling")
|
|
|
444 {
|
|
|
445 dwnsmpl <- downSample(dataX,dataY)
|
|
|
446 rawData <- dwnsmpl[,1:length(dwnsmpl)-1]
|
|
|
447 rawData$outcome <- dwnsmpl[,length(dwnsmpl)]
|
|
|
448 remove("dwnsmpl")
|
|
|
449 remove("dataX")
|
|
|
450 remove("dataY")
|
|
|
451 }else if(Smpling=="upsampling"){
|
|
|
452 upsmpl <- upSample(dataX,dataY)
|
|
|
453 rawData <- upsmpl[,1:length(upsmpl)-1]
|
|
|
454 rawData$outcome <- upsmpl[,length(upsmpl)]
|
|
|
455 remove("upsmpl")
|
|
|
456 remove("dataX")
|
|
|
457 remove("dataY")
|
|
|
458 }else{remove("dataX")
|
|
|
459 remove("dataY")
|
|
|
460 }
|
|
|
461
|
|
|
462
|
|
|
463
|
|
|
464 numSamples <- nrow(rawData)
|
|
|
465
|
|
|
466 predictorNames <- names(rawData)[names(rawData) != "outcome"]
|
|
|
467 numPredictors <- length(predictorNames)
|
|
|
468
|
|
|
469
|
|
|
470 classDist1 <- table(rawData$outcome)
|
|
|
471 classDistString1 <- paste("``",
|
|
|
472 names(classDist1),
|
|
|
473 "'' ($n$=",
|
|
|
474 classDist1,
|
|
|
475 ")",
|
|
|
476 sep = "")
|
|
|
477 classDistString1 <- listString(classDistString1)
|
|
|
478
|
|
|
479 pctTrain <- $PERCENT
|
|
|
480
|
|
|
481 if(pctTrain < 1)
|
|
|
482 {
|
|
|
483 ## OPTION: seed number can be changed
|
|
|
484 set.seed(1)
|
|
|
485 inTrain <- createDataPartition(rawData$$outcome,
|
|
|
486 p = pctTrain,
|
|
|
487 list = FALSE)
|
|
|
488 trainX <- rawData[ inTrain, predictorNames]
|
|
|
489 testX <- rawData[-inTrain, predictorNames]
|
|
|
490 trainY <- rawData[ inTrain, "outcome"]
|
|
|
491 testY <- rawData[-inTrain, "outcome"]
|
|
|
492 splitText <- paste("The original data were split into ",
|
|
|
493 "a training set ($$n$$=",
|
|
|
494 nrow(trainX),
|
|
|
495 ") and a test set ($$n$$=",
|
|
|
496 nrow(testX),
|
|
|
497 ") in a manner that preserved the ",
|
|
|
498 "distribution of the classes.",
|
|
|
499 sep = "")
|
|
|
500 isZeroVar <- apply(trainX, 2,
|
|
|
501 function(x) length(unique(x)) < 2)
|
|
|
502 if(any(isZeroVar))
|
|
|
503 {
|
|
|
504 trainX <- trainX[, !isZeroVar, drop = FALSE]
|
|
|
505 testX <- testX[, !isZeroVar, drop = FALSE]
|
|
|
506 }
|
|
|
507
|
|
|
508 } else {
|
|
|
509 trainX <- rawData[, predictorNames]
|
|
|
510 testX <- NULL
|
|
|
511 trainY <- rawData[, "outcome"]
|
|
|
512 testY <- NULL
|
|
|
513 splitText <- "The entire data set was used as the training set."
|
|
|
514 }
|
|
|
515 trainDist <- table(trainY)
|
|
|
516 nir <- max(trainDist)/length(trainY)*100
|
|
|
517 niClass <- names(trainDist)[which.max(trainDist)]
|
|
|
518 nirText <- paste("The non--information rate is the accuracy that can be ",
|
|
|
519 "achieved by predicting all samples using the most ",
|
|
|
520 "dominant class. For these data, the rate is ",
|
|
|
521 round(nir, 2), "\\\\% using the ``",
|
|
|
522 niClass,
|
|
|
523 "'' class.",
|
|
|
524 sep = "")
|
|
|
525
|
|
|
526 remove("rawData")
|
|
|
527
|
|
|
528 if((!is.null(testX)) && (!is.null(testY))){
|
|
|
529 save(trainX,trainY,testX,testY,file="datasets.RData")
|
|
|
530 } else {
|
|
|
531 save(trainX,trainY,file="datasets.RData")
|
|
|
532 }
|
|
|
533
|
|
|
534 @
|
|
|
535
|
|
|
536 \Sexpr{splitText}
|
|
|
537
|
|
|
538 \Sexpr{nirText}
|
|
|
539
|
|
|
540 The data set for model building consisted of \Sexpr{numSamples} samples and
|
|
|
541 \Sexpr{numPredictors} predictor variables. The breakdown of the
|
|
|
542 outcome data classes were: \Sexpr{classDistString1}.
|
|
|
543
|
|
|
544 <<nzv, eval= $NZVEVAL, results = $NZVRESULT, echo = $NZVECHO>>=
|
|
|
545 ## OPTION: other pre-processing steps can be used
|
|
|
546 ppSteps <- caret:::suggestions(modName)
|
|
|
547
|
|
|
548 set.seed(2)
|
|
|
549 if(ppSteps["nzv"])
|
|
|
550 {
|
|
|
551 nzv <- nearZeroVar(trainX)
|
|
|
552 if(length(nzv) > 0)
|
|
|
553 {
|
|
|
554 nzvVars <- names(trainX)[nzv]
|
|
|
555 trainX <- trainX[, -nzv]
|
|
|
556 nzvText <- paste("There were ",
|
|
|
557 length(nzv),
|
|
|
558 " predictors that were removed from train set due to",
|
|
|
559 " severely unbalanced distributions that",
|
|
|
560 " could negatively affect the model",
|
|
|
561 ifelse(length(nzv) > 10,
|
|
|
562 ".",
|
|
|
563 paste(": ",
|
|
|
564 listString(nzvVars),
|
|
|
565 ".",
|
|
|
566 sep = "")),
|
|
|
567 sep = "")
|
|
|
568 testX <- testX[, -nzv]
|
|
|
569 } else nzvText <- ""
|
|
|
570 } else nzvText <- ""
|
|
|
571 @
|
|
|
572
|
|
|
573 \Sexpr{nzvText}
|
|
|
574
|
|
|
575
|
|
|
576 <<corrFilter, eval = $CORRFILTEREVAL, results = $CORRFILTERRESULT, echo = $CORRFILTERECHO>>=
|
|
|
577 if(ppSteps["corr"])
|
|
|
578 {
|
|
|
579 ## OPTION:
|
|
|
580 corrThresh <- $THRESHHOLDCOR
|
|
|
581 highCorr <- findCorrelation(cor(trainX, use = "pairwise.complete.obs"),
|
|
|
582 corrThresh)
|
|
|
583 if(length(highCorr) > 0)
|
|
|
584 {
|
|
|
585 corrVars <- names(trainX)[highCorr]
|
|
|
586 trainX <- trainX[, -highCorr]
|
|
|
587 corrText <- paste("There were ",
|
|
|
588 length(highCorr),
|
|
|
589 " predictors that were removed due to",
|
|
|
590 " large between--predictor correlations that",
|
|
|
591 " could negatively affect the model fit",
|
|
|
592 ifelse(length(highCorr) > 10,
|
|
|
593 ".",
|
|
|
594 paste(": ",
|
|
|
595 listString(highCorr),
|
|
|
596 ".",
|
|
|
597 sep = "")),
|
|
|
598 " Removing these predictors forced",
|
|
|
599 " all pair--wise correlations to be",
|
|
|
600 " less than ",
|
|
|
601 corrThresh,
|
|
|
602 ".",
|
|
|
603 sep = "")
|
|
|
604 testX <- testX[, -highCorr]
|
|
|
605 } else corrText <- "No correlation among data on given threshold"
|
|
|
606 }else corrText <- ""
|
|
|
607 @
|
|
|
608
|
|
|
609 \Sexpr{corrText}
|
|
|
610
|
|
|
611 <<preProc, eval = $PREPROCEVAL, echo = $PREPROCECHO, results = $PREPROCRESULT>>=
|
|
|
612 ppMethods <- NULL
|
|
|
613 if(ppSteps["center"]) ppMethods <- c(ppMethods, "center")
|
|
|
614 if(ppSteps["scale"]) ppMethods <- c(ppMethods, "scale")
|
|
|
615 if(any(hasMissing) > 0) ppMethods <- c(ppMethods, "knnImpute")
|
|
|
616 ##OPTION other methods, such as spatial sign, can be added to this list
|
|
|
617
|
|
|
618 if(length(ppMethods) > 0)
|
|
|
619 {
|
|
|
620 ppInfo <- preProcess(trainX, method = ppMethods)
|
|
|
621 trainX <- predict(ppInfo, trainX)
|
|
|
622 if(pctTrain < 1) testX <- predict(ppInfo, testX)
|
|
|
623 ppText <- paste("The following pre--processing methods were",
|
|
|
624 " applied to the training",
|
|
|
625 ifelse(pctTrain < 1, " and test", ""),
|
|
|
626 " data: ",
|
|
|
627 listString(ppMethods),
|
|
|
628 ".",
|
|
|
629 sep = "")
|
|
|
630 ppText <- gsub("center", "mean centering", ppText)
|
|
|
631 ppText <- gsub("scale", "scaling to unit variance", ppText)
|
|
|
632 ppText <- gsub("knnImpute",
|
|
|
633 paste(ppInfo$$k, "--nearest neighbor imputation", sep = ""),
|
|
|
634 ppText)
|
|
|
635 ppText <- gsub("spatialSign", "the spatial sign transformation", ppText)
|
|
|
636 ppText <- gsub("pca", "principal component feature extraction", ppText)
|
|
|
637 ppText <- gsub("ica", "independent component feature extraction", ppText)
|
|
|
638 } else {
|
|
|
639 ppInfo <- NULL
|
|
|
640 ppText <- ""
|
|
|
641 }
|
|
|
642
|
|
|
643 predictorNames <- names(trainX)
|
|
|
644 if(nzvText != "" | corrText != "" | ppText != "")
|
|
|
645 {
|
|
|
646 varText <- paste("After pre--processing, ",
|
|
|
647 ncol(trainX),
|
|
|
648 "predictors remained for modeling.")
|
|
|
649 } else varText <- ""
|
|
|
650
|
|
|
651 @
|
|
|
652
|
|
|
653 \Sexpr{ppText}
|
|
|
654 \Sexpr{varText}
|
|
|
655
|
|
|
656 \clearpage
|
|
|
657
|
|
|
658 \section*{Model Building}
|
|
|
659
|
|
|
660 <<setupWorkers, eval = TRUE, echo = $SETUPWORKERSECHO, results = $SETUPWORKERSRESULT>>=
|
|
|
661 numWorkers <- $NUMWORKERS
|
|
|
662 ##OPTION: turn up numWorkers to use MPI
|
|
|
663 if(numWorkers > 1)
|
|
|
664 {
|
|
|
665 mpiCalcs <- function(X, FUN, ...)
|
|
|
666 {
|
|
|
667 theDots <- list(...)
|
|
|
668 parLapply(theDots$$cl, X, FUN)
|
|
|
669 }
|
|
|
670
|
|
|
671 library(snow)
|
|
|
672 cl <- makeCluster(numWorkers, "MPI")
|
|
|
673 }
|
|
|
674 @
|
|
|
675
|
|
|
676 <<setupResampling, echo = $SETUPRESAMPLINGECHO, results = $SETUPRESAMPLINGRESULT>>=
|
|
|
677 ##OPTION: the resampling options can be changed. See
|
|
|
678 ## ?trainControl for details
|
|
|
679
|
|
|
680 resampName <- "$RESAMPNAME"
|
|
|
681 resampNumber <- $RESAMPLENUMBER
|
|
|
682 numRepeat <- $NUMREPEAT
|
|
|
683 resampP <- $RESAMPLENUMBERPERCENT
|
|
|
684
|
|
|
685 modelInfo <- modelLookup(modName)
|
|
|
686
|
|
|
687 if(numClasses == 2)
|
|
|
688 {
|
|
|
689 foo <- if(any(modelInfo$$probModel)) twoClassSummary else twoClassNoProbs
|
|
|
690 } else foo <- defaultSummary
|
|
|
691
|
|
|
692 set.seed(3)
|
|
|
693 ctlObj <- trainControl(method = resampName,
|
|
|
694 number = resampNumber,
|
|
|
695 repeats = numRepeat,
|
|
|
696 p = resampP,
|
|
|
697 classProbs = any(modelInfo$$probModel),
|
|
|
698 summaryFunction = foo)
|
|
|
699
|
|
|
700
|
|
|
701 ##OPTION select other performance metrics as needed
|
|
|
702 optMetric <- if(numClasses == 2 & any(modelInfo$$probModel)) "ROC" else "Kappa"
|
|
|
703
|
|
|
704 if(numWorkers > 1)
|
|
|
705 {
|
|
|
706 ctlObj$$workers <- numWorkers
|
|
|
707 ctlObj$$computeFunction <- mpiCalcs
|
|
|
708 ctlObj$$computeArgs <- list(cl = cl)
|
|
|
709 }
|
|
|
710 @
|
|
|
711
|
|
|
712 <<setupGrid, results = $SETUPGRIDRESULT, echo = $SETUPGRIDECHO>>=
|
|
|
713 ##OPTION expand or contract these grids as needed (or
|
|
|
714 ## add more models
|
|
|
715
|
|
|
716 gridSize <- $SETUPGRIDSIZE
|
|
|
717
|
|
|
718 if(modName %in% c("svmPoly", "svmRadial", "svmLinear", "lvq", "ctree2", "ctree")) gridSize <- 5
|
|
|
719 if(modName %in% c("earth", "fda")) gridSize <- 7
|
|
|
720 if(modName %in% c("knn", "rocc", "glmboost", "rf", "nodeHarvest")) gridSize <- 10
|
|
|
721
|
|
|
722 if(modName %in% c("nb")) gridSize <- 2
|
|
|
723 if(modName %in% c("pam", "rpart")) gridSize <- 15
|
|
|
724 if(modName %in% c("pls")) gridSize <- min(20, ncol(trainX))
|
|
|
725
|
|
|
726 if(modName == "gbm")
|
|
|
727 {
|
|
|
728 tGrid <- expand.grid(.interaction.depth = -1 + (1:5)*2 ,
|
|
|
729 .n.trees = (1:10)*20,
|
|
|
730 .shrinkage = .1)
|
|
|
731 }
|
|
|
732
|
|
|
733 if(modName == "nnet")
|
|
|
734 {
|
|
|
735 tGrid <- expand.grid(.size = -1 + (1:5)*2 ,
|
|
|
736 .decay = c(0, .001, .01, .1))
|
|
|
737 }
|
|
|
738
|
|
|
739 if(modName == "ada")
|
|
|
740 {
|
|
|
741 tGrid <- expand.grid(.maxdepth = 1, .iter = c(100,200,300,400), .nu = 1 )
|
|
|
742
|
|
|
743 }
|
|
|
744
|
|
|
745
|
|
|
746 @
|
|
|
747
|
|
|
748 <<fitModel, results = $FITMODELRESULT, echo = $FITMODELECHO, eval = $FITMODELEVAL>>=
|
|
|
749 ##OPTION alter as needed
|
|
|
750
|
|
|
751 set.seed(4)
|
|
|
752 modelFit <- switch(modName,
|
|
|
753 gbm =
|
|
|
754 {
|
|
|
755 mix <- sample(seq(along = trainY))
|
|
|
756 train(
|
|
|
757 trainX[mix,], trainY[mix], modName,
|
|
|
758 verbose = FALSE,
|
|
|
759 bag.fraction = .9,
|
|
|
760 metric = optMetric,
|
|
|
761 trControl = ctlObj,
|
|
|
762 tuneGrid = tGrid)
|
|
|
763 },
|
|
|
764
|
|
|
765 multinom =
|
|
|
766 {
|
|
|
767 train(
|
|
|
768 trainX, trainY, modName,
|
|
|
769 trace = FALSE,
|
|
|
770 metric = optMetric,
|
|
|
771 maxiter = 1000,
|
|
|
772 MaxNWts = 5000,
|
|
|
773 trControl = ctlObj,
|
|
|
774 tuneLength = gridSize)
|
|
|
775 },
|
|
|
776
|
|
|
777 nnet =
|
|
|
778 {
|
|
|
779 train(
|
|
|
780 trainX, trainY, modName,
|
|
|
781 metric = optMetric,
|
|
|
782 linout = FALSE,
|
|
|
783 trace = FALSE,
|
|
|
784 maxiter = 1000,
|
|
|
785 MaxNWts = 5000,
|
|
|
786 trControl = ctlObj,
|
|
|
787 tuneGrid = tGrid)
|
|
|
788
|
|
|
789 },
|
|
|
790
|
|
|
791 svmRadial =, svmPoly =, svmLinear =
|
|
|
792 {
|
|
|
793 train(
|
|
|
794 trainX, trainY, modName,
|
|
|
795 metric = optMetric,
|
|
|
796 scaled = TRUE,
|
|
|
797 trControl = ctlObj,
|
|
|
798 tuneLength = gridSize)
|
|
|
799 },
|
|
|
800 {
|
|
|
801 train(trainX, trainY, modName,
|
|
|
802 trControl = ctlObj,
|
|
|
803 metric = optMetric,
|
|
|
804 tuneLength = gridSize)
|
|
|
805 })
|
|
|
806
|
|
|
807 @
|
|
|
808
|
|
|
809 <<modelDescr, echo = $MODELDESCRECHO, results = $MODELDESCRRESULT>>=
|
|
|
810 summaryText <- ""
|
|
|
811
|
|
|
812 resampleName <- switch(tolower(modelFit$$control$$method),
|
|
|
813 boot = paste("the bootstrap (", length(modelFit$$control$$index), " reps)", sep = ""),
|
|
|
814 boot632 = paste("the bootstrap 632 rule (", length(modelFit$$control$$index), " reps)", sep = ""),
|
|
|
815 cv = paste("cross-validation (", modelFit$$control$$number, " fold)", sep = ""),
|
|
|
816 repeatedcv = paste("cross-validation (", modelFit$$control$$number, " fold, repeated ",
|
|
|
817 modelFit$$control$$repeats, " times)", sep = ""),
|
|
|
818 lgocv = paste("repeated train/test splits (", length(modelFit$$control$$index), " reps, ",
|
|
|
819 round(modelFit$$control$$p, 2), "$$\\%$$)", sep = ""))
|
|
|
820
|
|
|
821 tuneVars <- latexTranslate(tolower(modelInfo$$label))
|
|
|
822 tuneVars <- gsub("\\#", "the number of ", tuneVars, fixed = TRUE)
|
|
|
823 if(ncol(modelFit$$bestTune) == 1 && colnames(modelFit$$bestTune) == ".parameter")
|
|
|
824 {
|
|
|
825 summaryText <- paste(summaryText,
|
|
|
826 "\n\n",
|
|
|
827 "There are no tuning parameters associated with this model.",
|
|
|
828 "To characterize the model performance on the training set,",
|
|
|
829 resampleName,
|
|
|
830 "was used.",
|
|
|
831 "Table \\\\ref{T:resamps} and Figure \\\\ref{F:profile}",
|
|
|
832 "show summaries of the resampling results. ")
|
|
|
833
|
|
|
834 } else {
|
|
|
835 summaryText <- paste("There",
|
|
|
836 ifelse(nrow(modelInfo) > 1, "are", "is"),
|
|
|
837 nrow(modelInfo),
|
|
|
838 ifelse(nrow(modelInfo) > 1, "tuning parameters", "tuning parameter"),
|
|
|
839 "associated with this model:",
|
|
|
840 listString(tuneVars, period = TRUE))
|
|
|
841
|
|
|
842
|
|
|
843
|
|
|
844 paramNames <- gsub(".", "", names(modelFit$$bestTune), fixed = TRUE)
|
|
|
845 for(i in seq(along = paramNames))
|
|
|
846 {
|
|
|
847 check <- modelInfo$$parameter %in% paramNames[i]
|
|
|
848 if(any(check))
|
|
|
849 {
|
|
|
850 paramNames[i] <- modelInfo$$label[which(check)]
|
|
|
851 }
|
|
|
852 }
|
|
|
853
|
|
|
854 paramNames <- gsub("#", "the number of ", paramNames, fixed = TRUE)
|
|
|
855 ## Check to see if there was only one combination fit
|
|
|
856 summaryText <- paste(summaryText,
|
|
|
857 "To choose",
|
|
|
858 ifelse(nrow(modelInfo) > 1,
|
|
|
859 "appropriate values of the tuning parameters,",
|
|
|
860 "an appropriate value of the tuning parameter,"),
|
|
|
861 resampleName,
|
|
|
862 "was used to generated a profile of performance across the",
|
|
|
863 nrow(modelFit$$results),
|
|
|
864 ifelse(nrow(modelInfo) > 1,
|
|
|
865 "combinations of the tuning parameters.",
|
|
|
866 "candidate values."),
|
|
|
867
|
|
|
868 "Table \\\\ref{T:resamps} and Figure \\\\ref{F:profile} show",
|
|
|
869 "summaries of the resampling profile. ", "The final model fitted to the entire training set was:",
|
|
|
870 listString(paste(latexTranslate(tolower(paramNames)), "=", modelFit$$bestTune[1,]), period = TRUE))
|
|
|
871
|
|
|
872 }
|
|
|
873 @
|
|
|
874
|
|
|
875 \Sexpr{summaryText}
|
|
|
876
|
|
|
877 <<resampTable, echo = $RESAMPTABLEECHO, results = $RESAMPTABLERESULT>>=
|
|
|
878 tableData <- modelFit$$results
|
|
|
879
|
|
|
880 if(all(modelInfo$$parameter == "parameter") && resampName == "boot632")
|
|
|
881 {
|
|
|
882 tableData <- tableData[,-1, drop = FALSE]
|
|
|
883 colNums <- c( length(modelFit$$perfNames), length(modelFit$$perfNames), length(modelFit$$perfNames))
|
|
|
884 colLabels <- c("Mean", "Standard Deviation","Apparant")
|
|
|
885 constString <- ""
|
|
|
886 isConst <- NULL
|
|
|
887 } else if (all(modelInfo$$parameter == "parameter") && (resampName == "boot" | resampName == "cv" | resampName == "repeatedcv" )){
|
|
|
888 tableData <- tableData[,-1, drop = FALSE]
|
|
|
889 colNums <- c(length(modelFit$$perfNames), length(modelFit$$perfNames))
|
|
|
890 colLabels <- c("Mean", "Standard Deviation")
|
|
|
891 constString <- ""
|
|
|
892 isConst <- NULL
|
|
|
893 } else if (all(modelInfo$$parameter == "parameter") && resampName == "LOOCV" ){
|
|
|
894 tableData <- tableData[,-1, drop = FALSE]
|
|
|
895 colNums <- length(modelFit$$perfNames)
|
|
|
896 colLabels <- c("Measures")
|
|
|
897 constString <- ""
|
|
|
898 isConst <- NULL
|
|
|
899 } else {
|
|
|
900 if (all(modelInfo$$parameter != "parameter") && resampName == "boot632" ){
|
|
|
901 isConst <- apply(tableData[, modelInfo$$parameter, drop = FALSE],
|
|
|
902 2,
|
|
|
903 function(x) length(unique(x)) == 1)
|
|
|
904
|
|
|
905 numParamInTable <- sum(!isConst)
|
|
|
906
|
|
|
907 if(any(isConst))
|
|
|
908 {
|
|
|
909 constParam <- modelInfo$$parameter[isConst]
|
|
|
910 constValues <- format(tableData[, constParam, drop = FALSE], digits = 4)[1,,drop = FALSE]
|
|
|
911 tableData <- tableData[, !(names(tableData) %in% constParam), drop = FALSE]
|
|
|
912 constString <- paste("The tuning",
|
|
|
913 ifelse(sum(isConst) > 1,
|
|
|
914 "parmeters",
|
|
|
915 "parameter"),
|
|
|
916 listString(paste("``", names(constValues), "''", sep = "")),
|
|
|
917 ifelse(sum(isConst) > 1,
|
|
|
918 "were",
|
|
|
919 "was"),
|
|
|
920 "held constant at",
|
|
|
921 ifelse(sum(isConst) > 1,
|
|
|
922 "a value of",
|
|
|
923 "values of"),
|
|
|
924 listString(constValues[1,]))
|
|
|
925
|
|
|
926 } else constString <- ""
|
|
|
927
|
|
|
928 cn <- colnames(tableData)
|
|
|
929 for(i in seq(along = cn))
|
|
|
930 {
|
|
|
931 check <- modelInfo$$parameter %in% cn[i]
|
|
|
932 if(any(check))
|
|
|
933 {
|
|
|
934 cn[i] <- modelInfo$$label[which(check)]
|
|
|
935 }
|
|
|
936 }
|
|
|
937 colnames(tableData) <- cn
|
|
|
938
|
|
|
939 colNums <- c(numParamInTable,
|
|
|
940 length(modelFit$$perfNames),
|
|
|
941 length(modelFit$$perfNames),
|
|
|
942 length(modelFit$$perfNames))
|
|
|
943 colLabels <- c("", "Mean", "Standard Deviation", "Apparant")
|
|
|
944
|
|
|
945 }else if (all(modelInfo$$parameter != "parameter") && (resampName == "boot" | resampName == "repeatedcv" | resampName == "cv") ){
|
|
|
946 isConst <- apply(tableData[, modelInfo$$parameter, drop = FALSE],
|
|
|
947 2,
|
|
|
948 function(x) length(unique(x)) == 1)
|
|
|
949
|
|
|
950 numParamInTable <- sum(!isConst)
|
|
|
951
|
|
|
952 if(any(isConst))
|
|
|
953 {
|
|
|
954 constParam <- modelInfo$$parameter[isConst]
|
|
|
955 constValues <- format(tableData[, constParam, drop = FALSE], digits = 4)[1,,drop = FALSE]
|
|
|
956 tableData <- tableData[, !(names(tableData) %in% constParam), drop = FALSE]
|
|
|
957 constString <- paste("The tuning",
|
|
|
958 ifelse(sum(isConst) > 1,
|
|
|
959 "parmeters",
|
|
|
960 "parameter"),
|
|
|
961 listString(paste("``", names(constValues), "''", sep = "")),
|
|
|
962 ifelse(sum(isConst) > 1,
|
|
|
963 "were",
|
|
|
964 "was"),
|
|
|
965 "held constant at",
|
|
|
966 ifelse(sum(isConst) > 1,
|
|
|
967 "a value of",
|
|
|
968 "values of"),
|
|
|
969 listString(constValues[1,]))
|
|
|
970
|
|
|
971 } else constString <- ""
|
|
|
972
|
|
|
973 cn <- colnames(tableData)
|
|
|
974 for(i in seq(along = cn))
|
|
|
975 {
|
|
|
976 check <- modelInfo$$parameter %in% cn[i]
|
|
|
977 if(any(check))
|
|
|
978 {
|
|
|
979 cn[i] <- modelInfo$$label[which(check)]
|
|
|
980 }
|
|
|
981 }
|
|
|
982 colnames(tableData) <- cn
|
|
|
983
|
|
|
984 colNums <- c(numParamInTable,
|
|
|
985 length(modelFit$$perfNames),
|
|
|
986 length(modelFit$$perfNames))
|
|
|
987 colLabels <- c("", "Mean", "Standard Deviation")
|
|
|
988
|
|
|
989 }
|
|
|
990 else if (all(modelInfo$$parameter != "parameter") && resampName == "LOOCV"){
|
|
|
991 isConst <- apply(tableData[, modelInfo$$parameter, drop = FALSE],
|
|
|
992 2,
|
|
|
993 function(x) length(unique(x)) == 1)
|
|
|
994
|
|
|
995 numParamInTable <- sum(!isConst)
|
|
|
996
|
|
|
997 if(any(isConst))
|
|
|
998 {
|
|
|
999 constParam <- modelInfo$$parameter[isConst]
|
|
|
1000 constValues <- format(tableData[, constParam, drop = FALSE], digits = 4)[1,,drop = FALSE]
|
|
|
1001 tableData <- tableData[, !(names(tableData) %in% constParam), drop = FALSE]
|
|
|
1002 constString <- paste("The tuning",
|
|
|
1003 ifelse(sum(isConst) > 1,
|
|
|
1004 "parmeters",
|
|
|
1005 "parameter"),
|
|
|
1006 listString(paste("``", names(constValues), "''", sep = "")),
|
|
|
1007 ifelse(sum(isConst) > 1,
|
|
|
1008 "were",
|
|
|
1009 "was"),
|
|
|
1010 "held constant at",
|
|
|
1011 ifelse(sum(isConst) > 1,
|
|
|
1012 "a value of",
|
|
|
1013 "values of"),
|
|
|
1014 listString(constValues[1,]))
|
|
|
1015
|
|
|
1016 } else constString <- ""
|
|
|
1017
|
|
|
1018 cn <- colnames(tableData)
|
|
|
1019 for(i in seq(along = cn))
|
|
|
1020 {
|
|
|
1021 check <- modelInfo$$parameter %in% cn[i]
|
|
|
1022 if(any(check))
|
|
|
1023 {
|
|
|
1024 cn[i] <- modelInfo$$label[which(check)]
|
|
|
1025 }
|
|
|
1026 }
|
|
|
1027 colnames(tableData) <- cn
|
|
|
1028
|
|
|
1029 colNums <- c(numParamInTable,
|
|
|
1030 length(modelFit$$perfNames))
|
|
|
1031 colLabels <- c("", "Measures")
|
|
|
1032
|
|
|
1033 }
|
|
|
1034
|
|
|
1035 }
|
|
|
1036
|
|
|
1037
|
|
|
1038
|
|
|
1039 colnames(tableData) <- gsub("SD$$", "", colnames(tableData))
|
|
|
1040 colnames(tableData) <- gsub("Apparent$$", "", colnames(tableData))
|
|
|
1041 colnames(tableData) <- latexTranslate(colnames(tableData))
|
|
|
1042 rownames(tableData) <- latexTranslate(rownames(tableData))
|
|
|
1043
|
|
|
1044 latex(tableData,
|
|
|
1045 rowname = NULL,
|
|
|
1046 file = "",
|
|
|
1047 cgroup = colLabels,
|
|
|
1048 n.cgroup = colNums,
|
|
|
1049 where = "h!",
|
|
|
1050 digits = 4,
|
|
|
1051 longtable = nrow(tableData) > 30,
|
|
|
1052 caption = paste(resampleName, "results from the model fit.", constString),
|
|
|
1053 label = "T:resamps")
|
|
|
1054 @
|
|
|
1055
|
|
|
1056 \setkeys{Gin}{ width = 0.9\textwidth}
|
|
|
1057 \begin{figure}[b]
|
|
|
1058 \begin{center}
|
|
|
1059
|
|
|
1060 <<profilePlot, echo = $PROFILEPLOTECHO, fig = $PROFILEPLOTFIG, width = 8, height = 6>>=
|
|
|
1061 trellis.par.set(caretTheme(), warn = TRUE)
|
|
|
1062 if(all(modelInfo$$parameter == "parameter") | all(isConst) | modName == "nb")
|
|
|
1063 {
|
|
|
1064 resultsPlot <- resampleHist(modelFit)
|
|
|
1065 plotCaption <- paste("Distributions of model performance from the ",
|
|
|
1066 "training set estimated using ",
|
|
|
1067 resampleName)
|
|
|
1068 } else {
|
|
|
1069 if(modName %in% c("svmPoly", "svmRadial", "svmLinear"))
|
|
|
1070 {
|
|
|
1071 resultsPlot <- plot(modelFit,
|
|
|
1072 metric = optMetric,
|
|
|
1073 xTrans = function(x) log10(x))
|
|
|
1074 resultsPlot <- update(resultsPlot,
|
|
|
1075 type = c("g", "p", "l"),
|
|
|
1076 ylab = paste(optMetric, " (", resampleName, ")", sep = ""))
|
|
|
1077
|
|
|
1078 } else {
|
|
|
1079 resultsPlot <- plot(modelFit,
|
|
|
1080 metric = optMetric)
|
|
|
1081 resultsPlot <- update(resultsPlot,
|
|
|
1082 type = c("g", "p", "l"),
|
|
|
1083 ylab = paste(optMetric, " (", resampleName, ")", sep = ""))
|
|
|
1084 }
|
|
|
1085 plotCaption <- paste("A plot of the estimates of the",
|
|
|
1086 optMetric,
|
|
|
1087 "values calculated using",
|
|
|
1088 resampleName)
|
|
|
1089 }
|
|
|
1090 print(resultsPlot)
|
|
|
1091 @
|
|
|
1092 \caption[Performance Plot]{\Sexpr{plotCaption}.}
|
|
|
1093 \label{F:profile}
|
|
|
1094 \end{center}
|
|
|
1095 \end{figure}
|
|
|
1096
|
|
|
1097
|
|
|
1098 <<stopWorkers, echo = $STOPWORKERSECHO, results = $STOPWORKERSRESULT>>=
|
|
|
1099 if(numWorkers > 1) stopCluster(cl)
|
|
|
1100 @
|
|
|
1101
|
|
|
1102 <<testPred, results = $TESTPREDRESULT, echo = $TESTPREDECHO>>=
|
|
|
1103 if(pctTrain < 1)
|
|
|
1104 {
|
|
|
1105 cat("\\clearpage\n\\section*{Test Set Results}\n\n")
|
|
|
1106 classPred <- predict(modelFit, testX)
|
|
|
1107 cm <- confusionMatrix(classPred, testY)
|
|
|
1108 values <- cm$$overall[c("Accuracy", "Kappa", "AccuracyPValue", "McnemarPValue")]
|
|
|
1109
|
|
|
1110 values <- values[!is.na(values) & !is.nan(values)]
|
|
|
1111 values <- c(format(values[1:2], digits = 3),
|
|
|
1112 format.pval(values[-(1:2)], digits = 5))
|
|
|
1113 nms <- c("the overall accuracy", "the Kappa statistic",
|
|
|
1114 "the $$p$$--value that accuracy is greater than the no--information rate",
|
|
|
1115 "the $$p$$--value of concordance from McNemar's test")
|
|
|
1116 nms <- nms[seq(along = values)]
|
|
|
1117 names(values) <- nms
|
|
|
1118
|
|
|
1119 if(any(modelInfo$$probModel))
|
|
|
1120 {
|
|
|
1121 classProbs <- extractProb(list(fit = modelFit),
|
|
|
1122 testX = testX,
|
|
|
1123 testY = testY)
|
|
|
1124 classProbs <- subset(classProbs, dataType == "Test")
|
|
|
1125 if(numClasses == 2)
|
|
|
1126 {
|
|
|
1127 tmp <- twoClassSummary(classProbs, lev = levels(classProbs$$obs))
|
|
|
1128 tmp <- c(format(tmp, digits = 3))
|
|
|
1129 names(tmp) <- c("the sensitivity", "the specificity",
|
|
|
1130 "the area under the ROC curve")
|
|
|
1131 values <- c(values, tmp)
|
|
|
1132
|
|
|
1133 }
|
|
|
1134 probPlot <- plotClassProbs(classProbs)
|
|
|
1135 }
|
|
|
1136 testString <- paste("Based on the test set of",
|
|
|
1137 nrow(testX),
|
|
|
1138 "samples,",
|
|
|
1139 listString(paste(names(values), "was", values), period = TRUE),
|
|
|
1140 "The confusion matrix for the test set is shown in Table",
|
|
|
1141 "\\\\ref{T:cm}.")
|
|
|
1142 testString <- paste(testString,
|
|
|
1143 " Using ", resampleName,
|
|
|
1144 ", the training set estimates were ",
|
|
|
1145 resampleStats(modelFit),
|
|
|
1146 ".",
|
|
|
1147 sep = "")
|
|
|
1148
|
|
|
1149 if(any(modelInfo$$probModel)) testString <- paste(testString,
|
|
|
1150 "Histograms of the class probabilities",
|
|
|
1151 "for the test set samples are shown in",
|
|
|
1152 "Figure \\\\ref{F:probs}",
|
|
|
1153 ifelse(numClasses == 2,
|
|
|
1154 " and the test set ROC curve is in Figure \\\\ref{F:roc}.",
|
|
|
1155 "."))
|
|
|
1156
|
|
|
1157
|
|
|
1158
|
|
|
1159 latex(cm$$table,
|
|
|
1160 title = "",
|
|
|
1161 file = "",
|
|
|
1162 where = "h",
|
|
|
1163 cgroup = "Observed Values",
|
|
|
1164 n.cgroup = numClasses,
|
|
|
1165 caption = "The confusion matrix for the test set",
|
|
|
1166 label = "T:cm")
|
|
|
1167
|
|
|
1168 } else testString <- ""
|
|
|
1169 @
|
|
|
1170 \Sexpr{testString}
|
|
|
1171
|
|
|
1172
|
|
|
1173 <<classProbsTex, results = $CLASSPROBSTEXRESULT, echo = $CLASSPROBSTEXECHO>>=
|
|
|
1174 if(any(modelInfo$probModel) && pctTrain < 1 ) {
|
|
|
1175 cat(
|
|
|
1176 paste("\\begin{figure}[p]\n",
|
|
|
1177 "\\begin{center}\n",
|
|
|
1178 "\\includegraphics{classProbs}",
|
|
|
1179 "\\caption[PCA Plot]{Class probabilities",
|
|
|
1180 "for the test set. Each panel contains ",
|
|
|
1181 "separate classes}\n",
|
|
|
1182 "\\label{F:probs}\n",
|
|
|
1183 "\\end{center}\n",
|
|
|
1184 "\\end{figure}"))
|
|
|
1185 }
|
|
|
1186 if(any(modelInfo$$probModel) & numClasses == 2 & pctTrain < 1 )
|
|
|
1187 {
|
|
|
1188 cat(
|
|
|
1189 paste("\\begin{figure}[p]\n",
|
|
|
1190 "\\begin{center}\n",
|
|
|
1191 "\\includegraphics[clip, width = .8\\textwidth]{roc}",
|
|
|
1192 "\\caption[ROC Plot]{ROC Curve",
|
|
|
1193 "for the test set.}\n",
|
|
|
1194 "\\label{F:roc}\n",
|
|
|
1195 "\\end{center}\n",
|
|
|
1196 "\\end{figure}"))
|
|
|
1197 } else {
|
|
|
1198 cat (paste(""))
|
|
|
1199 }
|
|
|
1200
|
|
|
1201 @
|
|
|
1202 <<classProbsTex, results = $CLASSPROBSTEXRESULT1, echo = $CLASSPROBSTEXECHO1 >>=
|
|
|
1203 if(any(modelInfo$probModel) && pctTrain < 1) {
|
|
|
1204 pdf("classProbs.pdf", height = 7, width = 7)
|
|
|
1205 trellis.par.set(caretTheme(), warn = FALSE)
|
|
|
1206 print(probPlot)
|
|
|
1207 dev.off()
|
|
|
1208 }
|
|
|
1209 if(any(modelInfo$probModel) & numClasses == 2 & pctTrain < 1) {
|
|
|
1210 resPonse<-testY
|
|
|
1211 preDictor<-classProbs[, levels(trainY)[1]]
|
|
|
1212 pdf("roc.pdf", height = 8, width = 8)
|
|
|
1213 # from pROC example at http://web.expasy.org/pROC/screenshots.htm
|
|
|
1214 plot.roc(resPonse, preDictor, # data
|
|
|
1215 percent=TRUE, # show all values in percent
|
|
|
1216 partial.auc=c(100, 90), partial.auc.correct=TRUE, # define a partial AUC (pAUC)
|
|
|
1217 print.auc=TRUE, #display pAUC value on the plot with following options:
|
|
|
1218 print.auc.pattern="Corrected pAUC (100-90%% SP):\n%.1f%%", print.auc.col="#1c61b6",
|
|
|
1219 auc.polygon=TRUE, auc.polygon.col="#1c61b6", # show pAUC as a polygon
|
|
|
1220 max.auc.polygon=TRUE, max.auc.polygon.col="#1c61b622", # also show the 100% polygon
|
|
|
1221 main="Partial AUC (pAUC)")
|
|
|
1222 plot.roc(resPonse, preDictor,
|
|
|
1223 percent=TRUE, add=TRUE, type="n", # add to plot, but don't re-add the ROC itself (useless)
|
|
|
1224 partial.auc=c(100, 90), partial.auc.correct=TRUE,
|
|
|
1225 partial.auc.focus="se", # focus pAUC on the sensitivity
|
|
|
1226 print.auc=TRUE, print.auc.pattern="Corrected pAUC (100-90%% SE):\n%.1f%%", print.auc.col="#008600",
|
|
|
1227 print.auc.y=40, # do not print auc over the previous one
|
|
|
1228 auc.polygon=TRUE, auc.polygon.col="#008600",
|
|
|
1229 max.auc.polygon=TRUE, max.auc.polygon.col="#00860022")
|
|
|
1230 dev.off()
|
|
|
1231 } else {
|
|
|
1232 cat("")
|
|
|
1233 }
|
|
|
1234
|
|
|
1235 @
|
|
|
1236
|
|
|
1237 \section*{Versions}
|
|
|
1238
|
|
|
1239 <<versions, echo = FALSE, results = tex>>=
|
|
|
1240 toLatex(sessionInfo())
|
|
|
1241
|
|
|
1242 @
|
|
|
1243
|
|
|
1244 <<save-data, echo = $SAVEDATAECHO, results = $SAVEDATARESULT>>=
|
|
|
1245 ## change this to the name of modName....
|
|
|
1246 Fit<-modelFit
|
|
|
1247 save(Fit,file="$METHOD-Fit.RData")
|
|
|
1248 @
|
|
|
1249 The model was built using $METHOD and is saved as $METHOD Model for reuse. This contains the variable Fit.
|
|
|
1250
|
|
|
1251 \end{document}'''
|
|
|
1252
|
|
|
1253 return template4Rnw
|
