R语言循环回归遇variable lengths differ错误及BMI结果提取求助

问题背景 #

需在R中对940个代谢物变量（对应数据框QBB_clean的第1653至2592列）分别拟合线性回归，自变量包含核心暴露因素bmi（连续）、连续变量Age及控制变量sex、lpa2c、smoking，目标是检验bmi对每个代谢物的影响。现有代码运行报错，同时需要提取每个回归中bmi的系数、p值、标准误、置信区间，并筛选显著结果。

报错代码 #

y<- c(1653:2592) # response 
x1<- c("bmi","Age", "sex","lpa2c", "smoking") # predictor 

for (i in x1){ 
  model <- lm(paste("y ~", i[[1]]), data= QBB_clean) 
  print(summary(model)) 
} 

错误信息 #

Error in model.frame.default(formula = paste("y ~", i[[1]]), data= QBB_clean, :
variable lengths differ (found for 'bmi').

数据示例 #

错误原因 #

1. 循环对象错误：代码循环的是自变量列表x1，但实际需要循环的是代谢物因变量的列索引/列名，且每个回归需包含所有自变量，而非逐个自变量单独拟合。
2. 因变量定义错误：y是列索引向量，直接写入公式会被当成长度940的独立向量，与数据框的行长度不匹配，触发"变量长度不一致"错误。

修正后的代码 #

方法1：基础循环实现 #

# 定义包含所有自变量的公式字符串
predictor_formula <- "bmi + Age + sex + lpa2c + smoking"

# 初始化结果存储数据框
results <- data.frame(
  metabolite = character(),
  coef_bmi = numeric(),
  se_bmi = numeric(),
  p_value = numeric(),
  ci_low = numeric(),
  ci_high = numeric(),
  stringsAsFactors = FALSE
)

# 循环每个代谢物列
for (col_idx in 1653:2592) {
  # 获取当前代谢物的列名
  met_name <- colnames(QBB_clean)[col_idx]
  # 构建完整回归公式
  formula_str <- paste(met_name, "~", predictor_formula)
  # 拟合线性模型
  model <- lm(formula_str, data = QBB_clean)
  # 提取bmi的统计结果
  model_summary <- summary(model)
  bmi_stats <- model_summary$coefficients["bmi", ]
  # 提取bmi的95%置信区间
  bmi_ci <- confint(model, "bmi")
  # 将结果存入数据框
  results <- rbind(results, data.frame(
    metabolite = met_name,
    coef_bmi = bmi_stats["Estimate"],
    se_bmi = bmi_stats["Std. Error"],
    p_value = bmi_stats["Pr(>|t|)"],
    ci_low = bmi_ci[1],
    ci_high = bmi_ci[2]
  ))
}

方法2：高效批量处理（purrr包） #

针对大样本量，推荐用purrr包避免循环rbind的性能损耗：

library(purrr)
library(dplyr)

# 提取所有代谢物的列名
metabolite_cols <- colnames(QBB_clean)[1653:2592]

# 批量拟合模型并提取结果
results <- map_dfr(metabolite_cols, function(met) {
  formula_str <- paste(met, "~ bmi + Age + sex + lpa2c + smoking")
  model <- lm(formula_str, data = QBB_clean)
  model_summary <- summary(model)
  bmi_stats <- model_summary$coefficients["bmi", ]
  bmi_ci <- confint(model, "bmi")
  
  tibble(
    metabolite = met,
    coef_bmi = bmi_stats["Estimate"],
    se_bmi = bmi_stats["Std. Error"],
    p_value = bmi_stats["Pr(>|t|)"],
    ci_low = bmi_ci[1],
    ci_high = bmi_ci[2]
  )
})

筛选显著结果 #

可按原始p值或校正后p值筛选（以Bonferroni校正为例）：

# 按原始p值<0.05筛选
significant_results <- results %>% filter(p_value < 0.05)

# 按Bonferroni校正后p值<0.05筛选（校正多重检验）
significant_results_adj <- results %>% 
  mutate(p_adj = p.adjust(p_value, method = "bonferroni")) %>%
  filter(p_adj < 0.05)

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