线性回归是最基础的统计模型，也应该是应用最广泛的应用模型。线性回归的因变量必须为连续变量；自变量可以为连续变量，也可以为分类变量（转化为哑变量/虚拟变量）

这里是我通过看资料对线性回归的总结。

参考资料

An Introduction to Statistical Learning，下文简称 ISL

线性回归的假设

1 自变量和因变量线性相关

线性相关的意思是，无论自变量多大，自变量每增加一个单位，因变量都会产生固定的变化。

2 多个自变量之间具有可加性

无论 X1 多大，X2 每增加一个单位，因变量都会产生固定的变化。如果满足这个条件，说明 X1和X2 之间具有可加性。否则说明 X1 和 X2 之间存在互作效应。用图表示如下（左图两条线平行说明不存在互作，右图不平行说明存在互作）：

3 多个自变量之间彼此独立

多个自变量之间不能存在相关关系。相关关系分为两种，一种是简单的两两之间的相关，这个通过相关系数矩阵就能看出来是否存在两个变量之间的明显的相关（看绝对值）；第二种是多重线性相关（Multicollinearity），例如 X1 = a + bX2 + cX3，一个自变量可以用其他的多个自变量来解释，就说明存在多重线性相关。

4 观测值的误差项之间互相独立

误差项即实际值 - 预测值。如果观测值之间存在我们没有添加到模型中的其他影响因素，那么误差项之间就可能存在相关关系，比如有些样本搜集自同一个家庭的成员，而我们没有把家庭这个因素加入到模型中。

5 误差项（ε）呈正态分布，期望为0，方差为定值

这里其实分为两个假设。第一个假设：误差项服从均值为0的正态分布。第二个假设：误差项的方差为定值（不变）。这两个假设是为了保证回归模型在小样本下能够顺利进行假设检验。正态分布假设仅在小样本的情况下需要，大样本的情况下则不需要，因为有中心极限定理做正态性的支撑。而方差齐性则保证最小二乘法估计出来的统计量具有最小的方差。如果违反了这个假设，置信区间会变宽，这称之为异方差性（heteroscedasticity）。当异方差性出现的时候，如果仍采用最小二乘法估计参数，会导致参数的t检验值被高估，可能造成本来不显著的某些参数变为显著，使假设检验失去意义。

https://www.cnblogs.com/HuZihu/p/10142737.html

6 自变量和误差项之间互相独立

模型中一个或多个自变量与随机误差项存在相关关系，这称之为内生性（endogeneity）。内生性通常由于遗漏变量而导致的，因此是一个普遍存在的问题。内生性会导致模型参数估计不准确。

https://www.cnblogs.com/HuZihu/p/10142737.html

简单线性回归

β0, β1 参数估计 - 最小二乘法

最小二乘的思路是，查找使得误差平方和或均方误差最小的参数。

我们首先知道，误差 = 实际值 - 预测值，我们的目的实际是想要所有观察值的总体的误差最小，因为误差有正有负，如果采用误差项直接相加会互相抵消，因为我们采用误差的平方进行相加，得到误差平方和 (RSS)。

或者写成

除以总共的样本数 n，便得到均方误差，使得 RSS/MSE 最小的参数便是最小二乘的最优解。写成公式如下（缺推导过程）：

画成等高图如下（作图横纵坐标为β0和β1，图中一圈表示相同的RSS，数字即为 RSS 大小），我们可以看到这是向下凹的图形，只有一个点RSS 最小。

为什么不用残差的绝对值？

我上统计课的时候就有这个问题，就一直有这个疑问，如果因为残差有方向，不能直接求和，那为什么不直接用绝对值呢？这不是更直接的更贴近残差本质的思路吗？

使用残差绝对值作为代价函数，在统计的历史上是存在的，称为最小一乘法。查看了很多网页的回答，感觉这个回答比较全 为什么是“最小二乘”而不是“最小一乘”？

最小一乘法的解不一定唯一

举个简单的例子，根据下图中的四个点，可以构建出无数条回归直线，只要这些直线同时与 AB, CD 相交即可。

最小二乘求解方便

后来上了吴恩达老师的机器学习课，知道了代价函数的概念，寻找代价函数最小值的最简单的方法就是求导。最小二乘法的代价函数就是 RSS/MSE，由于这里采用平方和（二次函数），求导会直接变成一次函数。而如果使用残差的绝对值的和，根本无法进行简单的求导，就无法简单地得到最优解的公式（而且绝对值在 x = 0 处不可导）。

