线性回归进行预测建模-R(结果解读)

预测模型对于预测未来结果和估算实际难以量化的指标非常有用。例如，数据科学家可以使用预测模型基于降雨和温度预测作物产量，或确定具有某些特征的患者是否更可能对新药物反应严重。

在我们具体讨论线性回归之前，让我们自己回顾一下典型的数据科学工作流程。很多时候，我们将从我们想要回答的问题开始，并执行以下操作：

• 收集与问题相关的一些数据（更多数据几乎总是更好）。

• 如果需要，可以将数据清理，扩充和预处理为方便的形式。

• 对数据进行探索性分析以更好地了解数据。

• 使用您的发现作为指南，构建数据某些方面的模型。

• 使用该模型回答您开始的问题，并验证您的结果。

线性回归是数据科学家用于预测建模的最简单和最常见的监督机器学习算法之一。 在这篇文章中，我们将使用线性回归来构建一个模型，该模型从人们易于研究的树木测量指标来预测樱桃树的体积。

我们将在此博客文章中使用R来探索此数据集并学习线性回归的基础知识。 如果您不熟悉R，我们建议从R Data Analyst路径获取R基础知识和R编程：中级课程。 它也有助于获得一些非常基本的统计知识，但如果你知道平均值和标准偏差是什么，你就可以跟随它。 如果您想自己练习构建模型和可视化，我们将使用以下R包：

• datasets
  该软件包包含各种练习数据集。 我们将使用其中一个“tree”来学习构建线性回归模型。

• GGPLOT2
  我们将使用这个流行的数据可视化包来构建我们的模型图。

• GGally
  该包扩展了ggplot2的功能。 我们将使用它来创建绘图矩阵，作为我们最初的探索性数据可视化的一部分。

• scatterplot3d
  我们将使用此包来显示具有多个预测变量的更复杂的线性回归模型。

无论如何，他们如何测量树木的体积？

树数据集包含在基础R的数据集包中，它将帮助我们回答这个问题。 由于我们正在使用现有的（干净的）数据集，上面的步骤1和2已经完成，因此我们可以跳到第3步中的一些初步探索性分析。这个数据集是什么样的？

data(trees) # 从datasets包中获取tree数据集
head(trees) # 查看数据集前面几行

str(trees) # 查看变量的结构

该数据集由31个描述黑樱桃树的3个数值变量组成：

• 树干周长（in）

• 高度（英尺）

• 体积（ft3）

这些指标对于研究树木生态学的林务员和科学家来说是有用的信息。 使用基本的林业工具测量树木高度和周长相当简单，但测量树木体积要困难得多。 如果您不想实际减少和拆除树木，您必须采取一些技术上具有挑战性和耗时的活动，如爬树和进行精确测量。 能够从高度和/或周长准确预测树木体积将是有用的。

为了确定我们是否可以建立预测模型，第一步是查看我们的预测变量和响应变量（在这种情况下是周长，高度和体积）之间是否存在关系。 让我们做一些探索性的数据可视化。 我们将使用GGally包中的ggpairs（）函数创建一个绘图矩阵，以查看变量如何相互关联。

ggpairs(data = trees, columns = 1:3, title = "trees data")

ggpairs（）函数给出了每个变量组合的散点图，以及每个变量的密度图和变量之间的相关强度。 如果您之前使用过ggplot2，这个符号可能看起来很熟悉：GGally是ggplot2的扩展，它提供了一个简单的界面来创建一些像这样复杂的图。

当我们查看这些图时，我们可以开始了解数据并提出问题。 相关系数提供有关变量与关系有多接近的信息; 相关系数越接近1，关系越强。 散点图让我们可视化变量对之间的关系。 散点图中点具有清晰的视觉模式（与看起来像无形的云相反）表明关系更强。

我们的问题：

哪些预测变量似乎与响应变量相关？
从ggpairs（）输出看，周长似乎与体积有关：相关系数接近于1，并且这些点似乎具有线性模式。 高度和体积之间可能存在关系，但看起来较弱：相关系数较小，散点图中的点更加分散。

变量之间关系的形状是什么？
这种关系似乎是线性的; 从散点图中，我们可以看到随着树周长的增加树木体积不断增加。

这种关系是强大的，还是淹没信号的数据中的噪音？
高度和体积之间的关系不是那么清楚，但是周长和体积之间的关系似乎很强。

现在我们对数据有了一个很好的整体把握，我们可以继续进行第4步并做一些预测建模。

形成一个假设

