R语言多站点分组分段回归断点提取与结果汇总方法

以下代码可直接跑通从分组拟合到结果输出的全流程，最终输出结构化断点表，支持直接计算平均恢复阈值。

依赖包加载 #

library(dplyr)
library(tidyr)
library(segmented)

分组拟合与断点提取 #

核心设置：固定分段回归仅识别1个断点，匹配你寻找"下降转增长/稳定"唯一临界位置的需求；npsi=1为强制单断点的核心参数，请勿删除；同时自动捕获拟合错误，提取斜率、拟合优度等辅助指标筛除无效拐点。

result_df <- df2 %>%
  # 按站点分组嵌套数据
  nest_by(site) %>%
  mutate(
    model_res = list(
      tryCatch({
        # 拟合基础线性模型
        lm_base <- lm(lam ~ ysw, data = data)
        # 拟合单断点分段回归
        seg_model <- segmented(lm_base, seg.Z = ~ysw, npsi = 1)
        # 提取核心指标
        bp_est <- seg_model$psi[2] # 断点ysw估计值
        bp_ci <- confint(seg_model) # 断点95%置信区间
        seg_slopes <- slope(seg_model)$ysw # 两段斜率
        model_r2 <- summary(seg_model)$r.squared # 模型拟合优度
        
        # 返回单站点结果
        tibble(
          ysw_breakpoint = bp_est,
          bp_ci_low = bp_ci[1],
          bp_ci_high = bp_ci[2],
          slope_before_bp = seg_slopes[1,1],
          slope_after_bp = seg_slopes[2,1],
          r_squared = model_r2,
          fit_status = "success"
        )
      }, error = function(e) {
        # 拟合失败时返回NA占位，记录错误原因
        tibble(
          ysw_breakpoint = NA_real_,
          bp_ci_low = NA_real_,
          bp_ci_high = NA_real_,
          slope_before_bp = NA_real_,
          slope_after_bp = NA_real_,
          r_squared = NA_real_,
          fit_status = paste0("failed: ", e$message)
        )
      })
    )
  ) %>%
  # 移除原始嵌套数据，解包拟合结果
  select(-data) %>%
  unnest(model_res) %>%
  ungroup()

有效断点筛选与平均阈值计算 #

分段回归可能识别出不符合"种群恢复"逻辑的假拐点，需要先过滤再计算平均阈值：

# 筛选有效恢复断点：拟合成功+断点前种群下降（斜率<0）+断点后种群增长/稳定（斜率≥0）+拟合质量达标
valid_breakpoints <- result_df %>%
  filter(
    fit_status == "success",
    slope_before_bp < 0,
    slope_after_bp >= 0,
    r_squared >= 0.1 # 可根据你的数据实际情况调整R²阈值
  )

# 计算多数站点种群恢复对应的平均ysw取值
mean_ysw_threshold <- mean(valid_breakpoints$ysw_breakpoint, na.rm = T)
# 可选：计算阈值的95%区间
threshold_95ci <- quantile(valid_breakpoints$ysw_breakpoint, c(0.025, 0.975), na.rm = T)

实操注意事项 #

• 提前筛除样本量不足的站点：单站点ysw维度的有效观测数少于6条时，分段回归结果可靠性极低，可在拟合前加summarise(n = n()) %>% filter(n >=6)步骤提前剔除
• 斜率阈值可灵活调整：如果你对"种群稳定"的定义是斜率大于-0.005而非严格≥0，直接修改slope_after_bp的过滤条件即可
• 拟合失败的站点可查看fit_status列的报错信息，针对性排查数据问题（比如ysw无梯度、lam无趋势等）

内容的提问来源于stack exchange，提问作者Amanda Goldberg