最小二乘估计结果与正态分布假设下的最大似然估计结果相同

可以证明，当误差服从正态分布时，最大似然法估计结果和最小二乘相同（注意最小二乘与最大似然是截然不同的两种思想，两种方法。最小二乘没有对观测值的分布做假设）。下方证明参考最小二乘法为什么使用平方误差 。

当你用线性回归拟合数据时，其实模型包括了两部分，回归项和残差项：

之后我们用最大似然方法来估计参数，假设样本数目为n，似然函数表达式如下 (给定一组参数下得到这些观测值的似然值，等于得到这些残差的似然值)：

我们需要找到使得似然函数最大的一组参数。

假设残差独立同分布于正态分布 N(0, σ²) ，单个残差的概率密度函数如下：

多个残差服从多元正态分布，联合概率密度函数如下：

因此这个问题就转换成了求使得**联合概率密度函数f(ε)最大 **的一组参数，这组参数就是最大似然法的估计参数。仔细查看 **联合概率密度函数f(ε) **的式子，不难发现该式最大时，即误差平方和最小。

因此当残差符合正态分布时，最小二乘法的参数估计值与最大似然法一致，最小二乘法的参数估计值就是最可能的参数估计值。

为什么不用三次方，平方和？

为什么最小二乘法要用最小误差平方和来拟合（为什么不用立方和，或者开方和）? 首先你要保证误差项永远是正数，因此三次方和其他奇数次方不成立。对于四次项及更高的偶数次项，效果还没有最小二乘好，因为会放大极端值的影响。开方更不可能了，因为残差可能是负数，无法求开方。

参数估计准确性

我们使用的模型为

这里的残差项包含了所有我们模型中没有考虑到的地方，比如Y 和 X 的关系不是线性的，还有其他的影响因素，存在测量误差等等。

为了估计参数的估计值与真实值之间的误差，我们需要计算两个参数估计值的标准误，使用下列公式（缺证明）

这里 σ2 是残差方差，是未知的。残差标准差 σ 可以用 RSE (residual standard error) 来估计（缺证明）。

知道了标准误，我们可以计算β0、β1的置信区间，例如 β1 的 95%的置信区间为：

如果 β1 的 95%的置信区间 包含 0，那说明不能拒绝 β1 = 0 的原假设，我们认为不存在线性相关。

我们也可以通过 t 检验得到 P值。构建 t 统计量如下（自由度 n-2）：

同理可以检验β0 估计值的P值（零假设：β0 = 0），但实际上截距的显著性没啥大的意义，主要还是看 β1。如果 β1 的 P值小于 0.05，我们认为 β1 ≠ 0。

模型准确性

上面是得到估计的两个参数的准确性，这里我们想知道，我们构建的简单线性模型与数据的契合程度（which the model fits the data）。

线性模型一般用两个指标， residual standard error(RSE) 和 R² 。

RSE

RSE 是残差的标准误的估计值，它表示预测值与真实值的平均偏差量（这句话需要理解标准误的概念：标准误是样本统计量的标准差，均值为真值，结合标准差的公式，可知标准误就是预测值与真实值的平均偏差量（绝对值），所以一般提供参数估计值的时候，都会写成 参数估计值±标准误的形式）。

RSE 是最直接说明模型准确性的统计量，它的缺点是有单位，是一个绝对值，仅仅是RSE你不知道到底是大还是小，所以你还得和因变量的均值啥的进行比较才知道相对大小。

R²

R² 指的是因变量方差能够被线性模型解释的比例，始终在0和1之间，是一个相对值，越大越好。公式如下：

其中TSS 计算如下，表示因变量的总的方差。

R² 比 RSE 更容易解释，因此更常用，但是还是面临一个问题，R² 多大才算结果很好呢？这可能得根据实际情况来定，在我们一般的实验中，实际上往往还有很多未测量的影响因素，所以 RSS 可能比较大，这时可能 R² 哪怕比较小可能也说明有用。

多重线性回归

简单线性回归只有一个自变量，但在实际情况中，我们往往有多个自变量（即生活中其实几乎不存在简单线性回归的应用场景）。例如，如果我们有三个影响因素，我们如何同时分析这三个影响因素呢？

