R语言做物种分布模型四：适生区预测与Boyce曲线评估

m0_51187291

335人浏览 · 2026-06-04 19:43:10

m0_51187291 · 2026-06-04 19:43:10 发布

简介

📖 上期回顾：在上一期推文中，我们使用 ENMeval 包进行了网格搜索调参，通过空间分块交叉验证（block partition）和 AICc 准则，从多种特征组合（fc）和正则化倍数（rm）中自动筛选出了最优模型。今天，我们继续深入——用调好的最优模型进行适生分布预测和模型评估。

本期以拟南芥（Arabidopsis thaliana）为例，使用上期选出的最优 maxnet 模型，完成以下工作：1）预测当前气候下的适生分布；2）绘制 Boyce 曲线定量评估模型预测能力；3）加载未来气候数据（SSP5-8.5 情景），预测未来适生分布。

🎯 为什么需要 Boyce 曲线？

仅靠 AUC 值评估模型是不够的——AUC 对"不存在"数据的选取方式敏感，且无法反映模型在适生梯度上的预测偏差（Lobo et al. 2008）。连续 Boyce 指数（Continuous Boyce Index, CBI）由 Boyce et al.（2002）提出，专门针对仅存在数据（presence-only）模型的评估：它计算预测适生度（HS）与观测/期望比（Predicted/Expected ratio）之间的 Spearman 秩相关系数。CBI 越接近 1，说明模型预测的适生度与物种实际分布越一致；CBI > 0 表示模型优于随机，CBI < 0 则表示模型表现不及随机预测（Hirzel et al. 2006）。

✉️ 关注本公众号，后台发送关键词"R语言做物种分布模型四"，

🎯 即可领取本期完整示例源代码和数据!

📦 一、加载所需R包

library(terra)
library(sp)
library(tictoc)
library(ENMeval)
library(maxnet)
library(dplyr)
library(ecospat)
library(ggplot2)

🗺️ 二、准备物种分布点与环境数据

本期沿用与上期相同的拟南芥分布数据和 WorldClim 气候因子。首先重新读取分布点并进行经纬度清洗，然后加载 19 个生物气候因子栅格。

大家可以参照第一期的内容，R语言做物种分布模型一：从GBIF获取物种分布数据，获取自己的物种分布数据进行分析。

气候数据可以参考R语言做物种分布模型二：根据经纬度进行气候数据提取，进行下载。

occ <- read.table(
"Arabidopsis thaliana/Arabidopsis thaliana-site.xls"
,

                  header = 
T
,

                  sep = 
'\t'
)


#清洗经纬度

occ <- occ %>%

  filter(

    Long >= -
180
 & Long <= 
180
,

    Lat  >= -
90
  & Lat  <= 
90

  ) %>%

  na.omit()


coordinates(occ) <- ~ Long + Lat

proj4string(occ) <- CRS(
"+proj=longlat +datum=WGS84"
)

occs <- as.data.frame(occ)[, c(
"Long"
, 
"Lat"
)]

colnames(occs) <- c(
"x"
, 
"y"
)


# 重新读取环境栅格（与之前相同路径，确保一致性）

env_files <- list.files(

"D:/桌面/jianzhi/code/气候因子抓取/worldclim/climate/wc2.1_2.5m"
,

  pattern = 
"\\.tif$"
,

  full.names = 
TRUE

)

env <- rast(env_files)

⚙️ 三、加载 ENMeval 结果并提取最优模型

上期的内容中，R语言做物种分布模型三：ENMeval高效调参，优化MaxEnt模型，我们已将 ENMevaluate() 的完整结果保存为 .rds 文件，现在直接加载。最优模型的选择标准是 AICc 最小，即 delta.AICc = 0 的参数组合。该组合在拟合优度与模型复杂度之间达到了最佳平衡。

eval <- readRDS(
"ENMevaluate_full_result.rds"
)

res <- eval@results


best_id <- which.min(res$AICc)

best_model <- eval@models[[best_id]]


# 保存最优模型，便于后续复用

