线性混合模型交互项缺失及显著性解读技术求助

1. 缺失参考组交互项的原因及解决方法 #

默认情况下，R对因子使用处理编码（treatment coding）：以第一个水平为参考组，模型中的交互项只会输出非参考组与参考组的斜率差异，参考组的交互效应被吸收进time主效应中，因此不会单独显示time:grup_int1。而nlme示例数据集能输出所有交互项，是因为它用了总和编码（sum coding）（如contr.sum），这种编码下每个组的交互项代表该组斜率与所有组平均斜率的差异，因此会显示所有水平的交互项。

两种解决方式： #

方式一：直接修改因子对比编码

无需循环relevel，直接设置总和编码即可输出所有交互项：

library(nlme)
# 设置总和编码（组数为3，所以用contr.sum(3)）
contrasts(genes_wide$grup_int) <- contr.sum(3)
# 重新拟合模型
lme_model <- lme(value ~ time * grup_int, random = ~1|id, data = genes_wide, na.action = na.exclude)
summary(lme_model)

方式二：循环relevel获取每组作为参考的对比结果

如果需要逐一将每个组设为参考，获取其他组与该组的斜率差异，可以用循环实现：

grup_levels <- levels(genes_wide$grup_int)
result_list <- list()

for (ref in grup_levels) {
  # 临时重新设置参考组
  genes_wide$grup_temp <- relevel(genes_wide$grup_int, ref = ref)
  # 拟合模型
  temp_model <- lme(value ~ time * grup_temp, random = ~1|id, data = genes_wide, na.action = na.exclude)
  # 提取包含time的系数（参考组斜率为time主效应，其他为组间斜率差异）
  temp_coef <- summary(temp_model)$tTable[grepl("time", rownames(summary(temp_model)$tTable)), ]
  result_list[[ref]] <- temp_coef
}

# 查看所有结果
print(result_list)

2. 交互项显著性的含义 #

交互项的显著性是当前组与参考组的time斜率差异的统计学检验结果：

• 当grup_int1为参考组时，time:grup_int2显著 → 组2的time斜率与组1存在显著差异；不显著 → 无足够证据证明两组斜率不同。
• 若你relevel后将其他组设为参考（比如组2），此时time:grup_int1显著 → 组1的斜率与组2存在显著差异，这和以组1为参考时time:grup_int2显著是同一结论（仅系数符号相反）。

如果出现time:grup_int1显著但time:grup_int2不显著，说明在当前参考组下，组1与参考组的斜率差异显著，而组2与参考组的斜率差异无统计学意义。

3. 获取组间斜率比较的整体p值 #

要检验所有组的time斜率是否存在整体差异（即交互效应的整体显著性），有三种常用方法：

方法1：似然比检验（模型比较） #

拟合包含/不包含交互项的两个模型，通过anova比较：

# 无交互项的空模型
model_null <- lme(value ~ time + grup_int, random = ~1|id, data = genes_wide, na.action = na.exclude)
# 包含交互项的全模型
model_full <- lme(value ~ time * grup_int, random = ~1|id, data = genes_wide, na.action = na.exclude)
# 似然比检验
anova(model_null, model_full)

输出的p-value即为交互效应的整体p值，代表是否存在至少两组的time斜率有显著差异。

方法2：直接提取模型的固定效应检验结果 #

对全模型使用anova()，可直接得到交互项的整体F检验p值：

anova(model_full)

结果中time:grup_int对应的p-value就是组间斜率差异的整体检验结果。

方法3：用emmeans包进行整体检验 #

若需要同时做事后两两比较，可使用emmeans的joint_tests()函数：

library(emmeans)
# 输出所有固定效应的整体检验结果
joint_tests(model_full)

该结果中time:grup_int对应的p值即为整体组间斜率差异的检验值。

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