机器学习的数据清理和探索-三-

机器学习的数据清理和探索（三）

原文：annas-archive.org/md5/25ad0fee8d118820a9d79ad1484952bd

译者：飞龙

协议：CC BY-NC-SA 4.0

第三部分 - 使用监督学习建模连续目标

本书最后十章介绍了广泛的各种机器学习算法，用于预测连续或分类目标，或者在没有目标的情况下。我们在本章探索连续目标的模型。

这些章节的一个持续主题是，找到最佳模型部分是关于平衡方差和偏差。当我们的模型对训练数据拟合得太好时，它们可能没有我们需要的那么具有可推广性。在这种情况下，它们可能具有低偏差但高方差。在这些章节中，我们检查的每个算法，我们都讨论了实现这种平衡的策略。这些策略从线性回归和支持向量回归模型的正则化，到 k 近邻中的 k 值，再到决策树的最大深度。

我们也有机会练习在第六章**，准备模型评估中使用的预处理、特征选择和模型评估策略。本节中我们讨论的每个算法都需要不同的预处理以获得最佳结果。例如，特征缩放对于支持向量回归很重要，但对于决策树回归通常不是。我们可能会使用多项式变换与线性回归模型，但在决策树中这也会是不必要的。我们在这部分的每一章中考虑这些选择。

本节包括以下章节：

• 第七章, 线性回归模型

• 第八章, 支持向量回归

• 第九章, K 近邻、决策树、随机森林和梯度提升回归

第七章：第七章: 线性回归模型

线性回归可能是最著名的机器学习算法，其起源至少可以追溯到 200 年前的统计学习。如果你在大学里学习了统计学、计量经济学或心理测量学课程，你很可能已经接触到了线性回归，即使你在机器学习在本科课程中教授之前就已经上了这门课。实际上，许多社会和物理现象都可以成功地被建模为预测变量线性组合的函数。尽管如此，线性回归对机器学习来说仍然非常有用，就像这些年来对统计学习一样，不过，在机器学习中，我们更关心预测而不是参数值。

假设我们的特征和目标具有某些特性，线性回归是建模连续目标的一个非常好的选择。在本章中，我们将讨论线性回归模型的假设，并使用与这些假设大部分一致的数据构建模型。然而，我们还将探索替代方法，例如非线性回归，当这些假设不成立时我们会使用它。我们将通过查看解决过拟合可能性的技术来结束本章，例如 lasso 回归。

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

• 关键概念

• 线性回归与梯度下降

• 使用经典线性回归

• 使用 lasso 回归

• 使用非线性回归

• 使用梯度下降进行回归

技术要求

在本章中，我们将坚持使用大多数 Python 科学发行版中可用的库——NumPy、pandas 和 scikit-learn。本章的代码可以在本书的 GitHub 仓库中找到，网址为github.com/PacktPublishing/Data-Cleaning-and-Exploration-with-Machine-Learning。

关键概念

经常从事预测建模的分析师通常会构建数十个，甚至数百个线性回归模型。如果你像我一样，在 20 世纪 80 年代末为一家大型会计师事务所工作，并且从事预测工作，你可能每天都会花整天时间指定线性模型。你会运行所有可能的独立变量排列和因变量变换，并勤奋地寻找异方差性（残差中的非恒定方差）或多重共线性（高度相关的特征）的证据。但最重要的是，你努力识别关键预测变量，并解决任何参数估计（你的系数或权重）中的偏差。

线性回归模型的关键假设

那些努力的大部分至今仍然适用，尽管现在对预测准确性的重视程度超过了参数估计。我们现在担心过度拟合，而 30 年前我们并没有这样做。当线性回归模型的假设被违反时，我们也更有可能寻求替代方案。这些假设如下：

• 特征（自变量）与目标（因变量）之间存在线性关系

• 偶然误差（实际值与预测值之间的差异）是正态分布的

• 偶然误差在观测之间是独立的

• 偶然误差的方差是恒定的

在现实世界中，这些假设中有一个或多个被违反并不罕见。特征与目标之间的关系通常是线性的。特征的影响可能在该特征的范围内变化。任何熟悉“厨房里的人太多”这个表达的人可能都会欣赏到，第五个厨师带来的边际生产力的增加可能不如第二个或第三个厨师那么大。

我们的偶然误差有时不是正态分布的。这表明我们的模型在目标的一些范围内可能不够准确。例如，在目标范围的中间部分（比如第 25 到第 75 百分位数）可能会有较小的偶然误差，而在两端可能会有较大的偶然误差。这可能是由于与目标的关系是非线性的。

偶然误差可能不独立的原因有几个。这在时间序列数据中通常是这种情况。对于一个日股价模型，相邻天的偶然误差可能是相关的。这被称为自相关。这也可能是纵向或重复测量数据的问题。例如，我们可能有 20 个不同教室中 600 名学生的考试成绩，或者 100 人的年度工资收入。如果我们的模型未能考虑到某些特征在群体中不存在变化——在这些例子中是教室决定和个人决定的特征——那么我们的偶然误差就不会是独立的。

最后，我们的偶然误差在不同特征的不同范围内可能会有更大的变异性。如果我们正在预测全球气象站的温度，并且纬度是我们使用的特征之一，那么在较高的纬度值上可能会有更大的偶然误差。这被称为异方差性。这也可能表明我们的模型遗漏了重要的预测因子。

除了这四个关键假设之外，线性回归的另一个常见挑战是特征之间的高度相关性。这被称为多重共线性。正如我们在第五章中讨论的[特征选择]，当我们的模型难以隔离特定特征的独立影响，因为它与另一个特征变化很大时，我们可能会增加过拟合的风险。对于那些花费数周时间构建模型的人来说，这将是熟悉的，其中系数随着每个新规格的提出而大幅变动。

当违反这些假设中的一个或多个时，我们仍然可以使用传统的回归模型。然而，我们可能需要以某种方式转换数据。我们将在本章中讨论识别这些假设违反的技术、这些违反对模型性能的影响以及解决这些问题的可能方法。

线性回归和普通最小二乘法

线性回归中最常见的估计技术是普通最小二乘法（OLS）。OLS 选择系数，使得实际目标值与预测值之间平方距离之和最小。更精确地说，OLS 最小化以下：

在这里，是第 i 个观测的实际值，是预测值。正如我们讨论过的，实际目标值与预测目标值之间的差异，，被称为残差。

从图形上看，普通最小二乘法（OLS）通过我们的数据拟合一条线，使得数据点到这条线的垂直距离最小。以下图表展示了一个具有一个特征（称为简单线性回归）的模型，使用的是虚构的数据点。每个数据点到回归线的垂直距离是残差，可以是正数或负数：

图 7.1 – 普通最小二乘法回归线

这条线，，给出了每个x值的y预测值。它等于估计的截距，，加上特征估计系数乘以特征值，。这就是 OLS 线。任何其他穿过数据的直线都会导致平方残差之和更高。这可以扩展到多线性回归模型 – 即具有一个以上特征的模型：

在这里，y 是目标，每个 x 是一个特征，每个 是一个系数（或截距），n 是特征的数量，ɛ 是一个误差项。每个系数是从关联特征 1 单位变化引起的目标估计变化。这是一个很好的地方来注意到系数在整个特征的范围上是恒定的；也就是说，特征从 0 到 1 的增加被认为对目标的影响与从 999 到 1000 的影响相同。然而，这并不总是有道理。在本章的后面，我们将讨论当特征与目标之间的关系不是线性的情况下如何使用变换。

线性回归的一个重要优点是它不像其他监督回归算法那样计算量大。当线性回归表现良好，基于我们之前章节讨论的指标时，它是一个好的选择。这尤其适用于你有大量数据需要训练或你的业务流程不允许大量时间用于模型训练的情况。算法的效率还可以使其使用更资源密集的特征选择技术成为可能，例如我们讨论过的包装方法，这在第五章 特征选择中提到。正如我们看到的，你可能不想使用决策树回归器的穷举特征选择。然而，对于线性回归模型来说，这可能是完全可行的。

线性回归和梯度下降

我们可以使用梯度下降而不是普通最小二乘法来估计我们的线性回归参数。梯度下降通过迭代可能的系数值来找到那些使残差平方和最小的值。它从随机的系数值开始，并计算该迭代的平方误差总和。然后，它为系数生成新的值，这些值比上一步的残差更小。当我们使用梯度下降时，我们指定一个学习率。学习率决定了每一步残差改进的量。

当处理非常大的数据集时，梯度下降通常是一个不错的选择。如果整个数据集无法装入你的机器内存，它可能就是唯一的选择。在下一节中，我们将使用 OLS 和梯度下降来估计我们的参数。

使用经典线性回归

在本节中，我们将指定一个相当直接的线性模型。我们将使用它来根据几个国家的经济和政治指标预测隐含的汽油税。但在我们指定模型之前，我们需要完成本书前几章讨论的预处理任务。

对我们的回归模型进行数据预处理

在本章中，我们将使用管道来预处理我们的数据，并在本书的其余部分也是如此。我们需要在数据缺失的地方填充值，识别和处理异常值，以及编码和缩放我们的数据。我们还需要以避免数据泄露并清理训练数据而不提前查看测试数据的方式来做这些。正如我们在第六章中看到的，准备模型评估，scikit-learn 的管道可以帮助我们完成这些任务。

我们将使用的数据集包含每个国家的隐含汽油税和一些可能的预测因子，包括人均国民收入、政府债务、燃料收入依赖性、汽车使用范围以及民主过程和政府有效性的衡量指标。

注意

这个关于各国隐含汽油税的数据集可以在哈佛数据共享平台dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/RX4JGK上供公众使用。它由Paasha Mahdavi、Cesar B. Martinez-Alvarez和Michael L. Ross编制。隐含汽油税是根据每升汽油的世界基准价格和当地价格之间的差异来计算的。当地价格高于基准价格表示征税。当基准价格更高时，它可以被视为补贴。我们将使用每个国家的 2014 年数据进行分析。

让我们先对数据进行预处理：

1. 首先，我们加载了许多我们在上一章中使用的库。然而，我们还需要两个新的库来构建我们的数据管道——ColumnTransformer和TransformedTargetRegressor。这些库允许我们构建一个对数值和分类特征进行不同预处理的管道，并且还可以转换我们的目标：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.linear_model import LinearRegression
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.compose import ColumnTransformer
   from sklearn.compose import TransformedTargetRegressor
   from sklearn.feature_selection import RFE
   from sklearn.impute import KNNImputer
   from sklearn.model_selection import cross_validate, KFold
   import sklearn.metrics as skmet
   import matplotlib.pyplot as plt
   

2. 我们可以通过添加自己的类来扩展 scikit-learn 管道的功能。让我们添加一个名为OutlierTrans的类来处理极端值。

要将这个类包含在管道中，它必须继承自BaseEstimator类。我们还必须继承自TransformerMixin，尽管还有其他可能性。我们的类需要fit和transform方法。我们可以在transform方法中放置将极端值赋为缺失值的代码。

但在我们能够使用我们的类之前，我们需要导入它。为了导入它，我们需要追加helperfunctions子文件夹，因为那里是我们放置包含我们的类的preprocfunc模块的地方：

import os
import sys
sys.path.append(os.getcwd() + “/helperfunctions”)
from preprocfunc import OutlierTrans

这导入了OutlierTrans类，我们可以将其添加到我们创建的管道中：

class OutlierTrans(BaseEstimator,TransformerMixin):
  def __init__(self,threshold=1.5):
    self.threshold = threshold

  def fit(self,X,y=None):
    return self

  def transform(self,X,y=None):
    Xnew = X.copy()
    for col in Xnew.columns:
      thirdq, firstq = Xnew[col].quantile(0.75),\
        Xnew[col].quantile(0.25)
      interquartilerange = self.threshold*(thirdq-firstq)
      outlierhigh, outlierlow = interquartilerange+thirdq,\
        firstq-interquartilerange
      Xnew.loc[(Xnew[col]>outlierhigh) | \
        (Xnew[col]<outlierlow),col] = np.nan
    return Xnew.values

OutlierTrans类使用一个相当标准的单变量方法来识别异常值。它计算每个特征的四分位数范围（IQR），然后将任何高于第三四分位数 1.5 倍或低于第一四分位数 1.5 倍以上的值设置为缺失。如果我们想更保守一些，可以将阈值改为其他值，例如 2.0。（我们在第一章，检查特征和目标的分布中讨论了识别异常值的技术。）

1. 接下来，我们加载 2014 年的汽油税数据。有 154 行——数据框中每个国家一行。一些特征有一些缺失值，但只有motorization_rate有两位数的缺失。motorization_rate是每人的汽车数量：

   fftaxrate14 = pd.read_csv(“data/fossilfueltaxrate14.csv”)
   fftaxrate14.set_index(‘countrycode’, inplace=True)
   fftaxrate14.info()
   <class ‘pandas.core.frame.DataFrame’>
   Index: 154 entries, AFG to ZWE
   Data columns (total 19 columns):
   #   Column                     Non-Null Count  Dtype
                                  --------------  -----  
   0   country                    154 non-null    object
   1   region                     154 non-null    object
   2   region_wb                  154 non-null    object
   3   year                       154 non-null    int64
   4   gas_tax_imp                154 non-null    float64
   5   bmgap_diesel_spotprice_la  146 non-null    float64
   6   fuel_income_dependence     152 non-null    float64
   7   national_income_per_cap    152 non-null    float64
   8   VAT_Rate                   151 non-null    float64
   9   gov_debt_per_gdp           139 non-null    float64
   10  polity                     151 non-null    float64
   11  democracy_polity           151 non-null    float64
   12  autocracy_polity           151 non-null    float64
   13  goveffect                  154 non-null    float64
   14  democracy_index            152 non-null    float64
   15  democracy                  154 non-null    int64
   16  nat_oil_comp               152 non-null    float64
   17  nat_oil_comp_state         152 non-null    float64
   18  motorization_rate          127 non-null    float64
   dtypes: float64(14), int64(2), object(3)
   memory usage: 24.1+ KB
   

2. 让我们将特征分为数值特征和二进制特征。我们将motorization_rate放入一个特殊类别，因为我们预计需要比其他特征做更多的工作：

   num_cols = [‘fuel_income_dependence’, 
     ’national_income_per_cap’, ‘VAT_Rate’,  
     ‘gov_debt_per_gdp’, ’polity’, ’goveffect’,
     ‘democracy_index’]
   dummy_cols = ‘democracy_polity’,’autocracy_polity’,
     ‘democracy’,’nat_oil_comp’,’nat_oil_comp_state’]
   spec_cols = [‘motorization_rate’]
   

3. 我们应该查看一些数值特征和目标的描述性统计。我们的目标gas_tax_imp的中位数为 0.52。请注意，一些特征的范围非常不同。超过一半的国家polity得分为 7 或更高；10 是可能的最高polity得分，意味着最民主。大多数国家的政府有效性值为负。democracy_index与polity是非常相似的一个度量，尽管有更多的变化：

   fftaxrate14[[‘gas_tax_imp’] + num_cols + spec_cols].\
     agg([‘count’,’min’,’median’,’max’]).T
                         count min    median   max
   gas_tax_imp             154 -0.80  0.52     1.73
   fuel_income_dependence  152 0.00   0.14     34.43
   national_income_per_cap 152 260.00 6,050.00 104,540.00
   VAT_Rate                151 0.00   16.50    27.00
   gov_debt_per_gdp        139 0.55   39.30    194.76
   polity                  151 -10.00 7.00     10.00
   goveffect               154 -2.04  -0.15    2.18
   democracy_index         152 0.03   0.57     0.93
   motorization_rate       127 0.00   0.20     0.81
   

4. 让我们再看看二进制特征的分布。我们必须将normalize设置为True以生成比率而不是计数。democracy_polity和autocracy_polity特征是polity特征的二进制版本；非常高的polity得分得到democracy_polity值为 1，而非常低的polity得分得到autocracy_polity值为 1。同样，democracy是一个虚拟特征，用于那些democracy_index值较高的国家。有趣的是，近一半的国家（0.46）拥有国家石油公司，而几乎四分之一（0.23）拥有国有国家石油公司：

   fftaxrate14[dummy_cols].apply(pd.value_counts, normalize=True).T
                                  0               1
   democracy_polity               0.41            0.59
   autocracy_polity               0.89            0.11
   democracy                      0.42            0.58
   nat_oil_comp                   0.54            0.46
   nat_oil_comp_state             0.77            0.23
   

所有这些都看起来相当不错。然而，我们需要对几个特征的缺失值做一些工作。我们还需要做一些缩放，但不需要进行任何编码，因为我们可以直接使用二进制特征。一些特征是相关的，因此我们需要进行一些特征消除。

1. 我们通过创建训练和测试数据框开始预处理。我们将只为测试保留 20%：

   target = fftaxrate14[[‘gas_tax_imp’]]
   features = fftaxrate14[num_cols + dummy_cols + spec_cols]
   X_train, X_test, y_train, y_test =  \
     train_test_split(features,\
     target, test_size=0.2, random_state=0)
   

2. 我们需要构建一个包含列转换的管道，以便我们可以对数值数据和分类数据进行不同的预处理。我们将为num_cols中的所有数值列构建一个管道，命名为standtrans。首先，我们想要将异常值设置为缺失值。我们将异常值定义为超过第三四分位数两倍以上的值，或者低于第一四分位数的值。我们将使用SimpleImputer将缺失值设置为该特征的中间值。

我们不希望对dummy_cols中的二元特征进行缩放，但我们确实想要使用SimpleImputer将每个分类特征的缺失值设置为最频繁的值。

我们不会为motorization_rate使用SimpleImputer。记住，motorization_rate不在num_cols列表中——它在spec_cols列表中。我们为spec_cols中的特征设置了一个特殊的管道，spectrans。我们将使用motorization_rate值：

standtrans = make_pipeline(OutlierTrans(2), 
  SimpleImputer(strategy=”median”), StandardScaler())
cattrans = make_pipeline(
  SimpleImputer(strategy=”most_frequent”))
spectrans = make_pipeline(OutlierTrans(2), 
  StandardScaler())
coltrans = ColumnTransformer(
  transformers=[
    (“stand”, standtrans, num_cols),
    (“cat”, cattrans, dummy_cols),
    (“spec”, spectrans, spec_cols)
  ]
)

这设置了我们在汽油税数据上想要进行的所有预处理。要进行转换，我们只需要调用列转换器的fit方法。然而，我们目前不会这样做，因为我们还想要将特征选择添加到管道中，并使其运行线性回归。我们将在接下来的几个步骤中完成这些操作。

运行和评估我们的线性模型

我们将使用递归特征消除（RFE）来选择模型的特征。RFE 具有包装特征选择方法的优点——它根据选定的算法评估特征，并在评估中考虑多元关系。然而，它也可能在计算上很昂贵。由于我们没有很多特征或观测值，在这种情况下这并不是一个大问题。

在选择特征后，我们运行线性回归模型并查看我们的预测。让我们开始吧：

1. 首先，我们创建线性回归和递归特征消除实例，并将它们添加到管道中。我们还创建了一个TransformedTargetRegressor对象，因为我们仍然需要转换目标。我们将我们的管道传递给TransformedTargetRegressor的回归器参数。

现在，我们可以调用目标回归器的fit方法。之后，管道的rfe步骤的support_属性将给我们提供选定的特征。同样，我们可以通过获取linearregression步骤的coef_值来获取系数。关键在于通过引用ttr.regressor我们能够访问到管道：

lr = LinearRegression()
rfe = RFE(estimator=lr, n_features_to_select=7)
pipe1 = make_pipeline(coltrans, 
  KNNImputer(n_neighbors=5), rfe, lr)
ttr=TransformedTargetRegressor(regressor=pipe1,
  transformer=StandardScaler())
ttr.fit(X_train, y_train)
selcols = X_train.columns[ttr.regressor_.named_steps[‘rfe’].support_]
coefs = ttr.regressor_.named_steps[‘linearregression’].coef_
np.column_stack((coefs.ravel(),selcols))
array([[0.44753064726665703, ‘VAT_Rate’],
       [0.12368913577287821, ‘gov_debt_per_gdp’],
       [0.17926454403985687, ‘goveffect’],
       [-0.22100930246392841, ‘autocracy_polity’],
       [-0.15726572731003752, ‘nat_oil_comp’],
       [-0.7013454686632653, ‘nat_oil_comp_state’],
       [0.13855012574945422, ‘motorization_rate’]], dtype=object)

我们的特征选择确定了增值税率、政府债务、政府有效性指标（goveffect）、国家是否属于威权主义类别、是否有国家石油公司以及是否为国有企业，以及机动化率作为前七个特征。特征的数目是一个超参数的例子，我们在这里选择七个是相当随意的。我们将在下一节讨论超参数调整的技术。

注意到在数据集中的几个威权/民主措施中，似乎最重要的是威权虚拟变量，对于polity得分非常低的国家的值为 1。它被估计对汽油税有负面影响；也就是说，会减少它们。

1. 让我们来看看预测值和残差。我们可以将测试数据中的特征传递给 transformer 的/pipeline 的predict方法来生成预测值。存在一点正偏斜和整体偏差；残差总体上是负的：

   pred = ttr.predict(X_test)
   preddf = pd.DataFrame(pred, columns=[‘prediction’],
     index=X_test.index).join(X_test).join(y_test)
   preddf[‘resid’] = preddf.gas_tax_imp-preddf.prediction
   preddf.resid.agg([‘mean’,’median’,’skew’,’kurtosis’])
   mean                -0.09
   median              -0.13
   skew                 0.61
   kurtosis             0.04
   Name: resid, dtype: float64
   

