### R code from vignette source 'caretTrain.Rnw'

###################################################
### code chunk number 1: loadLibs
###################################################
library(caret)
library(kernlab)
library(gbm)
library(ipred)
library(grid)
library(Hmisc)
library(randomForest)
data(BloodBrain)
data(mdrr)

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


###################################################
### code chunk number 2: preProc
###################################################
print(ncol(mdrrDescr))
nzv <- nearZeroVar(mdrrDescr)
filteredDescr <- mdrrDescr[, -nzv]

print(ncol(filteredDescr))

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

print(ncol(filteredDescr))

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

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

print(length(trainMDRR))
print(length(testMDRR))


###################################################
### code chunk number 3: setup
###################################################
fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 10,
                           ## repeated three times
                           repeats = 3,
                           ## Save all the resampling results
                           returnResamp = "all")


###################################################
### code chunk number 4: mdrrModel1
###################################################
gbmFit1 <- train(trainDescr, trainMDRR, 
                 method = "gbm", 
                 trControl = fitControl,
                 ## This last option is actually one
                 ## for gbm() that passes through
                 verbose = FALSE)
gbmFit1


###################################################
### code chunk number 5: getSeqMods
###################################################
seqModList <- paste(paste("\\\\pkg{", unique(subset(modelLookup(), seq)$model), "}", sep = ""), collapse = ", ")


###################################################
### code chunk number 6: getModelData
###################################################

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

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

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

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

  }

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

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




###################################################
### code chunk number 7: dualTable
###################################################
dualCats <- unique(as.character(dual$Class))
dualReps <- table(as.character(dual$Class))
dualReps <- dualReps[dualCats]

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

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


###################################################
### code chunk number 8: regTable
###################################################
regCats <- unique(as.character(reg$Class))
regReps <- table(as.character(reg$Class))
regReps <- regReps[regCats]

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

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


###################################################
### code chunk number 9: clasTable
###################################################
clasCats <- unique(as.character(clas$Class))
clasReps <- table(as.character(clas$Class))
clasReps <- clasReps[clasCats]

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

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


###################################################
### code chunk number 10: mdrrGrid
###################################################
gbmGrid <-  expand.grid(.interaction.depth = c(1, 3), 
                        .n.trees = c(10, 50, 100, 150, 200, 250, 300), 
                        .shrinkage = 0.1)


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


###################################################
### code chunk number 12: summaryFunc
###################################################
twoClassSummary


###################################################
### code chunk number 13: setup
###################################################
fitControl <- trainControl(method = "repeatedcv",
                           number = 10,
                           repeats = 3,
                           returnResamp = "all",
                           ## Estimate class probabilities
                           classProbs = TRUE,
                           ## Evaluate performance using 
                           ## the following function
                           summaryFunction = twoClassSummary)


###################################################
### code chunk number 14: reTune
###################################################
set.seed(3)
gbmFit3 <- train(trainDescr, trainMDRR, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 tuneGrid = gbmGrid,
                 ## Specify which metric to optimize
                 metric = "ROC")
gbmFit3


###################################################
### code chunk number 15: bestGBM
###################################################
printSelected <- function(x)
  {
    tmp <- x$bestTune
    names(tmp) <- gsub(".", " ", names(tmp), fixed = TRUE)
    tmp <- paste(names(tmp), "=", tmp)
    paste(tmp, collapse = ", ")
  }
getTrainPerf <- function(x)
  {
    bst <- x$bestTune
    names(bst) <- substring(names(bst), 2)
    merge(bst, x$results)
  }


###################################################
### code chunk number 16: tolerance
###################################################
whichTwoPct <- tolerance(gbmFit3$results, "ROC", 2, TRUE)  
cat("best model within 2 pct of best:\n")
gbmFit3$results[whichTwoPct,]


