机器学习的数据清理和探索-二-机器学习的数据清理和探索（二）第二部分 – 预处理、特征选择和采样任何进行过目标对数变

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

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第二部分 – 预处理、特征选择和采样

任何进行过目标对数变换或特征缩放的人都会深刻体会到分析预处理的重要性。请举手，如果你曾经自信地认为你的模型近似于真理，但后来尝试了一个相当明显的变换，并意识到你的原始模型与真理相差有多远。编码、变换和缩放数据并非花招，尽管有时人们会有这种印象。我们应用这种预处理是因为 1)它使我们更接近捕捉现实世界的过程，2)因为许多机器学习算法在缩放数据上表现更好。

特征选择同样重要。一句好的格言是：永远不要用 N 个特征来构建模型，当 N - 1 个特征就能达到同样的效果。值得记住的是，这比拥有太多特征要复杂得多。有时拥有 3 个特征就太多，而有时拥有 103 个特征则恰到好处。问题实际上在于特征是否过于相关，以至于它们对目标的影响难以分离。当这种情况不成立时，过拟合和不稳定结果的风险会大幅增加。我们在本书的这一部分和随后的部分都特别注意特征选择。

本部分关于模型评估的章节为我们准备在本书其余部分将要进行的工作。我们详细探讨了回归和分类模型的评估，分别针对连续和分类目标。我们还学习了如何构建管道和进行交叉验证。最重要的是，我们学习了数据泄露及其如何避免。

本节包括以下章节：

第四章，特征编码、变换和缩放
第五章，特征选择
第六章，为模型评估做准备

第四章：第四章：编码、转换和缩放特征

本书的前三章重点介绍了数据清洗、探索以及如何识别缺失值和异常值。接下来的几章将深入探讨特征工程，本章将从编码、转换和缩放数据以提高机器学习模型性能的技术开始。

通常，机器学习算法需要以某种形式对变量进行编码。此外，我们的模型在缩放后通常表现更好，这样具有更高变异性的特征就不会压倒优化过程。我们将向您展示如何在使用特征范围差异很大的情况下使用不同的缩放技术。

具体来说，在本章中，我们将探讨以下主要主题：

创建训练数据集并避免数据泄露
识别要删除的不相关或冗余观察结果
编码分类特征
使用中等或高基数编码特征
转换特征
分箱特征
特征缩放

技术要求

在本章中，我们将与feature-engine和category_encoders包以及sklearn库进行大量工作。您可以使用pip安装这些包，命令为pip install feature-engine，pip install category_encoders，以及pip install scikit-learn。本章中的代码使用了sklearn的 0.24.2 版本，feature-engine的 1.1.2 版本，以及category_encoders的 2.2.2 版本。请注意，无论是pip install feature-engine还是pip install feature_engine都可以工作。

本章的所有代码都可以在 GitHub 上找到，链接为github.com/PacktPublishing/Data-Cleaning-and-Exploration-with-Machine-Learning/tree/main/4.%20PruningEncodingandRescalingFeatures。

创建训练数据集并避免数据泄露

我们模型性能的最大威胁之一是数据泄露。数据泄露发生在我们的模型被训练数据集中没有的数据所告知的情况下。有时，我们无意中用无法仅从训练数据中获取的信息帮助我们的模型训练，最终导致我们对模型准确性的评估过于乐观。

数据科学家并不真的希望这种情况发生，因此有了“泄露”这个术语。这不是一种“不要这样做”的讨论。我们都知道不要这样做。这更像是一种“我应该采取哪些步骤来避免这个问题？”的讨论。实际上，除非我们制定预防措施，否则很容易出现数据泄露。

例如，如果我们有一个特征的缺失值，我们可能会使用整个数据集的平均值来插补这些值。然而，为了验证我们的模型，我们随后将数据分为训练和测试数据集。这样我们就会意外地将来自完整数据集（即全局平均值）的信息引入到训练数据集中。

数据科学家为了避免这种情况采取的一种做法是在分析开始尽可能早地建立单独的训练和测试数据集。在交叉验证等验证技术中，这可能会变得稍微复杂一些，但在接下来的章节中，我们将介绍如何在各种情况下避免数据泄露。

我们可以使用 scikit-learn 为国家纵向青年调查数据创建训练和测试 DataFrame。

注意

国家纵向青年调查（NLS）由美国劳工统计局进行。这项调查始于 1997 年，调查对象为 1980 年至 1985 年间出生的一批人，每年进行一次年度跟踪调查，直至 2017 年。在本节中，我从调查中的数百个数据项中提取了 89 个关于成绩、就业、收入和对政府态度的变量。可以从存储库下载 SPSS、Stata 和 SAS 的单独文件。NLS 数据可以从www.nlsinfo.org/investigator/pages/search下载供公众使用。

让我们开始创建 DataFrame：

首先，我们从sklearn导入train_test_split模块并加载 NLS 数据：

import pandas as pd
from sklearn.model_selection import train_test_split
nls97 = pd.read_csv("data/nls97b.csv")
nls97.set_index("personid", inplace=True)

然后，我们可以为特征（X_train和X_test）和目标（y_train和y_test）创建训练和测试 DataFrame。在本例中，wageincome是目标变量。我们将test_size参数设置为0.3，以保留 30%的观测值用于测试。请注意，我们只将使用 NLS 中的学术能力评估测试（SAT）和平均成绩点（GPA）数据：
```
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97[feature_cols],\
  nls97[['wageincome']], test_size=0.3, \
  random_state=0)
```

让我们看看使用train_test_split创建的训练 DataFrame。我们得到了预期的观测数，6,288，这是 NLS DataFrame 中 8,984 个观测总数的 70%：

nls97.shape[0]
8984
X_train.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 6288 entries, 574974 to 370933
Data columns (total 6 columns):
 #   Column        Non-Null Count     Dtype
---  ------        --------------   -------
 0   satverbal      1001 non-null   float64
 1   satmath      1001 non-null   float64
 2   gpascience     3998 non-null   float64
 3   gpaenglish     4078 non-null   float64
 4   gpamath        4056 non-null   float64
 5   gpaoverall     4223 non-null   float64
dtypes: float64(6)
memory usage: 343.9 KB
y_train.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 6288 entries, 574974 to 370933
Data columns (total 1 columns):
 #   Column        Non-Null Count    Dtype
---  ------        --------------  -------
 0   wageincome    3599 non-null   float64
dtypes: float64(1)
memory usage: 98.2 KB

此外，让我们看看测试 DataFrame。我们得到了预期的 30%的观测总数：

X_test.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 2696 entries, 363170 to 629736
Data columns (total 6 columns):
 #   Column        Non-Null Count    Dtype
---  ------        --------------  -------  
 0   satverbal      405 non-null   float64
 1   satmath        406 non-null   float64
 2   gpascience    1686 non-null   float64
 3   gpaenglish    1720 non-null   float64
 4   gpamath       1710 non-null   float64
 5   gpaoverall    1781 non-null   float64
dtypes: float64(6)
memory usage: 147.4 KB
y_test.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 2696 entries, 363170 to 629736
Data columns (total 1 columns):
 #   Column          Non-Null Count    Dtype
---  ------          --------------  -------  
 0   wageincome      1492 non-null   float64
dtypes: float64(1)
memory usage: 42.1 KB

我们将在本章的其余部分使用 scikit-learn 的test_train_split来创建单独的训练和测试 DataFrame。我们将在第六章 准备模型评估中介绍构建验证测试数据集的更复杂策略。

接下来，我们开始我们的特征工程工作，通过移除明显无用的特征。这是因为它们与另一个特征具有相同的数据，或者响应中没有变化。

移除冗余或无用的特征

在数据清洗和处理的过程中，我们经常会得到不再有意义的数据。也许我们根据单个特征值对数据进行子集划分，尽管现在所有观测值都具有相同的值，我们仍然保留了该特征。或者，对于我们所使用的数据子集，两个特征具有相同的值。理想情况下，我们在数据清洗过程中捕捉到这些冗余。然而，如果我们在这个过程中没有捕捉到它们，我们可以使用开源的feature-engine包来帮助我们。

此外，可能存在高度相关的特征，我们几乎不可能构建一个能够有效使用所有这些特征的模型。feature-engine有一个名为DropCorrelatedFeatures的方法，它使得在特征高度相关时移除特征变得容易。

在本节中，我们将处理陆地温度数据，以及 NLS 数据。请注意，我们在这里只加载波兰的温度数据。

数据备注

陆地温度数据集包含了 2019 年来自全球超过 12,000 个站点的平均温度读数（以摄氏度为单位），尽管大多数站点位于美国。原始数据是从全球历史气候学网络集成数据库中检索的。它已经由美国国家海洋和大气管理局在www.ncdc.noaa.gov/data-access/land-based-station-data/land-based-datasets/global-historical-climatology-network-monthly-version-4上提供给公众使用。

让我们开始移除冗余和无用的特征：

让我们从feature_engine和sklearn模块中导入所需的模块，并加载波兰的 NLS 数据和温度数据。波兰的数据是从全球 12,000 个气象站的大数据集中提取的。我们使用dropna来删除任何缺失数据的观测值：

import pandas as pd
import feature_engine.selection as fesel
from sklearn.model_selection import train_test_split
nls97 = pd.read_csv("data/nls97b.csv")
nls97.set_index("personid", inplace=True)
ltpoland = pd.read_csv("data/ltpoland.csv")
ltpoland.set_index("station", inplace=True)
ltpoland.dropna(inplace=True)

接下来，我们创建训练和测试 DataFrame，就像我们在上一节中所做的那样：

feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97[feature_cols],\
  nls97[['wageincome']], test_size=0.3, \
  random_state=0)

我们可以使用 pandas 的corr方法来查看这些特征之间的相关性：

X_train.corr()
          satverbal  satmath  gpascience  gpaenglish \
satverbal     1.000    0.729       0.439      0.444  
satmath       0.729    1.000       0.480      0.430
gpascience    0.439    0.480       1.000      0.672
gpaenglish    0.444    0.430       0.672      1.000
gpamath       0.375    0.518       0.606      0.600
gpaoverall    0.421    0.485       0.793      0.844
            gpamath   gpaoverall  
satverbal     0.375        0.421  
satmath       0.518        0.485  
gpascience    0.606        0.793  
gpaenglish    0.600        0.844  
gpamath       1.000        0.750  
gpaoverall    0.750        1.000

在这里，gpaoverall与gpascience、gpaenglish和gpamath高度相关。corr方法默认返回皮尔逊相关系数。当我们假设特征之间存在线性关系时，这是可以的。然而，当这个假设没有意义时，我们应该考虑请求 Spearman 相关系数。我们可以通过将spearman传递给corr方法的参数来实现这一点。

让我们删除与另一个特征相关性高于 0.75 的特征。我们将 0.75 传递给DropCorrelatedFeatures的threshold参数，表示我们想要使用皮尔逊相关系数，并且我们想要通过将变量设置为None来评估所有特征。我们在训练数据上使用fit方法，然后转换训练和测试数据。info方法显示，结果训练 DataFrame（X_train_tr）除了gpaoverall以外的所有特征，gpaoverall与gpascience和gpaenglish的相关性分别为 0.793 和 0.844（DropCorrelatedFeatures将从左到右进行评估，因此如果gpamath和gpaoverall高度相关，它将删除gpaoverall。如果gpaoverall在gpamath左侧，它将删除gpamath）：

tr = fesel.DropCorrelatedFeatures(variables=None, method='pearson', threshold=0.75)
tr.fit(X_train)
X_train_tr = tr.transform(X_train)
X_test_tr = tr.transform(X_test)
X_train_tr.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 6288 entries, 574974 to 370933
Data columns (total 5 columns):
 #   Column       Non-Null Count     Dtype
---  ------        --------------  -------  
 0   satverbal     1001 non-null   float64
 1   satmath       1001 non-null   float64
 2   gpascience    3998 non-null   float64
 3   gpaenglish    4078 non-null   float64
 4   gpamath       4056 non-null   float64
dtypes: float64(5)
memory usage: 294.8 KB

通常，我们在决定删除特征之前会更仔细地评估特征。然而，有时特征选择是管道的一部分，我们需要自动化这个过程。这可以通过DropCorrelatedFeatures来实现，因为所有的feature_engine方法都可以被纳入 scikit-learn 管道。

现在，让我们从波兰的陆地温度数据中创建训练和测试 DataFrame。year的值对所有观测值都是相同的，country的值也是如此。此外，对于每个观测值，latabs的值与latitude相同：

feature_cols = ['year','month','latabs', 
  'latitude','elevation', 'longitude','country']
X_train, X_test, y_train, y_test =  \
  train_test_split(ltpoland[feature_cols],\
  ltpoland[['temperature']], test_size=0.3, \
  random_state=0)
X_train.sample(5, random_state=99)
         year  month  latabs  latitude  elevation  longitude country
station  
SIEDLCE   2019   11    52    52    152    22    Poland
OKECIE    2019    6    52    52    110    21    Poland
BALICE    2019    1    50    50    241    20    Poland
BALICE    2019    7    50    50    241    20    Poland
BIALYSTOK 2019   11    53    53    151    23    Poland
X_train.year.value_counts()
2019    84
Name: year, dtype: int64
X_train.country.value_counts()
Poland    84
Name: country, dtype: int64
(X_train.latitude!=X_train.latabs).sum()
0

让我们删除在整个训练数据集中具有相同值的特征。注意，在转换后删除了year和country：

tr = fesel.DropConstantFeatures()
tr.fit(X_train)
X_train_tr = tr.transform(X_train)
X_test_tr = tr.transform(X_test)
X_train_tr.head()

         month  latabs  latitude  elevation  longitude
station                                                 
OKECIE       1      52        52        110         21
LAWICA       8      52        52         94         17
LEBA        11      55        55          2         18
SIEDLCE     10      52        52        152         22
BIALYSTOK   11      53        53        151         23

让我们删除具有与其他特征相同值的特征。在这种情况下，转换删除了latitude，因为它与latabs具有相同的值：

tr = fesel.DropDuplicateFeatures()
tr.fit(X_train_tr)
X_train_tr = tr.transform(X_train_tr)
X_train_tr.head()
           month   latabs   elevation   longitude
station                                       
OKECIE         1       52         110          21
LAWICA         8       52          94          17
LEBA          11       55           2          18
SIEDLCE       10       52         152          22
BIALYSTOK     11       53         151          23

这解决了 NLS 数据中我们的特征和波兰陆地温度数据中的一些明显问题。我们从包含其他 GPA 特征的 DataFrame 中删除了gpaoverall，因为它与它们高度相关。此外，我们删除了冗余数据，删除了在整个 DataFrame 中具有相同值的特征以及重复另一个特征值的特征。

本章的其余部分探讨了某些较为混乱的特征工程挑战：编码、转换、分箱和缩放。

对分类特征进行编码

我们可能需要在使用大多数机器学习算法之前对特征进行编码的几个原因。首先，这些算法通常需要数值数据。其次，当一个分类特征是用数字表示时，例如，女性为 1，男性为 2，我们需要对这些值进行编码，以便它们被识别为分类数据。第三，该特征实际上可能是有序的，具有代表某些有意义排名的离散数值。我们的模型需要捕捉这种排名。最后，一个分类特征可能具有大量值（称为高基数），我们可能希望我们的编码能够合并类别。

我们可以使用独热编码来处理具有有限值的特征，例如 15 或更少。在本节中，我们首先将介绍独热编码，然后讨论顺序编码。在下一节中，我们将探讨处理具有高基数类别特征的战略。

独热编码

对特征进行独热编码会为该特征的每个值创建一个二进制向量。因此，如果一个名为 letter 的特征有三个唯一值，A，B 和 C，独热编码会创建三个二进制向量来表示这些值。第一个二进制向量，我们可以称之为 letter_A，当 letter 的值为 A 时为 1，而当它是 B 或 C 时为 0。letter_B 和 letter_C 的编码方式类似。转换后的特征，letter_A，letter_B 和 letter_C，通常被称为虚拟变量。图 4.1 展示了独热编码：

图 4.1 – 类别特征的独热编码

NLS 数据中的许多特征适合进行独热编码。在下面的代码块中，我们将对其中一些特征进行编码：

让我们从导入 feature_engine 中的 OneHotEncoder 模块并加载数据开始。此外，我们还从 scikit-learn 中导入 OrdinalEncoder 模块，因为我们稍后会使用它：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.preprocessing import OrdinalEncoder
from sklearn.model_selection import train_test_split
nls97 = pd.read_csv("data/nls97b.csv")
nls97.set_index("personid", inplace=True)

接下来，我们为 NLS 数据创建训练和测试 DataFrame：

feature_cols =['gender','maritalstatus','colenroct99']
nls97demo = nls97[['wageincome'] + feature_cols].dropna()
X_demo_train, X_demo_test, y_demo_train, y_demo_test=\
  train_test_split(nls97demo[feature_cols],\
  nls97demo[['wageincome']], test_size=0.3, \
  random_state=0)