2. 让我们也生成一些整体模型评估统计信息。我们得到平均绝对误差为0.23。考虑到汽油税价格的中间值为0.52，这不是一个很好的平均误差。然而，r-squared 值还不错：

   print(“Mean Absolute Error: %.2f, R-squared: %.2f” % 
     (skmet.mean_absolute_error(y_test, pred),
     skmet.r2_score(y_test, pred)))
   Mean Absolute Error: 0.23, R-squared: 0.75
   

3. 通常查看残差的图表是有帮助的。让我们也画一条代表残差平均值的红色虚线：

   plt.hist(preddf.resid, color=”blue”, bins=np.arange(-0.5,1.0,0.25))
   plt.axvline(preddf.resid.mean(), color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Histogram of Residuals for Gax Tax Model”)
   plt.xlabel(“Residuals”)
   plt.ylabel(“Frequency”)
   plt.xlim()
   plt.show()
   

这产生了以下图表：

图 7.2 – 汽油税模型残差

这个图表显示了正偏斜。此外，我们的模型更有可能高估汽油税而不是低估它。（当预测值大于实际目标值时，残差为负。）

1. 让我们也看看预测值与残差的散点图。让我们在Y轴上画一条代表 0 的红色虚线：

   plt.scatter(preddf.prediction, preddf.resid, color=”blue”)
   plt.axhline(0, color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Scatterplot of Predictions and Residuals”)
   plt.xlabel(“Predicted Gax Tax”)
   plt.ylabel(“Residuals”)
   plt.show()
   

这产生了以下图表：

图 7.3 – 预测值与残差的散点图

在这里，高估发生在预测值的整个范围内，但没有低估（正残差）出现在预测值低于 0 或高于 1 的情况下。这应该让我们对我们的线性假设产生一些怀疑。

提高我们的模型评估

我们到目前为止评估模型的一个问题是，我们没有充分利用数据。我们只在大约 80%的数据上进行了训练。我们的指标也相当依赖于测试数据是否代表我们想要预测的真实世界。然而，可能并不总是如此。正如前一章所讨论的，我们可以通过 k 折交叉验证来提高我们的机会。

由于我们一直在使用 pipelines 进行我们的分析，我们已经为 k 折交叉验证做了很多工作。回想一下前一章，k 折模型评估将我们的数据分成 k 个相等的部分。其中一部分被指定为测试，其余部分用于训练。这重复 k 次，每次使用不同的折进行测试。

让我们尝试使用我们的线性回归模型进行 k 折交叉验证：

1. 我们将首先创建新的训练和测试 DataFrames，留下 10%用于后续验证。虽然保留数据用于验证不是必须的，但保留一小部分数据总是一个好主意：

   X_train, X_test, y_train, y_test =  \
     train_test_split(features,\
     target, test_size=0.1, random_state=1)
   

2. 我们还需要实例化KFold和LinearRegression对象：

   kf = KFold(n_splits=3, shuffle=True, random_state=0)
   

3. 现在，我们已经准备好运行我们的 k 折交叉验证。我们表示我们想要每个分割的 r-squared 和平均绝对误差。"cross_validate"自动为我们提供每个折叠的拟合和评分时间：

   scores = cross_validate(ttr, X=X_train, y=y_train,
     cv=kf, scoring=(‘r2’, ‘neg_mean_absolute_error’), n_jobs=1)
   print(“Mean Absolute Error: %.2f, R-squared: %.2f” % 
     (scores[‘test_neg_mean_absolute_error’].mean(),
     scores[‘test_r2’].mean()))
   Mean Absolute Error: -0.25, R-squared: 0.62
   

这些分数并不十分令人印象深刻。我们没有解释掉我们希望解释的那么多方差。R-squared 分数在三个折叠中平均约为 0.62。这部分的理由是每个折叠的测试 DataFrame 相当小，每个大约有 40 个观测值。尽管如此，我们应该探索对经典线性回归方法的修改，例如正则化和非线性回归。

正则化的一个优点是，我们可能在不需要经过计算成本高昂的特征选择过程的情况下获得类似的结果。正则化还可以帮助我们避免过拟合。在下一节中，我们将使用相同的数据探索 lasso 回归。我们还将研究非线性回归策略。

使用 lasso 回归

OLS 的关键特征是它以最小的偏差产生参数估计。然而，OLS 估计可能比我们想要的方差更高。当我们使用经典线性回归模型时，我们需要小心过拟合。减少过拟合可能性的一个策略是使用正则化。正则化还可能允许我们将特征选择和模型训练结合起来。这对于具有大量特征或观测值的数据集可能很重要。

与 OLS 最小化均方误差不同，正则化技术寻求最小误差和减少特征数量。我们在本节中探讨的 lasso 回归使用 L1 正则化，它惩罚系数的绝对值。岭回归类似。它使用 L2 正则化，惩罚系数的平方值。弹性网络回归使用 L1 和 L2 正则化。

再次强调，我们将使用上一节中的汽油税数据：

1. 我们将首先导入与上一节相同的库，除了我们将导入Lasso模块而不是linearregression模块：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.linear_model import Lasso
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.compose import ColumnTransformer
   from sklearn.compose import TransformedTargetRegressor
   from sklearn.model_selection import cross_validate, KFold
   import sklearn.metrics as skmet
   import matplotlib.pyplot as plt
   

2. 我们还需要我们创建的OutlierTrans类：

   import os
   import sys
   sys.path.append(os.getcwd() + “/helperfunctions”)
   from preprocfunc import OutlierTrans
   

3. 现在，让我们加载汽油税数据并创建测试和训练 DataFrame：

   fftaxrate14 = pd.read_csv(“data/fossilfueltaxrate14.csv”)
   fftaxrate14.set_index(‘countrycode’, inplace=True)
   num_cols = [‘fuel_income_dependence’,’national_income_per_cap’,
     ‘VAT_Rate’,  ‘gov_debt_per_gdp’,’polity’,’goveffect’,
     ‘democracy_index’]
   dummy_cols = [‘democracy_polity’,’autocracy_polity’,
   ‘democracy’,’nat_oil_comp’,’nat_oil_comp_state’]
   spec_cols = [‘motorization_rate’]
   target = fftaxrate14[[‘gas_tax_imp’]]
   features = fftaxrate14[num_cols + dummy_cols + spec_cols]
   X_train, X_test, y_train, y_test =  \
     train_test_split(features,\
     target, test_size=0.2, random_state=0)
   

4. 我们还需要设置列转换：

   standtrans = make_pipeline(
     OutlierTrans(2), SimpleImputer(strategy=”median”),
     StandardScaler())
   cattrans = make_pipeline(SimpleImputer(strategy=”most_frequent”))
   spectrans = make_pipeline(OutlierTrans(2), StandardScaler())
   coltrans = ColumnTransformer(
     transformers=[
       (“stand”, standtrans, num_cols),
       (“cat”, cattrans, dummy_cols),
       (“spec”, spectrans, spec_cols)
     ]
   )
   

5. 现在，我们已经准备好拟合我们的模型。我们将从一个相当保守的 alpha 值 0.1 开始。alpha 值越高，我们的系数惩罚就越大。在 0 时，我们得到与线性回归相同的结果。除了列转换和 lasso 回归之外，我们的管道还使用 KNN 插补缺失值。我们还将使用目标转换器来缩放汽油税目标。在我们拟合之前，我们将刚刚创建的管道传递给目标转换器的回归器参数：

   lasso = Lasso(alpha=0.1,fit_intercept=False)
   pipe1 = make_pipeline(coltrans, KNNImputer(n_neighbors=5), lasso)
   ttr=TransformedTargetRegressor(regressor=pipe1,transformer=StandardScaler())
   ttr.fit(X_train, y_train)
   

6. 让我们看看 lasso 回归的系数。如果我们将它们与上一节中线性回归的系数进行比较，我们会注意到我们最终选择了相同的特征。那些在递归特征选择中被消除的特征，在很大程度上与 lasso 回归中得到接近零值的特征相同：

   coefs = ttr.regressor_[‘lasso’].coef_
   np.column_stack((coefs.ravel(), num_cols + dummy_cols + spec_cols))
   array([[‘-0.0026505240129231175’, ‘fuel_income_dependence’],
          [‘0.0’, ‘national_income_per_cap’],
          [‘0.43472262042825915’, ‘VAT_Rate’],
          [‘0.10927136643326674’, ‘gov_debt_per_gdp’],
          [‘0.006825858127837494’, ‘polity’],
          [‘0.15823493727828816’, ‘goveffect’],
          [‘0.09622123660935211’, ‘democracy_index’],
          [‘0.0’, ‘democracy_polity’],
          [‘-0.0’, ‘autocracy_polity’],
          [‘0.0’, ‘democracy’],
          [‘-0.0’, ‘nat_oil_comp’],
          [‘-0.2199638245781246’, ‘nat_oil_comp_state’],
          [‘0.016680304258453165’, ‘motorization_rate’]], dtype=’<U32’)
   

7. 让我们看看这个模型的预测值和残差。残差看起来相当不错，几乎没有偏差，也没有很大的偏斜：

   pred = ttr.predict(X_test)
   preddf = pd.DataFrame(pred, columns=[‘prediction’],
     index=X_test.index).join(X_test).join(y_test)
   preddf[‘resid’] = preddf.gas_tax_imp-preddf.prediction
   preddf.resid.agg([‘mean’,’median’,’skew’,’kurtosis’])
   mean                 -0.06
   median               -0.07
   skew                  0.33
   kurtosis              0.10
   Name: resid, dtype: float64
   

8. 让我们也生成平均绝对误差和 r 平方。这些分数并不令人印象深刻。r 平方低于线性回归，但平均绝对误差大致相同：

   print(“Mean Absolute Error: %.2f, R-squared: %.2f” % 
     (skmet.mean_absolute_error(y_test, pred),
     skmet.r2_score(y_test, pred)))
   Mean Absolute Error: 0.24, R-squared: 0.68
   

9. 我们应该查看残差直方图。残差的分布与线性回归模型相当相似：

   plt.hist(preddf.resid, color=”blue”, bins=np.arange(-0.5,1.0,0.25))
   plt.axvline(preddf.resid.mean(), color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Histogram of Residuals for Gax Tax Model”)
   plt.xlabel(“Residuals”)
   plt.ylabel(“Frequency”)
   plt.show()
   

这产生了以下图表：

图 7.4 – 汽油税模型残差

1. 让我们也看看预测值与残差的散点图。我们的模型可能在较低范围内高估，在较高范围内低估。这与线性模型不同，线性模型在两个极端都持续高估：

   plt.scatter(preddf.prediction, preddf.resid, color=”blue”)
   plt.axhline(0, color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Scatterplot of Predictions and Residuals”)
   plt.xlabel(“Predicted Gax Tax”)
   plt.ylabel(“Residuals”)
   plt.show()
   

这产生了以下图表：

图 7.5 – 预测值与残差的散点图

1. 我们将通过在模型上执行 k 折交叉验证来结束。这些分数低于线性回归模型的分数，但接近：

   X_train, X_test, y_train, y_test =  \
     train_test_split(features,\
     target, test_size=0.1, random_state=22)
   kf = KFold(n_splits=4, shuffle=True, random_state=0)
   scores = cross_validate(ttr, X=X_train, y=y_train,
     cv=kf, scoring=(‘r2’, ‘neg_mean_absolute_error’), n_jobs=1)
   print(“Mean Absolute Error: %.2f, R-squared: %.2f” % 
     (scores[‘test_neg_mean_absolute_error’].mean(),
     scores[‘test_r2’].mean()))
   Mean Absolute Error: -0.27, R-squared: 0.57
   

这给我们提供了一个模型，它并不比我们的原始模型更好，但它至少更有效地处理了特征选择过程。也有可能如果我们尝试不同的 alpha 超参数值，我们可能会得到更好的结果。为什么不试试 0.05 或 1.0 呢？我们将在接下来的两个步骤中尝试回答这个问题。

使用网格搜索调整超参数

确定超参数的最佳值，例如前一个例子中的 alpha 值，被称为超参数调整。scikit-learn 中用于超参数调整的一个工具是 GridSearchCV。CV 后缀表示交叉验证。

使用 GridSearchCV 非常简单。如果我们已经有了管道，就像我们在这个例子中做的那样，我们将它传递给一个 GridSearchCV 对象，以及一个参数字典。GridSearchCV 将尝试所有参数组合，并返回最佳组合。让我们在我们的 lasso 回归模型上试一试：

1. 首先，我们将实例化一个 lasso 对象，并创建一个包含要调整的超参数的字典。这个字典 lasso_params 表示我们想要尝试 0.05 和 0.9 之间的所有 alpha 值，以 0.5 的间隔。我们无法为字典键选择任何想要的名称。regressor__lasso__alpha 是基于管道中步骤的名称。注意，我们正在使用双下划线。单下划线将返回错误：

   lasso = Lasso()
   lasso_params = {‘regressor__lasso__alpha’: np.arange(0.05, 1, 0.05)}
   

2. 现在，我们可以运行网格搜索。我们将传递管道，在这个案例中是一个 TransformedTargetRegressor，以及字典到 GridSearchCV。best_params_ 属性表明最佳 alpha 值为 0.05。当我们使用该值时，我们得到一个 r-squared 值为 0.60：

   gs = GridSearchCV(ttr,param_grid=lasso_params, cv=5)
   gs.fit(X_train, y_train)
   gs.best_params_
   {‘regressor__lasso__alpha’: 0.05}
   gs.best_score_
   0.6028804486340877
   

Lasso 回归模型在平均绝对误差和 r-squared 方面接近线性模型，但并不完全一样。Lasso 回归的一个优点是，在训练我们的模型之前，我们不需要进行单独的特征选择步骤。（回想一下，对于包装特征选择方法，模型需要在特征选择期间以及之后进行训练，正如我们在 第五章 中讨论的，特征选择。）

使用非线性回归

线性回归假设特征与目标之间的关系在特征的范围内是恒定的。你可能还记得，我们在本章开头讨论的简单线性回归方程为每个特征有一个斜率估计：

![图片 B17978_07_010.jpg]

在这里，y 是目标，每个 x 是一个特征，每个 β 是一个系数（或截距）。如果特征与目标之间的真实关系是非线性的，我们的模型可能表现不佳。

幸运的是，当我们无法假设特征与目标之间存在线性关系时，我们仍然可以很好地利用 OLS。我们可以使用与上一节相同的线性回归算法，但使用特征的多项式变换。这被称为多项式回归。

我们给特征添加一个幂来运行多项式回归。这给我们以下方程：

![图片 B17978_07_011.jpg]

下面的图比较了线性回归与多项式回归的预测值。多项式曲线似乎比线性回归线更好地拟合了虚构的数据点：

图 7.6 – 多项式方程曲线和线性方程线的示意图

在本节中，我们将对全球气象站平均年温度的线性模型和非线性模型进行实验。我们将使用纬度和海拔作为特征。首先，我们将使用多元线性回归来预测温度，然后尝试使用多项式变换的模型。按照以下步骤进行：

1. 我们将首先导入必要的库。如果你一直在本章中工作，这些库将很熟悉：

   import pandas as pd
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler, PolynomialFeatures
   from sklearn.linear_model import LinearRegression
   from sklearn.pipeline import make_pipeline
   from sklearn.model_selection import cross_validate
   from sklearn.model_selection import KFold
   from sklearn.impute import KNNImputer
   import matplotlib.pyplot as plt
   

2. 我们还需要导入包含我们用于识别异常值类别的模块：

   import os
   import sys
   sys.path.append(os.getcwd() + “/helperfunctions”)
   from preprocfunc import OutlierTrans
   

3. 我们加载陆地温度数据，识别我们想要的特征，并生成一些描述性统计。对于海拔高度有一些缺失值，对于平均年温度有一些极端的负值。目标和特征值的范围差异很大，所以我们可能需要缩放我们的数据：

   landtemps = pd.read_csv(“data/landtempsb2019avgs.csv”)
   landtemps.set_index(‘locationid’, inplace=True)
   feature_cols = [‘latabs’,’elevation’]
   landtemps[[‘avgtemp’] + feature_cols].\
     agg([‘count’,’min’,’median’,’max’]).T
                count      min      median    max
   avgtemp      12,095    -60.82    10.45     33.93
   latabs       12,095     0.02     40.67     90.00
   elevation    12,088    -350.00   271.30    4,701.00
   

4. 接下来，我们创建训练和测试数据框：

   X_train, X_test, y_train, y_test =  \
     train_test_split(landtemps[feature_cols],\
     landtemps[[‘avgtemp’]], test_size=0.1, random_state=0)
   

5. 现在，我们构建一个管道来处理我们的预处理 – 将异常值设置为缺失，对所有缺失值进行 KNN 插补，并对特征进行缩放 – 然后运行线性模型。我们进行 10 折交叉验证，得到平均 r-squared 为 0.79，平均绝对误差为-2.8：

   lr = LinearRegression()
   knnimp = KNNImputer(n_neighbors=45)
   pipe1 = make_pipeline(OutlierTrans(3),knnimp,
     StandardScaler(), lr)
   ttr=TransformedTargetRegressor(regressor=pipe1,
     transformer=StandardScaler())
   kf = KFold(n_splits=10, shuffle=True, random_state=0)
   scores = cross_validate(ttr, X=X_train, y=y_train,
     cv=kf, scoring=(‘r2’, ‘neg_mean_absolute_error’), n_jobs=1)
   scores[‘test_r2’].mean(), scores[‘test_neg_mean_absolute_error’].mean()
   (0.7933471824738406, -2.8047627785750913)
   

注意，我们在识别异常值方面非常保守。我们传递了一个阈值为 3，这意味着一个值需要比四分位数范围高或低三倍以上。显然，我们通常会更多地考虑异常值的识别。在这里，我们只演示了在决定这样做是有意义之后，如何在管道中处理异常值。

1. 让我们看看预测值和残差。总体上几乎没有偏差（残差的平均值是 0），但有一些负偏斜：

   ttr.fit(X_train, y_train)
   pred = ttr.predict(X_test)
   preddf = pd.DataFrame(pred, columns=[‘prediction’],
     index=X_test.index).join(X_test).join(y_test)
   preddf.resid.agg([‘mean’,’median’,’skew’,’kurtosis’])
   mean                 0.00
   median               0.50
   skew                -1.13
   kurtosis             3.48
   Name: resid, dtype: float64
   

2. 如果我们创建残差的直方图，很容易看到这种偏斜。有一些极端的负残差 – 即我们过度预测平均温度的次数：

   plt.hist(preddf.resid, color=”blue”)
   plt.axvline(preddf.resid.mean(), color=’red’, 
     linestyle=’dashed’, linewidth=1)
   plt.title(“Histogram of Residuals for Linear Model of Temperature”)
   plt.xlabel(“Residuals”)
   plt.ylabel(“Frequency”)
   plt.show()
   

这会产生以下图表：

图 7.7 – 温度模型残差

1. 也可以通过绘制预测值与残差的关系图来有所帮助：

   plt.scatter(preddf.prediction, preddf.resid, color=”blue”)
   plt.axhline(0, color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Scatterplot of Predictions and Residuals”)
   plt.xlabel(“Predicted Temperature”)
   plt.ylabel(“Residuals”)
   plt.xlim(-20,40)
   plt.ylim(-27,10)
   plt.show()
   

这会产生以下图表：

图 7.8 – 预测值与残差的散点图

我们的模型在预测大约 28 摄氏度以上的所有预测值时都会过度预测。它也可能会在预测 18 到 28 之间的值时低估。

让我们看看多项式回归是否能得到更好的结果：

1. 首先，我们将创建一个PolynomialFeatures对象，其degree为4，并进行拟合。我们可以通过传递原始特征名称给get_feature_names方法来获取变换后返回的列名。每个特征的二次、三次和四次幂值被创建，以及变量之间的交互效应（例如latabs * elevation）。在这里我们不需要运行拟合，因为那将在管道中发生。我们只是在这里做这个来了解它是如何工作的：

   polytrans = PolynomialFeatures(degree=4, include_bias=False)
   polytrans.fit(X_train.dropna())
   featurenames = polytrans.get_feature_names(feature_cols)
   featurenames
   [‘latabs’,
    ‘elevation’,
    ‘latabs²’,
    ‘latabs elevation’,
    ‘elevation²’,
    ‘latabs³’,
    ‘latabs² elevation’,
    ‘latabs elevation²’,
    ‘elevation³’,
    ‘latabs⁴’,
    ‘latabs³ elevation’,
    ‘latabs² elevation²’,
    ‘latabs elevation³’,
    ‘elevation⁴’]
   

2. 接下来，我们将为多项式回归创建一个管道。这个管道基本上与线性回归相同，只是在 KNN 插补之后插入多项式变换步骤：

   pipe2 = make_pipeline(OutlierTrans(3), knnimp,
     polytrans, StandardScaler(), lr)
   ttr2 = TransformedTargetRegressor(regressor=pipe2,\
     transformer=StandardScaler())
   

