R使用nlme包逐步构建线性混合效应模型及文献复现问题

R公式语法基础说明 #

首先明确两个符号的含义：

• : 仅代表两个变量的交互项，不包含对应主效应
• * 代表两个变量的主效应+交互项，例如MAT*MAP完全等价于MAT + MAP + MAT:MAP
• R公式会自动对重复输入的变量去重，你写的包含重复MAP的模型没有实际意义，R会自动忽略重复项。

问题1：变量重复、LRT检验逻辑错误 #

你当前的检验逻辑核心错误是：所有模型都和空模型对比，仅能证明「当前模型比空模型拟合效果好」，但无法证明「新加入的项能在已有变量的基础上带来额外的解释力」。
正确的嵌套检验逻辑应该是递进式的：

1. 以当前筛选得到的最优模型为基准，每次仅加入1个候选项
2. 对加入新项后的模型和基准模型做LRT检验，只有同时满足「AIC低于基准模型」且「LRT p值小于预设阈值（通常0.05）」才保留该新项
3. 保留新项后更新基准模型，再检验下一个候选项，直到当前组所有候选项检验完成
   按照这个逻辑，不会出现重复加变量的问题：比如已经在模型中保留了MAP，后续检验加MAT的时候，基准是已经有MAP的模型，LRT检验的就是MAT的额外贡献，不存在重复检验MAP的问题。

问题2：type变量的加入方法 #

根据原文献的分组递进逻辑，type属于功能控制变量，需要放在气候类变量筛选完成之后再加入：

1. 先完成第一组气候变量（主效应、二次项、交互项）的筛选，得到当前最优气候模型
2. 以最优气候模型为基准，加入type的主效应，做LRT检验，符合条件就保留
3. 若原方法要求检验类型和气候的交互，可在保留type主效应的基础上，逐一检验type和已保留气候变量的交互项，检验逻辑同上

修正后的筛选流程代码示例 #

library(nlme)

# 步骤1：初始化空模型作为初始基准
current_mod <- lme(C50_0.45 ~ 1, random = ~ 1|sitecode, data = practice, method = "ML")

# 步骤2：筛选第一组气候变量的一次主效应
candidate_main <- c("MAT", "MAP", "MGSL", "MGDD")
for (var in candidate_main) {
  test_mod <- update(current_mod, as.formula(paste0(".~. + ", var)))
  lrt_res <- anova(current_mod, test_mod)
  # 满足AIC降低且LRT显著则保留
  if (lrt_res$`p-value`[2] < 0.05 & lrt_res$AIC[2] < lrt_res$AIC[1]) {
    current_mod <- test_mod
  }
}

# 步骤3：筛选已保留主效应的二次项
retained_main <- labels(terms(current_mod))[-1] # 提取已保留的主效应，排除截距
for (var in retained_main) {
  test_mod <- update(current_mod, as.formula(paste0(".~. + I(", var, "^2)")))
  lrt_res <- anova(current_mod, test_mod)
  if (lrt_res$`p-value`[2] < 0.05 & lrt_res$AIC[2] < lrt_res$AIC[1]) {
    current_mod <- test_mod
  }
}

# 步骤4：筛选已保留主效应的一阶交互项
if (length(retained_main) >= 2) {
  var_combinations <- combn(retained_main, 2, simplify = FALSE)
  for (comb in var_combinations) {
    test_mod <- update(current_mod, as.formula(paste0(".~. + ", comb[1], "*", comb[2])))
    lrt_res <- anova(current_mod, test_mod)
    if (lrt_res$`p-value`[2] < 0.05 & lrt_res$AIC[2] < lrt_res$AIC[1]) {
      current_mod <- test_mod
    }
  }
}

# 步骤5：加入第二组控制变量（如type），重复检验逻辑
test_mod <- update(current_mod, .~. + type)
lrt_res <- anova(current_mod, test_mod)
if (lrt_res$`p-value`[2] < 0.05 & lrt_res$AIC[2] < lrt_res$AIC[1]) {
  current_mod <- test_mod
  # 如需检验type和气候变量的交互，可在此处继续按上述逻辑循环检验
}

# 最终得到最优模型后，换成REML方法拟合用于系数报告
final_mod <- update(current_mod, method = "REML")

注意事项 #

• 连续变量加入二次项前建议做标准化，降低多重共线性影响
• 交互项必须在对应两个主效应都保留在模型中的前提下才做检验
• 所有变量筛选阶段用ML方法拟合，最终报告系数用REML方法拟合

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