R语言逐行拟合线性回归计算每行对应p值的实现问题

需求合理性判断 #

直接为单行数据拟合线性回归模型的需求不符合统计逻辑：线性回归至少需要2个独立观测才能估计自变量系数、计算残差和对应的显著性p值，仅用1行数据拟合时自由度为0，无法输出有效的统计结果。

可行实现方案 #

根据组学分析的常见场景，这里提供两种符合实际需求的实现方式： #

方案1：按基因分组拟合回归，返回组层面p值 #

如果你是想评估每个基因的染色质开放区域倍数变化和表达量倍数变化的相关性，可以按ENSEMBL基因ID分组，每组内所有Peak观测拟合回归，将该基因的回归p值赋值到组内所有行：

# 加载依赖包
library(dplyr)
library(broom)

# 构造你的数据集
a <- structure(list(ENSEMBL = structure(c(1L, 2L, 3L, 3L, 3L, 4L), .Label = c("ENSG00000005187", 
"ENSG00000006740", "ENSG00000008277", "ENSG00000013810"), class = "factor"), 
    log2FoldChange_Expression = c(-2.2756549273843, -1.76655532051033, 
    -1.58489726654531, -1.58489726654531, -1.58489726654531, 
    -2.04282868170093), log2FoldChange_Region = c(-2.11261476936419, 
    -2.37119008459253, -1.59565539803813, -2.4954310786834, -2.11050911441613, 
    -1.81996408306615), Peak_Region = structure(c(5L, 6L, 4L, 
    2L, 3L, 1L), .Label = c("Peak147010", "Peak194531", "Peak194535", 
    "Peak194536", "Peak75759", "Peak81940"), class = "factor")), class = "data.frame", row.names = c(NA, 
-6L))

# 分组计算p值合并回原数据
result <- a %>%
  group_by(ENSEMBL) %>%
  # 仅当组内样本数≥2时计算p值，否则返回NA
  mutate(pvalue = ifelse(n() >= 2,
                         glance(lm(log2FoldChange_Expression ~ log2FoldChange_Region))$p.value,
                         NA_real_)) %>%
  ungroup()

该方案下，只有观测数≥2的基因会返回有效p值，你的示例数据中仅ENSG00000008277有3个观测，因此仅该基因对应的3行会有p值，其余基因返回NA。 #

方案2：拟合全局回归，返回每行的离群显著性p值 #

如果你确实需要每行对应一个p值，用来评估单个Peak-基因对的变化趋势是否偏离整体规律，可以先对全量数据拟合全局回归，再计算每行的学生化残差对应的p值：

# 拟合全局回归模型
full_lm <- lm(log2FoldChange_Expression ~ log2FoldChange_Region, data = a)
# 计算每行的学生化残差p值
a$pvalue <- 2 * pt(-abs(rstudent(full_lm)), df = full_lm$df.residual)

该方案返回的p值越小，说明对应行的观测越偏离整体回归的拟合趋势，可用于筛选异常的Peak-基因关联对。

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