多因变量与自变量的多种回归模型循环实现求助

需求说明 #

处理包含多个因变量（flux系列）和自变量（expl系列）的数据集，为每个因变量搭配每个自变量遍历多种回归模型，通过broom::glance()输出结果，筛选显著p值并比较BIC。同时解决flux变量负值的对数转换问题，后续需加入management因子变量。

数据集 #

library(tidyverse)
library(broom)

dat <- structure(list(management = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", "c"
), class = "factor"), flux1 = c(61.51, 73.3, 66.74, 117.75, 74.69, 
74.12, 67.28, 78.67, 63.11, 82.74, 101.46, 80.48, 81.47, 56.76, 
74.34), flux2 = c(0.72, 8.93, -0.2, 3.9, 0.05, 1.39, -0.69, 2.19, 
3.31, 1.43, 8.57, 2.42, 2.56, 1.91, 10.41), flux3 = c(-52.97, 
-59, -52.75, -44.95, -52.96, -62.68, -73.13, -120.33, -82.07, 
-53.11, -36.01, -75.93, -84.01, -80.23, -127.61), expl1 = c(32.98, 
76.13, 70.19, 106.29, 118.66, 102.61, 70.77, 45.93, 74.63, 90.48, 
43.34, 98.7, 92.5, 73.68, 83.62), expl2 = c(8.9, 9, 9.03, 10.03, 
9.75, 9.85, 9.15, 9.28, 8.97, 9.25, 9.72, 9.32, 9.65, 8.82, 9.1
 ), expl3 = c(10.12, 11.62, 11.18, 12.9, 14.05, 13.93, 10.8, 10.85, 
 10.62, 11.97, 11.03, 7.38, 11.55, 12.1, 7.9)), row.names = c(NA, 
-15L), class = "data.frame")

现有代码 #

遍历单个因变量与所有自变量的单模型：

a <- dat %>% 
  gather(key = "expl", value = "expl_value", c(expl1:expl3)) %>%
  group_split(expl) %>% 
  map_df(~lm(log(flux1) ~ expl_value, data = .x) %>% 
           glance() %>% 
           mutate(expl_group = unique(.x$expl)))

多模型定义（未成功映射）：

models <- function(.x) {
  linear <- lm(flux1 ~ expl_value, data = .x)
  exponential <- lm(log(flux1 ) ~ expl_value, data = .x)
  logarithmic <- lm(flux1 ~ log(expl_value), data = .x)
  power <- lm(log(flux1) ~ log(expl_value), data = .x)
  lst(linear, exponential, logarithmic, power)
}

完整解决方案 #

1. 预处理：解决flux负值的对数转换问题

通过统一偏移量将所有flux值转为正数，避免对数转换报错：

# 计算所有flux变量的最小负值，生成偏移量（加0.01避免0值）
min_flux <- min(select(dat, starts_with("flux")))
offset <- abs(min_flux) + 0.01

# 生成可用于对数转换的flux列
dat_transformed <- dat %>%
  mutate(across(starts_with("flux"), ~ .x + offset, .names = "{.col}_log"))

2. 构建所有因变量-自变量组合

生成所有flux（含转换后）与expl的配对，整理为长格式：

# 生成变量配对表
var_pairs <- crossing(
  dep_var = str_subset(names(dat_transformed), "^flux\\d+$"),
  dep_var_log = str_subset(names(dat_transformed), "^flux\\d+_log$"),
  ind_var = str_subset(names(dat_transformed), "^expl\\d+$")
)

# 整理为长格式数据集，每行对应一组因变量-自变量
dat_long <- dat_transformed %>%
  pivot_longer(cols = starts_with("flux"), names_to = "dep_var", values_to = "dep_value") %>%
  pivot_longer(cols = starts_with("flux_"), names_to = "dep_var_log", values_to = "dep_value_log") %>%
  pivot_longer(cols = starts_with("expl"), names_to = "ind_var", values_to = "ind_value") %>%
  inner_join(var_pairs, by = c("dep_var", "dep_var_log", "ind_var")) %>%
  select(management, dep_var, dep_value, dep_value_log, ind_var, ind_value)

3. 定义多模型拟合函数

包含4种常用模型，支持可选加入management因子：

fit_models <- function(data, include_management = FALSE) {
  # 根据是否加入management生成公式列表
  base_formula <- if (include_management) {
    list(
      linear = dep_value ~ ind_value + management,
      exponential = dep_value_log ~ ind_value + management,
      logarithmic = dep_value ~ log(ind_value) + management,
      power = dep_value_log ~ log(ind_value) + management
    )
  } else {
    list(
      linear = dep_value ~ ind_value,
      exponential = dep_value_log ~ ind_value,
      logarithmic = dep_value ~ log(ind_value),
      power = dep_value_log ~ log(ind_value)
    )
  }
  
  # 拟合模型并提取关键统计量
  map_dfr(base_formula, ~glance(lm(.x, data = data)), .id = "model_type")
}

4. 批量运行模型并整理结果

按因变量-自变量分组，批量拟合所有模型：

model_results <- dat_long %>%
  group_by(dep_var, ind_var) %>%
  nest() %>%
  # 这里设置include_management = TRUE即可加入因子变量
  mutate(model_stats = map(data, fit_models, include_management = FALSE)) %>%
  unnest(model_stats) %>%
  ungroup() %>%
  # 筛选关键指标
  select(dep_var, ind_var, model_type, r.squared, adj.r.squared, p.value, bic)

5. 结果筛选与对比

• 筛选显著模型（p<0.05）并按BIC排序（BIC越小模型越优）：

significant_models <- model_results %>%
  filter(p.value < 0.05) %>%
  arrange(bic)

• 可视化不同模型的BIC差异：

ggplot(model_results, aes(x = model_type, y = bic, fill = dep_var)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~ind_var) +
  theme_minimal() +
  labs(title = "各模型BIC值对比", x = "模型类型", y = "贝叶斯信息准则(BIC)")

额外建议 #

• 加入management时，可通过relevel(management, ref = "A")指定基准因子水平，避免结果歧义
• 若flux负值较多，可尝试Box-Cox转换（MASS::boxcox()）选择最优转换方式，比固定偏移更灵活
• 模型评估除BIC外，建议结合残差图（broom::augment()）检查拟合效果，避免过拟合
• 可使用交叉验证（如caret包）评估模型泛化能力，提升结果可靠性

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