栅格和随机 Forest::predict() 中的 NA

Question

这里是新人，如果您需要更多信息，请告诉我。

我的目标：我正在使用 Rehfeldt 气候数据和 eBird 存在/不存在数据来生成使用随机森林模型的利基模型。

我的问题：我想预测整个北美的利基模型。 Rehfeldt 气候栅格具有大陆上每个单元格的数据值，但这些数据值被“海洋单元格”中的 NA 包围。请参阅情节 here，我将 NA 着色为深绿色。如果独立数据集包含 NA，则 randomForest::predict() 不会运行。因此，我想裁剪我的气候栅格（或设置工作范围？），以便 predict() 函数仅对包含数据的单元格进行操作。

疑难解答：

1. 我使用较小范围运行随机森林模型，其中不包括栅格的“NA 海洋”，并且模型运行良好。所以，我知道 NA 是问题所在。但是，我不想仅针对北美的一个矩形区域预测我的利基模型。

2. 我使用 flowla 的方法 here 使用北美的多边形 shapefile 来裁剪和遮罩栅格。我希望这会删除 NA，但事实并非如此。我可以做类似的事情来删除 NA 吗？

3. 我已经阅读了一些内容，但无法找到调整随机森林代码本身的方法，以便 predict() 忽略 NA。 This post 看起来很相关，但我不确定它对我的情况是否有帮助。

数据

我的栅格、输入存在/不存在文本文件和其他功能的代码是here。与下面的主要代码一起使用以获得可重现的示例。

代码

require(sp)
require(rgdal)
require(raster)
library(maptools)
library(mapproj)
library(dismo)
library(maps)
library(proj4)
data(stateMapEnv)

# This source code has all of the functions necessary for running the Random Forest models, as well as the code for the function detecting multi-collinearity
source("Functions.R")

# Read in Rehfeldt climate rasters
# these rasters were converted to .img and given WGS 84 projection in ArcGIS

d100 <- raster("d100.img")
dd0 <- raster("dd0.img")
dd5 <- raster("dd5.img")
fday <- raster("fday.img")
ffp <- raster("ffp.img")
gsdd5 <- raster("gsdd5.img")
gsp <- raster("gsp.img")
map <- raster("map.img")
mat <- raster("mat_tenths.img")
mmax <- raster("mmax_tenths.img")
mmin <- raster("mmin_tenths.img")
mmindd0 <- raster("mmindd0.img")
mtcm <- raster("mtcm_tenths.img")
mtwm <- raster("mtwm_tenths.img")
sday <- raster("sday.img")
smrpb <- raster("smrpb.img")

# add separate raster files into one big raster, with each file being a different layer.
rehfeldt <- addLayer(d100, dd0, dd5, fday, ffp, gsdd5, gsp, map, mat, mmax, mmin, mmindd0, mtcm, mtwm, sday, smrpb)

# plot some rasters to make sure everything worked
plot(d100)
plot(rehfeldt)

# read in presence/absence data
LAZB.INBUtemp <- read.table("LAZB.INBU.txt", header=T, sep = "\t")
colnames(LAZB.INBUtemp) <- c("Lat", "Long", "LAZB", "INBU")
LAZB.INBUtemp <- LAZB.INBUtemp[c(2,1,3,4)]
LAZB.INBU <- LAZB.INBUtemp

latpr <- (LAZB.INBU$Lat)
lonpr <- (LAZB.INBU$Long)
sites <- SpatialPoints(cbind(lonpr, latpr))
LAZB.INBU.spatial <- SpatialPointsDataFrame(sites, LAZB.INBU, match.ID=TRUE)

# The below function extracts raster values for each of the different layers for each of the eBird locations
pred <- raster::extract(rehfeldt, LAZB.INBU.spatial)
LAZB.INBU.spatial@data = data.frame(LAZB.INBU.spatial@data, pred)
LAZB.INBU.spatial@data <- na.omit(LAZB.INBU.spatial@data)

# ITERATIVE TEST FOR MULTI-COLINEARITY
# Determines which variables show multicolinearity
cl <- MultiColinear(LAZB.INBU.spatial@data[,7:ncol(LAZB.INBU.spatial@data)], p=0.05)
xdata <- LAZB.INBU.spatial@data[,7:ncol(LAZB.INBU.spatial@data)]  
for(l in cl) {
  cl.test <- xdata[,-which(names(xdata)==l)]
  print(paste("REMOVE VARIABLE", l, sep=": "))
  MultiColinear(cl.test, p=0.05)    
}

# REMOVE MULTI-COLINEAR VARIABLES
for(l in cl) { LAZB.INBU.spatial@data <- LAZB.INBU.spatial@data[,-which(names(LAZB.INBU.spatial@data)==l)] }


################################################################################################

# FOR LAZB
# RANDOM FOREST MODEL AND RASTER PREDICTION

