caret::train()构建含交互项线性模型调用predict.lm()报错如何解决

该问题实际是caret::包的特性导致的。此前的方案仅能完成常规估计，因为可以让S3自动匹配predict()的对应实现，但如果要为线性模型计算95%预测区间，必须调用predict.lm()版本的预测函数。使用lm()直接构建的模型调用该函数无异常，但使用caret::train()构建的模型调用时会抛出如下错误：

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::包构建并交叉验证线性模型，用于从预测变量的rasterstack中预测属性值。实际数据中需要使用交互项并进行幂次变换，模型训练过程正常，在训练数据构成的dataframe上调用自定义的链接函数也正常，但调用raster::predict()时会抛出如下错误：

Error in eval(predvars, data, env) : object 'X3:X4' not found 

该错误表明函数尝试在rasterstack中查找名为'X3:X4'的图层，而非自动对rasterstack中独立存在的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")

内容的提问来源于stack exchange，提问作者Sean McKenzie