嘿，这个需求我之前做物种分布模型时也碰到过！核心思路很明确：要单独看每个环境变量对物种存在概率的影响，得把其他变量固定在平均水平（你用了scale()标准化，所以均值就是0），然后生成目标变量的全范围序列，再用模型预测概率和置信区间。下面给你两种实用方法，一种帮你理解原理，另一种高效省心：

方法一：手动构造预测数据（适合理解底层逻辑） #

这种方法一步步拆解，能让你清楚每一步在做什么：

1. 加载必备包 #

library(lme4)
library(dplyr)
library(ggplot2)

2. 定义要绘图的变量列表 #

vars_to_plot <- c("lotic_distance", "lentic_distance", "percent_impervious", "road_distance")

3. 写一个生成预测数据的函数 #

这个函数会自动处理单个变量的序列生成、概率预测和置信区间计算：

get_pred_data <- function(var_name, model, data) {
  # 提取目标变量的范围，生成100个均匀分布的点（保证曲线平滑）
  var_range <- range(data[[var_name]], na.rm = TRUE)
  var_seq <- seq(var_range[1], var_range[2], length.out = 100)
  
  # 构造预测数据集：目标变量取序列，其他变量设为0（标准化后的均值）
  # 随便选一个Municipality，后面用re.form=NA忽略随机效应的影响
  pred_df <- data.frame(
    lotic_distance = 0,
    lentic_distance = 0,
    percent_impervious = 0,
    road_distance = 0,
    Municipality = unique(data$Municipality)[1]
  )
  pred_df[[var_name]] <- var_seq
  
  # 在logit链接尺度预测（更适合计算置信区间）
  link_preds <- predict(model, newdata = pred_df, type = "link", re.form = NA, se.fit = TRUE)
  # 计算95%置信区间
  link_ci_lower <- link_preds$fit - 1.96 * link_preds$se.fit
  link_ci_upper <- link_preds$fit + 1.96 * link_preds$se.fit
  
  # 转换回概率尺度（用plogis函数）
  pred_df$prob <- plogis(link_preds$fit)
  pred_df$ci_lower <- plogis(link_ci_lower)
  pred_df$ci_upper <- plogis(link_ci_upper)
  
  # 标记变量名，方便后续批量绘图
  pred_df$variable <- var_name
  
  return(pred_df)
}

4. 批量生成所有变量的预测数据 #

all_pred_data <- bind_rows(lapply(vars_to_plot, get_pred_data, model = model, data = df))

5. 绘制带置信区间的曲线 #

用分面图展示每个变量的效应：

ggplot(all_pred_data, aes(x = .data[[variable]], y = prob)) +
  geom_line(size = 1, color = "#2c3e50") +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "#3498db") +
  facet_wrap(~variable, scales = "free_x") +  # 每个面板用独立的x轴范围
  labs(y = "物种存在概率", x = "标准化环境变量") +
  theme_bw() +
  theme(panel.grid = element_blank())

方法二：用emmeans包（高效省心） #

emmeans专门用来处理混合模型的边际效应计算，能省去手动构造数据的麻烦：

1. 加载包并计算边际效应 #

library(emmeans)

# 遍历每个变量，计算边际均值并转换到概率尺度
em_list <- lapply(vars_to_plot, function(var) {
  # 生成目标变量的序列，计算边际效应
  emmeans(model, specs = var, 
          at = setNames(list(seq(range(df[[var]], na.rm=T)[1], 
                                range(df[[var]], na.rm=T)[2], 
                                length.out=100)), var)) %>%
    regrid("response") %>%  # 从logit尺度转换为概率尺度
    as.data.frame() %>%
    mutate(variable = var)
})

all_em_data <- bind_rows(em_list)

2. 绘图（和方法一几乎一致） #

ggplot(all_em_data, aes(x = .data[[variable]], y = prob)) +
  geom_line(size = 1, color = "#2c3e50") +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = 0.2, fill = "#3498db") +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "物种存在概率", x = "标准化环境变量") +
  theme_bw() +
  theme(panel.grid = element_blank())

一些关键说明 #

• 为什么固定其他变量为0？：因为你用了scale()标准化，0代表该变量的均值，这样我们看的是当其他环境变量处于平均水平时，目标变量对物种存在概率的独立影响，这是单变量效应展示的标准方式。
• 关于随机效应：两种方法里的re.form=NA（手动方法）和emmeans默认行为，都是展示总体平均效应（整合了Municipality的随机差异）。如果想展示某个特定区域的效应，去掉re.form=NA并在预测数据里指定对应的Municipality即可。
• 置信区间的准确性：我们先在logit链接尺度计算CI再转概率，是因为二项式模型的误差在链接尺度更接近正态分布，这样得到的置信区间比直接在概率尺度计算更准确。

你可以根据自己的需求调整，比如修改序列长度、置信水平（把1.96换成2.576就是99%置信区间），或者调整绘图的颜色样式~