0

请原谅我的长问题，但我真的希望有人可以尝试帮助我改进我的代码。基本上这就是我想做的：用不同的输入重复相同的模型（例如随机森林）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)

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

4

1 回答 1

0

这至少是您使用和迭代列表的解决方案的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作为有用的工具。

编辑： PS请对任何不起作用的问题或部分发表评论，我会看看我是否可以在上面添加任何澄清或注释。

于 2019-12-16T21:16:32.007 回答