第2章 线 性 回 归

我们从第1章了解到回归问题会预测一个数值型的输出。最简单和最常见的回归类型就是线性回归。本章要探讨为什么线性回归如此常用，以及它的局限性和扩展问题。

2.1　线性回归入门

在线性回归（linear regression）中，输出变量是通过输入特征的一个线性加权组合来预测的。下面是简单线性模型的一个示例：
上述模型实质上表达的是，我们要估算一个用表示的输出，它是由字母标记的一个预测变量（也就是特征）的一个线性函数。由希腊字母表示的项是模型的参数，被称为回归系数（regression coefficient）。一旦我们训练好模型并给这些参数选定了值，我们就可以通过在方程里直接代入的任何值对输出变量进行预测。另一个线性模型的示例由下面的方程给出，这次的模型带有三个特征，回归系数也已经赋值了：

和前面那个方程一样，在该方程中我们也可以看到系数比特征数多1个。这个附加的系数0.49被称为截距（intercept），它是当所有输入特征值都为0时模型的期望值。其他的参数可以解释为某特征每增加1个单位时输出值的预期变化。例如，在前面的方程中，如果特征的值增加1个单位，输出的期望值就会增加1.91个单位。类似地，特征的值增加1个单位会让输出减少7.56个单位。在简单的一维回归问题中，我们可以把输出绘制在图的轴，把输入特征绘制在轴。在这个示例中，模型预测了两个轴之间的一种直线关系，其中0代表该直线和轴交叉或相交的交点，而1则代表该直线的斜率。我们通常把这种一个特征（因而有两个回归系数）的情况称为简单线性回归（simple linear regression），把两个或多个特征的情况称为多元线性回归（multiple linear regression）。

线性回归的假设

在深入钻研线性回归模型的训练方法及其性能的细节之前，我们要来看一下模型的假设。模型假设主要描述了模型对于我们尝试预测的输出变量y的判断。具体来说，线性回归模型会假设输出变量是一组特征变量的加权线性函数。此外，该模型还假设对于特征变量的固定值，输出是具有常数方差的正态分布。这就等于说，以两个输入特征为例，该模型假设的是，真实的输出变量y可以用如下的一个方程表示：

这里，代表一个误差项，它是一个均值为0且方差为常数2的正态分布：

我们可能听说过更正式地描述常数方差概念的术语—同方差性（homoscedasticity）。说到同方差性或者常数方差（constant variance），我们指的是这样的事实：误差部分的方差不会随着输入特征的值或水平（level）而变化。在下面的绘图中，我们可视化了一个具有异方差（heteroskedastic）误差（也就是不具有常数方差的误差）的线性关系虚构示例。其中，输入特征值较小的数据点离直线较近，因为在绘图的这个区域方差比较低，而较大的输入特征值就离直线更远，因为它的方差更大。

异方差误差

#上图的作图代码
set.seed(5427395)
xplt <- runif(500, 0, 200)
eplt <- rnorm(500, mean = 0, sd = 2.0)
yplt <- 1.67 * xplt - 2.93 + (xplt ^ 1/20) * eplt #注意误差项的设置
dfplt <- data.frame(yplt, xplt)
ggplot(dfplt, aes(x = xplt, y = yplt))+
  geom_point(shape = 19, alpha = 0.5)+
  xlab('x1') + ylab('y')+
  geom_smooth(method = lm)

项是实际函数的一个不能化简的误差分量，可以用来表示随机误差，例如在特征值里的测量误差。在训练一个线性回归模型时，即使得到了所有合适的特征和充足的数据，并且要建模的系统也确实是线性的，但我们还是会预期在对输出的估算中观测到一定数量的误差。换言之，即使是线性的实际函数，我们还是会预期，一旦找到了对训练样本拟合最佳的一条直线，但由于该误差分量体现的这个固有方差，这条直线也不会正好穿过所有（或者甚至不会穿过任何一个）数据点。不过，需要记住的关键一点是，在这个理想化的场景里，因为误差分量的均值为0且方差为常数，所以在具有足够大的样本的前提下，训练准则让我们能接近回归系数的真实值，因为这些误差会互相抵消。

