通过biomod2整合多个模型对物种进行预测——白芷为例

Riceneeder Lv2
  2025-10-30 17:54:35 2025-10-30 17:54:35 创建   2025-10-30 18:05:10 2025-10-30 18:05:10 更新    1.8k 字  8 分钟  

前准备

项目文件结构

1
2
3
4
5
6
7
8
9
10
11
12
13
Angelica_BIOMOD2/
├── data/
│   └── Angelica_points.csv
├── env/                    # 环境bio文件
│   ├── bio1.tif
│   ├── bio2.tif
│   ├── soil_ph.tif
│   ├── soil_organic.tif
│   └── ...
├── maxent.jar
├── outputs/
│   └── Angelica_dahurica/
└── BIOMOD2_angelica.R
  • renv是 R 的项目环境管理工具，用于固定项目依赖包的版本，确保代码可复现。
  • renv::init()用于初始化项目环境
1
2
3
install.packages("renv")
library(renv)
renv::init()

1. 环境准备

1
2
3
4
5
6
7
8
9
10
renv::install(c("biomod2", "terra", "sp", "dismo", "ggplot2", "viridis",'tidyterra'))

library(terra)
library(sp)
library(dismo)
library(ggplot2)
library(viridis)
library(biomod2)

setwd("E:\\甘政坤\\Angelica_BIOMOD2")  # ← 改成你的项目路径

2. 数据导入

白芷分布点数据（presence only） CSV格式示例：

1
2
3
species,lon,lat
Angelica_dahurica,104.08,30.67
Angelica_dahurica,105.24,31.48
  • 读取 CSV 格式的分布点数据
  • 将数据转换为空间点对象（sp包格式），指定经度和纬度为坐标
  • 设置坐标参考系（CRS）为 WGS84（经纬度坐标系，全球通用）
1
2
3
occ <- read.csv("data/Angelica_points.csv")
coordinates(occ) <- ~lon + lat
proj4string(occ) <- CRS("+proj=longlat +datum=WGS84")
  • 读取env文件夹下所有.tif格式的环境变量文件（如 WorldClim 的 19 个生物气候变量、HWSD2 土壤变量等）
  • 将多个环境变量合并为一个raster栈（env_stack），方便后续统一处理
  • 打印环境变量名称，确认数据正确导入
  • 绘制第一个环境变量（如 Bio1，年平均温度）的空间分布，并叠加白芷分布点（红色小点），直观展示分布点与环境的关系
1
2
3
4
5
6
env_files <- list.files("env/", pattern = ".tif$", full.names = TRUE)
env_stack <- stack(env_files)
print(names(env_stack))

plot(env_stack[[1]], main="示例环境变量（Bio1）")
points(occ, col="red", pch=20)

3. 格式化输入数据

  • BIOMOD_FormatingData是 BIOMOD2 的核心函数，用于将分布数据和环境数据格式化为模型输入格式
  • 因原始数据只有存在点（presence-only），需生成伪缺失点（Pseudo-Absences, PA）作为模型的 “缺失” 样本（模型训练需要存在 / 缺失对比）
  • 这里设置生成 3 组伪缺失点，每组 10000 个，随机分布在研究区域
1
2
3
4
5
6
7
8
9
biomod_data <- BIOMOD_FormatingData(
  resp.var = rep(1, nrow(occ)),
  resp.xy = occ@coords,
  resp.name = "Angelica.dahurica",
  expl.var = env_stack,
  PA.nb.rep = 3,
  PA.nb.absences = 10000,
  PA.strategy = "random"
)

4. 模型配置

  • 定义 5 种常用物种分布模型的参数：
    • GLM（广义线性模型）：使用二次项，无交互项，以 AIC 准则选择变量
    • GAM（广义可加模型）：平滑参数 k=4，二项分布链接函数
    • RF（随机森林）：500 棵决策树
    • GBM（梯度提升机）：3000 棵树，学习率 0.005，控制过拟合
    • MAXENT（最大熵模型）：包含线性项、二次项、 hinge 项等特征
  • user.val整合所有模型的参数，allModels指定最终运行的模型列表
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# 定义模型参数
user.GLM <- list('_allData_allRun' = list(type = 'quadratic',
                                          interaction.level = 0,
                                          test = 'AIC'))
