R语言lm类对象调用vcov()函数的计算原理及手动实现方法

线性回归系数的方差-协方差矩阵的标准计算公式为：

$Var(\hat{\beta}) = \hat{\sigma}^2 \times (X^TX){-1}$

其中：

• $\hat{\sigma}^2$ 是回归模型的残差均方误差，等于残差平方和除以模型自由度（样本量 - 待估参数个数）
• X 是模型的设计矩阵，默认包含截距项对应的全1列
• $(X^TX){-1}$ 是设计矩阵叉积的逆矩阵

该公式同时覆盖对角项（系数方差，即标准误的平方）和非对角项（不同系数的协方差）的计算需求，可运行的验证代码如下：

library(data.table)
# 固定随机种子保证结果可复现
set.seed(123)
df <- data.table(y = runif(100, 0, 100),
                 x = runif(100, 0, 100),
                 z = runif(100, 0, 100))

# 先拟合模型，取官方vcov作为对比基准
fit <- lm(y ~ x + z, df)
official_vcov <- vcov(fit)
cat("官方vcov输出结果：\n")
print(official_vcov)

# ------------------- 手动计算逻辑 -------------------
# 1. 提取包含截距项的设计矩阵
X <- model.matrix(fit)
# 2. 计算残差均方误差
n <- nrow(df) # 样本量
k <- length(coef(fit)) # 待估参数个数
res <- residuals(fit)
sigma_hat_sq <- sum(res^2)/(n - k)
# 3. 计算设计矩阵叉积的逆
xtx_inv <- solve(t(X) %*% X)
# 4. 得到完整方差-协方差矩阵
manual_vcov <- sigma_hat_sq * xtx_inv
cat("\n手动计算vcov结果：\n")
print(manual_vcov)

# 数值误差范围内验证二者完全相等
cat("\n二者一致性验证结果：\n")
print(all.equal(official_vcov, manual_vcov))

代码运行后会输出一致的方差-协方差矩阵，你需要的非对角项会和官方vcov()返回的结果完全匹配。

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