我们用于编码的一个选项是 pandas 的 get_dummies 方法。我们可以用它来指示我们想要转换 gender 和 maritalstatus 特征。get_dummies 为 gender 和 maritalstatus 的每个值提供一个虚拟变量。例如，gender 有 Female 和 Male 的值。get_dummies 创建一个特征，gender_Female，当 gender 为 Female 时为 1，而当 gender 为 Male 时为 0。当 gender 为 Male 时，gender_Male 为 1 而 gender_Female 为 0。这是一个经过验证的方法来进行此类编码，并且多年来一直为统计学家提供了良好的服务：

pd.get_dummies(X_demo_train, \
  columns=['gender','maritalstatus']).head(2).T
personid                    736081          832734
colenroct99                 1.Not enrolled  1.Not enrolled
gender_Female               1               0
gender_Male                 0               1
maritalstatus_Divorced      0               0
maritalstatus_Married       1               0
maritalstatus_Never-married 0               1
maritalstatus_Separated     0               0
maritalstatus_Widowed       0               0

我们没有保存由 get_dummies 创建的 DataFrame，因为在本节的后面部分，我们将使用不同的技术进行编码。

通常，我们为特征的 k 个唯一值创建 k-1 个虚拟变量。因此，如果 gender 在我们的数据中有两个值，我们只需要创建一个虚拟变量。如果我们知道 gender_Female 的值，我们也知道 gender_Male 的值；因此，后者变量是冗余的。同样，如果我们知道其他 maritalstatus 虚拟变量的值，我们也知道 maritalstatus_Divorced 的值。以这种方式创建冗余被不优雅地称为虚拟变量陷阱。为了避免这个问题，我们从每个组中删除一个虚拟变量。

注意

对于某些机器学习算法，例如线性回归，实际上需要删除一个虚拟变量。在估计线性模型的参数时，矩阵会被求逆。如果我们的模型有截距，并且所有虚拟变量都被包含在内，那么矩阵就无法求逆。

我们可以将get_dummies的drop_first参数设置为True以从每个组中删除第一个虚拟变量：

pd.get_dummies(X_demo_train, \
  columns=['gender','maritalstatus'],
  drop_first=True).head(2).T
personid                      736081          832734
colenroct99                  1\. Not enrolled  1\. Not enrolled
gender_Male                  0                1
maritalstatus_Married        1                0
maritalstatus_Never-married  0                1
maritalstatus_Separated      0                0
maritalstatus_Widowed        0                0

get_dummies的一个替代方案是sklearn或feature_engine中的 one-hot 编码器。这些 one-hot 编码器有优势，它们可以轻松地集成到机器学习流程中，并且可以将从训练数据集中收集到的信息持久化到测试数据集中。

让我们使用feature_engine中的OneHotEncoder模块来进行编码。我们将drop_last设置为True以从每个组中删除一个虚拟变量。我们将编码拟合到训练数据，然后转换训练数据和测试数据：

ohe = OneHotEncoder(drop_last=True,
  variables=['gender','maritalstatus'])
ohe.fit(X_demo_train)
X_demo_train_ohe = ohe.transform(X_demo_train)
X_demo_test_ohe = ohe.transform(X_demo_test)
X_demo_train_ohe.filter(regex='gen|mar', axis="columns").head(2).T
personid                     736081          832734
gender_Female                1               0
maritalstatus_Married        1               0
maritalstatus_Never-married  0               1
maritalstatus_Divorced       0               0
maritalstatus_Separated      0               0

这表明 one-hot 编码是准备名义数据供机器学习算法使用的一种相当直接的方法。但如果我们分类特征是有序的而不是名义的，那会怎样？在这种情况下，我们需要使用有序编码。

有序编码

如同在第一章中讨论的，分类特征可以是名义的或有序的。性别和婚姻状况是名义的。它们的值不表示顺序。例如，“未婚”并不比“离婚”的值高。

然而，当一个分类特征是有序的时，我们希望编码能够捕捉到值的排序。例如，如果我们有一个具有低、中、高值的特征，one-hot 编码会丢失这个排序。相反，一个具有低、中、高分别为 1、2、3 的转换特征会更好。我们可以通过有序编码来实现这一点。

NLS 数据集中的大学入学特征可以被认为是有序特征。其值范围从1. 未入学到3. 四年制大学。我们应该使用有序编码来为建模做准备。我们将在下一步做这件事：

我们可以使用sklearn的OrdinalEncoder模块来对 1999 年的大学入学特征进行编码。首先，让我们看一下编码前的colenroct99的值。这些值是字符串，但存在隐含的顺序：

X_demo_train.colenroct99.unique()
array(['1\. Not enrolled', '2\. 2-year college ', 
       '3\. 4-year college'], dtype=object)
X_demo_train.head()
            gender   maritalstatus   colenroct99
personid                                           
736081      Female   Married         1\. Not enrolled
832734      Male     Never-married   1\. Not enrolled
453537      Male     Married         1\. Not enrolled
322059      Female   Divorced        1\. Not enrolled
324323      Female   Married         2\. 2-year college

我们可以通过将前面的数组传递给categories参数来告诉OrdinalEncoder模块按相同的顺序对值进行排序。然后，我们可以使用fit_transform来转换大学入学字段colenroct99。（sklearn的OrdinalEncoder模块的fit_transform方法返回一个 NumPy 数组，因此我们需要使用 pandas DataFrame 方法来创建一个 DataFrame。）最后，我们将编码后的特征与训练数据中的其他特征合并：
```
oe = OrdinalEncoder(categories=\
  [X_demo_train.colenroct99.unique()])
colenr_enc = \
  pd.DataFrame(oe.fit_transform(X_demo_train[['colenroct99']]),
    columns=['colenroct99'], index=X_demo_train.index)
X_demo_train_enc = \
  X_demo_train[['gender','maritalstatus']].\
  join(colenr_enc)
```

让我们看看结果 DataFrame 的几个观察结果。此外，我们还应该比较原始大学入学特征计数与转换特征计数：

X_demo_train_enc.head()
             gender       maritalstatus    colenroct99
personid                                    
736081       Female       Married          0
832734       Male         Never-married    0
453537       Male         Married          0
322059       Female       Divorced         0
324323       Female       Married          1
X_demo_train.colenroct99.value_counts().sort_index()
1\. Not enrolled        3050
2\. 2-year college       142
3\. 4-year college       350
Name: colenroct99, dtype: int64
X_demo_train_enc.colenroct99.value_counts().sort_index()
0       3050
1       142
2       350
Name: colenroct99, dtype: int64

序列编码将 colenroct99 的初始值替换为从 0 到 2 的数字。现在它以许多机器学习模型可消费的形式存在，并且我们保留了有意义的排名信息。

注意

序列编码适用于非线性模型，如决策树。在线性回归模型中可能没有意义，因为这会假设在整个分布中值之间的距离具有同等意义。在本例中，这会假设从 0 到 1 的增加（即从无入学到 2 年入学）与从 1 到 2 的增加（即从 2 年入学到 4 年入学）是同一件事。

One-hot 编码和序列编码是工程化分类特征的相对直接的方法。当有更多唯一值时，处理分类特征可能更复杂。在下一节中，我们将介绍处理这些特征的一些技术。

对中等或高基数的分类特征进行编码

当我们处理具有许多唯一值的分类特征时，例如 10 个或更多，为每个值创建虚拟变量可能是不切实际的。当基数高，即具有非常多的唯一值时，某些值可能观察到的样本太少，无法为我们提供很多信息。在极端情况下，对于 ID 变量，每个值只有一个观察值。

处理中等或高基数的方法有几个。一种方法是为前 k 个类别创建虚拟变量，并将剩余的值组合到一个其他类别中。另一种方法是使用特征哈希，也称为哈希技巧。在本节中，我们将探讨这两种策略。我们将使用 COVID-19 数据集作为示例：

让我们从 COVID-19 数据中创建训练和测试 DataFrame，并导入 feature_engine 和 category_encoders 库：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from category_encoders.hashing import HashingEncoder
from sklearn.model_selection import train_test_split
covidtotals = pd.read_csv("data/covidtotals.csv")
feature_cols = ['location','population',
    'aged_65_older','diabetes_prevalence','region']
covidtotals = covidtotals[['total_cases'] + feature_cols].dropna()
X_train, X_test, y_train, y_test =  \
  train_test_split(covidtotals[feature_cols],\
  covidtotals[['total_cases']], test_size=0.3, 
  random_state=0)

特征区域有 16 个唯一值，其中前 6 个的计数为 10 或更多：

X_train.region.value_counts()
Eastern Europe  16
East Asia  12
Western Europe  12
West Africa  11
West Asia  10
East Africa  10
South America  7
South Asia  7
Central Africa  7
Southern Africa  7
Oceania / Aus  6
Caribbean  6
Central Asia  5
North Africa  4
North America  3
Central America  3
Name: region, dtype: int64

我们可以再次使用 feature_engine 中的 OneHotEncoder 模块来编码 region 特征。这次，我们使用 top_categories 参数来指示我们只想为前六个类别值创建虚拟变量。任何不属于前六个的值都将为所有虚拟变量设置 0：

ohe = OneHotEncoder(top_categories=6, variables=['region'])
covidtotals_ohe = ohe.fit_transform(covidtotals)
covidtotals_ohe.filter(regex='location|region',
  axis="columns").sample(5, random_state=99).T
          97      173      92         187        104
Location  Israel  Senegal  Indonesia  Sri Lanka  Kenya
region_Eastern Europe  0      0      0      0      0
region_Western Europe  0      0      0      0      0
region_West Africa     0      1      0      0      0
region_East Asia       0      0      1      0      0
region_West Asia       1      0      0      0      0
region_East Africa     0      0      0      0      1

当分类特征具有许多唯一值时，一种替代 one-hot 编码的方法是使用 特征哈希。

特征哈希

特征哈希将大量的唯一特征值映射到更少的虚拟变量。我们可以指定要创建的虚拟变量的数量。然而，可能出现冲突；也就是说，一些特征值可能映射到相同的虚拟变量组合。随着我们减少请求的虚拟变量数量，冲突的数量会增加。

我们可以使用 category_encoders 中的 HashingEncoder 进行特征哈希。我们使用 n_components 来表示我们想要六个虚拟变量（我们在变换之前复制了 region 特征，这样我们就可以将原始值与新的虚拟变量进行比较）：

X_train['region2'] = X_train.region

he = HashingEncoder(cols=['region'], n_components=6)

X_train_enc = he.fit_transform(X_train)

X_train_enc.\

 groupby(['col_0','col_1','col_2','col_3','col_4',

   'col_5','region2']).\

    size().reset_index().rename(columns={0:'count'})

  col_0 col_1 col_2 col_3 col_4 col_5 region2         count

0   0     0     0     0     0     1   Caribbean       6

1   0     0     0     0     0     1   Central Africa  7

2   0     0     0     0     0     1   East Africa     10

3   0     0     0     0     0     1   North Africa    4

4   0     0     0     0     1     0   Central America 3

5   0     0     0     0     1     0   Eastern Europe  16

6   0     0     0     0     1     0   North America   3

7   0     0     0     0     1     0   Oceania / Aus   6

8   0     0     0     0     1     0   Southern Africa 7

9   0     0     0     0     1     0   West Asia       10

10  0     0     0     0     1     0   Western Europe  12

11  0     0     0     1     0     0   Central Asia    5

12  0     0     0     1     0     0   East Asia       12

13  0     0     0     1     0     0   South Asia      7

14  0     0     1     0     0     0   West Africa     11

15  1     0     0     0     0     0   South America   7

不幸的是，这给我们带来了大量的冲突。例如，加勒比海、中非、东非和北非都得到了相同的虚拟变量值。在这种情况下，至少使用独热编码并指定类别数量，就像我们在上一节中所做的那样，是一个更好的解决方案。

在前两节中，我们介绍了常见的编码策略：独热编码、顺序编码和特征哈希。我们的大部分分类特征在使用模型之前都需要进行某种形式的编码。然而，有时我们需要以其他方式修改我们的特征，包括变换、分箱和缩放。在接下来的三个部分中，我们将考虑我们可能需要以这种方式修改特征的原因，并探讨实现这些修改的工具。

使用数学变换

有时，我们希望使用不具有高斯分布的特征，而机器学习算法假设我们的特征是以这种方式分布的。当这种情况发生时，我们可能需要改变我们关于使用哪种算法的想法（例如，我们可以选择 KNN 而不是线性回归）或者变换我们的特征，使它们近似于高斯分布。在本节中，我们将介绍几种实现后者的策略：

我们首先从 feature_engine 导入变换模块，从 sklearn 导入 train_test_split，从 scipy 导入 stats。此外，我们使用 COVID-19 数据创建训练和测试 DataFrame：

import pandas as pd
from feature_engine import transformation as vt
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
from scipy import stats
covidtotals = pd.read_csv("data/covidtotals.csv")
feature_cols = ['location','population',
    'aged_65_older','diabetes_prevalence','region']
covidtotals = covidtotals[['total_cases'] + feature_cols].dropna()
X_train, X_test, y_train, y_test =  \
  train_test_split(covidtotals[feature_cols],\
  covidtotals[['total_cases']], test_size=0.3, \
  random_state=0)

让我们看看按国家划分的病例总数是如何分布的。我们还应该计算偏斜：

y_train.total_cases.skew()
6.313169268923333
plt.hist(y_train.total_cases)
plt.title("Total COVID Cases  (in millions)")
plt.xlabel('Cases')
plt.ylabel("Number of Countries")
plt.show()

这产生了以下直方图：

图 4.2 – COVID 病例总数的直方图

这说明了病例总数的极高偏斜。实际上，它看起来是对数正态分布的，考虑到有大量非常低的值和几个非常高的值，这并不令人惊讶。

注意

有关偏斜和峰度的度量方法更多信息，请参阅第一章，检查特征和目标分布。

让我们尝试对数变换。我们只需要调用LogTranformer并传递我们想要变换的特征或特征即可：

tf = vt.LogTransformer(variables = ['total_cases'])
y_train_tf = tf.fit_transform(y_train)
y_train_tf.total_cases.skew()
-1.3872728024141519
plt.hist(y_train_tf.total_cases)
plt.title("Total COVID Cases (log transformation)")
plt.xlabel('Cases')
plt.ylabel("Number of Countries")
plt.show()

这会产生以下直方图：

图 4.3 – 对数变换后的总 COVID 病例数直方图

实际上，对数变换会增加分布下端的变异性，并减少上端的变异性。这会产生一个更对称的分布。这是因为对数函数的斜率对于较小的值比较大的值更陡峭。

这确实是一个很大的改进，但现在有一些负偏斜。也许 Box-Cox 变换会产生更好的结果。让我们试试：

tf = vt.BoxCoxTransformer(variables = ['total_cases'])
y_train_tf = tf.fit_transform(y_train)
y_train_tf.total_cases.skew()
0.07333475786753735
plt.hist(y_train_tf.total_cases)
plt.title("Total COVID Cases (Box-Cox transformation)")
plt.xlabel('Cases')
plt.ylabel("Number of Countries")
plt.show()

这会产生以下图表：

图 4.4 – Box-Cox 变换后的总 COVID 病例数直方图

Box-Cox 变换确定一个介于-5 和 5 之间的 lambda 值，该值生成一个与正态分布最接近的分布。它使用以下方程进行变换：

或者

在这里，是我们的变换特征。为了好玩，让我们看看用于变换total_cases的 lambda 值：

stats.boxcox(y_train.total_cases)[1]

0.10435377585681517

Box-Cox 变换的 lambda 值为0.104。相比之下，具有高斯分布的特征的 lambda 值为 1.000，这意味着不需要进行变换。

现在我们转换后的总病例特征看起来很好，我们可以用它作为目标来构建模型。此外，我们可以在预测时设置我们的管道以将值恢复到原始缩放。feature_engine有其他一些变换，它们的实现方式类似于对数和 Box-Cox 变换。

特征分箱

有时，我们可能希望将一个连续特征转换为分类特征。从分布的最小值到最大值创建k个等间隔区间的过程称为分箱，或者，不那么友好的术语，离散化。分箱可以解决特征的一些重要问题：偏斜、过度峰度和异常值的存在。

等宽和等频分箱

在 COVID 病例数据中，分箱可能是一个不错的选择。让我们试试（这可能在数据集中的其他变量中也很有用，包括总死亡人数和人口，但我们现在只处理总病例数。total_cases是以下代码中的目标变量，因此它是一个列——y_train DataFrame 上的唯一列）：

首先，我们需要从feature_engine导入EqualFrequencyDiscretiser和EqualWidthDiscretiser。此外，我们还需要从 COVID 数据创建训练集和测试集 DataFrame：

import pandas as pd
from feature_engine.discretisation import EqualFrequencyDiscretiser as efd
from feature_engine.discretisation import EqualWidthDiscretiser as ewd
from sklearn.preprocessing import KBinsDiscretizer
from sklearn.model_selection import train_test_split
covidtotals = pd.read_csv("data/covidtotals.csv")
feature_cols = ['location','population',
    'aged_65_older','diabetes_prevalence','region']
covidtotals = covidtotals[['total_cases'] + feature_cols].dropna()
X_train, X_test, y_train, y_test =  \
  train_test_split(covidtotals[feature_cols],\
  covidtotals[['total_cases']], test_size=0.3, random_state=0)