另一个重要的假设和误差项的独立性有关。这意味着我们不希望和某条具体的观测数据相关的残差（residual）或误差项以某种方式也和另一条观测数据具有相关性。作为测量过程有误的典型结果，当观测数据相互之间具有函数关系时，这个假设就不成立了。如果取出训练数据的一部分，把所有的特征值和输出值乘以2，并把这些新数据点加入训练数据，我们就会产生一种拥有更大数据集的错觉，不过，这样做的结果是导致存在误差项互相有依赖关系的配对观测数据，从而使我们的模型假设不能成立。附带再说一句，对于任何模型，这种人工扩大数据集的方式都永远是不可取的。类似地，如果观测数据以某种方式通过某个未测量的变量相互关联，也会产生互相关联的误差项。例如，如果我们在测量某个生产线零件的故障率，那么来自同一个工厂的零件可能会在误差里具有相关性，比如，来源于在组装过程的不同标准和规程。因此，如果不把工厂作为一个特征加进来，我们就会在样本里看到，来自同一个工厂的零件的观测数据之间会存在相互关联的误差。实验设计（experimental design）的研究内容就是关于识别并消除误差项里的相关性的，但这个主题超出了本书的范围。

最后还有一个重要的假设，它表达了这样一个概念：特征本身在统计学上是相互独立的。这里有必要阐明的是，在线性模型中，虽然输入特征必须是线性加权的，但是它们本身可以是另一个函数的输出。举例说明，如下方程式是一个具有三个特征sin(z1), ln(z2)和exp(z3)的线性模型，这看起来也许有点令人吃惊：

通过对输入特征进行一些变换并替换到模型中，我们就可以看出它其实是一个线性模型：

线性模型的分解

现在，我们有了一个更易于识别为线性回归模型的方程式。如果前面的示例让我们相信几乎任何方程都可以变换为一个线性模型，那么下面两个示例会有力地说服我们，实际上并非如此：

非线性模型

由于第一个回归系数（）的存在，上面两个模型都不是线性模型。第一个模型之所以不是线性模型，是因为成为第一个输入特征的指数。在第二个模型里，在一个正弦函数的内部。从这些示例要吸取的重要教训是，有些情况下，为了让数据符合线性模型，我们可以对输入特征进行变换，但是，需要小心从事，让回归系数始终是变换后新特征的线性权值。

2.2　简单线性回归

在着眼于某些真实环境的数据集之前，尝试在人造数据上训练模型是非常有帮助的。在这样的人造场景里，我们事先就知道了实际输出函数是什么，而这对于真实环境的数据来说通常是不成立的。进行这种练习的好处是，它会让我们对自己的模型在所有假设都完全成立的理想场景下的工作情况有清楚的了解，而且它有助于对具备理想的线性拟合时发生的情况进行可视化。我们先模拟一个简单线性回归模型。后面的R语言代码片段会为下面这个只有1个输入特征的线性模型创建一个带有100条模拟观测数据的数据框：

下面是用于这个简单线性回归模型的代码：

set.seed(5427395)
nObs <- 100
x1minrange <- 5
x1maxrange <- 25
x1 <- runif(nObs, x1minrange, x1maxrange)
e <- rnorm(nObs, mean = 0, sd = 2.0)
y <- 1.67 * x1 - 2.93 + e
df <- data.frame(y, x1)

对于输入特征，我们要从一个均匀分布中随机抽样一些点。我们用均匀分布让数据点理想地分散开。注意，最后那个df数据框的用途是模拟我们在实践中会获得的一个数据框，因此，其中并没有包含误差项e，因为它在实际环境里是无法得到的。

当利用类似上述数据框里的某些数据训练一个线性模型的时候，实际上希望的是能产生一个和数据背后的模型具有相同系数的线性模型。换言之，原始的系数就定义了一个总体回归线（population regression line）。在这种情况下，总体回归线代表了数据实际对应的模型。一般来说，我们会发现自己在试图对一个未必线性的函数进行建模。在这种情况下，仍然可以把总体回归线定义为最有可能的线性回归线，但一个线性回归模型显然不会有实际模型那样良好的表现。

估算回归系数

对于简单回归模型来说，训练模型的过程相当于根据数据集对两个回归系数进行估算。正如从之前构建的数据框里所看到的，我们的数据实际上是一系列的观测数据，其中每一个数据是一对值()，其中第一个元素是输入特征值，第二个元素是它对应的输出标签。其实，对于简单线性回归，有可能写出两个方程，用来计算这两个回归系数。我们不会直接给出这些方程，而是先要用简短篇幅来回顾一些读者在之前很可能遇到过的非常基础的统计量，因为它们很快就会起重要作用。

