基于Cohort对数据框多列做线性回归时遇NA/NaN/Inf错误求助

我在StackOverflow搜过相关答案但还是搞不懂，若问题重复先说声抱歉。我有如下样例数据框：

df <- data.frame(Cohort = c('con', 'con', 'dis', 'dis', 'con', 'dis'),
                 Sex = c('M', 'F', 'M', 'F', 'M', 'M'),
                 P1 = c(50, 40, 70, 80, 45, 75),
                 P2 = c(10, 9, 15, 13, 10, 8))

我想以Cohort为预测变量，对所有数值列做线性回归（之后还要加Sex这类特征）。我先剔除无关列：

new_df <- df[,-c(Sex)]

再执行回归：

fit <- lapply(new_df[-1], function(y){summary(lm(y ~ Cohort, data=new_df))})

小数据集（约5列）测试没问题，但真实数据集（约7300列）报错：

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
NA/NaN/Inf in 'y'

我以为是缺失值，但sum(is.na(new_df))返回0，加na.action=na.omit也没用。我的目标是提取p值（anova(fit)$'Pr(>F)'）和R平方（summary(fit)$r.squared）。请问怎么修正错误？有没有更优方法？另外后续加特征时，能不能不用子集化数据框就执行回归？

附真实数据框片段（dput输出）：

structure(list(Cohort = c("Disease", "Disease", "Control", "Control"), 
    seq.10010.10 = c(8.33449676839042, 8.39959836912012, 8.34385193344212, 
    8.43546191447928), seq.10011.65 = c(11.5222872738433, 11.7652860987237, 
    11.1661630826461, 11.008848763327), seq.10012.5 = c(10.5414838640543, 
    10.6862378767518, 10.5408061105915, 10.726558779105)), class = c("soma_adat", 
"data.frame"), row.names = c("258633854330_1", "258633854330_3", 
"258633854330_5", "258633854330_6"))

1. 修正「NA/NaN/Inf in 'y'」错误 #

sum(is.na(new_df))只能检测缺失值，但无法识别Inf或-Inf。可以检查数值列中是否存在无穷值：

# 检查所有数值列是否有Inf/-Inf
any_inf <- sapply(new_df[sapply(new_df, is.numeric)], function(col) any(is.infinite(col)))
# 输出有问题的列名
names(which(any_inf))

如果找到这类列，可按以下方式处理：

# 将Inf/-Inf替换为NA，再删除对应行（或根据需求替换为合理值）
new_df[sapply(new_df, is.numeric)] <- lapply(new_df[sapply(new_df, is.numeric)], function(col) {
  col[is.infinite(col)] <- NA
  col
})
# 删除含NA的行
new_df <- na.omit(new_df)

另外，原始剔除列的写法有误，Sex是列名，正确写法应为：

# 正确剔除Sex列的方式
new_df <- df[, !colnames(df) %in% c("Sex")]

2. 更高效的批量回归方法 #

针对7300列的大数据集，用purrr+broom组合比lapply更简洁，且方便结果整理：

library(purrr)
library(broom)

# 批量执行回归并提取所需统计量
results <- map(new_df[, -which(colnames(new_df) == "Cohort")], function(y_col) {
  fit <- lm(y_col ~ Cohort, data = new_df)
  # 一次性提取R平方和p值
  glance(fit) %>% 
    select(r.squared, p.value)
})

# 合并为结构化数据框
results_df <- bind_rows(results, .id = "Variable")

broom包的glance()函数可以直接提取关键统计量，避免手动从summary()或anova()中提取，大幅提升效率。

3. 无需子集化添加特征的方法 #

不用提前剔除无关列，直接动态构建回归公式即可加入特征。比如要加入Sex：

# 筛选出所有需要作为因变量的数值列
numeric_cols <- setdiff(colnames(df), c("Cohort", "Sex"))

results_with_sex <- map(numeric_cols, function(col) {
  # 动态生成公式：col ~ Cohort + Sex
  formula <- as.formula(paste(col, "~ Cohort + Sex"))
  fit <- lm(formula, data = df)
  glance(fit) %>% 
    select(r.squared, p.value) %>% 
    mutate(Variable = col)
})

# 合并结果
results_with_sex_df <- bind_rows(results_with_sex)

这种方法无需修改原始数据框，直接通过公式指定预测变量即可。

4. 额外优化建议 #

• 针对7300列的数据集，批量回归会消耗较多内存，可分批次处理，或用furrr包实现并行计算加速。
• 若Cohort是二分类变量，线性回归结果与t检验等价，用t.test()批量执行速度更快：

# 批量t检验示例（二分类Cohort）
t_results <- map(numeric_cols, function(col) {
  test <- t.test(df[[col]] ~ df$Cohort)
  tibble(Variable = col, p.value = test$p.value, mean_diff = test$estimate[1] - test$estimate[2])
})
t_results_df <- bind_rows(t_results)

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