一个假设是对我们认为我们的数据发生了什么的有根据的猜测。 在这种情况下，让我们假设樱桃树的周长和体积是相关的。 我们形成的每个假设都有相反的结论：“零假设”（H0）。 在这里，我们的零假设是周长和体积无关。

在统计学中，零假设是我们使用我们的数据来支持或拒绝的假设; 我们不能说我们“证明”一个假设。 我们称之为周长和体积与我们的“替代”假设（Ha）相关。

总结一下：

• H0：周长和体积之间没有关系

• Ha：周长和体积之间有一些关系

我们的线性回归模型将用于检验我们的假设。 如果我们找到足够强的证据来拒绝H0，那么我们就可以使用该模型来预测来自周长的樱桃树体积。

构建线性回归模型

线性回归描述了感兴趣的响应变量（或因变量）与一个或多个预测变量（或独立变量）之间的关系。 它有助于我们将信号（我们可以从预测变量中学习关于响应变量的信息）与噪声分离（我们无法从预测变量中了解响应变量）。 随着我们的前进，我们将深入探讨模型如何做到这一点。

让我们直接进入并建立一个将树体积与周长相关联的线性模型。 R使用基本函数lm（）简化了这一点。

fit_1 <- lm(Volume ~ Girth, data = trees)

lm（）函数拟合我们数据的一条线，尽可能接近我们所有的31个观测值。 更具体地说，它以这样的方式拟合线，即点和线之间的平方差的总和最小化; 这种方法被称为“最小化最小二乘法”。

即使线性回归模型非常适合数据，拟合也并不完美。 我们的观测值与模型预测值之间的距离称为残差。

在数学上，我们可以将线性回归的等式写成：

• Y和X变量是我们与彼此相关的数据的响应和预测变量

• β0是表示模型截距或与y轴交叉的模型系数

• β1是表示模型斜率的模型系数，该斜率是给出关于线的陡度及其方向（正面或负面）的信息的数字。

• ε是包含我们无法在模型中捕获的可变性的误差项（X不能告诉我们关于Y的信息）

就我们的例子而言:
Tree Volume ≈ Intercept + Slope(Tree Girth) + Error

lm（）函数估计它适合我们数据的线性模型的截距和斜率系数。 有了模型，我们可以继续执行第5步，记住我们还有一些工作要做，以验证这个模型实际上适合数据的想法。

我们可以使用这个模型进行预测吗？

我们是否可以使用我们的模型进行预测取决于：

• 我们是否可以拒绝零假设，即我们的变量之间没有关系。

• 该模型是否适合我们的数据。

让我们使用summary（）调用模型的输出。 模型输出将为我们提供测试我们的假设所需的信息，并评估模型与我们的数据的匹配程度。

让我们通过输出来回答每个问题。

该假设是否得到支持？

系数估计和标准误差：

如果周长值为零，我们示例中的截距是预期的树体积。当然，我们不能拥有负数量的树，但稍后会分析。
我们的例子中的斜率是树周长对树体积的影响。我们看到，对于每增加一英寸的周长，树木体积增加5.0659英尺3。
系数标准误差告诉我们估计系数与我们的响应变量的实际平均值的平均变化。

t值：

这是一个测试统计量，用于衡量估计系数与零的标准偏差。

PR（> | T |）：
该数字是p值，定义为如果H0为真，则任何观察值的概率大于或者等于t值。 t统计量越大，p值越小。一般来说，我们使用0.05作为显着性的截止值;当p值小于0.05时，我们拒绝H0。
我们可以拒绝零假设，而是相信树宽和体积之间存在关系。

模型与数据的拟合程度如何？

残差：

输出的这一部分为我们提供了残差的总结（回想一下这些是我们的观察与模型之间的距离），它告诉我们一些关于我们的模型如何适合我们的数据。残差应该在零附近具有非常对称的分布。一般来说，我们正在寻找正态分布在零附近的残差（即钟形曲线分布），但重要的是它们没有视觉上明显的模式，这表明线性模型不适合数据。

我们可以使用ggplot2制作直方图以使其可视化。

ggplot(data=trees, aes(fit_1$residuals)) + 
  geom_histogram(binwidth = 1, color = "black", fill = "purple4") +
  theme(panel.background = element_rect(fill = "white"),
        axis.line.x=element_line(),
        axis.line.y=element_line()) +
  ggtitle("Histogram for Model Residuals")

我们的残差在0附近看起来非常对称，这表明我们的模型很好地拟合了数据。

Residual standard error:：

该术语表示我们的响应变量测量值偏离拟合线性模型（模型误差项）的平均量。

Degrees of freedom (DoF)：

对自由度的讨论可能变得相当技术性。出于本文的目的，将它们视为用于计算估计的独立信息的数量就足够了。 DoF与测量数量相关但不相同。

Multiple R-squared：

R2值衡量我们的数据与线性回归模型的接近程度。 R2值始终在0和1之间;接近1的数字代表合适的模型。 R2总是随着模型中包含更多变量而增加，因此包含调整后的R2以计算用于制作模型的独立变量的数量。

F statistic：

此测试统计信息告诉我们正在测试的依赖变量和自变量之间是否存在关系。通常，大F表示更强的关系。

p-value：

该p值与F统计量相关联，用于解释适合我们数据的整个模型的显着性。

让我们来看看我们的模型是否适合我们的周长和体积数据。 我们可以通过使用ggplot（）将线性模型拟合到我们数据的散点图中来实现：

ggplot(data = trees, aes(x = Girth, y = Volume)) + geom_point()  +
  stat_smooth(method = "lm", col = "dodgerblue3") +
  theme(panel.background = element_rect(fill = "white"),
        axis.line.x=element_line(),
        axis.line.y=element_line()) +
  ggtitle("Linear Model Fitted to Data")

线周围的灰色阴影表示0.95的置信区间，这是stat_smooth（）函数的默认值，它可以平滑数据以使模式更易于显示。该0.95置信区间是所有黑樱桃树的周长和体积的真实线性模型将在符合我们数据的回归模型的置信区间内的概率。

即使这个模型非常适合我们的数据，我们的观察结果仍然存在差异。这是因为世界通常不整洁。在我们的模型中，树木体积不仅仅是树木周长的函数，而且还有我们不一定需要数据量化的东西（树干形状之间的个体差异，林务员树干周长测量技术的微小差异）。

有时，这种可变性会掩盖响应和预测变量之间可能存在的任何关系。但在这里，我们数据中的信号足够强大，可以让我们开发一个有用的模型来进行预测。

使用我们的简单线性模型进行预测

我们的模型适合做出预测！ 各地的树科学家都很高兴。 假设我们有一个遗漏在数据集之外的树的周长，高度和体积数据。 我们可以使用这棵树来测试我们的模型。

我们的模型在预测树木的体积方面有多好？

我们将使用predict（）函数，一个通用的R函数，用于根据模型拟合函数的变量进行预测。 predict（）将我们的线性回归模型和我们想要响应变量值的预测变量的值作为参数。

predict(fit_1, data.frame(Girth = 18.2))

我们的体积预测为55.2 ft3。

这与我们的实际值接近，但是将高度（我们的另一个预测变量）添加到我们的模型中可能可以让我们做出更好的预测。

添加更多预测变量：多元线性回归

如果我们使用我们可用的所有信息（宽度和高度）来预测树木体积，也许我们可以提高模型的预测能力。重要的是，从帖子开始的五步过程实际上是一个迭代过程 - 在现实世界中，您将获得一些数据，构建模型，根据需要调整模型以改进它，然后可能添加更多数据并建立一个新的模型，等等，直到你对结果感到满意和/或确信你不能做得更好。

我们可以构建两个单独的回归模型并对其进行评估，但这种方法存在一些问题。首先，想象一下，如果我们有5个，10个甚至50个预测变量，那将是多么繁琐。其次，两个预测模型会给我们预测两个单独的体积，而不是我们所追求的单个预测。也许最重要的是，建立两个单独的模型并不能让我们在估算模型系数时考虑预测变量之间的关系。

在我们的数据集中，我们怀疑树高和周长是基于我们的初始数据探索相关联的。由于我们将在本文中进一步明确地看到这一点，忽略预测变量之间的这种相关性可能会导致关于它们与树体积的关系的误导性结论。

更好的解决方案是构建包含多个预测变量的线性模型。我们可以通过为我们的模型中每个额外的自变量添加斜率系数来做到这一点。

这很容易用lm（）函数：我们只需要添加其他预测变量。

fit_2 <- lm(Volume ~ Girth + Height, data = trees)
summary(fit_2)

我们可以从模型输出中看到，周长和高度都与体积显着相关，并且模型很好地适合我们的数据。 我们调整后的R2值也略高于模型fit_1的调整后的R2值。

由于我们在此模型中有两个预测变量，因此我们需要第三个维度来对其进行可视化。 我们可以使用包scatterplot3d创建一个漂亮的三维散点图：

首先，我们为预测变量（在我们的数据范围内）制作一个值网格。 expand.grid（）函数从因子变量的所有组合创建数据框。

Girth <- seq(9,21, by=0.5) ## make a girth vector
Height <- seq(60,90, by=0.5) ## make a height vector
pred_grid <- expand.grid(Girth = Girth, Height = Height)

接下来，我们基于预测变量网格对体积进行预测：

pred_grid$Volume2 <-predict(fit_2, new = pred_grid)

现在我们可以从预测器网格和预测的体积中制作一个三维散点图：

fit_2_sp <- scatterplot3d(pred_grid$Girth, pred_grid$Height, pred_grid$Volume2, angle = 60, color = "dodgerblue", pch = 1, 
                          ylab = "Hight (ft)", xlab = "Girth (in)", zlab = "Volume (ft3)" )

让我们看看这个模型如何预测树的体积。 这次，我们包括树的高度，因为我们的模型使用Height作为预测变量：

predict(fit_2, data.frame(Girth = 18.2, Height = 72))

这次，我们得到的预测体积为52.13 ft3。 这个预测比我们只有周长作为预测器的简单模型得到的树体积更接近真实值。但是，正如我们将要看到的，我们可能能够改进。

考虑交互

虽然我们已经进行了改进，但我们刚刚建造的模型仍然没有说明整个故事。 它假设周长对体积的影响与树高对体积的影响无关。 事实显然不是这样，因为树高和周长是相关的; 更高的树往往更宽，我们的探索性数据可视化表明了更多。 换句话说，当高度的斜率增加时，周长的斜率应该增加。

为了解释我们模型中预测变量的这种非独立性，我们可以指定一个交互项，它是作为预测变量的乘积计算的。

再一次，使用lm（）很容易构建这个模型：

fit_3 <- lm(Volume ~ Girth * Height, data = trees)
summary(fit_3)

请注意，“Girth Height”术语是我们模型中“Girth + Height + Girth Height”的简写。

正如我们所怀疑的那样，周长和高度的相互作用是显着的，这表明我们应该在我们用来预测树木体积的模型中包含相互作用项。 调整后的R2值接近1，F的大值和p的小值也支持这一决定，这表明我们的模型非常适合数据。

让我们看一下散点图，使用该模型可视化树体积的预测值。 我们可以使用为fit_2可视化生成的相同网格预测值：

Girth <- seq(9,21, by=0.5)
Height <- seq(60,90, by=0.5)
pred_grid <- expand.grid(Girth = Girth, Height = Height)

与我们如何可视化fit_2模型类似，我们将使用fit_3模型和交互项来预测预测变量网格中的体积值：

pred_grid$Volume3 <-predict(fit_3, new = pred_grid)

现在我们制作预测器网格和预测体积的散点图：

fit_3_sp <- scatterplot3d(pred_grid$Girth, pred_grid$Height, pred_grid$Volume3, angle = 60, color = "dodgerblue", pch = 1, ylab = "Hight (ft)", xlab = "Girth (in)", zlab = "Volume (ft3)")

最后，我们覆盖了观察到的数据：

fit_3_sp$points3d(trees$Girth, trees$Height, trees$Volume, pch=16)

在这张图片中有点难以看到，但这次我们的预测是在一些曲面而不是平面上。 现在是时候了：让我们用这个模型来预测树木的体积。

predict(fit_3, data.frame(Girth = 18.2, Height = 72))

我们使用第三个模型的预测值是45.89，最接近我们的真实值46.2 ft3。

关于预测模型的一些注意事项

记住您的数据范围
使用模型进行预测时，最好避免尝试推断远远超出用于构建模型的值范围。 为了说明这一点，我们试着估计一棵小树苗（一棵幼树）的体积：

predict(fit_3, data.frame(Girth = 0.25, Height = 4))

我们得到的预测体积为62.88 ft3，比我们数据集中的高大树木更大。 当然这没有意义。 请记住，我们进行准确预测的能力受到我们用于构建模型的数据范围的限制。

避免制作对您的数据集过于具体的模型

在我们在这篇文章中调查的简单示例数据集中，向我们的模型添加第二个变量似乎可以提高我们的预测能力。然而，当尝试具有许多差异变量的各种多元线性回归模型时，选择最佳模型变得更具挑战性。如果添加了太多不能提高模型预测能力的术语，我们就会冒险“过度拟合”我们的模型到我们的特定数据集。过度拟合特定数据集的模型失去了预测未来事件或拟合不同数据集的功能，因此不是非常有用。

虽然我们在本文中用于评估模型有效性的方法（调整后的R2，残差分布）对于了解模型与数据的拟合程度非常有用，但将模型应用于数据集的不同子集可以提供有关模型执行情况的信息。在实践中。这种被称为“交叉验证”的方法通常用于测试预测模型。在我们的示例中，我们使用三个模型中的每一个来预测单个树的体积。但是，如果我们要构建更复杂的模型，我们希望将数据的子集用于交叉验证。

下一步

我们使用线性回归建立模型来预测来自两个连续预测变量的连续响应变量，但线性回归是许多其他常见场景的有用预测建模工具。

下一步，尝试构建线性回归模型，以预测来自两个以上预测变量的响应变量。想一想如何决定在回归模型中包含哪些变量;你怎么知道哪些是重要的预测因子？预测变量之间的关系如何干扰这一决定？ R中的数据集可用于处理多个线性回归问题，包括：airquality，iris和mtcars。

从数据构建模型的另一个重要概念是使用从现有数据计算的新预测变量来增加数据。这称为特征工程，您可以使用自己的专业知识来了解与问题相关的其他内容。例如，如果您正在查看具有时间戳作为变量之一的银行交易数据库，那么一周中的某一天可能与您想要回答的问题相关，因此您可以从时间戳计算并添加它将数据库作为新变量。这是一个复杂的主题，添加更多预测变量并不总是一个好主意，但是当您了解有关建模的更多信息时，您应该记住这一点。

在这篇文章中使用的树数据集中，你能想到你可以从周长和高度计算的任何额外数量来帮助你预测体积吗？ （提示：回想一下你学习各种几何形状体积的公式时，想一想树的样子。）

最后，虽然我们专注于连续数据，但也可以扩展线性回归以从分类变量进行预测。尝试使用线性回归模型来预测分类和连续预测变量的响应变量。 R中有一些数据集特别适合这个练习：ToothGrowth，PlantGrowth和npk。

附录：完整代码

library(ggplot2)
library(GGally)
library(scatterplot3d)

# 线性回归模型
data(trees) # 从datasets包中获取tree数据集
head(trees) # 查看数据集前面几行

str(trees) # 查看变量的结构

ggpairs(data = trees, columns = 1:3, title = "trees data")

# 构建线性回归模型
fit_1 <- lm(Volume ~ Girth, data = trees)

summary(fit_1)

# 残差分布的可视化
ggplot(data=trees, aes(fit_1$residuals)) + 
  geom_histogram(binwidth = 1, color = "black", fill = "purple4") +
  theme(panel.background = element_rect(fill = "white"),
        axis.line.x=element_line(),
        axis.line.y=element_line()) +
  ggtitle("Histogram for Model Residuals") 

ggplot(data = trees, aes(x = Girth, y = Volume)) + geom_point()  +
  stat_smooth(method = "lm", col = "dodgerblue3") +
  theme(panel.background = element_rect(fill = "white"),
        axis.line.x=element_line(),
        axis.line.y=element_line()) +
  ggtitle("Linear Model Fitted to Data")

predict(fit_1, data.frame(Girth = 18.2))

fit_2 <- lm(Volume ~ Girth + Height, data = trees)
summary(fit_2)

Girth <- seq(9,21, by=0.5) ## make a girth vector
Height <- seq(60,90, by=0.5) ## make a height vector
pred_grid <- expand.grid(Girth = Girth, Height = Height) 

pred_grid$Volume2 <-predict(fit_2, new = pred_grid)

fit_2_sp <- scatterplot3d(pred_grid$Girth, pred_grid$Height, pred_grid$Volume2, angle = 60, color = "dodgerblue", pch = 1, 
                          ylab = "Hight (ft)", xlab = "Girth (in)", zlab = "Volume (ft3)" )


predict(fit_2, data.frame(Girth = 18.2, Height = 72))


fit_3 <- lm(Volume ~ Girth * Height, data = trees)
summary(fit_3)

Girth <- seq(9,21, by=0.5)
Height <- seq(60,90, by=0.5)
pred_grid <- expand.grid(Girth = Girth, Height = Height)

pred_grid$Volume3 <-predict(fit_3, new = pred_grid)

fit_3_sp <- scatterplot3d(pred_grid$Girth, pred_grid$Height, pred_grid$Volume3, angle = 60, color = "dodgerblue", pch = 1, ylab = "Hight (ft)", xlab = "Girth (in)", zlab = "Volume (ft3)")

fit_3_sp$points3d(trees$Girth, trees$Height, trees$Volume, pch=16)

predict(fit_3, data.frame(Girth = 18.2, Height = 72))

predict(fit_3, data.frame(Girth = 0.25, Height = 4))

作者：Rose Martin