###################################################
### code chunk number 17: makePlots
###################################################
pdf(paste(getwd(), "/gbm1.pdf", sep = ""), width = 5, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3))
dev.off()
pdf(paste(getwd(), "/gbm2.pdf", sep = ""), width = 5, height = 5)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3, metric = "Spec"))
dev.off()
pdf(paste(getwd(), "/gbm3.pdf", sep = ""), width = 8, height = 4)
   trellis.par.set(caretTheme(), warn = FALSE)
   print(plot(gbmFit3, meric = "ROC", plotType = "level"))
dev.off()



###################################################
### code chunk number 18: bhExample
###################################################
library(mlbench)
data(BostonHousing)
# we could use the formula interface too
bhDesignMatrix <-  model.matrix(medv ~. - 1, BostonHousing)


###################################################
### code chunk number 19: bhSplit
###################################################
set.seed(4)
inTrain <- createDataPartition(BostonHousing$medv, p = .8, list = FALSE, times = 1)
trainBH <- bhDesignMatrix[inTrain,]
testBH <- bhDesignMatrix[-inTrain,]

trainMedv <- BostonHousing$medv[inTrain]
testMedv <- BostonHousing$medv[-inTrain]


###################################################
### code chunk number 20: bhModels
###################################################
set.seed(5)
plsFit <- train(trainBH, trainMedv, 
                "pls", 
                preProcess = c("center", "scale"),
                tuneLength = 10, 
                trControl = trainControl(verboseIter = FALSE,
                                         returnResamp = "all"))
set.seed(5)
marsFit <- train(trainBH, trainMedv, 
                 "earth", 
                 tuneLength = 10, 
                 trControl = trainControl(verboseIter = FALSE,
                                          returnResamp = "all"))



###################################################
### code chunk number 21: marPrediction1
###################################################
marsPred1 <- predict(marsFit$finalModel, newdata = testBH)
head(marsPred1)


###################################################
### code chunk number 22: marsPrediction2
###################################################
marsPred2 <- predict(marsFit, newdata = testBH)
head(marsPred2)


###################################################
### code chunk number 23: bhPrediction1
###################################################
bhModels <- list(pls = plsFit, mars = marsFit)

bhPred1 <- predict(bhModels, newdata = testBH)
str(bhPred1)


###################################################
### code chunk number 24: bhPrediction1
###################################################
allPred <- extractPrediction(bhModels,
                             testX = testBH,
                             testY = testMedv)
testPred <- subset(allPred, dataType == "Test")
head(testPred)

ddply(testPred, .(model), defaultSummary)



###################################################
### code chunk number 25: mdrConfusion
###################################################
mdrrPredictions <- extractPrediction(list(gbmFit3), testX = testDescr, testY = testMDRR)
mdrrPredictions <- mdrrPredictions[mdrrPredictions$dataType == "Test",]
sensitivity(mdrrPredictions$pred, mdrrPredictions$obs)
confusionMatrix(mdrrPredictions$pred, mdrrPredictions$obs)


###################################################
### code chunk number 26: resampCM
###################################################
confusionMatrix(gbmFit3)


###################################################
### code chunk number 27: mdrrPlots
###################################################
mdrrProbs <- extractProb(list(gbmFit3), testX = testDescr, testY = testMDRR)
mdrrProbs <- mdrrProbs[mdrrProbs$dataType == "Test",]

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


###################################################
### code chunk number 28: bhPredPlot
###################################################
pdf("bhPredPlot.pdf", width = 8, height = 5)
trellis.par.set(caretTheme(), warn = FALSE)
print(plotObsVsPred(testPred))
dev.off()


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


###################################################
### code chunk number 30: bbbModelFits (eval = FALSE)
###################################################
## data(BloodBrain)
## set.seed(1)
## tmp <- createDataPartition(logBBB, p = 0.8, times = 100)
## bbbControl <- trainControl(method = "LGOCV", index = tmp, timingSamps = 50)
## 
## rpartFit <- train(bbbDescr, logBBB,
##                    "rpart", 
##                    tuneLength = 16,
##                    trControl = bbbControl)
## 
## ctreeFit <- train(bbbDescr, logBBB,
##                  "ctree2",
##                   tuneLength = 10,
##                   trControl = bbbControl)
## 
## earthFit <- train(bbbDescr, logBBB,
##                    "earth",
##                    tuneLength = 20,
##                    trControl = bbbControl)
## m5Fit <- train(bbbDescr, logBBB,
##                 "M5Rules",
##                 trControl = bbbControl)  


