从 caret::train 获取预测的置信区间

Question

我试图弄清楚如何从插入符::训练线性模型中获取置信区间。

我的第一次尝试只是使用通常的 lm 置信区间参数来运行预测：

m <- caret::train(mpg ~ poly(hp,2), data=mtcars, method="lm")
predict(m, newdata=mtcars, interval="confidence", level=0.95)

但看起来从 caret::train 返回的对象没有实现这个。

我的第二次尝试是提取最终模型并对此进行预测：

m <- caret::train(mpg ~ poly(hp,2), data=mtcars, method="lm")
fm <- m$finalModel
predict(fm, newdata=mtcars, interval="confidence", level=0.95)

但我得到了错误

Error in eval(predvars, data, env) : object 'poly(hp, 2)1' not found

深入挖掘，最终模型似乎对公式有一些奇怪的表示，并且正在我的新数据中搜索“poly(hp, 2)1”列，而不是评估公式。m$finalModel 看起来像这样：

Call:
lm(formula = .outcome ~ ., data = dat)

Coefficients:
   (Intercept)  `poly(hp, 2)1`  `poly(hp, 2)2`  
         20.09          -26.05           13.15

我应该补充一点，我不只是使用lm，因为我使用插入符号通过交叉验证来拟合模型。

如何通过 caret::train 从线性模型拟合中获取置信区间？

Answer 1

免责声明：

这是一个可怕的答案，或者也许这个caret包只是对这个特定问题有一个可怕的实现。在这两种情况下，如果尚不存在的话，似乎适合在他们的 github上提出问题或希望（要么希望有更多样化的predict功能，要么修复 中使用的命名object$finalModel）

问题（在第二次试验中出现）源于该caret包如何在内部处理不同的装配程序，基本上限制了预测功能，以达到清洁和标准化的目的。

问题：

问题有两个方面。

1. 不允许predict.train预测/置信区间
2. 的finalModel输出中包含的train(...)公式格式异常。

这两个问题似乎源于 的格式train和用法predict.train。首先关注后一个问题，通过查看输出可以明显看出这一点

formula(m$finalModel)
#`.outcome ~ `poly(hp, 2)1` + `poly(hp, 2)2`)

显然，在运行时执行了一些格式化train，正如预期的输出一样mpg ~ poly(hp, 2)，而输出扩展了右侧（并添加了引号/标签）并更改了左侧。因此，最好修复公式，或者能够使用该公式。

研究caret包如何在函数中使用它揭示了下面的输入predict.train代码片段newdata

predict.formula
#output
--more code
if (!is.null(newdata)) {
    if (inherits(object, "train.formula")) {
        newdata <- as.data.frame(newdata)
        rn <- row.names(newdata)
        Terms <- delete.response(object$terms)
        m <- model.frame(Terms, newdata, na.action = na.action, 
            xlev = object$xlevels)
        if (!is.null(cl <- attr(Terms, "dataClasses"))) 
            .checkMFClasses(cl, m)
        keep <- match(row.names(m), rn)
        newdata <- model.matrix(Terms, m, contrasts = object$contrasts)
        xint <- match("(Intercept)", colnames(newdata), 
            nomatch = 0)
        if (xint > 0) 
            newdata <- newdata[, -xint, drop = FALSE]
    }
}
--more code
    out <- predictionFunction(method = object$modelInfo, 
                modelFit = object$finalModel, newdata = newdata, 
                preProc = object$preProcess)

对于经验较少的R用户，我们基本上看到， amodel.matrix是从头开始构造的，而不使用 的输出formula(m$finalModel)（我们可以使用这个！），然后调用一些函数来基于 进行预测m$finalModel。从同一个包中查看predictionFunction发现该函数只是调用m$modelInfo$predict(m$finalModel, newdata)（对于我们的示例）

最后查看m$modelInfo$predict下面的代码片段

m$modelInfo$predict
#output
function(modelFit, newdata, submodels = NULL) {
                    if(!is.data.frame(newdata)) 
                        newdata <- as.data.frame(newdata)
                    predict(modelFit, newdata)
                  }

请注意，modelFit = m$finalModel和newdata是由上面的输出生成的。另请注意，对 的调用predict不允许指定interval = "confidence"，这是第一个问题的原因。

解决问题（有点）：

有多种方法可以解决这个问题。一是使用lm(...)而不是train(...). 另一种方法是利用函数的内部结构来创建一个数据对象，该对象符合奇怪的模型规范，因此我们可以predict(m$finalModel, newdata = newdata, interval = "confidence")按照预期的方式使用。

我选择做后者。

caretNewdata <- caretTrainNewdata(m, mtcars)
preds <- predict(m$finalModel, caretNewdata, interval = "confidence")
head(preds, 3)
#output
                         fit      lwr      upr
Mazda RX4           22.03708 20.74297 23.33119
Mazda RX4 Wag       22.03708 20.74297 23.33119
Datsun 710          24.21108 22.77257 25.64960

下面提供了该函数。对于书呆子来说，我基本上model.matrix从predict.train、predictionFunction和中提取了构建过程m$modelInfo$predict。我不会保证这个函数适用于每个caret模型的一般情况使用，但它是一个起点。

caretTrainNewdata功能：

caretTrainNewdata <- function(object, newdata, na.action = na.omit){
    if (!is.null(object$modelInfo$library)) 
        for (i in object$modelInfo$library) do.call("requireNamespaceQuietStop", 
                                                    list(package = i))
    if (!is.null(newdata)) {
        if (inherits(object, "train.formula")) {
            newdata <- as.data.frame(newdata)
            rn <- row.names(newdata)
            Terms <- delete.response(object$terms)
            m <- model.frame(Terms, newdata, na.action = na.action, 
                             xlev = object$xlevels)
            if (!is.null(cl <- attr(Terms, "dataClasses"))) 
                .checkMFClasses(cl, m)
            keep <- match(row.names(m), rn)
            newdata <- model.matrix(Terms, m, contrasts = object$contrasts)
            xint <- match("(Intercept)", colnames(newdata), 
                          nomatch = 0)
            if (xint > 0) 
                newdata <- newdata[, -xint, drop = FALSE]
        }
    }
    else if (object$control$method != "oob") {
        if (!is.null(object$trainingData)) {
            if (object$method == "pam") {
                newdata <- object$finalModel$xData
            }
            else {
                newdata <- object$trainingData
                newdata$.outcome <- NULL
                if ("train.formula" %in% class(object) && 
                    any(unlist(lapply(newdata, is.factor)))) {
                    newdata <- model.matrix(~., data = newdata)[, 
                                                                -1]
                    newdata <- as.data.frame(newdata)
                }
            }
        }
        else stop("please specify data via newdata")
    } else
        stop("please specify data data via newdata")
    if ("xNames" %in% names(object$finalModel) & is.null(object$preProcess$method$pca) & 
        is.null(object$preProcess$method$ica)) 
        newdata <- newdata[, colnames(newdata) %in% object$finalModel$xNames, 
                           drop = FALSE]
    if(!is.null(object$preProcess))
       newdata <- predict(preProc, newdata)
    if(!is.data.frame(newdata) && 
      !is.null(object$modelInfo$predict) && 
      any(grepl("as.data.frame", as.character(body(object$modelInfo$predict)))))
           newdata <- as.data.frame(newdata)
    newdata
}