一组值的均值（mean）就是这些值的平均值，它经常被描述为位置的衡量指标，用来大致表示这些值在衡量它们的标尺上居中的位置。在统计学专业里，一个随机变量的平均值经常被称为期望值（expectation），因此我们经常可以看到一个随机变量X的均值被标记为。另一个普遍使用的记号是上横杠标记（bar notation），即可以通过给一个变量顶部加上一个横杠表示对一个变量取平均值的概念。要对此进行图示，下面的两个方程就表示了输出变量y和输入特征x的均值：

第二个很常见的统计量应该也是大家所熟悉的，它就是变量的方差（variance）。方差衡量的是每个变量值和均值之间距离的平方的平均值。从这个意义上说，它是一个离散度的衡量指标，因此低方差意味着大部分值是靠近均值进行聚集的，而更大的方差则会导致变量值是扩散开的。注意，方差的定义中包含了均值的定义，因此，我们会看到在下列显示输入特征x的方差的方程里用到了上部带一个横杠的x变量：

最后，用下列方程来定义两个随机变量x和y之间的协方差（covariance）：

从前面这个方程可以清楚地看到，前面定义的方差实际上就是协方差在两个变量相同时的一种特殊情况。协方差衡量的是两个变量相互之间关联的紧密程度，可以是正值或负值。正的协方差表明了正相关性，也就是说，当一个变量增大时，另一个也会增大。负的协方差则代表相反的含义，当一个变量增大时，另一个往往会减小。当两个变量从统计学角度是互相独立因而不相关的，它们的协方差就是0（不过需要说明的是：反过来，协方差为0并不一定能得出统计学角度的相互独立关系）。

掌握了这些基本概念之后，现在我们就可以给出简单线性回归里两个回归系数的近似计算公式了：

第一个回归系数等于输出与输入特征的协方差除以输入特征的方差所得到的比例。注意，如果输出变量和输入特征相互独立，则协方差为0，因而线性模型就是一条没有斜率的水平直线。在实践中要注意的是，即使两个变量从统计学角度是互相独立的，通常还是会得到少量的协方差，这是由于误差的随机性所致，因此如果要训练一个线性回归模型来表示它们之间的关系，第一个回归系数一般会是一个非零值。后面我们会看到，如何运用显著性检验（significance test）来检测出我们不应该纳入模型中的特征。

要在R语言中实现线性回归，进行这些计算并不是必需的，因为R语言提供了lm()函数，它会为我们构建一个线性回归模型。下面的代码样例会用到之前创建的df数据框，并计算出回归系数：

myfit <- lm(y ~ x1, df)
myfit

lm函数

在第一行我们可以看到lm()函数的用法：首先指定一个公式，然后接着给出data参数，它在示例中就是df数据框。对于简单线性回归，给lm()函数指定公式的语法是：输出变量的名字后面跟一个波浪符（~），然后再跟一个输入特征的名字。在本章后续内容里，当我们学习多元线性回归时，会看到如何指定更复杂的公式。最后，输出内容给出的是两个回归系数的值。注意，0系数被标记为截距（Intercept），而1系数是用线性模型方程里对应特征的名字（在这个示例中是x1）来标记的。
下图在同一张图上显示了总体回归线和近似线：

近似线和总体回归线

#上图代码
myfit <- lm(y ~ x1, df)
ggplot(df, aes(x = x1, y = y))+
  geom_point(shape = 19, alpha = 0.7)+
  xlab('x1') + ylab('y')+
  geom_smooth(method = lm)

正如我们所见，这两条线互相非常贴近，以至于几乎难以区分，这表明该模型对于实际总体回归线的近似是非常接近的。根据第1章的内容，我们知道可以用均方差作为规范来表示该模型拟合当前数据集的紧密程度，以及该模型将来会拟合类似测试集的紧密程度。我们在本章要了解它和其他几种模型性能及质量的衡量指标，但我们首先要把这个回归模型广义化，用来处理超过1个的输入特征。

2.3　多元线性回归

只要有多于一个输入特征，并且要构建一个线性回归模型，我们就处于多元线性回归的领域了。具有k个输入特征的多元线性回归模型的一般方程如下所示：

