R中如何通过循环按Octagonplot分组拟合线性模型计算预测值

方案1：tidyverse 分组批量计算（推荐） #

这是R生态中处理分组建模最简洁的实现方式，无需手动编写循环，鲁棒性更强：

library(tidyverse)

result <- dat %>%
  # 按Octagonplot分组
  group_by(Octagonplot) %>%
  # 每组数据打包为嵌套列表
  nest() %>%
  mutate(
    # 对每组数据拟合线性模型，加异常处理避免拟合失败中断运行
    mod = map(data, ~possibly(lm, otherwise = NA)(d13C ~ midpoint, data = .x)),
    # 直接预测midpoint=5对应的d13C值
    pred_d13C = map_dbl(mod, ~ifelse(
      inherits(.x, "lm"), 
      predict(.x, newdata = data.frame(midpoint = 5)), 
      NA
    ))
  ) %>%
  # 只保留需要的字段
  select(Octagonplot, pred_d13C) %>%
  ungroup()

运行后result就是你需要的、和Octagonplot一一对应的结果数据框。

方案2：基础for循环实现（符合你要求的循环写法） #

如果你需要用显式循环实现，可参考以下代码：

# 提取所有唯一的Octagonplot编号
plot_ids <- unique(dat$Octagonplot)
# 初始化结果数据框
result <- data.frame(Octagonplot = integer(), pred_d13C = numeric())

# 遍历所有分组
for (id in plot_ids) {
  # 筛选当前分组的所有数据
  sub_dat <- dat[dat$Octagonplot == id, ]
  # 清除d13C为缺失值的行
  sub_dat_clean <- sub_dat[!is.na(sub_dat$d13C), ]
  
  # 有效观测数不足2条时无法拟合线性回归，结果赋值为NA
  if (nrow(sub_dat_clean) < 2) {
    pred_val <- NA
  } else {
    # 拟合线性模型
    mod <- lm(d13C ~ midpoint, data = sub_dat_clean)
    # 计算midpoint=5时的预测值
    pred_val <- predict(mod, newdata = data.frame(midpoint = 5))
  }
  
  # 将当前分组的结果写入结果表
  result <- rbind(result, data.frame(Octagonplot = id, pred_d13C = pred_val))
}

注意事项 #

• 部分分组的d13C存在缺失值，有效样本量不足2时无法拟合模型，结果会返回NA，你可以根据需求调整缺失值的预处理逻辑
• 直接用predict函数计算预测值比手动提取回归系数计算更稳妥，可避免系数顺序搞反的错误

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