R-统计与机器学习研讨会-五-

R 统计与机器学习研讨会（五）

原文：annas-archive.org/md5/e61f21209b8b5ad343780cefdc55a102

译者：飞龙

协议：CC BY-NC-SA 4.0

第十二章：R 中的线性回归

在本章中，我们将介绍线性回归，这是一种基本的统计方法，用于建模目标变量与多个解释变量（也称为独立变量）之间的关系。我们将从简单线性回归开始，然后扩展到多元线性回归的概念。我们将学习如何估计模型系数，评估拟合优度，并使用假设检验测试系数的显著性。此外，我们还将讨论线性回归背后的假设，并探讨解决潜在问题（如非线性、交互效应、多重共线性异方差性）的技术。我们还将介绍两种广泛使用的正则化技术：岭回归和最小绝对收缩和选择算子（lasso）惩罚。

到本章结束时，您将学习线性回归的核心原则，其扩展到正则化线性回归，以及涉及的实现细节。

在本章中，我们将涵盖以下主题：

• 介绍线性回归

• 介绍惩罚线性回归

• 使用岭回归进行工作

• 使用 lasso 回归进行工作

要运行本章中的代码，您需要拥有以下包的最新版本：

• ggplot2, 3.4.0

• tidyr, 1.2.1

• dplyr, 1.0.10

• car, 3.1.1

• lmtest, 0.9.40

• glmnet, 4.1.7

请注意，前面提到的包的版本是在编写本章时的最新版本。本章的所有代码和数据均可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_12/working.R找到。

介绍线性回归

线性回归的核心是拟合一条直线——或者更一般地说，一个超平面——到数据点。这种拟合旨在最小化观察值和预测值之间的偏差。对于简单线性回归，一个目标变量由一个预测变量回归，目标是拟合一条最佳模拟两个变量之间关系的直线。对于多元线性回归，存在多个预测变量，目标是拟合一个最佳描述变量之间关系的超平面。这两个任务都可以通过最小化预测值和相应目标之间的偏差度量来实现。

在线性回归中，获得最佳模型意味着识别定义目标变量和输入预测变量之间关系的最佳系数。这些系数代表与相关预测变量的单位变化相关的目标变量的变化，假设其他所有变量保持不变。这使我们能够量化变量之间关系的大小（系数的大小）和方向（系数的符号），这些可以用于推断（强调可解释性）和预测。

当涉及到推断时，我们通常会观察输入变量单位变化对目标变量的相对影响。此类解释建模的例子包括营销支出如何影响季度销售额、吸烟状态如何影响保险费率、以及教育如何影响收入。另一方面，预测建模侧重于预测目标量。例如，根据营销支出预测季度销售额、根据投保人的个人资料信息（如年龄和性别）预测保险费率，以及根据某人的教育、年龄、工作经验和行业预测收入。

在线性回归中，预期的结果被建模为所有输入变量的加权总和。它还假设输出的变化与任何输入变量的变化成线性比例。这是回归方法的最简单形式。

让我们从简单线性 回归（SLR）开始。

理解简单线性回归

SLR（简单线性回归）是一种强大且广泛使用的统计模型，它规定了两个连续变量之间的关系，包括一个输入和一个输出。它使我们能够理解响应变量（也称为因变量或目标变量）如何随着解释变量（也称为自变量或输入变量）的变化而变化。通过将一条直线拟合到观察到的数据，SLR 量化了两个变量之间线性关联的强度和方向。这条直线被称为 SLR 模型。它使我们能够进行预测并推断预测变量对目标变量的影响。

具体来说，在 SLR 模型中，我们假设目标变量（y）和输入变量（x）之间存在线性关系。该模型可以用以下数学公式表示：

y = β 0 + β 1 x + ϵ

在这里，y 被称为响应变量、因变量、解释变量、预测变量、目标变量或回归量。x 被称为解释变量、自变量、控制变量、预测变量或回归变量。β 0 是线性线的截距，表示当 x 为 0 时 y 的期望值。β 1 是斜率，表示 x 增加一个单位时 y 的变化。最后，ϵ是随机误差项，它解释了目标变量 y 中预测变量 x 无法解释的变异性。

线性回归模型（SLR）的主要目标是估计 β₀ 和 β₁ 参数。一组最优的 β₀ 和 β₁ 参数将最小化观测目标值 y 和预测值 ˆy 之间的总平方偏差。这被称为最小二乘法，我们寻求最优的 β₀ 和 β₁ 参数，它们对应于最小平方误差和（SSR）：

minSSR = min∑i=1nui² = min(yi − ˆyi)²

在这里，每个残差 ui 是观测值 yi 和其拟合值 ˆyi 之间的差异。简单来说，目标是找到最接近数据点的直线。图 12*.1* 展示了一组数据点（蓝色）和线性模型（红色）：

图 12.1 – 线性回归模型，其中线性模型以线的形式出现，并通过最小化 SSR 进行训练

一旦我们估计了模型系数 β₀ 和 β₁，我们就可以使用该模型对变量之间线性关系的强度进行预测和推断。这种线性关系表明了拟合优度，通常使用确定系数，或 R² 来衡量。R² 的范围从 0 到 1，并量化了 y 的总变异中可以由 x 解释的部分。它定义如下：

R² = 1 − ∑i(yi − ˆyi)² / ∑i(yi − _y)²

这里，_y 表示观测目标变量 y 的平均值。

此外，我们还可以使用假设检验来检验得到的系数 β₀ 和 β₁ 的显著性，从而帮助我们确定变量之间的观测关系是否具有统计学意义。

让我们通过一个使用模拟数据集构建简单线性模型的例子来进行分析。

练习 12.1 – 构建线性回归模型

在这个练习中，我们将演示在 R 中实现线性回归模型。我们将使用内置函数和包的组合，使用模拟数据集来完成这项任务：

1. 模拟一个数据集，使得响应变量 Y 线性依赖于解释变量 X，并添加了一些噪声：

   
   # Set seed for reproducibility
   set.seed(123)
   # Generate independent variable X
   X = runif(100, min = 1, max = 100) # 100 random uniform numbers between 1 and 100
   # Generate some noise
   noise = rnorm(100, mean = 0, sd = 10) # 100 random normal numbers with mean 0 and standard deviation 10
   # Generate dependent variable Y
   Y = 5 + 0.5 * X + noise
   # Combine X and Y into a data frame
   data = data.frame(X, Y)
   

   在这里，我们使用 runif() 函数生成自变量 X，它是一个随机均匀数的向量。然后，我们在因变量 Y 上添加一些“噪声”，使观测数据更真实，且不那么完美地线性。这是通过使用 rnorm() 函数实现的，它创建了一个随机正态数的向量。然后，目标变量 Y 被创建为 X 的函数，加上噪声。

   此外，我们在开始时使用了一个种子（set.seed(123)）来确保可重复性。这意味着每次运行此代码时，我们都会得到相同的一组随机数。如果我们不设置种子，每次运行都会产生不同的随机数列表。

   在这个模拟中，真实的截距（β 0）是 5，真实的斜率（β 1）是 0.5，噪声是均值为 0、标准差为 10 的正态分布。

2. 使用lm()函数基于模拟数据集训练线性回归模型：

   
   # Fit a simple linear regression model
   model = lm(Y ~ X, data = data)
   # Print the model summary
   >>> summary(model)
   Call:
   lm(formula = Y ~ X, data = data)
   Residuals:
        Min       1Q   Median       3Q      Max
   -22.3797  -6.1323  -0.1973   5.9633  22.1723
   Coefficients:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept)  4.91948    1.99064   2.471   0.0152 *
   X            0.49093    0.03453  14.218   <2e-16 ***
   ---
   Signif. codes:
   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 9.693 on 98 degrees of freedom
   Multiple R-squared:  0.6735,  Adjusted R-squared:  0.6702
   F-statistic: 202.2 on 1 and 98 DF,  p-value: < 2.2e-16
   

   在这里，我们使用lm()函数来拟合数据，其中lm代表“线性模型”。此函数创建我们的简单线性回归（SLR）模型。Y ~ X语法是我们指定模型的方式：它告诉函数Y被建模为X的函数。

   summary()函数提供了模型的全面概述，包括估计系数、标准误差、t 值和 p 值等统计量。由于得到的 p 值极低，我们可以得出结论，输入变量具有强烈的统计显著性。

3. 使用plot()和abline()函数可视化数据和拟合的回归线：

   
   # Plot the data
   plot(data$X, data$Y, main = "Simple Linear Regression", xlab = "X", ylab = "Y")
   # Add the fitted regression line
   abline(model, col = "red")
   

   运行此代码生成图 12.2：

图 12.2 – 可视化数据和拟合的回归线

在这里，plot()函数创建了我们数据的散点图，而abline()函数将回归线添加到该图上。这种视觉表示对于理解拟合的质量非常有用。

在下一节中，我们将继续介绍多元线性回归（MLR）模型。

介绍多元线性回归

MLR 将 SLR 中的单个预测因子扩展到基于多个预测变量预测目标结果。在这里，MLR 中的“多个”一词指的是模型中的多个预测因子，其中每个特征都被赋予一个系数。特定的系数β代表在相关预测变量单位变化的情况下，假设所有其他预测变量保持不变，结果变量的变化。

MLR（多元线性回归）的一个显著优点是它能够包含多个预测因子，从而允许对现实世界进行更复杂和真实的（线性）表示。它可以提供对目标变量与所有输入变量之间关系的整体视角。这在结果变量可能受多个预测变量影响的领域中特别有用。它通过以下公式进行建模：

y = β 0 + β 1 x 1 + β 2 x 2 + … + β p x p + ϵ

在这里，我们总共有 p 个特征，因此由于截距项，有(p + 1)个系数。ϵ是表示未解释部分的常规噪声项。换句话说，我们使用 MLR 的预测如下：

ˆ y  = β 0 + β 1 x 1 + β 2 x 2 + … + β p x p

我们可以使用此公式进行 ceteris paribus 分析，这是一种拉丁语表达方式，意思是所有其他事物都相等，我们只改变一个输入变量以评估其对结果变量的影响。换句话说，MLR 允许我们明确控制（即保持不变）许多同时影响目标变量的其他因素，并观察仅一个因素的影响。

例如，假设我们向特征 x_j 添加一个小的增量Δx_j，并保持所有其他特征不变。新的预测ˆy_new 是通过以下公式获得的：

ˆy_new = β_0 + β_1 x_1 + … + β_j( x_j + Δx_j) + … + β_p x_p

我们知道原始预测如下：

ˆy_old = β_0 + β_1 x_1 + … + β_j x_j + … + β_p x_p

这两个之间的差异给我们带来了输出变量的变化：

Δˆy = ˆy_new − ˆy_old = β_j Δx_j

我们在这里做的实质上是控制所有其他输入变量，但只提高 x_j 以观察对预测ˆy 的影响。因此，系数β_j 衡量结果对特定特征的敏感性。当我们有一个单位变化，Δx_j = 1 时，变化正好是系数本身，给我们Δˆy = β_j。

下一节讨论了 MLR 模型的预测性度量。

寻求更高的确定系数

MLR 由于模型中使用的多个预测因子（例如更高的确定系数 R²）而往往比 SLR 表现更好。然而，具有更多输入变量和更高 R²的回归模型并不一定意味着模型拟合得更好，并且可以更好地预测测试集。

由于更多输入特征，更高的 R²可能是因为过拟合。过拟合发生在模型过于复杂，包括过多的预测因子甚至预测因子之间的交互项。在这种情况下，模型可能很好地拟合观察到的数据（从而导致高 R²），但当应用于新的、未见过的测试数据时，它可能表现不佳。这是因为模型可能不仅学会了训练数据的潜在结构，还学会了特定于数据集的随机噪声。

让我们更仔细地看看 R²的度量。虽然 R²衡量模型解释结果变量方差的好坏，但它有一个主要限制：随着更多预测因子的进入，它往往会变得更大，即使这些预测因子是不相关的。作为一种补救措施，我们可以使用调整后的 R²。与 R²不同，调整后的 R²明确考虑了预测因子的数量，并相应地调整了结果统计量。如果一个预测因子显著提高了模型，调整后的 R²将增加，但如果一个预测因子没有通过显著的数量提高模型，调整后的 R²甚至可能减少。

在构建统计模型时，当简单模型与更复杂的模型表现相似时，通常更倾向于选择简单模型。这个简约原则，也称为奥卡姆剃刀，表明在具有相似预测能力的模型中，应该选择最简单的一个。换句话说，向模型中添加更多的预测器会使模型更加复杂，更难以解释，并且更有可能过拟合。

更多关于调整后的 R²

调整后的 R²通过调整所选模型中的特征数量来改进 R²。具体来说，只有当添加这个特征所带来的价值超过添加随机特征所期望的价值时，调整后的 R²的值才会增加。本质上，添加到模型中的额外预测器必须是具有意义和预测性的，以导致调整后的 R²更高。然而，这些额外的预测器在添加到模型中时，总是会提高 R²。

调整后的 R²通过包含模型的自由度来解决这一问题。在这里，自由度指的是在统计计算中可以自由变化的值的数量。在回归模型的背景下，这通常意味着预测器的数量。调整后的 R²可以表示如下：

调整后的 R² = 1 − (1 − R²) × (n − 1) × (n − p − 1)

在这里，n 表示观测值的数量，p 表示模型中的特征数量。

该公式通过根据观测值和预测器的数量调整 R²的尺度来工作。术语(n − 1) × (n − p − 1)是一个反映模型中自由度的比率，其中(n − 1)代表模型中的总自由度。我们减去 1 是因为我们是从数据中估计因变量的均值。而(n − p − 1)代表误差的自由度，它反过来代表在估计模型参数后剩余的观测值数量。整个术语(1 − R²) × (n − 1) × (n − p − 1)表示调整了预测器数量的误差方差。

从 1 中减去这个值，得到的是模型解释的总方差的比例，这是在调整了模型中预测器的数量之后的结果。换句话说，它是一种 R²的版本，惩罚了不必要的预测器的添加。这有助于防止过拟合，使得调整后的 R²在比较具有不同预测器数量的模型时，成为一个更平衡的解释力度量。

让我们看看如何在 R 中开发一个多元线性回归（MLR）模型。

开发一个 MLR 模型

在本节中，我们将使用 R 中的相同lm()函数，基于预加载的mtcars数据集来开发一个 MLR 模型，该数据集在之前的练习中使用过。再次强调，mtcars数据集包含了 1974 年《汽车趋势》杂志中 32 辆汽车的测量数据。这些测量包括诸如每加仑英里数（mpg）、汽缸数（cyl）、马力（hp）和重量（wt）等属性。

练习 12.2 – 构建一个 MLR 模型

在这个练习中，我们将开发一个 MLR 模型，使用cyl、hp和wt预测mpg。然后我们将解释模型结果：

1. 加载mtcars数据集并构建一个基于cyl、hp和wt预测mpg的 MLR 模型：

   
   # Load the data
   data(mtcars)
   # Build the model
   model = lm(mpg ~ cyl + hp + wt, data = mtcars)
   

   在这里，我们首先使用data()函数加载mtcars数据集，然后使用lm()函数构建 MLR 模型。mpg ~ cyl + hp + wt公式用于指定模型。这个公式告诉 R，我们想要将mpg建模为cyl、hp和wt的函数。data = mtcars参数告诉 R 在mtcars数据集中寻找这些变量。lm()函数将模型拟合到数据，并返回一个模型对象，我们将其存储在变量 model 中。

2. 查看模型的摘要：

   
   # Print the summary of the model
   >>> summary(model)
   Call:
   lm(formula = mpg ~ cyl + hp + wt, data = mtcars)
   Residuals:
       Min      1Q  Median      3Q     Max
   -3.9290 -1.5598 -0.5311  1.1850  5.8986
   Coefficients:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
   cyl         -0.94162    0.55092  -1.709 0.098480 .
   hp          -0.01804    0.01188  -1.519 0.140015
   wt          -3.16697    0.74058  -4.276 0.000199 ***
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 2.512 on 28 degrees of freedom
   Multiple R-squared:  0.8431,  Adjusted R-squared:  0.8263
   F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11
   

   摘要包括模型的系数（每个预测变量的截距和斜率）、残差（实际观察值与目标预测值之间的差异），以及一些统计量，告诉我们模型如何拟合数据，包括 R²和调整 R²。

   让我们解释输出结果。每个系数代表在相关预测变量增加一个单位时，mpg的预期变化，假设其他所有预测变量保持不变。R²值，即0.8431，表示mpg中可由预测变量共同解释的方差比例（超过 84%）。再次强调，调整 R²值，即0.8263，是一个经过修改的 R²，它考虑了模型中的特征数量。

   此外，每个预测变量的 p 值测试了系数真实值为零的零假设。如果一个预测变量的 p 值小于预设的显著性水平（如 0.05），我们会拒绝这个零假设，并得出结论说该预测变量具有统计学意义。在这种情况下，与使用 5%显著性水平比较的其他因素相比，wt是唯一具有统计学意义的因素。

在 MLR 模型中，所有系数都是负数，表明输入变量和目标变量之间存在反向的旅行方向。然而，我们不能得出所有预测变量都与目标变量负相关的结论。在 SLR 中，单个预测变量与目标变量的相关性可能是正的或负的。

下一节将提供更多关于这一现象的背景信息。

介绍辛普森悖论

辛普森悖论指出，趋势在不同数据组中显现，但在合并时消失或改变。在回归分析的情况下，当控制其他变量时，一个看似与结果正相关的变量可能具有负相关性，辛普森悖论就可能出现。