我们可以使用 pandas 的qcut方法和其q参数来创建 10 个相对等频的箱子：

y_train['total_cases_group'] = pd.qcut(y_train.total_cases, q=10, labels=[0,1,2,3,4,5,6,7,8,9])
y_train.total_cases_group.value_counts().sort_index()
0   13
1   13
2   12
3   13
4   12
5   13
6   12
7   13
8   12
9   13
Name: total_cases_group, dtype: int64

我们可以使用EqualFrequencyDiscretiser实现相同的功能。首先，我们定义一个函数来运行转换。该函数接受一个feature_engine转换和训练集和测试集 DataFrame。它返回转换后的 DataFrame（定义函数不是必需的，但在这里这样做是有意义的，因为我们稍后会重复这些步骤）：
```
def runtransform(bt, dftrain, dftest):
  bt.fit(dftrain)
  train_bins = bt.transform(dftrain)
  test_bins = bt.transform(dftest)
  return train_bins, test_bins
```

接下来，我们创建一个EqualFrequencyDiscretiser转换器，并调用我们刚刚创建的runtransform函数：

y_train.drop(['total_cases_group'], axis=1, inplace=True)
bintransformer = efd(q=10, variables=['total_cases'])
y_train_bins, y_test_bins = runtransform(bintransformer, y_train, y_test)
y_train_bins.total_cases.value_counts().sort_index()
0  13
1  13
2  12
3  13
4  12
5  13
6  12
7  13
8  12
9  13
Name: total_cases, dtype: int64

这给我们带来了与qcut相同的结果，但它有一个优点，即更容易将其引入机器学习流程，因为我们使用feature_engine来生成它。等频分箱解决了偏斜和异常值问题。

注意

我们将在本书中详细探讨机器学习流程，从第六章，准备模型评估开始。在这里，关键点是特征引擎转换器可以是包含其他sklearn兼容转换器的流程的一部分，甚至包括我们自己构建的。

EqualWidthDiscretiser的工作方式类似：

bintransformer = ewd(bins=10, variables=['total_cases'])
y_train_bins, y_test_bins = runtransform(bintransformer, y_train, y_test)
y_train_bins.total_cases.value_counts().sort_index()
0  119
1  4
5  1
9  2
Name: total_cases, dtype: int64

这是一个远不如成功的转换。在分箱之前的数据中，几乎所有值都位于分布的底部，因此等宽分箱会产生相同的问题并不令人惊讶。它只产生了 4 个箱子，尽管我们请求了 10 个。

让我们检查每个箱子的范围。在这里，我们可以看到由于分布顶部的观察值数量很少，等宽分箱器甚至无法构建等宽箱子：

pd.options.display.float_format = '{:,.0f}'.format
y_train_bins = y_train_bins.\
  rename(columns={'total_cases':'total_cases_group'}).\
  join(y_train)
y_train_bins.groupby("total_cases_group")["total_cases"].agg(['min','max'])
  min  max
total_cases_group       
0  1           3,304,135
1  3,740,567   5,856,682
5  18,909,037  18,909,037
9  30,709,557  33,770,444

尽管在这种情况下，等宽分箱是一个糟糕的选择，但很多时候它是有意义的。当数据分布更均匀或等宽在实质上有意义时，它可能很有用。

K-means 分箱

另一个选项是使用 k-means 聚类来确定箱子。k-means 算法随机选择 k 个数据点作为聚类的中心，然后将其他数据点分配到最近的聚类。计算每个聚类的平均值，并将数据点重新分配到最近的新的聚类。这个过程重复进行，直到找到最佳中心。

当使用 k-means 进行分箱时，同一聚类中的所有数据点将具有相同的序数值：

我们可以使用 scikit-learn 的KBinsDiscretizer使用 COVID 病例数据创建箱子：

kbins = KBinsDiscretizer(n_bins=10, encode='ordinal', strategy='kmeans')
y_train_bins = \
  pd.DataFrame(kbins.fit_transform(y_train),
  columns=['total_cases'])
y_train_bins.total_cases.value_counts().sort_index()
0  49
1  24
2  23
3  11
4  6
5  6
6  4
7  1
8  1
9  1
Name: total_cases, dtype: int64

让我们比较原始总病例变量的偏斜和峰度与分箱变量的偏斜和峰度。回想一下，对于一个具有高斯分布的变量，我们预计偏斜为 0，峰度接近 3。分箱变量的分布与高斯分布非常接近：

y_train.total_cases.agg(['skew','kurtosis'])
skew          6.313
kurtosis     41.553
Name: total_cases, dtype: float64
y_train_bins.total_cases.agg(['skew','kurtosis'])
skew            1.439
kurtosis        1.923
Name: total_cases, dtype: float64

分箱可以帮助我们解决数据中的偏斜、峰度和异常值。然而，它确实掩盖了特征中的大部分变化，并减少了其解释潜力。通常，某种形式的缩放，如最小-最大或 z 分数，是一个更好的选择。让我们接下来检查特征缩放。

特征缩放

通常，我们想在模型中使用的特征在非常不同的尺度上。简单来说，最小值和最大值之间的距离，或者说范围，在可能的特征之间有很大的变化。例如，在 COVID-19 数据中，总病例特征从 1 到近 3400 万，而 65 岁及以上的人口从 9% 到 27%（数字代表人口百分比）。

特征尺度差异很大会影响许多机器学习算法。例如，KNN 模型通常使用欧几里得距离，范围更大的特征将对模型产生更大的影响。缩放可以解决这个问题。

在本节中，我们将介绍两种流行的缩放方法：最小-最大缩放和标准（或z 分数）缩放。最小-最大缩放将每个值替换为其在范围内的位置。更确切地说，以下情况发生：

在这里，是最小-最大分数，是观测的特征的值，而和是该特征的最小值和最大值。

标准缩放将特征值标准化到均值为 0 的周围。那些学习过本科统计学的人会将其识别为 z 分数。具体来说，如下所示：

在这里，是特征的观测的值，是特征的均值，而是该特征的标准差。

我们可以使用 scikit-learn 的预处理模块来获取最小-最大和标准缩放器：

我们首先导入预处理模块，并从 COVID-19 数据中创建训练和测试 DataFrame：

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler, StandardScaler, RobustScaler
covidtotals = pd.read_csv("data/covidtotals.csv")
feature_cols = ['population','total_deaths',
    'aged_65_older','diabetes_prevalence']
covidtotals = covidtotals[['total_cases'] + feature_cols].dropna()
X_train, X_test, y_train, y_test =  \
  train_test_split(covidtotals[feature_cols],\
  covidtotals[['total_cases']], test_size=0.3, random_state=0)

现在，我们可以运行最小-最大缩放器。sklearn 的 fit_transform 方法将返回一个 numpy 数组。我们使用训练 DataFrame 的列和索引将其转换为 pandas DataFrame。注意，现在所有特征现在都有介于 0 和 1 之间的值：

scaler = MinMaxScaler()
X_train_mms = pd.DataFrame(scaler.fit_transform(X_train),
  columns=X_train.columns, index=X_train.index)
X_train_mms.describe()
     population total_deaths aged_65_older diabetes_prevalence
count  123.00      123.00        123.00        123.00
mean   0.04        0.04          0.30          0.41
std    0.13        0.14          0.24          0.23
min    0.00        0.00          0.00          0.00
25%    0.00        0.00          0.10          0.26
50%    0.01        0.00          0.22          0.37
75%    0.02        0.02          0.51          0.54
max    1.00        1.00          1.00          1.00

我们以相同的方式运行标准缩放器：

scaler = StandardScaler()
X_train_ss = pd.DataFrame(scaler.fit_transform(X_train),
  columns=X_train.columns, index=X_train.index)
X_train_ss.describe()
       population  total_deaths  aged_65_older  diabetes_prevalence
count  123.00      123.00        123.00       123.00
mean  -0.00       -0.00         -0.00        -0.00
std    1.00        1.00          1.00         1.00
min   -0.29       -0.32         -1.24        -1.84
25%   -0.27       -0.31         -0.84        -0.69
50%   -0.24       -0.29         -0.34        -0.18
75%   -0.11       -0.18          0.87         0.59
max    7.58        6.75          2.93         2.63

如果我们的数据中有异常值，鲁棒缩放可能是一个不错的选择。鲁棒缩放从变量的每个值中减去中位数，并将该值除以四分位距。因此，每个值如下所示：

在这里，是特征的值，而、和分别是特征的均值、第三四分位数和第一四分位数。由于鲁棒缩放不使用均值或方差，因此它对极端值不太敏感。

我们可以使用 scikit-learn 的RobustScaler模块来进行鲁棒缩放：

scaler = RobustScaler()
X_train_rs = pd.DataFrame(
  scaler.fit_transform(X_train),
  columns=X_train.columns, index=X_train.index)
X_train_rs.describe()
     population total_deaths aged_65_older diabetes_prevalence
count  123.00      123.00      123.00      123.00
mean   1.47        2.22        0.20        0.14
std    6.24        7.65        0.59        0.79
min   -0.35       -0.19       -0.53       -1.30
25%   -0.24       -0.15       -0.30       -0.40
50%    0.00        0.00        0.00        0.00
75%    0.76        0.85        0.70        0.60
max    48.59       53.64       1.91        2.20

我们在大多数机器学习算法中使用特征缩放。尽管它不是经常必需的，但它会产生明显更好的结果。最小-最大缩放和标准缩放是流行的缩放技术，但在某些情况下，鲁棒缩放可能是更好的选择。

摘要

在本章中，我们涵盖了广泛的特征工程技术。我们使用了工具来删除冗余或高度相关的特征。我们探讨了最常见的编码类型——独热编码、顺序编码和哈希编码。在此之后，我们使用了转换来改善我们特征的分布。最后，我们使用了常见的分箱和缩放方法来解决偏斜、峰度和异常值，以及调整具有广泛不同范围的特性。

本章中我们讨论的一些技术对于大多数机器学习模型是必需的。我们几乎总是需要为算法编码我们的特征以便正确理解它们。例如，大多数算法无法理解女性或男性值，或者不知道不要将邮编视为有序值。虽然通常不是必需的，但当我们的特征具有非常不同的范围时，缩放通常是一个非常不错的想法。当我们使用假设特征具有高斯分布的算法时，可能需要对特征进行某种形式的转换，以便与该假设保持一致。

现在，我们对特征的分布有了很好的了解，已经填充了缺失值，并在必要时进行了一些特征工程。我们现在准备开始模型构建过程中最有趣和最有意义的一部分——特征选择。

在下一章中，我们将检查关键的特征选择任务，这些任务建立在到目前为止我们所做的特征清洗、探索和工程工作之上。

第五章：第五章: 特征选择

根据你开始数据分析工作和你的个人智力兴趣的不同，你可能会对特征选择这个话题有不同的看法。你可能认为，“嗯，嗯，这是一个重要的主题，但我真的想开始模型构建。”或者，在另一个极端，你可能会认为特征选择是模型构建的核心，并相信一旦你选择了特征，你就已经完成了模型构建的 90%。现在，让我们先达成共识，在我们进行任何严肃的模型指定之前，我们应该花一些时间来理解特征之间的关系——如果我们正在构建监督模型，那么它们与目标之间的关系。

以“少即是多”的态度来处理我们的特征选择工作是有帮助的。如果我们能用更少的特征达到几乎相同的准确度或解释更多的方差，我们应该选择更简单的模型。有时，我们实际上可以用更少的特征获得更好的准确度。这可能会很难理解，甚至对我们这些从构建讲述丰富和复杂故事模型的实践中成长起来的人来说有些令人失望。

但我们在拟合机器学习模型时，对参数估计的关注不如对预测准确性的关注。不必要的特征可能导致过拟合并消耗硬件资源。

有时，我们可能需要花费数月时间来指定模型的特征，即使数据中列的数量有限。例如，在第二章“检查特征与目标之间的双变量和多变量关系”中创建的双变量相关性，给我们一些预期的感觉，但一旦引入其他可能的解释特征，特征的重要性可能会显著变化。该特征可能不再显著，或者相反，只有在包含其他特征时才显著。两个特征可能高度相关，以至于包含两个特征与只包含一个特征相比，提供的额外信息非常有限。

本章将深入探讨适用于各种预测建模任务的特征选择技术。具体来说，我们将探讨以下主题：

为分类模型选择特征
为回归模型选择特征
使用正向和反向特征选择
使用穷举特征选择
在回归模型中递归消除特征
在分类模型中递归消除特征
使用 Boruta 进行特征选择
使用正则化和其他嵌入式方法
使用主成分分析

技术要求

本章中，我们将使用feature_engine、mlxtend和boruta包，以及scikit-learn库。您可以使用pip安装这些包。我选择了一个观测值数量较少的数据集用于本章的工作，因此代码即使在次优工作站上也能正常运行。

注意

在本章中，我们将专门使用美国劳工统计局进行的《青年纵向调查》数据。这项调查始于 1997 年，调查对象为 1980 年至 1985 年间出生的一代人，每年进行一次年度跟踪调查，直至 2017 年。我们将使用教育成就、家庭人口统计、工作周数和工资收入数据。工资收入列代表 2016 年赚取的工资。NLS 数据集可以下载供公众使用，网址为www.nlsinfo.org/investigator/pages/search。

为分类模型选择特征

最直接的特征选择方法是基于每个特征与目标变量的关系。接下来的两个部分将探讨基于特征与目标变量之间的线性或非线性关系来确定最佳k个特征的技术。这些被称为过滤方法。它们有时也被称为单变量方法，因为它们评估特征与目标变量之间的关系，而不考虑其他特征的影响。

当目标变量为分类变量时，我们使用的策略与目标变量为连续变量时有所不同。在本节中，我们将介绍前者，在下一节中介绍后者。

基于分类目标的互信息特征选择

当目标变量为分类变量时，我们可以使用互信息分类或方差分析（ANOVA）测试来选择特征。我们将首先尝试互信息分类，然后进行 ANOVA 比较。

互信息是衡量通过知道另一个变量的值可以获得多少关于变量的信息的度量。在极端情况下，当特征完全独立时，互信息分数为 0。

我们可以使用scikit-learn的SelectKBest类根据互信息分类或其他适当的度量选择具有最高预测强度的k个特征。我们可以使用超参数调整来选择k的值。我们还可以检查所有特征的分数，无论它们是否被识别为k个最佳特征之一，正如我们将在本节中看到的。

让我们先尝试互信息分类来识别与完成学士学位相关的特征。稍后，我们将将其与使用 ANOVA F 值作为选择依据进行比较：

我们首先从feature_engine导入OneHotEncoder来编码一些数据，并从scikit-learn导入train_test_split来创建训练和测试数据。我们还需要scikit-learn的SelectKBest、mutual_info_classif和f_classif模块来进行特征选择：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import SelectKBest,\
  mutual_info_classif, f_classif

我们加载了具有完成学士学位的二进制变量和可能与学位获得相关的特征的数据集：gender特征，并对其他数据进行缩放：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['gender','satverbal','satmath',
  'gpascience', 'gpaenglish','gpamath','gpaoverall',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
X_train_enc = ohe.fit_transform(X_train)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
X_train_enc = \
  pd.DataFrame(scaler.fit_transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])

注意

在本章中，我们将对 NLS 数据进行完整案例分析；也就是说，我们将删除任何特征缺失的观测值。这通常不是一个好的方法，尤其是在数据不是随机缺失或一个或多个特征有大量缺失值时尤其有问题。在这种情况下，最好使用我们在第三章中使用的某些方法，识别和修复缺失值。我们将在本章中进行完整案例分析，以使示例尽可能简单。

现在，我们已经准备好为我们的学士学位完成模型选择特征。一种方法是用互信息分类。为此，我们将SelectKBest的score_func值设置为mutual_info_classif，并指出我们想要五个最佳特征。然后，我们调用fit并使用get_support方法来获取五个最佳特征：
```
ksel = SelectKBest(score_func=mutual_info_classif, k=5)
ksel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[ksel.get_support()]
selcols
Index(['satverbal', 'satmath', 'gpascience', 'gpaenglish', 'gpaoverall'], dtype='object')
```

如果我们还想看到每个特征的得分，我们可以使用scores_属性，尽管我们需要做一些工作来将得分与特定的特征名称关联起来，并按降序排序：

pd.DataFrame({'score': ksel.scores_,
  'feature': X_train_enc.columns},
   columns=['feature','score']).\
   sort_values(['score'], ascending=False)
        feature              score
5       gpaoverall           0.108
1       satmath              0.074
3       gpaenglish           0.072
0       satverbal            0.069
2       gpascience           0.047
4       gpamath              0.038
8       parentincome         0.024
7       fatherhighgrade      0.022
6       motherhighgrade      0.022
9       gender_Female        0.015

注意

这是一个随机过程，所以每次运行它时我们都会得到不同的结果。

为了每次都能得到相同的结果，你可以将一个部分函数传递给score_func：

from functools import partial
SelectKBest(score_func=partial(mutual_info_classif, 
                               random_state=0), k=5)

我们可以使用使用get_support创建的selcols数组来创建仅包含重要特征的 DataFrame。（我们也可以使用SelectKBest的transform方法。这将返回所选特征的值作为 NumPy 数组。）

