线性模型系数解释中的常见陷阱#
在線性模型中，目標值被建模為特徵的線性組合（有關 scikit-learn 中可用的一組線性模型的描述，請參閱線性模型用戶指南部分）。多元線性模型中的係數表示給定特徵\(X_i\)與目標\(y\)之間的關係，假設所有其他特徵保持不變（條件依賴性）。這與繪製\(X_i\)與\(y\)並擬合線性關係不同：在這種情況下，所有其他特徵的可能值都將在估計中考慮在內（邊緣依賴性）。
本示例將提供一些關於解釋線性模型中係數的提示，指出當線性模型不適合描述數據集，或者當特徵相關時出現的問題。
注意
請記住，特徵\(X\)和結果\(y\)通常是我們未知的數據生成過程的結果。機器學習模型經過訓練，可以從樣本數據中逼近將\(X\)與\(y\)聯繫起來的未觀察到的數學函數。因此，對模型做出的任何解釋不一定能推廣到真實的數據生成過程。當模型質量較差或樣本數據不能代表總體時，尤其如此。
我們將使用來自“當前人口調查”的數據，從 1985 年開始，根據經驗、年齡或教育等各種特徵來預測工資。
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 DataFrame 的形式檢索數據。
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")
仔細觀察 WAGE 分佈表明它有一個長尾。出於這個原因，我們應該取其對數，以將其近似轉換為正態分佈（嶺迴歸或套索等線性模型最適合誤差的正態分佈）。
當 EDUCATION 增加時，WAGE 正在增加。請注意，這裡表示的 WAGE 和 EDUCATION 之間的依賴關係是邊緣依賴關係，即它描述了特定變量的行為，而沒有保持其他變量固定。
此外，EXPERIENCE 和 AGE 呈高度線性相關。
機器學習管道#
為了設計我們的機器學習管道，我們首先手動檢查我們正在處理的數據類型
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
)
為了將數據集描述為線性模型，我們使用帶有非常小的正則化的嶺迴歸器來模擬 WAGE 的對數。
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)
現在係數已經縮放，我們可以安全地比較它們。
警告
為什麼上面的圖表表明年齡增加導致工資減少？為什麼初始對圖顯示相反的結果？
上面的圖表告訴我們，當所有其他特徵保持不變時，特定特徵與目標之間的依賴關係，即**條件依賴關係**。當所有其他特徵保持不變時，AGE 的增加將導致 WAGE 的減少。相反，當所有其他特徵保持不變時，EXPERIENCE 的增加將導致 WAGE 的增加。此外，AGE、EXPERIENCE 和 EDUCATION 是影響模型的三個最重要的變量。
解釋係數：謹慎看待因果關係#
线性模型是衡量统计关联性的一个好工具，但我们在进行因果关系陈述时应该谨慎，毕竟相关性并不总是意味着因果关系。这在社会科学中尤其困难，因为我们观察到的变量只是作为潜在因果过程的代理。
在我们的特定情况下，我们可以将个人的教育视为其专业能力的代理，而专业能力是我们感兴趣但无法观察到的真实变量。我们当然希望认为，在学校学习的时间越长，技术能力就会越强，但也有可能因果关系是相反的。也就是说，那些技术能力强的人往往在学校学习的时间更长。
雇主不太可能关心是哪种情况（或者是否是两者兼而有之），只要他们仍然相信受教育程度更高的人更适合这份工作，他们就会很乐意支付更高的工资。
当考虑某种干预措施时，这种混杂效应就会变得有问题，例如政府对大学学位的补贴或鼓励个人接受高等教育的宣传材料。这些措施的效用最终可能会被夸大，尤其是在混杂程度很强的情况下。我们的模型预测，每增加一年的教育，每小时工资将增加 \(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)
相关变量的问题#
年龄和经验系数受到强变异性的影响，这可能是由于这两个特征之间的共线性：由于年龄和经验在数据中一起变化，因此很难将它们的影响区分开来。
为了验证这种解释，我们绘制了年龄和经验系数的变异性。
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")
有两个区域被填充：当经验系数为正时，年龄系数为负，反之亦然。
为了更进一步，我们删除了这两个特征中的一个，并检查对模型稳定性的影响。
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)
经验系数的估计现在显示出变异性大大降低。经验对于交叉验证期间训练的所有模型仍然很重要。
预处理数值变量#
如上所述（参见“机器学习管道”），我们也可以选择在训练模型之前对数值进行缩放。当我们在岭回归中对所有数值应用相同程度的正则化时，这很有用。预处理器被重新定义，以便减去均值并将变量缩放到单位方差。
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_
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()
正则化模型再现数据的性能与非正则化模型相似。
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)
系数有显著差异。年龄和经验系数均为正，但它们现在对预测的影响较小。
正则化减少了相关变量对模型的影响，因为权重在两个预测变量之间共享，因此单独的权重不会很大。
另一方面，使用正则化获得的权重更稳定（参见 岭回归和分类 用户指南部分）。这种增加的稳定性从交叉验证中获得的数据扰动图中可以看出。该图可以与 之前的图 进行比较。
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")
具有稀疏系数的线性模型#
考虑数据集中相关变量的另一种方法是估计稀疏系数。在某种程度上，我们在之前手动删除年龄列时已经做到了这一点。
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_
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 分钟 36.093 秒）
相关示例