3. 现在，让我们基于多项式模型创建预测值和残差。与线性模型相比，残差中的偏斜度要小一些：

   ttr2.fit(X_train, y_train)
   pred = ttr2.predict(X_test)
   preddf = pd.DataFrame(pred, columns=[‘prediction’],
     index=X_test.index).join(X_test).join(y_test)
   preddf[‘resid’] = preddf.avgtemp-preddf.prediction
   preddf.resid.agg([‘mean’,’median’,’skew’,’kurtosis’])
   mean                 0.01
   median               0.20
   skew                -0.98
   kurtosis             3.34
   Name: resid, dtype: float64
   

4. 我们应该看一下残差的直方图：

   plt.hist(preddf.resid, color=”blue”)
   plt.axvline(preddf.resid.mean(), color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Histogram of Residuals for Temperature Model”)
   plt.xlabel(“Residuals”)
   plt.ylabel(“Frequency”)
   plt.show()
   

这产生了以下图表：

图 7.9 – 温度模型残差

1. 让我们也做一个预测值与残差的散点图。这些看起来比线性模型的残差要好一些，尤其是在预测的上限范围内：

   plt.scatter(preddf.prediction, preddf.resid, color=”blue”)
   plt.axhline(0, color=’red’, linestyle=’dashed’, linewidth=1)
   plt.title(“Scatterplot of Predictions and Residuals”)
   plt.xlabel(“Predicted Temperature”)
   plt.ylabel(“Residuals”)
   plt.xlim(-20,40)
   plt.ylim(-27,10)
   plt.show()
   

这产生了以下图表：

图 7.10 – 预测值与残差的散点图

1. 让我们再次进行 k 折交叉验证，并取各折的平均 r-squared 值。与线性模型相比，r-squared 和平均绝对误差都有所提高：

   scores = cross_validate(ttr2, X=X_train, y=y_train,
     cv=kf, scoring=(‘r2’, ‘neg_mean_absolute_error’),
     n_jobs=1)
   scores[‘test_r2’].mean(), scores[‘test_neg_mean_absolute_error’].mean()
   (0.8323274036342788, -2.4035803290965507)
   

多项式变换提高了我们的整体结果，尤其是在预测变量的某些范围内。高温下的残差明显较低。当我们的残差表明特征和目标之间可能存在非线性关系时，尝试多项式变换通常是一个好主意。

梯度下降回归

梯度下降可以是有序最小二乘法优化线性模型损失函数的一个很好的替代方案。这在处理非常大的数据集时尤其正确。在本节中，我们将使用梯度下降和陆地温度数据集，主要为了演示如何使用它，并给我们另一个机会探索穷举网格搜索。让我们开始吧：

1. 首先，我们将加载我们迄今为止一直在使用的相同库，以及来自 scikit-learn 的随机梯度下降回归器：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.linear_model import SGDRegressor
   from sklearn.compose import TransformedTargetRegressor
   from sklearn.pipeline import make_pipeline
   from sklearn.impute import KNNImputer
   from sklearn.model_selection import GridSearchCV
   import os
   import sys
   sys.path.append(os.getcwd() + “/helperfunctions”)
   from preprocfunc import OutlierTrans
   

2. 然后，我们将再次加载陆地温度数据集并创建训练和测试数据框：

   landtemps = pd.read_csv(“data/landtempsb2019avgs.csv”)
   landtemps.set_index(‘locationid’, inplace=True)
   feature_cols = [‘latabs’,’elevation’]
   X_train, X_test, y_train, y_test =  \
     train_test_split(landtemps[feature_cols],\
     landtemps[[‘avgtemp’]], test_size=0.1, random_state=0)
   

3. 接下来，我们将设置一个管道来处理异常值，在运行梯度下降之前填充缺失值并缩放数据：

   knnimp = KNNImputer(n_neighbors=45)
   sgdr = SGDRegressor()
   pipe1 = make_pipeline(OutlierTrans(3),knnimp,StandardScaler(), sgdr)
   ttr=TransformedTargetRegressor(regressor=pipe1,transformer=StandardScaler())
   

4. 现在，我们需要创建一个字典来指示我们想要调整的超参数和要尝试的值。我们想要尝试 alpha、损失函数、epsilon 和惩罚的值。我们将为字典中的每个键添加前缀regressor__sgdregressor__，因为随机梯度下降回归器可以在管道中找到。

alpha参数决定了惩罚的大小。默认值是0.0001。我们可以选择 L1、L2 或弹性网络正则化。我们将选择huber和epsilon_insensitive作为要包含在搜索中的损失函数。默认损失函数是squared_error，但这只会给我们普通的有序最小二乘法。

huber损失函数对异常值不如 OLS 敏感。它的敏感度取决于我们指定的 epsilon 的值。使用epsilon_insensitive损失函数，给定范围（epsilon）内的错误不会受到惩罚（我们将在下一章构建具有 epsilon-insensitive 管的模型，我们将检查支持向量回归）：

sgdr_params = {
 ‘regressor__sgdregressor__alpha’: 10.0 ** -np.arange(1, 7),
 ‘regressor__sgdregressor__loss’: [‘huber’,’epsilon_insensitive’],
 ‘regressor__sgdregressor__penalty’: [‘l2’, ‘l1’, ‘elasticnet’],
 ‘regressor__sgdregressor__epsilon’: np.arange(0.1, 1.6, 0.1)
}

1. 现在，我们已经准备好运行一个全面的网格搜索。最佳参数是 alpha 为0.001，epsilon 为1.3，损失函数为huber，正则化为弹性网络：

   gs = GridSearchCV(ttr,param_grid=sgdr_params, cv=5, scoring=”r2”)
   gs.fit(X_train, y_train)
   gs.best_params_
   {‘regressor__sgdregressor__alpha’: 0.001,
    ‘regressor__sgdregressor__epsilon’: 1.3000000000000003,
    ‘regressor__sgdregressor__loss’: ‘huber’,
    ‘regressor__sgdregressor__penalty’: ‘elasticnet’}
   gs.best_score_
   0.7941051735846133
   

2. 我通常发现查看一些其他网格搜索迭代的超参数很有帮助。具有弹性网络或 L2 正则化的 Huber 损失模型表现最佳：

   Results = \
     pd.DataFrame(gs.cv_results_[‘mean_test_score’], \
       columns=[‘meanscore’]).\
     join(pd.DataFrame(gs.cv_results_[‘params’])).\
     sort_values([‘meanscore’], ascending=False)
   results.head(3).T
                                            254      252      534
   meanscore                           0.794105 0.794011 0.794009
   regressor__sgdregressor__alpha      0.001000 0.001000 0.000001
   regressor__sgdregressor__epsilon    1.300000 1.300000 1.500000
   regressor__sgdregressor__loss          huber    huber    huber
   regressor__sgdregressor__penalty  elasticnet       l2       l2
   

随机梯度下降是一种通用的优化方法，可以应用于许多机器学习问题。它通常非常高效，正如我们在这里所看到的。我们能够相当快地运行一个关于惩罚、惩罚大小、epsilon 和损失函数的全面网格搜索。

摘要

本章使我们能够探索一个非常著名的机器学习算法：线性回归。我们考察了使特征空间成为线性模型良好候选者的特性。我们还探讨了在必要时如何通过正则化和变换来改进线性模型。然后，我们研究了随机梯度下降作为 OLS 优化的替代方案。我们还学习了如何将我们自己的类添加到管道中以及如何进行超参数调整。

在下一章中，我们将探索支持向量回归。

第八章：第八章：支持向量回归

支持向量回归（SVR）在线性回归模型的假设不成立时可以是一个极佳的选择，例如当我们的特征与目标之间的关系过于复杂，无法用权重的线性组合来描述时。甚至更好，SVR 允许我们无需扩展特征空间来建模这种复杂性。

支持向量机识别出最大化两个类别之间边界的超平面。支持向量是距离边界最近的数据点，它们支持这个边界，如果可以这样表达的话。这证明它在回归建模中与在分类中一样有用。SVR 找到包含最多数据点的超平面。我们将在本章的第一节中讨论它是如何工作的。

与普通最小二乘回归不同，SVR 不是最小化平方残差的和，而是在一个可接受的误差范围内最小化系数。像岭回归和 Lasso 回归一样，这可以减少模型方差和过拟合的风险。SVR 在处理中小型数据集时效果最佳。

该算法也非常灵活，允许我们指定可接受的误差范围，使用核来建模非线性关系，并调整超参数以获得最佳的偏差-方差权衡。我们将在本章中展示这一点。

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

• SVR 的关键概念

• 基于线性模型的支持向量回归

• 使用核进行非线性 SVR

技术要求

在本章中，我们将使用 scikit-learn 和matplotlib库。您可以使用pip安装这些包。

SVR 的关键概念

我们将从这个部分开始讨论支持向量机在分类中的应用。在这里我们不会深入细节，将支持向量分类的详细讨论留给第十三章，支持向量机分类。但首先从支持向量机在分类中的应用开始，将很好地过渡到 SVR 的解释。

正如我在本章开头讨论的那样，支持向量机找到最大化类别之间边界的超平面。当只有两个特征存在时，这个超平面仅仅是一条线。考虑以下示例图：

图 8.1 – 基于两个特征的支持向量机分类

图中的两个类别，用红色圆圈和蓝色正方形表示，可以使用两个特征，x1 和 x2，进行线性可分。粗体线是决策边界。它是每个类别离边界数据点最远的线，或者说是最大边界。这些点被称为支持向量。

由于前一个图中的数据是线性可分的，我们可以无问题地使用所谓的硬间隔分类；也就是说，我们可以对每个类别的所有观测值位于决策边界正确一侧的要求非常严格。但如果我们数据点的样子像下面这个图所示的呢？

图 8.2 – 具有软间隔的支持向量机分类

这些数据点不是线性可分的。在这种情况下，我们可以选择软间隔分类并忽略异常值红色圆圈。

我们将在第十三章“支持向量机分类”中更详细地讨论支持向量分类，但这也说明了支持向量机的一些关键概念。这些概念可以很好地应用于涉及连续目标值的模型。这被称为支持向量回归或SVR。

在构建 SVR 模型时，我们决定可接受的预测误差量，ɛ。在一个特征模型中，预测值 在 ɛ 范围内的误差不会被惩罚。这有时被称为 epsilon-insensitive tube。SVR 最小化系数，使得所有数据点都落在该范围内。这在下图中得到说明：

图 8.3 – 具有可接受误差范围的 SVR

更精确地说，SVR 在满足误差 ε 不超过给定量的约束条件下，最小化系数的平方。

它在满足 的约束条件下最小化 ，其中 是权重（或系数）向量， 是实际目标值与预测值的差， 是可接受的误差量。

当然，期望所有数据点都落在期望范围内是不合理的。但我们仍然可以寻求最小化这种偏差。让我们用 ξ 表示偏离边缘的距离。这给我们一个新的目标函数。

我们在满足 的约束条件下最小化 ，其中 C 是一个超参数，表示模型对边缘外误差的容忍度。C 的值为 0 表示它根本不容忍那些大误差。这等价于原始目标函数：

图 8.4 – 具有超出可接受范围的点的 SVR

在这里，我们可以看到 SVR 的几个优点。有时，我们的误差不超过一定量比选择一个具有最低绝对误差的模型更重要。如果我们经常略微偏离但很少大幅偏离，可能比我们经常准确但偶尔严重偏离更重要。由于这种方法也最小化了我们的权重，它具有与正则化相同的优点，我们减少了过拟合的可能性。

非线性 SVR 和核技巧

我们尚未完全解决 SVR 中线性可分性的问题。为了简单起见，我们将回到一个涉及两个特征的分类问题。让我们看看两个特征与分类目标的关系图。目标有两个可能的值，由点和正方形表示。x1 和 x2 是数值，具有负值：

图 8.5 – 使用两个特征无法线性分离的类别标签

在这种情况下，我们该如何识别类别之间的边界呢？通常情况下，在更高的维度上可以识别出边界。在这个例子中，我们可以使用多项式变换，如下面的图表所示：

图 8.6 – 使用多项式变换建立边界

现在有一个第三维度，它是 x1 和 x2 平方的和。所有的点都高于平方。这与我们在上一章中使用多项式变换来指定非线性回归模型的方式相似。

这种方法的缺点之一是我们可能会迅速拥有太多特征，导致模型无法良好地表现。这时，核技巧就非常实用了。SVR 可以通过使用核函数隐式地扩展特征空间，而不实际创建更多特征。这是通过创建一个值向量来完成的，这些值可以用来拟合非线性边界。

虽然这允许我们拟合如前述图表中所示的一个假设的多项式变换，但 SVR 中最常用的核函数是径向基函数（RBF）。RBF 之所以受欢迎，是因为它比其他常见的核函数更快，并且因为它的伽马参数使其非常灵活。我们将在本章的最后部分探讨如何使用它。

但现在，让我们从一个相对简单的线性模型开始，看看 SVR 的实际应用。

线性模型的 SVR

我们通常有足够的领域知识，可以采取比仅仅最小化训练数据中的预测误差更细致的方法。利用这些知识，我们可能允许模型接受更多的偏差，当少量的偏差在实质上并不重要时，以减少方差。在使用 SVR 时，我们可以调整超参数，如 epsilon（可接受的误差范围）和C（调整该范围之外错误的容忍度），以改善模型的表现。

如果线性模型可以在你的数据上表现良好，线性 SVR 可能是一个不错的选择。我们可以使用 scikit-learn 的LinearSVR类构建一个线性 SVR 模型。让我们尝试使用我们在上一章中使用的汽油税数据创建一个线性 SVR 模型：

1. 我们需要使用与上一章相同的许多库来创建训练和测试 DataFrame，并预处理数据。我们还需要从 scikit-learn 和 scipy 中分别导入LinearSVR和uniform模块：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.svm import LinearSVR
   from scipy.stats import uniform
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.compose import ColumnTransformer
   from sklearn.compose import TransformedTargetRegressor
   from sklearn.impute import KNNImputer
   from sklearn.model_selection import cross_validate, \
     KFold, GridSearchCV, RandomizedSearchCV
   import sklearn.metrics as skmet
   import matplotlib.pyplot as plt
   

2. 我们还需要导入OutlierTrans类，这是我们首先在第七章中讨论的，线性回归模型，用于处理异常值：

   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

3. 接下来，我们加载汽油税数据并创建训练和测试 DataFrame。我们创建了数值和二进制特征的列表，以及一个单独的motorization_rate列表。正如我们在上一章查看数据时所见，我们需要对motorization_rate进行一些额外的预处理。

这个数据集包含了 2014 年每个国家的汽油税数据，以及燃料收入依赖性和衡量民主制度强度的指标：polity、democracy_polity和autocracy_polity。democracy_polity是一个二进制polity变量，对于polity得分高的国家取值为 1。autocracy_polity对于polity得分低的国家取值为 1。polity特征是衡量一个国家民主程度的一个指标：

fftaxrate14 = pd.read_csv("data/fossilfueltaxrate14.csv")
fftaxrate14.set_index('countrycode', inplace=True)
num_cols = ['fuel_income_dependence',
  'national_income_per_cap', 'VAT_Rate',  
  'gov_debt_per_gdp', 'polity','goveffect',
  'democracy_index']
dummy_cols = 'democracy_polity','autocracy_polity',
  'democracy', 'nat_oil_comp','nat_oil_comp_state']
spec_cols = ['motorization_rate']
target = fftaxrate14[['gas_tax_imp']]
features = fftaxrate14[num_cols + dummy_cols + spec_cols]
X_train, X_test, y_train, y_test =  \
  train_test_split(features,\
    target, test_size=0.2, random_state=0)

1. 让我们看看训练数据的摘要统计。我们需要标准化数据，因为数据范围差异很大，SVR 在标准化数据上表现更好。注意，motorization_rate有很多缺失值。我们可能想在这个特征上做得比简单的插补更好。对于虚拟列，我们有相当多的非缺失计数：

   X_train.shape
   (123, 13)
   X_train[num_cols + spec_cols].\
     agg(['count','min','median','max']).T
                         count min    median   max
   fuel_income_dependence  121 0.00   0.10     34.23
   national_income_per_cap 121 260.00 6,110.00 104,540.00
   VAT_Rate                121 0.00   16.00    27.00
   gov_debt_per_gdp        112 1.56   38.45    194.76
   polity                  121 -10.00 6.00     10.00
   goveffect               123 -2.04  -0.10    2.18
   democracy_index         121 0.03   0.54     0.93
   motorization_rate       100 0.00   0.20     0.81
   X_train[dummy_cols].apply(pd.value_counts, normalize=True).T
                                       0.00         1.00
   democracy_polity                    0.42         0.58
   autocracy_polity                    0.88         0.12
   democracy                           0.41         0.59
   nat_oil_comp                        0.54         0.46
   nat_oil_comp_state                  0.76         0.24
   X_train[dummy_cols].count()
   democracy_polity           121
   autocracy_polity           121
   democracy                  123
   nat_oil_comp               121
   nat_oil_comp_state         121
   

2. 我们需要构建一个列转换器来处理不同的数据类型。我们可以使用SimpleImputer来处理分类特征和数值特征，除了motorization_rate。我们将在稍后使用 KNN 插补来处理motorization_rate特征：

   standtrans = make_pipeline(OutlierTrans(2), 
    SimpleImputer(strategy="median"), StandardScaler())
   cattrans = make_pipeline(SimpleImputer(strategy="most_frequent"))
   spectrans = make_pipeline(OutlierTrans(2), StandardScaler())
   coltrans = ColumnTransformer(
     transformers=[
       ("stand", standtrans, num_cols),
       ("cat", cattrans, dummy_cols),
       ("spec", spectrans, spec_cols)
     ]
   )
   

3. 现在，我们已经准备好拟合我们的线性 SVR 模型。我们将为epsilon选择0.2的值。这意味着我们对实际值 0.2 标准差范围内的任何误差都感到满意（我们使用TransformedTargetRegressor来标准化目标）。我们将把C——决定模型对 epsilon 之外值容忍度的超参数——保留在其默认值 1.0。

在我们拟合模型之前，我们仍然需要处理motorization_rate的缺失值。我们将在列转换之后将 KNN 填充器添加到管道中。由于motorization_rate将是列转换后唯一的具有缺失值的特征，KNN 填充器只会改变该特征的价值。

我们需要使用目标转换器，因为列转换器只会改变特征，而不会改变目标。我们将把刚刚创建的管道传递给目标转换器的regressor参数以进行特征转换，并指出我们只想对目标进行标准缩放。

注意，线性 SVR 的默认损失函数是 L1，但我们可以选择 L2：

svr = LinearSVR(epsilon=0.2, max_iter=10000, 
  random_state=0)
pipe1 = make_pipeline(coltrans, 
  KNNImputer(n_neighbors=5), svr)
ttr=TransformedTargetRegressor(regressor=pipe1,
  transformer=StandardScaler())
ttr.fit(X_train, y_train)

1. 我们可以使用ttr.regressor_来访问管道中的所有元素，包括linearsvr对象。这就是我们如何获得coef_属性。与 0 显著不同的系数是VAT_Rate以及专制和国家级石油公司虚拟变量。我们的模型估计，在其他条件相同的情况下，增值税率和汽油税之间存在正相关关系。它估计拥有专制或国家级石油公司与汽油税之间存在负相关关系：

   coefs = ttr.regressor_['linearsvr'].coef_
   np.column_stack((coefs.ravel(), num_cols + dummy_cols + spec_cols))
   array([['-0.03040694175014407', 'fuel_income_dependence'],
          ['0.10549935644031803', 'national_income_per_cap'],
          ['0.49519936241642026', 'VAT_Rate'],
          ['0.0857845735264331', 'gov_debt_per_gdp'],
          ['0.018198547504343885', 'polity'],
          ['0.12656984468734492', 'goveffect'],
          ['-0.09889163752261303', 'democracy_index'],
          ['-0.036584519840546594', 'democracy_polity'],
          ['-0.5446613604546718', 'autocracy_polity'],
          ['0.033234557366924815', 'democracy'],
          ['-0.2048732386478349', 'nat_oil_comp'],
          ['-0.6142887840649164', 'nat_oil_comp_state'],
          ['0.14488410358761755', 'motorization_rate']], dtype='<U32')
   

注意，我们在这里没有进行任何特征选择。相反，我们依赖于 L1 正则化将特征系数推向接近 0。如果我们有更多的特征，或者我们更关心计算时间，那么仔细思考我们的特征选择策略就很重要了。

1. 让我们在该模型上进行一些交叉验证。平均绝对误差和 r-squared 值并不理想，这当然受到了小样本量的影响：

   kf = KFold(n_splits=3, shuffle=True, random_state=0)
   ttr.fit(X_train, y_train)
   scores = cross_validate(ttr, X=X_train, y=y_train,
     cv=kf, scoring=('r2', 'neg_mean_absolute_error'),
       n_jobs=1)
   print("Mean Absolute Error: %.2f, R-squared: %.2f" %
     (scores['test_neg_mean_absolute_error'].mean(),
     scores['test_r2'].mean()))
   Mean Absolute Error: -0.26, R-squared: 0.57
   

我们还没有进行任何超参数调整。我们不知道我们的epsilon和C的值是否是模型的最佳值。因此，我们需要进行网格搜索来尝试不同的超参数值。我们将从穷举网格搜索开始，这通常不切实际（除非你有一个性能相当高的机器，否则我建议不要运行接下来的几个步骤）。在穷举网格搜索之后，我们将进行随机网格搜索，这通常对系统资源的影响要小得多。

1. 我们将首先创建一个没有指定epsilon超参数的LinearSVR对象，并重新创建管道。然后，我们将创建一个字典svr_params，其中包含用于检查epsilon和C的值，分别称为regressor_linearsvr_epsilon和regressor_linearsvr_C。