X_train_analysis = X_train_enc[selcols] 
X_train_analysis.dtypes
satverbal       float64
satmath         float64
gpascience      float64
gpaenglish      float64
gpaoverall      float64
dtype: object

这就是我们使用互信息来选择模型中最佳 k 个特征所需做的所有事情。

使用分类目标的特征选择的 ANOVA F 值

或者，我们可以使用方差分析（ANOVA）而不是互信息。方差分析评估每个目标类中特征的平均值差异。当我们假设特征和目标之间存在线性关系，并且我们的特征是正态分布时，这是一个很好的单变量特征选择指标。如果这些假设不成立，互信息分类是一个更好的选择。

让我们尝试使用 ANOVA 进行特征选择。我们可以将SelectKBest的score_func参数设置为f_classif，以便基于 ANOVA 进行选择：

ksel = SelectKBest(score_func=f_classif, k=5)

ksel.fit(X_train_enc, y_train.values.ravel())

selcols = X_train_enc.columns[ksel.get_support()]

selcols

Index(['satverbal', 'satmath', 'gpascience', 'gpaenglish', 'gpaoverall'], dtype='object')

pd.DataFrame({'score': ksel.scores_,

  'feature': X_train_enc.columns},

   columns=['feature','score']).\

   sort_values(['score'], ascending=False)

       feature                score

5      gpaoverall           119.471

3      gpaenglish           108.006

2      gpascience            96.824

1      satmath               84.901

0      satverbal             77.363

4      gpamath               60.930

7      fatherhighgrade       37.481

6      motherhighgrade       29.377

8      parentincome          22.266

9      gender_Female         15.098

这选择了与我们使用互信息时选择的相同特征。显示得分给我们一些关于所选的k值是否合理的指示。例如，第五到第六个最佳特征的得分下降（77-61）比第四到第五个（85-77）的下降更大。然而，从第六到第七个的下降更大（61-37），这表明我们至少应该考虑k的值为 6。

ANOVA 测试和之前我们做的互信息分类没有考虑在多元分析中仅重要的特征。例如，fatherhighgrade可能在具有相似 GPA 或 SAT 分数的个人中很重要。我们将在本章后面使用多元特征选择方法。在下一节中，我们将进行更多单变量特征选择，以探索适合连续目标的特征选择技术。

选择回归模型的特征

scikit-learn的选择模块在构建回归模型时提供了几个选择特征的选择。在这里，我不指线性回归模型。我只是在指具有连续目标的模型）。两个好的选择是基于 F 检验的选择和基于回归的互信息选择。让我们从 F 检验开始。

基于连续目标的特征选择的 F 检验

F 统计量是目标与单个回归器之间线性相关强度的度量。Scikit-learn有一个f_regression评分函数，它返回 F 统计量。我们可以使用它与SelectKBest一起选择基于该统计量的特征。

让我们使用 F 统计量来选择工资模型的特征。我们将在下一节中使用互信息来选择相同目标的特征：

我们首先从feature_engine导入 one-hot 编码器，从scikit-learn导入train_test_split和SelectKBest。我们还导入f_regression以获取后续的 F 统计量：

import pandas as pd
import numpy as np
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import SelectKBest, f_regression

接下来，我们加载 NLS 数据，包括教育成就、家庭收入和工资收入数据：

nls97wages = pd.read_csv("data/nls97wages.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome',
  'completedba']

然后，我们创建训练和测试数据框，对gender特征进行编码，并对训练数据进行缩放。在这种情况下，我们需要对目标进行缩放，因为它是有连续性的：

X_train, X_test, y_train, y_test =  \
  train_test_split(nls97wages[feature_cols],\
  nls97wages[['wageincome']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
X_train_enc = ohe.fit_transform(X_train)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
X_train_enc = \
  pd.DataFrame(scaler.fit_transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Male']])
y_train = \
  pd.DataFrame(scaler.fit_transform(y_train),
  columns=['wageincome'], index=y_train.index)

注意

你可能已经注意到我们没有对测试数据进行编码或缩放。我们最终需要这样做以验证我们的模型。我们将在本章后面介绍验证，并在下一章中详细介绍。

现在，我们已准备好选择特征。我们将SelectKBest的score_func设置为f_regression，并指出我们想要五个最佳特征。SelectKBest的get_support方法对每个被选中的特征返回True：

ksel = SelectKBest(score_func=f_regression, k=5)
ksel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[ksel.get_support()]
selcols
Index(['satmath', 'gpascience', 'parentincome',
 'completedba','gender_Male'],
      dtype='object')

我们可以使用scores_属性来查看每个特征的得分：

pd.DataFrame({'score': ksel.scores_,
  'feature': X_train_enc.columns},
   columns=['feature','score']).\
   sort_values(['score'], ascending=False)

              feature              score
1             satmath              45
9             completedba          38
10            gender_Male          26
8             parentincome         24
2             gpascience           21
0             satverbal            19
5             gpaoverall           17
4             gpamath              13
3             gpaenglish           10
6             motherhighgrade       9
7             fatherhighgrade       8

F 统计量的缺点是它假设每个特征与目标之间存在线性关系。当这个假设不合理时，我们可以使用互信息进行回归。

对于具有连续目标的特征选择中的互信息

我们还可以使用SelectKBest通过回归中的互信息来选择特征：

我们需要将SelectKBest的score_func参数设置为mutual_info_regression，但存在一个小问题。为了每次运行特征选择时都能得到相同的结果，我们需要设置一个random_state值。正如我们在前一小节中讨论的，我们可以使用一个部分函数来做到这一点。我们将partial(mutual_info_regression, random_state=0)传递给评分函数。

我们可以运行fit方法，并使用get_support来获取选定的特征。我们可以使用scores_属性来为每个特征给出分数：

from functools import partial
ksel = SelectKBest(score_func=\
  partial(mutual_info_regression, random_state=0),
  k=5)
ksel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[ksel.get_support()]
selcols
Index(['satmath', 'gpascience', 'fatherhighgrade', 'completedba','gender_Male'],dtype='object')
pd.DataFrame({'score': ksel.scores_,
  'feature': X_train_enc.columns},
   columns=['feature','score']).\
   sort_values(['score'], ascending=False)
           feature               score
1          satmath               0.101
10         gender_Male           0.074
7          fatherhighgrade       0.047
2          gpascience            0.044
9          completedba           0.044
4          gpamath               0.016
8          parentincome          0.015
6          motherhighgrade       0.012
0          satverbal             0.000
3          gpaenglish            0.000
5          gpaoverall            0.000

我们在回归中的互信息得到了与 F 检验相当相似的结果。parentincome通过 F 检验被选中，而fatherhighgrade通过互信息被选中。否则，选中的特征是相同的。

与 F 检验相比，互信息在回归中的关键优势是它不假设特征与目标之间存在线性关系。如果这个假设被证明是不合理的，互信息是一个更好的方法。（再次强调，评分过程中也存在一些随机性，每个特征的分数可能会在一定范围内波动。）

注意

我们选择k=5以获取五个最佳特征是非常随意的。我们可以通过一些超参数调整使其更加科学。我们将在下一章中介绍调整。

我们迄今为止使用的特征选择方法被称为过滤器方法。它们检查每个特征与目标之间的单变量关系。它们是一个好的起点。类似于我们在前几章中讨论的，在开始检查多元关系之前，拥有相关性的有用性，至少探索过滤器方法是有帮助的。然而，通常我们的模型拟合需要考虑当其他特征也被包含时，哪些特征是重要的，哪些不是。为了做到这一点，我们需要使用包装器或嵌入式方法进行特征选择。我们将在下一节中探讨包装器方法，从前向和后向特征选择开始。

使用前向和后向特征选择

前向和后向特征选择，正如其名称所暗示的，通过逐个添加（或对于后向选择，逐个减去）特征来选择特征，并在每次迭代后评估对模型性能的影响。由于这两种方法都是基于给定的算法来评估性能，因此它们被认为是包装器选择方法。

包装特征选择方法相对于我们之前探索的过滤方法有两个优点。首先，它们在包含其他特征时评估特征的重要性。其次，由于特征是根据其对特定算法性能的贡献来评估的，因此我们能够更好地了解哪些特征最终会起作用。例如，根据我们上一节的结果，satmath似乎是一个重要的特征。但有可能satmath只有在使用特定模型时才重要，比如线性回归，而不是决策树回归等其他模型。包装选择方法可以帮助我们发现这一点。

包装方法的缺点主要在于它们在每次迭代后都会重新训练模型，因此在计算上可能相当昂贵。在本节中，我们将探讨前向和后向特征选择。

使用前向特征选择

前向特征选择首先识别出与目标有显著关系的特征子集，这与过滤方法类似。但它随后评估所有可能的选择特征的组合，以确定与所选算法表现最佳的组合。

我们可以使用前向特征选择来开发一个完成学士学位的模型。由于包装方法要求我们选择一个算法，而这是一个二元目标，因此让我们使用scikit-learn的mlxtend模块中的feature_selection来进行选择特征的迭代：

我们首先导入必要的库：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from mlxtend.feature_selection import SequentialFeatureSelector

然后，我们再次加载 NLS 数据。我们还创建了一个训练 DataFrame，对gender特征进行编码，并对剩余特征进行标准化：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
X_train_enc = ohe.fit_transform(X_train)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
X_train_enc = \
  pd.DataFrame(scaler.fit_transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])

我们创建一个随机森林分类器对象，然后将该对象传递给mlxtend的特征选择器。我们指出我们想要选择五个特征，并且应该进行前向选择。（我们也可以使用顺序特征选择器进行后向选择。）运行fit后，我们可以使用k_feature_idx_属性来获取所选特征的列表：

rfc = RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=0)
sfs = SequentialFeatureSelector(rfc, k_features=5,
  forward=True, floating=False, verbose=2,
  scoring='accuracy', cv=5)
sfs.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[list(sfs.k_feature_idx_)]
selcols
Index(['satverbal', 'satmath', 'gpaoverall',
'parentincome', 'gender_Female'], dtype='object')

你可能还记得本章的第一节，我们针对完成学士学位目标的多变量特征选择给出了不同的结果：

Index(['satverbal', 'satmath', 'gpascience',
 'gpaenglish', 'gpaoverall'], dtype='object')

有三个特征——satmath、satverbal和gpaoverall——是相同的。但我们的前向特征选择已经将parentincome和gender_Female识别为比在单变量分析中选择的gpascience和gpaenglish更重要的特征。实际上，gender_Female在早期分析中的得分最低。这些差异可能反映了包装特征选择方法的优点。我们可以识别出除非包含其他特征，否则不重要的特征，并且我们正在评估对特定算法（在这种情况下是随机森林分类）性能的影响。

前向选择的缺点之一是，一旦选择了特征，它就不会被移除，即使随着更多特征的添加，它的重要性可能会下降。（回想一下，前向特征选择是基于该特征对模型的贡献迭代添加特征的。）

让我们看看我们的结果是否随着反向特征选择而变化。

使用反向特征选择

反向特征选择从所有特征开始，并消除最不重要的特征。然后，它使用剩余的特征重复此过程。我们可以使用 mlxtend 的 SequentialFeatureSelector 以与正向选择相同的方式用于反向选择。

我们从 scikit-learn 库实例化了一个 RandomForestClassifier 对象，然后将其传递给 mlxtend 的顺序特征选择器：

rfc = RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=0)

sfs = SequentialFeatureSelector(rfc, k_features=5,

  forward=False, floating=False, verbose=2,

  scoring='accuracy', cv=5)

sfs.fit(X_train_enc, y_train.values.ravel())

selcols = X_train_enc.columns[list(sfs.k_feature_idx_)]

selcols

Index(['satverbal', 'gpascience', 'gpaenglish',

 'gpaoverall', 'gender_Female'], dtype='object')

也许并不令人惊讶，我们在特征选择上得到了不同的结果。satmath 和 parentincome 不再被选中，而 gpascience 和 gpaenglish 被选中。

反向特征选择与前向特征选择的缺点相反。一旦移除了特征，它就不会被重新评估，即使其重要性可能会随着不同的特征组合而改变。让我们尝试使用穷举特征选择。

使用穷举特征选择

如果你的正向和反向选择的结果没有说服力，而且你不在意在喝咖啡或吃午餐的时候运行模型，你可以尝试穷举特征选择。穷举特征选择会在所有可能的特征组合上训练给定的模型，并选择最佳的特征子集。但这也需要付出代价。正如其名所示，这个过程可能会耗尽系统资源和你的耐心。

让我们为学士学位完成情况的模型使用穷举特征选择：

我们首先加载所需的库，包括来自 scikit-learn 的 RandomForestClassifier 和 LogisticRegression 模块，以及来自 mlxtend 的 ExhaustiveFeatureSelector。我们还导入了 accuracy_score 模块，这样我们就可以使用选定的特征来评估模型：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from mlxtend.feature_selection import ExhaustiveFeatureSelector
from sklearn.metrics import accuracy_score

接下来，我们加载 NLS 教育达成度数据，并创建训练和测试 DataFrame：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)

然后，我们对训练和测试数据进行编码和缩放：

ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

我们创建了一个随机森林分类器对象，并将其传递给 mlxtend 的 ExhaustiveFeatureSelector。我们告诉特征选择器评估所有一至五个特征的组合，并返回预测学位达成度最高的组合。运行 fit 后，我们可以使用 best_feature_names_ 属性来获取选定的特征：

rfc = RandomForestClassifier(n_estimators=100, max_depth=2,n_jobs=-1, random_state=0)
efs = ExhaustiveFeatureSelector(rfc, max_features=5,
  min_features=1, scoring='accuracy', 
  print_progress=True, cv=5)
efs.fit(X_train_enc, y_train.values.ravel())
efs.best_feature_names_
('satverbal', 'gpascience', 'gpamath', 'gender_Female')

让我们评估这个模型的准确性。我们首先需要将训练和测试数据转换为只包含四个选定的特征。然后，我们可以仅使用这些特征再次拟合随机森林分类器，并生成学士学位完成情况的预测值。然后，我们可以计算我们正确预测目标的时间百分比，这是 67%：

X_train_efs = efs.transform(X_train)
X_test_efs = efs.transform(X_test)
rfc.fit(X_train_efs, y_train.values.ravel())
y_pred = rfc.predict(X_test_efs)
confusion = pd.DataFrame(y_pred, columns=['pred'],
  index=y_test.index).\
  join(y_test)
confusion.loc[confusion.pred==confusion.completedba].shape[0]\
  /confusion.shape[0]
0.6703296703296703

如果我们只使用 scikit-learn 的accuracy score，我们也会得到相同的答案。（我们在上一步计算它，因为它相当直接，并且让我们更好地理解在这种情况下准确率的含义。）
```
accuracy_score(y_test, y_pred)
0.6703296703296703
```
注意

准确率分数通常用于评估分类模型的性能。在本章中，我们将依赖它，但根据您模型的目的，其他指标可能同样重要或更重要。例如，我们有时更关心灵敏度，即我们的正确阳性预测与实际阳性数量的比率。我们在第六章中详细探讨了分类模型的评估，准备模型评估。

现在我们尝试使用逻辑模型进行全面特征选择：

lr = LogisticRegression(solver='liblinear')
efs = ExhaustiveFeatureSelector(lr, max_features=5,
  min_features=1, scoring='accuracy', 
  print_progress=True, cv=5)
efs.fit(X_train_enc, y_train.values.ravel())
efs.best_feature_names_
('satmath', 'gpascience', 'gpaenglish', 'motherhighgrade', 'gender_Female')

让我们看看逻辑模型的准确率。我们得到了相当相似的准确率分数：

X_train_efs = efs.transform(X_train_enc)
X_test_efs = efs.transform(X_test_enc)
lr.fit(X_train_efs, y_train.values.ravel())
y_pred = lr.predict(X_test_efs)
accuracy_score(y_test, y_pred)
0.6923076923076923

逻辑模型的一个关键优势是它训练得更快，这对于全面特征选择来说确实有很大影响。如果我们为每个模型计时（除非你的电脑相当高端或者你不在乎离开电脑一会儿，否则这通常不是一个好主意），我们会看到平均训练时间有显著差异——从随机森林的惊人的 5 分钟到逻辑回归的 4 秒。（当然，这些绝对数字取决于机器。）

rfc = RandomForestClassifier(n_estimators=100, max_depth=2, 
  n_jobs=-1, random_state=0)
efs = ExhaustiveFeatureSelector(rfc, max_features=5,
  min_features=1, scoring='accuracy', 
  print_progress=True, cv=5)
%timeit efs.fit(X_train_enc, y_train.values.ravel())
5min 8s ± 3 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
lr = LogisticRegression(solver='liblinear')
efs = ExhaustiveFeatureSelector(lr, max_features=5,
  min_features=1, scoring='accuracy', 
  print_progress=True, cv=5)
%timeit efs.fit(X_train_enc, y_train.values.ravel())
4.29 s ± 45.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

如我所述，全面特征选择可以提供关于要选择哪些特征的非常清晰的指导，但这可能对许多项目来说代价太高。实际上，它可能更适合于诊断工作而不是用于机器学习管道。如果一个线性模型是合适的，它可以显著降低计算成本。

