注意
前往结尾 下载完整的示例代码。或通过 JupyterLite 或 Binder 在浏览器中运行此示例。
线性模型系数解释中的常见陷阱#
在线性模型中，目标值被建模为特征的线性组合（有关 scikit-learn 中可用的一组线性模型的描述，请参见线性模型用户指南部分）。多元线性模型中的系数表示给定特征 \(X_i\) 和目标 \(y\) 之间的关系，假设所有其他特征保持不变（条件依赖性）。这与绘制 \(X_i\) 与 \(y\) 并拟合线性关系不同：在这种情况下，所有其他特征的可能值都考虑在估计中（边缘依赖性）。
此示例将提供一些解释线性模型系数的提示，指出当线性模型不适合描述数据集，或当特征相关时出现的问题。
注意
请记住，特征 \(X\) 和结果 \(y\) 通常是未知于我们的数据生成过程的结果。机器学习模型经过训练，可以根据样本数据逼近连接 \(X\) 和 \(y\) 的未观察到的数学函数。因此，对模型作出的任何解释不一定能推广到真实的数据生成过程。当模型质量差或样本数据不能代表总体时，尤其如此。
我们将使用来自 1985 年“当前人口调查”的数据来预测工资作为经验、年龄或教育等各种特征的函数。
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
数据集：工资#
我们从OpenML获取数据。请注意，将参数as_frame设置为 True 将检索数据作为 pandas 数据框。
from sklearn.datasets import fetch_openml
survey = fetch_openml(data_id=534, as_frame=True)
然后，我们识别特征X和目标y：WAGE 列是我们的目标变量（即我们要预测的变量）。
X = survey.data[survey.feature_names]
X.describe(include="all")
请注意，数据集包含分类变量和数值变量。之后预处理数据集时，我们需要考虑这一点。
X.head()
我们的预测目标：工资。工资以美元/小时的浮点数表示。
y = survey.target.values.ravel()
survey.target.head()
0 5.10
1 4.95
2 6.67
3 4.00
4 7.50
Name: WAGE, dtype: float64
我们将样本分成训练集和测试集。以下探索性分析中只使用训练集。这是一种模拟实际情况的方法，在这种情况下，预测是在未知目标上进行的，我们不希望我们的分析和决策受到我们对测试数据的了解的偏差。
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
首先，让我们通过查看变量分布及其之间的成对关系来获得一些见解。只使用数值变量。在下图中，每个点代表一个样本。
train_dataset = X_train.copy()
train_dataset.insert(0, "WAGE", y_train)
_ = sns.pairplot(train_dataset, kind="reg", diag_kind="kde")
仔细观察工资分布表明它有一个长尾。因此，我们应该取其对数，使其近似于正态分布（岭回归或 lasso 等线性模型最适合误差的正态分布）。
当教育程度增加时，工资也在增加。请注意，此处表示的工资与教育程度之间的依赖性是边缘依赖性，即它描述了特定变量的行为，而没有保持其他变量不变。
此外，经验和年龄高度线性相关。
机器学习管道#
为了设计我们的机器学习管道，我们首先手动检查我们正在处理的数据类型
survey.data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 534 entries, 0 to 533
Data columns (total 10 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 EDUCATION 534 non-null int64
1 SOUTH 534 non-null category
2 SEX 534 non-null category
3 EXPERIENCE 534 non-null int64
4 UNION 534 non-null category
5 AGE 534 non-null int64
6 RACE 534 non-null category
7 OCCUPATION 534 non-null category
8 SECTOR 534 non-null category
9 MARR 534 non-null category
dtypes: category(7), int64(3)
memory usage: 17.2 KB
如前所述，数据集包含不同数据类型的列，我们需要为每种数据类型应用特定的预处理。特别是，如果分类变量没有首先编码为整数，则不能包含在线性模型中。此外，为了避免将分类特征视为有序值，我们需要对其进行独热编码。我们的预处理器将
对分类列进行独热编码（即按类别生成一列），仅针对非二元分类变量；
作为第一种方法（我们稍后将看到数值的标准化将如何影响我们的讨论），保持数值不变。
from sklearn.compose import make_column_transformer
from sklearn.preprocessing import OneHotEncoder
categorical_columns = ["RACE", "OCCUPATION", "SECTOR", "MARR", "UNION", "SEX", "SOUTH"]
numerical_columns = ["EDUCATION", "EXPERIENCE", "AGE"]
preprocessor = make_column_transformer(
(OneHotEncoder(drop="if_binary"), categorical_columns),
remainder="passthrough",
verbose_feature_names_out=False, # avoid to prepend the preprocessor names
)
为了将数据集描述为线性模型，我们使用具有非常小的正则化的岭回归器来模拟工资的对数。
from sklearn.compose import TransformedTargetRegressor
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
model = make_pipeline(
preprocessor,
TransformedTargetRegressor(
regressor=Ridge(alpha=1e-10), func=np.log10, inverse_func=sp.special.exp10
),
)
处理数据集#
首先，我们拟合模型。
model.fit(X_train, y_train)
然后我们检查计算模型的性能，将其预测绘制在测试集上，并计算例如模型的中位数绝对误差。
from sklearn.metrics import PredictionErrorDisplay, median_absolute_error
mae_train = median_absolute_error(y_train, model.predict(X_train))
y_pred = model.predict(X_test)
mae_test = median_absolute_error(y_test, y_pred)
scores = {
"MedAE on training set": f"{mae_train:.2f} $/hour",
"MedAE on testing set": f"{mae_test:.2f} $/hour",
}
_, ax = plt.subplots(figsize=(5, 5))
display = PredictionErrorDisplay.from_predictions(
y_test, y_pred, kind="actual_vs_predicted", ax=ax, scatter_kwargs={"alpha": 0.5}
)
ax.set_title("Ridge model, small regularization")
for name, score in scores.items():
ax.plot([], [], " ", label=f"{name}: {score}")
ax.legend(loc="upper left")
plt.tight_layout()
学习到的模型远非一个能够进行准确预测的好模型：从上图可以看出这一点，在该图中，好的预测应该位于虚线线上。
在下一节中，我们将解释模型的系数。在这样做的同时，我们应该记住，我们得出的任何结论都是关于我们构建的模型，而不是关于数据的真实（现实世界）生成过程。
解释系数：规模很重要#
首先，我们可以查看我们已拟合的回归器的系数值。
feature_names = model[:-1].get_feature_names_out()
coefs = pd.DataFrame(
model[-1].regressor_.coef_,
columns=["Coefficients"],
index=feature_names,
)
coefs
AGE 系数以“美元/小时/年”表示，而 EDUCATION 系数以“美元/小时/年教育”表示。这种系数表示方法的好处在于可以清楚地说明模型的实际预测：AGE 增加\(1\)年意味着每小时工资减少\(0.030867\)美元，而 EDUCATION 增加\(1\)年意味着每小时工资增加\(0.054699\)美元。另一方面，分类变量（如 UNION 或 SEX）是无量纲数，取值要么为 0 要么为 1。它们的系数以美元/小时表示。因此，由于特征具有不同的自然尺度，因此值范围不同，因为它们的计量单位不同，我们无法比较不同系数的大小。如果我们绘制系数图，这一点就更明显了。
coefs.plot.barh(figsize=(9, 7))
plt.title("Ridge model, small regularization")
plt.axvline(x=0, color=".5")
plt.xlabel("Raw coefficient values")
plt.subplots_adjust(left=0.3)
事实上，从上图可以看出，决定 WAGE 的最重要因素似乎是 UNION 变量，即使我们的直觉可能会告诉我们，像 EXPERIENCE 这样的变量应该有更大的影响。
查看系数图来衡量特征重要性可能会产生误导，因为其中一些特征的变化范围很小，而另一些特征，如 AGE，变化范围则大得多，跨越几十年。
如果我们比较不同特征的标准差，这一点就显而易见了。
X_train_preprocessed = pd.DataFrame(
model[:-1].transform(X_train), columns=feature_names
)
X_train_preprocessed.std(axis=0).plot.barh(figsize=(9, 7))
plt.title("Feature ranges")
plt.xlabel("Std. dev. of feature values")
plt.subplots_adjust(left=0.3)
将系数乘以相关特征的标准差，可以将所有系数都转换为相同的计量单位。正如我们将在之后看到的，这等效于将数值变量归一化为其标准差，如\(y = \sum{coef_i \times X_i} = \sum{(coef_i \times std_i) \times (X_i / std_i)}\)。
这样，我们强调，在其他条件相同的情况下，特征的方差越大，相应系数对输出的权重就越大。
coefs = pd.DataFrame(
model[-1].regressor_.coef_ * X_train_preprocessed.std(axis=0),
columns=["Coefficient importance"],
index=feature_names,
)
coefs.plot(kind="barh", figsize=(9, 7))
plt.xlabel("Coefficient values corrected by the feature's std. dev.")
plt.title("Ridge model, small regularization")
plt.axvline(x=0, color=".5")
plt.subplots_adjust(left=0.3)
现在系数已经进行了缩放，我们可以安全地进行比较了。
警告
为什么上图表明年龄增加会导致工资减少？为什么初始 pairplot显示的结果却相反？
上图告诉我们的是，在所有其他特征保持不变的情况下，特定特征与目标变量之间的依赖关系，即**条件依赖关系**。在所有其他特征保持不变的情况下，AGE 的增加会导致 WAGE 的减少。相反，在所有其他特征保持不变的情况下，EXPERIENCE 的增加会导致 WAGE 的增加。此外，AGE、EXPERIENCE 和 EDUCATION 是三个最影响模型的变量。
解释系数：谨慎对待因果关系#
线性模型是衡量统计关联的强大工具，但在进行因果关系陈述时应谨慎，毕竟相关性并不总是意味着因果关系。这在社会科学中尤其困难，因为我们观察到的变量只是潜在因果过程的代理变量。
在我们的特定案例中，我们可以将个人的 EDUCATION 视为其专业能力的代理变量，这是我们感兴趣但无法观察到的真实变量。我们当然希望认为在学校学习的时间越长，技术能力就会越高，但因果关系也可能相反。也就是说，那些技术能力强的人往往会在学校学习更长时间。
雇主不太可能关心哪种情况（或者是否是两者兼而有之），只要他们确信受教育程度更高的人更适合这份工作，他们就会很乐意支付更高的 WAGE。
这种混杂效应在考虑某种干预措施时会成为问题，例如政府对大学学位的补贴或鼓励个人接受高等教育的宣传材料。这些措施的效用最终可能会被夸大，特别是如果混杂程度很强的话。我们的模型预测，每增加一年的教育，每小时工资就会增加\(0.054699\)美元。由于这种混杂效应，实际的因果效应可能会更低。
检查系数的可变性#
我们可以通过交叉验证来检查系数的可变性：这是一种数据扰动形式（与重采样相关）。
如果在改变输入数据集时系数变化很大，则它们的稳健性无法保证，并且可能应该谨慎解释。
from sklearn.model_selection import RepeatedKFold, cross_validate
cv = RepeatedKFold(n_splits=5, n_repeats=5, random_state=0)
cv_model = cross_validate(
model,
X,
y,
cv=cv,
return_estimator=True,
n_jobs=2,
)
coefs = pd.DataFrame(
[
est[-1].regressor_.coef_ * est[:-1].transform(X.iloc[train_idx]).std(axis=0)
for est, (train_idx, _) in zip(cv_model["estimator"], cv.split(X, y))
],
columns=feature_names,
)
plt.figure(figsize=(9, 7))
sns.stripplot(data=coefs, orient="h", palette="dark:k", alpha=0.5)
sns.boxplot(data=coefs, orient="h", color="cyan", saturation=0.5, whis=10)
plt.axvline(x=0, color=".5")
plt.xlabel("Coefficient importance")
plt.title("Coefficient importance and its variability")
plt.suptitle("Ridge model, small regularization")
plt.subplots_adjust(left=0.3)
相关变量的问题#
AGE 和 EXPERIENCE 系数受到强烈可变性的影响，这可能是由于这两个特征之间的共线性造成的：由于 AGE 和 EXPERIENCE 在数据中一起变化，因此很难区分它们的影响。
为了验证这种解释，我们绘制了 AGE 和 EXPERIENCE 系数的可变性图。
plt.ylabel("Age coefficient")
plt.xlabel("Experience coefficient")
plt.grid(True)
plt.xlim(-0.4, 0.5)
plt.ylim(-0.4, 0.5)
plt.scatter(coefs["AGE"], coefs["EXPERIENCE"])
_ = plt.title("Co-variations of coefficients for AGE and EXPERIENCE across folds")
有两个区域的数据点比较密集：当 EXPERIENCE 系数为正时，AGE 系数为负，反之亦然。
为了进一步研究，我们删除了这两个特征中的一个，并检查了对模型稳定性的影响。
column_to_drop = ["AGE"]
cv_model = cross_validate(
model,
X.drop(columns=column_to_drop),
y,
cv=cv,
return_estimator=True,
n_jobs=2,
)
coefs = pd.DataFrame(
[
est[-1].regressor_.coef_
* est[:-1].transform(X.drop(columns=column_to_drop).iloc[train_idx]).std(axis=0)
for est, (train_idx, _) in zip(cv_model["estimator"], cv.split(X, y))
],
columns=feature_names[:-1],
)
plt.figure(figsize=(9, 7))
sns.stripplot(data=coefs, orient="h", palette="dark:k", alpha=0.5)
sns.boxplot(data=coefs, orient="h", color="cyan", saturation=0.5)
plt.axvline(x=0, color=".5")
plt.title("Coefficient importance and its variability")
plt.xlabel("Coefficient importance")
plt.suptitle("Ridge model, small regularization, AGE dropped")
plt.subplots_adjust(left=0.3)
EXPERIENCE 系数的估计现在显示出可变性大大降低。在交叉验证期间训练的所有模型中，EXPERIENCE 仍然很重要。
数值变量预处理#
如上所述（参见“机器学习管道”），我们也可以选择在训练模型之前缩放数值。当我们在岭回归中对所有变量应用大致相同的正则化时，这很有用。预处理器被重新定义为减去均值并将变量缩放为单位方差。
from sklearn.preprocessing import StandardScaler
preprocessor = make_column_transformer(
(OneHotEncoder(drop="if_binary"), categorical_columns),
(StandardScaler(), numerical_columns),
)
模型将保持不变。
model = make_pipeline(
preprocessor,
TransformedTargetRegressor(
regressor=Ridge(alpha=1e-10), func=np.log10, inverse_func=sp.special.exp10
),
)
model.fit(X_train, y_train)
同样，我们使用例如模型的中位数绝对误差和 R 平方系数来检查计算模型的性能。
mae_train = median_absolute_error(y_train, model.predict(X_train))
y_pred = model.predict(X_test)
mae_test = median_absolute_error(y_test, y_pred)
scores = {
"MedAE on training set": f"{mae_train:.2f} $/hour",
"MedAE on testing set": f"{mae_test:.2f} $/hour",
}
_, ax = plt.subplots(figsize=(5, 5))
display = PredictionErrorDisplay.from_predictions(
y_test, y_pred, kind="actual_vs_predicted", ax=ax, scatter_kwargs={"alpha": 0.5}
)
ax.set_title("Ridge model, small regularization")
for name, score in scores.items():
ax.plot([], [], " ", label=f"{name}: {score}")
ax.legend(loc="upper left")
plt.tight_layout()
对于系数分析，这次不需要缩放，因为缩放是在预处理步骤中进行的。
coefs = pd.DataFrame(
model[-1].regressor_.coef_,
columns=["Coefficients importance"],
index=feature_names,
)
coefs.plot.barh(figsize=(9, 7))
plt.title("Ridge model, small regularization, normalized variables")
plt.xlabel("Raw coefficient values")
plt.axvline(x=0, color=".5")
plt.subplots_adjust(left=0.3)
我们现在检查多个交叉验证折中的系数。与上面的例子一样，我们不需要按特征值的标准差缩放系数，因为这种缩放已经在管道的预处理步骤中完成了。
cv_model = cross_validate(
model,
X,
y,
cv=cv,
return_estimator=True,
n_jobs=2,
)
coefs = pd.DataFrame(
[est[-1].regressor_.coef_ for est in cv_model["estimator"]], columns=feature_names
)
plt.figure(figsize=(9, 7))
sns.stripplot(data=coefs, orient="h", palette="dark:k", alpha=0.5)
sns.boxplot(data=coefs, orient="h", color="cyan", saturation=0.5, whis=10)
plt.axvline(x=0, color=".5")
plt.title("Coefficient variability")
plt.subplots_adjust(left=0.3)
结果与非归一化情况非常相似。
带有正则化的线性模型#
在机器学习实践中，岭回归更常与非可忽略的正则化一起使用。
上面，我们将这种正则化限制在一个非常小的范围内。正则化提高了问题的条件数并降低了估计值的方差。RidgeCV 应用交叉验证以确定哪个正则化参数值（alpha）最适合预测。
from sklearn.linear_model import RidgeCV
alphas = np.logspace(-10, 10, 21) # alpha values to be chosen from by cross-validation
model = make_pipeline(
preprocessor,
TransformedTargetRegressor(
regressor=RidgeCV(alphas=alphas),
func=np.log10,
inverse_func=sp.special.exp10,
),
)
model.fit(X_train, y_train)
首先我们检查已选择哪个\(\alpha\)值。
model[-1].regressor_.alpha_
np.float64(10.0)
然后我们检查预测的质量。
mae_train = median_absolute_error(y_train, model.predict(X_train))
y_pred = model.predict(X_test)
mae_test = median_absolute_error(y_test, y_pred)
scores = {
"MedAE on training set": f"{mae_train:.2f} $/hour",
"MedAE on testing set": f"{mae_test:.2f} $/hour",
}
_, ax = plt.subplots(figsize=(5, 5))
display = PredictionErrorDisplay.from_predictions(
y_test, y_pred, kind="actual_vs_predicted", ax=ax, scatter_kwargs={"alpha": 0.5}
)
ax.set_title("Ridge model, optimum regularization")
for name, score in scores.items():
ax.plot([], [], " ", label=f"{name}: {score}")
ax.legend(loc="upper left")
plt.tight_layout()
正则化模型再现数据的 kemampuan 与非正则化模型相似。
coefs = pd.DataFrame(
model[-1].regressor_.coef_,
columns=["Coefficients importance"],
index=feature_names,
)
coefs.plot.barh(figsize=(9, 7))
plt.title("Ridge model, with regularization, normalized variables")
plt.xlabel("Raw coefficient values")
plt.axvline(x=0, color=".5")
plt.subplots_adjust(left=0.3)
系数存在显著差异。“年龄”和“经验”系数均为正值，但它们对预测的影响较小了。
正则化降低了相关变量对模型的影响，因为权重在两个预测变量之间共享，因此任何一个变量单独都不会具有很强的权重。
另一方面，使用正则化获得的权重更加稳定（参见Ridge回归和分类用户指南部分）。通过交叉验证中数据扰动获得的图可以看出这种稳定性的提高。该图可以与之前的图进行比较。
cv_model = cross_validate(
model,
X,
y,
cv=cv,
return_estimator=True,
n_jobs=2,
)
coefs = pd.DataFrame(
[est[-1].regressor_.coef_ for est in cv_model["estimator"]], columns=feature_names
)
plt.ylabel("Age coefficient")
plt.xlabel("Experience coefficient")
plt.grid(True)
plt.xlim(-0.4, 0.5)
plt.ylim(-0.4, 0.5)
plt.scatter(coefs["AGE"], coefs["EXPERIENCE"])
_ = plt.title("Co-variations of coefficients for AGE and EXPERIENCE across folds")
具有稀疏系数的线性模型#
考虑数据集中的相关变量的另一种方法是估计稀疏系数。在某种程度上，当我们在之前的Ridge估计中删除“年龄”列时，我们已经手动进行了这种操作。
Lasso模型（参见Lasso用户指南部分）估计稀疏系数。LassoCV应用交叉验证来确定哪个正则化参数值（alpha）最适合模型估计。
from sklearn.linear_model import LassoCV
alphas = np.logspace(-10, 10, 21) # alpha values to be chosen from by cross-validation
model = make_pipeline(
preprocessor,
TransformedTargetRegressor(
regressor=LassoCV(alphas=alphas, max_iter=100_000),
func=np.log10,
inverse_func=sp.special.exp10,
),
)
_ = model.fit(X_train, y_train)
首先，我们验证选择哪个\(\alpha\)值。
model[-1].regressor_.alpha_
np.float64(0.001)
然后我们检查预测的质量。
mae_train = median_absolute_error(y_train, model.predict(X_train))
y_pred = model.predict(X_test)
mae_test = median_absolute_error(y_test, y_pred)
scores = {
"MedAE on training set": f"{mae_train:.2f} $/hour",
"MedAE on testing set": f"{mae_test:.2f} $/hour",
}
_, ax = plt.subplots(figsize=(6, 6))
display = PredictionErrorDisplay.from_predictions(
y_test, y_pred, kind="actual_vs_predicted", ax=ax, scatter_kwargs={"alpha": 0.5}
)
ax.set_title("Lasso model, optimum regularization")
for name, score in scores.items():
ax.plot([], [], " ", label=f"{name}: {score}")
ax.legend(loc="upper left")
plt.tight_layout()
对于我们的数据集，模型的预测能力仍然不强。
coefs = pd.DataFrame(
model[-1].regressor_.coef_,
columns=["Coefficients importance"],
index=feature_names,
)
coefs.plot(kind="barh", figsize=(9, 7))
plt.title("Lasso model, optimum regularization, normalized variables")
plt.axvline(x=0, color=".5")
plt.subplots_adjust(left=0.3)
Lasso模型识别“年龄”和“经验”之间的相关性，并为了预测而抑制其中一个。
需要注意的是，被舍弃的系数本身可能仍然与结果相关：模型选择舍弃它们是因为它们相对于其他特征几乎没有提供额外的信息。此外，对于相关特征，这种选择是不稳定的，应谨慎解释。
事实上，我们可以检查各折中系数的变异性。
cv_model = cross_validate(
model,
X,
y,
cv=cv,
return_estimator=True,
n_jobs=2,
)
coefs = pd.DataFrame(
[est[-1].regressor_.coef_ for est in cv_model["estimator"]], columns=feature_names
)
plt.figure(figsize=(9, 7))
sns.stripplot(data=coefs, orient="h", palette="dark:k", alpha=0.5)
sns.boxplot(data=coefs, orient="h", color="cyan", saturation=0.5, whis=100)
plt.axvline(x=0, color=".5")
plt.title("Coefficient variability")
plt.subplots_adjust(left=0.3)
我们观察到，“年龄”和“经验”系数根据不同的折而变化很大。
错误的因果解释#
政策制定者可能想知道教育对工资的影响，以评估旨在鼓励人们接受更多教育的某种政策是否具有经济意义。虽然机器学习模型非常擅长测量统计关联，但它们通常无法推断因果关系。
我们可能很想查看我们最后一个模型（或任何模型）中教育对工资的系数，并得出结论认为它捕捉到了标准化教育变量变化对工资的真实影响。
不幸的是，可能存在未观察到的混杂变量，这些变量会夸大或缩小该系数。混杂变量是指同时导致教育和工资的变量。此类变量的一个例子是能力。据推测，能力较强的人更有可能接受教育，同时在任何教育水平上都更有可能获得更高的每小时工资。在这种情况下，能力会在教育系数上产生正的遗漏变量偏差（OVB），从而夸大教育对工资的影响。
参见机器学习无法推断因果关系，了解能力OVB的模拟案例。
经验教训#
必须将系数缩放至相同的度量单位才能检索特征重要性。使用特征的标准差对其进行缩放是一个有用的近似值。
当存在混杂效应时，解释因果关系很困难。如果两个变量之间的关系也受某些未观察到的因素影响，那么在得出关于因果关系的结论时，我们应该谨慎。
多元线性模型中的系数表示给定特征与目标之间的依赖关系，**以其他特征为条件**。
相关特征会导致线性模型系数的不稳定性，并且它们的影响难以很好地区分开来。
不同的线性模型对特征相关性的响应不同，系数可能彼此之间存在显著差异。
检查交叉验证循环各折中的系数，可以了解其稳定性。
系数不太可能具有任何因果意义。它们往往会受到未观察到的混杂因素的影响。
检查工具不一定能提供关于真实数据生成过程的见解。
脚本总运行时间：（0 分钟 13.554 秒）
相关示例
