问题描述
我正在使用 Biomod2 包在 R 中运行一系列物种分布模型。我使用的建模技术之一是使用 rpart 包的分类树分析 (CTA)。
这些模型中的响应是植物物种的存在/不存在，并且预测变量包含在 rasterStack 中。 rasterStack 中的大多数变量都是连续的数值变量，除了一个土地覆盖变量，geology，这是一个因素。我将每个单独的 rasterLayer 堆叠起来，然后使用 as.factor() 将地质层转换为一个因子。
我在尝试从 CTA 进行预测时遇到错误消息。 CTA 模型是用一个数据框构建的，其中“地质”是一个因素（见下文），并在 predict 上使用了栅格 rasterStack 函数（“地质”是一个因素，见下文）。但是，在运行 predict 函数时，我收到一条错误消息，指出我提供的是数字而不是因子。我已经检查了所有可能的点，看看是否以某种方式将“地质”转换回数字，但它似乎是我所看到的任何地方的一个因素（应该如此）。
编辑：更改数据以使其可重现。
library(raster)
library(rpart)
set.seed(123)
# Create sample rasterStack
data.rast <- stack(system.file("external/rlogo.Grd",package = "raster"))
# Create one layer as a factor
data.rast$geology <- as.factor(sampleInt(7,length(data.rast$red),replace = TRUE))
# Create sample presence/absence data by randomly selecting cells of raster
data <- as.data.frame(data.rast)
data <- data[sample(nrow(data),300,replace = FALSE),]
data$pa <- as.factor(sample(0:1,nrow(data),replace = TRUE))
names(data)[4] <- "geology"
head(data)
# red green blue geology pa
#2463 251 255 255 7 1
#1944 191 190 186 5 0
#5016 162 174 226 7 0
#5771 255 255 253 4 1
#3739 204 205 199 7 0
#5483 131 133 122 3 0
# Build CTA model using presence/absence dataframe
# Parameters set as the defaults in Biomod2 modeling options
cta <- rpart(pa ~ .,data = data,na.action = na.omit,method = "class",control = list(xval = 5,minbucket = 5,minsplit = 5,cp = 0.001,maxdepth = 25))
# Confirm classes of data before running predict function
data.frame(ctaClass = attr(terms(cta),"dataClasses")[2:5],rasterFactor = is.factor(data.rast))
# ctaClass rasterFactor
#red numeric FALSE
#green numeric FALSE
#blue numeric FALSE
#geology factor TRUE
# Once again confirming this rasterLayer is a factor
levels(data.rast$geology)
#[[1]]
# ID VALUE
#1 1 1
#2 2 2
#3 3 3
#4 4 4
#5 5 5
#6 6 6
#7 7 7
# Run predict function on rasterStack
cta.predict <- predict(object = data.rast,model = cta,type = "class")
#Error: variable 'geology' was fitted with type "factor" but type "numeric" was #supplied
#In addition: Warning message:
#In model.frame.default(Terms,newdata,na.action = na.action,xlev = #attr(object,:
# variable 'geology' is not a factor
编辑：添加证明它适用于 randomForests 模型
library(randomForest)
rf <- randomForest(pa ~ .,na.action = na.omit)
rf.predict <- predict(data.rast,rf)
rf.predict
#class : RasterLayer
#dimensions : 77,101,7777 (nrow,ncol,ncell)
#resolution : 1,1 (x,y)
#extent : 0,77 (xmin,xmax,ymin,ymax)
#crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs
#source : memory
#names : layer
#values : 0,1 (min,max)
#attributes :
# ID value
# 1 0
# 2 1
解决方法
在这种情况下，您需要通过提供因子名称和级别来帮助 predict
data$geology <- as.factor(data$geology)
cta.predict <- predict(data.rast,cta,type="class",factors=list(geology=levels(data$geology)))
还要注意 type= 中的 type=class，您不能只执行 class（除非您希望 filename 为 class.grd）
使用 terra 这会更好一些，我认为（希望）
library(terra)
x <- rast(data.rast*1)
x$geology <- as.factor(x$geology)
cta.predict <- predict(x,type="class")