从 caret::train 获取预测的置信区间
Getting confidence intervals on prediction from caret::train
我正在尝试找出如何从 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 获得线性模型拟合的置信区间？
免责声明：
这是一个糟糕的答案，或者 caret 包可能只是对这个特定问题的实现很糟糕。在任何一种情况下，如果不存在，它似乎适合打开一个问题或希望 their github（希望更多样化的 predict 功能或修复 object$finalModel 中使用的命名）
问题（出现在第二次试验中）源于 caret 包如何在内部处理不同的拟合程序，基本上限制了似乎是清洁和标准化目的的预测功能。
问题：
问题有两个。
predict.train 不允许 prediction/confidence 间隔train(...) 的输出中包含的 finalModel 包含格式异常的公式。
这两个问题似乎源于 train 的格式和 predict.train 中的用法。首先关注后一个问题，通过查看
的输出可以明显看出这一点
formula(m$finalModel)
#`.outcome ~ `poly(hp, 2)1` + `poly(hp, 2)2`)
显然在 运行 train 时执行了一些格式化，因为预期输出将是 mpg ~ poly(hp, 2)，而输出扩展了 RHS（并添加了 quotes/tags）并更改了 LHS。因此，修复公式或能够使用公式会很好。
查看 caret 包如何在 predict.train 函数中使用它揭示了下面的代码片段 newdata input
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
函数如下。对于书呆子，我基本上从 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
}
我正在尝试找出如何从 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 获得线性模型拟合的置信区间？
免责声明：
这是一个糟糕的答案，或者 caret 包可能只是对这个特定问题的实现很糟糕。在任何一种情况下，如果不存在，它似乎适合打开一个问题或希望 their github（希望更多样化的 predict 功能或修复 object$finalModel 中使用的命名）
问题（出现在第二次试验中）源于 caret 包如何在内部处理不同的拟合程序，基本上限制了似乎是清洁和标准化目的的预测功能。
问题：
问题有两个。
predict.train不允许 prediction/confidence 间隔train(...)的输出中包含的finalModel包含格式异常的公式。
这两个问题似乎源于 train 的格式和 predict.train 中的用法。首先关注后一个问题，通过查看
formula(m$finalModel)
#`.outcome ~ `poly(hp, 2)1` + `poly(hp, 2)2`)
显然在 运行 train 时执行了一些格式化，因为预期输出将是 mpg ~ poly(hp, 2)，而输出扩展了 RHS（并添加了 quotes/tags）并更改了 LHS。因此，修复公式或能够使用公式会很好。
查看 caret 包如何在 predict.train 函数中使用它揭示了下面的代码片段 newdata input
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
函数如下。对于书呆子，我基本上从 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
}