saveRDS(

  best_model,

  file = paste0(
"best_maxnet_model_"
, best_id, 
"_cv.rds"
)

)

🔍 AICc（小样本校正赤池信息准则）同时考虑了模型似然度和参数数量，相比 AUC 等纯预测指标，更能避免选择过度参数化的模型（Warren & Seifert 2011）。在 SDM 研究中，基于 AICc 选择最优模型已成为标准做法（Muscarella et al. 2014）。

🌍 四、当前气候适生分布预测

有了最优模型和环境栅格，使用 terra::predict() 即可对整个研究区进行空间预测。type = "cloglog" 表示输出为互补双对数（complementary log-log）变换后的值，范围为 0–1，可解释为物种出现的相对概率。

💡 cloglog 与 logistic 的区别：MaxEnt 默认输出两种格式——logistic 和 cloglog。cloglog 是 Phillips et al.（2017）推荐的形式，它能更好地处理极低适生度区域的估计，且在生态学意义上更接近"出现概率"的直觉。

这里我们在32G的运行内存情况下，仅运行1分钟就得到了结果，相比之前从原数据调参将近30分钟的时间，是大大缩短了的。

tic(
"当前气候预测："
)

pred <- terra::predict(env, best_model, type = 
"cloglog"
, na.rm = 
T
)

toc()


plot(pred, main = 
"Current Suitability"
)


terra::writeRaster(

  pred,

  filename = paste0(
"current_suitability_best_model_"
, best_id, 
".tif"
),

  overwrite = 
TRUE

)

以下就是基于最优模型预测的当前气候情况下，拟南芥的适生区分布。颜色越偏向黄色，则表示出现概率越大，即适生度越高。

📊 五、Boyce 曲线评估模型预测能力

5.1 提取分布点的预测值

首先从当前预测栅格中提取已知分布点位置的适生度值：

pres_vals <- terra::extract(pred, occs)[, 
2
]

pres_vals <- pres_vals[is.finite(pres_vals)]  
# 去掉 NA 或 Inf


fit_vals <- values(pred)

fit_vals <- fit_vals[is.finite(fit_vals)]   
# 去掉 NA / Inf

5.2 计算连续 Boyce 指数

ecospat.boyce() 来自 ecospat 包（Di Cola et al. 2017），它通过将适生度划分为多个区间，计算每个区间内观测点与期望值的比值，再求比值与适生度均值的 Spearman 秩相关系数。

参数解释：

fit：整个研究区的预测适生度值（背景）。
obs：物种分布点处的预测适生度值。
nclass = 0：自动选择最优区间数。
window.w = "default"：移动窗口宽度，用于平滑比值曲线。
res = 100：分辨率参数，控制计算精度。

boyce <- ecospat.boyce(fit = fit_vals, obs = pres_vals,

                       nclass = 
0
, window.w = 
"default"
, res = 
100
)

5.3 绘制 Boyce 曲线

boyce_df <- data.frame(

  HS = boyce$HS,

  Fratio = boyce$F.ratio

)


p_boyce <- ggplot(boyce_df, aes(x = HS, y = Fratio)) +

  geom_line(linewidth = 
1.1
, color = 
"#2C7BB6"
) +

  geom_point(size = 
2
, color = 
"#2C7BB6"
) +

  geom_hline(yintercept = 
1
, linetype = 
"dashed"
, color = 
"grey40"
, size = 
1.2
) +

  labs(

    x = 
"Habitat suitability"
,

    y = 
"Predicted / Expected ratio"
,

    title = 
""

  ) +

  annotate(

    
"text"
,

    x = 
0.25
,

    y = max(boyce_df$Fratio, na.rm = 
TRUE
),

    label = paste0(
"Boyce = "
, round(boyce$cor, 
3
)),

    hjust = 
1
,

    vjust = 
1
,

    size = 
5

  ) +

  theme_bw() +

  theme(

    panel.grid = element_blank(),

    panel.border = element_rect(size = 
1.2
),

    legend.key.spacing.y = unit(
0.2
, 
'cm'
),

    legend.title = element_text(size = 
20
, face = 
'bold'
),

    axis.title = element_text(size = 
20
),

    legend.text = element_text(size = 
15
),

    axis.ticks = element_line(size = 
1
),

    axis.text = element_text(size = 
15
, color = 
'black'
))