###################################################
### code chunk number 31: caretTrain.Rnw:1394-1395
###################################################
print(getwd())


###################################################
### code chunk number 32: resamps (eval = FALSE)
###################################################
## resamps <- resamples(list(CART = rpartFit,
##                           CondInfTree = ctreeFit,
##                           MARS = earthFit,
##                           M5 = m5Fit))


###################################################
### code chunk number 33: resampLoad
###################################################
load("resamps.RData")


###################################################
### code chunk number 34: resampFuncs
###################################################
resamps
summary(resamps)


###################################################
### code chunk number 35: resamplePlots
###################################################
bwplot(resamps, metric = "RMSE")

densityplot(resamps, metric = "RMSE")

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

splom(resamps, metric = "RMSE")


###################################################
### code chunk number 36: resamplePlots2
###################################################
pdf("resampleScatter.pdf", width = 6, height = 6)
trellis.par.set(caretTheme())
print(xyplot(resamps, models = c("CART", "MARS")))
dev.off()
pdf("resampleDens.pdf", width = 7, height = 4.5)
print(
      densityplot(resamps, 
                  scales = list(x = list(relation = "free")), 
                  adjust = 1.2, 
                  plot.points = FALSE, 
                  auto.key = list(columns = 2)))

dev.off() 


###################################################
### code chunk number 37: diffs
###################################################
difValues <- diff(resamps)

difValues

summary(difValues)


###################################################
### code chunk number 38: diffPlots
###################################################
dotplot(difValues)

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

bwplot(difValues,
       metric = "RMSE")

levelplot(difValues, what = "differences")


###################################################
### code chunk number 39: diffPlots2
###################################################
pdf("diffLevel.pdf", width = 7.5, height = 6)
trellis.par.set(caretTheme())
print(levelplot(difValues, what = "differences"))
dev.off()

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



###################################################
### code chunk number 40: getCp1
###################################################
## rpart requires a formula method
tmpData <- bbbDescr
tmpData$logBBB <- logBBB

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


###################################################
### code chunk number 41: cpGrid
###################################################
rpartGrid <- data.frame(.cp = c(cpValues10, cpValues30),
                        .minsplit = 
                        c(rep(10, length(cpValues10)), 
                          rep(30, length(cpValues30))))



###################################################
### code chunk number 42: modelFoo
###################################################
modelFunc <- function(data, parameter, levels, last, ...)
  {
    library(rpart)
    ctrl <- rpart.control(cp = parameter$.cp,
                          minsplit = parameter$.minsplit)
    
    list(fit = rpart(.outcome ~ ., data = data, control = ctrl))                             
  }


###################################################
### code chunk number 43: predFoo
###################################################
predFunc <- function(object, newdata)
  {
    library(rpart)
    predict(object$fit, newdata)                            
  }


###################################################
### code chunk number 44: sortFoo
###################################################
sortFunc <- function(x)  x[order(x$cp, x$minsplit),]                           


###################################################
### code chunk number 45: trControl (eval = FALSE)
###################################################
## ctrl <- trainControl(custom = list(
##                        parameters = rpartGrid,
##                        model = modelFunc,
##                        prediction = predFunc,
##                        probability = NULL,
##                        sort = sortFunc),
##                      method = "cv")
## 
## set.seed(581)
## customRpart <- train(bbbDescr, logBBB, "custom", trControl = ctrl)


###################################################
### code chunk number 46: plotResults (eval = FALSE)
###################################################
## rpartPlot <- plot(customRpart, xTrans = log10)
## rpartPlot <- update(rpartPlot, xlab = "Complexity Parameter")


###################################################
### code chunk number 47: session
###################################################
toLatex(sessionInfo())