前向、后向和全面特征选择等包装方法会消耗系统资源，因为它们需要每次迭代时都进行训练，而选择的算法越难实现，这个问题就越严重。递归特征消除（RFE）在过滤方法的简单性和包装方法提供的信息之间是一种折衷。它与后向特征选择类似，但它在每次迭代中通过基于模型的整体性能而不是重新评估每个特征来简化特征的移除。我们将在下一节中探讨递归特征选择。

在回归模型中递归消除特征

一个流行的包装方法是 RFE。这种方法从所有特征开始，移除权重最低的一个（基于系数或特征重要性度量），然后重复此过程，直到确定最佳拟合模型。当移除一个特征时，它会得到一个反映其移除点的排名。

RFE 可以用于回归模型和分类模型。我们将从在回归模型中使用它开始：

我们导入必要的库，其中三个我们尚未使用：来自 scikit-learn 的 RFE、RandomForestRegressor 和 LinearRegression 模块：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression

接下来，我们加载工资的 NLS 数据并创建训练和测试 DataFrame：

nls97wages = pd.read_csv("data/nls97wages.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','motherhighgrade',
  'fatherhighgrade','parentincome','gender','completedba']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97wages[feature_cols],\
  nls97wages[['weeklywage']], test_size=0.3, random_state=0)

我们需要编码 gender 特征并标准化其他特征以及目标（wageincome）。我们不编码或缩放二进制特征 completedba：

ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = feature_cols[:-2]
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Male','completedba']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Male','completedba']])
scaler.fit(y_train)
y_train, y_test = \
  pd.DataFrame(scaler.transform(y_train),
  columns=['weeklywage'], index=y_train.index),\
  pd.DataFrame(scaler.transform(y_test),
  columns=['weeklywage'], index=y_test.index)

现在，我们准备进行一些递归特征选择。由于 RFE 是一种包装方法，我们需要选择一个算法，该算法将围绕选择进行包装。在这种情况下，回归的随机森林是有意义的。我们正在模拟一个连续的目标，并且不希望假设特征和目标之间存在线性关系。

使用 scikit-learn 实现 RFE 比较简单。我们实例化一个 RFE 对象，在过程中指定我们想要的估计器。我们指示 RandomForestRegressor。然后我们拟合模型并使用 get_support 获取选定的特征。我们将 max_depth 限制为 2 以避免过拟合：

rfr = RandomForestRegressor(max_depth=2)
treesel = RFE(estimator=rfr, n_features_to_select=5)
treesel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[treesel.get_support()]
selcols
 Index(['satmath', 'gpaoverall', 'parentincome', 'gender_Male', 'completedba'], dtype='object')

注意，这与使用带有 F 检验的滤波方法（针对工资收入目标）得到的特征列表略有不同。在这里选择了 gpaoverall 和 motherhighgrade，而不是 gender 标志或 gpascience。

我们可以使用 ranking_ 属性来查看每个被消除的特征何时被移除：

pd.DataFrame({'ranking': treesel.ranking_,
  'feature': X_train_enc.columns},
   columns=['feature','ranking']).\
   sort_values(['ranking'], ascending=True)
           feature                ranking
1          satmath                1
5          gpaoverall             1
8          parentincome           1
9          gender_Male            1
10         completedba            1
6          motherhighgrade        2
2          gpascience             3
0          satverbal              4
3          gpaenglish             5
4          gpamath                6
7          fatherhighgrade        7

在第一次交互后移除了 fatherhighgrade，在第二次交互后移除了 gpamath。

让我们运行一些测试统计量。我们仅在随机森林回归器模型上拟合选定的特征。RFE 选择器的 transform 方法给我们的是 treesel.transform(X_train_enc) 中选定的特征。我们可以使用 score 方法来获取 r 平方值，也称为确定系数。r 平方是我们模型解释的总变异百分比的度量。我们得到了一个非常低的分数，表明我们的模型只解释了很少的变异。（请注意，这是一个随机过程，所以我们每次拟合模型时可能会得到不同的结果。）
```
rfr.fit(treesel.transform(X_train_enc), y_train.values.ravel())
rfr.score(treesel.transform(X_test_enc), y_test)
0.13612629794428466
```

让我们看看使用带有线性回归模型的 RFE 是否能得到更好的结果。此模型返回与随机森林回归器相同的特征：

lr = LinearRegression()
lrsel = RFE(estimator=lr, n_features_to_select=5)
lrsel.fit(X_train_enc, y_train)
selcols = X_train_enc.columns[lrsel.get_support()]
selcols
Index(['satmath', 'gpaoverall', 'parentincome', 'gender_Male', 'completedba'], dtype='object')

让我们评估线性模型：

lr.fit(lrsel.transform(X_train_enc), y_train)
lr.score(lrsel.transform(X_test_enc), y_test)
0.17773742846314056

线性模型实际上并不比随机森林模型好多少。这可能是这样一个迹象，即我们可用的特征总体上只捕捉到每周工资变异的一小部分。这是一个重要的提醒，即我们可以识别出几个显著的特征，但仍然有一个解释力有限的模型。（也许这也是一个好消息，即我们的标准化测试分数，甚至我们的学位获得，虽然重要但不是多年后我们工资的决定性因素。）

让我们尝试使用分类模型进行 RFE。

在分类模型中递归消除特征

RFE 也可以是分类问题的一个很好的选择。我们可以使用 RFE 来选择完成学士学位模型的特征。你可能还记得，我们在本章前面使用穷举特征选择来选择该模型的特征。让我们看看使用 RFE 是否能获得更高的准确率或更容易训练的模型：

我们导入本章迄今为止一直在使用的相同库：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFE
from sklearn.metrics import accuracy_score

接下来，我们从 NLS 教育成就数据中创建训练和测试数据：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, 
  random_state=0)

然后，我们编码和缩放训练和测试数据：

ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

我们实例化一个随机森林分类器并将其传递给 RFE 选择方法。然后我们可以拟合模型并获取所选特征。

rfc = RandomForestClassifier(n_estimators=100, max_depth=2, 
  n_jobs=-1, random_state=0)
treesel = RFE(estimator=rfc, n_features_to_select=5)
treesel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[treesel.get_support()]
selcols
Index(['satverbal', 'satmath', 'gpascience', 'gpaenglish', 'gpaoverall'], dtype='object')

我们还可以使用 RFE 的ranking_属性来展示特征的排名：

pd.DataFrame({'ranking': treesel.ranking_,
  'feature': X_train_enc.columns},
   columns=['feature','ranking']).\
   sort_values(['ranking'], ascending=True)

           feature                 ranking
0          satverbal               1
1          satmath                 1
2          gpascience              1
3          gpaenglish              1
5          gpaoverall              1
4          gpamath                 2
8          parentincome            3
7          fatherhighgrade         4
6          motherhighgrade         5
9          gender_Female           6

让我们看看使用与我们的基线模型相同的随机森林分类器，使用所选特征的模型的准确率：

rfc.fit(treesel.transform(X_train_enc), y_train.values.ravel())
y_pred = rfc.predict(treesel.transform(X_test_enc))
accuracy_score(y_test, y_pred)
0.684981684981685

回想一下，我们使用穷举特征选择获得了 67%的准确率。这里我们得到的准确率大致相同。然而，RFE 的好处是它比穷举特征选择更容易训练。

包装和类似包装特征选择方法中的另一种选择是scikit-learn集成方法。我们将在下一节中使用scikit-learn的随机森林分类器来使用它。

使用 Boruta 进行特征选择

Boruta 包在特征选择方面采用独特的方法，尽管它与包装方法有一些相似之处。对于每个特征，Boruta 创建一个影子特征，它与原始特征具有相同的值范围，但具有打乱后的值。然后它评估原始特征是否比影子特征提供更多信息，逐渐移除提供最少信息的特征。Boruta 在每个迭代中输出已确认、尝试和拒绝的特征。

让我们使用 Boruta 包来选择完成学士学位分类模型的特征（如果你还没有安装 Boruta 包，可以使用pip安装）：

我们首先加载必要的库：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from boruta import BorutaPy
from sklearn.metrics import accuracy_score

我们再次加载 NLS 教育成就数据并创建训练和测试 DataFrame：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)

接下来，我们对训练和测试数据进行编码和缩放：

ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

我们以与运行 RFE 特征选择相同的方式运行 Boruta 特征选择。我们再次使用随机森林作为基线方法。我们实例化一个随机森林分类器并将其传递给 Boruta 的特征选择器。然后我们拟合模型，该模型在100次迭代后停止，识别出提供信息的9个特征：

rfc = RandomForestClassifier(n_estimators=100, 
  max_depth=2, n_jobs=-1, random_state=0)
borsel = BorutaPy(rfc, random_state=0, verbose=2)
borsel.fit(X_train_enc.values, y_train.values.ravel())
BorutaPy finished running.
Iteration:            100 / 100
Confirmed:            9
Tentative:            1
Rejected:             0
selcols = X_train_enc.columns[borsel.support_]
selcols
Index(['satverbal', 'satmath', 'gpascience', 'gpaenglish', 'gpamath', 'gpaoverall', 'motherhighgrade', 'fatherhighgrade', 'parentincome', 'gender_Female'], dtype='object')

我们可以使用ranking_属性来查看特征的排名：

pd.DataFrame({'ranking': borsel.ranking_,
  'feature': X_train_enc.columns},
   columns=['feature','ranking']).\
   sort_values(['ranking'], ascending=True)
           feature               ranking
0          satverbal             1
1          satmath               1
2          gpascience            1
3          gpaenglish            1
4          gpamath               1
5          gpaoverall            1
6          motherhighgrade       1
7          fatherhighgrade       1
8          parentincome          1
9          gender_Female         2

为了评估模型的准确率，我们仅使用所选特征来拟合随机森林分类器模型。然后我们可以对测试数据进行预测并计算准确率：

rfc.fit(borsel.transform(X_train_enc.values), y_train.values.ravel())
y_pred = rfc.predict(borsel.transform(X_test_enc.values))
accuracy_score(y_test, y_pred)
0.684981684981685

Boruta 的吸引力之一在于其对每个特征选择的说服力。如果一个特征被选中，那么它很可能提供了信息，这些信息不是通过排除它的特征组合所捕获的。然而，它在计算上相当昂贵，与穷举特征选择不相上下。它可以帮助我们区分哪些特征是重要的，但可能并不总是适合那些训练速度很重要的流水线。

最后几节展示了包装特征选择方法的某些优点和缺点。在下一节中，我们将探讨嵌入式选择方法。这些方法比过滤器方法提供更多信息，但又不具备包装方法的计算成本。它们通过将特征选择嵌入到训练过程中来实现这一点。我们将使用我们迄今为止所使用的数据来探讨嵌入式方法。

使用正则化和其他嵌入式方法

正则化方法是嵌入式方法。与包装方法一样，嵌入式方法根据给定的算法评估特征。但它们的计算成本并不高。这是因为特征选择已经嵌入到算法中，所以随着模型的训练而发生。

嵌入式模型使用以下过程：

训练一个模型。
估计每个特征对模型预测的重要性。
移除重要性低的特征。

正则化通过向任何模型添加惩罚来约束参数来实现这一点。L1 正则化，也称为lasso 正则化，将回归模型中的某些系数缩小到 0，从而有效地消除了这些特征。

使用 L1 正则化

我们将使用 L1 正则化和逻辑回归来选择学士学位达成模型的特征：我们需要首先导入所需的库，包括我们将首次使用的模块，scikit-learn中的SelectFromModel：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectFromModel
from sklearn.metrics import accuracy_score

接下来，我们加载关于教育成就的 NLS 数据：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade','fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, 
  random_state=0)

然后，我们对训练数据和测试数据进行编码和缩放：

ohe = OneHotEncoder(drop_last=True, 
                    variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

现在，我们准备根据逻辑回归和 L1 惩罚进行特征选择：

lr = LogisticRegression(C=1, penalty="l1", 
                        solver='liblinear')
regsel = SelectFromModel(lr, max_features=5)
regsel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[regsel.get_support()]
selcols
Index(['satmath', 'gpascience', 'gpaoverall', 
'fatherhighgrade', 'gender_Female'], dtype='object')

让我们来评估模型的准确性。我们得到了一个准确率分数为0.68：

lr.fit(regsel.transform(X_train_enc), 
       y_train.values.ravel())
y_pred = lr.predict(regsel.transform(X_test_enc))
accuracy_score(y_test, y_pred)
0.684981684981685

这给我们带来了与学士学位完成的前向特征选择相当相似的结果。在那个例子中，我们使用随机森林分类器作为包装方法。

在这种情况下，Lasso 正则化是特征选择的一个好选择，尤其是当性能是一个关键关注点时。然而，它确实假设特征与目标之间存在线性关系，这可能并不合适。幸运的是，有一些嵌入式特征选择方法不做出这种假设。对于嵌入式模型来说，逻辑回归的一个好替代品是随机森林分类器。我们将使用相同的数据尝试这种方法。

使用随机森林分类器

在本节中，我们将使用随机森林分类器：

我们可以使用SelectFromModel来使用随机森林分类器而不是逻辑回归：

rfc = RandomForestClassifier(n_estimators=100, 
  max_depth=2, n_jobs=-1, random_state=0)
rfcsel = SelectFromModel(rfc, max_features=5)
rfcsel.fit(X_train_enc, y_train.values.ravel())
selcols = X_train_enc.columns[rfcsel.get_support()]
selcols
Index(['satverbal', 'gpascience', 'gpaenglish', 
  'gpaoverall'], dtype='object')

这实际上选择与 lasso 回归非常不同的特征。satmath、fatherhighgrade和gender_Female不再被选中，而satverbal和gpaenglish被选中。这很可能部分是由于线性假设的放宽。

让我们评估随机森林分类器模型的准确性。我们得到了0.67的准确率。这几乎与我们在 lasso 回归中得到的分数相同：

rfc.fit(rfcsel.transform(X_train_enc), 
        y_train.values.ravel())
y_pred = rfc.predict(rfcsel.transform(X_test_enc))
accuracy_score(y_test, y_pred)
0.673992673992674

嵌入式方法通常比包装方法 CPU-/GPU 密集度低，但仍然可以产生良好的结果。在本节的学士学位完成模型中，我们得到了与基于穷举特征选择模型相同的准确率。

我们之前讨论的每种方法都有重要的应用场景，正如我们所讨论的。然而，我们还没有真正讨论一个非常具有挑战性的特征选择问题。如果你简单地有太多的特征，其中许多特征在你的模型中独立地解释了某些重要内容，你会怎么做？在这里，“太多”意味着有如此多的特征，以至于模型无法高效地运行，无论是训练还是预测目标值。我们如何在不牺牲模型部分预测能力的情况下减少特征集？在这种情况下，主成分分析（PCA）可能是一个好的方法。我们将在下一节中讨论 PCA。

使用主成分分析

PCA 是一种与之前讨论的任何方法都截然不同的特征选择方法。PCA 允许我们用有限数量的组件替换现有的特征集，每个组件都解释了重要数量的方差。它是通过找到一个捕获最大方差量的组件，然后是一个捕获剩余最大方差量的第二个组件，然后是一个第三个组件，以此类推来做到这一点的。这种方法的一个关键优势是，这些被称为主成分的组件是不相关的。我们在第十五章，主成分分析中详细讨论 PCA。

虽然我在这里将主成分分析（PCA）视为一种特征选择方法，但可能更合适将其视为降维工具。当我们需要限制维度数量而又不希望牺牲太多解释力时，我们使用它来进行特征选择。

让我们再次使用 NLS 数据，并使用 PCA 为学士学位完成模型选择特征：

我们首先加载必要的库。在本章中，我们还没有使用过的模块是scikit-learn的PCA：

import pandas as pd
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

接下来，我们再次创建训练和测试 DataFrame：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpascience',
  'gpaenglish','gpamath','gpaoverall','gender',
  'motherhighgrade', 'fatherhighgrade','parentincome']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3,
  random_state=0)

我们需要对数据进行缩放和编码。在 PCA 中，缩放尤其重要：