p_boyce


ggsave(

  filename = 
"Boyce_curve_current.png"
,

  plot = p_boyce,

  width = 
6
,

  height = 
4
,

  dpi = 
300

)

📈 如何解读 Boyce 曲线？

连续Boyce指数衡量的是模型的预测与物种真实分布之间的一致性，值越接近1，说明模型性能越好。
横轴为适生度（Habitat suitability, HS），从低到高排列。
纵轴为 Predicted/Expected ratio：>1 表示该适生度区间内观测点多于随机期望，模型在该区间预测偏保守；<1 则预测偏高。
水平虚线（ratio = 1）为随机期望基线。

整体来讲，本模型的连续Boyce指数高达 0.975，预测/期望比随适宜性单调急剧上升（从 0.01 到 133），说明模型具有极强的排序能力，能几乎完美地分离极适宜与不适宜生境。

🔮 六、未来气候情景适生分布预测

接下来使用未来气候数据对最优模型进行投影。本示例采用 ACCESS-CM2 气候模式在 SSP5-8.5 情景下 2061–2080 年的预测数据。

future_files <- list.files(

  
"D:/桌面/jianzhi/code/气候因子抓取/future_climate/wc2.1_2.5m_ACCESS-CM2_ssp585_2061-2080/"
,

  pattern = 
"\\.tif$"
,

  full.names = 
TRUE

)


env_future_all <- rast(future_files)


# 统一变量名，使其与当前气候数据一致

future_bio_id <- gsub(
".*_(\\d+)$"
, 
"\\1"
, names(env_future_all))

std_names <- paste0(
"wc2.1_2.5m_bio_"
, future_bio_id)

names(env_future_all) <- std_names

⚠️ 跨时间投影的关键前提是变量名一致。未来气候栅格的波段名必须与模型训练时使用的变量名完全匹配，否则 predict() 会因找不到对应变量而报错。

tic(
"未来气候预测："
)

future_pred_ssp585 <- terra::predict(

  env_future_all,

  best_model,

  type = 
"cloglog"
,

  na.rm = 
TRUE

)

toc()


plot(future_pred_ssp585, main = 
"Future Suitability"
)


terra::writeRaster(

  future_pred_ssp585,

  filename = paste0(
"future_ssp585_suitability_best_model_"
, best_id, 
".tif"
),

  overwrite = 
TRUE

)

🔜 下期预告

模型预测做好了，可整张适生图该怎么解读？哪些区域是"高适生区"、哪些又是"不适生区"？

下期我们进行适生区划分与时空对比：

✅ 基于 5% FCV（Field Confirmed Value）阈值划分适生等级

✅ 将适生区分为不适生、低适生、中适生、高适生四级

✅ 绘制当前 vs 未来（SSP5-8.5）适生区对比图

✅ 统计各等级面积占比，直观展示气候变化对物种分布的影响

👉 关注更新，不迷路

📚 参考文献

Boyce, M. S., Vernier, P. R., Nielsen, S. E., & Schmiegelow, F. K. A. (2002). Evaluating resource selection functions. Ecological Modelling, 157(2–3), 281–300.
Di Cola, V., et al. (2017). ecospat: An R package to support spatial analyses and modeling of species niches and distributions. Ecography, 40(6), 774–787.
Hirzel, A. H., et al. (2006). Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling, 199(2), 142–152.
Lobo, J. M., Jiménez-Valverde, A., & Real, R. (2008). AUC: A misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography, 17(2), 145–151.
Muscarella, R., et al. (2014). ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution, 5(11), 1198–1205.
Phillips, S. J., et al. (2017). Opening the black box: An open-source release of Maxent. Ecography, 40(7), 887–893.
Warren, D. L., & Seifert, S. N. (2011). Ecological niche modeling in Maxent: The importance of model complexity and the performance of model selection criteria. Ecological Applications, 21(2), 335–342.