一种思路是我们分别做上次简单的线性回归，但是这种思路无法令人满意。首先，如果给定一组三个自变量的值，我们无法给出因变量的预测值；其次，每一次简单回归都忽视了其他两个影响因素的作用，这可能会带来严重的误判。更好的解决思路是同时将多个自变量均放到线性回归的模型中。

ISL 书中提到了一个例子，销量和三种媒体预算的回归分析。在简单线性回归下，报纸的作用是显著的，如下图。

但是，同时考虑三种媒体后，报纸的P值是不显著的。这说明简单线性回归和多重线性回归的结果是不一样的。这种差距的出现的原因是，简单回归得到的斜率是不考虑其他因素（TV 和 radio）下得到的，而多重回归得到的斜率是保持 TV 和 radio 固定不变的基础上计算得到的。

这种情况说明了什么呢，说明这三个自变量之间并不是互相独立的，通过下面的相关矩阵我们可以看到，radio 和 newspaper 的相关系数为0.35，表示当投资人在 radio advertising 花更多的钱的同时，存在一种趋势在 newspaper advertising 上也花更多的钱。因此，我们对 newspaper 进行简单回归时得到的显著性，实际不是 newspaper 的作用，而是 radio 因素的作用。

这样的例子实际上比比皆是。ISL 提了一个离谱的例子来说明这一点。如果我们计算沙滩冰淇淋销量对鲨鱼伤人事件的数量求回归，我们应该能得到一个正向的关系，类似于上面的 newpaper 的例子。但是这实际是很荒谬的，这两个事件不存在因果关系，你无法通过禁止销售冰淇淋来减少鲨鱼伤人的事件。实际上是高温天气促使了人们去海边，从而导致更多的冰淇淋销量和鲨鱼伤人事件。如果你同时考虑冰淇淋销量和天气，对鲨鱼伤人事件做多重回归，你就会发现冰淇淋销量不再是一个显著的影响因素。

决定重要的影响因素

当我们有多个影响因素时，很有可能所有的影响因素都和因变量是有联系的，但是因变量可能仅仅是由其中的一小部分自变量决定的（相关不等于因果，比如上面的冰淇淋销量的例子）。我们把如何挑选自变量的过程称为 variable selection。

理想情况下，我们可以通过穷举所有的模型，每种模型均包含一部分的自变量。通过对比所有模型的表现，我们就可以得到最佳的自变量的集合。例如假如我们手上总共有2个可选的自变量X1 X2，那么总共就有4种可能的模型，（1）不包含任何自变量 ；（2）包含X1 ；（3）包含X2；（4）包含X1和X2。但是我们用什么指标来判断不同模型的优劣呢？这里有很多的判定模型优劣的统计量，比如AIC， BIC , 校正后R² 等。

为什么使用校正后的R²

我们看R² 的计算公式，可以看到 TSS 是固定值。如果你增加了变量后（比如原来只包含X1，到包含X1和X2），残差平方和只会减少（这就是overfitting的原理，哪怕仅仅是增加了无关变量，也会造成与本次数据的拟合程度“变好”）。这导致增加了变量后，R²只会进一步增加。

校正后的R²的公式如下，考虑了自变量个数的影响。可以看到，校正后的R² 随着自变量个数的增加而降低，并且始终低于未校正的R²。

ISL 还提到了了使用 RSE 来校正自变量个数的影响。

三种经典的挑选变量的思路

• Forward selection: 先从null model (没有任何变量的模型)出发，对所有的候选因素（假设 p 个）各做一次简单回归，然后挑选结果中 RSS 最低的变量加入到模型中。然后我们将剩下的 p-1 个候选因素，挨个加入到目前的模型中，再挑选结果中 RSS 最低的变量加入到模型中。一直循环下去，知道达到了某个终止条件。
• Backward selection: 首先将所有的变量均加入到模型中，然后移除其中 P值最大的变量，我们就得到了一个 p-1个变量的模型。循环往复直到达到终止条件（比如所有剩下的变量的P值均低于某个值）。
• Mixed selection: 结合了上述的两个方法，我们先从空的模型出发，就像第一种方法的示例，不停地往里加变量，直到最佳。在加变量的过程中，我们观察有没有变量的P值变大并且超出了我们设定的阈值，如果存在这种现象，我们就从模型中移除这个变量（backward）。我们一直重复这种添加变量和删除变量的操作，直到最终实现所有剩下的变量的P值均很低，所有剔除的变量只要加入到模型中就会有一个很大的P值。

