如何通过循环拟合含所有双向交互项的多个glm模型

1. 生成所有不重复的双向交互组合 #

先把23个变量的名字存成向量，方便后续生成交互组合：

vars <- paste0("x", 1:23)

用combn()函数生成所有两两不重复的变量对，返回的矩阵中每一列就是一组交互项：

interactions <- combn(vars, 2)

2. 定义主效应基础公式 #

先构建包含所有主效应的基础公式，后续所有模型都基于这个公式添加单个交互项：

base_formula <- as.formula(paste("y ~", paste(vars, collapse = " + ")))

3. 循环拟合模型并提取结果 #

用for循环实现（直观易读） #

# 初始化列表存储所有拟合好的模型
model_list <- list()
# 初始化数据框存储交互项的显著性结果
result_df <- data.frame(
  交互项 = character(),
  系数估计值 = numeric(),
  标准误 = numeric(),
  Z值 = numeric(),
  P值 = numeric(),
  stringsAsFactors = FALSE
)

for (i in 1:ncol(interactions)) {
  # 取出当前要添加的交互项
  int_pair <- interactions[, i]
  int_term <- paste(int_pair, collapse = ":") # 用:表示纯交互，避免重复添加主效应
  # 更新公式，把交互项加入基础模型
  full_formula <- update(base_formula, paste(". ~ . +", int_term))
  # 拟合模型，记得替换成你的数据集名称
  current_model <- glm(full_formula, data = your_data)
  # 将模型存入列表，用交互项命名方便后续查找
  model_list[[int_term]] <- current_model
  # 提取交互项的系数及显著性信息
  coef_info <- summary(current_model)$coefficients[int_term, ]
  # 将结果写入数据框
  result_df[i, ] <- list(
    交互项 = int_term,
    系数估计值 = coef_info[1],
    标准误 = coef_info[2],
    Z值 = coef_info[3],
    P值 = coef_info[4]
  )
}

用lapply简化写法（代码更简洁） #

如果习惯函数式编程，可以用lapply批量处理：

# 先生成所有带单个交互项的完整公式
all_formulas <- apply(interactions, 2, function(pair) {
  update(base_formula, paste(". ~ . +", paste(pair, collapse = ":")))
})

# 批量拟合模型
model_list <- lapply(all_formulas, glm, data = your_data)

# 批量提取结果到数据框
result_df <- do.call(rbind, lapply(names(model_list), function(term) {
  coef_info <- summary(model_list[[term]])$coefficients[term, ]
  data.frame(
    交互项 = term,
    系数估计值 = coef_info[1],
    标准误 = coef_info[2],
    Z值 = coef_info[3],
    P值 = coef_info[4],
    stringsAsFactors = FALSE
  )
}))

关键说明 #

• 用:而非*表示交互项：基础模型已经包含所有主效应，x1*x2等价于x1 + x2 + x1:x2，会重复加入主效应，虽然不影响结果，但用:更高效。
• 结果数据框result_df可直接筛选P值小于显著性水平（如0.05）的交互项，无需逐个查看模型摘要。
• 所有拟合好的模型都存在model_list中，如需查看某交互项的详细结果，直接调用summary(model_list[["x1:x2"]])即可。

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