记住我们从前一章的网格搜索中提到的，键的名称必须与我们的管道步骤相对应。在这个例子中，我们可以通过转换目标对象的regressor属性来访问我们的管道。管道中有一个linearsvr对象，具有epsilon和C的属性。

我们将svr_params字典传递给GridSearchCV对象，并指示我们希望评分基于 r-squared（如果我们想基于平均绝对误差进行评分，我们也可以这样做）。

然后，我们将运行网格搜索对象的fit方法。我应该重复之前提到的警告，你可能不希望运行穷举网格搜索，除非你使用的是高性能机器，或者你不在乎在去拿一杯咖啡的时候让它运行。请注意，在我的机器上每次运行大约需要 26 秒：

svr = LinearSVR(max_iter=100000, random_state=0)
pipe1 = make_pipeline(coltrans, 
  KNNImputer(n_neighbors=5), svr)
ttr=TransformedTargetRegressor(regressor=pipe1,
  transformer=StandardScaler())
svr_params = {
  'regressor__linearsvr__epsilon': np.arange(0.1, 1.6, 0.1),
  'regressor__linearsvr__C': np.arange(0.1, 1.6, 0.1)
}
gs = GridSearchCV(ttr,param_grid=svr_params, cv=3, 
  scoring='r2')
%timeit gs.fit(X_train, y_train)
26.2 s ± 50.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

1. 现在，我们可以使用网格搜索的best_params_属性来获取与最高分数相关的超参数。我们可以通过best_scores_属性查看这些参数的分数。这告诉我们，我们以 0.1 的C值和 0.2 的epsilon值获得了最高的 r-squared 值，即 0.6：

   gs.best_params_
   {'regressor__linearsvr__C': 0.1, 'regressor__linearsvr__epsilon': 0.2}
   gs.best_score_
   0.599751107082899
   

了解为我们的超参数选择哪些值是很好的。然而，穷举网格搜索在计算上相当昂贵。让我们尝试随机搜索。

1. 我们将指示epsilon和C的随机值应来自介于 0 和 1.5 之间的均匀分布。然后，我们将该字典传递给RandomizedSearchCV对象。这比穷举网格搜索快得多——每次迭代略超过 1 秒。这给我们提供了比穷举网格搜索更高的epsilon和C值——即，分别为 0.23 和 0.7。然而，r-squared 值略低：

   svr_params = {
    'regressor__linearsvr__epsilon': uniform(loc=0, scale=1.5),
    'regressor__linearsvr__C': uniform(loc=0, scale=1.5)
   }
   rs = RandomizedSearchCV(ttr, svr_params, cv=3, scoring='r2')
   %timeit rs.fit(X_train, y_train)
   1.21 s ± 24.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
   rs.best_params_
   {'regressor__linearsvr__C': 0.23062453444814285,
    'regressor__linearsvr__epsilon': 0.6976844872643301}
   rs.best_score_
   0.5785452537781279
   

2. 让我们查看基于随机网格搜索最佳模型的预测。随机网格搜索对象的predict方法可以为我们生成这些预测：

   pred = rs.predict(X_test)
   preddf = pd.DataFrame(pred, columns=['prediction'],
     index=X_test.index).join(X_test).join(y_test)
   preddf['resid'] = preddf.gas_tax_imp-preddf.prediction
   

3. 现在，让我们看看我们残差的分布：

   plt.hist(preddf.resid, color="blue", bins=np.arange(-0.5,1.0,0.25))
   plt.axvline(preddf.resid.mean(), color='red', linestyle='dashed', linewidth=1)
   plt.title("Histogram of Residuals for Gas Tax Model")
   plt.xlabel("Residuals")
   plt.ylabel("Frequency")
   plt.xlim()
   plt.show()
   

这产生了以下图表：

图 8.7 – 加油税线性 SVR 模型的残差分布

在这里，存在一点偏差（整体上有些过度预测）和一些正偏斜。

1. 让我们再查看一下预测值与残差之间的散点图：

   plt.scatter(preddf.prediction, preddf.resid, color="blue")
   plt.axhline(0, color='red', linestyle='dashed', linewidth=1)
   plt.title("Scatterplot of Predictions and Residuals")
   plt.xlabel("Predicted Gas Tax")
   plt.ylabel("Residuals")
   plt.show()
   

这产生了以下图表：

图 8.8 – 加油税线性 SVR 模型的预测值与残差的散点图

这些残差是有问题的。我们在预测值的低和高范围内总是过度预测（预测值高于实际值）。这不是我们想要的，也许在提醒我们存在未考虑的非线性关系。

当我们的数据是线性可分时，线性支持向量回归（SVR）可以是一个高效的选择。它可以用在许多我们本会使用线性回归或带有正则化的线性回归的相同情况下。它的相对效率意味着我们对于使用包含超过 10,000 个观测值的数据集并不像使用非线性 SVR 那样担忧。然而，当线性可分性不可能时，我们应该探索非线性模型。

使用核函数进行非线性 SVR

回想一下本章开头我们讨论的内容，我们可以使用核函数来拟合一个非线性ε敏感管。在本节中，我们将使用我们在上一章中使用过的土地温度数据运行非线性 SVR。但首先，我们将使用相同的数据构建一个线性 SVR 以进行比较。

我们将把气象站的平均温度建模为纬度和海拔的函数。按照以下步骤进行：

1. 我们将首先加载熟悉的库。唯一的新类是来自 scikit-learn 的SVR：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.svm import LinearSVR, SVR
   from scipy.stats import uniform
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.compose import TransformedTargetRegressor
   from sklearn.impute import KNNImputer
   from sklearn.model_selection import RandomizedSearchCV
   import sklearn.metrics as skmet
   import matplotlib.pyplot as plt
   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

2. 接下来，我们将加载土地温度数据并创建训练和测试数据框。我们还将查看一些描述性统计。海拔高度有多个缺失值，两个特征的范围差异很大。还有一些异常低的平均温度：

   landtemps = pd.read_csv("data/landtempsb2019avgs.csv")
   landtemps.set_index('locationid', inplace=True)
   feature_cols = ['latabs','elevation']
   landtemps[['avgtemp'] + feature_cols].\
     agg(['count','min','median','max']).T
                  count       min     median     max
   avgtemp        12,095     -61      10         34
   latabs         12,095      0       41         90
   elevation      12,088     -350     271        4,701
   X_train, X_test, y_train, y_test =  \
     train_test_split(landtemps[feature_cols],\
     landtemps[['avgtemp']], test_size=0.1, random_state=0)
   

3. 让我们从平均温度的线性 SVR 模型开始。我们可以对处理异常值保持相当保守，只有当四分位数范围超过三倍时，才将它们设置为缺失值。（我们在第七章，线性回归模型）中创建了OutlierTrans类。）我们将使用 KNN 插补缺失的海拔值并缩放数据。记住，我们需要使用目标转换器来缩放目标变量。

正如我们在上一节中所做的那样，我们将使用一个字典，svr_params，来表示我们想要从均匀分布中采样超参数的值——即epsilon和C。我们将把这个字典传递给RandomizedSearchCV对象。

在运行fit之后，我们可以得到epsilon和C的最佳参数，以及最佳模型的平均绝对误差。平均绝对误差相当不错，大约为 2.8 度：

svr = LinearSVR(epsilon=1.0, max_iter=100000)
knnimp = KNNImputer(n_neighbors=45)
pipe1 = make_pipeline(OutlierTrans(3), knnimp, StandardScaler(), svr)
ttr=TransformedTargetRegressor(regressor=pipe1,
  transformer=StandardScaler())
svr_params = {
 'regressor__linearsvr__epsilon': uniform(loc=0, scale=1.5),
 'regressor__linearsvr__C': uniform(loc=0, scale=20)
}
rs = RandomizedSearchCV(ttr, svr_params, cv=10, scoring='neg_mean_absolute_error')
rs.fit(X_train, y_train)
rs.best_params_
{'regressor__linearsvr__C': 15.07662849482442,
 'regressor__linearsvr__epsilon': 0.06750238486004034}
rs.best_score_
-2.769283402595076

1. 让我们看看预测结果：

   pred = rs.predict(X_test)
   preddf = pd.DataFrame(pred, columns=['prediction'],
     index=X_test.index).join(X_test).join(y_test)
   preddf['resid'] = preddf.avgtemp-preddf.prediction
   plt.scatter(preddf.prediction, preddf.resid, color="blue")
   plt.axhline(0, color='red', linestyle='dashed', linewidth=1)
   plt.title("Scatterplot of Predictions and Residuals")
   plt.xlabel("Predicted Gas Tax")
   plt.ylabel("Residuals")
   plt.show()
   

这产生了以下图表：

图 8.9 – 土地温度线性 SVR 模型的预测值和残差散点图

在预测值的上限范围内有相当多的过度预测。我们通常在预测的汽油税值 15 到 25 度之间低估值。也许我们可以通过非线性模型来提高拟合度。

1. 运行非线性 SVR 不需要做太多改变。我们只需要创建一个SVR对象并选择一个核函数。通常选择rbf。（除非你使用的是良好的硬件，或者你不在乎暂时做其他事情然后回来获取结果，否则你可能不想在你的机器上拟合这个模型。）看看下面的代码：

   svr = SVR(kernel='rbf')
   pipe1 = make_pipeline(OutlierTrans(3), knnimp, StandardScaler(), svr)
   ttr=TransformedTargetRegressor(regressor=pipe1,
     transformer=StandardScaler())
   svr_params = {
    'regressor__svr__epsilon': uniform(loc=0, scale=5),
    'regressor__svr__C': uniform(loc=0, scale=20),
    'regressor__svr__gamma': uniform(loc=0, scale=100)
    }
   rs = RandomizedSearchCV(ttr, svr_params, cv=10, scoring='neg_mean_absolute_error')
   rs.fit(X_train, y_train)
   rs.best_params_
   {'regressor__svr__C': 5.3715128489311255,
    'regressor__svr__epsilon': 0.03997496426101643,
    'regressor__svr__gamma': 53.867632383007994}
   rs.best_score_
   -2.1319240416548775
   

在平均绝对误差方面有明显的改进。在这里，我们可以看到gamma和 C 超参数为我们做了很多工作。如果我们对平均偏差 2 度左右没有异议，这个模型就能达到这个目标。

在第十三章**，支持向量机分类中，我们详细讨论了 gamma 和 C 超参数。我们还探讨了除了 rbf 核以外的其他核函数。

1. 让我们再次查看残差，看看我们的误差分布是否有问题，就像我们的线性模型那样：

   pred = rs.predict(X_test)
   preddf = pd.DataFrame(pred, columns=['prediction'],
     index=X_test.index).join(X_test).join(y_test)
   preddf['resid'] = preddf.avgtemp-preddf.prediction
   plt.scatter(preddf.prediction, preddf.resid, color="blue")
   plt.axhline(0, color='red', linestyle='dashed', linewidth=1)
   plt.title("Scatterplot of Predictions and Residuals")
   plt.xlabel("Predicted Gas Tax")
   plt.ylabel("Residuals")
   plt.show()
   

这产生了以下图表：

图 8.10 – 非线性 SVR 模型预测值和残差的散点图

这些残差看起来比线性模型的残差要好得多。

这说明了使用核函数如何在不增加特征空间的情况下增加我们模型的复杂性。通过使用rbf核并调整 C 和gamma超参数，我们解决了线性模型中看到的欠拟合问题。这是非线性 SVR 的一个巨大优点。缺点，正如我们之前看到的，是对系统资源的需求很大。包含 12,000 个观测值的数据库是非线性 SVR 可以轻松处理的极限，尤其是在进行最佳超参数的网格搜索时。

摘要

本章中的示例展示了 SVR 的一些优点。该算法允许我们调整超参数来解决欠拟合或过拟合问题。这可以在不增加特征数量的情况下完成。与线性回归等方法相比，SVR 对异常值也不那么敏感。

当我们可以用线性 SVR 构建一个好的模型时，这是一个完全合理的选择。它比非线性模型训练得快得多。然而，我们通常可以通过非线性 SVR 来提高性能，就像我们在本章的最后部分看到的那样。

这一讨论引出了我们在下一章将要探讨的内容，我们将探讨两种流行的非参数回归算法：k 近邻和决策树回归。这两个算法对特征和目标分布几乎没有假设。与 SVR 类似，它们可以在不增加特征空间的情况下捕捉数据中的复杂关系。

第九章：第九章：K 近邻、决策树、随机森林和梯度提升回归

对于支持向量机而言，K 近邻和决策树模型最著名的用途是分类模型。然而，它们也可以用于回归，并且相对于经典的线性回归，它们具有一些优势。K 近邻和决策树可以很好地处理非线性，并且不需要对特征的高斯分布做出任何假设。此外，通过调整K 近邻（KNN）的k值或决策树的最大深度，我们可以避免对训练数据拟合得太精确。

这使我们回到了前两章的主题——如何在不过度拟合的情况下增加模型复杂度，包括考虑非线性。我们已经看到，允许一些偏差可以减少方差，并给我们提供更可靠的模型性能估计。我们将继续在本章中探索这种平衡。

具体来说，我们将涵盖以下主要内容：

• K 近邻回归的关键概念

• K 近邻回归

• 决策树和随机森林回归的关键概念

• 决策树和随机森林回归

• 使用梯度提升回归

技术要求

在本章中，我们将使用 scikit-learn 和matplotlib库。我们还将使用 XGBoost。你可以使用pip来安装这些包。

K 近邻回归的关键概念

KNN 算法的部分吸引力在于它相当直接且易于解释。对于每个需要预测目标的观测，KNN 会找到k个训练观测，其特征与该观测最相似。当目标是分类时，KNN 会选择k个训练观测中最频繁的目标值。（我们通常在分类问题中选择一个奇数作为k，以避免平局。）

当目标是数值时，KNN 会给出k个训练观测的目标平均值。通过训练观测，我指的是那些具有已知目标值的观测。KNN 实际上并不进行真正的训练，因为它被称为懒惰学习器。我将在本节稍后详细讨论这一点。

图 9.1展示了使用 K 近邻进行分类，其中k的值为 1 和 3：当k为 1 时，我们的新观测将被分配红色标签。当k为 3 时，它将被分配蓝色标签：

图 9.1 – K 近邻，k 值为 1 和 3

但我们所说的相似或最近的观测意味着什么？有几种方法可以衡量相似度，但一个常见的衡量标准是欧几里得距离。欧几里得距离是两点之间平方差的和。这可能会让你想起勾股定理。从点 a 到点 b 的欧几里得距离如下：

曼哈顿距离是欧几里得距离的一个合理替代方案。从点 a 到点 b 的曼哈顿距离如下：

曼哈顿距离有时也被称为出租车距离。这是因为它反映了两个点在网格路径上的距离。图 9.2 展示了曼哈顿距离，并将其与欧几里得距离进行了比较：

图 9.2 – 欧几里得和曼哈顿距离度量

当特征在类型或尺度上非常不同时，使用曼哈顿距离可以获得更好的结果。然而，我们可以将距离度量的选择视为一个经验问题；也就是说，我们可以尝试两者（或其他距离度量），并看看哪个给我们提供了性能最好的模型。我们将在下一节通过网格搜索来演示这一点。

如你所怀疑的，KNN 模型的敏感性取决于 k 的选择。k 的值较低时，模型会试图识别观测值之间的细微差别。当然，在 k 值非常低的情况下，过拟合的风险很大。但在 k 值较高时，我们的模型可能不够灵活。我们再次面临方差-偏差权衡。较低的 k 值导致偏差较小而方差较大，而较高的 k 值则相反。

对于 k 的选择没有明确的答案。一个很好的经验法则是从观测值的平方根开始。然而，正如我们对距离度量的处理一样，我们应该测试模型在不同 k 值下的性能。

如我之前提到的，K 近邻是一种懒惰学习算法。在训练时间不进行任何计算。学习主要发生在测试期间。这有一些缺点。当数据中有许多实例或维度时，KNN 可能不是一个好的选择，而且预测速度很重要。它也往往在稀疏数据（即具有许多 0 值的数据集）的情况下表现不佳。

K 近邻是非参数算法。不对底层数据的属性做出假设，例如线性或正态分布的特征。当线性模型不起作用时，它通常能给出相当不错的结果。我们将在下一节构建 KNN 回归模型。

K 近邻回归

如前所述，当普通最小二乘法的假设不成立，且观测值和维度数量较少时，K 近邻算法可以成为线性回归的一个很好的替代方案。它也非常容易指定，因此即使我们最终模型中不使用它，它也可以用于诊断目的。

在本节中，我们将使用 KNN 构建一个国家层面的女性与男性收入比模型。我们将基于劳动力参与率、教育成就、青少年出生频率和最高级别的女性政治参与率。这是一个很好的数据集进行实验，因为小样本量和特征空间意味着它不太可能消耗你的系统资源。特征数量较少也使得它更容易解释。缺点是可能难以找到显著的结果。话虽如此，让我们看看我们发现了什么。

注意

我们将在本章中一直使用收入差距数据集。该数据集由联合国开发计划署在www.kaggle.com/datasets/undp/human-development上提供供公众使用。每个国家都有一个记录，包含按性别划分的 2015 年的综合就业、收入和教育数据。

让我们开始构建我们的模型：

1. 首先，我们必须导入我们在前两章中使用的一些相同的sklearn库。我们还必须导入KNeighborsRegressor和来自第五章的老朋友——特征选择——即SelectFromModel。我们将使用SelectFromModel来将特征选择添加到我们将构建的管道中：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.model_selection import RandomizedSearchCV
   from sklearn.neighbors import KNeighborsRegressor
   from sklearn.linear_model import LinearRegression
   from sklearn.feature_selection import SelectFromModel
   import seaborn as sns
   import matplotlib.pyplot as plt
   

2. 我们还需要在第七章的线性回归模型中创建的OutlierTrans类。我们将使用它根据四分位数范围来识别异常值，正如我们在第三章的识别和修复缺失值中首先讨论的那样：

   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

3. 接下来，我们必须加载收入数据。我们还需要构建一个女性与男性收入比、教育年限比、劳动力参与比和人类发展指数比的序列。这些指标中的任何一项的值较低都表明男性可能具有优势，假设这些特征与女性与男性收入比之间存在正相关关系。例如，我们预计收入比会提高——也就是说，接近 1.0——当劳动力参与比接近 1.0 时——也就是说，当女性的劳动力参与率等于男性的时候。

4. 我们必须删除目标incomeratio缺失的行：

   un_income_gap = pd.read_csv("data/un_income_gap.csv")
   un_income_gap.set_index('country', inplace=True)
   un_income_gap['incomeratio'] = \
     un_income_gap.femaleincomepercapita / \
       un_income_gap.maleincomepercapita
   un_income_gap['educratio'] = \
     un_income_gap.femaleyearseducation / \
        un_income_gap.maleyearseducation
   un_income_gap['laborforcepartratio'] = \
     un_income_gap.femalelaborforceparticipation / \
        un_income_gap.malelaborforceparticipation
   un_income_gap['humandevratio'] = \
     un_income_gap.femalehumandevelopment / \
        un_income_gap.malehumandevelopment
   un_income_gap.dropna(subset=['incomeratio'], inplace=True)
   

5. 让我们看看一些数据行：

   num_cols = ['educratio','laborforcepartratio',
   'humandevratio','genderinequality','maternalmortality',
     'adolescentbirthrate','femaleperparliament',
   'incomepercapita']
   gap_sub = un_income_gap[['incomeratio'] + num_cols]
   gap_sub.head()
   incomeratio  educratio  laborforcepartratio  humandevratio\
   country
   Norway         0.78    1.02    0.89    1.00
   Australia      0.66    1.02    0.82    0.98
   Switzerland    0.64    0.88    0.83    0.95
   Denmark        0.70    1.01    0.88    0.98
   Netherlands    0.48    0.95    0.83    0.95
   genderinequality  maternalmortality  adolescentbirthrate\
   country
   Norway        0.07    4.00    7.80
   Australia     0.11    6.00    12.10
   Switzerland   0.03    6.00    1.90
   Denmark       0.05    5.00    5.10
   Netherlands   0.06    6.00    6.20
                femaleperparliament  incomepercapita  
   country                                            
   Norway       39.60    64992
   Australia    30.50    42261
   Switzerland  28.50    56431
   Denmark      38.00    44025
   Netherlands  36.90    45435
   

6. 让我们也看看一些描述性统计：

   gap_sub.\
     agg(['count','min','median','max']).T
                       count  min    median    max
   incomeratio         177.00 0.16   0.60      0.93
   educratio           169.00 0.24   0.93      1.35
   laborforcepartratio 177.00 0.19   0.75      1.04
   humandevratio       161.00 0.60   0.95      1.03
   genderinequality    155.00 0.02   0.39      0.74
   maternalmortality   174.00 1.00   60.00     1,100.00
   adolescentbirthrate 177.00 0.60   40.90     204.80
   femaleperparliament 174.00 0.00   19.35     57.50
   incomepercapita     177.00 581.00 10,512.00 123,124.00
   

我们有 177 个观测值和目标变量incomeratio。一些特征，如humandevratio和genderinequality，有超过 15 个缺失值。我们将在那里需要填充一些合理的值。我们还需要进行一些缩放，因为一些特征的范围与其他特征非常不同，从一端的incomeratio和incomepercapita到另一端的educratio和humandevratio。

