【发布时间】:2021-12-19 20:30:03
【问题描述】:
2021 年 11 月 9 日更新
由于这个问题,情节有点复杂。我的问题显然是 caret:: 包的一个怪癖。我之前的答案仅适用于运行估计,因为我可以让 S3 确定要使用 predict() 的哪个变体。但是,我想为我的线性模型计算 95% 的预测区间。为此,我必须使用predict() 函数的predict.lm() 变体。如果我使用lm() 函数创建模型,它工作正常,但是当我从使用caret::train() 构建的模型运行predict.lm() 时出现以下错误
Error in eval(predvars, data, env) :
object 'Sepal.Length:Sepal.Width' not found
这是一个更简短的可重现示例来说明问题。
##Loading Packages and Data##
library(caret)
data(iris)
####Model with an Interaction Term####
##Building a model Using lm()##
mod.lm<-lm(Petal.Length~Petal.Width+Sepal.Length*Sepal.Width, data=iris)
print(mod.lm)#Notice the interaction term is not quoted
class(mod.lm)
##Building a model using caret::train()##
trCtrl<-trainControl(method="LOOCV", savePredictions = TRUE)
mod.caret<-train(Petal.Length~Petal.Width+Sepal.Length*Sepal.Width, data=iris, method="lm", trControl=trCtrl)
print(mod.caret$finalModel)#Notice the interaction term is now quoted
class(mod.caret$finalModel)#Notice this is also of class lm
##Getting Prediction Intervals##
PI95.lm<-predict.lm(mod.lm, iris, interval="prediction") #No Error
PI95.caret<-predict.lm(mod.caret$finalModel, iris, interval="prediction")#This will throw an error##
####Model with only additive terms Term####
##Building a model Using lm()##
mod.lm.add<-lm(Petal.Length~Petal.Width+Sepal.Length+Sepal.Width, data=iris)
print(mod.lm.add)#Notice the interaction term is not quoted
class(mod.lm.add)
##Building a model using caret::train()##
trCtrl<-trainControl(method="LOOCV", savePredictions = TRUE)
mod.caret.add<-train(Petal.Length~Petal.Width+Sepal.Length+Sepal.Width, data=iris, method="lm", trControl=trCtrl)
print(mod.caret.add$finalModel)#Notice there are now no quotes
class(mod.caret.add$finalModel)#Notice this is also of class lm
##Getting Prediction Intervals##
PI95.lm.add<-predict.lm(mod.lm.add, iris, interval="prediction") #No Error
PI95.caret.add<-predict.lm(mod.caret.add$finalModel, iris, interval="prediction")#No more error because interaction term is gone
原帖
我正在使用caret:: 包来创建和交叉验证一个线性模型,然后我将使用该模型从预测值的栅格堆栈中预测一个属性。在我的真实数据集中,我需要使用交互项并执行幂变换。我在底部用我的假数据复制了这个问题。为此,我创建了一个链接函数,然后将其提供给raster::predict() 函数。建模工作正常。当我在训练数据的数据帧上运行链接函数时,它工作正常。但是,当我运行 raster:predict() 时,我收到以下错误消息
Error in eval(predvars, data, env) : object 'X3:X4' not found
此错误消息表明该函数正在尝试在我的栅格堆栈中查找名为 'X3:X4' 的图层,而不是创建变量 X3 和 X4 的交互,这两个变量是栅格堆栈中的单独图层。这个问题在下面使用我创建的一些假数据的可重现示例中得到了例证:
##Loading Necessary Packages##
library(caret)#For modeling and crossvalidation
library(raster)#To handle raster data
library(snow)#For Parallel processing
set.seed(88)#For Reproducibility
##Creating some fake data for model training##
X1<-runif(100, 20, 50)
X2<-runif(100, 100, 220)
X3<-runif(100, 80, 150)
X4<-runif(100, 1000, 15000)
Y<-rnorm(100, 2000, 50)
df<-data.frame(X1=X1, X2=X2, X3=X3, X4=X4, Y=Y)
##Creating a fake model##
trCtrl<-trainControl(method="LOOCV", savePredictions = TRUE)
mod<-train(Y^0.33~X1+X2+X3*X4, data=df, method="lm", trControl=trCtrl)
##Creating rasters with fake predictor data##
WGS84<-crs("+init=epsg:4326")
RAST<-raster(xmn=-122.0, xmx=-121.5, ymn = 45.0, ymx=45.5, crs=WGS84, ncol=100, nrow=100)
R1<-RAST
values(R1)<-runif(10000, 20, 50)
R2<-RAST
values(R2)<-runif(10000, 100, 220)
R3<-RAST
values(R3)<-runif(10000, 80, 150)
R4<-RAST
values(R4)<-runif(10000, 1000, 15000)
##Combining into a single raster stack and renaming with the same predictor variable names as the training dataset##
PRED<-stack(R1, R2, R3, R4)
names(PRED)<-c("X1", "X2", "X3", "X4")
##Creating a linking function to undo the power transformation in the model##
linkfun<-function(mod, x){
out<-predict(mod, x)^(1/0.33)
return(out)
}
##Starting Parallel Processing##
beginCluster()
##Attempting to predict a new raster of the response variable from the fake model##
raster::predict(PRED, mod$finalModel, fun=linkfun, filename=paste(tempdir(), "/Example.tif", sep=""), datatype="FLT4S", format="GTiff")
【问题讨论】: