泊松GLM纳入situatie后三组klasse的agressiekrab事后检验咨询

此前两种检验方法不适配/报错的核心原因 #

1. TukeyHSD运行报错的问题 #

• 表层报错原因：调用aov(agressiekrab~klasse*situatie)时，klasse和situatie以纯数值型传入模型，未被定义为分类因子，aov无法识别可开展比较的分组变量，因此抛出no factors in the fitted model错误。
• 方法本质不适配的深层原因：哪怕将变量转为因子后顺利运行TukeyHSD，该方法是专门为高斯正态线性模型（ANOVA）设计的事后比较逻辑，完全不适用于计数类非正态数据的泊松GLM分析框架，结果存在系统性偏差。

2. pairwise.wilcox.test结果无法解读的问题 #

• 分组逻辑完全偏离研究目标：将分组变量设为klasse*situatie的交互组合时，函数会自动拆分所有理论上的3个klasse×3个situatie细分交互组，输出表格里的1、2、3、4、6等编号是细分交互组的标签（受样本量限制，拆分后的交互组存在部分空组，编号不连续），并非计划比较的3个klasse主分组，自然无法解读。
• 分析框架不满足控制要求：威尔科克森秩和检验是无协变量控制的非参数检验，完全没有纳入situatie作为混杂因素，相当于把不同情境下的攻击数据混在一起直接比较组间差，结果会被情境效应污染，和“将situatie作为控制变量纳入分析”的核心目标完全背离。
• 细节参数问题：传入的p值校正参数P.adj="Bonj"是无效写法，函数无法识别后自动切换为默认的holm校正，不影响流程但属于写法错误。

适配当前分析需求的事后检验实现流程 #

前置关键修正（初始模型存在拟合偏差） #

初始拟合泊松模型时，klasse和situatie均为数值型变量，模型会自动将其识别为连续线性预测变量，估计的系数是“变量每升高1个单位的对数攻击数变化”，并非3个组别之间的分类效应，会导致后续所有解读完全偏离目标，必须先完成变量类型修正：

# 读取原始数据并将分类变量转为因子
dat <- data.frame(
  agressiekrab = c(8,5,1,10,6,12,4,17,1,1,5,9,11,2,3,0,21,17,4,1,10,4,14,15,22,8,19,0,6,16,4,10),
  klasse = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3),
  situatie = c(1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,1,1,1,2,2,2,1,1,1,1,1,2,2,2,2,2)
)
dat <- transform(
  dat,
  klasse = factor(klasse),
  situatie = factor(situatie)
)

完成类型转换后重新拟合正确的泊松GLM：

agmodel_correct <- glm(
  agressiekrab ~ klasse * situatie,
  family = poisson(link = "log"),
  data = dat
)
summary(agmodel_correct)

事后两两比较的R实现 #

使用emmeans包开展边际均值的两两比较，该方法完全适配泊松GLM框架，可在控制situatie效应的基础上完成组间检验，是目前计数类GLM事后比较的标准方案：

# 加载emmeans包
install.packages("emmeans")
library(emmeans)

# 计算边际均值（因模型存在显著的klasse:situatie交互，需分情境展示组间差）
em_res <- emmeans(agmodel_correct, ~ klasse | situatie)

# 对所有klasse分组做两两比较，自动采用Tukey方法校正多重检验的p值
pair_compare <- pairs(em_res, adjust = "tukey")
print(pair_compare)

结果解读逻辑 #

• 因模型中交互项显著（p=0.0054），说明不同情境下三组的攻击性差异规律并不统一，需按situatie的每个水平分别解读两两比较结果：输出会列出每个情境下，1组vs2组、1组vs3组、2组vs3组的对数计数差值、标准误、z值、校正后p值。
• 若需将对数尺度的差值转换为更易理解的发生率比（IRR），只需将比较结果中的估计值取指数运算（exp(pair_compare$estimate)），即代表某一组的攻击性事件计数是另一组的多少倍。
• 不建议跨情境计算无拆分的平均组间差，否则会掩盖交互带来的效应异质性，导致解读偏差。若需获得控制了situatie分布的总体组间比较，可将emmeans的公式修改为~ klasse，但需在报告中明确提及交互显著的前提。

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