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

Question

我正在尝试弄清楚如何从 caret::train 线性模型中获取置信区间。

我的第一次尝试只是使用通常的 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 包只是对这个特定问题有一个糟糕的实现。在任何一种情况下，它似乎都适合在 their github 上提出问题或希望如果还不存在(希望更多样化的 predict 函数或修复 object$finalModel 中使用的命名)

问题(出现在第二次试验中)源于 caret 包如何在内部处理不同的拟合过程，基本上限制了似乎是清洁和标准化目的的预测功能。

问题:

问题是双重的。

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

这两个问题似乎源于train 的格式和predict.train 中的用法。首先关注后一个问题，通过查看

的输出可以明显看出这一点

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

显然，在运行 train 时执行了一些格式化，因为预期输出将是 mpg ~ poly(hp, 2)，而输出扩展了 RHS(和添加引号/标签)并更改了 LHS。因此，修复公式或能够使用公式会很好。

查看 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 用户，我们基本上看到，model.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

函数如下。对于 Nerd ，我基本上从 predict.train、predictionFunction 和 m$modelInfo$predict< 中提取了 model.matrix 构建过程。我不保证此函数适用于每个 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
}