user.GAM <- list('_allData_allRun' = list(k = 4,
                                          family = binomial(link = 'logit')))
user.RF <- list('_allData_allRun' = list(ntree = 500))
user.GBM <- list('_allData_allRun' = list(n.trees = 3000,        # Increase trees
                                          interaction.depth = 3, # Keep this
                                          shrinkage = 0.005,     # Decrease learning rate
                                          n.minobsinnode = 10,     # Increase this slightly
                                          distribution = 'bernoulli'))
user.MAXENT <- list('_allData_allRun' = list(
  memory_allocated = NULL,
  initial_heap_size = NULL,
  max_heap_size = NULL,
  background_data_dir = NULL,
  visible = FALSE,
  linear = TRUE,
  quadratic = TRUE,
  product = FALSE,
  threshold = FALSE,
  hinge = TRUE,
  lq2lqptthreshold = 80,
  l2lqthreshold = 10,
  hingethreshold = 15,
  beta_threshold = -1,
  beta_categorical = 1,
  beta_lqp = -1,
  beta_hinge = -1,
  betamultiplier = 1,
  defaultprevalence = 0.5
))

user.val <- list(
  MAXENT.binary.MAXENT.MAXENT = user.MAXENT,
  GLM.binary.stats.glm = user.GLM,
  GAM.binary.mgcv.gam = user.GAM,
  RF.binary.randomForest.randomForest = user.RF,
  GBM.binary.gbm.gbm = user.GBM
)

allModels <- c( 'MAXENT','GLM', 'GAM', 'RF', 'GBM')

  • 用bm_ModelingOptions统一配置建模参数，指定数据类型、模型列表、参数来源等，为后续建模做准备。
1
2
3
4
5
6
7
8
myOptions <- bm_ModelingOptions(
  data.type = 'binary',
  models = allModels,
  strategy = 'user.defined',
  user.val = user.val,
  bm.format = biomod_data, 
  user.base = "bigboss"
)

5. 运行建模

  • BIOMOD_Modeling是模型训练的核心函数，基于输入数据和配置运行所有模型
  • 交叉验证（10 次重复，75% 训练 / 25% 验证）用于评估模型稳定性
  • 记录 TSS 和 ROC 作为模型性能指标，并计算变量重要性（评估环境变量对分布的影响程度）
1
2
3
4
5
6
7
8
9
10
11
biomod_model_out <- BIOMOD_Modeling(
  bm.format = biomod_data,
  modeling.id = 'Angelica.dahurica',
  models = allModels,
  OPT.user = myOptions,
  CV.strategy = 'random',
  CV.nb.rep = 10,
  CV.perc = 0.75,
  metric.eval = c('TSS', 'ROC'),
  var.import = 3
)

6. 模型评估

  • 提取模型评估结果（如各模型在训练集和验证集上的 TSS、ROC 值），用于判断模型表现（TSS>0.7、ROC>0.8 通常认为模型较好）
1
2
3
evals <- get_evaluations(biomod_model_out)
evals_summary <- BIOMOD_LoadModels(biomod_model_out)
print(evals)
  • 可视化模型评估结果：
    • 平均评估图：展示各模型在训练 / 验证集上的平均性能；
    • 箱线图：展示不同模型在多次重复中的性能分布，反映稳定性。
1
2
3
bm_PlotEvalMean(bm.out = biomod_model_out, dataset = 'calibration')
bm_PlotEvalMean(bm.out = biomod_model_out, dataset = 'validation')
bm_PlotEvalBoxplot(bm.out = biomod_model_out, group.by = c('algo', 'run'))

7. 变量重要性

  • 绘制变量重要性箱线图，展示不同环境变量在各模型（algo）和多次重复（run）中的重要性分布
  • 帮助识别对白芷分布影响最大的环境变量（如某个气候因子或土壤因子）