如果变量的个数超过了记录数，那么我们无法用 Backward selection (因为将所有变量均加入到模型中，会发生维数灾难)，此时 Forward selection 是可用的。但是Forward selection 是一种贪婪算法，早期加入的变量可能是多余的，而 Mixed selection 算法修正了这一点。

变量之间是否存在互作

还是采用上面三种媒介的例子，我们已经确定了 newpaper 没有作用， TV 和 radio 对 sales 有作用，那是不是模型就是下面的形式呢？

1

Y = β0 + β1*TV + β2*radio + e

我们可以通过这个模型画出一个三维图，图中的平面为预测值。我们可以看到当主要的钱都投向一种媒体时，销量会低于预测值；而如果我们均衡投资这两种媒体，销量会高于预测值（虽然我没有很好的看懂这个三维图）。这说明TV 和 radio 这两个因素之间存在互作。

预测的不确定性

一旦我们得到了最终的多重线性回归模型，我们就可以基于一套自变量的值来预测反应变量的值。但是，我们所做的预测至少存在三种不确定性。

1. 我们得到的仅仅是 β0, β1, β2……等参数的估计值。这些参数的估计值与真实值存在差距，这种差距是可消除误差(reducible error) 的一种。我们统计计算一个 confidence interval 来决定估计参数得到的预测值与真实参数得到的预测值的差距。

2. 第二种不确定性是由于线性模型仅仅是对实际数据的一个近似，这种带来的错误我们称之为model bias 。

3. 即使我们知道真实的模型 f(x)（我们知道所有参数的真实值），我们也不能完全精确地预测因变量，因为还存在随机误差(random error)，这种属于不可消除误差(irreducible error)。真实模型的预测值与真实值之间的差距，用prediction intervals 来评价。prediction intervals 通常比 confidence intervals 更宽， 因为 prediction intervals 不仅包括了可消除误差，也包含了不可预测误差（ prediction intervals 包含了 confidence interval ）。下面这段话解释地更明确一点：

   

线性回归潜在的问题

1 自变量与因变量之间为非线性关系

线性回归假设因变量与自变量之间存在直接的线性关系。如果因变量与自变量之间不是线性关系，那么线性回归的准确性会显著降低。

我们一般通过残差图来确定非线性关系。如果是简单线性回归，横坐标可以是自变量X，纵坐标是残差e；但是如果是多重线性回归，由于存在多个自变量，横坐标改为因变量的预测值。

下边左图为简单线性回归的残差图，可以看到残差不是均匀地分布在均值为0这条线地上下两侧，说明存在非线性关系。右图通过多项式回归，残差图的结果变好了很多，虽然看上去还是有一点非随机分布。

2 残差项的相关

线性回归的一个重要的假设是残差项，e₁,e₂,e₃…… e_n 之间是不相关的。这句话的含义是，如果残差项是不相关的，那么某一个残差项的值不受其他残差项的影响。ISL 原文如下：

For instance, if the errors are uncorrelated, then the fact that ε_i is positive provides little or no information about the sign of ε_i+1.

预测参数或预测值的标准误都是基于残差不相关的假设计算得到的。如果实际上残差项之间存在相关，那么标准误倾向于被低估（缺解释），因此置信区间也会比正常情况下更窄，因此 p值会更低，容易出现假阳性。

下面这段话举了个例子，但是我没理解。这里假如我们不小心把数据重复了一遍（本来是n条数据，现在是2n条数据），参数估计的置信区间会缩减 √2 。

**什么时候会出现残差项的相关？**这种情况频繁出现在 time series 数据，也有在不同的时间点去观测得到的数据。在很多情况下，相邻时间观测到的记录之间会有正相关的残差项。为了确定这种情况，我们可以画一个横坐标是时间的残差图，如果残差项之间是不相关的，你就不会观测到明显的模式（残差应该是忽上忽下）；相反，如果残差项之间是正相反的，你可以会发现相邻时间的残差会倾向于有相似的值。

下图就体现了这一点，最上面是没有残差相关的图，越往下残差相关越高，我们可以看到很清晰的模式 - 相邻的点之间倾向于有相似的值。