实质上，这个悖论说明了考虑混杂变量的重要性，以及在未理解背景的情况下，不要从汇总数据中得出结论。混杂变量是那些不在考虑的解释变量中，但与目标变量和预测变量都相关的变量。

让我们通过以下练习来考虑一个简单的例子。

练习 12.3 – 说明辛普森悖论

在这个练习中，我们将观察两个场景，在 SLR 和 MLR 中，相同特征的系数值符号相反：

1. 创建一个包含两个预测变量和一个输出变量的虚拟数据集：

   
   set.seed(123)
   x1 = rnorm(100)
   x2 = -3 * x1 + rnorm(100)
   y = 2 + x1 + x2 + rnorm(100)
   df = data.frame(y = y, x1 = x1, x2 = x2)
   

   在这里，x1是一组从标准正态分布中随机生成的 100 个数字。x2是x1的线性函数，但具有负相关性，并添加了一些随机噪声（通过rnorm(100)）。然后，y作为x1和x2的线性函数生成，同样添加了一些随机噪声。所有三个变量都编译到一个 DataFrame，df中。

2. 使用y作为结果变量和x1作为输入特征来训练一个简单线性回归（SLR）模型。检查模型的摘要：

   
   # Single linear regression
   single_reg = lm(y ~ x1, data = df)
   >>> summary(single_reg)
   Call:
   lm(formula = y ~ x1, data = df)
   Residuals:
       Min      1Q  Median      3Q     Max
   -2.7595 -0.8365 -0.0564  0.8597  4.3211
   Coefficients:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept)   2.0298     0.1379   14.72   <2e-16 ***
   x1           -2.1869     0.1511  -14.47   <2e-16 ***
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 1.372 on 98 degrees of freedom
   Multiple R-squared:  0.6813,  Adjusted R-squared:  0.678
   F-statistic: 209.5 on 1 and 98 DF,  p-value: < 2.2e-16
   

   结果显示，由于系数为-2.1869，x1与y呈负相关。

3. 使用y作为目标变量，x1和x2作为输入特征来训练一个 SLR 模型。检查模型的摘要：

   
   # Multiple linear regression
   multi_reg = lm(y ~ x1 + x2, data = df)
   >>> summary(multi_reg)
   Call:
   lm(formula = y ~ x1 + x2, data = df)
   Residuals:
       Min      1Q  Median      3Q     Max
   -1.8730 -0.6607 -0.1245  0.6214  2.0798
   Coefficients:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept)  2.13507    0.09614  22.208  < 2e-16 ***
   x1           0.93826    0.31982   2.934  0.00418 **
   x2           1.02381    0.09899  10.342  < 2e-16 ***
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 0.9513 on 97 degrees of freedom
   Multiple R-squared:  0.8484,  Adjusted R-squared:  0.8453
   F-statistic: 271.4 on 2 and 97 DF,  p-value: < 2.2e-16
   

   结果显示，x1的估计系数现在是一个正数。这表明x1突然与y呈正相关吗？不，因为可能存在其他导致系数为正的混杂变量。

关键的启示是，我们只能在 SLR 设置中（正或负）对预测变量和目标结果之间的相关性进行推断。例如，如果我们构建一个 SLR 模型来对y进行回归，如果得到的系数是正的（β > 0），我们可以得出结论，x和y是正相关的。同样，如果β > 0，我们可以得出结论，x和y是正相关的。同样适用于负相关的情况。

然而，在多元线性回归（MLR）设置中，这样的推断会失效 – 也就是说，如果我们有β > 0，我们不能得出正相关的结论，反之亦然。

让我们利用这个机会来解释结果。Estimate列显示了估计的回归系数。这些值表示当相应的预测变量增加一个单位，同时保持所有其他特征不变时，y变量预期会增加多少。在这种情况下，对于x1每增加一个单位，y预计会增加大约0.93826个单位，对于x2每增加一个单位，y预计会增加大约1.02381个单位。(Intercept)行显示了当模型中的所有预测变量都为零时y的估计值。在这个模型中，估计的截距是2.13507。

标准误差表示估计的标准误差。这里的值越小，表示估计越精确。t 值列显示了假设检验的 t 统计量，即给定其他预测变量在模型中，对应的总体回归系数为零。t 统计量绝对值越大，表明对零假设的证据越强。Pr(>|t|)列给出了假设检验的 p 值。在这种情况下，x1 和 x2 的 p 值都低于 0.05，表明在 5% 的显著性水平下，两者都是 y 的统计显著预测变量。

最后，多重 R 平方和调整 R 平方值提供了模型拟合数据的度量。多重 R 平方值为 0.8484，表明该模型解释了 y 的变异性约 84.84%。调整 R 平方值对此度量进行了模型中特征数量的调整。如前所述，当比较具有不同数量预测变量的模型时，这是一个更好的度量。这里的调整 R 平方值为 0.8453。F 统计量及其相关的 p 值用于检验所有总体回归系数为零的假设。小的 p 值（小于 0.05）表明我们可以拒绝这个假设，并得出结论，至少有一个预测变量在预测 y 时是有用的。

下一个部分将探讨在 MLR 模型中有一个分类预测变量的情况。

与分类变量一起工作

在多元线性回归（MLR）中，包含一个二元预测变量的过程与包含一个数值预测变量的过程类似。然而，解释不同。考虑一个数据集，其中 y 是目标变量，x 1 是数值预测变量，x 2 是二元预测变量：

y = β 0 + β 1 x 1 + β 2 x 2 + ϵ

在这个模型中，x 2 被编码为 0 和 1，其对应的系数 β 2 代表了由 x 2 识别的两个组之间 y 的均值差异。

例如，如果 x 2 是一个表示性别的二元变量（男性为 0，女性为 1），y 是工资，那么 β 2 代表在考虑 x 1 的值后，女性和男性之间的平均工资差异。

注意，二元预测变量的系数解释依赖于模型中的其他变量。因此，在前面的例子中，β 2 是在给定 x 1 的值时，女性和男性之间的工资差异。

在实现方面，当在回归模型中使用因子时，R 会自动创建虚拟变量。因此，如果 x 2 是一个具有“男性”和“女性”级别的因子，R 将在拟合模型时内部处理将其转换为 0 和 1。

让我们看看一个具体的例子。在下面的代码中，我们正在构建一个 MLR 模型，使用 qsec 和 am 来预测 mpg：

# Fit the model
model <- lm(mpg ~ qsec + am, data = mtcars)
# Display the summary of the model
>>> summary(model)
Call:
lm(formula = mpg ~ qsec + am, data = mtcars)
Residuals:
    Min      1Q  Median      3Q     Max
-6.3447 -2.7699  0.2938  2.0947  6.9194
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -18.8893     6.5970  -2.863  0.00771 **
qsec          1.9819     0.3601   5.503 6.27e-06 ***
am            8.8763     1.2897   6.883 1.46e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.487 on 29 degrees of freedom
Multiple R-squared:  0.6868,  Adjusted R-squared:  0.6652
F-statistic:  31.8 on 2 and 29 DF,  p-value: 4.882e-08

注意，am 变量被当作数值变量处理。因为它代表汽车中的传动类型（0 = 自动，1 = 手动），它本应被当作分类变量处理。这可以通过将其转换为因子来实现，如下所示：

# Convert am to categorical var
mtcars$am_cat = as.factor(mtcars$am)
# Fit the model
model <- lm(mpg ~ qsec + am_cat, data = mtcars)
# Display the summary of the model
>>> summary(model)
Call:
lm(formula = mpg ~ qsec + am_cat, data = mtcars)
Residuals:
    Min      1Q  Median      3Q     Max
-6.3447 -2.7699  0.2938  2.0947  6.9194
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -18.8893     6.5970  -2.863  0.00771 **
qsec          1.9819     0.3601   5.503 6.27e-06 ***
am_cat1       8.8763     1.2897   6.883 1.46e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.487 on 29 degrees of freedom
Multiple R-squared:  0.6868,  Adjusted R-squared:  0.6652
F-statistic:  31.8 on 2 and 29 DF,  p-value: 4.882e-08

注意，只为分类变量 am_cat 创建了一个变量 am_cat1。这是因为 am_cat 是二元的，因此我们只需要一个虚拟列（在这种情况下保留 am_cat1 并删除 am_cat0）来表示原始分类变量。一般来说，对于一个有 k 个分类的分类变量，R 将在模型中自动创建 (k − 1) 个虚拟变量。

这个过程被称为 am 等于相应的水平，否则为 0。这实际上创建了一组指标，捕捉每个类别的存在或不存在。最后，由于我们可以根据前一个（k-1）虚拟变量的值推断最后一个虚拟变量，因此我们可以从结果中删除最后一个虚拟变量。

使用分类变量会给模型估计引入一个垂直位移，如以下章节所述。为了看到这一点，让我们更仔细地看看分类变量 am_cat1 的影响。我们的 MLR 模型现在假设以下形式：

ˆ y  = β 0 + β 1 x qsec + β 2 x am_cat

我们知道 x am_cat 是一个二元变量。当 x am_cat = 0 时，预测如下：

ˆ y  = β 0 + β 1 x qsec

当 x am_cat = 1 时，预测如下：

ˆ y  = β 0 + β 1 x qsec + β 2 = (β 0 + β 2) + β 1 x qsec

通过比较这两个量，我们可以看出它们是两个相互平行的线性模型，因为斜率相同，唯一的区别是截距项中的 β 2。

一个视觉说明在这里很有帮助。在下面的代码片段中，我们首先创建一个新的 DataFrame，newdata，它覆盖了原始数据中 qsec 值的范围，对于每个 am_cat 值（0 和 1）。然后，我们使用 predict() 函数从模型中获取新数据的预测 mpg 值。接下来，我们使用 geom_point() 绘制原始数据点，并使用 geom_line() 添加两条回归线，其中线基于 newdata 中的预测值。color = am_cat 美学设置为不同的 am_cat 值添加不同的颜色，并在 scale_color_discrete() 中调整标签，以便 0 对应于“自动”，1 对应于“手动”：

# Load required library
library(ggplot2)
# Create new data frame for the predictions
newdata = data.frame(qsec = seq(min(mtcars$qsec), max(mtcars$qsec), length.out = 100),
                      am_cat = c(rep(0, 100), rep(1, 100)))
newdata$am_cat = as.factor(newdata$am_cat)
# Get predictions
newdata$mpg_pred = predict(model, newdata)
# Plot the data and the regression lines
ggplot(data = mtcars, aes(x = qsec, y = mpg, color = am_cat)) +
  geom_point() +
  geom_line(data = newdata, aes(y = mpg_pred)) +
  labs(title = "mpg vs qsec by Transmission Type",
       x = "Quarter Mile Time (qsec)",
       y = "Miles per Gallon (mpg)",
       color = "Transmission Type") +
  scale_color_discrete(labels = c("Automatic", "Manual")) +
  theme(text = element_text(size = 16),  # Default text size
        title = element_text(size = 15),  # Title size
        axis.title = element_text(size = 18),  # Axis title size
        legend.title = element_text(size = 16),  # Legend title size
        legend.text = element_text(size = 16), # Legend text size
        legend.position = "bottom")  # Legend position

运行此代码生成 图 12*.3*：

图 12.3 – 基于不同传动类型的两个线性回归模型的可视化。由于截距项的位移，这两条线是平行的

这个图表明，在相同的四分之一英里时间（qsec）下，手动变速汽车比自动变速汽车的每加仑英里数（mpg）相同。然而，在实际情况中，由于不同类型的汽车（手动与自动）可能具有不同的四分之一英里时间，这不太可能。换句话说，这两个变量之间存在交互作用。

下一个部分将介绍交互项作为这种情况的补救措施。

引入交互项

在回归分析中，当一个预测变量的影响因另一个预测变量的水平而不同时，就会发生交互作用。在我们的运行示例中，我们实际上是在查看mpg和qsec之间的关系是否因am的不同值而不同。换句话说，我们正在测试将mpg和qsec联系起来的线的斜率是否因手动（am=1）和自动（am=0）变速而不同。

例如，如果没有交互作用，那么qsec对mpg的影响是相同的，无论汽车是手动还是自动。这意味着描述手动和自动汽车之间mpg与qsec关系的线条将是平行的。

如果存在交互作用，那么qsec对mpg的影响在手动和自动汽车之间是不同的。这意味着描述手动和自动汽车之间mpg与qsec关系的线条不会平行。它们可能相交，或者更常见的是，具有不同的斜率。

为了描绘这些关系中的差异，我们可以在模型中添加一个交互项，这可以通过*运算符来完成。例如，具有qsec和am_cat之间交互的回归模型的公式将是mpg ~ qsec * am_cat。这相当于mpg ~ qsec + am``_cat + qsec:am_cat，其中qsec:am_cat代表交互项。以下代码显示了详细信息：

# Adding interaction term
model_interaction <- lm(mpg ~ qsec * am_cat, data = mtcars)
# Print model summary
>>> summary(model_interaction)
Call:
lm(formula = mpg ~ qsec * am_cat, data = mtcars)
Residuals:
    Min      1Q  Median      3Q     Max
-6.4551 -1.4331  0.1918  2.2493  7.2773
Coefficient s:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   -9.0099     8.2179  -1.096  0.28226
qsec           1.4385     0.4500   3.197  0.00343 **
am_cat1      -14.5107    12.4812  -1.163  0.25481
qsec:am_cat1   1.3214     0.7017   1.883  0.07012 .
Signif. code s:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.343 on 28 degrees of freedom
Multiple R-squared:  0.722,  Adjusted R-squared:  0.6923
F-statistic: 24.24 on 3 and 28 DF,  p-value: 6.129e-08

让我们再绘制一个更新后的模型，该模型由由于交互作用而产生的两条相交线组成。在以下代码片段中，使用geom_smooth(method =""l"", se = FALSE)为每组点（自动和手动汽车）拟合不同的线性线。as.factor(am_cat)用于将am_cat视为一个因子（分类）变量，以便为每个类别拟合一条单独的线：

# Create scatter plot with two intersecting lines
ggplot(mtcars, aes(x = qsec, y = mpg, color = as.factor(am_cat))) +
  geom_point() +
  geom_smooth(method =""l"", se = FALSE) + # fit separate regression lines per group
  scale_color_discrete(name =""Transmission Typ"",
                       labels = c""Automati"",""Manua"")) +
  labs(x =""Quarter mile time (seconds"",
       y =""Miles per gallo"",
       title =""Separate regression lines fit for automatic and manual car"") +
  theme(text = element_text(size = 16),
        title = element_text(size = 15),
        axis.title = element_text(size = 20),
        legend.title = element_text(size = 16),
        legend.text = element_text(size = 16))

运行此代码生成图 12.4：

图 12.4 – 由于四分之一英里时间和变速类型之间的交叉项而产生的两条相交线

下一个部分将关注另一个相关主题：非线性项。

处理非线性项

线性回归是理解响应变量和解释变量之间线性关系的一个广泛使用的模型。然而，数据中的所有潜在关系并不都是线性的。在许多情况下，特征和响应变量可能没有直线关系，因此需要在线性回归模型中处理非线性项以增加其灵活性。

将非线性引入的最简单方法，即构建曲线而不是直线，是在回归模型中包含预测变量的多项式项。在多项式回归中，一些或所有预测变量被提升到特定的多项式项 – 例如，将特征 x 转换为 x²或 x³。

让我们通过一个练习来了解向线性回归模型添加多项式特征的影响。

练习 12.4 – 向线性回归模型添加多项式特征

在这个练习中，我们将创建一个简单的数据集，其中输入特征 x 和目标变量 y 之间的关系不是线性的，而是二次的。首先，我们将拟合一个简单线性回归（SLR）模型，然后添加二次项并比较结果：

1. 生成一个从-10 到 10 的x值序列。对于每个x，计算相应的y作为x的平方，并添加一些随机噪声以显示（有噪声的）二次关系。将x和y放入 DataFrame 中：

   
   # Create a quadratic dataset
   set.seed(1)
   x = seq(-10, 10, by = 0.5)
   y = x² + rnorm(length(x), sd = 5)
   # Put it in a dataframe
   df = data.frame(x = x, y = y)
   Fit a simple linear regression to the data and print the summary.
   lm1 <- lm(y ~ x, data = df)
   >>> summary(lm1)
   Call:
   lm(formula = y ~ x, data = df)
   Residuals:
       Min      1Q  Median      3Q     Max
   -43.060 -29.350  -5.451  19.075  64.187
   Coefficient s:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept) 35.42884    5.04627   7.021 2.01e-08 ***
   x           -0.04389    0.85298  -0.051    0.959
   Signif. code s:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 32.31 on 39 degrees of freedom
   Multiple R-squared:  6.787e-05,  Adjusted R-squared:  -0.02557
   F-statistic: 0.002647 on 1 and 39 DF,  p-value: 0.9592
   

   结果表明，模型对数据的拟合度不佳，数据具有非线性关系。

2. 使用I()函数将 x²作为预测变量包括在内，对数据进行二次模型拟合。打印模型的摘要：

   
   lm2 <- lm(y ~ x + I(x²), data = df)
   >>> summary(lm2)
   Call:
   lm(formula = y ~ x + I(x²), data = df)
   Residuals:
       Min      1Q  Median      3Q     Max
   -11.700  -2.134  -0.078   2.992   7.247
   Coefficients:
               Estimate Std. Error t value Pr(>|t|)
   (Intercept)  0.49663    1.05176   0.472    0.639
   x           -0.04389    0.11846  -0.370    0.713
   I(x²)       0.99806    0.02241  44.543   <2e-16 ***
   Signif. code s:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   Residual standard error: 4.487 on 38 degrees of freedom
   Multiple R-squared:  0.9812,  Adjusted R-squared:  0.9802
   F-statistic: 992.1 on 2 and 38 DF,  p-value: < 2.2e-16
   

   结果显示，多项式模型比简单线性模型更好地拟合数据。因此，当预测变量和响应变量之间的关系不是严格线性时，添加非线性项可以提高模型拟合度。

3. 将线性模型和二次模型与数据一起绘制：

   
   ggplot(df, aes(x = x, y = y)) +
     geom_point() +
     geom_line(aes(y = linear_pred), color =""blu"", linetype =""dashe"") +
     geom_line(aes(y = quadratic_pred), color =""re"") +
     labs(title =""Scatter plot with linear and quadratic fit"",
          x ="""",
          y ="""") +
     theme(text = element_text(size = 15)) +
     scale_color_discrete(name =""Mode"",
                          labels = c"«Linear Mode"»,"«Quadratic Mode"»)) +
     annotate""tex"", x = 0, y = 40, label =""Linear Mode"", color =""blu"") +
     annotate""tex"", x = 6, y = 80, label =""Quadratic Mode"", color =""re"")
   

   在这里，我们首先计算两个模型的预测值并将它们添加到 DataFrame 中。然后，我们创建一个散点图来表示数据，并添加两条线来表示线性模型（蓝色）和二次模型（红色）的预测值。

   运行此代码生成图 12.5。结果表明，添加多项式特征可以扩展线性模型的灵活性：

图 12.5 – 非线性数据的线性与二次拟合可视化

引入非线性的一些其他常见方法包括对数变换（logx）或平方根变换（√x）。这些变换也可以应用于目标变量 y，我们可以在同一个线性模型中有多个转换特征。

注意，具有转换特征的模型仍然是一个线性模型。如果系数有非线性转换，则模型将是非线性模型。

下一节将更详细地介绍一种广泛使用的变换类型：对数变换。

更多关于对数变换的内容

对数变换，或称对数转换，基于对数函数将输入映射到相应的输出，即 y = logx。使用这种变换的一个流行原因是向线性回归模型中引入非线性。当输入特征与目标输出之间的关系是非线性时，应用变换有时可以使关系线性化，从而可以使用线性回归模型来建模这种关系。对于对数变换，当预测变量的值增加时，结果变量的变化率增加或减少时，它可以帮助实现这一点。

具体来说，当输入变得更为极端时，变化率会降低。这种变换的自然结果是，输入数据中的潜在异常值被压缩，使得它们在变换后的列中看起来不那么极端。换句话说，由此产生的线性回归模型将由于对数变换而对原始异常值不那么敏感。

使用对数变换的另一个副作用是它能够处理异方差性。异方差性是指回归模型中误差项的变异性在预测变量的所有水平上并不恒定。这违反了线性回归模型的一个假设，可能导致效率低下和有偏的估计。在这种情况下，通过对数变换可以缩小潜在的大的误差项，使得误差项的方差在不同预测变量的水平上更加恒定。

最后，当预测变量与结果之间的关系是乘法而非加法时，对预测变量和/或结果变量取对数可以将关系转换为加法关系，这可以使用线性回归进行建模。

让我们考虑一个例子，其中我们预测每加仑英里数（mpg）从马力（hp）。我们将比较直接从hp预测mpg的模型，以及从hp预测mpg的对数的模型，如下面的代码片段所示：

# Fit the original model
model_original = lm(mpg ~ hp, data = mtcars)
# Fit the log-transformed model
mtcars$log_mpg = log(mtcars$mpg)
model_log = lm(log_mpg ~ hp, data = mtcars)
# Predictions from the original model
mtcars$pred_original = predict(model_original, newdata = mtcars)
# Predictions from the log-transformed model (back-transformed to the original scale using exp)
mtcars$pred_log = exp(predict(model_log, newdata = mtcars))
library(tidyr)
library(dplyr)
# Reshape data to long format
df_long <- mtcars %>%
  gather(key =""Mode"", value =""Predictio"", pred_original, pred_log)
# Create plot
ggplot(df_long, aes(x = hp, y = mpg)) +
  geom_point(data = mtcars, aes(x = hp, y = mpg)) +
  geom_line(aes(y = Prediction, color = Model)) +
  labs(
    x =""Horsepower (hp"",
    y =""Miles per gallon (mpg"",
    color =""Mode""
  ) +
  scale_color_manual(values = c""pred_origina"" =""blu"",""pred_lo"" =""re"")) +
  theme(
    legend.position =""botto"",
    text = element_text(size = 16),
    legend.title = element_text(size = 16),
    axis.text = element_text(size = 16),  # control the font size of axis labels
    legend.text = element_text(size = 16)  # control the font size of legend text
  )

运行此代码生成图 12.6，其中我们可以看到蓝色线有轻微的弯曲：

图 12.6 – 可视化原始和对数变换后的模型

注意，对数变换只能应用于正数数据。在mtcars$mpg的情况下，所有值都是正数，因此我们可以安全地应用对数变换。如果变量包含零或负值，我们需要考虑不同的变换或方法。

下一节将重点介绍推导和使用线性回归模型的闭式解。

使用闭式解

在开发线性回归模型时，已知的训练集（X, y）被给出，唯一未知的参数是系数，β。在这里，粗体小写字母表示向量（例如 β 和 y），而粗体大写字母表示矩阵（例如 X）。结果证明，线性回归模型的闭式解可以通过使用普通最小二乘法（OLS）估计的概念来推导，该估计法旨在最小化模型中残差的平方和。拥有闭式解意味着我们可以简单地插入所需元素（在这种情况下，X 和 y）并执行计算以获得解，而无需求助于任何优化过程。

具体来说，给定一个数据矩阵，X（其中包含一个用于截距项的列，并用粗体表示表示有多个特征），其维度为 n × p（其中 n 是观测数，p 是预测数）和一个长度为 n 的响应向量，y，OLS 对系数向量 β 的估计由以下公式给出：

β = (X T X) −1 X T y

这个解假设项 (X T X) 是可逆的，这意味着它应该是一个满秩矩阵。如果不是这种情况，解可能不存在或不唯一。

现在，让我们看看如何推导这个解。我们从最小二乘法的最小化问题开始：最小化 (y − Xβ) T(y − Xβ) 关于 β。这个二次型可以展开为 y T y − β T X T y − y T Xβ + β T X T Xβ。注意 β T X T y = y T Xβ，因为这两个项都是标量，因此在转置操作后它们相等。我们可以将残差平方和（RSS）表达式写为 y T y − 2 β T X T y + β T X T Xβ。

在这里，我们应用一阶条件来求解使此表达式最小化的 β 值（回想一下，图形上最小化或最大化的点具有导数为 0）。这意味着我们会将其一阶导数设为零，从而得到以下公式：

∂ ( y T y − 2 β T X T y + β T X T Xβ)  ____________________  ∂ β  = − 2 X T y + 2 X T Xβ = 0

X T Xβ = X T y

β = (X T X) −1 X T y

因此，我们已经推导出 β 的闭式解，该解最小化了残差平方和。让我们通过一个示例来看看它是如何实现的。

实现闭式解

让我们看看如何在 R 中实现线性回归模型的 OLS 估计。以下代码片段展示了使用合成数据的示例：

# Set seed for reproducibility
set.seed(123)
# Generate synthetic data
n = 100 # number of observations
x = runif(n, -10, 10) # predictors
beta0 = 2 # intercept
beta1 = 3 # slope
epsilon = rnorm(n, 0, 2) # random error term
y = beta0 + beta1*x + epsilon # response variable
# Design matrix X
X = cbind(1, x)

在这里，我们生成具有单个输入特征的 100 个观测值，观测值受到噪声扰动，并遵循 y = β 0 + β 1 x + ϵ 的过程。误差项假设为具有均值为 0 和标准差为 2 的正态分布。

在进行估计之前，请注意，我们在输入特征x的左侧也添加了一个 1 的列，以形成一个矩阵X。这个 1 的列用于表示截距项，通常被称为偏差技巧。也就是说，截距项的系数β0 将是系数向量的一部分，因此不需要为截距创建一个单独的系数。

让我们使用闭式解来计算结果：

beta_hat = solve(t(X) %*% X) %*% t(X) %*% y
>>> print(beta_hat)
       [,1]
  1.985344
x 3.053152

在这里，%*%用于矩阵乘法，t(X)是X的转置，solve()用于计算矩阵的逆。

我们也可以使用lm()函数运行线性回归过程进行比较：

# Fit linear regression model for comparison
model = lm(y ~ x)
>>> print(coef(model))
(Intercept)           x
   1.985344    3.053152

结果与通过手动方法获得的结果相同。

下两节涵盖了线性回归设置中的两个常见问题：共线性和非同方差性。

处理共线性

共线性指的是在多元回归模型中，两个（或更多）预测变量高度相关的情况。这意味着一个自变量可以以很高的精度从其他变量线性预测出来。这是我们不想陷入的情况。换句话说，我们希望看到预测变量与目标变量之间有高度的相关性，而预测变量之间本身的相关性较低。

面对线性回归模型中的共线性，得到的模型往往会生成不可靠和不稳定的回归系数估计。它可能会膨胀参数的系数，使它们在统计上不显著，尽管它们可能具有实质性的重要性。此外，共线性使得难以评估每个自变量对因变量的影响，因为影响是交织在一起的。然而，它不会影响模型的预测能力或可解释性；相反，它只改变了单个特征的计算。

检测预测变量之间是否存在潜在的共线性可以通过检查成对的相关性来完成。或者，我们可以求助于一个特定的统计量，称为方差膨胀因子（VIF），它量化了由于共线性导致方差增加的程度。VIF 值为 1 表示两个变量不相关，而 VIF 值大于 5（在许多领域）则表明存在问题的共线性程度。

当线性回归模型中存在共线性时，我们可以选择只保留一个预测变量，并移除所有其他高度相关的预测变量。我们还可以通过主成分分析（PCA），这是一种广泛用于降维的技术，将这些相关变量组合成几个不相关的变量。此外，我们可以求助于岭回归来控制系数的大小；这将在本章后面介绍。

要使用 VIF 检查多重共线性，我们可以使用car包中的vif()函数，如下面的代码片段所示：

# install the package if not already installed
if(!require(car)) install.packages('car')
# load the package
library(car)
# fit a linear model
model = lm(mpg ~ hp + wt + disp, data = mtcars)
# calculate VIF
vif_values = vif(model)
>>> print(vif_values)
      hp       wt     disp
2.736633 4.844618 7.324517

从结果来看，disp似乎具有高度的多重共线性（VIF = 7.32 > 5），这表明它与hp和wt有很强的相关性。这意味着disp并没有提供其他两个预测变量中已经包含的信息。

为了处理这里的多重共线性问题，我们可以考虑从模型中移除disp，因为它具有最高的 VIF 值，应用 PCA 来合并三个预测变量，或者使用岭回归或 Lasso 回归（关于这一点，请参阅本章的最后两部分）。

下一节将重点讨论异方差性的问题。

处理异方差性

异方差性（或异方差性）指的是误差项或残差的变异性在不同独立变量的所有水平上并不相同的情况。这违反了 OLS 线性回归的一个关键假设，即假设残差具有恒定的方差——换句话说，残差是同方差的。违反这个假设可能导致对系数的统计显著性推断不正确，因为回归系数的估计标准误差可能比应有的更大或更小。

处理异方差性的方法有几种。首先，我们可以使用之前介绍的对数函数对结果变量进行转换。其他函数，如取原始结果变量的平方根或倒数，也可能有助于减少异方差性。探索高级回归模型，如加权最小二乘法（WLS）或广义最小二乘法（GLS），也可能有助于减少异方差性的影响。

为了正式检验异方差性，我们可以使用lmtest包中的bptest()函数进行 Breusch-Pagan 检验。在下面的代码片段中，我们拟合了一个 MLR 模型，使用wt和hp预测mpg，然后执行 Breusch-Pagan 检验：

# Load library
library(lmtest)
# Fit a simple linear regression model on mtcars dataset
model = lm(mpg ~ wt + hp, data = mtcars)
# Perform a Breusch-Pagan test to formally check for heteroskedasticity
>>> bptest(model)
  studentized Breusch-Pagan test
data:  model
BP = 0.88072, df = 2, p-value = 0.6438

由于 p 值（0.6438）大于 0.05，我们不能拒绝 Breusch-Pagan 检验的零假设。这表明没有足够的证据表明回归模型中存在异方差性。因此，我们可以得出结论，残差的方差与恒定的方差没有显著差异，或者说，残差是同方差的。

下一节将转向查看正则化线性回归模型和岭回归和 Lasso 惩罚。

引入惩罚线性回归

惩罚回归模型，如岭回归和 Lasso，是用于处理多重共线性、减少过拟合甚至进行变量选择的技术，尤其是在处理具有多个输入特征的高维数据时。

岭回归（也称为 L2 正则化）是一种方法，它添加了一个与系数幅度的平方相当的惩罚。我们会在加权一个额外的超参数后，通常表示为 λ，将其添加到损失函数中，以控制惩罚项的强度。

另一方面，Lasso 回归（L1 正则化）是一种方法，类似于岭回归，它对非零系数添加惩罚，但与岭回归不同，当惩罚调整参数足够大时，它可以强制某些系数精确等于零。超参数 λ 的值越大，收缩量就越大。对系数大小的惩罚有助于减少模型复杂性和多重共线性，从而使得模型在未见数据上具有更好的泛化能力。然而，岭回归包括最终模型中的所有特征，因此不会产生任何稀疏性。因此，它对变量选择并不特别有用。

总结来说，岭回归和 Lasso 回归都是惩罚线性回归方法，它们在模型优化过程中添加了对估计系数幅度的约束，这有助于防止过度拟合、管理多重共线性并减少模型复杂度。然而，岭回归倾向于将所有预测变量包含在模型中并帮助减少它们的影响，而 Lasso 可以完全排除模型中的预测变量，从而得到更简单、更可解释的模型。

让我们从岭回归开始，更仔细地看看它的损失函数。

使用岭回归

岭回归，也称为 L2 正则化，是一种常用的技术，通过惩罚线性回归模型中估计系数的幅度来减轻过度拟合。

回想一下，在简单线性回归（SLR）模型中，我们试图最小化预测值和实际值之间平方差的和，这被称为最小二乘法。我们希望最小化的损失函数是均方误差（RSS）：

RSS = ∑_i=1^n (y_i - (β_0 + ∑_j=1^p β_j x_ij))²

在这里，y_i 是实际的目标值，β_0 是截距项，{β_j} 是每个预测变量 x_ij 的系数估计，求和是总的观测值和预测变量。

仅最小化 RSS 会导致过度拟合的模型，如结果系数的高幅度所示。作为补救措施，我们可以通过向这个损失函数添加惩罚项来应用岭回归。这个惩罚项是每个系数的平方和乘以一个调整参数 λ。因此，岭回归损失函数（也称为成本函数）如下：

L_ridge = RSS + λ∑_j=1^p β_j²

在这里，λ参数是一个用户定义的调整参数。较大的λ意味着更高的惩罚，较小的λ意味着更少的正则化效果。λ = 0 给出普通最小二乘回归结果，而当λ接近无穷大时，惩罚项的影响占主导地位，系数估计接近零。

通过添加这个惩罚项，岭回归倾向于减小系数的大小，这有助于缓解多重共线性问题（预测变量高度相关）。它是通过将相关预测变量的系数估计分散到彼此之间来实现的，这可能导致一个更稳定且可解释的模型。

然而，需要注意的是，岭回归通常不会产生稀疏解，并且不执行变量选择。换句话说，它不会导致某些系数恰好为零的模型（除非λ是无穷大），因此所有预测变量都包含在模型中。如果特征选择很重要，那么 lasso（L1 正则化）或弹性网络（L1 和 L2 正则化的组合）等方法可能更合适。

注意，我们通常会惩罚截距，β 0。

让我们通过一个练习来学习如何开发岭回归模型。

练习 12.5 – 实现岭回归

在这个练习中，我们将实现一个岭回归模型，并将估计系数与 OLS 模型进行比较。我们的实现将基于glmnet包：

1. 安装并加载glmnet包：

   
   # install the package if not already installed
   if(!require(glmnet)) install.packages('glmnet')
   library(glmnet)
   

   在这里，我们使用if-else语句来检测glmnet包是否已安装。

2. 将除mpg之外的所有列作为预测变量存储在X中，将mpg作为目标变量存储在y中：

   
   # Prepare data
   data(mtcars)
   X = as.matrix(mtcars[, -1]) # predictors
   y = mtcars[, 1] # response
   

3. 使用glmnet()函数拟合岭回归模型：

   
   # Fit ridge regression model
   set.seed(123) # for reproducibility
   ridge_model = glmnet(X, y, alpha = 0)
   

   在这里，alpha 参数控制着我们拟合的模型类型。alpha = 0 拟合岭回归模型，alpha = 1 拟合 lasso 模型，而介于两者之间的任何值都会拟合弹性网络模型。

4. 使用交叉验证选择最佳的lambda值：

   
   # Use cross-validation to find the optimal lambda
   cv_ridge = cv.glmnet(X, y, alpha = 0)
   best_lambda = cv_ridge$lambda.min
   >>> best_lambda
   2.746789
   

   在这里，我们使用交叉验证方法来识别平均误差最低的lambda值。所有重复的交叉验证步骤都是通过cv.glmnet()函数完成的。

5. 使用最优lambda拟合一个新的岭回归模型，并提取没有截距的系数：

   
   # Fit a new ridge regression model using the optimal lambda
   opt_ridge_model = glmnet(X, y, alpha = 0, lambda = best_lambda)
   # Get coefficients
   ridge_coefs = coef(opt_ridge_model)[-1]  # remove intercept
   >>> ridge_coefs
   [1] -0.371840170 -0.005260088 -0.011611491  1.054511975 -1.233657799  0.162231830
    [7]  0.771141047  1.623031037  0.544153807 -0.547436697
   

6. 拟合线性回归模型并提取其系数：

   
   # Ordinary least squares regression
   ols_model = lm(mpg ~ ., data = mtcars)
   # Get coefficients
   ols_coefs = coef(ols_model)[-1] # remove intercept
   >>> ols_coefs
           cyl        disp          hp        drat          wt        qsec          vs
   -0.11144048  0.01333524 -0.02148212  0.78711097 -3.71530393  0.82104075  0.31776281
            am        gear        carb
    2.52022689  0.65541302 -0.19941925
   

7. 在同一张图上绘制两个模型的系数：

   
   plot(1:length(ols_coefs), ols_coefs, type="b", col="blue", pch=19, xlab="Coefficient", ylab="Value", ylim=c(min(ols_coefs, ridge_coefs), max(ols_coefs, ridge_coefs)))
   lines(1:length(ridge_coefs), ridge_coefs, type="b", col="red", pch=19)
   legend("bottomright", legend=c("OLS", "Ridge"), col=c("blue", "red"), pch=19)
   

   运行此代码生成图 12**.7：

图 12.7 – 可视化岭回归和 OLS 模型的估计系数

此图显示，岭回归模型的估计系数通常小于 OLS 模型的估计系数。

下一个部分将重点介绍 lasso 回归。

使用 lasso 回归进行工作

Lasso 回归是另一种正则化线性回归。它与岭回归相似，但在计算系数大小具体过程上有所不同。具体来说，它使用系数的 L1 范数，即系数绝对值总和的总和，作为添加到 OLS 损失函数中的惩罚项。

Lasso 回归的成本函数可以写成如下形式：

L lasso = RSS + λ∑ j=1 p | β j|

Lasso 回归的关键特性是它可以精确地将一些系数减少到 0，从而有效地执行变量选择。这是 L1 惩罚项的结果，而岭回归只能将系数缩小到接近 0。因此，当我们认为只有预测变量的一部分对预测结果有影响时，Lasso 回归特别有用。

此外，与岭回归不同，岭回归无法执行变量选择，因此可能不太容易解释，而 Lasso 回归自动选择最重要的特征并丢弃其余特征，这使得最终模型更容易解释。

与岭回归一样，Lasso 回归的惩罚项也受一个调整参数λ的影响。最佳λ参数通常通过交叉验证或更智能的搜索策略，如贝叶斯优化来选择。

让我们通过一个练习来了解如何开发 Lasso 回归模型。

练习 12.6 – 实现 Lasso 回归

要实现 Lasso 回归模型，我们可以遵循与岭回归模型类似的过程：

1. 要拟合 Lasso 回归模型，在glmnet()函数中将alpha设置为1：

   
   lasso_model = glmnet(X, y, alpha = 1)
   

2. 使用相同的交叉验证过程来识别lambda的最佳值：

   
   # Use cross-validation to find the optimal lambda
   cv_lasso = cv.glmnet(X, y, alpha = 1)
   best_lambda = cv_lasso$lambda.min
   >>> best_lambda
   0.8007036
   

3. 使用最佳lambda值拟合一个新的 Lasso 回归模型：

   
   # Fit a new lasso regression model using the optimal lambda
   opt_lasso_model = glmnet(X, y, alpha = 1, lambda = best_lambda)
   

   也可以使用coef()函数提取结果系数，然后跟[-1]以移除截距项：

   # Get coefficients
   lasso_coefs = coef(opt_lasso_model)[-1]  # remove intercept
   >>> lasso_coefs
   [1] -0.88547684  0.00000000 -0.01169485  0.00000000 -2.70853300  0.00000000  0.00000000
    [8]  0.00000000  0.00000000  0.00000000
   

4. 将估计系数与前面两个模型一起绘制：

   
   plot(1:length(ols_coefs), ols_coefs, type="b", col="blue", pch=19, xlab="Coefficient", ylab="Value", ylim=c(min(ols_coefs, ridge_coefs), max(ols_coefs, ridge_coefs)))
   lines(1:length(ridge_coefs), ridge_coefs, type="b", col="red", pch=19)
   lines(1:length(lasso_coefs), lasso_coefs, type="b", col="green", pch=19)
   legend("bottomright", legend=c("OLS", "Ridge", "Lasso"), col=c("blue", "red", "green"), pch=19)
   

   运行此代码生成图 12.8，这表明结果模型中只保留了两个变量。因此，Lasso 回归模型可以通过将某些特征的系数设置为 0 来生成稀疏解：

图 12.8 – 可视化岭回归、Lasso 回归和 OLS 回归模型的估计系数

总结，Lasso 回归模型通过将非显著变量的系数设置为 0，从而实现特征选择和模型估计。

概述

在本章中，我们介绍了线性回归模型的基本原理。我们首先介绍了单变量线性回归模型（SLR），它只包含一个输入变量和一个目标变量，然后扩展到多变量线性回归模型（MLR），它包含两个或更多预测变量。这两个模型都可以使用 R²来评估，或者更偏好使用调整后的 R²指标。接下来，我们讨论了特定场景，例如处理分类变量和交互项，通过变换处理非线性项，使用闭式解，以及处理多重共线性异方差性。最后，我们介绍了广泛使用的正则化技术，如岭回归和 Lasso 惩罚，这些技术可以作为损失函数中的惩罚项被纳入，生成正则化模型，在 Lasso 回归的情况下，还可以生成稀疏解。

在下一章中，我们将介绍另一种广泛使用的线性模型：逻辑回归模型。

第十三章：R 中的逻辑回归

本章将介绍逻辑回归，包括其理论结构、与线性回归的联系以及实际应用。由于它是一个在解释性重要的领域（如信用风险建模）中广泛使用的分类模型，我们将关注其在不同情境下的建模过程，以及添加正则化到损失函数和预测多于两个类别的扩展。

到本章结束时，你将理解逻辑回归模型的基本原理及其与线性回归的比较，包括扩展概念，如sigmoid函数、优势比和交叉熵损失（CEL）。你还将掌握在分类设置中常用的评估指标，以及处理不平衡数据集和目标变量中多个类别的改进方法。

本章将涵盖以下内容：

• 介绍逻辑回归

• 比较逻辑回归与线性回归

• 更多关于对数优势比和优势比的内容

• 介绍交叉熵损失

• 评估逻辑回归模型

• 处理不平衡数据集

• 惩罚逻辑回归

• 扩展到多类分类

技术要求

要运行本章中的代码，你需要以下软件包的最新版本：

• caret – 6.0.94

• tibble – 3.2.1

• dplyr – 1.0.10

• pROC – 1.18.2

• nnet – 7.3.18

• glmnet – 4.1.7

上述提到的版本以及前述列表中的软件包都是我在写这本书时的最新版本。

本章的所有代码和数据均可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_13/working.R找到。

介绍逻辑回归

逻辑回归是一种二元分类模型。它仍然是一个线性模型，但现在输出被限制为二元变量，取值为0或1，而不是像线性回归那样模拟连续结果。换句话说，我们将观察和模拟结果 y = 1 或 y = 0。例如，在信用风险建模的情况下，y = 0 表示非违约贷款申请，而 y = 1 表示违约贷款。

然而，逻辑回归模型不是直接预测二元结果，而是预测 y 取特定值的概率，例如 P(y = 1)。假设其他类别的概率 P(y = 0) = 1 − P(y = 1)，因为总概率应该始终加起来为1。最终的预测将是两个中的胜者，如果 P(y = 1) > P(y = 0)，则取值为1，否则为0。在信用风险示例中，P(y = 1)将被解释为贷款违约的概率。

在逻辑回归中，术语 logistic 与 logit 有关，它指的是对数几率。几率是描述概率的另一种方式；不是指定个体 P(y = 1) 和 P(y = 0)，而是指 P(y = 1) 与 P(y = 0) 的比率。因此，对数几率是通过 logP(y = 1) - P(y = 0) 计算的。因此，我们可以简单地使用术语 odds 来描述事件发生的概率（y = 1）与事件不发生的概率（y = 0）之间的比率。

首先，让我们看看逻辑回归模型如何将连续输出（如线性回归中的情况）转换为概率分数，这是一个介于 0 和 1 之间的数字。

理解 sigmoid 函数

sigmoid 函数是关键成分，它将任何连续数字（从负无穷大到正无穷大）映射到概率。也称为逻辑函数，sigmoid 函数的特点是 S 形曲线，它将任何实数作为输入，并将其映射到 0 和 1 之间的分数，这恰好是有效概率分数的范围。

标准的 sigmoid 函数具有以下形式：

f(x) = 1 / (1 + e^(-x))

注意，这是一个非线性函数。也就是说，当通过这个函数进行转换时，输入值将得到不成比例的缩放。它也是一个连续函数（因此可导）和单调函数（f(x) 将随着 x 的增加而增加），因此作为二分类任务中典型神经网络模型最后一层的激活函数，它享有很高的流行度。

让我们尝试可视化这个函数。在下面的代码片段中，我们使用 seq(-10, 10, by = 0.1) 创建一个从 -10 到 10 的等间距数字序列，步长为 0.1。对于每个数字，我们使用 sigmoid 函数计算相应的输出。在这里，我们直接传递函数的所有数字，然后以并行模式（称为向量化）计算所有输出。在这里，向量化是指同时将一个操作应用于整个向量，而不是像 for 循环那样逐个元素循环。最后，我们绘制函数以显示 sigmoid 函数的特征 S 形曲线，并使用 grid() 函数添加网格线：

# Create a range of equally spaced numbers between -10 and 10
x = seq(-10, 10, by = 0.1)
# Calculate the output value for each number using sigmoid function
sigmoid = 1 / (1 + exp(-x))
# Plot the sigmoid function
plot(x, sigmoid, type = "l", lwd = 2,
     main = "Sigmoid Function",
     xlab = "x",
     ylab = "f(x)",
     col = "red")
# Add grid lines
grid()

运行前面的代码生成了 图 13.1 中所示的输出。该图显示了整个域内不同级别的陡峭程度，其中函数在中间区域更敏感，在两个极端处变得更加饱和。

图 13.1 – sigmoid 函数的可视化

现在我们已经理解了 sigmoid 函数，让我们看看逻辑回归模型的数学结构。

理解逻辑回归模型

逻辑回归模型本质上是一种线性回归模型，它被推广到依赖结果变量为二元的设置中。换句话说，它是一个线性回归模型，用于建模事件发生概率的对数几率。

为了理解这一点，让我们首先回顾以下线性回归模型，其中我们使用总共 p 个特征来建模目标输出变量 z：

z = β0 + β1 x1 + β2 x2 + … + βp xp

在这里，z 被解释为事件 y = 1 的对数几率，或 logit。我们感兴趣的是从 β0 到 βp 的参数估计。

现在，我们知道 z 变量是无界的，这意味着它可以从负无穷大到正无穷大变化。我们需要一种方法来限制这个输出并将其转换为介于 0 和 1 之间的概率分数。这是通过使用 sigmoid 函数的额外转换来实现的，它恰好满足我们的所有需求。从数学上讲，我们有以下公式：

P(y = 1) = 1 / (1 + e^(-z))

将 z 的定义代入，给出完整的逻辑回归模型：

P(y = 1) = 1 / (1 + e^(-β0 - β1x1 - β2x2 - … - βpxp))

在这里，P(y = 1) 指的是 y = 1 成功的概率（这是一个一般性陈述），相应地，P(y = 0) 表示失败的概率。

注意，我们可以等效地将模型表达如下：

log P(y = 1) / P(y = 0) = β0 + β1 x1 + β2 x2 + … + βp xp

在这里，术语 log P(y = 1) / P(y = 0) 表示对数几率。

这里一个关键的变化是引入了 sigmoid 函数。这使得预测变量与结果概率之间的关系不再是线性的，而是 sigmoid 形式的。为了观察这里的微妙之处，我们可以查看 sigmoid 函数在整个定义域中的不同区域。例如，当查看 0 附近的区域时，输入的微小变化会导致结果概率输出的相对较大变化。然而，当输入位于函数的两个极端一侧时，同样的输入变化将导致输出非常小的变化。此外，随着输入变得更加极端，无论是向负方向还是正方向，结果概率将逐渐接近 0 或 1。图 13.2 总结了逻辑回归模型的特点。

图 13.2 – 总结逻辑回归模型

注意，逻辑回归模型与线性回归模型具有相似的假设。具体来说，它假设观测值之间相互独立，并且目标结果遵循参数为 p 的伯努利分布，即 y ∼ Bernoulli(p)。鉴于这一点，我们并不假设输出变量和输入预测之间存在线性回归；相反，我们使用逻辑链接函数来转换并引入非线性到输入变量中。

以下部分进一步比较逻辑回归与线性回归。

比较逻辑回归与线性回归

在本节中，我们将专注于使用German Credit数据集的二进制信用分类任务，该数据集包含 1,000 个观测值和 20 列。每个观测值表示一位客户，该客户曾向银行申请贷款，并按信用风险被标记为好或坏。该数据集在 R 的caret包中可用。

对于我们的研究，我们将尝试根据Duration预测目标二元变量Class，并比较线性回归和逻辑回归在预测结果上的差异。我们特意选择一个预测变量，以便在二维图中可视化和比较结果模型的决策边界。

练习 13.1 – 比较线性回归与逻辑回归

在这个练习中，我们将展示使用逻辑回归模型产生概率输出的优势，与使用线性回归模型产生的无界输出相比：

1. 从caret包中加载German Credit数据集。将目标变量（Class）转换为数值：

   
   # install the caret package if you haven't done so
   install.packages("caret")
   # load the caret package
   library(caret)
   # load the German Credit dataset
   data(GermanCredit)
   GermanCredit$Class_num = ifelse(GermanCredit$Class == "Bad", 1, 0)
   

   在这里，我们创建一个新的目标变量Class_num，将原始的Class变量映射到1，如果它取值为"Bad"，否则为0。这是必要的，因为线性回归和逻辑回归模型都不能接受基于字符串的变量作为目标（或预测变量）。

2. 建立一个线性回归模型，将Class_num对Duration进行回归：

   
   lm_model = lm(Class_num ~ Duration, data=GermanCredit)
   coefs = coefficients(lm_model)
   intercept = coefs[1]
   slope = coefs[2]
   

   在这里，我们使用lm()函数构建线性回归模型，并使用coefficients()提取模型系数，包括截距和斜率。

3. 可视化预测和目标：

   
   ggplot(GermanCredit,
          aes(Duration, Class_num)) +
     geom_point() +
     geom_abline(intercept=intercept, slope=slope) +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   在这里，我们将观察到的目标变量作为散点图绘制，并使用geom_abline()函数根据估计的斜率和截距绘制模型为直线。

   运行前面的代码生成图 13.3。

图 13.3 – 可视化线性回归模型

由于所有目标值都是0或1，我们可以将预测视为介于0和1之间的概率值。然而，当我们放大视图时，无界概率的问题就会显现出来，如下面的步骤所示。

1. 通过将x轴的范围扩大到(-30, 120)和y轴的范围扩大到(-0.5, 1.5)来重新绘制图形：

   
   ggplot(GermanCredit,
          aes(Duration, Class_num)) +
     geom_point() +
     geom_abline(intercept=intercept, slope=slope) +
     xlim(-30, 120) +
     ylim(-0.5, 1.5) +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   在这里，我们使用xlim()和ylim()函数扩大了x轴和y轴可能值的范围。

   运行前面的代码生成 图 13.4 中所示的输出，这表明当 Duration 的值变得极端时，预测值超出了 [0,1] 的范围，这种情况称为超出观察值范围的外推。这意味着预测概率将小于 0 或大于 1，这显然是一个无效的输出。这需要一种称为逻辑回归的广义线性回归模型，其中响应将根据 sigmoid 函数的转换遵循逻辑、S 形曲线。

图 13.4 – 可视化扩展范围的线性回归模型

1. 使用 glm() 函数构建逻辑回归模型：

   
   glm_model = glm(Class_num ~ Duration, data=GermanCredit, family=binomial)
   >>> glm_model
   Call:  glm(formula = Class_num ~ Duration, family = binomial, data = GermanCredit)
   Coefficients:
   (Intercept)     Duration
      -1.66635      0.03754
   Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
   Null Deviance:     1222
   Residual Deviance: 1177   AIC: 1181
   

   结果显示了逻辑回归模型的估计截距和斜率，以及残差偏差，这是衡量模型拟合优度的一个指标。

2. 在之前的图形上绘制估计的逻辑曲线：

   
   ggplot(GermanCredit,
          aes(Duration, Class_num)) +
     geom_point() +
     geom_abline(intercept=intercept, slope=slope) +
     geom_smooth(
       method = "glm",
       se = FALSE,
       method.args = list(family=binomial)
     ) +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   运行前面的代码将生成 图 13.5 中所示的输出，这表明线性回归线和逻辑回归曲线之间存在轻微的偏差。

图 13.5 – 可视化逻辑回归曲线

再次，我们可以放大图形并关注当超出数据集中可能值的观察范围时的差异。

1. 在更宽的值范围内绘制逻辑曲线，超出观察值范围：

   
   # Get coefficients from logistic model
   intercept_glm = coef(glm_model)[1]
   slope_glm = coef(glm_model)[2]
   # Generate sequence of x-values
   x_values = seq(from = min(GermanCredit$Duration) - 150,
                   to = max(GermanCredit$Duration) + 150,
                   by = 0.1)
   # Compute probabilities using logistic function
   y_values = 1 / (1 + exp(-(intercept_glm + slope_glm * x_values)))
   # Data frame for plot
   plot_df = data.frame(x = x_values, y = y_values)
   # Plot
   ggplot() +
     geom_point(data = GermanCredit, aes(Duration, Class_num)) +
     geom_abline(intercept=intercept, slope=slope) +
     geom_line(data = plot_df, aes(x, y), color = "blue") +
     theme_minimal() +
     xlim(-30, 120) +
     ylim(-0.5, 1.5) +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   在这里，我们首先提取逻辑回归模型的系数，然后使用 sigmoid 函数转换生成输入值的序列及其相应的输出。最后，我们在同一图形中与观察数据一起绘制逻辑曲线和线性拟合。

   运行前面的代码生成 图 13.6 中所示的输出，这表明随着输入值的极端化，逻辑回归曲线逐渐饱和。此外，所有值现在都被限制在 [0,1] 的范围内，这使得它成为将解释为概率而不是无界值的有效候选。

图 13.6 – 可视化扩展范围的逻辑回归模型

下一节将探讨如何使用逻辑回归模型进行预测。

使用逻辑回归模型进行预测

如前所述，逻辑回归模型的直接预测形式是介于0和1之间的概率。要将它们转换为二元预测，我们可以通过使用0.5的阈值对概率进行四舍五入来实现。例如，如果预测概率是 P(y = 1) = 0.8，四舍五入操作将导致最终的二元预测 y = 1。另一方面，如果 P(y = 1) = 0.3，四舍五入将导致 y = 0。

让我们通过以下练习来了解如何使用逻辑回归模型进行预测。

练习 13.2 – 使用逻辑回归模型进行预测

我们已经看到了如何通过提取逻辑回归模型的斜率和截距来执行预测。在这个练习中，我们将探索使用predict()函数的更便捷的方法：

1. 生成一个从5到80的Duration值序列，步长为2，并使用之前的逻辑回归模型通过predict()函数预测该序列的对应概率：

   
   library(tibble)
   library(dplyr)
   # making predictions
   pred_df = tibble(
     Duration = seq(5, 80, 2)
   )
   pred_df = pred_df %>%
     mutate(
       pred_prob = predict(glm_model, pred_df, type="response")
     )
   

   这里，我们使用seq()函数创建等间隔向量，并将其存储在名为pred_df的tibble对象中。然后我们使用predict()通过指定type="response"来预测相应的概率。

2. 将预测概率与原始数据一起可视化：

   
   ggplot() +
     geom_point(data = GermanCredit, aes(Duration, Class_num)) +
     geom_point(data = pred_df, aes(Duration, pred_prob) , color="blue") +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   上述代码将生成以下输出：

图 13.7 – 可视化预测的概率

1. 使用round()函数将概率转换为二元结果：

   
   # getting the most likely outcome
   pred_df = pred_df %>%
     mutate(
       most_likely_outcome = round(pred_prob)
     )
   

   这里，我们使用默认阈值0.5对预测概率进行四舍五入。

2. 将二元结果作为绿色点添加到之前的图中：

   
   ggplot() +
     geom_point(data = GermanCredit, aes(Duration, Class_num)) +
     geom_point(data = pred_df, aes(Duration, pred_prob) , color="blue") +
     geom_point(data = pred_df, aes(Duration, most_likely_outcome) , color="green") +
     theme(axis.title.x = element_text(size = 18),
           axis.title.y = element_text(size = 18))
   

   运行上述代码将生成以下输出，表明所有预测概率高于0.5的都被转换为1，而低于0.5的都被转换为0。

图 13.8 – 可视化预测的二元结果

下一个部分将进一步讨论对数赔率。

关于对数赔率和赔率比率的更多内容

记住，赔率是指事件发生的概率与其补数的比率：

odds =  事件发生的概率   ________________________   事件不发生的概率  =  p _ 1 − p  =  P(y = 1) _ P(y = 0)

这里，概率的计算如下：

p = P(y = 1) =  1 _ 1 + e −z

1 − p = 1 −  1 _ 1 + e −z  =  e −z _ 1 + e −z

将 p 和 1 − p 的定义代入，我们得到以下结果：

odds =  p _ 1 − p  = e z

我们通常不直接处理赔率，而是使用对数赔率或 logit。这个术语通常通过以下方式在逻辑回归模型中作为预测器的线性组合进行建模：

log P(y = 1) _ P(y = 0)  = z = β 0 + β 1 x 1 + … + β p x p

在这里，我们可以将每个系数 β j 解释为在保持所有其他预测变量不变的情况下，jth 预测变量 x j 单位增加时对对数优势的预期变化。这个方程本质上表明，目标值 y 为 1 的对数优势与输入变量呈线性关系。

现在假设 x i 是一个二元输入变量，使得 x i = 1 或 0。我们可以如下计算 x i = 1 的优势：

p 1 _ 1 − p 1  = e z 1

这衡量了 x i = 1 事件发生的概率与不发生事件的概率。同样，我们可以如下计算 x i = 0 的优势：

p 0 _ 1 − p 0  = e z 0

这衡量了 x i = 0 事件发生的概率与不发生事件的概率。

然后，我们可以计算 x i 的优势比，这是 x i = 1 的优势与 x i = 0 的优势之比：

p 1 _ 1 − p 1 _  p 0 _ 1 − p 0  =  e β 0+β 1x 1+…β i1+…+β px p ___________ e β 0+β 1x 1+…β i0+…+β px p  = e β i

在这里，e β i 衡量了二元输入变量 x i 对结果 y 为 1 的优势的量化影响，同时保持所有其他输入变量不变。这为我们提供了一种衡量逻辑回归模型中任何预测变量影响的方法，包括分类和数值输入变量。

对于分类输入变量，我们可以使用性别（男性为 0，女性为 1）来预测是否购买保险（购买为 1，不购买为 0）。我们将男性的基准类别设置为 0。如果性别估计系数 β gender = 0.2，其优势比计算为 e 0.2 ≈ 1.22。因此，女性客户购买保险的优势是男性客户购买保险优势的 1.22 倍，假设所有其他变量保持不变。

对于数值输入变量，我们可以使用年龄来预测是否购买保险。在这种情况下，无需设置基准类别。如果估计的系数 β age = 0.3，相应的优势比计算为 e 0.3 ≈ 1.35。这意味着客户的购买概率是比一年年轻的人高 1.35 倍，假设其他条件相同（ceteris paribus）。

注意，我们可以使用预测概率来计算对数优势，如下面的代码片段所示：

# calculate log odds using predicted probabilities
pred_df = pred_df %>%
  mutate(
    log_odds = log(pred_prob / (1 - pred_prob))
  )

下一节将介绍更多关于逻辑回归模型损失函数的内容。

引入交叉熵损失

二元交叉熵损失，也称为 log loss，通常用作逻辑回归中的成本函数。这是逻辑回归模型将通过移动参数来尝试最小化的损失。这个函数将预测概率和相应的目标作为输入，并输出一个标量分数，表示拟合优度。对于具有目标 y i 和预测概率 p i 的单个观察值，损失计算如下：

Q i(y i, p i) = − [ y i logp i + (1 − y i)log(1 − p i)]

将所有个别损失相加得到总二元交叉熵损失：

Q(y, p) =  1 _ N  ∑ i N Q i =  1 _ N  ∑ i=1 N − [ y i logp i + (1 − y i)log(1 − p i)]

二元 CEL 函数是二元分类问题的合适选择，因为它对自信但错误的预测进行重罚。例如，当 p i 接近 0 或 1 时，如果预测错误，则 resulting CEL 将趋于无穷大。因此，这一特性鼓励学习过程输出接近目标真实概率的概率。

更一般地，我们使用 CEL 来模拟目标变量中有两个或更多类别的分类问题。对于第 i 个观察 x i，分类函数将产生一个概率输出，表示为 p i,k = f( x i; w)，以指示属于第 k 个类别的可能性。当我们有一个包含总共 C 个类别的分类任务时，x i 的 CEL 定义为 Q i(w) = − ∑ k=1 C  y i,k log(p i,k)，这实际上是对所有类别的求和。再次强调，如果第 i 个观察的目标标签属于第 k 个类别，则 y i,k = 1，否则 y i,k = 0。

求和意味着对每个类别执行类内评估（即，项 y i,k log(p i,k)），并将它们相加以产生第 i 个观察的总成本。对于每个观察，总共有 C 个预测与属于每个类别的相应概率相对应。因此，CEL 通过将它们求和成一个单一的数字来聚合预测概率矩阵。此外，结果取反以产生一个正数，因为当 x 是介于 0 和 1 之间的概率时，log(x) 是负数。

注意，CEL 计算中的目标标签需要是一维编码。这意味着单个分类标签被转换成一个包含 1 的二进制数字数组，表示类别标签，其他为 0。例如，对于一个数字 8 的图像，而不是传递原始类别作为目标输出，将使用结果的一维编码目标 [0, 0, 0, 0, 0, 0, 0, 1, 0, 0]，其中第 8 位被激活（即，热），其余被禁用。目标数组也将与概率输出数组具有相同的长度，其元素对应于每个类别。

直观上，我们希望预测的正确类别的概率接近 1，而错误类别的概率为 0。也就是说，随着预测概率偏离实际类别标签，损失应该增加。CEL 设计是为了遵循这种直觉。具体来说，我们可以查看第 i 个观察的以下四种情况：

• 当目标标签属于第 k 个类别（即，y i,k = 1）且第 k 个类别的预测概率非常强（即，p i,k ≈ 1）时，成本应该很低。

• 当目标标签属于第 k 个类别（即，y i,k = 1）且第 k 个类别的预测概率非常弱（即，p i,k ≈ 0）时，成本应该很高。

• 当目标标签不属于第 k 个类别（即 y_i,k = 0）且第 k 个类别的预测概率非常强（即 p_i,k ≈ 1）时，成本应该高。

• 当目标标签不属于第 k 个类别（即 y_i,k = 0）且第 k 个类别的预测概率非常弱（即 p_i,k ≈ 0）时，成本应该低。

要计算 CEL，我们首先计算每个观察值在所有类别中二进制标签向量与预测概率向量之间的加权求和。将所有观察值的结果相加，然后加一个负号以将成本转换为正数。CEL 的设计是为了匹配成本的感觉：当预测和目标非常接近时，成本应该低，否则高。换句话说，为了计算总成本，我们只需将单个成本相加，导致 Q(w) = − ∑_i=1^N ∑_k=1^C y_i,k log(p_i,k)。图 13.9 总结了关于 CEL 的前述讨论。

图 13.9 – 说明 CEL

注意，最后两种情况对损失计算没有任何贡献，因为目标值等于 0；任何乘以 0 的值都变为 0。此外，观察这些两个类别的特定观察值的预测概率需要加起来等于 1，这使得我们只需要关注一个类别的预测概率（通常是类别 1）。

下一节将介绍如何评估逻辑回归模型。

评估逻辑回归模型

我们可以使用多个指标来评估逻辑回归模型。这些是我们用来确定拟合优度（在测试集上）的指标，需要与我们在训练模型时使用的 CEL（在训练集上）区分开来。

以下列表提供了常用的指标：

• 准确率：这是模型对观察值进行正确预测的数量与观察值总数的比例。由于正确的预测可以是真正的阳性或真正的阴性，准确率是通过将真正的阳性和真正的阴性相加，然后除以观察值总数来计算的。

• 错误率：这是模型在总观察值中做出错误预测的比例。错误的预测可以是假阳性或假阴性。它被计算为 1 - 准确率；也就是说，错误率是准确率的补数。换句话说，它被计算为 *(假阳性 + 假阴性) / 总观察值。

• 精确度：精确度是所有正预测中正确正预测的比例。这个指标实际上告诉我们，在所有预测为正的实例中，有多少是正确的。因此，它表明模型在预测正观察方面的准确性，并计算为真正例 / (真正例 + 假正例)。在分母中，我们注意到在所有正实例中，一些是真正例，其余是假正例。

• 召回率：召回率指的是模型正确预测的实际正实例的比例。也称为敏感性或真正例率（TPR），召回率衡量模型检测正观察的能力。它计算为真正例 / (真正例 + 假反例)，其中分母的公式与精确度不同。

• 特异性：也称为真正例率（TNR），特异性衡量模型正确预测的实际负实例的比例。这与敏感性相反，敏感性关注模型捕捉真正例的能力。换句话说，特异性衡量模型正确识别负实例或非事件的能力。它计算为真正例 / (真正例 + 假正例)。

• 0 和 1 表示两个类别之间分离的程度。一个完美模型，如果预测正确率达到 100%，其 AUC（曲线下面积）为1，而一个完全错误的模型，其 AUC 为0。一个随机猜测（以 50%的概率选择0或1）的模型对应于 AUC 为0.5，这表明没有类别分离能力。

注意，在评估两个具有同等良好评估指标的模型时，我们将遵循简约性原则。简约性原则指出，如果两个竞争模型对数据的拟合程度相似，那么应该选择输入变量更少的模型，从而更倾向于简单而非复杂。其基本假设是最准确的模型不一定是最好的模型。

图 13*.10* 描述了构建混淆矩阵的过程，该矩阵捕捉了不同场景下的预测结果，以及计算上述评估指标的具体细节。

图 13.10 – 展示二元分类任务中的混淆矩阵和常见评估指标

注意，我们只能在选择一个阈值来截断预测概率之后才能计算这些指标。具体来说，预测概率大于截止阈值的观测值将被分类为正类，否则为负类。

精确率和召回率通常与分类阈值的调整呈反向关系。在目标变量值存在巨大不平衡的情况下，同时审查精确率和召回率是有用的。由于精确率和召回率是两个相关但不同的指标，我们应该优化哪一个？

要回答这个问题，我们需要评估做出错误阴性预测的相对影响。这通过错误阴性率来衡量，它是召回率的相反（或补数）。如图 13.11 所示，未能识别垃圾邮件的风险低于错过欺诈交易或阳性癌症患者。我们应该旨在优化精确度（以便模型的预测更加精确和有针对性）对于第一种情况，以及召回率（以便最小化错过潜在阳性案例的机会）对于第二种情况。

图 13.11 – 具有不同错误阴性预测影响的三个案例

至于 AUC，或 ROC 曲线下方的面积，不需要选择特定的阈值，因为它通过评估从 0 到 1 的阈值序列来计算。ROC 曲线在 y 轴上绘制灵敏度，在 x 轴上绘制 1 - 特异性。这也对应于绘制 TPR 与 1-TNR，或 TPR 与 FPR。随着分类阈值的上升，FPR 下降，导致曲线向左移动。

完美的二元分类器具有 AUC 得分为 1。这意味着 FPR，或 1 – 特异性，为 0。也就是说，没有假阳性，所有负例都没有被预测为阳性。此外，灵敏度，或 TPR，为 1，意味着所有正例都被正确预测为阳性。

图 13.12 展示了三种不同的 AUC 曲线。最上面的曲线（绿色）对应于更好的模型，因为其 AUC 值最高。两种模型，由绿色和红色曲线表示，表现优于随机猜测，如蓝色直线非对角线所示。

图 13.12 – 展示三种不同的 AUC 曲线

继续进行之前的练习，我们现在可以计算相应的评估指标。首先，我们使用训练好的逻辑回归模型对训练集中的所有观测值进行评分，并使用 predict() 函数和设置 type="response" 获得相应的概率，如下面的代码片段所示。请注意，我们需要传递一个包含相应特征名称的 DataFrame 作为模型的输入：

# Create new data frame with all durations
new_data = data.frame(Duration = GermanCredit$Duration)
# Calculate predicted classes based on predicted probabilities
predicted_probs = predict(glm_model, new_data, type="response")

接下来，我们设置一个单一的截止阈值（在这种情况下为 0.5）使用 ifelse() 函数将预测概率转换为相应的二元结果：

# Convert to binary outcomes
predicted_classes = ifelse(predicted_probs > 0.5, 1, 0)

对于二元结果和真实目标标签，我们可以通过以下方式获得混淆矩阵：

# Create confusion matrix
conf_matrix = table(predicted = predicted_classes, actual = GermanCredit$Class_num)
>>> conf_matrix
         actual
predicted   0   1
        0 670 260
        1  30  40

在这里，混淆矩阵提供了模型正确和错误分类的细分。在混淆矩阵中，左上角的单元格表示真正的负例，右上角的单元格表示错误的正例，左下角的单元格表示错误的负例，而右下角的单元格表示真正的正例。

根据混淆矩阵，我们可以计算评估指标如下：

# Accuracy
accuracy = sum(diag(conf_matrix)) / sum(conf_matrix)
>>> print(paste("Accuracy: ", accuracy))
"Accuracy:  0.71"
# Error rate
error_rate = 1 - accuracy
>>> print(paste("Error rate: ", error_rate))
"Error rate:  0.29"
# Precision
precision = conf_matrix[2,2] / sum(conf_matrix[2,])
print(paste("Precision: ", precision))
"Precision:  0.571428571428571"
# Recall / Sensitivity
recall = conf_matrix[2,2] / sum(conf_matrix[,2])
print(paste("Recall: ", recall))
>>> "Recall:  0.133333333333333"
# Specificity
specificity = conf_matrix[1,1] / sum(conf_matrix[,1])
print(paste("Specificity: ", specificity))
>>> "Specificity:  0.957142857142857"

在这里，我们提取混淆矩阵的相应项，并将不同评估指标的定义代入以完成计算。

我们也可以计算 AUC，从使用pROC包计算 ROC 曲线开始：

library(pROC)
# Calculate ROC curve
roc_obj = roc(GermanCredit$Class_num, predicted_probs)
# Plot ROC curve
>>> plot(roc_obj)

运行代码生成图 13.13，这表明模型的表现略好于随机猜测。

图 13.13 – 可视化 ROC 曲线

要计算 AUC，我们可以调用auc()函数：

# Calculate AUC
auc = auc(roc_obj)
>>> print(paste(«AUC: «, auc))
«AUC:  0.628592857142857»

注意，当我们对精确度或召回率没有偏好时，我们可以使用 F1 分数，其定义如下：

F1 分数 = 2(精确度 × 召回率) / (精确度 + 召回率)

下一个部分讨论了一个具有挑战性的建模情况，即我们一开始就有一个不平衡的数据集。

处理不平衡数据集

当使用目标为二元结果的数据库构建逻辑回归模型时，目标值可能不是均匀分布的。这意味着我们会观察到比事件（y = 1）更多的非事件（y = 0），这在银行欺诈交易、企业员工的垃圾邮件/钓鱼邮件、癌症等疾病的识别以及地震等自然灾害等应用中通常是这种情况。在这些情况下，分类性能可能会被多数类主导。

这样的主导地位可能导致误导性的高准确率得分，这与较差的预测性能相对应。为了理解这一点，假设我们正在使用一个包含 1,000 个观察值的数据库来开发一个默认预测模型，其中只有 10 个（或 1%）是违约案例。一个简单的模型会简单地预测每个观察值都不是违约，从而得到 99%的准确率。

当我们遇到不平衡数据集时，我们通常对少数类更感兴趣，这代表了分类问题中要检测的结果。由于少数类的信号相对较弱，我们需要依赖良好的建模技术来识别良好的模式，以便正确检测信号。

我们可以使用多种技术来解决不平衡数据集的挑战。我们将介绍一种流行的称为数据重采样的方法，它需要原始数据集的过采样和/或欠采样，以使整体分布不那么不平衡。重采样包括对少数类进行过采样、对多数类进行欠采样，或使用它们的组合，如**合成少数过采样技术（SMOTE）**所表示。然而，这种补救措施并非没有风险。在这里，过采样可能导致过拟合，因为更多的样本被添加到原始数据集中，而欠采样可能导致信息丢失，因为一些多数观测值被从原始数据集中移除。

图 13.14 展示了过采样少数类（左侧面板）和欠采样多数类（右侧面板）的过程。请注意，一旦基于平衡数据集建立模型，我们还需要在新测试集上对其进行校准，以便它在新的真实数据集上也能表现良好。

图 13.14 – 展示了过采样少数类和欠采样多数类的过程

让我们通过一个练习来了解如何在逻辑回归设置中执行欠采样或过采样。

练习 13.3 – 执行欠采样和过采样

在这个练习中，我们将基于欠采样和过采样创建两个人工数据集。然后，我们将使用混淆矩阵评估所得到的逻辑回归模型的表现：

1. 将原始数据集划分为训练集（70%）和测试集（30%）：

   
   set.seed(2)
   index = sample(1:nrow(GermanCredit), nrow(GermanCredit)*0.7)
   train = GermanCredit[index, ]
   test = GermanCredit[-index, ]
   

   在这里，我们随机采样一组索引，用于选择训练集的观测值，并将剩余的分配到测试集。

2. 检查训练集中的类比例：

   
   >>> table(train$Class_num)
     0   1
   504 196
   

   结果显示，多数类（类别 0）的大小是少数类（类别 1）的两倍多。

3. 根据类别标签将训练集划分为多数集和少数集：

   
   # separate the minority and majority classes
   table(train$Class_num)
   minority_data = train[train$Class_num == 1,]
   majority_data = train[train$Class_num == 0,]
   

   然后，我们将使用这两个数据集来执行欠采样和过采样。

4. 对多数类进行欠采样，并将欠采样的多数类与少数类合并。检查结果类的比例：

   
   # undersample the majority class
   undersampled_majority = majority_data[sample(1:nrow(majority_data), nrow(minority_data)),]
   # combine undersampled majority class and minority class
   undersampled_data = rbind(minority_data, undersampled_majority)
   >>> table(undersampled_data$Class_num)
     0   1
   196 196
   

   类比例现在平衡了。让我们对少数类执行过采样。

5. 对少数类进行过采样，并将过采样的少数类与多数类合并。检查结果类的比例：

   
   # oversample the minority class
   oversampled_minority = minority_data[sample(1:nrow(minority_data), nrow(majority_data), replace = TRUE),]
   # combine majority class and oversampled minority class
   oversampled_data = rbind(majority_data, oversampled_minority)
   >>> table(oversampled_data$Class_num)
     0   1
   504 504
   

6. 在欠采样和过采样的数据集上拟合逻辑回归模型：

   
   # fit logistic regression models on undersampled and oversampled data
   undersampled_model = glm(Class_num ~ Duration, family = binomial(link = 'logit'), data = undersampled_data)
   oversampled_model = glm(Class_num ~ Duration, family = binomial(link = 'logit'), data = oversampled_data)
   

7. 在测试集上获取预测概率：

   
   # get the predicted probabilities on the test set
   undersampled_pred = predict(undersampled_model, newdata = test, type = "response")
   oversampled_pred = predict(oversampled_model, newdata = test, type = "response")
   

8. 应用 0.5 的阈值将概率转换为二进制类别：

   
   # apply threshold to convert the probabilities into binary classes
   undersampled_pred_class = ifelse(undersampled_pred > 0.5, 1, 0)
   oversampled_pred_class = ifelse(oversampled_pred > 0.5, 1, 0)
   

9. 计算混淆矩阵：

   
   # calculate the confusion matrix
   undersampled_cm = table(predicted = undersampled_pred_class, actual = test$Class_num)
   oversampled_cm = table(predicted = oversampled_pred_class, actual = test$Class_num)
   >>> undersampled_cm
            actual
   predicted   0   1
           0 117  59
           1  79  45
   >>> oversampled_cm
            actual
   predicted   0   1
           0 115  59
           1  81  45
   

   结果表明，使用欠采样或过采样训练数据集，这两种模型都表现出相似的性能。再次，我们会使用另一个验证数据集，它可以来自原始训练集，以进一步校准模型参数，以便在测试集上表现更好。

结果表明，我们还可以在逻辑回归模型中添加 lasso 或岭惩罚，如下一节所述。

惩罚逻辑回归

如其名所示，惩罚逻辑回归模型在通常逻辑回归模型的损失函数中包含一个额外的惩罚项。回想一下，标准逻辑回归模型试图最小化负对数似然函数（或等价地，最大化对数似然函数），其定义如下：

Q(𝜷) =  1 _ N  ∑ i=1 N − [ y i logp i + (1 − y i)log(1 − p i)]

在这里，p i =  1 _ 1 + e −(β 0+β 1x 1 (i)+β 2x 2 (i)+…+β px p (i)) 是输入 x (i) 的预测概率，y i 是相应的目标标签，𝜷 = { β 0, β 1, … , β p} 是要估计的模型参数。请注意，我们现在将损失表示为系数向量的函数，因为它直接由模型中使用的参数集确定。

由于惩罚项旨在缩小估计系数的幅度，因此我们会将其添加到损失函数中，以便惩罚项将相对较小（受一个调整超参数λ的限制）。对于岭惩罚的情况，我们会添加平方系数的总和，从而得到以下惩罚负对数似然函数：

Q ridge(𝜷) = Q(𝜷) + λ ||𝜷|| 2 2 = Q(𝜷) + λ∑ j=1 p β j 2

相应地，使用 lasso 正则化项的惩罚负对数似然函数具有以下形式：

Q lasso(𝜷) = Q(𝜷) + λ|𝜷| = Q(𝜷) + λ∑ j=1 p | β j|

在这两种情况下，惩罚项都有可能将系数的幅度缩小到相对于原始最大似然估计的 0。这可以通过控制模型估计过程的复杂性来帮助防止过拟合。调整超参数λ控制缩放的程度。特别是，较大的λ会给惩罚项添加更多的权重，从而导致更多的缩放效应，而较小的λ则对估计系数的整体幅度赋予较少的权重。

让我们说明惩罚逻辑回归模型的发展过程。可以使用支持 lasso 和岭惩罚的glmnet包来完成此任务。在下面的代码片段中，我们使用前九列作为预测变量来建模二元结果：

# Create a matrix of predictors and a response vector
# For glmnet, we need to provide our data as matrices/vectors
X = GermanCredit[1:nrow(GermanCredit), 1:9]
y = GermanCredit$Class_num
# Define an alpha value: 0 for ridge, 1 for lasso, between 0 and 1 for elastic net
alpha_value = 1 # for lasso
# Run the glmnet model
fit = glmnet(X, y, family = "binomial", alpha = alpha_value)

在这里，我们设置 alpha=1 以启用 lasso 惩罚。设置 alpha=0 启用岭惩罚，将 alpha 设置为0和1之间对应于弹性网络惩罚。

注意，此过程将评估λ超参数的一系列值，从而让我们了解基于惩罚对结果系数的影响程度。特别是，我们可以绘制出系数路径，如下所示，表示不同λ值下的结果系数：

# plot coefficient paths
>>> plot(fit, xvar = "lambda", label = TRUE)

运行代码生成图 13.15，它表明当λ值变大时，更多参数会缩小到 0。

图 13.15 – 使用 lasso 惩罚逻辑回归可视化系数路径

下一节讨论一个更一般的设置：多项式逻辑回归，这是一个多类分类的模型类别。

扩展到多类分类

许多问题具有超过两个类别。例如，AAA、AA、A以及更多类似的情况。银行的企业客户账户被分类为良好信用、逾期、拖欠、可疑或损失。这种设置需要多项式逻辑回归模型，这是在多类分类背景下对二元逻辑回归模型的推广。本质上，目标变量 y 可以取超过两个可能的离散结果，并允许超过两个分类值。

假设目标变量可以取三个值，即 y ∈ {0,1, 2}。让我们选择类别0作为基线值。我们将建模其他类别（类别1和2）相对于这个基线的概率。换句话说，我们有以下：

p(y = 1) _ p(y = 0)  = e z 1

p(y = 2) _ p(y = 0)  = e z 2

因此，每个类别的预测概率的相对比例如下：

p(y = 2) : p(y = 1) : p(y = 0) = e z 2 : e z 1 : 1

我们知道以下内容：

p(y = 2) + p(y = 1) + p(y = 0) = 1

因此，我们有以下：

p(y = 2) =  e z 2 _ e z 2 + e z 1 + 1

p(y = 1) =  e z 1 _ e z 2 + e z 1 + 1

p(y = 0) =  1 _ e z 2 + e z 1 + 1

再次强调，多项式逻辑回归模型中的一个主要假设是，对数几率是预测变量的线性组合。这与二元逻辑回归中的假设相同。在多项式逻辑回归中，系数的解释也会略有变化。特别是，每个系数现在代表对应预测变量单位变化时，相对于基线类别的对数几率的变化，同时保持所有其他预测变量不变。

我们可以依赖nnet包中的multinom()函数来创建多项式逻辑回归模型。在下面的代码片段中，我们使用之前的mtcars数据集，并将gear变量转换为因子，这将作为目标变量：

library(nnet)
# convert gear to factor
mtcars$gear = as.factor(mtcars$gear)
>>> table(mtcars$gear)
3  4  5
15 12  5

频率计数显示变量中有三个类别。接下来，我们使用mpg、hp和disp在多项式逻辑回归模型中预测gear：

# fit the model
multinom_model = multinom(gear ~ mpg + hp + disp, data = mtcars)
# weights:  15 (8 variable)
initial  value 35.155593
iter  10 value 10.945783
iter  20 value 9.011992
iter  30 value 8.827997
iter  40 value 8.805003
iter  50 value 8.759821
iter  60 value 8.742738
iter  70 value 8.737492
iter  80 value 8.736569
final  value 8.735812
converged

输出信息表明，模型中共有八个变量。这是有道理的，因为我们有四个变量（截距、mpg、hp和disp）来模型一个子模型中四档和三档汽车之间的差异，另外四个变量来模型另一个子模型中五档和三档汽车之间的差异。

让我们看一下模型的摘要：

# view summary of the model
>>> summary(multinom_model)
Call:
multinom(formula = gear ~ mpg + hp + disp, data = mtcars)
Coefficients:
  (Intercept)       mpg         hp        disp
4   0.3892548 0.2707320 0.02227133 -0.04428756
5 -17.6837050 0.6115097 0.15511207 -0.08815984
Std. Errors:
  (Intercept)       mpg         hp       disp
4    17.30456 0.5917790 0.05813736 0.02735148
5    15.46373 0.5754793 0.08651377 0.06060359
Residual Deviance: 17.47162
AIC: 33.47162

如预期，摘要包括两个子模型（分别以4和5索引）的两组系数。

最后，让我们使用多项式逻辑回归模型进行预测，并计算混淆矩阵：

# make prediction
predicted_gears = predict(multinom_model, newdata = mtcars)
# view the confusion matrix
>>> table(Predicted = predicted_gears, Actual = mtcars$gear)
         Actual
Predicted  3  4  5
        3 14  0  0
        4  1 12  1
        5  0  0  4

结果表明，只有两个分类错误，分类性能相当不错。

摘要

在本章中，我们深入探讨了逻辑回归的世界，其理论基础及其实际应用。我们首先探讨了逻辑回归的基本结构及其与线性回归的比较。然后，我们介绍了 Sigmoid 变换的概念，这是逻辑回归中的一个关键元素，确保我们模型的输出被限制在0和1之间。这一部分帮助我们更好地理解逻辑回归在二元分类任务中的优势。

接下来，我们深入探讨了对数几率和几率比的概念，这两个是逻辑回归模型的关键组成部分。理解这些概念使我们能够理解模型预测的现实世界影响，并有效地解释其参数。随后章节介绍了 CEL，这是逻辑回归中使用的成本函数。具体来说，我们讨论了该损失函数如何确保我们的模型学会预测二元类别的准确概率。

当涉及到评估逻辑回归时，我们学习了各种指标，包括准确率、错误率、精确度、召回率、灵敏度、特异性和 AUC。这种理解将使我们能够准确评估我们的逻辑回归模型的表现。

重要的讨论围绕着处理不平衡数据集的问题，这是现实世界数据中常见的场景。我们理解了数据不平衡对我们模型的影响，并学习了如何通过重采样技术等策略有效地处理这种情况。此外，我们还讨论了惩罚逻辑回归，其中我们将 L1（lasso）或 L2（ridge）正则化纳入我们的逻辑回归模型。这种惩罚技术通过保持模型权重的大小，并在处理高维数据时创建更简单的模型，帮助我们防止过拟合。

最后，我们简要提到了多项逻辑回归，这是逻辑回归的扩展，用于多类分类问题。这部分内容为我们提供了处理目标变量包含两个以上类别的情况的见解。

在本章结束时，我们广泛了解了逻辑回归、其实现及其细微差别。这些知识为我们深入研究更复杂的分类方法和策略奠定了基础。

在下一章中，我们将介绍贝叶斯统计，这是统计建模的另一个主要分支。

第十四章：贝叶斯统计学

在本章中，我们将介绍贝叶斯推断框架，涵盖其核心组件和实现细节。贝叶斯推断引入了一个有用的框架，它提供了对目标结果预测的合理猜测以及量化的不确定性估计。从包含领域专业知识的先验分布开始，贝叶斯推断方法允许我们从数据中持续学习更新信息，并更新后验分布，以形成对潜在参数的更现实的观点。

在本章结束时，你将掌握使用贝叶斯推断框架的基本技能。你将学习贝叶斯定理的核心理论和它在贝叶斯线性回归模型中的应用。

本章我们将涵盖以下主要内容：

• 介绍贝叶斯统计学

• 深入了解贝叶斯推断

• 完整的贝叶斯推断过程

• 带有分类变量的贝叶斯线性回归

技术要求

要运行本章中的代码，你需要拥有以下软件包的最新版本：

• ggplot2, 3.4.0

• ggridges, 0.5.4

• rjags, 4.13

• coda, 0.19.4

请注意，前面提到的软件包版本是在我撰写本章时的最新版本。

本章的所有代码和数据均可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_14/working.R找到。

介绍贝叶斯统计学

贝叶斯方法在统计学和机器学习（ML）中提供了一个逻辑清晰、透明且可解释的框架。这是一个统一的框架，可以为统计推断和预测构建特定问题的模型。特别是，贝叶斯推断提供了一种方法，在已知事实（观察数据）的情况下，通过概率描述未知量的可能值的不确定性——即感兴趣的随机变量。

使用贝叶斯统计学，我们能够表达我们对未知量的先验假设，并根据观察数据调整这一假设。它提供了常见的统计程序（如假设检验和线性回归）的贝叶斯版本，这些内容在第十一章，统计估计和第十二章，R 中的线性回归中有所涉及。与迄今为止我们采用的频率主义方法相比，贝叶斯方法还允许我们构建特定问题的模型，以连续学习的方式（通过贝叶斯后验更新，将在下一节中介绍）充分利用数据。

例如，未知量对应于我们在线性或逻辑回归模型中试图估计的参数。我们不是将它们视为固定量，并使用最大似然原理来估计它们的值，而是贝叶斯方法将它们视为具有各自可能值概率分布的移动变量。

让我们首先了解一下著名的贝叶斯定理。

初探贝叶斯定理

贝叶斯定理描述了统计量条件概率之间的关系。在线性回归的背景下，我们会将参数 β 视为一个随机变量，而不是像我们在第十二章中处理线性回归那样将其视为固定量。例如，在一个简单的线性回归模型 y = βx 中，我们不会通过最小化给定数据 (x, y) 的普通最小二乘法（OLS）来获得单个最佳参数 β，而是将 β 视为一个遵循特定分布的随机变量。

做这件事涉及到关于 β 的两个分布，即感兴趣的参数，如下所示：

• 第一个分布是先验分布 P(β)，这对应于我们在观察任何数据之前分配给 β 的主观分布。这个分布封装了我们在观察任何实际数据之前对 β 可能值概率的先验信念。

• 第二个分布是后验分布 P(β|x, y)，这对应于观察数据后我们对该分布更新的信念。这是我们希望通过更新来估计的分布。这种更新是为了使我们的先验信念与我们在现实中实际观察到的相符合。

自然地，我们希望随着训练规模的增大，后验分布 P(β|x, y) 越来越接近数据反映的内容，并且相应地越来越远离先验信念。在这里，数据指的是 (x, y)。为了进行更新，数据将以我们所说的似然函数 P(y|x, β) 的形式进入，也被称为生成模型。也就是说，P(y|x, β) 表示在给定输入特征 x 和参数值 β 的情况下观察目标 y 的似然（类似于概率，尽管是非归一化的方式），我们将 β 视为一个特定的值而不是一个随机变量。换句话说，我们首先从分布 P(β) 中采样以获得 β 的具体值，然后根据特定的观察模型，在给定输入特征 x 的情况下获得实际数据点 y。例如，如果我们假设误差遵循正态分布，我们会假设观察模型为正态分布。

现在，我们准备通过以下贝叶斯定理将这些三个量一起编译：

P(β | x, y) =  P(y|x, β)P(β) _ P(y|x)

在这里，P(y|x) 被称为证据，它作为归一化常数，确保后验分布是一个有效的概率分布，这意味着每个概率都是非负的，并且总和（或积分）为 1。

图 14.1 展示了贝叶斯定理：

图 14.1 – 说明贝叶斯定理计算后验分布 P(𝜷 | x, y)

注意，贝叶斯线性回归中的先验分布 P(β) 可以选择来模拟我们对参数 β 的先验信念，这在使用基于 OLS 的线性回归的频率派框架时是不可用的。在实践中，我们通常采用正态先验分布，但它可以是任何在观察任何数据之前捕捉 β 可能值概率的先验信念的分布。因此，贝叶斯框架允许我们以原则性的方式将先验知识纳入建模中。例如，如果我们认为所有特征应该有相似的影响，我们就可以配置系数的先验分布以围绕相同的值为中心。

由于参数 β 遵循后验分布 P(β | x, y)，给定新的输入数据 x 的结果预测 y 将不会是一个单一的数字，正如频率派方法中的情况。相反，我们将获得一系列可能的 y 值，这些值遵循后验预测分布 P( y*| x*, β)。也就是说，由于参数 β 的随机性，预测 y 被视为一个随机变量。然后我们可以使用这个分布来理解结果预测的不确定性。例如，如果后验预测分布 P( y*| x*, β) 较宽，那么从 P( y*| x*, β) 中采样的结果预测将包含更高的不确定性。另一方面，如果分布较窄，结果预测将更加集中，因此更加自信。随着新数据的出现，后验分布 P(β | x, y) 也将继续演变。

下一个部分将介绍更多关于生成模型的内容。

理解生成模型

在贝叶斯推理中，生成模型指定了控制数据如何生成的概率分布。例如，当可用的目标数据是二进制时，我们可以假设它是按照参数 p 的伯努利分布生成的，其中 p 代表成功的概率。为了得到一系列的二进制结果，我们首先为 p 分配一个概率值，然后使用这个伯努利分布通过从这个分布中重复采样来生成二进制标签。

让我们通过一个练习来理解生成过程。

练习 14.1 – 生成二进制结果

在这个练习中，我们将根据伯努利分布生成一系列二进制结果。这涉及到将来自均匀分布的 0 到 1 之间的随机样本与预设的成功概率进行比较。遵循以下步骤：

1. 我们首先使用均匀分布生成一个介于 0 和 1 范围内的随机数：

   
   set.seed(1)
   random_prob = runif(1, min = 0, max = 1)
   >>> random_prob
   0.2655087
   

   再次提醒，为了可重复性，请设置随机种子。

2. 接下来，我们将数字与预设的概率 0.2 进行比较：

   
   prop_success = 0.2
   >>> random_prob < prop_success
   FALSE
   

   这完成了单个二进制数的生成。现在，让我们将其扩展到 10 个数字。

3. 在以下代码的帮助下，我们将生成 10 个二进制数：

   
   n_samples = 10
   data = c()
   for(sample_idx in 1:n_samples) {
     data[sample_idx] <- runif(1, min = 0, max = 1) < prop_success
   }
   >>> data
   FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
   

   在这里，我们使用 for 循环重复生成均匀随机数，将其与预设的成功概率进行比较，然后将结果存储在向量中。

4. 将向量转换为 0 和 1 数字，如下所示：

   
   data = as.numeric(data)
   >>> data
   0 0 0 0 0 0 0 0 1 0
   

   注意，尽管成功的概率为 20%，但我们观察到 10 次抽取中有 1 次成功。随着样本量的增加，我们预计成功的经验概率（本例中为 10%）将接近理论值（20%）。

结果表明，这个生成模型对应于二项过程或二项分布，这允许我们一次性生成二进制结果。具体来说，我们可以使用 rbinom() 函数来模拟二项分布的数据，如下面的代码片段所示：

set.seed(1)
>>> rbinom(n = n_samples, size = 1, prob = prop_success)

在这里，n 是从生成模型中生成的样本数量，size 是要运行的试验次数，prob 是介于 0.0 和 1.0 之间的成功基础概率。

注意，我们本质上是在处理一个已知的参数 p，即成功的概率。在实践中，这将是未知参数，是我们从数据中感兴趣估计的东西。贝叶斯推理将允许我们这样做，在下一节中介绍的先验分布和似然函数的帮助下。

理解先验分布

先验分布，贝叶斯推理的一个基本组成部分，表示在我们观察实际数据之前对潜在参数的先验知识或信念。它本质上基于特定领域的偏好或专业知识指定参数的概率分布。如果我们有合理的理由相信某些参数值更有可能，我们可以选择一个反映这种偏好的先验分布。

先验分布，表示为 P(β)，将β视为一个随机变量并指定其概率分布。也就是说，它告诉我们哪些β的值比其他值更有可能。在我们关于贝叶斯线性回归的运行示例中，β的先验分布通常选择为多元高斯分布。这主要是出于数学上的便利，因为高斯分布具有一些良好的性质，使得它更容易处理。然而，如果我们没有先验偏好，均匀分布（对所有可能的选择给予相同的概率，因此是无信息的）可能是一个好的候选者。

注意，我们也可以使用先验对模型施加一种正则化的形式。通过选择一个偏好较小参数值的先验分布（例如，以0为中心的高斯分布），我们可以阻止模型找到具有大系数的解，从而防止过拟合。

在下面的代码片段中，我们随机生成 10 个 P(β)样本，这些样本遵循在0和0.2之间的均匀分布：

set.seed(1)
prop_successes = runif(n_samples, min = 0.0, max = 0.2)
>>> prop_successes
0.05310173 0.07442478 0.11457067 0.18164156 0.04033639 0.17967794 0.18893505 0.13215956 0.12582281 0.01235725

当我们将成功的概率 p ∈ [0,1]建模为一个随机变量时，我们通常将贝塔分布作为先验分布。例如，在下面的代码片段中，我们使用rbeta()函数从贝塔分布 p ∼ Beta(35,55)中生成 1,000 个 p 的样本，然后将其转换为 DataFrame 格式后展示其密度图：

library(ggplot2)
# Sample 1000 draws from Beta(35,55) prior
prior_A = rbeta(n = 1000, shape1 = 35, shape2 = 55)
# Store the results in a data frame
prior_sim = data.frame(prior_A)
# Construct a density plot of the prior sample
ggplot(prior_sim, aes(x = prior_A)) +
  geom_density()

运行前面的代码将生成图 14.2所示的输出。它显示了成功概率在0和1范围内的0.4附近的先验集中：

图 14.2 – 可视化先验分布的密度图

Beta(a, b)分布定义在从0到1的区间上，因此为概率随机变量提供了一个自然且灵活的先验。我们可以调整Beta形状参数 a 和 b 以产生不同的先验模型。在下面的代码片段中，我们比较了原始Beta(35,55)先验分布与两种替代方案：Beta(1,1)和Beta(100,100)。然后我们将所有三个先验分布一起绘制出来：

# Sample draws from the Beta(1,1) prior
prior_B = rbeta(n = 1000, shape1 = 1, shape2 = 1)
# Sample draws from the Beta(100,100) prior
prior_C = rbeta(n = 1000, shape1 = 100, shape2 = 100)
# Combine the results in a single data frame
prior_all = data.frame(samples = c(prior_A, prior_B, prior_C),
                        priors = rep(c("A","B","C"), each = 1000))
# Plot the 3 priors
ggplot(prior_all, aes(x = samples, fill = priors)) +
  geom_density(alpha = 0.5)

前面的代码返回图 14.3所示的输出。每个先验分布都有一个不同的偏好区域。例如，分布 B 接近均匀分布，没有特定的偏好，而分布 C 对0.5有强烈的偏好，如峰值在0.5附近和狭窄的分布所指示：

图 14.3 – 可视化三种不同的先验分布

在下一节中，我们将向您介绍似然函数。

介绍似然函数

似然函数描述了在给定一组固定模型参数的情况下，观察到的数据有多可能。在参数模型（假设一组特定参数的模型）中，似然是观察数据的概率作为参数的函数。似然函数的具体形式取决于数据所假设的分布（更具体地说，是观测模型）。例如，如果我们假设数据遵循正态分布，似然函数将采用正态概率密度函数的形式。

让我们来看一个具体的例子。假设我们正在开发一个具有标准正态误差的简单线性回归模型，表示为 y = βx + ϵ，其中 ϵ ∼ N(0,1)。在这里，我们忽略了截距项，只考虑了斜率。对于特定的数据点 (x_i, y_i)，我们可以将似然 l_i 表示为在误差项的概率密度函数上的概率评估：

l_i = 1 / √2π e^(-(y_i−βx_i)² / 2)

由于数据集由总共 n 个输入输出对组成，所有数据点的联合似然可以表示如下：

L = Π i=1 n l_i = (1 / √2π)^n e^−∑ i=1 n (y_i−βx_i)² / 2

在实践中，我们通常会引入对数变换后处理对数似然，如下所示：

logL = log(1 / √2π)^n e^−∑ i=1 n (y_i−βx_i)² / 2 = −0.5nlog2π − 0.5∑ i=1 n (y_i − β x_i)²

与基于 OLS 的线性回归中使用的目标函数相比，我们发现指数项 ∑ i=1 n (y_i − β x_i)² 正好是平方误差的总和。当使用最大似然估计过程时，这两个不同的目标函数变得彼此等价。换句话说，我们有以下：

logL ≈ − ∑ i=1 n (y_i − β x_i)²

在最后一步，我们忽略了常数项 (1 / √2π)^n。

让我们通过一个例子来了解如何计算一组观察数据的联合似然。在下面的代码列表中，我们创建了一个包含数据点 x 和 y 的列表，以及一个具有系数 b 的简单线性回归模型，以生成预测值 y_pred，以及残差项 residuals：

# observed data
x = c(1, 2, 3, 4, 5)
y = c(2, 3, 5, 6, 7)
# parameter value
b = 0.8
# calculate the predicted values
y_pred = b * x
# calculate the residuals
residuals = y - y_pred

然后，我们可以将联合似然的闭式表达式代入，并计算总对数似然，如下所示：

log_likelihood = -0.5 * length(y) * log(2 * pi) - 0.5 * sum(residuals²)
log_likelihood
-18.09469

让我们看看二项模型的另一个例子。如前所述，当潜在参数 p 代表成功的概率时，我们可以使用 rbinom() 函数来获得在总抽取次数中观察到特定结果（成功次数）的概率。在下面的代码片段中，我们首先创建一个概率向量来表示不同的成功概率，并计算从二项模型 Bin(p, 10)（其中 10 是总抽取次数）中抽取的 1,000 次试验中的对应似然。最后，我们使用 ggridges 包中的 geom_density_ridges() 函数通过堆叠密度图可视化所有似然函数：

library(ggridges)
# Define a vector of 1000 p values
p_grid = seq(from = 0, to = 1, length.out = 1000)
# Simulate 10 trials for each p in p_grid, each trial has 1000 samples
sim_result = rbinom(n = 1000, size = 10, prob = p_grid)
# Collect results in a data frame
likelihood_sim = data.frame(p_grid, sim_result)
# Density plots of p_grid grouped by sim_result
ggplot(likelihood_sim, aes(x = p_grid, y = sim_result, group = sim_result)) +
  geom_density_ridges()

上述代码产生以下输出：

图 14.4 – 可视化堆叠密度图作为来自二项分布不同抽样的似然函数

下一个部分介绍了后验模型。

引入后验模型

后验分布代表了观察可用数据后我们对未知参数的了解。它结合了通过先验分布表达的前验信念和数据通过似然函数提供的证据，形成一个新分布，覆盖可能的参数值，这些值可以是离散的或连续的。

基于贝叶斯定理，后验分布与先验分布和似然函数的乘积成比例：

P(β | x, y) ∝ P(y|x, β)P(β)

注意，在求解参数 β 的最优值时，我们不需要知道分母中的证据项，因为 P(y|x) 与 β 完全独立。随着我们收集更多数据并更新我们的信念，后验分布通常会围绕真实参数值变得更加尖锐，这表明置信度有所提高。

然而，当我们需要知道 P(y|x) 以计算 P(β | x, y) 时，任务并不那么直接，因为可能没有封闭形式的解，或者参数是多维的，这阻止了直接计算嵌套积分。在这种情况下，我们通常会求助于数值方法，如马尔可夫链蒙特卡洛（MCMC）或近似方法，如变分推断。

让我们通过一个练习来了解整体推理过程。我们将使用 rjags 包根据我们的运行示例进行计算。

练习 14.2 – 获取后验分布

在这个练习中，我们将使用 rjags 包进行贝叶斯推理，包括指定模型架构和获取潜在参数的后验分布。按照以下步骤进行：

1. 我们首先定义似然函数为一个二项分布（使用dbin），参数p代表成功概率，n代表样本总数。我们还将定义参数p的先验分布为一个贝塔分布（使用dbeta），参数为a和b：

   
   library(rjags)
   # define the model
   bayes_model = "model{
       # Likelihood model for X
       X ~ dbin(p, n)
       # Prior model for p
       p ~ dbeta(a, b)
   }"
   Compile the model using the bayes.model() function.
   # compile the model
   bayes_jags = jags.model(textConnection(bayes_model),
                           data = list(a = 1, b = 1, X = 3, n = 10))
   

   这里，我们指定textConnection(bayes_model)将模型规范字符串传递给jags。数据参数是一个列表，指定了观察到的数据和先验分布的参数。

2. 接下来，我们从后验分布中抽取样本：

   
   # simulate the posterior
   bayes_sim = coda.samples(model = bayes_jags, variable.names = c("p"), n.iter = 10000)
   

   这里，使用coda.samples()函数运行 MCMC 模拟并抽取参数样本。n.iter参数指定了 MCMC 模拟的迭代次数。

3. 最后，我们绘制后验分布图：

   
   # plot the posterior
   plot(bayes_sim, trace = FALSE, xlim = c(0,1), ylim = c(0,3))
   

   运行前面的代码将产生图 14.5中显示的输出：

图 14.5 – 可视化基础参数（成功概率）的后验分布

在下一节中，我们将更深入地探讨贝叶斯推理，首先介绍正态-正态模型，这是贝叶斯推理中常用的一种模型类型。

深入贝叶斯推理

贝叶斯推理是一种统计方法，它利用条件概率来更新关于统计模型参数的先验信念，基于观察到的数据。贝叶斯推理的输出是一个后验分布，它表示我们在观察数据后对参数更新的信念。

当计算精确的后验分布很困难时，我们通常会求助于 MCMC，这是一种估计随机变量分布的技术。它是贝叶斯推理中从后验分布中生成样本的常用方法，尤其是在模型参数的维度很高时，使得解析解变得不可行。

在下一节中，我们将介绍正态-正态模型，并使用 MCMC 来估计其后验分布。

介绍正态-正态模型

正态-正态模型是贝叶斯推理中的另一个基础模型。它指的是当似然函数和先验分布都是正态分布，且都遵循钟形曲线的情况。这类模型在贝叶斯统计中经常被使用，因为它对后验分布有封闭形式的解，而这个解恰好也是正态分布。

让我们通过一个基于 MCMC 的贝叶斯推理的 normal-normal 模型的具体练习来观察。

练习 14.3 – 使用正态-正态模型

在这个练习中，我们将定义一个正态似然函数，其均值遵循正态先验分布，而标准差遵循均匀先验分布。然后我们将使用rjags来获取均值和标准差的后验估计。请按照以下步骤操作：

1. 让我们从模拟 100 个遵循具有真实均值2和真实标准差1的正态分布的数据点开始：

   
   library(coda)
   set.seed(1)
   mu_true = 2
   sd_true = 1
   n = 100
   data = rnorm(n, mean = mu_true, sd = sd_true)
   

2. 现在，我们指定模型架构，包括数据的一个正态似然函数、均值变量的正态先验分布和标准差变量的均匀先验。正态先验由0和0.1分别参数化均值和标准差，而均匀先验的范围从0到10：

   
   model_string = "model {
       for (i in 1:n) {
           y[i] ~ dnorm(mu, prec)
       }
       mu ~ dnorm(0, 0.1)
       sigma ~ dunif(0, 10)
       prec <- pow(sigma, -2)
   }"
   

   前面的模型表明，我们数据中的每个观测值（y[i]）遵循一个具有均值mu和精度prec的正态分布（定义为方差的倒数，因此prec <- pow(sigma, -2)）。均值mu遵循一个具有均值0和接近0的精度的正态分布，这对应于相对较大的方差，因此对mu的先验信念较弱。标准差sigma遵循一个从0到10的均匀分布，这表达了在这两个界限之间对其值的完全不确定性。

3. 接下来，在jags中编译模型并烧尽马尔可夫链，如下所示：

   
   data_jags = list(y = data, n = n)
   model = jags.model(textConnection(model_string), data = data_jags)
   update(model, 1000)  # burn-in
   

   在这里，burn-in期是一系列初始迭代次数，我们在进行贝叶斯推理时将其丢弃，假设在此期间产生的链可能没有收敛。想法是让马尔可夫链烧尽，直到它达到一个稳定且能反映我们感兴趣的后续分布的分布。在这种情况下，我们将丢弃马尔可夫链的前 1,000 次迭代，这些迭代在后续分析中不被使用。

4. 在这里，我们从后验分布中生成样本：

   
   params = c("mu", "sigma")
   samples = coda.samples(model, params, n.iter = 10000)
   # print summary statistics for the posterior samples
   >>> summary(samples)
   Iterations = 2001:12000
   Thinning interval = 1
   Number of chains = 1
   Sample size per chain = 10000
   1\. Empirical mean and standard deviation for each variable,
      plus standard error of the mean:
           Mean      SD  Naive SE Time-series SE
   mu    2.1066 0.09108 0.0009108      0.0009108
   sigma 0.9097 0.06496 0.0006496      0.0008996
   2\. Quantiles for each variable:
          2.5%    25%    50%    75% 97.5%
   mu    1.929 2.0453 2.1064 2.1666 2.287
   sigma 0.792 0.8652 0.9055 0.9513 1.046
   

   在这里，summary()函数为mu和sigma的后验样本提供了有用的汇总统计信息，包括它们的均值、中位数和可信区间。

