基于最小二乘法的线性回归是你在学习数据科学和机器学习时最先遇到的模型之一，它不仅简单易懂，还能在很多问题中发挥作用，并且已经集成在了很多种编程语言之中。大部分用户对R语言中的lm()函数肯定不陌生，它让你能简易而快速地拟合一个线性回归模型。然而，这个函数并不现实参数估计和很多检验统计量的计算过程，所以本文就打算手把手地计算实现线性回归模型的参数估计。
本文中，我只会使用矩阵、向量和矩阵操作符来获得参数估计（基于最小二乘法的线性回归模型的参数估计有解析表达式）。在读完本文之后，你会发现自己动手也很简单，并且你可以把这个过程套用在任何数据集之上（尽管使用lm()函数会更快更稳健）。
我将使用MASS包中大名鼎鼎的Boston数据集作为例子，它包含了14个关于经济，地理和人口方面的变量，一共有506个观测值。每个观测代表1990年代波士顿某个地区的特征。我们的被解释变量将会是每个地区的房价的中位数，对应列名是'medv'，让我们先看看数据集长什么样：
library(MASS)
str(Boston)
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
你可以运行help(Boston)来看看每个变量的定义。
线性回归的表达式一般如下：
其中，y是响应变量（被解释变量），X是设计矩阵（特征变量），是需要估计的参数向量，
是残差项，它反映了哪些没被纳入模型的因素对y的影响。
万事俱备，第一步，我们要声明响应变量
y <- Boston$medv
然后，我们需要确定输入特征的设计矩阵，也就是公式中的X。我们可以使用as.matrix()函数完成这一操作，吧‘medv’以外的所有变量纳入。我们还需要在X中添加一列元素全为1的列来估计截距项。
# 设计矩阵
X <- as.matrix(Boston[-ncol(Boston)])
# 生成用来估计截距项的列
int <- rep(1, length(y))
# 添加用来估计截距项的列
X <- cbind(int, X)
现在被解释变量和设计矩阵都有了，我们可以通过如下的计算式来得到系数的估计值：
这个表达式的推导超出了本文的范畴，如果你有兴趣，可以参看Wikipedia 或者任何一本相关的教科书。
我们可以使用刚刚赋值好的X和y，结合%*%操作符来实现矩阵乘法。t()函数能将矩阵转置，solve()函数能够得到任意可逆矩阵的逆矩阵。
# 实现解析表达式
betas <- solve(t(X) %*% X) %*% t(X) %*% y
# 结果保留两位小数，更美观
betas <- round(betas, 2)
现在我们有系数向量了，它代表了y和解释变量之间的线性关系。比如，房间数每增长1单位，房价就增长3810刀。我们可以遵循这个方式来解释估计结果，但需要注意每个变量的量纲。注意，变量‘chas’是个只取0，1两个值的虚拟变量。1代表该区毗邻查理斯河，0代表相反情况。而对应的系数则表明，河边的房子的房价中位数比内陆的搞2690刀（在其他变量一定时）。
译者注：作者没有考虑变量的显著性就直接解释结果了，这是不对的。如果变量的系数并不显著，我们可以在统计意义上认为这个系数是0，也就是该变量和y之间没有明显的相关性。
print(betas)
## [,1]
## int 36.46
## crim -0.11
## zn 0.05
## indus 0.02
## chas 2.69
## nox -17.77
## rm 3.81
## age 0.00
## dis -1.48
## rad 0.31
## tax -0.01
## ptratio -0.95
## black 0.01
## lstat -0.52
现在，让我们来看看我们得到的结果是否正确，调用lm()建模并比较结果，我知道你们很想知道这个。
# 线性回归模型
lm.mod <- lm(medv ~ ., data=Boston)
# 保留两位小数
lm.betas <- round(lm.mod$coefficients, 2)
# 结果放到一个数据框内
results <- data.frame(our.results=betas, lm.results=lm.betas)
print(results)
## our.results lm.results
## int 36.46 36.46
## crim -0.11 -0.11
## zn 0.05 0.05
## indus 0.02 0.02
## chas 2.69 2.69
## nox -17.77 -17.77
## rm 3.81 3.81
## age 0.00 0.00
## dis -1.48 -1.48
## rad 0.31 0.31
## tax -0.01 -0.01
## ptratio -0.95 -0.95
## black 0.01 0.01
## lstat -0.52 -0.52
系数估计的结果一模一样，现在你知道怎么手把手得到系数估计值了。但是使用lm()函数无疑更快，并且它还提供了诸如拟合优度，t值，p值等检验统计量，还是用QR分解等技术使得计算更加稳健。可无论如何，我希望展示最为原始的解法。
线性回归的系数估计还可以通过其他方法得到，比如在特征数量很大的时候使用梯度下降法就是不错的选择。这个时候再利用解析表达式结算可能就很费时间。而随机梯度下降则在特征数和样本数都很大的时候适用，我会在后续文章中介绍它们。
注：原文刊载于Datascience+网站