如何复制一个模型10次，并从中提取几个对象(测试结果)，然后计算平均值？

Question

请原谅我的长问题，但我真的希望有人能帮助我改进我的代码。基本上，这就是我想要做的:用不同的输入重复相同的模型(例如随机森林) 10次。作为每次迭代的结果，我想从每个模型中提取几个参数，并在所有迭代之后对它们进行均值和标准差(例如，均值AUC，均值偏差)。我可能会上传输入文件，但我的问题是连接到一个不直接依赖于它们的步骤，我认为可以使用一些编码来解决它。下面是一个示例：

我正在使用来自vignette随附的"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)

在这个例子中，a只做了2次重复，但理想情况下我想做10或100。如何让它变得更优雅和自动化，而不是手动创建100个对象。

Answer 1

这至少是使用purrr和dplyr并迭代列表的解决方案的一部分。这将带来将样本和结果存储在一个数据帧中的优势。

在下面的示例中，我使用了一个随机生成的数据帧和一个非常简单的函数来演示。我将在最后指出如何将其应用到您自己的数据和流程中。我还没有在上面的代码和数据上尝试过，因为它相当长和复杂，并且需要一段时间才能理解您的方法。但希望您能够看到如何将其融入到您自己的流程中。

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作为一个有用的工具进行迭代的过程中有一些有用的步骤。

编辑： P.S.请对任何问题或部件进行评论，我会看看是否可以在上面添加任何澄清或注释。