关于模型和误差分量的假设还是和简单线性回归的一样，记住，因为现在有了超过1个的输入特征，我们假设它们是相互独立的。我们在讲解多元线性回归时不会再使用模拟数据，而是要分析两套实际环境下的数据集。

2.3.1　预测CPU性能

我们的第一个实际环境下的数据集由研究者Dennis F. Kibler,David W. Aha,和Marc K. Albert于1989年在一篇论文中所提出，该论文名为《基于实例的实值属性预测》（Instance-based prediction of real-valued attributes），发表于《计算智能杂志》（Journal of Computational Intelligence）。该数据包含了不同CPU型号的特征，例如周期时间和缓存数量。在选择不同的处理器时，要考虑到所有这些因素，但理想情况下，我们希望在单一的数值标尺上对处理器进行比较。因此，我们经常会开发程序，用标杆来比较某个CPU的相对性能。数据集还带有已发布的CPU相对性能数据，我们的目标就是把这些可以获得的CPU特性作为预测所采用的特征。该数据集可以从UCI机器学习资源库（UCIMachine Learning Repository）在线获取，它的链接是：http://archive.ics.uci.edu/ml/datasets/Computer+Hardware。

UCI机器学习资源库是一个非常好的在线资源，它存放了大量数据集，其中很多都经常被各种著作和教材的作者引用。花费精力去熟悉这个网站及其数据集是非常值得的。有一种学习预测分析技术的好方法，就是把你在本书中学习的技术实际运用到不同的数据集上，而UCI资源库恰恰为此提供了很多这样的数据集。

machine.data文件包含了我们所有的数据，它使用逗号分开的格式（类似于.CSV文件），每一行对应一个CPU型号。我们要把它导入R并标记所有的列。注意，里面一共有10列，但我们的分析不需要前两列，因为它们不过是CPU的品牌和型号名。类似地，最后一列是研究者们自己得出的相对性能预测估算，实际输出变量PRP在第9列。我们会把需要的数据保存在一个名为machine的数据框里：

machine <- read.csv('YuCeFenXiRawData/machine.data', header = F)
names(machine) <- c('VENDOR', 'MODEL', 'MYCT', 'MMIN', 'MMAX', 'CACH', 'CHMIN', 'CHMAX',
                    'PRP', 'ERP')
machine = machine[, 3:9]
head(machine)

machine数据框的内容

该数据集还带有各列数据的定义：

列　　名	定　　义
MYCT	机器周期时间，以纳秒为单位
MMIN	最小主内存，以KB为单位
MMAX	最大主内存，以KB为单位
CACH	缓存，以KB为单位
CHMIN	最小通道数，以个数为单位
CHMAX	最大通道数，以个数为单位
PRP	发布的相对性能（即输出变量）

该数据集不包含缺失值，所以就不需要删除或修改任何观测数据。我们会注意到的一个问题是我们只有大约200个数据点，一般来说这算是非常小的样本了。尽管如此，我们还是会继续把数据划分为训练集和测试集，划分比例是85% / 15%，如下所示：

library(caret)
set.seed(4352345)
machine_sampling_vector <- createDataPartition(machine$PRP, p = 0.85, list = F)
machine_train <- machine[machine_sampling_vector, ]
machine_train_features <- machine[, 1:6]
machine_train_labels <- machine[machine_sampling_vector, 'PRP']
machine_test <- machine[-machine_sampling_vector, ]
machine_test_labels <- machine[-machine_sampling_vector, 'PRP']

既然已经把数据设置好了，我们通常就要继续分析和检查针对线性回归的假设是否成立。例如，我们会希望了解是否存在任何高度相关的特征。为此，可以用caret包的cor()函数构建一个相关度矩阵，并利用findCorrelation()函数获得关于哪些特征应该去掉的建议：

machine_correlations <- cor(machine_train_features)
findCorrelation(machine_correlations)
findCorrelation(machine_correlations, cutoff = 0.75)
cor(machine_train$MMIN, machine_train$MMAX)

findCorrelation函数

在采用代表高度相关性的默认阈值0.9的情况下，我们发现没有哪个特征应该去掉。当把这个阈值减小到0.75的时候，我们会看到caret包推荐去掉第三个特征（MMAX）。正如前面代码的最后一行所显示的，该特征和MMIN之间的相关度为0.768。虽然这个值并不是特别大，但它还是大到让我们会在某种程度上担心它会影响我们的模型。当然，直观地说，如果我们看一下输入特征的定义，当然我们往往会预期某个具有相对较高的最小主内存的型号也会具有相对较高的最大主内存。线性回归有时候对于具有相关性的变量也能给出一个较好的模型，但如果变量是不相关的，我们会期望得到更好的结果。就目前而言，我们决定还是保留数据集里的所有特征。