注意

数据集有妇女和男人的单独人类发展指数。该指数是健康、知识获取和生活水平的衡量标准。我们之前计算的humandevratio特征将妇女的得分除以男人的得分。genderinequality特征是那些对妇女有不成比例影响的国家的健康和劳动力市场政策的指数。femaleperparliament是最高国家立法机构中女性所占的百分比。

1. 我们还应该查看特征与目标之间的相关性热图。在我们建模时，记住更高的相关性（无论是负的还是正的）是一个好主意。更高的正相关用较暖的颜色表示。laborforcepartratio、humandevratio和maternalmortality都与我们的目标正相关，后者有些令人惊讶。humandevratio和laborforcepartratio也是相关的，因此我们的模型可能难以区分每个特征的影响。一些特征选择将帮助我们弄清楚哪个特征更重要。（我们需要使用包装器或嵌入式特征选择方法来很好地揭示这一点。我们将在第五章，特征选择）中详细讨论这些方法。）看看以下代码：

   corrmatrix = gap_sub.corr(method="pearson")
   corrmatrix
   sns.heatmap(corrmatrix, xticklabels=corrmatrix.columns,
   yticklabels=corrmatrix.columns, cmap="coolwarm")
   plt.title('Heat Map of Correlation Matrix')
   plt.tight_layout()
   plt.show()
   

这会产生以下图表：

图 9.3 – 相关系数矩阵

1. 接下来，我们必须设置训练和测试数据框：

   X_train, X_test, y_train, y_test =  \
     train_test_split(gap_sub[num_cols],\
     gap_sub[['incomeratio']], test_size=0.2, random_state=0)
   

现在，我们已经准备好设置 KNN 回归模型。我们还将构建一个管道来处理异常值，基于每个特征的中间值进行插补，缩放特征，并使用 scikit-learn 的SelectFromModel进行一些特征选择：

1. 我们将使用线性回归进行特征选择，但我们可以选择任何会返回特征重要性值的算法。我们将特征重要性阈值设置为平均特征重要性的 80%。平均值是默认值。我们的选择在这里相当随意，但我喜欢保留那些低于平均特征重要性水平的特征的想法，当然还有那些重要性更高的特征：

   knnreg = KNeighborsRegressor()
   feature_sel = SelectFromModel(LinearRegression(), threshold="0.8*mean")
   pipe1 = make_pipeline(OutlierTrans(3), \
     SimpleImputer(strategy="median"), StandardScaler(), \
     feature_sel, knnreg)
   

2. 现在，我们已经准备好进行网格搜索以找到最佳参数。首先，我们将创建一个字典，knnreg_params，以表明我们希望 KNN 模型从 3 到 19 选择k的值，跳过偶数。我们还希望网格搜索找到最佳的距离度量 – 欧几里得、曼哈顿或闵可夫斯基：

   knnreg_params = {
    'kneighborsregressor__n_neighbors': \
        np.arange(3, 21, 2),
    'kneighborsregressor__metric': \
        ['euclidean','manhattan','minkowski']
   }
   

3. 我们将把这些参数传递给RandomizedSearchCV对象，然后拟合模型。我们可以使用RandomizedSearchCV的best_params_属性来查看特征选择和 KNN 回归选择的超参数。这些结果表明，对于 KNN 的k最佳超参数值是 11，对于距离度量是曼哈顿：

最佳模型具有-0.05 的负均方误差。考虑到样本量小，这相当不错。它小于incomeratio的中位数值的 10%，而incomeratio的中位数值为 0.6：

rs = RandomizedSearchCV(pipe1, knnreg_params, cv=4, n_iter=20, \
  scoring='neg_mean_absolute_error', random_state=1)
rs.fit(X_train, y_train)
rs.best_params_
{'kneighborsregressor__n_neighbors': 11,
 'kneighborsregressor__metric': 'manhattan'}
rs.best_score_
-0.05419731104389228

1. 让我们来看看在管道的特征选择步骤中选定的特征。只选择了两个特征 – laborforcepartratio 和 humandevratio。请注意，这一步不是运行我们的模型的必要步骤。它只是帮助我们解释它：

   selected = rs.best_estimator_['selectfrommodel'].get_support()
   np.array(num_cols)[selected]
   array(['laborforcepartratio', 'humandevratio'], dtype='<U19')
   

2. 如果你使用的是scikit-learn 1.0 或更高版本，这会容易一些。在这种情况下，你可以使用get_feature_names_out方法：

   rs.best_estimator_['selectfrommodel'].\
     get_feature_names_out(np.array(num_cols))
   array(['laborforcepartratio', 'humandevratio'], dtype=object)
   

3. 我们还应该看看一些其他的前几名结果。有一个使用euclidean距离的模型，其表现几乎与最佳模型一样好：

   results = \
     pd.DataFrame(rs.cv_results_['mean_test_score'], \
       columns=['meanscore']).\
     join(pd.DataFrame(rs.cv_results_['params'])).\
     sort_values(['meanscore'], ascending=False)
   results.head(3).T
   13       1      3
   Meanscore   -0.05   -0.05   -0.05
   regressor__kneighborsregressor__n_neighbors  11  13  9
   regressor__kneighborsregressor__metric  manhattan  manhattan  euclidean
   

4. 让我们看看这个模型的残差。我们可以使用RandomizedSearchCV对象的predict方法在测试数据上生成预测。残差在 0 附近很好地平衡。有一点负偏斜，但这也不算坏。峰度低，但我们对此没有太多尾巴的情况感到满意。这很可能反映了异常值残差不多：

   pred = rs.predict(X_test)
   preddf = pd.DataFrame(pred, columns=['prediction'],
     index=X_test.index).join(X_test).join(y_test)
   preddf['resid'] = preddf.incomeratio-preddf.prediction
   preddf.resid.agg(['mean','median','skew','kurtosis'])
   mean            -0.01
   median          -0.01
   skew            -0.61
   kurtosis         0.23
   Name: resid, dtype: float64
   

5. 让我们绘制残差图：

   plt.hist(preddf.resid, color="blue")
   plt.axvline(preddf.resid.mean(), color='red', linestyle='dashed', linewidth=1)
   plt.title("Histogram of Residuals for Gax Tax Model")
   plt.xlabel("Residuals")
   plt.ylabel("Frequency")
   plt.xlim()
   plt.show()
   

这产生了以下图表：

图 9.4 – 使用 KNN 回归的收入的比率模型的残差

当我们绘制残差时，它们看起来也很不错。然而，有几个国家，我们的预测误差超过了 0.1。在这两种情况下，我们都高估了。 （虚线红色线是平均残差量。）

1. 让我们也看看散点图。在这里，我们可以看到两个大的高估值位于预测范围的两端。总的来说，残差在预测的收入比率范围内相对恒定。我们可能只是想对两个异常值做些处理：

   plt.scatter(preddf.prediction, preddf.resid, color="blue")
   plt.axhline(0, color='red', linestyle='dashed', linewidth=1)
   plt.title("Scatterplot of Predictions and Residuals")
   plt.xlabel("Predicted Income Gap")
   plt.ylabel("Residuals")
   plt.show()
   

这产生了以下图表：

图 9.5 – 使用 KNN 回归的收入的比率模型的预测值和残差的散点图

我们应该更仔细地看看那些残差较高的国家。我们的模型在预测阿富汗或荷兰的收入比率方面做得不好，在两种情况下都高估了很多。回想一下，我们的特征选择步骤给我们提供了一个只有两个预测因子（laborforcepartratio 和 humandevratio）的模型。

对于阿富汗，劳动力参与率（女性相对于男性的参与率）非常接近最低的 0.19，而人类发展比率处于最低水平。这仍然不能使我们接近预测非常低的收入比率（女性相对于男性的收入），这个比率也是最低的。

对于荷兰来说，劳动力参与率 0.83 比中位数 0.75 高出不少，但人类发展比率正好是中位数。这就是为什么我们的模型预测的收入比率略高于 0.6 的中位数。荷兰的实际收入比率出人意料地低：

preddf.loc[np.abs(preddf.resid)>=0.1,
  ['incomeratio', 'prediction', 'resid', 
   'laborforcepartratio', 'humandevratio']].T
country                     Afghanistan    Netherlands
incomeratio                  0.16           0.48
prediction                   0.32           0.65
resid                       -0.16          -0.17
laborforcepartratio          0.20           0.83
humandevratio                0.60           0.95

在这里，我们可以看到 KNN 回归的一些优点。我们可以在不花费大量时间指定模型的情况下，对难以建模的数据获得不错的预测。除了某些插补和缩放之外，我们没有进行任何转换或创建交互效应。我们也不必担心非线性。KNN 回归可以很好地处理这一点。

但这种方法可能扩展得不是很好。在这个例子中，一个懒惰的学习者是可以的。然而，对于更工业级的工作，我们通常需要转向一个具有许多 KNN 优点但没有一些缺点的算法。我们将在本章的剩余部分探讨决策树和随机森林回归。

决策树和随机森林回归的关键概念

决策树是一个非常有用的机器学习工具。它们具有与 KNN 相似的一些优点——它们是非参数的，易于解释，并且可以处理广泛的数据——但没有一些限制。

决策树根据数据集中观测值的特征值对观测值进行分组。这是通过一系列二元决策来完成的，从根节点的一个初始分割开始，以每个分组的叶子节点结束。所有沿着从根节点到该叶子的分支具有相同值或相同值范围的观测值，都会得到相同的预测值。当目标是数值时，这就是该叶子节点训练观测值的平均目标值。图 9.6展示了这一点：

图 9.6 – 每晚睡眠小时数的决策树模型

这是一个基于每周工作时间、孩子数量、家庭中其他成年人数量以及个人是否在学校注册的每周睡眠小时数的个人模型。（这些结果基于假设数据。）根节点基于每周工作时间，将数据分为工作时间超过 30 小时和 30 小时或以下的观察结果。括号中的数字是该节点达到的训练数据百分比。60%的观察结果工作时间超过 30 小时。在树的左侧，我们的模型进一步通过孩子数量和家中其他成年人的数量来分割数据。在树的另一侧，代表工作时间少于或等于 30 小时的观察结果，唯一的额外分割是通过学校注册情况。

我现在意识到，并非所有读者都能看到这个颜色。我们可以从每个叶子节点向上导航树，描述树是如何分割数据的。15%的观测值在家中没有其他成年人，有超过 1 个孩子，每周工作时间超过 30 小时。这些观测值的平均每晚睡眠时间为 4.5 小时。这将是对具有相同特征的新观测值的预测值。

你可能想知道决策树算法是如何选择数值特征的阈值。例如，为什么每周工作时间大于 30 小时或孩子数量大于 1 小时？算法从根节点开始，在每个级别上选择分割，以最小化平方误差之和。更精确地说，选择的分割可以最小化：

你可能已经注意到了与线性回归优化的相似之处。但是，决策树回归相对于线性回归有几个优点。决策树可以用来模拟线性关系和非线性关系，而无需我们修改特征。我们还可以使用决策树避免特征缩放，因为算法可以处理我们特征中的非常不同的范围。

决策树的主要缺点是它们的方差很高。根据我们数据的特点，每次拟合决策树时，我们可能会得到一个非常不同的模型。我们可以使用集成方法，如袋装法或随机森林，来解决这个问题。

使用随机森林回归

随机森林，可能不出所料，是决策树的集合。但这并不能区分随机森林和自助聚合，通常称为袋装法。袋装法通常用于减少具有高方差的机器学习算法的方差，如决策树。在袋装法中，我们从数据集中生成随机样本。然后，我们在每个样本上运行我们的模型，例如决策树回归，并对预测进行平均。

然而，使用袋装法生成的样本可能相关，并且产生的决策树可能有很多相似之处。这种情况在只有少数特征可以解释大部分变化时更为可能。随机森林通过限制每个树可以选定的特征数量来解决此问题。一个很好的经验法则是将可用特征的数量除以 3，以确定每个决策树每个分割要使用的特征数量。例如，如果有 21 个特征，我们会在每个分割中使用 7 个。我们将在下一节中构建决策树和随机森林回归模型。

决策树和随机森林回归

在本节中，我们将使用决策树和随机森林构建一个回归模型，使用与本章前面相同的工作收入差距数据。 我们还将使用调整来识别给我们带来最佳性能模型的超参数，就像我们在 KNN 回归中所做的那样。 让我们开始吧：

1. 我们必须加载与 KNN 回归相同的许多库，以及来自 scikit-learn 的DecisionTreeRegressor和RandomForestRegressor：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.model_selection import RandomizedSearchCV
   from sklearn.tree import DecisionTreeRegressor, plot_tree
   from sklearn.ensemble import RandomForestRegressor
   from sklearn.linear_model import LinearRegression
   from sklearn.feature_selection import SelectFromModel
   

2. 我们还必须导入我们用于处理异常值的类：

   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

3. 我们必须加载之前使用过的相同收入差距数据，并创建测试和训练数据框：

   un_income_gap = pd.read_csv("data/un_income_gap.csv")
   un_income_gap.set_index('country', inplace=True)
   un_income_gap['incomeratio'] = \
     un_income_gap.femaleincomepercapita / \
       un_income_gap.maleincomepercapita
   un_income_gap['educratio'] = \
     un_income_gap.femaleyearseducation / \
        un_income_gap.maleyearseducation
   un_income_gap['laborforcepartratio'] = \
     un_income_gap.femalelaborforceparticipation / \
        un_income_gap.malelaborforceparticipation
   un_income_gap['humandevratio'] = \
     un_income_gap.femalehumandevelopment / \
        un_income_gap.malehumandevelopment
   un_income_gap.dropna(subset=['incomeratio'], 
     inplace=True)
   num_cols = ['educratio','laborforcepartratio',
     'humandevratio', 'genderinequality', 
     'maternalmortality', 'adolescentbirthrate', 
     'femaleperparliament', 'incomepercapita']
   gap_sub = un_income_gap[['incomeratio'] + num_cols]
   X_train, X_test, y_train, y_test =  \
     train_test_split(gap_sub[num_cols],\
     gap_sub[['incomeratio']], test_size=0.2, 
       random_state=0)
   

让我们从相对简单的决策树开始 – 一个没有太多层级的树。 一个简单的树可以很容易地在一页上展示。

带有解释的决策树示例

在我们构建决策树回归器之前，让我们先看看一个最大深度设置为低值的快速示例。 随着深度的增加，决策树解释和绘图变得更加困难。 让我们开始吧：

1. 我们首先实例化一个决策树回归器，限制深度为三，并要求每个叶节点至少有五个观测值。 我们创建一个仅预处理数据的管道，并将结果 NumPy 数组X_train_imp传递给决策树回归器的fit方法：

   dtreg_example = DecisionTreeRegressor(
     min_samples_leaf=5,
     max_depth=3)
   pipe0 = make_pipeline(OutlierTrans(3),
     SimpleImputer(strategy="median"))
   X_train_imp = pipe0.fit_transform(X_train)
   dtreg_example.fit(X_train_imp, y_train)
   plot_tree(dtreg_example, 
     feature_names=X_train.columns,
     label="root", fontsize=10)
   

这生成了以下图表：

图 9.7 – 最大深度为 3 的决策树示例

我们不会遍历这棵树上的所有节点。 通过描述通往几个叶节点的路径，我们可以了解如何解释决策树回归图的一般方法：

• 解释劳动力参与率 <= 0.307 的叶节点：

根节点分割是基于劳动力参与率小于或等于 0.601 的。 (回想一下，劳动力参与率是女性参与率与男性参与率的比率。) 有 34 个国家属于这一类别。 (分割测试的真实值在左侧。虚假值在右侧。) 之后还有另一个基于劳动力参与率的分割，这次分割点为 0.378。 有 13 个国家的值小于或等于这个值。 最后，我们到达了劳动力参与率小于或等于 0.307 的最左侧的叶节点。 有六个国家的劳动力参与率如此之低。 这六个国家的平均收入比率为 0.197。 然后，我们的决策树回归器将预测劳动力参与率小于或等于 0.307 的测试实例的收入比率为 0.197。

• 解释劳动力参与率在 0.601 到 0.811 之间，且 humandevratio <= 0.968 的叶节点：

有 107 个国家的劳动力参与率大于 0.601。这显示在树的右侧。当劳动力参与率小于或等于 0.811 时，还有一个二进制分割，进一步基于人类发展比率小于或等于 0.968 进行分割。这带我们到一个有 31 个国家的叶子节点，这些国家的人类发展比率小于或等于 0.968，劳动力参与率小于或等于 0.811，但大于 0.601。决策树回归器将预测这 31 个国家的收入比率平均值，为 0.556，对于所有测试实例，其人类发展比率和劳动力参与率的值都在这些范围内。

有趣的是，我们还没有进行任何特征选择，但这次构建决策树模型的初步尝试已经表明，仅用两个特征就可以预测收入比率：laborforcepartratio和humandevratio。

尽管这个模型的简单性使得它非常容易解释，但我们还没有完成寻找最佳超参数所需的工作。让我们接下来做这件事。

构建和解释我们的实际模型

按照以下步骤进行：

1. 首先，我们实例化一个新的决策树回归器并创建一个使用它的管道。我们还为一些超参数创建了一个字典——也就是说，对于最大树深度和每个叶子的最小样本数（观察值）。请注意，我们不需要对我们的特征或目标进行缩放，因为在决策树中这不是必要的：

   dtreg = DecisionTreeRegressor()
   feature_sel = SelectFromModel(LinearRegression(),
     threshold="0.8*mean")
   pipe1 = make_pipeline(OutlierTrans(3),
     SimpleImputer(strategy="median"),
     feature_sel, dtreg)
   dtreg_params={
    "decisiontreeregressor__max_depth": np.arange(2, 20),
    "decisiontreeregressor__min_samples_leaf": np.arange(5, 11)
   }
   

2. 接下来，我们必须根据上一步的字典设置一个随机搜索。我们的决策树的最佳参数是 5 个最小样本和 9 的最大深度：

   rs = RandomizedSearchCV(pipe1, dtreg_params, cv=4, n_iter=20,
     scoring='neg_mean_absolute_error', random_state=1)
   rs.fit(X_train, y_train.values.ravel())
   rs.best_params_
   {'decisiontreeregressor__min_samples_leaf': 5,
    'decisiontreeregressor__max_depth': 9}
   rs.best_score_
   -0.05268976358459662
   

正如我们在上一节中讨论的，决策树在回归方面具有许多 KNN 的优点。它们容易解释，并且对底层数据没有太多假设。然而，决策树仍然可以与大型数据集合理地工作。决策树的一个不那么重要但仍然有用的优点是，它们不需要特征缩放。

但决策树确实有高方差。通常值得牺牲决策树的可解释性以换取相关方法，例如随机森林，这可以显著减少方差。我们在上一节中概念性地讨论了随机森林算法。我们将在下一节中使用收入差距数据尝试它。

随机森林回归

回想一下，随机森林可以被视为具有袋装法的决策树；它们通过减少样本之间的相关性来改进袋装法。这听起来很复杂，但它的实现与决策树一样简单。让我们看看：

1. 我们将首先实例化一个随机森林回归器并为超参数创建一个字典。我们还将为预处理和回归器创建一个管道：

   rfreg = RandomForestRegressor()
   rfreg_params = {
    'randomforestregressor__max_depth': np.arange(2, 20),
    'randomforestregressor__max_features': ['auto', 'sqrt'],
    'randomforestregressor__min_samples_leaf':  np.arange(5, 11)
   }
   pipe2 = make_pipeline(OutlierTrans(3), 
     SimpleImputer(strategy="median"),
     feature_sel, rfreg)
   

2. 我们将把管道和超参数字典传递给RandomizedSearchCV对象以运行网格搜索。在得分方面有轻微的改进：

   rs = RandomizedSearchCV(pipe2, rfreg_params, cv=4, n_iter=20,
     scoring='neg_mean_absolute_error', random_state=1)
   rs.fit(X_train, y_train.values.ravel())
   rs.best_params_
   {'randomforestregressor__min_samples_leaf': 5,
    'randomforestregressor__max_features': 'auto',
    'randomforestregressor__max_depth': 9}
   rs.best_score_
   -0.04930503752638253
   

3. 让我们看看残差：

   pred = rs.predict(X_test)
   preddf = pd.DataFrame(pred, columns=['prediction'],
     index=X_test.index).join(X_test).join(y_test)
   preddf['resid'] = preddf.incomegap-preddf.prediction
   plt.hist(preddf.resid, color="blue", bins=5)
   plt.axvline(preddf.resid.mean(), color='red', linestyle='dashed', linewidth=1)
   plt.title("Histogram of Residuals for Income Gap")
   plt.xlabel("Residuals")
   plt.ylabel("Frequency")
   plt.xlim()
   plt.show()
   

这会产生以下图表：

图 9.8 – 随机森林模型在收入比率上的残差直方图

1. 让我们也看看残差与预测的散点图：

   plt.scatter(preddf.prediction, preddf.resid, color="blue")
   plt.axhline(0, color='red', linestyle='dashed', linewidth=1)
   plt.title("Scatterplot of Predictions and Residuals")
   plt.xlabel("Predicted Income Gap")
   plt.ylabel("Residuals")
   plt.show()
   

这会产生以下图表：

图 9.9 – 随机森林模型在收入比率上的预测与残差散点图

