本文将基于实际数据集，探讨一种先进的方法来解答哪些特征（数据集的列）阻碍了模型在新数据上的有效泛化。
机器学习的核心目标是在未见过的新数据上实现准确预测。
当模型在训练数据上表现良好，但在测试数据上表现不佳时，即出现“过拟合”。这意味着模型从训练数据中学习了过多的噪声模式，从而丧失了在新数据上的泛化能力。
那么，过拟合的根本原因是什么？具体来说，
哪些特征（数据集的列）阻碍了模型在新数据上的有效泛化？
本文将基于实际数据集，探讨一种先进的方法来解答这一问题。
特征重要性在此场景下不再适用
如果你的第一反应是“我会查看特征重要性”，那么请重新考虑。
特征重要性无法直接反映特征在新数据上的表现。
实际上，特征重要性仅是模型在训练阶段所学内容的表现。如果模型在训练过程中学习到关于“年龄”特征的复杂模式，那么该特征的特征重要性将会很高。但这并不意味着
这些
模式是准确的
（“准确”指的是一种具备泛化能力的模式，即在新的数据上依然成立）。
因此，我们需要采用不同的方法来解决这个问题。
案例
为了阐述该方法，我将使用一个包含1984年至1988年德国健康登记数据的数据集（该数据集可通过 Pydataset 库获得，并遵循 MIT 许可证）。
下载数据的方式非常简单：
import pydataset X = pydataset.data('rwm5yr') y = (X['hospvis'] > 1).astype(int) X = X.drop('hospvis', axis = 1)
该数据集包含19,609行，每行记录了一名患者在特定年份的一些信息。需要注意的是，患者在不同年份都可能被观测到，因此同一患者可能出现在数据框的多行中。
目标变量是：
我们拥有16个特征列：
-
**id**：患者ID（1-7028）；
-
**docvis**：年内就诊次数（0-121）；
-
**year**：年份（1984-1988）；
-
**edlevel**：教育水平（1-4）；
-
**age**：年龄（25-64）；
-
**outwork**：是否失业，1表示失业，0表示在职；
-
**female**：性别，1表示女性，0表示男性；
-
**married**：婚姻状况，1表示已婚，0表示未婚；
-
**kids**：是否有子女，1表示有，0表示无；
-
**hhninc**：家庭年收入（马克）；
-
**educ**：受教育年限（7-18）；
-
**self**：是否自雇，1表示自雇，0表示非自雇；
-
**edlevel1**：是否为高中以下学历，1表示是，0表示否；
-
**edlevel2**：是否为高中学历，1表示是，0表示否；
-
**edlevel3**：是否为大学/学院学历，1表示是，0表示否；
-
**edlevel4**：是否为研究生学历，1表示是，0表示否。
现在，将数据集划分为训练集和测试集。虽然存在更严谨的方法，如交叉验证，但为了保持简洁，我们采用简单的划分方法。但在本文中我们（简单地）将所有列视为数值特征。
from sklearn.model_selection import train_test_split from catboost import CatBoostClassifier
X_train, X_test, y_train, y_test = train_test_split( X, y, test_size = .2, stratify = y)
cat = CatBoostClassifier(silent = True).fit(X_train, y_train)
模型训练完成后，我们来分析特征重要性：
import pandas as pd
fimpo = pd.Series(cat.feature_importances_, index = X_train.columns)
已训练模型的特征重要性。
不出所料，docvis（就诊次数）在预测患者是否住院超过1天方面至关重要。age（年龄）和hhninc（收入）这两个特征的重要性也符合预期。但是，患者ID（id）的重要性排名第二则值得警惕，尤其是在我们将其视为数值特征的情况下。
接下来，计算模型在训练集和测试集上的性能指标（ROC曲线下面积，AUC）：
from sklearn.metrics import roc_auc_score roc_train = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1]) roc_test = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1])
模型在训练集和测试集上的性能。
结果显示，训练集和测试集之间的性能差距显著。
这表明存在明显的过拟合现象。
那么，究竟是哪些特征导致了过拟合？
什么是SHAP值？
我们有多种指标可以衡量
模型
在特定数据上的表现。但如何衡量
特征
在特定数据上的表现？
“SHAP值”是解决此问题的有力工具。
通常，你可以使用专门的Python库来高效计算任何预测模型的SHAP值。但这里为了简单，我们将利用Catboost的原生方法：
from catboost import Pool
shap_train = pd.DataFrame( data = cat.get_feature_importance( data = Pool(X_train), type = 'ShapValues')[:, :-1], index = X_train.index, columns = X_train.columns )
shap_test = pd.DataFrame( data = cat.get_feature_importance( data = Pool(X_test), type = 'ShapValues')[:, :-1], index = X_test.index, columns = X_test.columns )
观察shap_train和shap_test，你会发现它们的形状与各自的数据集相同。
SHAP值可以量化
每个特征对模型在单个或多个观测值上的最终预测的影响。
看几个例子：
原始数据及其对应的SHAP值。
第12071行的患者就诊次数为0，相应的SHAP值为-0.753，这意味着该信息将患者住院超过1天的概率（实际上是对数几率）降低了0.753。相反，第18650行的患者就诊4次，这使得她住院超过1天的对数几率提高了0.918。
认识ParShap
一个特征在数据集上的性能可以通过
该特征的SHAP值与目标变量之间的相关性
来近似表示。如果模型在某个特征上学习到有效的模式，那么该特征的SHAP值应与目标变量高度正相关。
例如，如果我们想计算docvis特征与测试集中观测数据的目标变量之间的相关性：
import numpy as np
np.corrcoef(shap_test['docvis'], y_test)
然而，SHAP值具有可加性，即最终预测是所有特征SHAP值的总和。因此在
计算相关性之前，先消除其他特征的影响
会更有意义。这正是
“偏相关”
的定义。偏相关的便捷实现方式可以在 Python 库 Pingouin 中找到：
import pingouin
pingouin.partial_corr( data = pd.concat([shap_test, y_test], axis = 1).astype(float), x = 'docvis', y = y_test.name, x_covar = [feature for feature in shap_test.columns if feature != 'docvis'] )
这段代码的含义是：“计算docvis特征的SHAP值与测试集观测数据的目标变量之间的相关性，同时消除所有其他特征的影响。”
为了方便将此公式称为
“ParShap”
（Partial correlation of Shap values，SHAP值的偏相关）。
我们可以对训练集和测试集中的每个特征重复此过程：
parshap_train = partial_correlation(shap_train, y_train) parshap_test = partial_correlation(shap_test, y_test)
注意：你可以在本文末尾找到 partial_correlation 函数的定义。
现在在 x 轴上绘制 parshap_train，在 y 轴上绘制 parshap_test。
import matplotlib.pyplot as plt
plt.scatter(parshap_train, parshap_test)
SHAP值与目标变量的偏相关性，分别在训练集和测试集上。注意：颜色条表示特征重要性。[作者提供的图片]
如果一个特征位于对角线上，则表示它在训练集和测试集上的表现完全一致。
这是理想情况，既没有过拟合也没有欠拟合。反之如果一个特征位于对角线下方，则表示它在测试集上的表现不如训练集。这属于过拟合区域。
通过视觉观察，我们可以立即发现哪些特征表现不佳：我已用蓝色圆圈标记出它们。
蓝色圆圈标记的特征是当前模型中最容易出现过拟合的特征。[作者提供的图片]
因此，parshap_test和parshap_train之间的算术差（等于每个特征与对角线之间的垂直距离）可以量化该特征对模型的过拟合程度。
parshap_diff = parshap_test - parshap_train
parshap_test和parshap_train之间的算术差。[作者提供的图片]
应该如何解读这个结果？基于以上分析，
该值越负，则该特征导致的过拟合程度越高。
验证
我们是否能找到一种方法来验证本文提出的观点的正确性？
从逻辑上讲，如果从数据集中移除“过拟合特征”，应该能够减少过拟合现象（即，缩小roc_train和roc_test之间的差距）。
因此尝试每次删除一个特征，并观察ROC曲线下面积的变化。
根据特征重要性（左）或ParShap（右）排序，依次删除一个特征时，模型在训练集和测试集上的性能。
在
左侧
的图中，每次移除一个特征，并按照
特征重要性
进行排序。首先移除最不重要的特征（edlevel4），然后移除两个最不重要的特征（edlevel4和edlevel1），以此类推。
在
右侧
的图中，执行相同的操作，但是移除的顺序由
ParShap
决定。首先移除ParShap值最小（最负）的特征（id），然后移除两个ParShap值最小的特征（id和year），以此类推。
正如预期的那样，移除ParShap值最负的特征显著减少了过拟合。实际上，roc_train的值逐渐接近roc_test的值。
需要注意的是，
这仅仅是用于验证我们推理过程的测试。
通常来说，ParShap
不应被用作特征选择的方法。
因为，
某些特征容易出现过拟合并不意味着这些特征完全没有用处！
（例如，本例中的收入和年龄）。
但是ParShap在为我们提供模型调试的线索方面非常有用。它可以帮助我们将注意力集中在那些需要更多特征工程或正则化的特征上。
最后以下是本文的完整，有兴趣的可以进行复现测试
本文使用的完整代码（由于随机种子，结果可能略有不同）：
import pandas as pd import pydataset from sklearn.model_selection import train_test_split from catboost import CatBoostClassifier, Pool from sklearn.metrics import roc_auc_score from pingouin import partial_corr import matplotlib.pyplot as plt print('################# Print docs') pydataset.data('rwm5yr', show_doc = True) X = pydataset.data('rwm5yr') y = (X['hospvis'] > 1).astype(int) X = X.drop('hospvis', axis = 1) X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = .2, stratify = y) cat = CatBoostClassifier(silent = True).fit(X_train, y_train) fimpo = pd.Series(cat.feature_importances_, index = X_train.columns) fig, ax = plt.subplots() fimpo.sort_values().plot.barh(ax = ax) fig.savefig('fimpo.png', dpi = 200, bbox_inches="tight") fig.show() roc_train = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1]) roc_test = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1]) print('\n################# Print roc') print('roc_auc train: {:.2f}'.format(roc_train)) print('roc_auc test: {:.2f}'.format(roc_test)) shap_train = pd.DataFrame( data = cat.get_feature_importance(data = Pool(X_train), type = 'ShapValues')[:, :-1], index = X_train.index, columns = X_train.columns ) shap_test = pd.DataFrame( data = cat.get_feature_importance(data = Pool(X_test), type = 'ShapValues')[:, :-1], index = X_test.index, columns = X_test.columns ) print('\n################# Print df shapes') print(f'X_train.shape: {X_train.shape}') print(f'X_test.shape: {X_test.shape}\n') print(f'shap_train.shape: {shap_train.shape}') print(f'shap_test.shape: {shap_test.shape}') print('\n################# Print data and SHAP') print('Original data:') display(X_test.head(3)) print('\nCorresponding SHAP values:') display(shap_test.head(3).round(3)) def partial_correlation(X, y): out = pd.Series(index = X.columns, dtype = float) for feature_name in X.columns: out[feature_name] = partial_corr( data = pd.concat([X, y], axis = 1).astype(float), x = feature_name, y = y.name, x_covar = [f for f in X.columns if f != feature_name] ).loc['pearson', 'r'] return out parshap_train = partial_correlation(shap_train, y_train) parshap_test = partial_correlation(shap_test, y_test) parshap_diff = pd.Series(parshap_test - parshap_train, name = 'parshap_diff') print('\n################# Print parshap_diff') print(parshap_diff.sort_values()) plotmin, plotmax = min(parshap_train.min(), parshap_test.min()), max(parshap_train.max(), parshap_test.max()) plotbuffer = .05 * (plotmax - plotmin) fig, ax = plt.subplots() if plotmin < 0: ax.vlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color = 'darkgrey', zorder = 0) ax.hlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color = 'darkgrey', zorder = 0) ax.plot( [plotmin - plotbuffer, plotmax + plotbuffer], [plotmin - plotbuffer, plotmax + plotbuffer], color = 'darkgrey', zorder = 0 ) sc = ax.scatter( parshap_train, parshap_test, edgecolor = 'grey', c = fimpo, s = 50, cmap = plt.cm.get_cmap('Reds'), vmin = 0, vmax = fimpo.max()) ax.set(title = 'Partial correlation bw SHAP and target...', xlabel = '... on Train data', ylabel = '... on Test data') cbar = fig.colorbar(sc) cbar.set_ticks([]) for txt in parshap_train.index: ax.annotate(txt, (parshap_train[txt], parshap_test[txt] + plotbuffer / 2), ha = 'center', va = 'bottom') fig.savefig('parshap.png', dpi = 300, bbox_inches="tight") fig.show() n_drop_max = 5 iterations = 4 features = {'parshap': parshap_diff, 'fimpo': fimpo} features_dropped = {} roc_auc_scores = { 'fimpo': {'train': pd.DataFrame(), 'test': pd.DataFrame()}, 'parshap': {'train': pd.DataFrame(), 'test': pd.DataFrame()} } for type_ in ['parshap', 'fimpo']: for n_drop in range(n_drop_max + 1): features_drop = features[type_].sort_values().head(n_drop).index.to_list() features_dropped[type_] = features_drop X_drop = X.drop(features_drop, axis = 1) for i in range(iterations): X_train, X_test, y_train, y_test = train_test_split(X_drop, y, test_size = .2, stratify = y) cat = CatBoostClassifier(silent = True).fit(X_train, y_train) roc_auc_scores[type_]['train'].loc[n_drop, i] = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1]) roc_auc_scores[type_]['test'].loc[n_drop, i] = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1]) fig, axs = plt.subplots(1, 2, sharey = True, figsize = (8, 3)) plt.subplots_adjust(wspace = .1) axs[0].plot(roc_auc_scores['fimpo']['train'].index, roc_auc_scores['fimpo']['train'].mean(axis = 1
