raster中的NAs和randomForest::predict()函数

3

新来的,请让我知道如果您需要更多信息。

我的目标:我正在使用Rehfeldt气候数据和eBird存在/不存在数据,使用Random Forest模型生成生态位模型。

我的问题:我想要预测整个北美洲的生态位模型。 Rehfeldt气候光栅具有整个大陆上每个单元格的数据值,但是这些被“海洋单元格”中的NAs所包围。请参见此处的图,其中我已将NAs涂成深绿色。如果独立数据集包含NAs,则randomForest :: predict()无法运行。因此,我希望裁剪我的气候光栅(或设置工作范围?),使predict()函数仅在包含数据的单元格上运行。

故障排除:

  1. 我已经使用不包括光栅的“NA oceans”的较小范围运行了Random Forest模型,并且该模型运行得非常好。因此,我知道NAs是问题。但是,我不想仅为北美洲的一个矩形块预测我的生态位模型。

  2. 我使用flowla的方法(此处)裁剪和遮罩了用于北美洲的多边形形状文件的光栅。我希望这将删除NAs,但它并没有。是否有类似的方法可以删除NAs?

  3. 我阅读了一些文章,但无法找到一种调整Random Forest代码本身的方法,以使predict()忽略NAs。这篇文章看起来很相关,但我不确定它是否有助于我的情况。

数据

我的光栅,输入存在/不存在的文本文件以及其他函数的代码在此处。与下面的主要代码一起使用以获得可重复的示例。

代码

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()

Thanks!!


你能添加一些数据以便重现你的问题吗? - Paulo E. Cardoso
嗨@Prophet60091,很抱歉这是一个非常老的问题,我认为它已经过时了。据我记得,函数的作者告诉我na.action = omit调用在当时不起作用,可能是基础R的问题?我不记得细节了。我同意na.action = omit 应该是解决方案,但我认为它对我来说失败了。我现在看到我没有清楚地表明我是否尝试了na.action参数;可能是我在写完这个问题后才尝试做到这一点...最终我采取了不同的分析方法,没有回到RF。 - LCM
1个回答

0

你尝试在调用随机森林函数 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:rfImputena.roughfix:这将填充您的NAs,以便您可以继续进行预测。注意,这可能会给您带来虚假的预测,因为NAs被插补/“修复”的位置。

选项3:从选项2开始,在获得预测结果后,将栅格图像导入您选择的GIS /图像处理软件中,并遮罩掉您不想要的区域。在您的情况下,遮罩水体将非常简单。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接