1. 让我们仔细看看一个显著的异常值，我们在这里严重高估了：

   preddf.loc[np.abs(preddf.resid)>=0.12,
     ['incomeratio','prediction','resid',
     'laborforcepartratio', 'humandevratio']].T
   country              Netherlands
   incomeratio                 0.48
   prediction                  0.66
   resid                      -0.18
   laborforcepartratio         0.83
   humandevratio               0.95
   

我们在荷兰仍然有困难，但残差的相对均匀分布表明这是异常的。实际上，这是一个好消息，从我们预测新实例收入比率的能力来看，表明我们的模型并没有过于努力地拟合这个不寻常的案例。

使用梯度提升回归

我们有时可以通过使用梯度提升来改进随机森林模型。与随机森林类似，梯度提升是一种集成方法，它结合了学习者，通常是树。但与随机森林不同，每棵树都是根据前一棵树的错误来学习的。这可以显著提高我们建模复杂性的能力。

虽然梯度提升不太容易过拟合，但我们必须比随机森林模型更加小心地进行超参数调整。我们可以减慢学习率，也称为收缩。我们还可以调整估计器的数量（树的数量）。学习率的选择影响所需估计器的数量。通常，如果我们减慢学习率，我们的模型将需要更多的估计器。

有几种工具可以用于实现梯度提升。我们将使用其中的两个：来自 scikit-learn 的梯度提升回归和 XGBoost。

在本节中，我们将处理房价数据。我们将尝试根据房屋及其附近房屋的特征，预测华盛顿州金斯县的房价。

注意

该数据集关于金斯县的房价可以由公众在www.kaggle.com/datasets/harlfoxem/housesalesprediction下载。它包含多个卧室、浴室和楼层，房屋和地块的平方英尺，房屋状况，15 个最近房屋的平方英尺，以及更多作为特征。

让我们开始构建模型：

1. 我们将首先导入所需的模块。其中两个新模块来自 XGBoost，分别是GradientBoostingRegressor和XGBRegressor：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.impute import SimpleImputer
   from sklearn.pipeline import make_pipeline
   from sklearn.preprocessing import OneHotEncoder
   from sklearn.preprocessing import MinMaxScaler
   from sklearn.compose import ColumnTransformer
   from sklearn.model_selection import RandomizedSearchCV
   from sklearn.ensemble import GradientBoostingRegressor
   from xgboost import XGBRegressor
   from sklearn.linear_model import LinearRegression
   from sklearn.feature_selection import SelectFromModel
   import matplotlib.pyplot as plt
   from scipy.stats import randint
   from scipy.stats import uniform
   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

2. 让我们加载房价数据并查看一些实例：

   housing = pd.read_csv("data/kc_house_data.csv")
   housing.set_index('id', inplace=True)
   num_cols = ['bedrooms', 'bathrooms', 'sqft_living', 
     'sqft_lot', 'floors', 'view', 'condition', 
     'sqft_above', 'sqft_basement', 'yr_built', 
     'yr_renovated', 'sqft_living15', 'sqft_lot15']
   cat_cols = ['waterfront']
   housing[['price'] + num_cols + cat_cols].\
     head(3).T
   id              7129300520  6414100192  5631500400
   price           221,900     538,000     180,000
   bedrooms        3           3           2
   bathrooms       1           2           1
   sqft_living     1,180       2,570       770
   sqft_lot        5,650       7,242       10,000
   floors          1           2           1
   view            0           0           0
   condition       3           3           3
   sqft_above      1,180       2,170       770
   sqft_basement   0           400         0
   yr_built        1,955       1,951       1,933
   yr_renovated    0           1,991       0
   sqft_living15   1,340       1,690       2,720
   sqft_lot15      5,650       7,639       8,062
   waterfront      0           0           0
   

3. 我们还应该查看一些描述性统计。我们没有缺失值。我们的目标变量price有一些极端值，不出所料。这可能会在建模中引起问题。我们还需要处理一些特征的极端值：

   housing[['price'] + num_cols].\
     agg(['count','min','median','max']).T
                   count   min      median    max
   price          21,613   75,000   450,000   7,700,000
   bedrooms       21,613   0        3         33
   bathrooms      21,613   0        2         8
   sqft_living    21,613   290      1,910     13,540
   sqft_lot       21,613   520      7,618     1,651,359
   floors         21,613   1        2         4
   view           21,613   0        0         4
   condition      21,613   1        3         5
   sqft_above     21,613   290      1,560     9,410
   sqft_basement  21,613   0        0         4,820
   yr_built       21,613   1,900    1,975     2,015
   yr_renovated   21,613   0        0         2,015
   sqft_living15  21,613   399      1,840     6,210
   sqft_lot15     21,613   651      7,620     871,200
   

4. 让我们创建一个房价直方图：

   plt.hist(housing.price/1000)
   plt.title("Housing Price (in thousands)")
   plt.xlabel('Price')
   plt.ylabel("Frequency")
   plt.show()
   

这生成了以下图表：

图 9.10 – 房价直方图

1. 如果我们使用目标变量的对数变换来进行建模，可能会更有运气，正如我们在第四章**，编码、转换和缩放特征中尝试的那样，使用 COVID 总病例数据。

   housing['price_log'] = np.log(housing['price'])
   plt.hist(housing.price_log)
   plt.title("Housing Price Log")
   plt.xlabel('Price Log')
   plt.ylabel("Frequency")
   plt.show()
   

这产生了以下图表：

图 9.11 – 房价对数直方图

1. 这看起来更好。让我们看看价格和价格对数的偏度和峰度。对数看起来有很大的改进：

   housing[['price','price_log']].agg(['kurtosis','skew'])
                    price       price_log
   kurtosis         34.59            0.69
   skew              4.02            0.43
   

2. 我们还应该查看一些相关性。居住面积平方英尺、地面以上面积平方英尺、最近 15 个家庭的居住面积平方英尺以及浴室数量是与价格最相关的特征。居住面积平方英尺和地面以上居住面积平方英尺高度相关。我们可能需要在模型中在这两者之间做出选择：

   corrmatrix = housing[['price_log'] + num_cols].\
      corr(method="pearson")
   sns.heatmap(corrmatrix, 
     xticklabels=corrmatrix.columns,
     yticklabels=corrmatrix.columns, cmap="coolwarm")
   plt.title('Heat Map of Correlation Matrix')
   plt.tight_layout()
   plt.show()
   

这产生了以下图表：

图 9.12 – 房地产特征的相关矩阵

1. 接下来，我们创建训练和测试数据框：

   target = housing[['price_log']]
   features = housing[num_cols + cat_cols]
   X_train, X_test, y_train, y_test =  \
     train_test_split(features,\
     target, test_size=0.2, random_state=0)
   

2. 我们还需要设置我们的列转换。对于所有数值特征，即除了waterfront之外的所有特征，我们将检查极端值，然后缩放数据：

   ohe = OneHotEncoder(drop='first', sparse=False)
   standtrans = make_pipeline(OutlierTrans(2),
     SimpleImputer(strategy="median"),
     MinMaxScaler())
   cattrans = make_pipeline(ohe)
   coltrans = ColumnTransformer(
     transformers=[
       ("stand", standtrans, num_cols),
       ("cat", cattrans, cat_cols)
     ]
   )
   

3. 现在，我们已准备好设置预处理和模型的管道。我们将实例化一个GradientBoostingRegressor对象并设置特征选择。我们还将创建一个超参数字典，用于我们在下一步中进行的随机网格搜索：

   gbr = GradientBoostingRegressor(random_state=0)
   feature_sel = SelectFromModel(LinearRegression(),
     threshold="0.6*mean")
   gbr_params = {
    'gradientboostingregressor__learning_rate': uniform(loc=0.01, scale=0.5),
    'gradientboostingregressor__n_estimators': randint(500, 2000),
    'gradientboostingregressor__max_depth': randint(2, 20),
    'gradientboostingregressor__min_samples_leaf': randint(5, 11)
   }
   pipe1 = make_pipeline(coltrans, feature_sel, gbr)
   

4. 现在，我们已准备好运行随机网格搜索。考虑到price_log的平均值约为 13，我们得到了相当不错的均方误差分数：

   rs1 = RandomizedSearchCV(pipe1, gbr_params, cv=5, n_iter=20,
     scoring='neg_mean_squared_error', random_state=0)
   rs1.fit(X_train, y_train.values.ravel())
   rs1.best_params_
   {'gradientboostingregressor__learning_rate': 0.118275177212,
    'gradientboostingregressor__max_depth': 2,
    'gradientboostingregressor__min_samples_leaf': 5,
    'gradientboostingregressor__n_estimators': 1577}
   rs1.best_score_
   -0.10695077555421204
   y_test.mean()
   price_log   13.03
   dtype: float64
   

5. 不幸的是，平均拟合时间相当长：

   print("fit time: %.3f, score time: %.3f"  %
     (np.mean(rs1.cv_results_['mean_fit_time']),\
     np.mean(rs1.cv_results_['mean_score_time'])))
   fit time: 35.695, score time: 0.152
   

6. 让我们尝试使用 XGBoost：

   xgb = XGBRegressor()
   xgb_params = {
    'xgbregressor__learning_rate': uniform(loc=0.01, scale=0.5),
    'xgbregressor__n_estimators': randint(500, 2000),
    'xgbregressor__max_depth': randint(2, 20)
   }
   pipe2 = make_pipeline(coltrans, feature_sel, xgb)
   

7. 我们没有获得更好的分数，但平均拟合时间和分数时间有了显著提高：

   rs2 = RandomizedSearchCV(pipe2, xgb_params, cv=5, n_iter=20,
     scoring='neg_mean_squared_error', random_state=0)
   rs2.fit(X_train, y_train.values.ravel())
   rs2.best_params_
   {'xgbregressor__learning_rate': 0.019394900218177573,
    'xgbregressor__max_depth': 7,
    'xgbregressor__n_estimators': 1256}
   rs2.best_score_
   -0.10574300757906044
   print("fit time: %.3f, score time: %.3f"  %
     (np.mean(rs2.cv_results_['mean_fit_time']),\
     np.mean(rs2.cv_results_['mean_score_time'])))
   fit time: 3.931, score time: 0.046
   

XGBoost 由于许多原因已经成为一个非常受欢迎的梯度提升工具，其中一些你已经在本例中看到了。它可以快速产生非常好的结果，而且模型指定很少。我们确实需要仔细调整我们的超参数以获得首选的方差-偏差权衡，但这也适用于其他梯度提升工具，正如我们所看到的。

摘要

在本章中，我们探讨了最受欢迎的一些非参数回归算法：K 最近邻算法、决策树和随机森林。使用这些算法构建的模型可以表现良好，但也存在一些限制。我们讨论了这些技术的优缺点，包括维度和观测限制，以及对于 KNN 模型所需训练时间的担忧。我们讨论了决策树的关键挑战，即高方差，以及随机森林模型如何解决这一问题。我们还探讨了梯度提升回归树，并讨论了它们的一些优点。我们继续提高我们在超参数调整方面的技能，因为每个算法都需要一种不同的策略。

在接下来的几章中，我们将讨论监督学习算法，其中目标变量是分类的。首先，我们将从可能最熟悉的分类算法——逻辑回归开始讨论。

第四部分 - 使用监督学习建模二分类和多类目标

对于预测分类目标，有许多高性能算法。在本部分，我们将考察最流行的分类算法。我们还将考虑为什么根据我们的数据和领域知识，我们可能会选择一种算法而不是其他算法。

我们对分类模型中的欠拟合和过拟合的关注程度与之前部分中对回归模型的关注程度一样。当特征与目标之间的关系复杂时，我们需要使用能够捕捉这种复杂性的算法。但往往存在非同小可的过拟合风险。在本部分章节中，我们将讨论建模复杂性而不出现过拟合的策略。这通常涉及对逻辑回归模型进行某种形式的正则化、对决策树的树深度进行限制，以及调整支持向量分类中边缘违规的容忍度。

如果我们试图建模复杂性而不出现过拟合，我们必须准备好花大量时间进行超参数调整。在这些章节中，我们肯定会在这方面花费相当多的时间。与此相关，我们也会在交叉验证、生成和解释评估指标方面变得非常熟练。在接下来的五个章节中，我们将讨论准确率、精确度、敏感度和特异性。我们还将非常习惯于查看混淆矩阵。

我们还将探讨如何将这些算法扩展到多类目标。对于 k-最近邻和决策树来说，这是直截了当的，但对于逻辑回归和支持向量回归算法则需要扩展。这些内容将在本章中介绍。

本节包括以下章节：

• 第十章, 逻辑回归

• 第十一章, 决策树和随机森林分类

• 第十二章, K-最近邻分类

• 第十三章, 支持向量机分类

• 第十四章, 朴素贝叶斯分类

第十章：第十章：逻辑回归

在本章和接下来的几章中，我们将探讨分类模型。这些模型涉及具有两个或多个类别值的目标，例如学生是否会通过一门课程，或者顾客在只有鸡肉、牛肉和豆腐三种选择的情况下会选择哪一种。对于这类分类问题，存在几种机器学习算法。在本章中，我们将查看其中一些最受欢迎的算法。

逻辑回归已经用于构建具有二元目标的模型数十年了。传统上，它被用来生成独立变量或变量对二元结果概率影响估计。由于我们的重点是预测，而不是每个特征的影响，我们还将探讨正则化技术，如 lasso 回归。这些技术可以提高我们的分类预测准确性。我们还将检查预测多类别目标（当存在超过两个可能的目标值时）的策略。

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

• 逻辑回归的关键概念

• 使用逻辑回归进行二元分类

• 使用逻辑回归进行正则化

• 多项式逻辑回归

技术要求

在本章中，我们将坚持使用在大多数 Python 科学发行版中可用的库：pandas、NumPy 和 scikit-learn。本章中的所有代码都可以在 scikit-learn 版本 0.24.2 和 1.0.2 上正常运行。

逻辑回归的关键概念

如果你熟悉线性回归，或者阅读了本书的第七章，线性回归模型，你可能会预见到我们将在本章讨论的一些问题——正则化、回归器的线性关系和正态分布的残差。如果你过去构建过监督机器学习模型，或者阅读过本书的最后一章，那么你很可能预见到我们将花一些时间讨论偏差-方差权衡以及它如何影响我们选择模型。

我记得 35 年前在一门大学课程中第一次接触到逻辑回归。在本科教科书中，它通常几乎被呈现为线性回归的一个特例；也就是说，具有二元因变量的线性回归，并伴随一些变换以保持预测值在 0 到 1 之间。

它确实与数值目标变量的线性回归有许多相似之处。逻辑回归相对容易训练和解释。线性回归和逻辑回归的优化技术都是高效的，可以生成低偏差的预测器。

同样地，与线性回归一样，逻辑回归也是基于分配给每个特征的权重来预测目标。但为了将预测概率约束在 0 到 1 之间，我们使用 sigmoid 函数。这个函数将任何值映射到 0 到 1 之间的值：

当x趋近于无穷大时，趋近于 1。当x趋近于负无穷大时，趋近于 0。

下面的图示说明了 sigmoid 函数：

图 10.1 – Sigmoid 函数

我们可以将线性回归的熟悉方程代入 sigmoid 函数来预测类成员的概率：

在这里，是二元情况下类成员的预测概率。系数（β）可以转换为优势比以进行解释，如下所示：

在这里，r是优势比，β是系数。特征值增加 1 个单位会乘以类成员的优势比。同样，对于二元特征，一个真值有倍类成员的优势比，而一个假值也有相同特征的优势比，其他条件相同。

逻辑回归作为分类问题的算法具有几个优点。特征可以是二元的、分类的或数值的，且不需要服从正态分布。目标变量可以具有超过两个的可能值，正如我们稍后将要讨论的，它可以是无序的或有序的。另一个关键优点是，特征与目标之间的关系不假设是线性的。

这里的命名有点令人困惑。为什么我们使用回归算法来解决分类问题？好吧，逻辑回归预测类成员的概率。我们应用决策规则来预测这些概率。对于二元目标，默认阈值通常是 0.5；预测概率大于或等于 0.5 的实例被赋予正类或 1 或 True；那些小于 0.5 的实例被分配 0 或 False。

逻辑回归的扩展

在本章中，我们将考虑逻辑回归的两个关键扩展。我们将探讨多类模型——即目标值超过两个的模型。我们还将检查逻辑模型的正则化以改善（减少）方差。

在构建多类模型时，多项式逻辑回归（MLR）是一个流行的选择。使用 MLR，预测概率分布是一个多项式概率分布。我们可以用 softmax 函数替换我们用于二元分类器的方程：

在这里，。这为每个类别标签j计算一个概率，其中k是类别的数量。

当我们拥有超过两个类别时，一对余（OVR）逻辑回归是多项式逻辑回归的一个替代方案。这种逻辑回归的扩展将多类别问题转化为二分类问题，估计类别成员资格相对于所有其他类别的成员资格的概率。这里的关键假设是每个类别的成员资格是独立的。在本章的例子中，我们将使用 MLR。它相较于 OVR 的一个优点是预测概率更加可靠。

如前所述，逻辑回归与线性回归有一些相同的挑战，包括我们的预测的低偏差伴随着高方差。当几个特征高度相关时，这更有可能成为一个问题。幸运的是，我们可以通过正则化来解决这个问题，就像我们在第七章，“线性回归模型”中看到的那样。

正则化会给损失函数添加一个惩罚项。我们仍然寻求最小化误差，但同时也约束参数的大小。L1正则化，也称为 lasso 回归，惩罚权重（或系数）的绝对值：

在这里，p是特征的数量，λ决定了正则化的强度。L2正则化，也称为岭回归，惩罚权重（或系数）的平方值：

L1 和 L2 正则化都将权重推向 0，尽管 L1 正则化更有可能导致稀疏模型。在 scikit-learn 中，我们使用C参数来调整λ的值，其中C只是λ的倒数：

我们可以通过弹性网络回归在 L1 和 L2 之间取得平衡。在弹性网络回归中，我们调整 L1 比率。0.5 的值表示 L1 和 L2 同等使用。我们可以使用超参数调整来选择 L1 比率的最佳值。

正则化可能导致具有更低方差模型的产生，当我们对预测的关注超过对系数的关注时，这是一个很好的权衡。

在构建具有正则化的模型之前，我们将构建一个相当简单的具有二元目标的逻辑模型。我们还将花大量时间评估该模型。这将是本书中我们将构建的第一个分类模型，并且模型评估对于这些模型与回归模型看起来非常不同。

逻辑回归的二分类

当目标为二元时，逻辑回归常用于建模健康结果，例如，一个人是否患有疾病。在本节中，我们将通过一个例子来展示这一点。我们将构建一个模型，根据个人的吸烟和饮酒习惯、健康特征（包括 BMI、哮喘、糖尿病和皮肤癌）以及年龄等个人特征来预测一个人是否会患有心脏病。

注意

在本章中，我们将专门使用可在www.kaggle.com/datasets/kamilpytlak/personal-key-indicators-of-heart-disease公开下载的心脏病数据。这个数据集来源于 2020 年美国疾病控制与预防中心超过 40 万个人的数据。数据列包括受访者是否曾经患有心脏病、体重指数、是否吸烟、大量饮酒、年龄、糖尿病和肾病。在本节中，我们将使用 30,000 个个体样本以加快处理速度，但完整的数据集可以在本书 GitHub 仓库的同一文件夹中找到。

在本章中，我们将比以前章节进行更多的预处理。我们将把大部分工作整合到我们的管道中。这将使将来重用此代码更容易，并减少数据泄露的可能性。请按照以下步骤操作：

1. 我们将首先导入我们在过去几章中使用的相同库。我们还将导入 LogisticRegression 和 metrics 模块。我们将使用 scikit-learn 的 metrics 模块来评估本书这一部分中的每个分类模型。除了 matplotlib 用于可视化外，我们还将使用 seaborn：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.preprocessing import OneHotEncoder
   from sklearn.pipeline import make_pipeline
   from sklearn.impute import SimpleImputer
   from sklearn.compose import ColumnTransformer
   from sklearn.model_selection import StratifiedKFold
   from sklearn.feature_selection import RFECV
   from sklearn.linear_model import LogisticRegression
   import sklearn.metrics as skmet
   import matplotlib.pyplot as plt
   import seaborn as sns
   

2. 我们还需要几个自定义类来处理预处理。我们已经看到了 OutlierTrans 类。在这里，我们添加了两个新的类——MakeOrdinal 和 ReplaceVals：

   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans,\
     MakeOrdinal, ReplaceVals
   

MakeOrdinal 类接受一个字符特征并根据字母数字排序分配数值。例如，一个有三个可能值（不好、一般、好）的特征将被转换为一个序数特征，其值分别为 0、1 和 2。

请记住，scikit-learn 管道转换器必须具有 fit 和 transform 方法，并且必须继承自 BaseEstimator。它们通常也继承自 TransformerMixin，尽管还有其他选项。

MakeOrdinal 类中的所有操作都在 transform 方法中完成。我们遍历通过列转换器传递给它的所有列。对于每一列，我们找到所有唯一的值并按字母数字顺序排序，将唯一的值存储在我们命名为 cats 的 NumPy 数组中。然后，我们使用 lambda 函数和 NumPy 的 where 方法来找到与每个特征值关联的 cats 的索引：

class MakeOrdinal(BaseEstimator,TransformerMixin):
  def fit(self,X,y=None):
    return self

  def transform(self,X,y=None):
    Xnew = X.copy()
    for col in Xnew.columns:
      cats = np.sort(Xnew[col].unique())
      Xnew[col] = Xnew.\
        apply(lambda x: int(np.where(cats==\
        x[col])[0]), axis=1)
    return Xnew.values

