如何在CARET中自定义模型来执行PLS-[Classifier]两步分类模型？

这个问题是同一线程的延续here https://stats.stackexchange.com/questions/81727/what-is-the-best-strategy-to-train-and-validate-classification-using-pls-classi。下面是从本书中摘取的一个最小的工作示例：

Wehrens R. 化学计量学与 R 多变量数据分析 自然科学和生命科学。第一版。海德堡；纽约： 施普林格。 2011 年。（第 250 页）。

该示例取自本书及其包ChemometricsWithR。它强调了使用交叉验证技术建模时的一些陷阱。

The Aim:
使用同一组重复 CV 来执行已知策略的交叉验证方法PLS通常紧随其后的是LDA或者类似逻​​辑回归、SVM、C5.0、CART 的精神caret包裹。所以每次调用等待分类器之前都需要PLS来对PLS进行分类score空间而不是观察本身。插入符包中最近的方法是PCA作为使用任何分类器建模之前的预处理步骤。下面是一个 PLS-LDA 程序，只有一次交叉验证来测试分类器的性能，没有 10 倍 CV 或任何重复。下面的代码取自上述书籍，但进行了一些更正，否则会引发错误：

library(ChemometricsWithR)
data(prostate)
prostate.clmat <- classvec2classmat(prostate.type) # convert Y to a dummy var

odd <- seq(1, length(prostate.type), by = 2) # training
even <- seq(2, length(prostate.type), by = 2) # holdout test

prostate.pls <- plsr(prostate.clmat ~ prostate, ncomp = 16, validation = "CV", subset=odd)

Xtst <- scale(prostate[even,], center = colMeans(prostate[odd,]), scale = apply(prostate[odd,],2,sd))

tst.scores <- Xtst %*% prostate.pls$projection # scores for the waiting trained LDA to test

prostate.ldapls <- lda(scores(prostate.pls)[,1:16],prostate.type[odd]) # LDA for scores
table(predict(prostate.ldapls, new = tst.scores[,1:16])$class, prostate.type[even])

predictionTest <- predict(prostate.ldapls, new = tst.scores[,1:16])$class)

library(caret)    
confusionMatrix(data = predictionTest, reference= prostate.type[even]) # from caret

输出：

Confusion Matrix and Statistics

          Reference
Prediction bph control pca
   bph       4       1   9
   control   1      35   7
   pca      34       4  68

Overall Statistics

               Accuracy : 0.6564          
                 95% CI : (0.5781, 0.7289)
    No Information Rate : 0.5153          
    P-Value [Acc > NIR] : 0.0001874       

                  Kappa : 0.4072          
 Mcnemar's Test P-Value : 0.0015385       

Statistics by Class:

                     Class: bph Class: control Class: pca
Sensitivity             0.10256         0.8750     0.8095
Specificity             0.91935         0.9350     0.5190
Pos Pred Value          0.28571         0.8140     0.6415
Neg Pred Value          0.76510         0.9583     0.7193
Prevalence              0.23926         0.2454     0.5153
Detection Rate          0.02454         0.2147     0.4172
Detection Prevalence    0.08589         0.2638     0.6503
Balanced Accuracy       0.51096         0.9050     0.6643

然而，混淆矩阵与书中的不匹配，无论如何，书中的代码确实被破坏了，但这里的这个对我有用！

Notes:
虽然这只是一份简历，但其目的是首先就这个方法论达成一致，sd and mean将训练集的数据应用于测试集，PLUS根据特定数量的PC转化为PLS分数ncomp。我希望这种情况在插入符号中的每一轮简历中都发生。如果代码方法在这里是正确的，那么它可以作为修改插入符包的代码时的最小工作示例的良好开端。

附注：
缩放和居中可能会非常混乱，我认为 R 中的一些 PLS 函数在内部进行缩放，有或没有居中，我不确定，所以在插入符中构建自定义模型应该小心处理，以避免缺少或多重缩放或居中（我对这些事情保持警惕）。

多重居中/缩放的危险
下面的代码只是为了展示多重居中/缩放如何更改数据，此处仅显示居中，但缩放也存在相同的问题。

