检验单因素对多分类分组对象效应的线性模型实现方法咨询

你原来用do的写法已经可以实现需求，不过目前tidyverse生态下有更易读、更易维护的实现方案，完全保留你需要的按species+gear分组检验的逻辑：

方案1：输出organization各水平的效应系数与显著性 #

library(tidyverse)
library(broom)

df %>%
  # 按物种、渔具分组，将分组内数据嵌套为列表列
  nest_by(species, gear) %>%
  # 分组拟合线性模型，整理效应结果
  mutate(
    model = list(lm(rate ~ organization, data = data)),
    effect_res = list(tidy(model))
  ) %>%
  # 展开效应结果，过滤所需内容
  unnest(effect_res) %>%
  mutate(p.value = round(p.value, 3)) %>%
  # 仅保留organization的效应，排除无意义的截距项，过滤显著结果
  filter(term != "(Intercept)", p.value < 0.05) %>%
  # 按需保留输出列
  select(species, gear, term, estimate, std.error, p.value)

优势说明：

• 替换了dplyr中已经不推荐使用的do函数，nest_by是目前官方推荐的分组嵌套建模方案，可读性更强
• 可以灵活保留完整的模型对象，后续需要提取拟合优度、残差等其他指标时不需要重构流程

方案2：仅检验organization整体显著性（更精简输出） #

如果你不需要看organization不同水平之间的效应差异，只需要判断每个species+gear分组下organization整体是否对渔获率有显著影响，可以用这个更精简的写法：

library(tidyverse)

df %>%
  nest_by(species, gear) %>%
  mutate(
    model = lm(rate ~ organization, data = data),
    # 提取organization的整体检验p值
    org_global_p = anova(model)$`Pr(>F)`[1]
  ) %>%
  # 过滤出organization效应显著的分组
  filter(org_global_p < 0.05) %>%
  select(species, gear, org_global_p)

额外说明 #

• 这套框架的兼容性极强，如果后续确认渔获率数据不符合线性模型假设（比如零值多、右偏），仅需要将lm替换为对应的模型函数（如glm、zeroinfl等）即可，整体流程不需要调整
• 所有输出天然保留species和gear的分组标识，不需要额外关联匹配

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