如何为各蜂箱线性模型批量应用segmented()函数获取断点？

你可以在已有的tidyverse管道中继续扩展，用map对每个蜂箱的线性模型批量应用segmented()函数，同时指定2个初始断点的猜测值。

完整代码实现 #

library(tidyverse)
library(segmented)

# 基于你已生成的lm.fit对象，添加分段模型拟合步骤
segmented_results <- lm.fit %>%
  mutate(
    # 对每个线性模型应用segmented，指定自变量和2个初始断点
    seg_model = map(lmMod, ~segmented(
      obj = ., 
      seg.Z = ~ndate, 
      psi = c(1620000000, 1630000000)  # 替换成你数据中合理的初始时间戳
    )),
    # 提取每个模型的最终断点估计值
    breakpoints = map(seg_model, ~.$psi[, "Est."]),
    # 提取分段模型的系数
    seg_coefficients = map(seg_model, coefficients)
  )

关键注意事项 #

• psi参数需要填入两个符合你ndate范围的初始时间戳，作为断点的初始猜测，segmented包会基于这个值迭代优化出最优断点。如果你的时间戳范围不同，一定要替换成对应区间的数值。
• 如果出现收敛报错，可以通过seg.control参数调整迭代设置，比如：

  seg_model = map(lmMod, ~segmented(
    obj = ., 
    seg.Z = ~ndate, 
    psi = c(1620000000, 1630000000),
    control = seg.control(maxit = 100, tol = 1e-4)
  ))
  

• 若要把断点整理成结构化数据框方便查看，可添加以下步骤：

  breakpoint_summary <- segmented_results %>%
    select(id, breakpoints) %>%
    unnest(breakpoints) %>%
    mutate(breakpoint_order = row_number(), .by = id)
  

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