5. 最后，通过运行以下命令来可视化结果：

   
   >>> plot(samples)
   

   运行前面的代码行将生成如图 14.6 所示的输出，该图显示了mu和sigma的后验样本的轨迹图和密度图：

图 14.6 – 可视化后验样本的轨迹图和密度图

在下一节中，我们将进一步讨论 MCMC。

介绍 MCMC

马尔可夫链是一个数学模型，在有限或可数的可能状态之间从一个状态过渡到另一个状态。它是一系列随机变量，其中未来的状态只取决于当前状态，而不取决于它之前的事件序列。这个特性被称为马尔可夫性质，或无记忆性。

在统计推断的背景下，我们可以从复杂的概率分布中进行采样，并通过蒙特卡洛模拟创建序列数据的模型，因此术语 MCMC。MCMC 算法构建参数值的马尔可夫链，其中链的平稳分布是潜在参数的后验分布。链是通过迭代提出新的参数值并根据预设规则接受或拒绝这些候选值来生成的，确保样本收敛到后验分布。

让我们分析之前 MCMC 链的细节。在下面的代码片段中，我们将链转换为 DataFrame 并打印其前几行：

# Store the chains in a data frame
mcmc_chains <- data.frame(samples[[1]], iter = 1:10000)
# Check out the head
>>> head(mcmc_chains)
        mu     sigma iter
1 2.159540 0.8678513    1
2 2.141280 0.8719263    2
3 1.975057 0.8568497    3
4 2.054670 0.9313297    4
5 2.144810 1.0349093    5
6 2.001104 1.0597861    6

这些是生成以近似 mu 和 sigma 的后验分布的 MCMC 样本。每个样本仅依赖于前一个样本，与其他先验样本无关。我们可以将这些样本绘制在称为迹图的线图中，如下代码片段所示：

# Use plot() to construct trace plots
>>> plot(samples, density = FALSE)

运行上述代码生成 图 14.7 中显示的输出。迹图是 MCMC 中常用的诊断工具。它是 MCMC 算法每次迭代或步骤中样本值的图形表示。在这种情况下，迹图没有明显的趋势，表明两条链都稳定收敛：

图 14.7 – 显示 MCMC 链的迹图

让我们使用 ggplot() 观察 mu 的前 100 个样本：

# Trace plot the first 100 iterations of the mu chain
>>> ggplot(mcmc_chains[1:100, ], aes(x = iter, y = mu)) +
  geom_line() +
  theme(axis.title.x = element_text(size = 20),  # Increase x-axis label size
        axis.title.y = element_text(size = 20))  # Increase y-axis label size

运行上述代码生成 图 14.8 中显示的输出。我们可以看到采样器在充分探索了前一个区域后移动到了另一个区域：

图 14.8 – 显示 mu 链的前 100 次迭代

注意，我们也可以仅通过设置 trace = FALSE 来显示密度图，如下所示：

# Use plot() to construct density plots
>>> plot(samples, trace = FALSE)

上述代码返回了 图 14.9 中显示的输出：

图 14.9 – 显示两条链的密度图

在实践中，我们通常会运行多个 MCMC 链。这允许我们检查所有链是否收敛到相同的分布，这是确保 MCMC 采样正确执行的关键步骤。由于 MCMC 中的每个链都是从不同的初始点开始的，运行多个链可以帮助检查后验分布是否依赖于起始值。当所有链无论初始点如何都收敛到相同的分布时，我们对 MCMC 过程的稳定性和鲁棒性有更高的信心。

以下代码片段在四个链上运行 MCMC：

model2 = jags.model(textConnection(model_string), data = data_jags, n.chains = 4)
# simulate the posterior
samples2 <- coda.samples(model = model2, variable.names = params, n.iter = 1000)

我们可以如下检查迹图：

# Construct trace plots
>>> plot(samples2, density = FALSE)

如 图 14.10 所示，所有四个链或多或少都收敛到了一个稳定状态：

图 14.10 – 可视化两个链的密度图

最后，我们可以检查 MCMC 样本的摘要，如下所示：

>>> summary(samples2)
Iterations = 1001:2000
Thinning interval = 1
Number of chains = 4
Sample size per chain = 1000
1\. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:
        Mean      SD Naive SE Time-series SE
mu    2.1052 0.09046  0.00143        0.00144
sigma 0.9089 0.06454  0.00102        0.00134
2\. Quantiles for each variable:
        2.5%    25%   50%    75% 97.5%
mu    1.9282 2.0456 2.104 2.1655 2.282
sigma 0.7952 0.8626 0.906 0.9522 1.041

在下一节中，我们将介绍完整的贝叶斯推理过程，包括量化后验不确定性和基于参数后验分布进行预测。

完整的贝叶斯推理过程

完整的贝叶斯推理首先需要指定模型架构，包括未知（未观测）参数的先验分布以及确定数据生成方式的似然函数。然后，我们可以执行 MCMC 来推断给定观测数据集的这些参数的后验分布。最后，我们可以使用后验分布来量化这些参数的不确定性，或者对新输入数据进行预测，并量化预测的不确定性。

以下练习使用 mtcars 数据集说明了这个过程。

练习 14.4 – 执行完整贝叶斯推理

在这个练习中，我们将使用单个特征和两个未知参数（intercept 和 slope）执行贝叶斯线性回归。该模型研究 mtcars 数据集中汽车重量 (wt) 和马力 (hp) 之间的关系。按照以下步骤进行：

1. 指定一个贝叶斯推理模型，其中每个目标 wt 被建模为服从正态分布的随机变量的实现。mean 参数是两个参数（a 为 intercept，b 为 slope）与相应输入特征的线性组合。a 和 b 都服从正态分布，而 variance 参数服从均匀分布。以下代码片段展示了代码：

   
   # load the necessary libraries
   library(rjags)
   library(coda)
   # define the model
   model = "model{
       # Define model for data Y[i]
       for(i in 1:length(Y)) {
         Y[i] ~ dnorm(m[i], s^(-2))
         m[i] <- a + b * X[i]
       }
       # Define the a, b, s priors
       a ~ dnorm(0, 0.5^(-2))
       b ~ dnorm(1, 0.5^(-2))
       s ~ dunif(0, 20)
   }"
   

2. 使用三个链编译模型，并控制随机种子以实现模型的可重复性：

   
   # compile the model
   model = jags.model(textConnection(model),
                       data = list(Y = mtcars$wt, X = mtcars$hp),
                       n.chains = 3,
                       inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 100))
   

3. 运行三次迭代的预热期：

   
   # burn-in
   update(model, 1000)
   

4. 参数后验分布的 MCMC 样本，如图所示：

   
   # generate MCMC samples
   samples = coda.samples(model, variable.names = c("a", "b", "s"), n.iter = 5000)
   

5. 创建 MCMC 样本的迹图以评估收敛性，如下所示：

   
   # check convergence using trace plot
   >>> plot(samples)
   

   运行前面的代码生成如图 14.11 所示的输出，表明所有三个参数都有良好的收敛性：

图 14.11 – 可视化所有三个参数的迹图

1. 计算参数 a 和 b 的后验分布的均值，并使用这些均值进行点预测，并在数据的散点图上绘制预测线，如下所示：

   
   # Get the posterior estimates
   posterior_estimates = summary(samples)
   # Calculate the mean for each parameter
   a_mean = posterior_estimates$statistics["a", "Mean"]
   b_mean = posterior_estimates$statistics["b", "Mean"]
   # Plot the prediction line
   ggplot(mtcars, aes(x = hp, y = wt)) +
     geom_point() +
     geom_abline(intercept = a_mean, slope = b_mean) +
     labs(title = "Bayesian Linear Regression",
          x = "Horsepower",
          y = "Weight") +
     theme(plot.title = element_text(hjust = 0.5))
   

   在这里，我们使用后验分布的均值作为参数值，对每个输入特征进行点预测。生成的输出如图 14.12 所示：

图 14.12 – 使用每个参数的后验分布的均值进行点预测

值得注意的是，与之前章节中介绍的频率主义线性回归相比，贝叶斯线性回归提供了量化的不确定性。这种不确定性以可信区间（credible intervals）的形式出现，与置信区间不同。具体来说，我们通过将参数视为随机变量，将数据视为固定量来获得可信区间，这在大多数情况下我们只能观察一次数据时更有意义。

1. 使用HPDinterval()函数计算a和b，并在直方图中绘制置信区间，如下所示：

   
   # Extract samples
   a_samples = as.matrix(samples[, "a"])
   b_samples = as.matrix(samples[, "b"])
   # Calculate credible intervals
   a_hpd = coda::HPDinterval(coda::as.mcmc(a_samples))
   b_hpd = coda::HPDinterval(coda::as.mcmc(b_samples))
   # Plot histograms and credible intervals
   par(mfrow=c(2,1))  # Create 2 subplots
   # Parameter a
   hist(a_samples, freq=FALSE, xlab="a", main="Posterior distribution of a", col="lightgray")
   abline(v=a_hpd[1,1], col="red", lwd=2)  # Lower limit of the credible interval
   abline(v=a_hpd[1,2], col="red", lwd=2)  # Upper limit of the credible interval
   # Parameter b
   hist(b_samples, freq=FALSE, xlab="b", main="Posterior distribution of b", col="lightgray")
   abline(v=b_hpd[1,1], col="red", lwd=2)  # Lower limit of the credible interval
   abline(v=b_hpd[1,2], col="red", lwd=2)  # Upper limit of the credible interval
   

   运行命令生成图 14.13所示的输出。在这里，我们默认计算 95%的可信区间：

图 14.13 – 可视化模型参数的后验分布和可信区间

1. 最后，对新输入值 120 进行后验预测，并基于参数的后验样本列表绘制后验预测分布，如下所示：

   
   # make posterior predictions
   # Obtain the mean values of the MCMC samples for each parameter
   a_mean = mean(samples[[1]][,"a"])
   b_mean = mean(samples[[1]][,"b"])
   # New input (e.g., horsepower = 120)
   new_input = 120
   # Prediction
   predicted_weight = a_mean + b_mean * new_input
   print(predicted_weight)
   # Predictive distribution
   predicted_weights = samples[[1]][,"a"] + samples[[1]][,"b"] * new_input
   # Plot the predictive distribution
   >>> hist(predicted_weights, breaks = 30, main = "Posterior predictive distribution", xlab = "Predicted weight")
   

   运行命令生成图 14.14所示的输出。此后的验预测分布捕捉了模型对预测的不确定性：

图 14.14 – 可视化新输入特征的后验预测分布

下一节将介绍使用分类输入变量的贝叶斯线性回归模型。

带有分类变量的贝叶斯线性回归

当预测变量是分类变量，例如二元特征时，我们会对每个对应的类别设置一个参数。以下练习演示了这样一个例子。

练习 14.5 – 使用分类变量进行贝叶斯推理

在这个练习中，我们将检查am（自动或手动变速，一个分类变量）与mpg（每加仑英里数，一个连续变量）之间的关系。我们将定义mpg的正态似然均值作为am的函数，每个am级别有不同的均值mu[i]。我们还将给mu一个正态先验，给标准差s一个均匀先验。按照以下步骤进行：

1. 按照以下方式指定上述模型架构：

   
   # define the model
   model = "model{
       # Define model for data Y[i]
       for(i in 1:length(Y)) {
         Y[i] ~ dnorm(mu[am[i]+1], s^(-2))
       }
       # Define the mu, s priors
       for(j in 1:2){
         mu[j] ~ dnorm(20, 10^(-2))
       }
       s ~ dunif(0, 20)
   }"
   

2. 编译模型，生成后验样本，并显示收敛图：

   
   # compile the model
   model = jags.model(textConnection(model),
                       data = list(Y = mtcars$mpg, am = mtcars$am),
                       n.chains = 3,
                       inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 100))
   # burn-in
   update(model, 1000)
   # generate MCMC samples
   samples = coda.samples(model, variable.names = c("mu", "s"), n.iter = 5000)
   # check convergence using trace plot
   >>> plot(samples)
   

   运行前面的代码生成图 14.15所示的输出，表明所有模型参数的收敛性良好：

图 14.15 – 可视化收敛图

1. 按照以下方式绘制数据集的分布以及分类变量两个级别的参数均值估计：

   
   # Get the posterior estimates
   posterior_estimates = summary(samples)
   # Calculate the mean for each parameter
   mu1_mean = posterior_estimates$statistics["mu[1]", "Mean"]
   mu2_mean = posterior_estimates$statistics["mu[2]", "Mean"]
   # Plot the prediction line
   ggplot(mtcars, aes(x = as.factor(am), y = mpg)) +
     geom_jitter(width = 0.2) +
     geom_hline(aes(yintercept = mu1_mean, color = "Automatic"), linetype = "dashed") +
     geom_hline(aes(yintercept = mu2_mean, color = "Manual"), linetype = "dashed") +
     scale_color_manual(name = "Transmission", values = c("Automatic" = "red", "Manual" = "blue")) +
     labs(title = "Bayesian Linear Regression",
          x = "Transmission (0 = automatic, 1 = manual)",
          y = "Miles Per Gallon (mpg)") +
     theme(plot.title = element_text(hjust = 0.5),
           legend.position = "bottom")
   

   上述代码将生成图 14.16中所示的结果：

图 14.16 – 可视化数据集的分布和均值估计

摘要

本章全面介绍了贝叶斯统计，从对贝叶斯定理的基本探讨开始。我们深入研究了其组成部分，从理解生成模型开始，这有助于我们模拟数据并检查参数变化如何影响数据生成过程。

然后，我们专注于理解先验分布，这是贝叶斯统计的一个基本组成部分，它代表了我们对于一个不确定参数的先验知识。随后，我们介绍了似然函数，这是一个统计函数，它决定了在给定的参数值下，一组观测值发生的可能性。

接下来，我们介绍了后验模型的概念。这结合了我们的先验分布和似然函数，给出一个新的概率分布，它代表了在看到数据后的更新信念。我们还探讨了更复杂的模型，例如正态-正态模型，其中似然和先验都是正态分布的。我们进一步通过 MCMC 方法研究了贝叶斯推理的机制，这是一种强大的工具，用于估计参数分布和进行预测。完整的贝叶斯推理过程的详细说明伴随着这一部分。

最后，我们讨论了具有分类变量的贝叶斯线性回归，这扩展了方法到包括分类预测因子的模型。

恭喜！您已成功到达本书的结尾，这是您奉献和努力的证明。这次旅程，希望对您来说富有成效，对我而言当然也是如此，标志着您在动态的统计学和机器学习领域持续探索中的一个重要里程碑。我很荣幸您选择这本书作为这次航行的伴侣，我相信它为您的未来追求奠定了坚实的基础。