2.3.2　预测二手汽车的价格

我们的第二个数据集是在caret包里包含的汽车数据框，它是由Shonda Kuiper于2008年从Kelly Blue Book网站www.kbb.com收集的。这是一个获取二手车可靠价格的在线资源。该数据集包括了804种通用汽车（GM）品牌的汽车，所有汽车都是2005年出厂的。它包含了一组汽车属性，例如里程数、发动机排量，以及建议的销售价格。其中很多特征是二元的指示变量，例如Buick（别克商标）特征，它代表具体的某辆车的商标是否是Buick。所有的二手车在标价时都处于完好状态，使用时间少于一年，因此车况没有作为一个特征被纳入。我们对这个数据集的目标是创建一个能根据这些属性预测二手车售价的模型。特征的定义如下所示：

列名	定义
Price	以美元为单位的建议零售价（输出变量）
Mileage	车辆已经行驶的英里数
Cylinder	车辆发动机里的气缸数量
Doors	车门的数量
Cruise	代表车辆是否具有定速巡航控制的指示变量
Sound	代表车辆是否具有升级的音响设备的指示变量
Leather	代表车辆是否具有真皮座椅的指示变量
Buick	代表车辆的商标是否是别克（Buick）的指示变量
Cadillac	代表车辆的商标是否是凯迪拉克（Cadillac）的指示变量
Chevy	代表车辆的商标是否是雪佛兰（Chevy）的指示变量
Pontiac	代表车辆的商标是否是庞蒂亚克（Pontiac）的指示变量
Saab	代表车辆的商标是否是萨博（Saab）的指示变量
Saturn	代表车辆的商标是否是土星（Saturn）的指示变量
convertible	代表车辆的类型是否是敞篷车的指示变量
coupe	代表车辆的类型是否是双门轿车的指示变量
hatchback	代表车辆的类型是否是斜背式汽车的指示变量
sedan	代表车辆的类型是否是四门轿车的指示变量
wagon	代表车辆的类型是否是旅行轿车的指示变量

cars数据集

和前面关于CPU数据的machine数据集的处理方式一样，我们要分析输入特征的相关系数：

library(caret)
data(cars)
cars_train_features <- cars[, 2:18]
cars_cor <- cor(cars_train_features)
findCorrelation(cars_cor)
findCorrelation(cars_cor, cutoff = 0.75)
cor(cars$Doors, cars$coupe)
table(cars$Doors, cars$coupe)

cars数据集输入特征的相关系数

类似对前面CPU数据的machine数据集的处理，我们把caret包中findCorrelation()函数的cutoff参数设定为0.75，就得到了特征之间的相关情况。通过直接查看相关矩阵，我们发现在Doors特征和coupe特征之间存在相对较高程度的相关性。通过对这两个特征进行交叉分析（cross-tabulating），可以看到为什么会存在这种情况。如果知道某个汽车的类型是双门轿车，那么车门的数量就总是等于2。如果该汽车不是双门轿车，那么它最有可能的车门数是4。

这些二手车数据的另一个有问题的表现是，某些特征是其他一些特征的完全线性组合。这个问题是通过使用caret包里的findLinearCombos()函数发现的：

findLinearCombos(cars)

共线性问题

在这里，我们获得的指示是去掉coupe列和wagon列，它们分别是第15和第18个特征，因为它们恰好是其他特征的组合。我们会从我们的数据框里去掉它们，这样还能消除之前看到的相关性问题。

下一步，我们要把数据划分为训练集和测试集：

cars <- cars[, c(-15, -18)]
set.seed(232455)
cars_sampling_vector <- createDataPartition(cars$Price, p = 0.85, list = F)
cars_train <- cars[cars_sampling_vector,]
cars_train_features <- cars[, -1]
cars_train_labels <- cars[cars_sampling_vector, 'Price']
cars_test <- cars[-cars_sampling_vector, ]
cars_test_labels <- cars[-cars_sampling_vector, 'Price']

既然已经准备好了数据，我们就要创建一些模型了。