ohe = OneHotEncoder(drop_last=True, 
                    variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

现在，我们实例化一个PCA对象并拟合模型：
```
pca = PCA(n_components=5)
pca.fit(X_train_enc)
```
PCA对象的components_属性返回了所有 10 个特征在每个 5 个成分上的得分。对第一个成分贡献最大的特征是得分绝对值最高的那些，在这种情况下，是gpaoverall、gpaenglish和gpascience。对于第二个成分，最重要的特征是motherhighgrade、fatherhighgrade和parentincome。satverbal和satmath驱动第三个成分。

在以下输出中，列0到4是五个主成分：

pd.DataFrame(pca.components_,
  columns=X_train_enc.columns).T
                   0       1      2       3       4
satverbal         -0.34   -0.16  -0.61   -0.02   -0.19
satmath           -0.37   -0.13  -0.56    0.10    0.11
gpascience        -0.40    0.21   0.18    0.03    0.02
gpaenglish        -0.40    0.22   0.18    0.08   -0.19
gpamath           -0.38    0.24   0.12    0.08    0.23
gpaoverall        -0.43    0.25   0.23   -0.04   -0.03
motherhighgrade   -0.19   -0.51   0.24   -0.43   -0.59
fatherhighgrade   -0.20   -0.51   0.18   -0.35    0.70
parentincome      -0.16   -0.46   0.28    0.82   -0.08
gender_Female     -0.02    0.08   0.12   -0.04   -0.11

另一种理解这些得分的方式是，它们表明每个特征对成分的贡献程度。（实际上，如果对每个成分，你将 10 个得分平方然后求和，你会得到一个总和为 1。）

让我们也检查每个成分解释了特征中多少方差。第一个成分单独解释了 46%的方差，第二个成分额外解释了 19%。我们可以使用 NumPy 的cumsum方法来查看五个成分累积解释了多少特征方差。我们可以用 5 个成分解释 10 个特征中的 87%的方差：
```
pca.explained_variance_ratio_
array([0.46073387, 0.19036089, 0.09295703, 0.07163009, 0.05328056])
np.cumsum(pca.explained_variance_ratio_)
array([0.46073387, 0.65109476, 0.74405179, 0.81568188, 0.86896244])
```

让我们根据这五个主成分来转换测试数据中的特征。这返回了一个只包含五个主成分的 NumPy 数组。我们查看前几行。我们还需要转换测试 DataFrame：

X_train_pca = pca.transform(X_train_enc)
X_train_pca.shape
(634, 5)
np.round(X_train_pca[0:6],2)
array([[ 2.79, -0.34,  0.41,  1.42, -0.11],
       [-1.29,  0.79,  1.79, -0.49, -0.01],
       [-1.04, -0.72, -0.62, -0.91,  0.27],
       [-0.22, -0.8 , -0.83, -0.75,  0.59],
       [ 0.11, -0.56,  1.4 ,  0.2 , -0.71],
       [ 0.93,  0.42, -0.68, -0.45, -0.89]])
X_test_pca = pca.transform(X_test_enc)

现在，我们可以使用这些主成分来拟合一个关于学士学位完成情况的模型。让我们运行一个随机森林分类。

我们首先创建一个随机森林分类器对象。然后，我们将带有主成分和目标值的训练数据传递给其fit方法。我们将带有成分的测试数据传递给分类器的predict方法，然后得到一个准确度分数：
```
rfc = RandomForestClassifier(n_estimators=100, 
  max_depth=2, n_jobs=-1, random_state=0)
rfc.fit(X_train_pca, y_train.values.ravel())
y_pred = rfc.predict(X_test_pca)
accuracy_score(y_test, y_pred)
0.7032967032967034
```

当特征选择挑战是我们有高度相关的特征，并且我们希望在不过度减少解释方差的情况下减少维度数量时，PCA 等降维技术可以是一个好的选择。在这个例子中，高中 GPA 特征一起移动，父母的教育水平和收入水平以及 SAT 特征也是如此。它们成为了我们前三个成分的关键特征。（可以认为我们的模型只需要那三个成分，因为它们共同解释了特征变异的 74%。）

根据你的数据和建模目标，PCA（主成分分析）有几种修改方式可能是有用的。这包括处理异常值和正则化的策略。通过使用核函数，PCA 还可以扩展到那些成分不能线性分离的情况。我们将在第十五章**，《主成分分析》*中详细讨论 PCA。

让我们总结一下本章所学的内容。

摘要

在本章中，我们讨论了从过滤方法到包装方法再到嵌入式方法的一系列特征选择方法。我们还看到了它们如何与分类和连续目标一起工作。对于包装和嵌入式方法，我们考虑了它们如何与不同的算法一起工作。

过滤方法运行和解释都非常简单，且对系统资源的影响较小。然而，它们在评估每个特征时并没有考虑其他特征。而且，它们也没有告诉我们这种评估可能会因所使用的算法而有所不同。包装方法没有这些限制，但计算成本较高。嵌入式方法通常是一个很好的折衷方案，它们根据多元关系和给定的算法选择特征，而不像包装方法那样对系统资源造成过多负担。我们还探讨了如何通过降维方法 PCA 来改进我们的特征选择。

你可能也注意到了，我在本章中稍微提到了一点模型验证。我们将在下一章更详细地介绍模型验证。

第六章：第六章: 准备模型评估

在开始运行模型之前，思考如何评估模型性能是一个好主意。一种常见的技术是将数据分为训练集和测试集。我们在早期阶段就做这件事，以避免所谓的数据泄露；也就是说，基于原本打算用于模型评估的数据进行分析。在本章中，我们将探讨创建训练集的方法，包括如何确保训练数据具有代表性。我们还将探讨交叉验证策略，如K 折，它解决了使用静态训练/测试分割的一些局限性。我们还将开始更仔细地评估模型性能。

你可能会想知道为什么我们在详细讨论任何算法之前讨论模型评估。这是因为有一个实际考虑。我们倾向于在具有相似目的的算法中使用相同的指标和评估技术。在评估分类模型时，我们检查准确率和敏感性，在检查回归模型时，我们检查平均绝对误差和 R 平方。我们对所有监督学习模型进行交叉验证。因此，我们将在以下章节中多次重复介绍这些策略。你甚至可能会在概念稍后重新引入时回到这些页面。

除了这些实际考虑因素之外，当我们不将数据提取、数据清洗、探索性分析、特征工程和预处理、模型指定和模型评估视为离散的、顺序的任务时，我们的建模工作会得到改善。如果你只构建了 6 个月的机器学习模型，或者超过 30 年，你可能会欣赏到这种严格的顺序与数据科学家的工作流程不一致。我们总是在准备模型验证，并且总是在清理数据。这是好事。当我们整合这些任务时，我们会做得更好；当我们选择特征时继续审查我们的数据清洗，以及当我们计算精确度或均方根误差后回顾双变量相关或散点图时。

我们还将花费相当多的时间构建这些概念的可视化。在处理分类问题时，养成查看混淆矩阵和累积准确率轮廓的习惯是一个好主意，而在处理连续目标时，则查看残差图。这同样会在后续章节中对我们大有裨益。

具体来说，在本章中，我们将涵盖以下主题：

测量二分类的准确率、敏感性、特异性和精确度
检查二分类的 CAP、ROC 和精确度-敏感性曲线
评估多分类模型
评估回归模型
使用 K 折交叉验证
使用管道预处理数据

技术要求

在本章中，我们将使用feature_engine和matplotlib库，以及 scikit-learn 库。您可以使用pip安装这些包。本章的代码文件可以在本书的 GitHub 仓库中找到，网址为github.com/PacktPublishing/Data-Cleaning-and-Exploration-with-Machine-Learning。

测量二元分类的准确率、灵敏度、特异度和精确度

在评估分类模型时，我们通常想知道我们正确的情况有多频繁。在二元目标的情况下——目标有两个可能的分类值——我们计算准确率为预测正确分类的次数与观察总数之比。

但是，根据分类问题，准确率可能不是最重要的性能指标。也许我们愿意接受更多的假阳性，以换取能够识别更多真正正面的模型，即使这意味着较低的准确率。这可能适用于预测患有乳腺癌、安全漏洞或桥梁结构损坏的可能性模型。在这些情况下，我们可能更强调灵敏度（识别正案例的倾向）而不是准确率。

另一方面，我们可能希望有一个模型能够以高可靠性识别出负面案例，即使这意味着它不能很好地识别正面案例。特异度是模型识别出的所有负面案例的百分比。

精确度，即预测为正的预测值实际上是正的百分比，是另一个重要的度量。对于某些应用，限制假阳性可能很重要，即使这意味着我们必须容忍较低的灵敏度。一个苹果种植者，使用图像识别来识别坏苹果，可能更倾向于选择高精确度的模型，而不是更灵敏的模型，不希望不必要地丢弃苹果。

这可以通过查看混淆矩阵来更清楚地说明：

图 6.1 – 混淆矩阵

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

灵敏度是指我们正确预测正面的次数除以正面的总数。回顾一下混淆矩阵，确认实际的正值可以是预测正值（TP）或预测负值（FN）。灵敏度也被称为召回率或真正率：

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

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

当存在类别不平衡时，如准确率和灵敏度这样的指标可能会给我们关于模型性能的非常不同的估计。一个极端的例子将说明这一点。黑猩猩有时会尝试白蚁捕捞，把一根棍子插进蚁丘里，希望能捕获几只白蚁。这只有偶尔会成功。我不是灵长类学家，但我们可以将成功的捕捞尝试建模为使用棍子的大小、年份和时间以及黑猩猩年龄的函数。在我们的测试数据中，捕捞尝试只有 2%的时间会成功。（这些数据是为了演示而编造的。）

让我们再假设我们构建了一个成功白蚁捕捞的分类模型，其灵敏度为 50%。因此，如果我们测试数据中有 100 次捕捞尝试，我们只会正确预测其中一次成功的尝试。还有一个假阳性，即我们的模型在捕捞失败时预测了成功的捕捞。这给我们以下混淆矩阵：

图 6.2 – 成功白蚁捕捞的混淆矩阵

注意，我们得到了一个非常高的准确率 98%——即（97+1）/ 100。我们得到了高准确率和低灵敏度，因为大部分捕捞尝试都是负的，这很容易预测。一个总是预测失败的模型也会有一个 98%的准确率。

现在，让我们用真实数据来查看这些模型评估指标。我们可以通过实验一个k 最近邻（KNN）模型来预测学士学位的获得，并评估其准确率、灵敏度、特异性和精确率：

我们将首先加载用于编码和标准化数据的库，以及用于创建训练和测试数据框（DataFrames）的库。我们还将加载 scikit-learn 的 KNN 分类器和metrics库：

import pandas as pd
import numpy as np
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
import sklearn.metrics as skmet
import matplotlib.pyplot as plt

现在，我们可以创建训练和测试数据框，并对数据进行编码和缩放：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpaoverall',
  'parentincome','gender']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

让我们创建一个 KNN 分类模型。我们不会太在意如何指定它，因为我们只想关注本节中的评估指标。我们将使用feature_cols中列出的所有特征。我们使用 KNN 分类器的预测方法从测试数据中生成预测：
```
knn = KNeighborsClassifier(n_neighbors = 5)
knn.fit(X_train_enc, y_train.values.ravel())
pred = knn.predict(X_test_enc)
```

我们可以使用 scikit-learn 绘制混淆矩阵。我们将传递测试数据中的实际值（y_test）和预测值到confusion_matrix方法：

cm = skmet.confusion_matrix(y_test, pred, labels=knn.classes_)
cmplot = skmet.ConfusionMatrixDisplay(
  confusion_matrix=cm, 
  display_labels=['Negative', 'Positive'])
cmplot.plot()
cmplot.ax_.set(title='Confusion Matrix', 
  xlabel='Predicted Value', ylabel='Actual Value')

这生成了以下图表：

图 6.3 – 实际值和预测值的混淆矩阵

我们也可以只返回真正的负值、假阳性、假阴性和真正阳性的计数：

tn, fp, fn, tp = skmet.confusion_matrix(
  y_test.values.ravel(), pred).ravel()
tn, fp, fn, tp
(53, 63, 31, 126)

我们现在有了计算准确率、灵敏度、特异性和精确率所需的所有信息：

accuracy = (tp + tn) / pred.shape[0]
accuracy
0.6556776556776557
sensitivity = tp / (tp + fn)
sensitivity
0.802547770700637
specificity = tn / (tn+fp)
specificity
0.45689655172413796
precision = tp / (tp + fp)
precision
0.6666666666666666

这个模型相对精度较低，但灵敏度略好；也就是说，它更好地识别了测试数据中完成学士学位的人，而不是正确识别所有学位完成者和未完成者的整体情况。如果我们回顾混淆矩阵，我们会看到有相当数量的假阳性，因为我们的模型预测测试数据中有 63 个人将会有学士学位，而实际上并没有。

我们还可以使用 scikit-learn 提供的便捷方法直接生成这些统计数据：

skmet.accuracy_score(y_test.values.ravel(), pred)
0.6556776556776557
skmet.recall_score(y_test.values.ravel(), pred)
0.802547770700637
skmet.precision_score(y_test.values.ravel(), pred)
0.6666666666666666

仅为了比较，让我们尝试使用随机森林分类器，看看是否能得到更好的结果。

让我们将随机森林分类器拟合到相同的数据，并再次调用confusion_matrix：

rfc = RandomForestClassifier(n_estimators=100, 
  max_depth=2, n_jobs=-1, random_state=0)
rfc.fit(X_train_enc, y_train.values.ravel())
pred = rfc.predict(X_test_enc)
tn, fp, fn, tp = skmet.confusion_matrix(
  y_test.values.ravel(), pred).ravel()
tn, fp, fn, tp
(49, 67, 17, 140)
accuracy = (tp + tn) / pred.shape[0]
accuracy
0.6923076923076923
sensitivity = tp / (tp + fn)
sensitivity
0.89171974522293
specificity = tn / (tn+fp)
specificity
0.4224137931034483
precision = tp / (tp + fp)
precision
0.6763285024154589

第二个模型比第一个模型显著减少了假阴性，并增加了真阳性。它不太可能预测测试数据中的个体没有完成学士学位，而当个人完成了学士学位时，更有可能预测这个人有学士学位。较低的 FP 和较高的 TP 的主要影响是灵敏度显著提高。第二个模型有 89%的时间识别出实际正值，而第一个模型只有 80%。

我们在本节中讨论的措施——准确性、灵敏度、特异性和精确度——在评估分类模型时都值得一看。但是，例如，在精度和灵敏度之间，我们有时可能会面临难以权衡的情况。数据科学家在构建分类模型时，会依赖几种标准的可视化方法来提高我们对这些权衡的认识。我们将在下一节中探讨这些可视化方法。

检查二元分类的 CAP、ROC 和精度-灵敏度曲线

可视化二元分类模型性能的方法有很多。一种相对直接的可视化方法是累积准确率曲线（CAP），它显示了我们的模型识别正类（或积极案例）的能力。它显示了 X 轴上的累积案例和 Y 轴上的累积积极结果。CAP 曲线是了解我们的模型在区分正类观察方面做得如何的好方法。（在讨论二元分类模型时，我将交替使用正类和积极这两个术语。）

接收者操作特征（ROC）曲线说明了在调整分类正值的阈值时，模型灵敏度（能够识别正值）与假阳性率之间的权衡。同样，精度-灵敏度曲线显示了在调整阈值时，我们积极预测的可靠性（它们的精度）与灵敏度（我们的模型识别实际正值的能）力之间的关系。

构建 CAP 曲线

让我们从学士学位完成 KNN 模型的 CAP 曲线开始。让我们也将其与决策树模型进行比较。同样，我们在这里不会进行太多的特征选择。上一章详细介绍了特征选择。

除了我们模型的曲线外，CAP 曲线还有用于比较的随机模型和完美模型的图表。随机模型除了提供正值的整体分布信息外，没有其他信息。完美模型精确地预测正值。为了说明这些图表是如何绘制的，我们将从一个假设的例子开始。想象一下，你从一副洗好的牌中抽取前六张牌。你创建一个表格，其中一列是累积牌总数，下一列是红牌的数量。它可能看起来像这样：

图 6.4 – 玩牌样本

我们可以根据我们对红牌数量的了解绘制一个随机模型。随机模型只有两个点，(0,0)和(6,3)，但这就足够了。

完美模型图表需要更多的解释。如果我们的模型完美预测红牌，并且按预测降序排列，我们会得到图 6.5。累积 in-class count 与牌的数量相匹配，直到红牌耗尽，在这个例子中是 3 张。使用完美模型的累积 in-class total 图表将有两个斜率；在达到 in-class total 之前等于 1，之后为 0：

图 6.5 – 玩牌样本

现在我们已经足够了解如何绘制随机模型和完美模型。完美模型将有三点：(0,0)，(in-class count, in-class count)，和(number of cards, in-class count)。在这种情况下，in-class count 是3，卡片数量是6：

numobs = 6

inclasscnt = 3

plt.yticks([1,2,3])

plt.plot([0, numobs], [0, inclasscnt], c = 'b', label = 'Random Model')

plt.plot([0, inclasscnt, numobs], [0, inclasscnt, inclasscnt], c = 'grey', linewidth = 2, label = 'Perfect Model')

plt.title("Cumulative Accuracy Profile")

plt.xlabel("Total Cards")

plt.ylabel("In-class (Red) Cards")

这会产生以下图表：

图 6.6 – 使用玩牌数据的 CAP

要理解完美模型相对于随机模型的改进，可以考虑随机模型在中间点预测多少红牌 – 那就是说，3 张牌。在那个点上，随机模型会预测 1.5 张红牌。然而，完美模型会预测 3 张。(记住，我们已经按预测顺序降序排列了牌。)

使用虚构数据构建了随机模型和完美模型的图表后，让我们用我们的学士学位完成数据试试：

首先，我们必须导入与上一节相同的模块：

import pandas as pd
import numpy as np
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
import sklearn.metrics as skmet
import matplotlib.pyplot as plt
import seaborn as sb

然后，我们加载、编码和缩放 NLS 学士学位数据：