1
2
3
bm_PlotVarImpBoxplot(bm.out = biomod_model_out, group.by = c('expl.var', 'algo', 'algo'))
bm_PlotVarImpBoxplot(bm.out = biomod_model_out, group.by = c('expl.var', 'algo', 'run'))
bm_PlotVarImpBoxplot(bm.out = biomod_model_out, group.by = c('algo', 'expl.var', 'run'))

8. 集成建模

  • 集成多个单个模型的结果，减少单一模型的不确定性，提高预测可靠性
  • 先筛选出性能较好的模型（TSS≥0.7），再通过 3 种算法（平均、交叉验证平均、加权平均）生成集成模型
1
2
3
4
5
6
7
8
9
10
biomod_em <- BIOMOD_EnsembleModeling(
  bm.mod = biomod_model_out,
  models.chosen = 'all',
  em.by = 'all',
  metric.select = c('TSS'),
  metric.select.thresh = c(0.7),
  metric.eval = c('TSS', 'ROC'),
  var.import = 3,
  em.algo = c('EMmean', 'EMcv', 'EMwmean')
)

9. 集成预测

  • BIOMOD_Projection：用单个模型基于当前环境变量（env_stack）预测白芷的适宜性分布；
  • BIOMOD_EnsembleForecasting：基于集成模型生成最终的适宜性预测结果（综合多个模型的优势）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
biomod_proj <- BIOMOD_Projection(
  bm.mod = biomod_model_out,
  new.env = env_stack,
  proj.name = "current",
  models.chosen = 'all',
  metric.binary = 'TSS'
)

biomod_em_proj <- BIOMOD_EnsembleForecasting(
  bm.em = biomod_em,
  bm.proj = biomod_proj,
  proj.name = "currentEM",
  models.chosen = 'all',
  metric.binary = 'TSS'
)
1
plot(biomod_em_proj) # 快速可视化集成预测结果

10. 结果导出

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
output_path <- paste0(getwd(), "/outputs/Angelica_dahurica/proj_current/")
dir.create(output_path, recursive = TRUE, showWarnings = FALSE)

unpacked_val <- terra::unwrap(biomod_em_proj@proj.out@val)

# 保存结果
writeRaster(
  unpacked_val,
  filename = paste0(output_path, "Angelica_EMmeanByTSS.tif"),
  overwrite = TRUE
)
writeRaster(
  unpacked_val[[1]],
  filename = paste0(output_path, "Angelica_EMmeanByTSS_mean.tif"),
  overwrite = TRUE
)
writeRaster(
  unpacked_val[[2]],
  filename = paste0(output_path, "Angelica_EMmeanByTSS_cv.tif"),
  overwrite = TRUE
)
writeRaster(
  unpacked_val[[3]],
  filename = paste0(output_path, "Angelica_EMmeanByTSS_wmean.tif"),
  overwrite = TRUE
)

11. 可视化预测图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
pred_raster <- raster(paste0(output_path, "Angelica_EMmeanByTSS.tif"))
pred_df <- as.data.frame(rasterToPoints(pred_raster))
names(pred_df) <- c("x", "y", "suitability")

ggplot(pred_df, aes(x = x, y = y, fill = suitability)) +
  geom_raster() +
  scale_fill_viridis(name = "Habitat suitability", option = "C") +
  coord_equal() +
  theme_minimal() +
  labs(
    title = "Potential Suitable Distribution of Angelica dahurica",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.text = element_text(size = 10)
  )

最终效果：

  • 标题: 通过biomod2整合多个模型对物种进行预测——白芷为例
  • 作者: Riceneeder
  • 创建于 : 2025-10-30 17:54:35
  • 更新于 : 2025-10-30 18:05:10
  • 链接: https://gankun.cn.lu/posts/2025-10-30/
  • 版权声明: 版权所有 © Riceneeder，禁止转载。
评论
目录
通过biomod2整合多个模型对物种进行预测——白芷为例
  1. 前准备
    1. 1. 环境准备
    2. 2. 数据导入
    3. 3. 格式化输入数据
    4. 4. 模型配置
    5. 5. 运行建模
    6. 6. 模型评估
    7. 7. 变量重要性
    8. 8. 集成建模
    9. 9. 集成预测
    10. 10. 结果导出
    11. 11. 可视化预测图