为何运行相同的R语言lm回归代码会得到不同结果？

问题描述 #

执行无随机性的二次项线性回归代码，每次运行得到不同的系数结果，该问题仅在本地Ubuntu 20.04环境中出现，其他机器运行结果一致。相关代码及环境信息如下：

运行代码及输出 #

d = read.csv("d.csv")
str(d)
## 'data.frame':    5000 obs. of  2 variables:
##  $ elevation_sd: num  -0.6123 1.3988 -0.111 -0.0763 0.4129 ...
##  $ pred_rich   : num  6.74 10.35 10.81 15.24 37.58 ...
coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
##       (Intercept)      elevation_sd I(elevation_sd^2) 
##       31.40967132      -14.51366622       -0.04686489
coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
##       (Intercept)      elevation_sd I(elevation_sd^2) 
##         24.791451        -22.377711          2.530645
coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
##       (Intercept)      elevation_sd I(elevation_sd^2) 
##         29.009348        -24.539297          2.184579

系统及R环境信息 #

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libmkl_rt.so
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     

排查与解决建议 #

• 检查MKL库的非确定性优化
  你的环境使用Intel MKL作为BLAS/LAPACK后端，MKL多线程模式下可能因浮点运算的并行顺序差异导致结果波动。先强制单线程运行MKL测试：

  Sys.setenv(MKL_NUM_THREADS=1)
  d <- read.csv("d.csv")
  coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
  

  若结果稳定，说明是多线程优化导致；若仍不稳定，可切换到系统默认的OpenBLAS库：

  # 终端安装OpenBLAS（未安装时执行）
  sudo apt install libopenblas-dev
  

  然后在R中设置环境变量并重启R：

  Sys.setenv(R_BLAS="/usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3")
  Sys.setenv(R_LAPACK="/usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3")
  

• 排查数据隐性修改问题
  排除CSV读取过程中的潜在歧义，将数据保存为二进制RDS格式，每次运行前重新读取：

  # 仅需执行一次保存操作
  saveRDS(d, "d.rds")
  # 每次回归前读取RDS文件
  d <- readRDS("d.rds")
  coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
  

• 解决共线性导致的系数不稳定
  二次项模型中，自变量与平方项易存在严重共线性，放大浮点计算误差。先对elevation_sd做中心化处理：

  d$elev_center <- d$elevation_sd - mean(d$elevation_sd)
  coef(lm(pred_rich ~ elev_center + I(elev_center^2), data = d))
  

  也可通过方差膨胀因子（VIF）验证共线性程度：

  library(car)
  model <- lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d)
  vif(model) # VIF>10说明共线性严重
  

• 排除环境干扰
  启动干净的R会话（不加载任何额外包），仅运行核心代码，排除其他包的潜在影响：

  # 关闭所有R会话后重新启动，仅执行以下代码
  d <- read.csv("d.csv")
  coef(lm(pred_rich ~ elevation_sd + I(elevation_sd^2), data = d))
  

• 检查系统浮点环境
  在终端执行以下命令，检查是否存在异常的浮点精度设置：

  printenv | grep -i FP
  

  若有异常设置，临时取消后再启动R测试。

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