如何复制模型10次并从中提取几个对象(测试结果),然后求平均值?

问题描述 投票:0回答:1

请问我一个很长的问题,但是我真的希望有人可以尝试帮助我改善我的代码。基本上,这就是我想做的事情:用不同的输入重复10次相同的模型(例如随机森林)。每次迭代的结果是,我想从每个模型中提取几个参数,并在所有迭代之后根据它们得出均值和标准差(例如均值AUC,均值偏差)。我可以上载输入文件,但是我的问题与不直接依赖于它们的步骤有关,我想可以使用一些编码来解决。这是一个例子:

我正在使用附带“ dismo”包的小插图中的数据来处理物种分布模型。所有代码都可以在这里找到:https://rspatial.org/raster/sdm/6_sdm_methods.html#random-forest首先,我要创建一个物种出现率(pb = 1)和伪出现率(pb = 0)的数据。它们在两列中伴随着经度和纬度坐标,随后的环境变量被连接到每个点。此处一切正常,因此我可以创建一个模型。但我想建立几个模型并平均其结果。

这些是我的初始步骤:

require(raster)
#that is my file with occurrence points:
points_herb <- read.csv("herbarium.csv",header=TRUE)
points_herb <- points_herb[,2:3]
points_herb <- SpatialPointsDataFrame(coords = points_herb, data = points_herb, proj4string + CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
> head(points_herb)
 lon_x    lat_y
1 19.62083 49.62917
2 19.64583 49.62917
3 20.23750 49.61250...

#Variables (I use variables from PCA ran on climate)
files <- list.files("D:/variables/",pattern='asc',full.names=TRUE)
predictors <- raster::stack(files)
> predictors
class      : RasterStack 
dimensions : 1026, 1401, 1437426, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 16.36667, 28.04167, 42.7, 51.25  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : PCA1, PCA2

#Assigning variables to points
presvals <- extract(predictors, points_herb)
reading background points (about 20000):
points_back <- read.csv("back.csv",header=TRUE,dec = ".",sep = ",")
points_back <- points_back[,2:3]
points_back <- SpatialPointsDataFrame(coords = points_back, data = points_back, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

assigning variables to background/pseudoabsence points
absvals <- extract(predictors, points_back)
absvals <- unique(absvals)

#**this is important!** Sampling 1000 random points from my entire dataset containing ca. 20000
absvals_1 <- absvals[sample(nrow(absvals), 1000), ]

#making an input file for the modeling
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals_1)))
sdmdata1 <- data.frame(cbind(pb, rbind(presvals, absvals_1)))
sdmdata1 <- na.omit(sdmdata1)```

> head(sdmdata1)
  pb   PCA1    PCA2 
1  1  9.985359 2.419048 
2  1  8.711462 2.229476 
...

我运行模型:

#Random Forest
library(dismo)
library(randomForest)
#rf1- first random forest model
model_rf1 <- pb ~ PCA1 + PCA2
bc <- randomForest(model_rf1, data=sdmdata1)
#the model is predicted over a geographic space
bc_mod <- predict(predictors, bc, progress='')

#let's test it using CalibratR
require(CalibratR)
#extracting model probabilities to presence and absence points (those are actually from a separate dataset)
points_pres1 <- extract(bc_mod, points_pres1, cellnumbers=TRUE)
points_abs1 <- extract(bc_mod, points_abs1, cellnumbers=TRUE)
#prepare those data to test the model 
testECE <- c(rep(1, nrow(points_pres1)), rep(0, nrow(points_abs1)))
testECE <- data.frame(cbind(testECE, rbind(points_pres1, points_abs1)))
testECE <- na.omit(testECE)
testECE <- subset(testECE, select = c(testECE, layer))
#make Expected Calibration Error
ECE <- getECE(testECE$testECE, testECE$layer, n_bins = 10)
#make Maximum Calibration Error
MCE <- getMCE(testECE$testECE, testECE$layer, n_bins = 10)
#some other test
require(Metrics)
#get RMSE values
RMSE <- rmse(testECE$testECE, testECE$layer)

random_forest_1 <- data.frame(mget(c('ECE', 'RMSE', 'MCE')))
rownames(random_forest_1) <- "random_forest1"

然后,我想运行相同的模型,但使用不同的背景点。因此,在这种情况下,我将创建另一个输入文件,并从整个数据集中获取另外1000个随机点:

absvals_2 <- absvals[sample(nrow(absvals), 1000), ]
pb <- c(rep(1, nrow(presvals_2)), rep(0, nrow(absvals_2)))
sdmdata2 <- data.frame(cbind(pb, rbind(presvals_2, absvals_2)))
sdmdata2 <- na.omit(sdmdata2)

model_rf2 <- pb ~ variable1 + variable2
bc <- randomForest(model_rf2, data=sdmdata2)
bc_mod <- predict(predictors, bc, progress='')

#again, let's test it using CalibratR
points_pres2 <- extract(bc_mod, points_pres2, cellnumbers=TRUE)
points_abs2 <- extract(bc_mod, points_abs2, cellnumbers=TRUE)
# everything just as above, the objects are overwritten
testECE <- c(rep(1, nrow(points_pres2)), rep(0, nrow(points_abs2)))
testECE <- data.frame(cbind(testECE, rbind(points_pres2, points_abs2)))
testECE <- na.omit(testECE)
testECE <- subset(testECE, select = c(testECE, layer))
ECE <- getECE(testECE$testECE, testECE$layer, n_bins = 10)
MCE <- getMCE(testECE$testECE, testECE$layer, n_bins = 10)
RMSE <- rmse(testECE$testECE, testECE$layer)

random_forest_2 <- data.frame(mget(c('ECE', 'RMSE', 'MCE')))
rownames(random_forest_2) <- "random_forest2"

#And finally let's make a mean from ECE, MCE, and RMSE
rf_results <- rbind(random_forest_1, random_forest_2)
rf_results_mean <- sapply(rf_results, 2, FUN=mean)
#and standard deviation
rf_results_sd <- sapply(rf_results, 2, FUN=sd)

result <- rbind(rf_results_mean, rf_results_sd)

在此示例中,仅进行了2次重复,但是理想情况下,我想进行10或100次重复。如何使其更优雅,更自动,而不是手动创建100个对象。?

r repeat random-forest spatial dismo
1个回答
0
投票

这里至少是使用purrrdplyr并遍历列表的解决方案的一部分。这将具有将样本和结果存储在一个数据框中的优势。

在下面的示例中,我使用了一个随机生成的数据框和一个非常简单的函数来演示。最后,我将指出如何将其适合您自己的数据和流程。我没有在上面的代码和数据上尝试过此操作,因为它很长而且很复杂,需要花费一些时间来了解您的方法。但是希望您能够看到如何将其应用于您自己的过程。

library(dplyr)
library(purrr)

# step 1: create a function that takes a dataframe and returns a dataframe
calculate_mean_sd <- function(df){
  tibble(
    mean_lat = mean(df$lat),
    sd_lad = sd(df$lat),
    mean_long = mean(df$long),
    sd_long = sd(df$long)
  )
}

# random dataframe with all values you'd want to use (i.e. your `absvals` above)
full_df <- tibble(
  id = 1:100000,
  lat = runif(100000, 0, 100),
  long = runif(100000, 0, 100)
)

# step 2: create an empty list with the number (100) of loops you want to do
df <- as.list(1:100) %>% 
  map(~ tibble(iteration = .x))  # makes the iteration number into a dataframe to use later

# step 3: for each of 100 rows take a sample of a specified size and add to list as a dataframe
samples <- df %>%
  map( ~ mutate(.x, sample = list(full_df[sample(nrow(full_df), 100),])))

# step 4: iterate over list and pass your dataframes to the function, add results to new column
results <- samples %>% 
  map_df( ~ mutate(.x, results = list(bind_cols(.x[1], calculate_mean_sd(.x$sample[[1]]))))) 

# final, optional step: output a dataframe with iteration labelled and results
results$results %>% 
  reduce(bind_rows)

在上面的数据中,您想使用absvals[sample(nrow(absvals), 1000), ]在步骤3中对数据进行采样,然后将其后的部分放入函数中,该函数将返回包含所需列的数据框。

这可能无法给出完整的答案,但是希望在使用purrr作为有用的工具进行迭代时需要一些有用的步骤。

编辑:附言请针对任何无法解决的问题或部分发表评论,我将看看是否可以在上面添加任何说明或注释。

© www.soinside.com 2019 - 2024. All rights reserved.