nls97compba = pd.read_csv("data/nls97compba.csv")
feature_cols = ['satverbal','satmath','gpaoverall',
  'parentincome','gender']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97compba[feature_cols],\
  nls97compba[['completedba']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

接下来，我们创建KNeighborsClassifier和RandomForestClassifier实例：

knn = KNeighborsClassifier(n_neighbors = 5)
rfc = RandomForestClassifier(n_estimators=100, max_depth=2, 
  n_jobs=-1, random_state=0)

我们现在可以开始绘制我们的 CAP 曲线了。我们将首先绘制一个随机模型，然后是一个完美模型。这些模型不使用任何信息（除了正值的整体分布）并且提供完美信息。

我们计算测试数据中的观测数和正值的数量。我们将使用(0,0)和(观测数，类内计数)来绘制随机模型线。对于完美模型，我们将从(0,0)到(类内计数，类内计数)绘制一条线，因为该模型可以完美地区分类内值（它永远不会出错）。在该点右侧是平的，因为再也没有更多的正值可以找到。

我们还将绘制一条垂直线在中间，以及与随机模型线相交的水平线。这将在以后更有用：

numobs = y_test.shape[0]
inclasscnt = y_test.iloc[:,0].sum()
plt.plot([0, numobs], [0, inclasscnt], c = 'b', label = 'Random Model')
plt.plot([0, inclasscnt, numobs], [0, inclasscnt, inclasscnt], c = 'grey', linewidth = 2, label = 'Perfect Model')
plt.axvline(numobs/2, color='black', linestyle='dashed', linewidth=1)
plt.axhline(numobs/2, color='black', linestyle='dashed', linewidth=1)
plt.title("Cumulative Accuracy Profile")
plt.xlabel("Total Observations")
plt.ylabel("In-class Observations")
plt.legend()

这产生了以下图表：

图 6.7 – 仅使用随机和完美模型的 CAP

接下来，我们定义一个函数来绘制我们传递给它的模型的 CAP 曲线。我们将使用predict_proba方法来获取一个数组，该数组包含测试数据中每个观测值在类内（在这种情况下，已完成学士学位）的概率。然后，我们将创建一个包含这些概率和实际目标值的 DataFrame，按概率降序排序，并计算正实际目标值的累计总和。

我们还将得到中间观测值的累计值，并在该点绘制一条水平线。最后，我们将绘制一条线，其 x 值为从 0 到观测数的数组，y 值为累计的类内总数：

def addplot(model, X, Xtest, y, modelname, linecolor):
  model.fit(X, y.values.ravel())
  probs = model.predict_proba(Xtest)[:, 1]

  probdf = pd.DataFrame(zip(probs, y_test.values.ravel()),
    columns=(['prob','inclass']))
  probdf.loc[-1] = [0,0]
  probdf = probdf.sort_values(['prob','inclass'],
    ascending=False).\
    assign(inclasscum = lambda x: x.inclass.cumsum())
  inclassmidpoint = \
    probdf.iloc[int(probdf.shape[0]/2)].inclasscum
  plt.axhline(inclassmidpoint, color=linecolor,
    linestyle='dashed', linewidth=1)
  plt.plot(np.arange(0, probdf.shape[0]),
    probdf.inclasscum, c = linecolor,
    label = modelname, linewidth = 4)

现在，让我们使用相同的数据运行 KNN 和随机森林分类器模型的函数：

addplot(knn, X_train_enc, X_test_enc, y_train,
  'KNN', 'red')
addplot(rfc, X_train_enc, X_test_enc, y_train,
  'Random Forest', 'green')
plt.legend()

这更新了我们的早期图表：

图 6.8 – 使用 KNN 和随机森林模型的 CAP 更新

毫不奇怪，CAP 曲线显示我们的 KNN 和随机森林模型比随机猜测要好，但不如完美模型好。问题是，分别有多好和多差。水平线给我们一些想法。一个完美模型会在 138 个观测值中正确识别出 138 个正值。（回想一下，观测值是按最高可能性为正的顺序排序的。）随机模型会识别出 70 个（线未显示），而 KNN 和随机森林模型分别会识别出 102 和 103 个。我们的两个模型在区分正值方面与完美模型一样好，分别是 74%和 75%。在 70%到 80%之间被认为是好的模型；高于这个百分比的模型非常好，而低于这个百分比的模型较差。

绘制接收者操作特征（ROC）曲线

ROC 曲线说明了在调整阈值时，假阳性率与真阳性率（也称为灵敏度）之间的权衡。在进一步讨论之前，我们应该讨论假阳性率。它是模型错误地将实际负例（真负例加上假阳性）识别为正例的百分比：

在这里，你可以看到假阳性率与特异性之间的关系，这是在本章开头讨论过的。差异是分子。特异性是我们模型正确地将实际负例识别为负例的百分比：

我们还可以将假阳性率与灵敏度进行比较，灵敏度是实际正例（真阳性加上假阴性）的百分比，我们模型正确地将它们识别为正例：

我们通常面临灵敏度和假阳性率之间的权衡。我们希望我们的模型能够识别大量实际正例，但我们不希望假阳性率过高。什么是“过高”取决于你的上下文。

当区分负例和正例变得更加困难时，灵敏度和假阳性率之间的权衡就变得更加复杂。当我们绘制预测概率时，我们可以通过我们的学士学位完成模型看到这一点：

首先，让我们再次调整我们的随机森林分类器并生成预测和预测概率。我们会看到，当预测概率大于 0.500 时，该模型预测这个人将完成学士学位：

rfc.fit(X_train_enc, y_train.values.ravel())
pred = rfc.predict(X_test_enc)
pred_probs = rfc.predict_proba(X_test_enc)[:, 1]
probdf = pd.DataFrame(zip(
  pred_probs, pred, y_test.values.ravel()),
  columns=(['prob','pred','actual']))
probdf.groupby(['pred'])['prob'].agg(['min','max'])
                min             max
pred             
0.000           0.305           0.500
1.000           0.502           0.883

将这些概率的分布与实际类别值进行比较是有帮助的。我们可以使用密度图来完成这项工作：

sb.kdeplot(probdf.loc[probdf.actual==1].prob, 
  shade=True, color='red',
  label="Completed BA")
sb.kdeplot(probdf.loc[probdf.actual==0].prob,  
  shade=True, color='green',
  label="Did Not Complete")
plt.axvline(0.5, color='black', linestyle='dashed', linewidth=1)
plt.axvline(0.65, color='black', linestyle='dashed', linewidth=1)
plt.title("Predicted Probability Distribution")
plt.legend(loc="upper left")

这会产生以下图表：

图 6.9 – 在课内和课外观察的密度图

在这里，我们可以看到我们的模型在区分实际正例和负例方面有些困难，因为课内和课外重叠的部分相当多。阈值为 0.500（左侧虚线）会产生很多假阳性，因为课外观察分布（那些没有完成学士学位的人）中有相当一部分预测概率大于 0.500。如果我们提高阈值，比如到 0.650，我们会得到更多的假阴性，因为许多课内观察的概率低于 0.65。

基于测试数据和随机森林模型，构建 ROC 曲线很容易。roc_curve方法返回不同阈值（ths）下的假阳性率（fpr）和灵敏度（真阳性率，tpr）。

首先，让我们通过阈值绘制单独的假阳性率和灵敏度线：

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(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()

这会产生以下图表：

图 6.10 – 假阳性率和灵敏度线

在这里，我们可以看到提高阈值将提高（降低）我们的假阳性率，但也会降低我们的灵敏度。

现在，让我们绘制相关的 ROC 曲线，该曲线绘制了每个阈值下的假阳性率与灵敏度：

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')

这会产生以下图表：

图 6.11 – 带有假阳性率和灵敏度的 ROC 曲线

ROC 曲线表明，假阳性率和灵敏度之间的权衡在假阳性率约为 0.5 或更高时相当陡峭。让我们看看这对随机森林模型预测中使用的 0.5 阈值意味着什么。

让我们从阈值数组中选择一个接近 0.5 的索引，以及一个接近 0.4 和 0.6 的索引进行比较。然后，我们将为那些索引绘制垂直线表示假阳性率，以及为那些索引绘制水平线表示灵敏度值：

tholdind = np.where((ths>0.499) & (ths<0.501))[0][0]
tholdindlow = np.where((ths>0.397) & (ths<0.404))[0][0]
tholdindhigh = np.where((ths>0.599) & (ths<0.601))[0][0]
plt.vlines((fpr[tholdindlow],fpr[tholdind],
  fpr[tholdindhigh]), 0, 1, linestyles ="dashed", 
  colors =["green","blue","purple"])
plt.hlines((tpr[tholdindlow],tpr[tholdind],
  tpr[tholdindhigh]), 0, 1, linestyles ="dashed", 
  colors =["green","blue","purple"])

这更新了我们的图表：

图 6.12 – 带有阈值线的 ROC 曲线

这说明了在 0.5 阈值（蓝色虚线）下，预测中假阳性率和灵敏度之间的权衡。ROC 曲线在 0.5 以上的阈值处斜率非常小，例如 0.6 阈值（绿色虚线）。因此，将阈值从 0.6 降低到 0.5 会导致假阳性率显著降低（从超过 0.8 降低到低于 0.6），但灵敏度降低不多。然而，通过将阈值从 0.5 降低到 0.4（从蓝色到紫色线），假阳性率（降低）将导致灵敏度显著下降。它从近 90%下降到略高于 70%。

绘制精度-灵敏度曲线

当调整阈值时，检查精度和灵敏度之间的关系通常很有帮助。记住，精度告诉我们当我们预测一个正值时，我们正确的时间百分比：

我们可以通过提高将值分类为正的阈值来提高精度。然而，这可能会意味着灵敏度的降低。随着我们提高预测正值时正确的时间（精度），我们将减少我们能够识别的正值数量（灵敏度）。精度-灵敏度曲线，通常称为精度-召回曲线，说明了这种权衡。

在绘制精度-灵敏度曲线之前，让我们先看看分别针对阈值绘制的精度和灵敏度线：

我们可以使用precision_recall_curve方法获得精度-灵敏度曲线的点。我们移除最高阈值值的一些不规则性，这有时可能发生：

prec, sens, ths = skmet.precision_recall_curve(y_test, pred_probs)
prec = prec[1:-10]
sens = sens[1:-10]
ths  = ths[:-10]
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.set_xlim(0.3,0.9)
ax.legend()

这会产生以下图表：

图 6.13 – 精确度和灵敏度线

在这里，我们可以看到，当阈值高于 0.5 时，灵敏度下降更为陡峭。这种下降并没有在 0.6 阈值以上带来多少精确度的改进。

现在，让我们绘制灵敏度与精确度之间的关系，以查看精确度-灵敏度曲线：

fig, ax = plt.subplots()
ax.plot(sens, prec)
ax.set_title('Precision-Sensitivity Curve')
ax.set_xlabel('Sensitivity')
ax.set_ylabel('Precision')
plt.yticks(np.arange(0.2, 0.9, 0.2))

这生成了以下图表：

图 6.14 – 精确度-灵敏度曲线

精确度-灵敏度曲线反映了这样一个事实：对于这个特定的模型，灵敏度对阈值的反应比精确度更敏感。这意味着我们可以将阈值降低到 0.5 以下，以获得更高的灵敏度，而不会显著降低精确度。

注意

阈值的选择部分是判断和领域知识的问题，并且在存在显著类别不平衡的情况下主要是一个问题。然而，在第十章**，逻辑回归中，我们将探讨如何计算一个最优阈值。

本节以及上一节演示了如何评估二元分类模型。它们表明模型评估不仅仅是一个简单的赞成或反对的过程。它更像是在做蛋糕时品尝面糊。我们对模型规格做出良好的初始假设，并使用模型评估过程进行改进。这通常涉及准确度、灵敏度、特异性和精确度之间的权衡，以及抵制一刀切建议的建模决策。这些决策在很大程度上取决于特定领域，并且是专业判断的问题。

本节中的讨论以及大多数技术同样适用于多类建模。我们将在下一节讨论如何评估多类模型。

评估多类模型

我们用来评估二元分类模型的所有相同原则也适用于多类模型评估。计算混淆矩阵同样重要，尽管它更难解释。我们还需要检查一些相互竞争的指标，如精确度和灵敏度。这也比二元分类更复杂。

再次，我们将使用 NLS 学位完成数据。在这种情况下，我们将目标从学士学位完成与否更改为高中完成、学士学位完成和研究生学位完成：

我们将首先加载必要的库。这些库与前面两节中使用的相同：

import pandas as pd
import numpy as np
from feature_engine.encoding import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
import sklearn.metrics as skmet
import matplotlib.pyplot as plt

接下来，我们将加载 NLS 学位达成数据，创建训练和测试数据框，并对数据进行编码和缩放：

nls97degreelevel = pd.read_csv("data/nls97degreelevel.csv")
feature_cols = ['satverbal','satmath','gpaoverall',
  'parentincome','gender']
X_train, X_test, y_train, y_test =  \
  train_test_split(nls97degreelevel[feature_cols],\
  nls97degreelevel[['degreelevel']], test_size=0.3, random_state=0)
ohe = OneHotEncoder(drop_last=True, variables=['gender'])
ohe.fit(X_train)
X_train_enc, X_test_enc = \
  ohe.transform(X_train), ohe.transform(X_test)
scaler = StandardScaler()
standcols = X_train_enc.iloc[:,:-1].columns
scaler.fit(X_train_enc[standcols])
X_train_enc = \
  pd.DataFrame(scaler.transform(X_train_enc[standcols]),
  columns=standcols, index=X_train_enc.index).\
  join(X_train_enc[['gender_Female']])
X_test_enc = \
  pd.DataFrame(scaler.transform(X_test_enc[standcols]),
  columns=standcols, index=X_test_enc.index).\
  join(X_test_enc[['gender_Female']])

现在，我们将运行一个 KNN 模型，并为每个学位级别类别预测值：

knn = KNeighborsClassifier(n_neighbors = 5)
knn.fit(X_train_enc, y_train.values.ravel())
pred = knn.predict(X_test_enc)
pred_probs = knn.predict_proba(X_test_enc)[:, 1]

我们可以使用这些预测来生成一个混淆矩阵：

cm = skmet.confusion_matrix(y_test, pred)
cmplot = skmet.ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=['High School', 'Bachelor','Post-Graduate'])
cmplot.plot()
cmplot.ax_.set(title='Confusion Matrix', 
  xlabel='Predicted Value', ylabel='Actual Value')

这生成了以下图表：

图 6.15 – 具有多类目标的混淆矩阵

可以手动计算评估指标。精度是我们预测中实际属于类内的百分比。因此，对于我们的高中预测，它是 48 / (48 + 38 + 8) = 0.51。高中类别的灵敏度——即我们模型预测的实际高中值的百分比——是 48 / (48 + 19 + 5) = 0.67。然而，这相当繁琐。幸运的是，scikit-learn 可以为我们完成这项工作。

我们可以通过调用 classification_report 方法来获取这些统计数据，传递实际和预测值（记住召回率和灵敏度是相同的度量）：

print(skmet.classification_report(y_test, pred,
  target_names=['High School', 'Bachelor', 'Post-Graduate']))
               precision    recall  f1-score   support
  High School       0.51      0.67      0.58        72
     Bachelor       0.51      0.49      0.50        92
Post-Graduate       0.42      0.24      0.30        42
     accuracy                           0.50       206
    macro avg       0.48      0.46      0.46       206
 weighted avg       0.49      0.50      0.49       206

除了按类别计算的精度和灵敏度比率之外，我们还会得到一些其他统计数据。F1 分数是精度和灵敏度的调和平均值。

在这里，p 代表精度，而 s 代表灵敏度。

要获得类别的平均精度、灵敏度和 F1 分数，我们可以使用简单的平均（宏平均）或调整类别大小的加权平均。使用加权平均，我们得到精度、灵敏度和 F1 分数分别为 0.49、0.50 和 0.49。（由于这里的类别相对平衡，宏平均和加权平均之间没有太大差异。）

这演示了如何将我们讨论的二分类模型评估指标扩展到多类评估。虽然实现起来更困难，但相同的概念和技术同样适用。

到目前为止，我们关注的是指标和可视化，以帮助我们评估分类模型。我们尚未检查评估回归模型的指标。这些指标可能比分类指标更为直接。我们将在下一节中讨论它们。

评估回归模型

回归模型评估的指标通常基于目标变量的实际值和模型预测值之间的距离。最常见的指标——均方误差、均方根误差、平均绝对误差和 R 平方——都追踪我们的预测如何成功地捕捉目标变量的变化。

实际值和我们的预测值之间的距离被称为残差或误差。均方误差（MSE）是残差平方的平均值：

在这里，是第 i 次观察的实际目标变量值，而是我们对目标值的预测。由于预测值高于实际值，残差被平方以处理负值。为了使我们的测量值返回到一个更有意义的尺度，我们通常使用均方误差（MSE）的平方根。这被称为均方根误差（RMSE）。

由于平方，均方误差（MSE）将对较大的残差进行更多的惩罚，而不是较小的残差。例如，如果我们对五个观测值进行预测，其中一个残差为 25，其他四个残差为 0，我们将得到均方误差为 (0+0+0+0+625)/5 = 125。然而，如果所有五个观测值的残差都是 5，均方误差将是 (25+25+25+25+25)/5 = 25。

将残差的平方作为替代方案是取它们的绝对值。这给我们带来了平均绝对误差：

R-squared，也称为确定系数，是我们模型捕获目标变量变化的估计比例。我们平方残差，就像我们在计算均方误差（MSE）时做的那样，并将其除以每个实际目标值与其样本均值之间的偏差。这给我们带来了仍然未解释的变异，我们从 1 中减去它以得到解释的变异：

幸运的是，scikit-learn 使得生成这些统计数据变得容易。在本节中，我们将构建一个关于陆地温度的线性回归模型，并使用这些统计数据来评估它。我们将使用来自美国国家海洋和大气管理局 2019 年气象站平均年度温度、海拔和纬度的数据。

注意

陆地温度数据集包含了 2019 年来自全球超过 12,000 个站点的平均温度读数（以摄氏度为单位），尽管大多数站点位于美国。原始数据是从全球历史气候学网络综合数据库中检索的。它已由美国国家海洋和大气管理局在www.ncdc.noaa.gov/data-access/land-based-station-data/land-based-datasets/global-historical-climatology-network-monthly-version-4上提供给公众使用。

让我们开始构建一个线性回归模型：

我们将首先加载所需的库和陆地温度数据。我们还将创建训练和测试数据框：

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import sklearn.metrics as skmet
import matplotlib.pyplot as plt
landtemps = pd.read_csv("data/landtemps2019avgs.csv")
feature_cols = ['latabs','elevation']
X_train, X_test, y_train, y_test =  \
  train_test_split(landtemps[feature_cols],\
  landtemps[['avgtemp']], test_size=0.3, random_state=0)

注意

latabs 特征是纬度的值，不带北或南指示符；因此，埃及开罗大约在北纬 30 度，巴西的波尔图阿雷格里大约在南纬 30 度，它们具有相同的值。

现在，我们缩放我们的数据：

scaler = StandardScaler()
scaler.fit(X_train)
X_train = \
  pd.DataFrame(scaler.transform(X_train),
  columns=feature_cols, index=X_train.index)
X_test = \
  pd.DataFrame(scaler.transform(X_test),
  columns=feature_cols, index=X_test.index)
scaler.fit(y_train)
y_train, y_test = \
  pd.DataFrame(scaler.transform(y_train),
  columns=['avgtemp'], index=y_train.index),\
  pd.DataFrame(scaler.transform(y_test),
  columns=['avgtemp'], index=y_test.index)

接下来，我们实例化一个 scikit-learn LinearRegression 对象，并在训练数据上拟合一个模型。我们的目标是年度平均温度 (avgtemp)，而特征是纬度 (latabs) 和 elevation。coef_ 属性给我们每个特征的系数：
```
lr = LinearRegression()
lr.fit(X_train, y_train)
np.column_stack((lr.coef_.ravel(),
  X_test.columns.values))
array([[-0.8538957537748768, 'latabs'],
       [-0.3058979822791853, 'elevation']], dtype=object)
```

latabs系数的解释是，标准化平均年温度每增加一个标准差，将下降 0.85。（LinearRegression模块不返回 p 值，这是系数估计的统计显著性的度量。你可以使用statsmodels来查看普通最小二乘模型的完整摘要。）

现在，我们可以获取预测值。让我们也将测试数据中的特征和目标与返回的 NumPy 数组结合起来。然后，我们可以通过从实际值（avgtemp）中减去预测值来计算残差。尽管存在轻微的负偏斜和过度的峰度，但残差看起来还不错：

pred = lr.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.021
median           0.032
skew              -0.641
kurtosis        6.816
Name: resid, dtype: float64

值得注意的是，在本书中，我们大多数时候在处理回归模型时都会以这种方式生成预测值和计算残差。如果你对前面代码块中我们刚刚做的事情感到有些不清楚，再次回顾一下可能是个好主意。

我们应该绘制残差图，以更好地了解它们的分布。

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")

这会产生以下图表：

图 6.16 – 线性回归模型的残差直方图

这看起来并不太糟糕，但我们有更多的正残差，在我们预测的测试数据中的温度低于实际温度的情况下，比负残差更多。

通过残差绘制我们的预测可能让我们更好地理解正在发生的情况：

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")

这会产生以下图表：

图 6.17 – 线性回归模型的预测残差散点图

这看起来并不糟糕。残差在 0 附近随机波动。然而，在 1 到 2 个标准差之间的预测值更有可能过低（具有正残差），而不是过高。超过 2，预测值总是过高（它们具有负残差）。这个模型线性假设可能并不合理。我们应该探索我们在第四章中讨论的一些转换，或者尝试一个非参数模型，如 KNN 回归。

也可能极端值在相当大的程度上拉动了我们的系数。一个不错的下一步可能是移除异常值，正如我们在第一章的识别极端值和异常值部分所讨论的，检查特征和目标的分布。然而，我们在这里不会这么做。

让我们看看一些评估指标。这可以通过 scikit-learn 的metrics库轻松完成。我们可以调用相同的函数来获取 RMSE 作为 MSE。我们只需要将平方参数设置为False：

mse = skmet.mean_squared_error(y_test, pred)
mse
0.18906346144036693
rmse = skmet.mean_squared_error(y_test, pred, squared=False)
rmse
0.4348142838504353
mae = skmet.mean_absolute_error(y_test, pred)
mae
0.318307379728143
r2 = skmet.r2_score(y_test, pred)
r2
0.8162525715296725

标准差以下 0.2 的均方误差（MSE）和标准差以下 0.3 的绝对误差（MAE）看起来相当不错，尤其是对于这样一个稀疏模型。R-squared 超过 80%也是相当有希望的。

让我们看看如果我们使用 KNN 模型会得到什么结果：

knn = KNeighborsRegressor(n_neighbors=5)
knn.fit(X_train, y_train)
pred = knn.predict(X_test)
mae = skmet.mean_absolute_error(y_test, pred)
mae
0.2501829988751876
r2 = skmet.r2_score(y_test, pred)
r2
0.8631113217183314

这个模型实际上在 MAE 和 R-squared 方面都有所改进。

我们也应该再次审视残差：

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 with KNN Model")
plt.xlabel("Predicted Temperature")
plt.ylabel("Residuals")
plt.show()

这会产生以下图表：

图 6.18 – KNN 模型的预测残差散点图

这个残差图表看起来也好多了。在目标分布的任何部分，我们都不太可能过度预测或低估。

本节介绍了评估回归模型的关键措施及其解释方法。同时，还展示了如何通过可视化，尤其是模型残差的可视化，来提高这种解释的准确性。

然而，到目前为止，我们在使用回归和分类度量时都受到了我们构建的训练和测试数据框的限制。如果，出于某种原因，测试数据在某些方面不寻常呢？更普遍地说，我们基于什么结论认为我们的评估措施是准确的？如果我们使用 K 折交叉验证，我们可以更有信心地使用这些措施，我们将在下一节中介绍。

使用 K 折交叉验证

到目前为止，我们已经保留了 30%的数据用于验证。这不是一个坏策略。它防止我们在训练模型时提前查看测试数据。然而，这种方法并没有充分利用所有可用数据，无论是用于训练还是测试。如果我们使用 K 折交叉验证，我们就可以使用所有数据，同时避免数据泄露。也许这听起来太好了，但事实并非如此，这并不是因为一个巧妙的小技巧。

K 折交叉验证在 K 个折（或部分）中的所有但一个上训练我们的模型，留出一个用于测试。这重复k次，每次排除一个不同的折用于测试。性能指标是基于 K 个折的平均分数。

在开始之前，我们还需要再次考虑数据泄露的可能性。如果我们对我们将用于训练模型的所有数据进行缩放，然后将其分成折，我们将在训练中使用所有折的信息。为了避免这种情况，我们需要在每个迭代中仅对训练折进行缩放，以及进行任何其他预处理。虽然我们可以手动完成这项工作，但 scikit-learn 的pipeline库可以为我们做很多这项工作。我们将在本节中介绍如何使用管道进行交叉验证。

让我们尝试使用 K 折交叉验证评估上一节中指定的两个模型。同时，我们也来看看随机森林回归器可能表现如何：

除了我们迄今为止已经使用的库之外，我们还需要 scikit-learn 的make_pipeline、cross_validate和Kfold库：

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_validate
from sklearn.model_selection import KFold

我们再次加载陆地温度数据并创建训练和测试 DataFrame。我们仍然想留出一些数据用于最终验证，但这次，我们只留出 10%。我们将使用剩余的 90%进行训练和测试：

landtemps = pd.read_csv("data/landtemps2019avgs.csv")
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)