除了 time series 数据以外，其他数据也可能出现残差相关。举个例子，假设我们有一个研究通过体重去预测一个人的身高。如果有一些研究个体来自于相同的家庭，或吃相同的饮食，又或者面临某些相同的环境因素，那么残差不相关的假设就可能被打破。

一般来说，残差不相关的假设对于线性回归和其他统计方法是非常重要的。良好的实验设计对于减少残差不相关的风险至关重要。

为什么会出现残差相关？

因为残差是个筐，啥都往里装。我们的模型中所有没有考虑到的影响因素都会包含在这个残差里。实际上我们的模型不可能包含所有的影响因素，我们也不可能收集到所有影响因素的数据，所以我感觉残差相关是不可能完全避免的，或多或少有一点，只要不是很严重。

所以实验设计很重要，你要事先把重要的影响因素考虑到，而且把这些数据收集到。如果没有收集到相应的数据，貌似都无法发现残差相关，比如上面的 time series 数据，如果你没有搜集每个数据的时间信息，你都没法画上面的残差图。你不知道就等于没有了。

3 残差方差不齐性

线性回归的另一个重要的假设是所有残差项均有一个固定的方差。我们计算标准误、置信区间等都依赖于这个假设。

但不幸地是，残差项的方差不是固定的，才是通常的情况。比如，残差方差可能会随着因变量值的提高而提高。我们可以通过残差图来发现这一点，比如下图的左图。当面临这种情况时，一个可能的解决思路是将反应变量Y 转换为 log(Y) 或 √Y ，这两个函数对于较大的值都会有一个很好的收缩效应，从而导致方差不齐性的改善。下图的右图就展示了转换后的残差图，我们可以看到方差不齐性的结果大大改善了（残差的宽度都差不多了），但是存在轻微的非线性关系。

4 Outliers

Outliers 的定义是这个点的真实值与预测值偏差很大。Outliers 出现的原因有很多，比如错误记录。

下图中的红点（观测值为20）就是一个典型的Outlier。左边图中的红线是包含了Outlier的回归线，蓝色虚线是剔除了Outlier的回归线，我们可以看到移除或不移除这个Outlier对参数估计的影响很小。

但是Outliers 还会造成另外的问题。比如，在这个例子中，如果包括Outlier 计算得到的 RSE 为 1.09 ，如果剔除Outlier 得到的 RSE 为 0.77 。因为 RSE 在计算置信区间和 p值时均会用到，这样因为一个点导致RSE剧增会影响到回归方差拟合效果的解释。相似地，加入这个 Outlier， R² 从 0.892 降低到 0.805。

我们可以通过残差图来找到Outlier。在本例中我们轻松地找到这个Outlier，但是在实际应用中，我们很难决定到底偏离多少才可以视为一个 Outlier。为了解决这个问题，我们可以画图studentized residuals 图，纵坐标时残差除以它的标准误估计值。如果学生残差超过3，我们一般就认为可能是Outlier。在下图的右图中，我们发现这个Outlier的学生残差超过了6，而其他点均在 -2 与 2 之间。

如果你相信Outlier的出现是由于数据记录或收集时出错，那么解决办法就是移除这个观测值。但是你需要担心，Outlier的出现也可能是说明你的模型不完善，比如缺了一个影响因子。

5 High Leverage Points

High Leverage Points 仅指拥有异常的 x_i 的观测值。例如，下图中的观测值 41 就是一个 High Leverage Point。 左图中的红线为包含这个点的回归直线，蓝色虚线是剔除了这个点的回归直线。我们可以看到， High Leverage Point 的影响很大，因此找 High Leverage Point 很重要。

对于简单线性回归， High Leverage Point 很容易找到，我们可以直接从 X-Y 图中直接看出。但是多重线性回归存在多个因变量，有可能一个观测值的所有的自变量单独看都在合理的范围，但是这一组自变量的组合是异常的。下图中间的图就显示了这一点，这是一个有两个自变量(X1, X2)的回归。很多观测值的两个自变量均落在蓝色虚线框定的椭圆中，但是这个红色的点不在这个范围内，虽然这个点的 X1 和 X2 值单独看都不是异常的值。**因此对于多重回归，如果仅仅看单个自变量是否存在 High Leverage Point，我们可能会漏掉一些 High Leverage Point。 ** 但是对于超过2个变量的多重回归模型，如果能够同时查看多个自变量组合的 High Leverage Point 就是一个问题，因为你没法针对2维以上的变量进行画图。

