R语言用lm做回归时，如何提取底层的d和S矩阵？

示例数据与模型拟合 #

首先生成模拟数据并拟合线性模型：

# 生成100个正态分布的模拟变量
Y <- rnorm(100)
x <- rnorm(100)
z <- rnorm(100)

# 组合成数据框
data <- data.frame(Y, x, z)

# 拟合包含x、z的线性回归模型
fit <- lm(Y ~ x + z, data)

正规方程组的矩阵定义 #

线性回归的正规方程组矩阵形式为 (S\beta = d)：

• (S) 矩阵：设计矩阵 (X) 的叉积，即 (S = X^TX)（(X) 包含模型的截距项和所有自变量）
• (d) 向量：设计矩阵 (X) 与响应变量 (Y) 的叉积，即 (d = X^TY)

提取S矩阵（(X^TX)） #

有两种常用方法：

• 方法1：从lm对象的QR分解结果中推导
  lm模型默认用QR分解求解，利用QR分解的R矩阵可直接计算 (X^TX)：

  # 提取QR分解的R矩阵，计算R'R得到S矩阵
  S <- t(qr.R(fit$qr)) %*% qr.R(fit$qr)
  

• 方法2：构造设计矩阵后直接计算
  先获取模型的设计矩阵，再计算叉积：

  # 生成模型对应的设计矩阵X
  X <- model.matrix(fit)
  # 计算X的叉积得到S矩阵
  S <- crossprod(X)  # 等价于 t(X) %*% X，crossprod在R中更高效
  

提取d矩阵（(X^TY)） #

同样有两种常用方法：

• 方法1：从lm对象的QR分解结果中推导
  利用QR分解的Q矩阵与响应变量计算：

  # 提取QR分解的Q矩阵，计算Q'Y得到d矩阵
  d <- t(qr.Q(fit$qr)) %*% fit$y
  

• 方法2：构造设计矩阵后直接计算
  用设计矩阵与响应变量计算叉积：

  # 计算X与Y的叉积得到d矩阵
  d <- crossprod(X, fit$y)  # 等价于 t(X) %*% fit$y
  

内容的提问来源于stack exchange，提问作者J. Doe