现在，我们创建一个KFold对象，并指出我们想要五个折，并且数据要打乱（如果数据尚未随机排序，打乱数据是一个好主意）：
```
kf = Kfold(n_splits=5, shuffle=True, random_state=0)
```

接下来，我们定义一个函数来创建一个管道。然后，该函数运行cross_validate，它接受管道和我们之前创建的KFold对象：

def getscores(model):
  pipeline = make_pipeline(StandardScaler(), model)
  scores = cross_validate(pipeline, X=X_train, 
    y=y_train, cv=kf, scoring=['r2'], n_jobs=1)
  scorelist.append(dict(model=str(model),
    fit_time=scores['fit_time'].mean(),
    r2=scores['test_r2'].mean()))

现在，我们准备调用线性回归、随机森林回归和 KNN 回归模型的getscores函数：

scorelist = []
getscores(LinearRegression())
getscores(RandomForestRegressor(max_depth=2))
getscores(KNeighborsRegressor(n_neighbors=5))

我们可以将scorelist列表打印出来查看我们的结果：

scorelist
[{'model': 'LinearRegression()',
  'fit_time': 0.004968833923339844,
  'r2': 0.8181125031214872},
 {'model': 'RandomForestRegressor(max_depth=2)',
  'fit_time': 0.28124608993530276,
  'r2': 0.7122492698889024},
 {'model': 'KNeighborsRegressor()',
  'fit_time': 0.006945991516113281,
  'r2': 0.8686733636724104}]

根据 R-squared 值，KNN 回归模型比线性回归或随机森林回归模型表现更好。随机森林回归模型也有一个显著的缺点，那就是它的拟合时间要长得多。

使用管道预处理数据

在上一节中，我们只是触及了 scikit-learn 管道可以做的事情的表面。我们经常需要将所有的预处理和特征工程折叠到一个管道中，包括缩放、编码以及处理异常值和缺失值。这可能很复杂，因为不同的特征可能需要不同的处理。我们可能需要用数值特征的中间值和分类特征的众数来填充缺失值。我们可能还需要转换我们的目标变量。我们将在本节中探讨如何做到这一点。

按照以下步骤：

我们将首先加载本章中已经使用过的库。然后，我们将添加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 feature_engine.encoding import OneHotEncoder
from sklearn.impute import KNNImputer
from sklearn.model_selection import cross_validate, KFold
import sklearn.metrics as skmet
from sklearn.compose import ColumnTransformer
from sklearn.compose import TransformedTargetRegressor

列转换器非常灵活。我们甚至可以使用我们自定义的预处理函数。以下代码块从helperfunctions子文件夹中的preprocfunc模块导入OutlierTrans类：
```
import os
import sys
sys.path.append(os.getcwd() + "/helperfunctions")
from preprocfunc import OutlierTrans
```
OutlierTrans类通过距离四分位数范围来识别极端值。这是我们演示过的技术，见第三章，识别和修复缺失值。

要在 scikit-learn 管道中工作，我们的类必须具有 fit 和 transform 方法。我们还需要继承BaseEstimator和TransformerMixin类。

在这个类中，几乎所有操作都在transform方法中完成。任何超过第三四分位数 1.5 倍或低于第一四分位数的值都被分配为缺失值：

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)
      inlierrange = self.threshold*(thirdq-firstq)
      outlierhigh, outlierlow = inlierrange+thirdq,\
        firstq-inlierrange
      Xnew.loc[(Xnew[col]>outlierhigh) | \
        (Xnew[col]<outlierlow),col] = np.nan
    return Xnew.values

我们的OutlierTrans类可以在我们的管道中以与我们在上一节中使用StandardScaler相同的方式使用。我们将在稍后这样做。

现在，我们已经准备好加载需要处理的数据。在本节中，我们将使用 NLS 每周工资数据。每周工资将是我们的目标，我们将使用高中 GPA、母亲和父亲最高学历、家庭收入、性别以及个人是否完成学士学位作为特征。

我们将创建一个特征列表，以不同的方式处理这些特征。这将在我们指导管道对数值、分类和二进制特征执行不同操作时很有帮助：

nls97wages = pd.read_csv("data/nls97wagesb.csv")
nls97wages.set_index("personid", inplace=True)
nls97wages.dropna(subset=['wageincome'], inplace=True)
nls97wages.loc[nls97wages.motherhighgrade==95,
  'motherhighgrade'] = np.nan
nls97wages.loc[nls97wages.fatherhighgrade==95,
  'fatherhighgrade'] = np.nan
num_cols = ['gpascience','gpaenglish','gpamath','gpaoverall',
  'motherhighgrade','fatherhighgrade','parentincome']
cat_cols = ['gender']
bin_cols = ['completedba']
target = nls97wages[['wageincome']]
features = nls97wages[num_cols + cat_cols + bin_cols]
X_train, X_test, y_train, y_test =  \
  train_test_split(features,\
  target, test_size=0.2, random_state=0)

让我们看看一些描述性统计。一些变量有超过一千个缺失值（gpascience、gpaenglish、gpamath、gpaoverall和parentincome）：

nls97wages[['wageincome'] + num_cols].agg(['count','min','median','max']).T
                 count    min      median    max
wageincome       5,091    0        40,000    235,884
gpascience       3,521    0        284       424
gpaenglish       3,558    0        288       418
gpamath          3,549    0        280       419
gpaoverall       3,653    42       292       411
motherhighgrade  4,734    1        12        20
fatherhighgrade  4,173    1        12        29
parentincome     3,803   -48,100   40,045    246,474

现在，我们可以设置一个列转换器。首先，我们将为处理数值数据（standtrans）、分类数据和二进制数据创建管道。

对于数值数据，我们希望将异常值视为缺失。在这里，我们将2这个值传递给OutlierTrans的阈值参数，表示我们希望将高于或低于该范围两倍四分位距的值设置为缺失。记住，默认值是 1.5，所以我们相对保守一些。

然后，我们将创建一个ColumnTransformer对象，将其传递给刚刚创建的三个管道，并指明使用哪个管道来处理哪些特征：

standtrans = make_pipeline(OutlierTrans(2),
  StandardScaler())
cattrans = make_pipeline(SimpleImputer(strategy="most_frequent"),
  OneHotEncoder(drop_last=True))
bintrans = make_pipeline(SimpleImputer(strategy="most_frequent"))
coltrans = ColumnTransformer(
  transformers=[
    ("stand", standtrans, num_cols),
    ("cat", cattrans, ['gender']),
    ("bin", bintrans, ['completedba'])
  ]
)

现在，我们可以将列转换器添加到包含我们想要运行的线性模型的管道中。我们将向管道中添加 KNN 插补来处理缺失值。

我们还需要对目标进行缩放，这在我们管道中无法完成。我们将使用 scikit-learn 的TransformedTargetRegressor来完成这个任务。我们将刚刚创建的管道传递给目标回归器的regressor参数：

lr = LinearRegression()
pipe1 = make_pipeline(coltrans,
  KNNImputer(n_neighbors=5), lr)
ttr=TransformedTargetRegressor(regressor=pipe1,
  transformer=StandardScaler())

让我们使用这个管道进行 K 折交叉验证。我们可以通过目标回归器ttr将我们的管道传递给cross_validate函数：

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)
print("Mean Absolute Error: %.2f, R-squared: %.2f" % 
  (scores['test_neg_mean_absolute_error'].mean(),
  scores['test_r2'].mean()))
Mean Absolute Error: -23781.32, R-squared: 0.20

虽然这些分数并不理想，但这并不是这次练习的真正目的。关键要点是我们通常希望将大部分预处理工作整合到管道中。这是避免数据泄露的最佳方式。列转换器是一个极其灵活的工具，允许我们对不同的特征应用不同的转换。

摘要

本章介绍了关键模型评估指标和技术，以便在本书的剩余章节中广泛使用和扩展时，它们将变得熟悉。我们研究了分类和回归模型评估的非常不同的方法。我们还探讨了如何使用可视化来改进我们对预测的分析。最后，我们使用管道和交叉验证来获取模型性能的可靠估计。

我希望这一章也给了你一个机会去适应这本书接下来的一般方法。尽管在接下来的章节中我们将讨论大量算法，但我们仍将继续探讨我们在前几章中讨论的预处理问题。当然，我们将讨论每个算法的核心概念。但是，以真正的动手实践方式，我们还将处理现实世界数据的杂乱无章。每一章将从相对原始的数据开始，到特征工程，再到模型指定和模型评估，高度依赖 scikit-learn 的管道来整合所有内容。

在接下来的几章中，我们将讨论回归算法——那些允许我们建模连续目标的算法。我们将探讨一些最受欢迎的回归算法——线性回归、支持向量回归、K 最近邻回归和决策树回归。我们还将考虑对回归模型进行修改，以解决欠拟合和过拟合问题，包括非线性变换和正则化。