如何基于分组数据对多变量执行分段回归并提取断点

问题背景 #

现有长格式数据，包含自变量x、分组变量Group、因变量标识Var和响应值y。需要对每个Group+Var的组合，以x为自变量、y为因变量执行分段回归，批量提取断点及其标准误，并整理成结构化数据框。

原手动处理代码如下：

library(segmented)

mod_lm <- lm(y ~ x, data = df) # 拟合线性模型
mod_seg <- segmented(mod_lm, seg.Z = ~ x) # 执行分段回归
mod_seg$psi # 提取断点及标准误

数据结构示例：

x   Group   Var     y
9   Group1  Var1    0.6901
6   Group1  Var1    0.6346
5   Group1  Var1    0.8089
5   Group1  Var1    0.1274
7   Group1  Var1    0.6426
1   Group1  Var2    0.1059
2   Group1  Var2    0.6989
4   Group1  Var2    0.1129
7   Group1  Var2    0.1458
7   Group1  Var2    0.8185
2   Group2  Var1    0.7950
0   Group2  Var1    0.0533
1   Group2  Var1    0.1866
3   Group2  Var1    0.3876
8   Group2  Var1    0.2788
2   Group2  Var2    0.1559
8   Group2  Var2    0.3382
1   Group2  Var2    0.6346
9   Group2  Var2    0.6038
8   Group2  Var2    0.2026

解决方案 #

结合dplyr分组功能与purrr映射函数，可实现批量处理：

步骤1：加载依赖包 #

library(segmented)
library(dplyr)
library(purrr)

步骤2：批量处理并提取结果 #

# 按Group和Var分组，嵌套数据后批量拟合模型
breakpoint_results <- df %>%
  group_by(Group, Var) %>%
  nest() %>%
  mutate(
    # 拟合线性模型
    lm_mod = map(data, ~lm(y ~ x, data = .x)),
    # 拟合分段回归
    seg_mod = map(lm_mod, ~segmented(.x, seg.Z = ~x)),
    # 提取断点及标准误，转换为数据框格式
    breakpoint_info = map(seg_mod, ~as.data.frame(.x$psi) %>%
                           rename(breakpoint = Est., se = St.Err.))
  ) %>%
  # 展开断点信息，合并分组标识
  unnest(breakpoint_info) %>%
  # 移除中间临时列
  select(-data, -lm_mod, -seg_mod)

# 查看最终结果
print(breakpoint_results)

结果说明 #

最终生成的breakpoint_results数据框包含4列：

• Group：分组标识
• Var：因变量标识
• breakpoint：分段回归的断点估计值
• se：断点估计的标准误

容错处理（可选） #

若部分分组因数据量不足或无分段趋势导致模型拟合失败，可添加错误捕获机制避免流程中断：

# 替换seg_mod的映射逻辑，遇到错误返回NA
seg_mod = map(lm_mod, ~possibly(segmented, otherwise = NA)(.x, seg.Z = ~x))

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