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

Question

这个问题是同一线程的延续 here .以下是本书中的一个最小工作示例:

Wehrens R. Chemometrics with R multivariate data analysis in the natural sciences and life sciences. 1st edition. Heidelberg; New York: Springer. 2011. (page 250).

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

宗旨:
使用相同的重复 CV 集执行已知策略 PLS 的交叉验证方法后面通常​​是 LDA或类似逻辑回归、SVM、C5.0、CART 的表亲，具有 caret 的精神包裹。因此每次在调用等待分类器之前都需要 PLS，以便对 PLS 分数空间而不是观察值本身进行分类。 caret 包中最近的方法是做 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

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

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

附注:
缩放和居中可能会非常困惑，我认为 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

如果我在本类(class)中遗漏了某些内容，请发表评论。谢谢。

Answer 1

如果你想用 caret 来适应这些类型的模型，您需要在 CRAN 上使用最新版本。创建最后更新以便人们可以使用 non-standard models他们认为合适。

我下面的方法是联合拟合 PLS 和其他模型(我在下面的示例中使用了随机森林)并同时调整它们。因此，对于每个折叠，ncomp 的 2D 网格和 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

和

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

最大限度