set.seed(1)
x <- rnorm(200, 2, 1)
xCentered1 <- scale(x, center=TRUE, scale=FALSE)
xCentered2 <- scale(xCentered1, center=TRUE, scale=FALSE)
xCentered3 <- scale(xCentered2, center=TRUE, scale=FALSE)
sapply (list(xNotCentered= x, xCentered1 = xCentered1, xCentered2 = xCentered2, xCentered3 = xCentered3), mean)

输出：

xNotCentered    xCentered1    xCentered2    xCentered3 
 2.035540e+00  1.897798e-16 -5.603699e-18 -5.332377e-18

如果我在本课程中遗漏了某些内容，请发表评论。谢谢。

如果您想将这些类型的模型与caret，您需要使用 CRAN 上的最新版本。创建最后一个更新是为了让人们可以使用非标型号 http://topepo.github.io/caret/custom_models.html他们认为合适的。

我下面的方法是联合拟合 PLS 和其他模型（我在下面的示例中使用随机森林）并同时调整它们。因此，对于每次折叠，都有一个 2D 网格ncomp and mtry用来。

“技巧”是将 PLS 载荷附加到随机森林对象，以便可以在预测期间使用它们。这是定义模型的代码（仅分类）：

 modelInfo <- list(label = "PLS-RF",
              library = c("pls", "randomForest"),
              type = "Classification",
              parameters = data.frame(parameter = c('ncomp', 'mtry'),
                                      class = c("numeric", 'numeric'),
                                      label = c('#Components', 
                                                '#Randomly Selected Predictors')),
              grid = function(x, y, len = NULL) {
                grid <- expand.grid(ncomp = seq(1, min(ncol(x) - 1, len), by = 1),
                            mtry = 1:len)
                grid <- subset(grid, mtry <= ncomp)
                },
              loop = NULL,
              fit = function(x, y, wts, param, lev, last, classProbs, ...) { 
                     ## First fit the pls model, generate the training set scores,
                     ## then attach what is needed to the random forest object to 
                     ## be used later
                     pre <- plsda(x, y, ncomp = param$ncomp)
                     scores <- pls:::predict.mvr(pre, x, type = "scores")
                     mod <- randomForest(scores, y, mtry = param$mtry, ...)
                     mod$projection <- pre$projection
                     mod
                   },
                   predict = function(modelFit, newdata, submodels = NULL) {       
                     scores <- as.matrix(newdata)  %*% modelFit$projection
                     predict(modelFit, scores)
                   },
                   prob = NULL,
                   varImp = NULL,
                   predictors = function(x, ...) rownames(x$projection),
                   levels = function(x) x$obsLevels,
                   sort = function(x) x[order(x[,1]),])

这是调用train:

 library(ChemometricsWithR)
 data(prostate)

 set.seed(1)
 inTrain <- createDataPartition(prostate.type, p = .90)
 trainX <-prostate[inTrain[[1]], ]
 trainY <- prostate.type[inTrain[[1]]]
 testX <-prostate[-inTrain[[1]], ]
 testY <- prostate.type[-inTrain[[1]]]

 ## These will take a while for these data
 set.seed(2)
 plsrf <- train(trainX, trainY, method = modelInfo,
                preProc = c("center", "scale"),
                tuneLength = 10,
                trControl = trainControl(method = "repeatedcv",
                                         repeats = 5))

 ## How does random forest do on its own?
 set.seed(2)
 rfOnly <- train(trainX, trainY, method = "rf",
                tuneLength = 10,
                trControl = trainControl(method = "repeatedcv",
                                         repeats = 5))

只是为了好玩，我得到了：

 > getTrainPerf(plsrf)
   TrainAccuracy TrainKappa method
 1     0.7940423    0.65879 custom
 > getTrainPerf(rfOnly)
   TrainAccuracy TrainKappa method
 1     0.7794082  0.6205322     rf

and

 > postResample(predict(plsrf, testX), testY)
  Accuracy     Kappa 
 0.7741935 0.6226087 
 > postResample(predict(rfOnly, testX), testY)
  Accuracy     Kappa 
 0.9032258 0.8353982 

Max