为了确定一个观测值的 leverage，我们可以计算一个 leverage statistic 。这个值越大说明 leverage 程度越高。对于一个简单线性回归，公式如下（多重回归的公式未给）：

这个统计量的值在 1/n 与 1 之间，而且所有的观测值的 leverage statistic 的均值是 (p+1)/n 。所以如果一个观测值的 leverage statistic 远超均值，我们就可以怀疑这是一个 High Leverage Point

上图 3.13 中的右图是 学生残差 vs h_i 的散点图。观测值41 的 h_i 和 学生残差均特别高，说明**这个点即是 High Leverage Point ，也是一个 outlier ** 。这种情况是非常危险的组合，会对回归结果造成很大的影响，相比之下观测值20 虽然是一个 outlier, 但是 leverage 很低，最起码对回归方程没什么影响。

6 Collinearity

Collinearity 表示两个或更多个自变量之间具有紧密的相关关系。Collinearity 可以用下图来表示，左图中的两个变量之间表现出来没有明显的关系，而右图的两个变量之间存在高度的相关关系（我不知道这个等高图为什么可以这么解释）。

Collinearity 的影响在于无非区分开多个自变量对因变量的影响程度。比如右图中，RSS最低的（β_Limit , β_Rating）有很多个，他们在一个狭长的椭圆当中，这导致了参数估计的不确定性。

Collinearity 会降低回归参数估计的准确性，导致斜率估计的标准误增加，因此 P 值会增加，会导致 power（正确地拒绝原假设的概率） 降低。下图就是两次回归得到的回归系数和 P值，我们可以看到由于存在 Collinearity ，导致 limit 因素第二次分析时变得不显著。

因此我们需要在进行回归分析前进行 Collinearity 的检测。一个简单的思路是看自变量之间的相关矩阵，看是否存在相关系数很高的情况。但是有可能存在一种情况，就是自变量之间不存在两两之间的强相关，但是存在三个及以上的变量之间的相关（例如 X1 = a + bX2 + cX3）。我们把这种情况称为 multicollinearity 。因此我们不采用自变量的相关矩阵，我们引入了一个新的统计量来检测 Collinearity，就是 VIF 。定义如下：

The VIF is the ratio of the variance of ˆβ_j when fitting the full model divided by the variance of ˆβ_j if fit on its own.

VIF 的最小值为1， 说明完全没有Collinearity ，一般来说自变量之间存在一定的 Collinearity , 一般我们认为 VIF 超过 5 或 10 就说明存在严重的 Collinearity 。对于每个变量的 VIF 计算公式为

在上面的例子中，对 age, rating 和 limit 这三个变量做 VIF 分析，计算结果分别为 1.01， 160.67 和 160.59 。

如果面临 Collinearity 的问题，有两个简单的思路：一是删除其中一个变量，这通常不会导致模型拟合程度的降低，因为变量之间的信息是冗余的；二是合并共显性的变量为一个变量。

R代码实现

简单线性回归

回归分析

1
2
3
4
5
6

# 简单回归
library(MASS)
library(ISLR)
## 采用MASS包的Boston数据集
lm.fit = lm(medv~lstat, data = Boston)
summary(lm.fit) #各种回归的信息

回归分析结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

