R语言大矩阵跨列迭代成对比较线性模型检验实现

核心优化逻辑 #

• 简单线性回归的自变量斜率p值与两变量的皮尔逊相关系数检验p值完全等价，避免调用lm、summary这类包含大量冗余计算的函数，直接通过相关系数推导p值，速度可提升数十倍
• 用向量化运算批量处理每一轮基准列对应的所有后续列计算，消除逐次拟合的循环开销
• 逐轮写出结果无需累积全量数据，内存占用极低

实现代码 #

library(data.table)

# 此处可替换为你自己的数据集
set.seed(123)
n_sample <- 200
n_feature <- 1500
df <- as.data.frame(matrix(runif(n_sample * n_feature, min=0, max=100), nrow = n_sample))
colnames(df) <- paste0("data_", 1:n_feature)

# 预计算固定参数
n <- nrow(df)
df_resid <- n - 2 # 等价于简单线性回归残差自由度
col_names <- colnames(df)
n_col <- length(col_names)

# 循环执行成对比较
for (i in 1:(n_col - 1)) {
  # 提取当前基准列
  x <- df[[i]]
  x_cent <- x - mean(x)
  x_sd <- sd(x)
  
  # 提取所有排在当前列之后的待比较列，转为矩阵批量计算
  y_cols <- (i + 1):n_col
  y_mat <- as.matrix(df[, y_cols, drop = FALSE])
  y_cent <- t(t(y_mat) - colMeans(y_mat))
  y_sd <- apply(y_mat, 2, sd)
  
  # 批量计算皮尔逊相关系数
  r <- colSums(x_cent * y_cent) / ((n - 1) * x_sd * y_sd)
  # 计算t值与双侧p值，和lm输出的p值完全一致
  t_val <- r * sqrt(df_resid / (1 - r^2))
  p_val <- 2 * pt(-abs(t_val), df = df_resid)
  
  # 构造结果表并写入文件，输出格式与你的要求完全匹配
  res <- data.table(
    id1 = col_names[i],
    id2 = col_names[y_cols],
    pval = p_val
  )
  fwrite(res, file = paste0("data", i, ".tsv"), sep = "\t")
}

结果验证 #

你可以用5列以内的小样本数据对比上述代码输出的p值和你原代码中lm提取的p值，二者结果完全一致。

性能说明 #

该实现完全消除了lm函数的冗余计算开销，1500列合计110万次成对比较仅需数秒即可完成，远快于combn+lm的实现方案。

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