多元线性回归的系数及其标准差估计

专注系列化、高质量的R语言教程

推文索引 | 联系小编 | 付费合集

线性回归是最基础的回归模型，但不知道有多少读者了解它的回归系数以及标准差是如何估计出来的。本篇就来介绍一下，目录如下：

1 符号说明
2 系数估计
3 系数标准差
4 相关函数和操作符
- 4.1 %*%
- 4.2 t函数
- 4.3 solve函数
- 4.4 diag函数
5 案例

1 符号说明

使用表示样本标识，表示样本的因变量取值，表示自变量表示（，其中为自变量个数），表示样本的一系列自变量取值，表示随机项。

线性回归的方程如下：

使用矩阵可以表示为如下形式：

其中，和都来自已有的样本数据。

为的满秩矩阵（为样本数，为自变量个数），行表示样本，列表示变量，也称设计矩阵：

是长度为的列向量：

为待估计的模型系数，是长度为的列向量：

为随机项，也是模型的残差，是长度为的列向量：

从方程组的角度看，和都属于未知数，共计个，而方程个数为，因此方程组是不可解的，它的自由度为未知数个数与方程个数之差，即。

2 系数估计

既然方程组是不可解的，我们可以使用优化的方法去估计出“最优”的系数组合。

众所周知，多元线性回归的优化目标是残差平方和最小。残差平方和为

复习一下，转置矩阵有如下运算性质：

因此，

从而，

问题转化为求取得最小值时的。可以看出，是一个二次型，它的最小值在偏导为0处取得。

使用矩阵直接求导有如下运算性质^[1]：

其中，、、表示列向量，表示方阵。

因此，

令，即

可得的估计值为

3 系数标准差

因为，

显然有。

的方差是下面矩阵的对角线元素：

线性回归假设所有样本的随机项都服从同一个均值为0的正态分布，即

因此，

并且不同样本之间的随机项相互独立。因此，

所以，

进而，

取上面矩阵的对角线元素即为系数估计值的方差：

标准差为：

在回归模型中，系数估计值的标准差一般称为标准误（standard error, se）。其中，的无偏估计为：

4 相关函数和操作符

上面推导过程中涉及到一些R语言中的函数和操作符。

4.1 %*%

*用于矩阵相乘表示同型矩阵对应位置的元素相乘，而%*%才表示矩阵真正的乘法。

A = matrix(1:12, nrow = 3)
B = matrix(1:12, nrow = 4) A %*% B
##      [,1] [,2] [,3]
## [1,]   70  158  246
## [2,]   80  184  288
## [3,]   90  210  330

4.2 t函数

t()函数表示矩阵的转置。

A
##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12t(A)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
## [4,]   10   11   12

在R语言中，向量是没有维度的：

a <- 1:3dim(a)
## NULL

但若取其转置，会默认其为列向量，转置结果为行向量（对应的数据结构实际上是矩阵）：

t(a)
##      [,1] [,2] [,3]
## [1,]    1    2    3

4.3 solve函数

solve()函数的本职功能是解形如a %*% x=b（其中a为方阵）的矩阵方程：

a = matrix(c(1,1, 1, -1), nrow = 2)
b = c(3,1)solve(a, b)
## [1] 2 1

但当b缺失时，会默认其为单位矩阵，进而方程组的解为，等价于求a的逆矩阵：

solve(a)
##      [,1] [,2]
## [1,]  0.5  0.5
## [2,]  0.5 -0.5

4.4 diag函数

diag()函数的功能主要有4个：

获取已知方阵的对角线元素：

A = matrix(1:9, nrow = 3)
A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9diag(A)
## [1] 1 5 9

将已知方阵的对角线元素修改为指定值：

diag(A) <- c(1,2,3)
A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    2    8
## [3,]    3    6    3

生成指定大小的单位矩阵：

diag(3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

生成指定元素和大小的对角矩阵：

diag(2, 3)
##      [,1] [,2] [,3]
## [1,]    2    0    0
## [2,]    0    2    0
## [3,]    0    0    2diag(c(1,2, 3), 3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3

第一个参数表示对角线元素，第二个参数表示矩阵行数。

5 案例

下面使用一个案例验证系数和标准误的估计结果。

使用R语言的lm()函数运行一个线性模型并输出结果：

model <- lm(mpg ~ wt + qsec, data = mtcars)
summary(model)## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  19.7462     5.2521   3.760 0.000765 ***
## wt           -5.0480     0.4840 -10.430 2.52e-11 ***
## qsec          0.9292     0.2650   3.506 0.001500 **

手动计算模型系数和标准误：

X <- mtcars[c("wt", "qsec")]
Y <- mtcars["mpg"]X <- as.matrix(X)
Y <- as.matrix(Y)n <- dim(X)[1]
m <- dim(X)[2]
X <- cbind(rep(1,n), X)  ## 系数估计值inv <- solve(t(X) %*% X)
inv %*% t(X) %*% Y
##            mpg
##      19.746223
## wt   -5.047982
## qsec  0.929198### 标准误估计值 sigma2 <- sum(residuals(model)^2)/(n-m-1)
sqrt(diag(sigma2 * inv))
##                  wt      qsec
## 5.2520617 0.4839974 0.2650173

比较可知，手动计算结果与lm()函数运行结果一致。

参考资料

[1]

矩阵求导公式的数学推导（矩阵求导——基础篇）: https://zhuanlan.zhihu.com/p/273729929