当字母数字顺序与一个有意义的顺序相匹配时，MakeOrdinal 将正常工作，就像前面的例子一样。当这不是真的时，我们可以使用 ReplaceVals 来分配适当的序数值。这个类根据传递给它的字典替换任何特征中的值。

我们本可以使用 pandas 的 replace 方法而不将其放入管道中，但这样更容易将我们的重新编码与其他管道步骤（如特征缩放）集成：

class ReplaceVals(BaseEstimator,TransformerMixin):
  def __init__(self,repdict):
    self.repdict = repdict
  def fit(self,X,y=None):
    return self

  def transform(self,X,y=None):
    Xnew = X.copy().replace(self.repdict)
return Xnew.values

如果你现在不完全理解我们将如何使用这些类，请不要担心。当我们将它们添加到列转换中时，一切都会变得清晰。

1. 接下来，我们将加载心脏病数据并查看几行。一些字符串特征在概念上是二元的，例如alcoholdrinkingheavy，当一个人是重度饮酒者时为Yes，否则为No。在运行模型之前，我们需要对这些特征进行编码。

agecategory特征是表示年龄区间的字符数据。我们需要将这个特征转换为数值：

healthinfo = pd.read_csv("data/healthinfo.csv")
healthinfo.set_index("personid", inplace=True)
healthinfo.head(2).T
personid                    299391       252786
heartdisease                Yes          No
bmi                         28.48        25.24
smoking                     Yes          Yes
alcoholdrinkingheavy        No           No
stroke                      No           No
physicalhealthbaddays       7            0
mentalhealthbaddays         0            2
walkingdifficult            No           No
gender                      Male         Female
agecategory                 70-74        65-69
ethnicity                   White        White
diabetic    No, borderline diabetes      No
physicalactivity            Yes          Yes
genhealth                   Good         Very good
sleeptimenightly            8            8
asthma                      No           No
kidneydisease               No           No
skincancer                  No           Yes

1. 让我们看看 DataFrame 的大小以及有多少缺失值。有 30,000 个实例，但 18 个数据列中没有任何缺失值。这太好了。当我们构建管道时，我们不必担心这一点：

   healthinfo.shape
   (30000, 18)
   healthinfo.isnull().sum()
   heartdisease                0
   bmi                         0
   smoking                     0
   alcoholdrinkingheavy        0
   stroke                      0
   physicalhealthbaddays       0
   mentalhealthbaddays         0
   walkingdifficult            0
   gender                      0
   agecategory                 0
   ethnicity                   0
   diabetic                    0
   physicalactivity            0
   genhealth                   0
   sleeptimenightly            0
   asthma                      0
   kidneydisease               0
   skincancer                  0
   dtype: int64
   

2. 让我们将heartdisease变量，也就是我们的目标变量，转换为0和1变量。这将减少我们以后需要担心的事情。立即要注意的一件事是，目标变量的值非常不平衡。我们观察结果中不到 10%的人患有心脏病。这当然是好消息，但也给我们的建模带来了一些挑战，我们需要处理这些挑战：

   healthinfo.heartdisease.value_counts()
   No        27467
   Yes       2533
   Name: heartdisease, dtype: int64
   healthinfo['heartdisease'] = \
     np.where(healthinfo.heartdisease=='No',0,1).\
        astype('int')
   healthinfo.heartdisease.value_counts()
   0        27467
   1        2533
   Name: heartdisease, dtype: int64
   

3. 我们应该根据我们将要对它们进行的预处理来组织我们的特征。我们将对数值特征进行缩放，并对分类特征进行独热编码。我们希望将当前为字符串的agecategory和genhealth特征转换为有序特征。

我们需要对diabetic特征进行特定的清理。有些人表示没有，但他们处于边缘状态。为了我们的目的，我们将它们视为no。有些人怀孕期间只有糖尿病。我们将它们视为yes。对于genhealth和diabetic，我们将设置一个字典，以指示特征值应该如何替换。我们将在管道的ReplaceVals转换器中使用该字典：

num_cols = ['bmi','physicalhealthbaddays',
   'mentalhealthbaddays','sleeptimenightly']
binary_cols = ['smoking','alcoholdrinkingheavy',
  'stroke','walkingdifficult','physicalactivity',
  'asthma','kidneydisease','skincancer']
cat_cols = ['gender','ethnicity']
spec_cols1 = ['agecategory']
spec_cols2 = ['genhealth']
spec_cols3 = ['diabetic']
rep_dict = {
  'genhealth': {'Poor':0,'Fair':1,'Good':2,
    'Very good':3,'Excellent':4},
  'diabetic': {'No':0,
    'No, borderline diabetes':0,'Yes':1,
    'Yes (during pregnancy)':1}           
}

1. 我们应该查看一些二元特征以及其他分类特征的频率。很大比例的个人（42%）报告说他们是吸烟者。14%的人报告说他们走路有困难：

   healthinfo[binary_cols].\
     apply(pd.value_counts, normalize=True).T
                           No       Yes
   smoking                 0.58     0.42
   alcoholdrinkingheavy    0.93     0.07
   stroke                  0.96     0.04
   walkingdifficult        0.86     0.14
   physicalactivity        0.23     0.77
   asthma                  0.87     0.13
   kidneydisease           0.96     0.04
   skincancer              0.91     0.09
   

2. 让我们也看看其他分类特征的频率。男性和女性的数量几乎相等。大多数人报告他们的健康状况非常好或非常好：

   for col in healthinfo[cat_cols + 
   ['genhealth','diabetic']].columns:
     print(col, "----------------------",
     healthinfo[col].value_counts(normalize=True).\
         sort_index(), sep="\n", end="\n\n")
   

这会产生以下输出：

gender
----------------------
Female   0.52
Male     0.48
Name: gender, dtype: float64
ethnicity
----------------------
American Indian/Alaskan Native   0.02
Asian                            0.03
Black                            0.07
Hispanic                         0.09
Other                            0.03
White                            0.77
Name: ethnicity, dtype: float64
genhealth
----------------------
Excellent   0.21
Fair        0.11
Good        0.29
Poor        0.04
Very good   0.36
Name: genhealth, dtype: float64
diabetic
----------------------
No                        0.84
No, borderline diabetes   0.02
Yes                       0.13
Yes (during pregnancy)    0.01
Name: diabetic, dtype: float64

1. 我们还应该查看一些数值特征的描述性统计。对于不良的身体健康和心理健康天数的中位数都是 0；也就是说，至少一半的观察结果报告没有不良的身体健康天数，至少一半报告没有不良的心理健康天数：

   healthinfo[num_cols].\
     agg(['count','min','median','max']).T
                          count    min    median  max
   bmi                    30,000   12     27      92
   physicalhealthbaddays  30,000   0      0       30
   mentalhealthbaddays    30,000   0      0       30
   sleeptimenightly       30,000   1      7       24
   

我们需要进行一些缩放。我们还需要对分类特征进行编码。数值特征也有一些极端值。sleeptimenightly的值为 24 似乎不太可能！处理它们可能是个好主意。

1. 现在，我们已经准备好构建我们的流水线。让我们创建训练和测试的 DataFrame：

   X_train, X_test, y_train, y_test =  \
     train_test_split(healthinfo[num_cols + 
       binary_cols + cat_cols + spec_cols1 +
       spec_cols2 + spec_cols3],\
     healthinfo[['heartdisease']], test_size=0.2,
       random_state=0)
   

2. 接下来，我们将设置列转换。我们将创建一个 one-hot 编码器实例，我们将使用它来处理所有分类特征。对于数值列，我们将使用OutlierTrans对象去除极端值，然后填充中位数。

我们将使用MakeOrdinal转换器将agecategory特征转换为有序特征，并使用ReplaceVals转换器对genhealth和diabetic特征进行编码。

我们将在下一步将列转换添加到我们的流水线中：

ohe = OneHotEncoder(drop='first', sparse=False)
standtrans = make_pipeline(OutlierTrans(3),
  SimpleImputer(strategy="median"),
  StandardScaler())
spectrans1 = make_pipeline(MakeOrdinal(),
  StandardScaler())
spectrans2 = make_pipeline(ReplaceVals(rep_dict),
  StandardScaler())
spectrans3 = make_pipeline(ReplaceVals(rep_dict))
bintrans = make_pipeline(ohe)
cattrans = make_pipeline(ohe)
coltrans = ColumnTransformer(
  transformers=[
    ("stand", standtrans, num_cols),
    ("spec1", spectrans1, spec_cols1),
    ("spec2", spectrans2, spec_cols2),
    ("spec3", spectrans3, spec_cols3),
    ("bin", bintrans, binary_cols),
    ("cat", cattrans, cat_cols),
  ]
)

1. 现在，我们已经准备好设置和调整我们的流水线。首先，我们将实例化逻辑回归和分层 k 折对象，我们将使用递归特征消除。回想一下，递归特征消除需要一个估计器。我们使用分层 k 折来确保每个折叠中目标值的分布大致相同。

现在，我们必须为我们的模型创建另一个逻辑回归实例。我们将class_weight参数设置为balanced。这应该会提高模型处理类别不平衡的能力。然后，我们将列转换、递归特征消除和逻辑回归实例添加到我们的流水线中，然后对其进行调整：

lrsel = LogisticRegression(random_state=1, 
  max_iter=1000)
kf = StratifiedKFold(n_splits=5, shuffle=True)
rfecv = RFECV(estimator=lrsel, cv=kf)
lr = LogisticRegression(random_state=1,
  class_weight='balanced', max_iter=1000)
pipe1 = make_pipeline(coltrans, rfecv, lr)
pipe1.fit(X_train, y_train.values.ravel())

1. 在调整后，我们需要做一些工作来恢复流水线中的列名。我们可以使用bin转换器的 one-hot 编码器的get_feature_names方法和cat转换器的get_feature_names方法。这为我们提供了编码后的二元和分类特征的列名。数值特征的名称保持不变。我们将在后面使用这些特征名：

   new_binary_cols = \
     pipe1.named_steps['columntransformer'].\
     named_transformers_['bin'].\
     named_steps['onehotencoder'].\
     get_feature_names(binary_cols)
   new_cat_cols = \
     pipe1.named_steps['columntransformer'].\
     named_transformers_['cat'].\
     named_steps['onehotencoder'].\
     get_feature_names(cat_cols)
   new_cols = np.concatenate((np.array(num_cols +
     spec_cols1 + spec_cols2 + spec_cols3),
     new_binary_cols, new_cat_cols))
   new_cols
   array(['bmi', 'physicalhealthbaddays',
          'mentalhealthbaddays', 'sleeptimenightly',
          'agecategory', 'genhealth', 'diabetic',
          'smoking_Yes', 'alcoholdrinkingheavy_Yes',
          'stroke_Yes', 'walkingdifficult_Yes',
          'physicalactivity_Yes', 'asthma_Yes',
          'kidneydisease_Yes', 'skincancer_Yes',
          'gender_Male', 'ethnicity_Asian',
          'ethnicity_Black', 'ethnicity_Hispanic',
          'ethnicity_Other', 'ethnicity_White'],
         dtype=object)
   

2. 现在，让我们看看递归特征消除的结果。我们可以使用rfecv对象的ranking_属性来获取每个特征的排名。那些排名为1的特征将被选入我们的模型。

如果我们使用rfecv对象的get_support方法或support_属性代替ranking_属性，我们只会得到那些将在我们的模型中使用的特点——即那些排名为 1 的特点。我们将在下一步做这个：

rankinglabs = \
 np.column_stack((pipe1.named_steps['rfecv'].ranking_,
 new_cols))
pd.DataFrame(rankinglabs,
 columns=['rank','feature']).\
 sort_values(['rank','feature']).\
 set_index("rank")
                       feature
rank                          
1                  agecategory
1     alcoholdrinkingheavy_Yes
1                   asthma_Yes
1                     diabetic
1              ethnicity_Asian
1              ethnicity_Other
1              ethnicity_White
1                  gender_Male
1                    genhealth
1            kidneydisease_Yes
1                  smoking_Yes
1                   stroke_Yes
1         walkingdifficult_Yes
2           ethnicity_Hispanic
3               skincancer_Yes
4                          bmi
5        physicalhealthbaddays
6             sleeptimenightly
7          mentalhealthbaddays
8         physicalactivity_Yes
9              ethnicity_Black

1. 我们可以从逻辑回归的系数中获取优势比。回想一下，优势比是指数化的系数。有 13 个系数，这是有意义的，因为我们之前学习到有 13 个特征得到了 1 的排名。

我们将使用rfecv步骤的get_support方法来获取所选特征的名称，并创建一个包含这些名称和优势比（oddswithlabs）的 NumPy 数组。然后我们创建一个 pandas DataFrame，并按优势比降序排序。

毫不奇怪，那些曾经中风的人和老年人患心脏病的可能性要大得多。如果个人曾经中风，他们在其他条件相同的情况下患心脏病的几率是三倍。另一方面，随着年龄类别的增加，患心脏病的几率增加 2.88 倍。另一方面，随着总体健康状况的提高，患心脏病的几率大约减少一半（57%）；比如说，从“一般”到“良好”。令人惊讶的是，在控制其他条件的情况下，大量饮酒与心脏病几率降低有关：

oddsratios = np.exp(pipe1.\
  named_steps['logisticregression'].coef_)
oddsratios.shape
(1, 13)
selcols = new_cols[pipe1.\
  named_steps['rfecv'].get_support()]
oddswithlabs = np.column_stack((oddsratios.\
  ravel(), selcols))
pd.DataFrame(oddswithlabs, 
  columns=['odds','feature']).\
  sort_values(['odds'], ascending=False).\
  set_index('odds')
                        feature
odds                          
3.01                stroke_Yes
2.88               agecategory
2.12               gender_Male
1.97         kidneydisease_Yes
1.75                  diabetic
1.55               smoking_Yes
1.52                asthma_Yes
1.30      walkingdifficult_Yes
1.27           ethnicity_Other
1.22           ethnicity_White
0.72           ethnicity_Asian
0.61  alcoholdrinkingheavy_Yes
0.57                 genhealth

现在我们已经拟合了逻辑回归模型，我们准备对其进行评估。在下一节中，我们将花时间探讨各种性能指标，包括准确率和灵敏度。我们将使用我们在第六章，“准备模型评估”中介绍的一些概念。

评估逻辑回归模型

一个分类模型性能的最直观的衡量标准是其准确率——也就是说，我们的预测有多正确。然而，在某些情况下，我们可能至少和准确率一样关心灵敏度——即我们正确预测的阳性案例的百分比；我们甚至可能愿意牺牲一点准确率来提高灵敏度。疾病预测模型通常属于这一类。但是，每当存在类别不平衡时，准确率和灵敏度等指标可能会给我们提供关于模型性能的非常不同的估计。

除了关注准确率或灵敏度之外，我们还可能担心我们的模型的特异性或精确度。我们可能希望有一个模型能够以高可靠性识别出负面案例，即使这意味着它不能很好地识别正面案例。特异性是模型识别出的所有负面案例所占的百分比。

精确度，即预测为正面的案例中实际为正面的比例，是另一个重要的衡量指标。对于某些应用来说，限制误报非常重要，即使这意味着我们必须容忍较低的灵敏度。一个使用图像识别来识别坏苹果的苹果种植者可能更倾向于选择一个高精确度的模型，而不是一个更灵敏的模型，不希望不必要地丢弃苹果。

通过查看混淆矩阵可以使这一点更加清晰：

图 10.2 – 二元目标按预测值预测的实际值混淆矩阵

混淆矩阵帮助我们理解准确率、灵敏度、特异性和精确度。准确率是我们预测正确的观察值的百分比。这可以更精确地表述如下：

敏感性是指我们正确预测正值的次数除以正值的总数。再次查看混淆矩阵可能会有所帮助，以确认实际正值可以是预测正值（TP）或预测负值（FN）。敏感性也被称为召回率或真正阳性率：

特异性是指我们正确预测负值的次数（TN）除以实际负值总数（TN + FP）。特异性也被称为真正阴性率：

精确度是指我们正确预测正值的次数（TP）除以预测的正值总数：

我们在第六章中更详细地介绍了这些概念，准备模型评估。在本节中，我们将检查心脏病逻辑回归模型的准确性、敏感性、特异性和精确度：

1. 我们可以使用上一节中拟合的管道的predict方法来从我们的逻辑回归中生成预测。然后，我们可以生成一个混淆矩阵：

   pred = pipe1.predict(X_test)
   cm = skmet.confusion_matrix(y_test, pred)
   cmplot = skmet.ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=['Negative', 'Positive'])
   cmplot.plot()
   cmplot.ax_.set(title='Heart Disease Prediction Confusion Matrix', 
     xlabel='Predicted Value', ylabel='Actual Value')
   

这产生了以下图表：

图 10.3 – 心脏病预测混淆矩阵

这里首先要注意的是，大部分动作都在左上角，我们在测试数据中正确预测了实际负值。这将大大有助于我们的准确性。尽管如此，我们还是有相当数量的假阳性。当没有心脏病时，我们预测心脏病 1,430 次（在 5,506 个负实例中）。我们似乎在识别正性心脏病实例方面做得还不错，正确分类了 392 个实例（在 494 个正实例中）。

1. 让我们计算准确性、敏感性、特异性和精确度。总体准确性并不高，为 74%。敏感性相当不错，为 79%。（当然，敏感性的好坏取决于领域和判断。对于像心脏病这样的领域，我们可能希望它更高。）这可以在以下代码中看到：

   tn, fp, fn, tp = skmet.confusion_matrix(y_test.values.ravel(), pred).ravel()
   tn, fp, fn, tp
   (4076, 1430, 102, 392)
   accuracy = (tp + tn) / pred.shape[0]
   accuracy
   0.7446666666666667
   sensitivity = tp / (tp + fn)
   sensitivity
   0.7935222672064778
   specificity = tn / (tn+fp)
   specificity
   0.7402833272793317
   precision = tp / (tp + fp)
   precision
   0.21514818880351264
   

2. 我们可以使用metrics模块以更直接的方式来进行这些计算（我在上一步中选择了更迂回的方法来展示计算过程）：

   print("accuracy: %.2f, sensitivity: %.2f, specificity: %.2f, precision: %.2f"  %
     (skmet.accuracy_score(y_test.values.ravel(), pred),
     skmet.recall_score(y_test.values.ravel(), pred),
     skmet.recall_score(y_test.values.ravel(), pred,
     pos_label=0),
     skmet.precision_score(y_test.values.ravel(), pred)))
   accuracy: 0.74, sensitivity: 0.79, specificity: 0.74, precision: 0.22
   

我们模型的最大问题是精确度非常低——即 22%。这是由于大量的假阳性。我们的模型预测正性的大多数情况下都是错误的。

除了我们已经计算出的四个指标之外，获取假阳性率也可能很有帮助。假阳性率是指我们的模型在实际值为负时预测为正的倾向：

1. 让我们计算假阳性率：

   falsepositiverate = fp / (tn + fp)
   falsepositiverate
   0.25971667272066834
   

因此，26%的时间，当一个人没有心脏病时，我们预测他有心脏病。虽然我们当然希望限制假阳性的数量，但这通常意味着牺牲一些敏感性。我们将在本节后面演示为什么这是真的。

1. 我们应该仔细查看模型生成的预测概率。在这里，正类预测的阈值是 0.5，这在逻辑回归中通常是默认值。（回想一下，逻辑回归预测的是类成员的概率。我们需要一个伴随的决策规则，比如 0.5 阈值，来预测类别。）这可以在以下代码中看到：

   pred_probs = pipe1.predict_proba(X_test)[:, 1]
   probdf = \
     pd.DataFrame(zip(pred_probs, pred,
     y_test.values.ravel()),
     columns=(['prob','pred','actual']))
   probdf.groupby(['pred'])['prob'].\
     agg(['min','max','count'])
           min        max        count
   pred                 
   0       0.01       0.50       4178
   1       0.50       0.99       1822
   

2. 我们可以使用核密度估计（KDE）图来可视化这些概率。我们还可以看到不同的决策规则可能如何影响我们的预测。例如，我们可以将阈值从 0.5 移动到 0.25。乍一看，这有一些优点。两个可能阈值之间的区域比没有心脏病病例的心脏病病例要多一些。我们将得到虚线之间的棕色区域，正确预测心脏病，而不会在 0.5 阈值下这样做。这比线之间的绿色区域更大，在 0.5 阈值下，我们将一些真正的阴性预测变成了 0.25 阈值下的假阳性：

   sns.kdeplot(probdf.loc[probdf.actual==1].prob,
     shade=True, color='red',label="Heart Disease")
   sns.kdeplot(probdf.loc[probdf.actual==0].prob,
     shade=True,color='green',label="No Heart Disease")
   plt.axvline(0.25, color='black', linestyle='dashed',
     linewidth=1)
   plt.axvline(0.5, color='black', linestyle='dashed',
     linewidth=1)
   plt.title("Predicted Probability Distribution")
   plt.legend(loc="upper left")
   

这会产生以下图表：

图 10.4 – 预测的心脏病概率分布

让我们比之前更仔细地考虑精确度和敏感性之间的权衡。记住，精确度是我们预测正类值时正确率的比率。敏感性，也称为召回率或真正阳性率，是我们识别实际正例为正例的比率。