Call:
lm(formula = medv ~ lstat, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.168  -3.990  -1.318   2.034  24.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.55384    0.56263   61.41   <2e-16 ***
lstat       -0.95005    0.03873  -24.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.216 on 504 degrees of freedom
Multiple R-squared:  0.5441,	Adjusted R-squared:  0.5432 
F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

只看系数

1
2
3

> coef(lm.fit)
(Intercept)       lstat 
 34.5538409  -0.9500494 

查看系数的置信区间（confidence interval ）

1
2
3
4

> confint(lm.fit)
                2.5 %     97.5 %
(Intercept) 33.448457 35.6592247
lstat       -1.026148 -0.8739505

画 x-y 散点图，并添加回归线

1
2

plot(Boston$lstat, Boston$medv)
abline(lm.fit)

我们可以查看检查线性回归是否正常的图

1
2

par(mfrow=c(2,2))
plot(lm.fit) # 直接画4副图

绘图如下：上方左图就是正常的残差图，可以查看非线性，残差异方差等，我们可以看到这里存在非线性关系；上方右图就是QQ图，用于检验残差是否符合正态分布；下方左图也是一个残差图，不过就是把纵坐标改成了 standardized residuals（缺公式）；下方右图用于判断 Outlier 和 High Leverage Point ，横坐标是 leverage， 通常大于 4/n (n 是样本点的数目)，就算 High Leverage Point ； 纵轴是residuals，通常小样本大于2，大样本大于4的算作outliers 。

多重线性回归

看一下数据结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

> 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 ...

两个回归因素，lstat 和 age

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

> lm.fit=lm(medv~lstat+age,data=Boston)
> summary(lm.fit)

Call:
lm(formula = medv ~ lstat + age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.981  -3.978  -1.283   1.968  23.158 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
age          0.03454    0.01223   2.826  0.00491 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.173 on 503 degrees of freedom
Multiple R-squared:  0.5513,	Adjusted R-squared:  0.5495 
F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

Boston 数据集总共有 13 个自变量，我们可以通过这种简写方式，直接使用所有的自变量（前提是除了因变量以外，其他全是自变量）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

> lm.fit=lm(medv~.,data=Boston)
> summary(lm.fit)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288    
chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229    
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,	Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

我们可以通过car包的 vif() 函数计算 各个因子的VIF，这里可以看到 rad 和 tax 这两个因素的 VIF 比较高。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

> library(car)
载入需要的程辑包：carData
Registered S3 methods overwritten by 'tibble':
  method     from  
  format.tbl pillar
  print.tbl  pillar
Warning messages:
1: 程辑包‘car’是用R版本4.0.3 来建造的 
2: 程辑包‘carData’是用R版本4.0.3 来建造的 
> vif(lm.fit)
    crim       zn    indus     chas      nox       rm      age      dis      rad      tax 
1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 7.484496 9.008554 
 ptratio    black    lstat 
1.799084 1.348521 2.941491 

使用除了某一个因子外的其他所有因子的回归（采用 - 号），例如剔除 age 因子。

1
2

>lm.fit1 = lm(medv∼.-age,data=Boston)
>summary(lm.fit1)

互作项

回归模型中 ， lstat:black 表示新增一个 lstat 和 black 的互作项。lstat*age 表示同时lstat, age 和他们之间的互作项，是 lstat+age+lstat:age 的简写。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

> summary(lm(medv~lstat*age,data=Boston))

Call:
lm(formula = medv ~ lstat * age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.806  -4.045  -1.333   2.085  27.552 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 36.0885359  1.4698355  24.553  < 2e-16 ***
lstat       -1.3921168  0.1674555  -8.313 8.78e-16 ***
age         -0.0007209  0.0198792  -0.036   0.9711    
lstat:age    0.0041560  0.0018518   2.244   0.0252 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.149 on 502 degrees of freedom
Multiple R-squared:  0.5557,	Adjusted R-squared:  0.5531 
F-statistic: 209.3 on 3 and 502 DF,  p-value: < 2.2e-16

自变量的非线性转换

我们可以创建一个新的自变量 X²，采用 I(X^2) 函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

> lm.fit2 = lm(medv~lstat+I(lstat^2), data = Boston)
> summary(lm.fit2)

Call:
lm(formula = medv ~ lstat + I(lstat^2), data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2834  -3.8313  -0.5295   2.3095  25.4148 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 42.862007   0.872084   49.15   <2e-16 ***
lstat       -2.332821   0.123803  -18.84   <2e-16 ***
I(lstat^2)   0.043547   0.003745   11.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.524 on 503 degrees of freedom
Multiple R-squared:  0.6407,	Adjusted R-squared:  0.6393 
F-statistic: 448.5 on 2 and 503 DF,  p-value: < 2.2e-16

我们可以看到二次项的P值也是显著的，说明模型拟合结果有所改善。我们可以使用 anova() 函数来比对两个模型。

这里 Model 1 表示原来的线性模型，Model 2 表示二次项回归模型。这个 anova() 函数执行了一个比对两个模型的假设检验，零假设是两个模型拟合程度一致，备择假设是更全的模型更好。这里F统计量是135， 相应的 P值也很低，说明二次项模型拟合程度更好。

1
2
3
4
5
6
7
8
9
10
11

> lm.fit = lm(medv~lstat, data=Boston)
> anova(lm.fit, lm.fit2)
Analysis of Variance Table

Model 1: medv ~ lstat
Model 2: medv ~ lstat + I(lstat^2)
  Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
1    504 19472                                 
2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

如果要建立三次型方程，我们可以新增一个变量 I(X^3) 。但是，这种写法比较啰嗦。一个更好的写法是使用poly() 函数来创建一个多项式。比如下面就会创建一个自由度为5的二次型拟合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

> lm.fit5 = lm(medv~poly(lstat,5),data=Boston)
> summary(lm.fit5)

Call:
lm(formula = medv ~ poly(lstat, 5), data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5433  -3.1039  -0.7052   2.0844  27.1153 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.215 on 500 degrees of freedom
Multiple R-squared:  0.6817,	Adjusted R-squared:  0.6785 
F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16

这里说明如果增加二次项，加到5， 可以增加模型的拟合程度。但是，进一步研究会发现，如果增加超过项数超过5的项，这些超过5的项的p值均不显著。

当然，除了使用二次项以外，我们还可以使用其他转换函数，比如log函数。

1

summary(lm(medv~log(rm), data=Boston))

分类因子

这里我们用 ISLR 包的 Carseats 数据集。我们想要预测的变量是 Sales

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

> library(ISLR)
> fix(Carseats)
> str(Carseats)
'data.frame':	400 obs. of  11 variables:
 $ Sales      : num  9.5 11.22 10.06 7.4 4.15 ...
 $ CompPrice  : num  138 111 113 117 141 124 115 136 132 132 ...
 $ Income     : num  73 48 35 100 64 113 105 81 110 113 ...
 $ Advertising: num  11 16 10 4 3 13 0 15 0 0 ...
 $ Population : num  276 260 269 466 340 501 45 425 108 131 ...
 $ Price      : num  120 83 80 97 128 72 108 120 124 124 ...
 $ ShelveLoc  : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
 $ Age        : num  42 65 59 55 38 78 71 67 76 76 ...
 $ Education  : num  17 10 12 14 13 16 15 10 10 17 ...
 $ Urban      : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
 $ US         : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...

这里包含了类似于 ShelveLoc 的分类变量。如果使用分类变量，R会自动转变为哑变量。下面我们进行包含了一些互作项的多重回归模型。

这里 ShelveLoc 生成了两个哑变量 ShelveLocGood 和 ShelveLocMedium ，Bad 水平视为基础水平。我们可以看到 Good 和 Medium 水平的斜率都是正数，且均显著，说明 Good 和 Medium 水平的因变量的均值均比 Bad 水平要高。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

> lm.fit=lm(Sales~.+Income:Advertising +Price:Age,data=Carseats)
> summary(lm.fit)

Call:
lm(formula = Sales ~ . + Income:Advertising + Price:Age, data = Carseats)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9208 -0.7503  0.0177  0.6754  3.3413 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         6.5755654  1.0087470   6.519 2.22e-10 ***
CompPrice           0.0929371  0.0041183  22.567  < 2e-16 ***
Income              0.0108940  0.0026044   4.183 3.57e-05 ***
Advertising         0.0702462  0.0226091   3.107 0.002030 ** 
Population          0.0001592  0.0003679   0.433 0.665330    
Price              -0.1008064  0.0074399 -13.549  < 2e-16 ***
ShelveLocGood       4.8486762  0.1528378  31.724  < 2e-16 ***
ShelveLocMedium     1.9532620  0.1257682  15.531  < 2e-16 ***
Age                -0.0579466  0.0159506  -3.633 0.000318 ***
Education          -0.0208525  0.0196131  -1.063 0.288361    
UrbanYes            0.1401597  0.1124019   1.247 0.213171    
USYes              -0.1575571  0.1489234  -1.058 0.290729    
Income:Advertising  0.0007510  0.0002784   2.698 0.007290 ** 
Price:Age           0.0001068  0.0001333   0.801 0.423812    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.011 on 386 degrees of freedom
Multiple R-squared:  0.8761,	Adjusted R-squared:  0.8719 
F-statistic:   210 on 13 and 386 DF,  p-value: < 2.2e-16