线性混合效应模型手动构建交互项及ANCOMBC2适配问题

核心问题原因 #

你之前的错误在于同时将treatment、age和合并后的treatment_age放入模型，这三个变量存在完全共线性：treatment_age的每个水平可由前两者的水平完全推导，导致模型矩阵秩亏，进而出现系数丢弃、结果偏差，甚至ANCOMBC2的估计失败报错。

正确的手动构建与建模步骤 #

1. 规范生成交互项变量 #

先确保treatment和age是因子类型，用interaction()函数生成交互项（比paste更规范，能保留因子属性）：

# 转换为因子（若数据中还不是）
df$treatment <- as.factor(df$treatment)
df$age <- as.factor(df$age)

# 生成交互项因子，分隔符可自定义
df$treatment_age <- interaction(df$treatment, df$age, sep = "_")

2. 调整模型公式，避免共线性 #

不要同时放入主效应和交互因子，二选一即可：

• 方案1：直接使用交互因子建模（最简便，等价于treatment*age的全交互模型）
  公式只保留交互因子和随机效应：

  formula <- value ~ treatment_age + (1|year)
  

  这个模型会将每个treatment-age组合作为单独的系数，和treatment*age的模型结论完全一致，只是参数化方式不同（前者是每个组合相对参考组的差异，后者是主效应+交互效应）。

• 方案2：保留单个主效应+交互项（适合需要单独看某主效应的场景）
  先提取交互项的虚拟变量（排除参考组，避免共线性），再建模：

  # 提取treatment:age的交互项虚拟变量（去掉第一列参考组）
  int_vars <- model.matrix(~ treatment:age, data = df)[, -1]
  # 将交互项变量合并回原数据框
  df <- cbind(df, int_vars)
  
  # 模型公式：保留treatment主效应 + 交互项 + 随机效应
  formula <- value ~ treatment + treated_old + treated_new + (1|year)
  

  注：这里的treated_old、treated_new是int_vars生成的列名，需根据你的实际因子水平调整。

3. 解决协变量估计失败报错 #

除了共线性，样本量不足也会导致这个报错，先检查每个交互水平的样本量：

table(df$treatment_age)

若存在样本量为0或极少的水平，需合并该水平或剔除对应样本，再重新运行模型。

结果一致性验证 #

运行上述模型后，可将ANCOMBC2的结果与lmer(value ~ treatment*age + (1|year), data = df)的结果对比，两者的交互效应结论应完全一致，仅参数呈现方式不同。

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