请问我一个很长的问题,但是我真的希望有人可以尝试帮助我改善我的代码。基本上,这就是我想做的事情:用不同的输入重复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个对象。?
这里至少是使用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
作为有用的工具进行迭代时需要一些有用的步骤。
编辑:附言请针对任何无法解决的问题或部分发表评论,我将看看是否可以在上面添加任何说明或注释。