多元线性回归的系数及其标准差估计
阿里云国内75折 回扣 微信号:monov8 |
阿里云国际,腾讯云国际,低至75折。AWS 93折 免费开户实名账号 代冲值 优惠多多 微信号:monov8 飞机:@monov6 |
专注系列化、高质量的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 12
t(A)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9
## [4,] 10 11 12
在R语言中向量是没有维度的
a <- 1:3
dim(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 9
diag(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 2
diag(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