require(randomForest)

# NUMBER OF BOOTSTRAP REPLICATES
b=1001

# CREATE X,Y DATA

# use column 3 for LAZB and 4 for INBU
ydata <- as.factor(LAZB.INBU.spatial@data[,3])
xdata <- LAZB.INBU.spatial@data[,7:ncol(LAZB.INBU.spatial@data)]

# PERCENT OF PRESENCE OBSERVATIONS
( dim(LAZB.INBU.spatial[LAZB.INBU.spatial$LAZB == 1, ])[1] / dim(LAZB.INBU.spatial)[1] ) * 100

# RUN RANDOM FORESTS MODEL SELECTION FUNCTION
# This model is using the model improvement ratio to select a final model.

pdf(file = "LAZB Random Forest Model Rehfeldt.pdf")
( rf.model <- rf.modelSel(x=xdata, y=ydata, imp.scale="mir", ntree=b) ) 
dev.off()

# RUN RANDOM FORESTS CLASS BALANCE BASED ON SELECTED VARIABLES
# This code would help in the case of imbalanced sample
mdata <- data.frame(y=ydata, xdata[,rf.model$SELVARS])
rf.BalModel <- rfClassBalance(mdata[,1], mdata[,2:ncol(mdata)], "y", ntree=b)

# CREATE NEW XDATA BASED ON SELECTED MODEL AND RUN FINAL RF MODEL
sel.vars <- rf.model$PARAMETERS[[3]]
rf.data <- data.frame(y=ydata, xdata[,sel.vars])  
write.table(rf.data, "rf.data.txt", sep = ",", row.names = F)

# This the code given to me; takes forever to run for my dataset (I haven't tried to let it finish)
# ( rf.final <- randomForest(y ~ ., data=rf.data, ntree=b, importance=TRUE, norm.votes=TRUE, proximity=TRUE) )

# I use this form because it's a lot faster
( rf.final <- randomForest(x = rf.data[2:6], y = rf.data$y, ntree=1000, importance=TRUE, norm.votes=TRUE, proximity=F) )

################################################################################################         
# MODEL VALIDATION 
# PREDICT TO VALIDATION DATA

# Determines the percent correctly classified
rf.pred <- predict(rf.final, rf.data[,2:ncol(rf.data)], type="response")
rf.prob <- as.data.frame(predict(rf.final, rf.data[,2:ncol(rf.data)], type="prob"))
ObsPred <- data.frame(cbind(Observed=as.numeric(as.character(ydata)), 
                            PRED=as.numeric(as.character(rf.pred)), Prob1=rf.prob[,2], 
                            Prob0=rf.prob[,1]) )
op <- (ObsPred$Observed == ObsPred$PRED)
( pcc <- (length(op[op == "TRUE"]) / length(op))*100 )

# PREDICT MODEL PROBABILITIES RASTER

# The first line of code says what directory I'm working, and then what folder in that directory has the raster files that I'm using to predict the range
# The second line defines the x variable, wich is my final Random Forest model

rpath=paste('~/YOURPATH', "example", sep="/")
xvars <- stack(paste(rpath, paste(rownames(rf.final$importance), "img", sep="."), sep="/"))
tr <-  blockSize(xvars)
s <- writeStart(xvars[[1]], filename=paste('~/YOURPATH', "prob_LAZB_Rehfeldt.img", sep="/"), overwrite=TRUE)                                           
for (i in 1:tr$n) {
  v <- getValuesBlock(xvars, row=tr$row[i], nrows=tr$nrows[i])
  v <- as.data.frame(v)         
  rf.pred <- predict(rf.final, v, type="prob")[,2]           
  writeValues(s, rf.pred, tr$row[i])
}
s <- writeStop(s)   

prob_LAZB <- raster("prob_LAZB_Rehfeldt.img")

# Write range prediction raster to .pdf
pdf(file="LAZB_range_pred.pdf")
plot(prob_LAZB)
map("state", add = TRUE)
dev.off()

谢谢！！

Answer 1

您是否尝试在对 RF 的调用中设置“na.action”？该选项在randomForest R manual 中明确标记。您对 RF 的调用如下所示：

rf.final <- randomForest(x = rf.data[2:6], y = rf.data$y, ntree=1000, importance=TRUE, norm.votes=TRUE, proximity=F, na.action = omit)

这将告诉 RF 忽略存在 NA 的行，从而丢弃这些观察结果。这不一定是最好的方法，但在您的情况下可能会很方便。

选项 2：rfImpute 或 na.roughfix：这将填写您的 NA，以便您继续进行预测。请注意，因为这可能会给您在估算/“修复” NA 的任何地方提供虚假预测。

选项 3：从选项 2 开始，在您获得预测后，将您的栅格带入您选择的 GIS/图像处理软件中，并屏蔽掉您不想要的区域。在您的情况下，掩盖水体将非常简单。