一般线性模型中p值与线性回归斜率的相关性问题咨询

我将参考fasta中600条序列的计数表（read counts）输入到一个使用R语言lm()函数实现一般线性模型的工具包中，得到每条序列的趋势线斜率（b.value）和p值。绘制火山图后，发现回归斜率与p值存在意外的相关性，正在排查工具包代码，想了解该情况是否可能由输入数据而非模型代码导致？

附火山图：图中显示回归斜率（b.value）与p值存在明显的相关性趋势。

用到的线性回归代码： #

do.lm <- function(dat) {
  count <- dat$count
  col_data <- dat$col_data
  tidy.count <- count %>% 
    cbind(name = rownames(count)) %>%
    gather(sample, count, -name) %>% 
      left_join(col_data %>% rownames_to_column("sample"), by = "sample")    
  # TIP: values in a column must be atomic, can't have a vector
  res <- tidy.count %>% group_by(name) %>% 
    summarise(count = list(count), 
              condition = list(condition)) %>%
    group_by(name) %>%
    mutate( lm = list(summary(lm(unlist(condition)~unlist(count)))),
            baseMean = mean(unlist(count))) %>%
    mutate( p.value = tryCatch({lm[[1]]$coefficients[2,4]}, error = function(e) NA),
            b.value = tryCatch({lm[[1]]$coefficients[2,3]}, error = function(e) NA)) %>%
    select(name, baseMean, b.value, p.value)
  res$padj <- p.adjust(res$p.value, method="fdr")
  res
}

一、核心问题：代码逻辑错误直接导致相关性 #

你看到的斜率与p值的相关性，首先是代码错误导致的，而非输入数据，关键错误点有两个：

1. 回归方向完全颠倒：当前公式是unlist(condition)~unlist(count)，把condition作为因变量、count作为自变量，这和常规的差异分析逻辑（用条件解释计数变化）完全相反，会导致统计量的解读彻底混乱。
2. 统计量提取错误：你提取的b.value是系数矩阵的第3列，这是t值而非斜率（斜率是系数矩阵的第1列）。而t值和p值本来就直接相关——t值的绝对值越大，对应的p值越小，这才是你看到两者关联的根本原因。

二、输入数据可能导致关联的场景（修正代码后再考虑） #

如果修正代码后仍存在类似相关性，输入数据层面的可能原因包括：

• 测序深度偏差：部分样本测序深度过高，count值异常大，回归时被视为极端值，既拉动斜率估计，又让t值变大、p值变小，形成关联。
• 条件变量编码问题：若condition是分组变量却用了不恰当的数值编码，或存在未控制的混杂变量（如批次效应），可能导致回归结果出现虚假关联。
• 低表达序列干扰：大量低表达（count接近0）的序列，回归残差波动大，p值不稳定，同时斜率估计易出现极端值，形成聚集性关联。

三、代码修正示例 #

修正回归方向和统计量提取逻辑，修正后的代码如下：

do.lm <- function(dat) {
  count <- dat$count
  col_data <- dat$col_data
  tidy.count <- count %>% 
    cbind(name = rownames(count)) %>%
    gather(sample, count, -name) %>% 
      left_join(col_data %>% rownames_to_column("sample"), by = "sample")    
  
  res <- tidy.count %>% group_by(name) %>% 
    summarise(count = list(count), 
              condition = list(condition)) %>%
    group_by(name) %>%
    mutate( lm = list(summary(lm(unlist(count) ~ unlist(condition)))), # 修正回归方向：count为因变量，condition为自变量
            baseMean = mean(unlist(count))) %>%
    mutate( p.value = tryCatch({lm[[1]]$coefficients[2,4]}, error = function(e) NA),
            b.value = tryCatch({lm[[1]]$coefficients[2,1]}, error = function(e) NA)) %>% # 修正斜率提取：取系数矩阵第1列（估计值）
    select(name, baseMean, b.value, p.value)
  res$padj <- p.adjust(res$p.value, method="fdr")
  res
}

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