【问题标题】:Error with `predict.lm()` on linear model built using `caret::train()` when interaction in model在模型中交互时,使用“caret::train()”构建的线性模型上的“predict.lm()”出错
【发布时间】: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")

【问题讨论】:

    标签: r r-caret


    【解决方案1】:

    正如我在给 Robert 的评论中所指出的,这只是将 raster::predict() 中的模型参数从 mod$finalModel 更改为 mod 的问题。我已经包含了 mod$finalModel,因为在我的工作流程的其他地方,我创建了一个 95% 预测区间的栅格,并且对 confint() 的调用需要一个 lm 对象,该对象在 @987654327 中作为 $finalModel 属性生成@。以下是更正后的代码:

    ##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+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, fun=linkfun, filename=paste(tempdir(), "/Example.tif", sep=""), datatype="FLT4S", format="GTiff")
    

    【讨论】:

      【解决方案2】:

      我认为这与 raster::predict 无关,因为我在使用 predict 和 data.frame 时也看到了该错误:

      p <- predict(mod$finalModel, data.frame(values(PRED)), fun=linkfun)
      #Error in eval(predvars, data, env) : object 'X3:X4' not found
      

      我可以使用这些值运行您的代码

      library(raster)
      library(caret)
      
      ## 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)
      
      ## model
      trCtrl <- trainControl(method="LOOCV", savePredictions = TRUE)
      mod <- train(Y^0.33~X1+X2+X3*X4, data=df, method="lm", rControl=trCtrl)
      
      ## predictor raster data
      PRED <- brick(xmn=-122.0, xmx=-121.5, ymn = 45.0, ymx=45.5, ncol=100, nrow=100, nl=4)
      values(PRED) <- runif(40000, 1000, 15000)
      names(PRED) <- c("X1", "X2", "X3", "X4")
      
      linkfun <- function(mod, x){
        predict(mod, x)^(1/0.33)
      }
      
      r <- raster::predict(PRED, mod, fun=linkfun)
      r
      #class      : RasterLayer 
      #dimensions : 100, 100, 10000  (nrow, ncol, ncell)
      #resolution : 0.005, 0.005  (x, y)
      #extent     : -122, -121.5, 45, 45.5  (xmin, xmax, ymin, ymax)
      #crs        : +proj=longlat +datum=WGS84 +no_defs 
      #source     : memory
      #names      : layer 
      #values     : 1.47405e-07, 2200.354  (min, max)
      

      【讨论】:

      • 谢谢罗伯特。顺便说一下,对于 raster 和 terra 包,这是一个重要的提示。它们非常有用,让我在笨重的 ESRI 环境中尝试进行地统计分析时节省了数小时的挫败感。事实证明这是一个“liveware”错误。我忘记了,如果你想创建预测区间,你只需要调用“mod$finalModel”。只需致电mod 即可解决此问题。我已经在下面发布了答案。
      猜你喜欢
      • 2021-08-30
      • 2016-04-24
      • 2021-06-16
      • 2012-05-16
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多