1. 我们可以如下绘制精确度和敏感性曲线：

   prec, sens, ths = skmet.precision_recall_curve(y_test, pred_probs)
   sens = sens[1:-20]
   prec = prec[1:-20]
   ths  = ths[:-20]
   fig, ax = plt.subplots()
   ax.plot(ths, prec, label='Precision')
   ax.plot(ths, sens, label='Sensitivity')
   ax.set_title('Precision and Sensitivity by Threshold')
   ax.set_xlabel('Threshold')
   ax.set_ylabel('Precision and Sensitivity')
   ax.legend()
   

这会产生以下图表：

图 10.5 – 阈值值下的精确度和敏感性

当阈值超过 0.2 时，敏感性的下降比精确度的增加更为明显。

1. 通常，查看假阳性率与敏感性率也是很有帮助的。假阳性率是我们模型在实际值为负时预测为正的倾向。了解这种关系的一种方法是通过 ROC 曲线：

   fpr, tpr, ths = skmet.roc_curve(y_test, pred_probs)
   ths = ths[1:]
   fpr = fpr[1:]
   tpr = tpr[1:]
   fig, ax = plt.subplots()
   ax.plot(fpr, tpr, linewidth=4, color="black")
   ax.set_title('ROC curve')
   ax.set_xlabel('False Positive Rate')
   ax.set_ylabel('Sensitivity')
   

这会产生以下图表：

图 10.6 – ROC 曲线

在这里，我们可以看到，随着假阳性率的增加，我们获得的敏感性增加越来越少。超过 0.5 的假阳性率，几乎没有回报。

1. 可能还有助于仅通过阈值绘制假阳性率和敏感性：

   fig, ax = plt.subplots()
   ax.plot(ths, fpr, label="False Positive Rate")
   ax.plot(ths, tpr, label="Sensitivity")
   ax.set_title('False Positive Rate and Sensitivity by Threshold')
   ax.set_xlabel('Threshold')
   ax.set_ylabel('False Positive Rate and Sensitivity')
   ax.legend()
   

这会产生以下图表：

图 10.7 – 灵敏度和假阳性率

这里，我们可以看到，当我们把阈值降低到 0.25 以下时，假阳性率比灵敏度增加得更快。

这最后两个可视化暗示了找到最佳阈值值——即在灵敏度和假阳性率之间有最佳权衡的值；至少在数学上，忽略领域知识。

1. 我们将计算argmax函数。我们想要这个索引处的阈值值。根据这个计算，最佳阈值是 0.46，这与默认值并不太不同：

   jthresh = ths[np.argmax(tpr – fpr)]
   jthresh
   0.45946882675453804
   

2. 我们可以根据这个替代阈值重新做混淆矩阵：

   pred2 = np.where(pred_probs>=jthresh,1,0)
   cm = skmet.confusion_matrix(y_test, pred2)
   cmplot = skmet.ConfusionMatrixDisplay(
     confusion_matrix=cm, 
     display_labels=['Negative', 'Positive'])
   cmplot.plot()
   cmplot.ax_.set(
     title='Heart Disease Prediction Confusion Matrix', 
     xlabel='Predicted Value', ylabel='Actual Value')
   

这产生了以下图表：

图 10.8 – 心脏病预测混淆矩阵

1. 这给我们带来了灵敏度的微小提升：

   skmet.recall_score(y_test.values.ravel(), pred)
   0.7935222672064778
   skmet.recall_score(y_test.values.ravel(), pred2)
   0.8380566801619433
   

这里要说明的不是我们应该随意更改阈值。这通常是一个坏主意。但我们应该记住两点。首先，当我们有一个高度不平衡的类别时，0.5 的阈值可能没有意义。其次，这是依赖领域知识的一个重要地方。对于某些分类问题，假阳性远不如假阴性重要。

在本节中，我们关注了灵敏度、精确度和假阳性率作为模型性能的度量。这部分是因为空间限制，也因为这个特定目标的问题——不平衡的类别和可能对灵敏度的偏好。在接下来的几章中，我们将强调其他度量，如准确性和特异性，在其他我们将构建的模型中。在本章的其余部分，我们将探讨逻辑回归的几个扩展，包括正则化和多项式逻辑回归。

使用逻辑回归进行正则化

如果你已经阅读了第七章《线性回归模型》并阅读了本章的第一节，你对正则化的工作原理已经有了很好的了解。我们向估计器添加一个惩罚，以最小化我们的参数估计。这个惩罚的大小通常基于模型性能的度量来调整。我们将在本节中探讨这一点。按照以下步骤进行：

1. 我们将加载与上一节中使用的相同模块，以及我们需要进行必要超参数调整的模块。我们将使用RandomizedSearchCV和uniform来找到我们惩罚强度的最佳值：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.preprocessing import OneHotEncoder
   from sklearn.pipeline import make_pipeline
   from sklearn.impute import SimpleImputer
   from sklearn.compose import ColumnTransformer
   from sklearn.model_selection import RepeatedStratifiedKFold
   from sklearn.linear_model import LogisticRegression
   from sklearn.model_selection import RandomizedSearchCV
   from scipy.stats import uniform
   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans,\
     MakeOrdinal, ReplaceVals
   

2. 接下来，我们将加载心脏病数据并进行一些处理：

   healthinfo = pd.read_csv("data/healthinfosample.csv")
   healthinfo.set_index("personid", inplace=True)
   healthinfo['heartdisease'] = \
     np.where(healthinfo.heartdisease=='No',0,1).\
     astype('int')
   

3. 接下来，我们将组织我们的特征，以便于我们在接下来的几个步骤中进行的列转换：

   num_cols = ['bmi','physicalhealthbaddays',
      'mentalhealthbaddays','sleeptimenightly']
   binary_cols = ['smoking','alcoholdrinkingheavy',
     'stroke','walkingdifficult','physicalactivity',
     'asthma','kidneydisease','skincancer']
   cat_cols = ['gender','ethnicity']
   spec_cols1 = ['agecategory']
   spec_cols2 = ['genhealth']
   spec_cols3 = ['diabetic']
   rep_dict = {
     'genhealth': {'Poor':0,'Fair':1,'Good':2,
       'Very good':3,'Excellent':4},
     'diabetic': {'No':0,
       'No, borderline diabetes':0,'Yes':1,
       'Yes (during pregnancy)':1}           
   }
   

4. 现在，我们必须创建测试和训练数据框：

   X_train, X_test, y_train, y_test =  \
     train_test_split(healthinfo[num_cols + 
       binary_cols + cat_cols + spec_cols1 +
       spec_cols2 + spec_cols3],\
     healthinfo[['heartdisease']], test_size=0.2,
       random_state=0)
   

5. 然后，我们必须设置列转换：

   ohe = OneHotEncoder(drop='first', sparse=False)
   standtrans = make_pipeline(OutlierTrans(3),
     SimpleImputer(strategy="median"),
     StandardScaler())
   spectrans1 = make_pipeline(MakeOrdinal(),
     StandardScaler())
   spectrans2 = make_pipeline(ReplaceVals(rep_dict),
     StandardScaler())
   spectrans3 = make_pipeline(ReplaceVals(rep_dict))
   bintrans = make_pipeline(ohe)
   cattrans = make_pipeline(ohe)
   coltrans = ColumnTransformer(
     transformers=[
       ("stand", standtrans, num_cols),
       ("spec1", spectrans1, spec_cols1),
       ("spec2", spectrans2, spec_cols2),
       ("spec3", spectrans3, spec_cols3),
       ("bin", bintrans, binary_cols),
       ("cat", cattrans, cat_cols),
     ]
   )
   

6. 现在，我们已经准备好运行我们的模型了。我们将实例化逻辑回归和重复分层 k 折对象。然后，我们将创建一个包含之前步骤中的列转换和逻辑回归的管道。

之后，我们将为超参数创建一个字典列表，而不是像在这本书中之前所做的那样只创建一个字典。这是因为并非所有超参数都能一起工作。例如，我们不能使用newton-cg求解器与 L1 惩罚一起使用。字典键名前缀的logisticregression__（注意双下划线）表示我们希望将这些值传递到管道的逻辑回归步骤。

我们将设置随机网格搜索的n_iter参数为20，以便它采样超参数 20 次。每次，网格搜索都将从列出的一个字典中选择超参数。我们将指示我们希望网格搜索的评分基于 ROC 曲线下的面积：

lr = LogisticRegression(random_state=1, class_weight='balanced', max_iter=1000)
kf = RepeatedStratifiedKFold(n_splits=7, n_repeats=3, random_state=0)
pipe1 = make_pipeline(coltrans, lr)
reg_params = [
  {
    'logisticregression__solver': ['liblinear'],
    'logisticregression__penalty': ['l1','l2'],
    'logisticregression__C': uniform(loc=0, scale=10)
  },
  {
    'logisticregression__solver': ['newton-cg'],
    'logisticregression__penalty': ['l2'],
    'logisticregression__C': uniform(loc=0, scale=10)
  },
  {
    'logisticregression__solver': ['saga'],
    'logisticregression__penalty': ['elasticnet'],
    'logisticregression__l1_ratio': uniform(loc=0, scale=1),   
    'logisticregression__C': uniform(loc=0, scale=10)
  }
]
rs = RandomizedSearchCV(pipe1, reg_params, cv=kf, 
  n_iter=20, scoring='roc_auc')
rs.fit(X_train, y_train.values.ravel())

1. 在拟合搜索后，best_params属性给出了与最高分数相关的参数。弹性网络回归，其 L1 比率更接近 L1 而不是 L2，表现最佳：

   rs.best_params_
   {'logisticregression__C': 0.6918282397356423,
    'logisticregression__l1_ratio': 0.758705704020254,
    'logisticregression__penalty': 'elasticnet',
    'logisticregression__solver': 'saga'}
   rs.best_score_
   0.8410275986723489
   

2. 让我们看看网格搜索的一些其他高分。最好的三个模型得分几乎相同。一个使用弹性网络回归，另一个使用 L1，另一个使用 L2。

网格搜索的cv_results_字典为我们提供了关于尝试过的 20 个模型的大量信息。该字典中的params列表结构有些复杂，因为某些键在某些迭代中不存在，例如L1_ratio。我们可以使用json_normalize来简化结构：

results = \
  pd.DataFrame(rs.cv_results_['mean_test_score'], \
    columns=['meanscore']).\
  join(pd.json_normalize(rs.cv_results_['params'])).\
  sort_values(['meanscore'], ascending=False)
results.head(3).T
                              15          4      12
meanscore                     0.841       0.841  0.841
logisticregression__C         0.692       1.235  0.914
logisticregression__l1_ratio  0.759       NaN    NaN
logisticregression__penalty   elasticnet  l1     l2
logisticregression__solver  saga  liblinear  liblinear

1. 让我们看看混淆矩阵：

   pred = rs.predict(X_test)
   cm = skmet.confusion_matrix(y_test, pred)
   cmplot = \
     skmet.ConfusionMatrixDisplay(confusion_matrix=cm,
     display_labels=['Negative', 'Positive'])
   cmplot.plot()
   cmplot.ax_.\
     set(title='Heart Disease Prediction Confusion Matrix', 
     xlabel='Predicted Value', ylabel='Actual Value')
   

这生成了以下图表：

图 10.9 – 心脏病预测混淆矩阵

1. 让我们也看看一些度量指标。我们的分数与没有正则化的模型基本没有变化：

   print("accuracy: %.2f, sensitivity: %.2f, specificity: %.2f, precision: %.2f"  %
     (skmet.accuracy_score(y_test.values.ravel(), pred),
     skmet.recall_score(y_test.values.ravel(), pred),
     skmet.recall_score(y_test.values.ravel(), pred,
       pos_label=0),
     skmet.precision_score(y_test.values.ravel(), pred)))
   accuracy: 0.74, sensitivity: 0.79, specificity: 0.74, precision: 0.21
   

尽管正则化并没有明显提高我们模型的表现，但很多时候它确实做到了。在使用 L1 正则化时，也不必过于担心特征选择，因为不太重要的特征的权重将会是 0。

尽管我们还没有解决如何处理目标值超过两个的可能值的模型，尽管上一两节的几乎所有讨论也适用于多类模型。在下一节中，我们将学习如何使用多项式逻辑回归来建模多类目标。

多项式逻辑回归

如果逻辑回归只适用于二元分类问题，那么它就不会那么有用。幸运的是，当我们的目标值超过两个时，我们可以使用多项式逻辑回归。

在本节中，我们将处理关于空气和工艺温度、扭矩和旋转速度作为机器故障函数的数据。

注意

这个关于机器故障的数据集可在www.kaggle.com/datasets/shivamb/machine-predictive-maintenance-classification公开使用。有 10,000 个观测值，12 个特征，以及两个可能的目标。一个是二元的——也就是说，机器故障或未故障。另一个是故障类型。这个数据集中的实例是合成的，由一个旨在模仿机器故障率和原因的过程生成。

让我们学习如何使用多项式逻辑回归来建模机器故障：

1. 首先，我们将导入现在熟悉的库。我们还将导入cross_validate，我们首次在第六章准备模型评估中使用它，以帮助我们评估我们的模型：

   import pandas as pd
   import numpy as np
   from sklearn.model_selection import train_test_split
   from sklearn.preprocessing import StandardScaler
   from sklearn.preprocessing import OneHotEncoder
   from sklearn.pipeline import make_pipeline
   from sklearn.impute import SimpleImputer
   from sklearn.compose import ColumnTransformer
   from sklearn.model_selection import RepeatedStratifiedKFold
   from sklearn.linear_model import LogisticRegression
   from sklearn.model_selection import cross_validate
   import os
   import sys
   sys.path.append(os.getcwd() + "/helperfunctions")
   from preprocfunc import OutlierTrans
   

2. 我们将加载机器故障数据并查看其结构。我们没有缺失数据。这是个好消息：

   machinefailuretype = pd.read_csv("data/machinefailuretype.csv")
   machinefailuretype.info()
   <class 'pandas.core.frame.DataFrame'>
   RangeIndex: 10000 entries, 0 to 9999
   Data columns (total 10 columns):
    #   Column              Non-Null Count      Dtype
   ---  ------             --------------       ----  
    0   udi                 10000 non-null      int64
    1   product             10000 non-null      object 
    2   machinetype         10000 non-null      object 
    3   airtemp             10000 non-null      float64
    4   processtemperature  10000 non-null      float64
    5   rotationalspeed     10000 non-null      int64
    6   torque              10000 non-null      float64
    7   toolwear            10000 non-null      int64
    8   fail                10000 non-null      int64
    9   failtype            10000 non-null      object 
   dtypes: float64(3), int64(4), object(3)
   memory usage: 781.4+ KB
   

3. 让我们查看几行。machinetype的值为L、M和H。这些值分别是低、中、高质机器的代理：

   machinefailuretype.head()
      udi product machinetype airtemp processtemperature\
   0  1   M14860  M           298     309 
   1  2   L47181  L           298     309 
   2  3   L47182  L           298     308 
   3  4   L47183  L           298     309 
   4  5   L47184  L           298     309 
      Rotationalspeed  torque  toolwear  fail  failtype  
   0   1551             43       0        0    No Failure
   1   1408             46       3        0    No Failure
   2   1498             49       5        0    No Failure
   3   1433             40       7        0    No Failure
   4   1408             40       9        0    No Failure
   

4. 我们还应该生成一些频率：

   machinefailuretype.failtype.value_counts(dropna=False).sort_index()
   Heat Dissipation Failure    112
   No Failure                  9652
   Overstrain Failure          78
   Power Failure               95
   Random Failures             18
   Tool Wear Failure           45
   Name: failtype, dtype: int64
   machinefailuretype.machinetype.\
     value_counts(dropna=False).sort_index()
   H    1003
   L    6000
   M    2997
   Name: machinetype, dtype: int64
   

5. 让我们将failtype值合并，并为它们创建数字代码。由于随机故障的计数很低，我们将随机故障和工具磨损故障合并：

   def setcode(typetext):
     if (typetext=="No Failure"):
       typecode = 1
     elif (typetext=="Heat Dissipation Failure"):
       typecode = 2
     elif (typetext=="Power Failure"):
       typecode = 3
     elif (typetext=="Overstrain Failure"):
       typecode = 4
     else:
       typecode = 5
     return typecode
   machinefailuretype["failtypecode"] = \
     machinefailuretype.apply(lambda x: setcode(x.failtype), axis=1)
   

6. 我们应该确认failtypecode是否按我们的意图工作：

   machinefailuretype.groupby(['failtypecode','failtype']).size().\
     reset_index()
     failtypecode   failtype                     0
   0     1          No Failure                   9652
   1     2          Heat Dissipation Failure     112
   2     3          Power Failure                95
   3     4          Overstrain Failure           78
   4     5          Random Failures              18
   5     5          Tool Wear Failure            45
   

7. 让我们也获取一些描述性统计信息：

   num_cols = ['airtemp','processtemperature','rotationalspeed',
     'torque','toolwear']
   cat_cols = ['machinetype']
   machinefailuretype[num_cols].agg(['min','median','max']).T
                         min      median    max
   airtemp               295      300       304
   processtemperature    306      310       314
   rotationalspeed       1,168    1,503     2,886
   torque                4        40        77
   toolwear              0        108       253
   

8. 现在，让我们创建测试和训练数据框。我们还将设置列转换：

   X_train, X_test, y_train, y_test =  \
     train_test_split(machinefailuretype[num_cols +
     cat_cols], machinefailuretype[['failtypecode']],
     test_size=0.2, random_state=0)
   ohe = OneHotEncoder(drop='first', sparse=False)
   standtrans = make_pipeline(OutlierTrans(3),
     SimpleImputer(strategy="median"),
     StandardScaler())
   cattrans = make_pipeline(ohe)
   coltrans = ColumnTransformer(
     transformers=[
       ("stand", standtrans, num_cols),
       ("cat", cattrans, cat_cols),
     ]
   )
   

9. 现在，让我们设置一个包含我们的列转换和多项逻辑回归模型的管道。当我们实例化逻辑回归时，只需将multi_class属性设置为多项式即可：

   lr = LogisticRegression(random_state=0, 
     multi_class='multinomial', solver='lbfgs',
     max_iter=1000)
   kf = RepeatedStratifiedKFold(n_splits=10,
     n_repeats=5, random_state=0)
   pipe1 = make_pipeline(coltrans, lr)
   

10. 现在，我们可以生成一个混淆矩阵：

    cm = skmet.confusion_matrix(y_test, 
       pipe1.fit(X_train, y_train.values.ravel()).\
       predict(X_test))
    cmplot = \
       skmet.ConfusionMatrixDisplay(confusion_matrix=cm,
       display_labels=['None', 'Heat','Power','Overstrain','Other'])
    cmplot.plot()
    cmplot.ax_.\
      set(title='Machine Failure Type Confusion Matrix', 
      xlabel='Predicted Value', ylabel='Actual Value')
    

这会产生以下图表：

图 10.10 – 预测机器故障类型的混淆矩阵

混淆矩阵显示，当存在故障时，我们的模型在预测故障类型方面做得并不好，尤其是电力故障或其他故障。

1. 我们可以使用cross_validate来评估这个模型。我们主要得到准确率、精确率和敏感度（召回率）的优异成绩。然而，这是误导性的。当类别如此不平衡（几乎所有实例都是无故障）时，加权分数会受到包含几乎所有值的类别的严重影响。我们的模型能够可靠地正确预测无故障。

如果我们查看f1_macro分数（回忆一下第六章，准备模型评估，其中f1是精确率和敏感性的调和平均值），我们会看到，除了无故障类别之外，我们的模型在其它类别上表现并不好。（macro分数只是一个简单的平均值。）

我们本可以使用分类报告在这里，就像我们在第六章，“准备模型评估”中做的那样，但我有时发现生成我需要的统计数据很有帮助：

scores = cross_validate(
  pipe1, X_train, y_train.values.ravel(), \
  scoring=['accuracy', 'precision_weighted',
           'recall_weighted', 'f1_macro',
           'f1_weighted'], 
  cv=kf, n_jobs=-1)
accuracy, precision, sensitivity, f1_macro, f1_weighted = \
  np.mean(scores['test_accuracy']),\
  np.mean(scores['test_precision_weighted']),\
  np.mean(scores['test_recall_weighted']),\
  np.mean(scores['test_f1_macro']),\
  np.mean(scores['test_f1_weighted'])
accuracy, precision, sensitivity, f1_macro, f1_weighted
(0.9716499999999999,
 0.9541025493784612,
 0.9716499999999999,
 0.3820938909478524,
 0.9611411229222823)

在本节中，我们探讨了如何构建多项逻辑回归模型。这种方法无论目标变量是名义变量还是有序变量都适用。在这种情况下，它是名义变量。我们还看到了如何将用于具有二元目标的逻辑回归模型中的模型评估方法进行扩展。我们回顾了当我们有超过两个类别时如何解释混淆矩阵和评分指标。

摘要

逻辑回归多年来一直是我预测分类目标的首选工具。它是一个高效且偏差低的算法。一些其缺点，如高方差和难以处理高度相关的预测因子，可以通过正则化和特征选择来解决。我们在本章中探讨了如何做到这一点。我们还考察了如何处理不平衡类别，以及这些目标对建模和结果解释的含义。

在下一章中，我们将探讨分类中一个非常流行的逻辑回归替代方案——决策树。我们会看到决策树具有许多优点，使得它们在需要建模复杂性时成为一个特别好的选择，而不必像使用逻辑回归那样过多地担心我们的特征是如何指定的。