复现 Nature 图表可视化——基于模型残差分析与显著性检验的模型解释
2025-01-06
背景
在数据科学和机器学习领域，模型的评估和结果的可视化是至关重要的环节，残差分析通过衡量预测值与真实值之间的差异，为我们提供了深入了解模型性能的关键手段，Nature文章中的b图展示了一种有效的残差可视化方式，不仅能够直观呈现预测误差，还通过对不同数据集的对比揭示模型的优缺点
接下来作者将尝试复现该残差可视化图表，并扩展其应用，让我们得以更加清晰地分析不同模型和数据集之间的预测效果差异，为模型优化提供了进一步的依据
代码实现
数据导入及筛选
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings("ignore")
df = pd.read_excel('实验数据.xlsx')
df_1 = df[df['AgeBin'] == '年轻'].reset_index(drop=True)
df_2 = df[df['AgeBin'] == '中年'].reset_index(drop=True)
df_3 = df[df['AgeBin'] == '老年'].reset_index(drop=True)
df
读取一个Excel文件中的数据，并根据特定的条件对数据进行筛选和分类，对不同的数据建立多个模型方便对该可视化的复现（后文会针对单一模型进行类似绘图）
数据集划分及模型建立
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
# 划分特征和目标变量
X_1 = df_1.drop(['price', 'AgeBin'], axis=1)
y_1 = df_1['price']
X_2 = df_2.drop(['price', 'AgeBin'], axis=1)
y_2 = df_2['price']
X_3 = df_3.drop(['price', 'AgeBin'], axis=1)
y_3 = df_3['price']
# 对三个数据集分别进行训练集、测试集划分
X_train_1, X_test_1, y_train_1, y_test_1 = train_test_split(X_1, y_1, test_size=0.3, random_state=42)
X_train_2, X_test_2, y_train_2, y_test_2 = train_test_split(X_2, y_2, test_size=0.3, random_state=42)
X_train_3, X_test_3, y_train_3, y_test_3 = train_test_split(X_3, y_3, test_size=0.3, random_state=42)
# 创建随机森林回归器实例
rf_regressor_1 = RandomForestRegressor(
n_estimators=100,
criterion='squared_error', # 回归模型的评估标准
max_depth=None,
min_samples_split=2,
min_samples_leaf=1,
min_weight_fraction_leaf=0.0,
random_state=42,
max_leaf_nodes=None,
min_impurity_decrease=0.0
)
rf_regressor_2 = RandomForestRegressor(
n_estimators=100,
criterion='squared_error',
max_depth=None,
min_samples_split=2,
min_samples_leaf=1,
min_weight_fraction_leaf=0.0,
random_state=42,
max_leaf_nodes=None,
min_impurity_decrease=0.0
)
rf_regressor_3 = RandomForestRegressor(
n_estimators=100,
criterion='squared_error',
max_depth=None,
min_samples_split=2,
min_samples_leaf=1,
min_weight_fraction_leaf=0.0,
random_state=42,
max_leaf_nodes=None,
min_impurity_decrease=0.0
)
# 分别训练三个随机森林回归模型
rf_regressor_1.fit(X_train_1, y_train_1)
rf_regressor_2.fit(X_train_2, y_train_2)
rf_regressor_3.fit(X_train_3, y_train_3)对三个不同年龄段的数据集分别创建和训练了随机森林回归模型，用于预测目标变量 price
模型评价指标输出
from sklearn import metrics
def evaluate_model(model, X_train, y_train, X_test, y_test, model_name):
# 预测
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)
# 计算训练集的指标
mse_train = metrics.mean_squared_error(y_train, y_pred_train)
rmse_train = np.sqrt(mse_train)
mae_train = metrics.mean_absolute_error(y_train, y_pred_train)
r2_train = metrics.r2_score(y_train, y_pred_train)
# 计算测试集的指标
mse_test = metrics.mean_squared_error(y_test, y_pred_test)
rmse_test = np.sqrt(mse_test)
mae_test = metrics.mean_absolute_error(y_test, y_pred_test)
r2_test = metrics.r2_score(y_test, y_pred_test)
# 输出结果
print(f"{model_name} 训练集评价指标:")
print("均方误差 (MSE):", mse_train)
print("均方根误差 (RMSE):", rmse_train)
print("平均绝对误差 (MAE):", mae_train)
print("拟合优度 (R-squared):", r2_train)
print("\n")
print(f"{model_name} 测试集评价指标:")
print("均方误差 (MSE):", mse_test)
print("均方根误差 (RMSE):", rmse_test)
print("平均绝对误差 (MAE):", mae_test)
print("拟合优度 (R-squared):", r2_test)
print("\n" + "-"*40 + "\n")
# 分别计算三个模型的指标
evaluate_model(rf_regressor_1, X_train_1, y_train_1, X_test_1, y_test_1, "模型 1 (年轻)")
evaluate_model(rf_regressor_2, X_train_2, y_train_2, X_test_2, y_test_2, "模型 2 (中年)")
evaluate_model(rf_regressor_3, X_train_3, y_train_3, X_test_3, y_test_3, "模型 3 (老年)")
输出模型各数据集的评价指标
模型残差计算
def calculate_residuals(model, X_test, y_test, model_name):
# 预测
y_pred_test = model.predict(X_test)
# 计算残差
residuals = y_test - y_pred_test
# 记录残差到 DataFrame
residuals_df = pd.DataFrame({
'True Value': y_test,
'Predicted Value': y_pred_test,
'Residual': residuals
})
residuals_df['Model'] = model_name
return residuals_df
# 分别计算三个模型的残差
residuals_1 = calculate_residuals(rf_regressor_1, X_test_1, y_test_1, "Model 1 (Young)")
residuals_2 = calculate_residuals(rf_regressor_2, X_test_2, y_test_2, "Model 2 (Middle Age)")
residuals_3 = calculate_residuals(rf_regressor_3, X_test_3, y_test_3, "Model 3 (Old)")
# 合并三个模型的残差数据
residuals_all = pd.concat([residuals_1, residuals_2, residuals_3], ignore_index=True)
residuals_all
计算三个随机森林回归模型在测试集上的预测残差，并将这些残差合并为一个包含所有模型的 DataFrame
残差 Mann-Whitney U 检验
from scipy.stats import mannwhitneyu
# 获取模型的唯一组合
models = residuals_all['Model'].unique()
comparisons = [(models[i], models[j]) for i in range(len(models)) for j in range(i+1, len(models))]
# 存储每个比较的结果
results = []
# 进行 Mann-Whitney U 检验
for group1, group2 in comparisons:
group1_data = residuals_all[residuals_all['Model'] == group1]['Residual']
group2_data = residuals_all[residuals_all['Model'] == group2]['Residual']
# 进行 Mann-Whitney U 检验
stat, p_value = mannwhitneyu(group1_data, group2_data, alternative='two-sided')
# 存储结果
results.append([group1, group2, p_value])
# 将结果转换为 DataFrame
mannwhitney_results = pd.DataFrame(results, columns=['group1', 'group2', 'p_value'])
# 输出结果
mannwhitney_results
对三个模型的残差进行两两比较，使用 Mann-Whitney U 检验判断它们的残差分布是否存在显著差异，并将每次比较的结果（包括 p 值）汇总为一个 DataFrame，其中p 值用于比较两个模型的残差分布差异，具体而言，它表示两个模型的预测误差（残差）是否来自相同的分布，通常情况下，将 p 值与显著性水平（例如 0.05）进行比较：若 p 值小于 0.05，则可以认为这两个模型的残差分布存在显著差异；若 p 值大于或等于 0.05，则无法拒绝两个模型残差分布相同的假设（即没有显著差异），显著性差异存在表示两个模型的预测误差分布不同，模型性能表现有差异；显著性差异不存在表示误差分布相似，模型性能表现相近（所以作者认为也可以用在单模型上针对训练集和测试集去做来判别是否存在显著性差异，进一步说明模型是否过拟合，后文会给出该可视化代码）
不同模型预测残差的可视化（不带p值）
import seaborn as sns
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
palette = ["C20", "#A184BC", "#ff7f0e"]
plt.figure(figsize=(12, 10), dpi=800)
# 创建 JointGrid 对象
g = sns.JointGrid(data=residuals_all, x="True Value", y="Predicted Value", height=10)
# 在主图中绘制散点图，使用指定的配色
g.plot_joint(sns.scatterplot, hue='Model', data=residuals_all, palette=palette, alpha=0.6)
# 分别绘制三个模型的堆叠边缘密度图，不显示图例
sns.kdeplot(data=residuals_all, x='True Value', hue='Model', ax=g.ax_marg_x, fill=True, common_norm=False, palette=palette, alpha=0.5, legend=False)
sns.kdeplot(data=residuals_all, y='Predicted Value', hue='Model', ax=g.ax_marg_y, fill=True, common_norm=False, palette=palette, alpha=0.5, legend=False)
# 添加 y=x 对角线
g.ax_joint.plot([residuals_all['True Value'].min(), residuals_all['True Value'].max()],
[residuals_all['True Value'].min(), residuals_all['True Value'].max()],
color='red', linestyle='--', label='y=x')
# 设置标签和标题
g.set_axis_labels('True Values', 'Predicted Values', fontsize=14)
# 仅在主图中显示图例
g.ax_joint.legend(title='Model', loc='upper left')
# 在右下角添加箱线图，调整位置以避免重叠
ax_inset = inset_axes(g.ax_joint, width="40%", height="20%", loc='lower right',
bbox_to_anchor=(0.2, 0.05, 0.8, 0.8), # 调整这里的坐标来平移
bbox_transform=g.ax_joint.transAxes)
# 绘制横向箱线图
sns.boxplot(data=residuals_all, y='Model', x='Residual', palette=palette, ax=ax_inset)
# 设置箱线图的标题和标签
ax_inset.set_title('Residuals', fontsize=10)
ax_inset.set_xlabel('')
ax_inset.set_ylabel('')
# 关闭箱线图 y 轴显示
ax_inset.yaxis.set_visible(False)
plt.savefig("表1.pdf", format='pdf', bbox_inches='tight')
# 显示图像
plt.show()
不同模型预测残差的可视化（带p值）
# 定义整个可视化的配色方案
palette = {"Model 1 (Young)": "#1f77b4", "Model 2 (Middle Age)": "#9467bd", "Model 3 (Old)": "#ff7f0e"}
plt.figure(figsize=(12, 10), dpi=800)
# 创建 JointGrid 对象
g = sns.JointGrid(data=residuals_all, x="True Value", y="Predicted Value", height=10)
# 在主图中绘制散点图，使用指定的配色
g.plot_joint(sns.scatterplot, hue='Model', data=residuals_all, palette=palette, alpha=0.6)
# 分别绘制三个模型的堆叠边缘密度图，不显示图例
sns.kdeplot(data=residuals_all, x='True Value', hue='Model', ax=g.ax_marg_x, fill=True, common_norm=False, palette=palette.values(), alpha=0.5, legend=False)
sns.kdeplot(data=residuals_all, y='Predicted Value', hue='Model', ax=g.ax_marg_y, fill=True, common_norm=False, palette=palette.values(), alpha=0.5, legend=False)
# 添加 y=x 对角线
g.ax_joint.plot([residuals_all['True Value'].min(), residuals_all['True Value'].max()],
[residuals_all['True Value'].min(), residuals_all['True Value'].max()],
color='red', linestyle='--', label='y=x')
# 设置标签和标题
g.set_axis_labels('True Values', 'Predicted Values', fontsize=14)
# 仅在主图中显示图例
g.ax_joint.legend(title='Model', loc='upper left')
# 在右下角添加箱线图，调整位置以避免重叠
ax_inset = inset_axes(g.ax_joint, width="40%", height="20%", loc='lower right',
bbox_to_anchor=(0.2, 0.05, 0.8, 0.8), # 调整这里的坐标来平移
bbox_transform=g.ax_joint.transAxes)
# 绘制横向箱线图
sns.boxplot(data=residuals_all, y='Model', x='Residual', palette=palette, ax=ax_inset)
# 设置箱线图的标题和标签
ax_inset.set_title('Residuals', fontsize=10)
ax_inset.set_xlabel('')
ax_inset.set_ylabel('')
# 关闭箱线图 y 轴显示
ax_inset.yaxis.set_visible(False)
# 准备文本内容并设置更紧凑的布局
y_offset = 0.2 # 调整整体垂直位置
x_offset_point1 = 0.44 # 第一个点的水平位置
x_offset_vs = 0.45 # 'vs' 文字的位置
x_offset_point2 = 0.47 # 第二个点的位置
x_offset_pvalue = 0.48 # p 值的位置
text_content = "" # 用于保存每一行的内容
for i, (group1, group2, p_value) in enumerate(results):
# 获取每个模型对应的颜色
color1 = palette[group1]
color2 = palette[group2]
# 在文本框中添加带颜色的点和 p 值
plt.gcf().text(x_offset_point1, y_offset, '●', fontsize=10, ha='right', va='top', color=color1)
plt.gcf().text(x_offset_vs, y_offset, 'vs', fontsize=8, ha='center', va='top')
plt.gcf().text(x_offset_point2, y_offset, '●', fontsize=10, ha='center', va='top', color=color2)
plt.gcf().text(x_offset_pvalue, y_offset, f'p = {p_value:.4f}', fontsize=8, ha='left', va='top')
y_offset -= 0.03 # 控制行间距，减小间距让内容更紧凑
plt.savefig("表2.pdf", format='pdf', bbox_inches='tight')
plt.show()
对比 Nature 文章中的 b 图，可以发现我们的 p 值可视化效果并不如原图那般美观。目前，作者暂时没有通过代码实现原图中这种 p 值可视化的更好方法。如果有读者能够实现，作者非常欢迎指教。鉴于此，作者认为可以采用上一节中不带 p 值的可视化结果，并通过 PPT 或 Photoshop 等软件手动添加 p 值标注，这种方式比通过代码实现更加便捷和灵活，在这里作者是通过在图中手动添加文本标注，使用彩色圆点和 p 值来表示不同模型残差之间的显著性差异
针对单一模型训练集、测试集可视化
def calculate_residuals(model, X, y, dataset_name):
# 预测
y_pred = model.predict(X)
# 计算残差
residuals = y - y_pred
# 记录残差到 DataFrame
residuals_df = pd.DataFrame({
'True Value': y,
'Predicted Value': y_pred,
'Residual': residuals
})
# 添加列标记是训练集或测试集
residuals_df['Dataset'] = dataset_name
return residuals_df
# 分别计算模型1在训练集和测试集的残差
residuals_train = calculate_residuals(rf_regressor_1, X_train_1, y_train_1, "Train")
residuals_test = calculate_residuals(rf_regressor_1, X_test_1, y_test_1, "Test")
# 合并训练集和测试集的残差数据
residuals_all = pd.concat([residuals_train, residuals_test], ignore_index=True)
# 定义整个可视化的配色方案
palette = ["#1f77b4", "#ff7f0e"] # 分别对应训练集和测试集
# 设置更大的画布
plt.figure(figsize=(12, 10), dpi=800)
# 创建 JointGrid 对象
g = sns.JointGrid(data=residuals_all, x="True Value", y="Predicted Value", height=10)
# 在主图中绘制散点图，区分训练集和测试集
g.plot_joint(sns.scatterplot, hue='Dataset', data=residuals_all, palette=palette, alpha=0.6)
# 分别绘制两个数据集的边缘密度图，不显示图例
sns.kdeplot(data=residuals_all, x='True Value', hue='Dataset', ax=g.ax_marg_x, fill=True, common_norm=False, palette=palette, alpha=0.5, legend=False)
sns.kdeplot(data=residuals_all, y='Predicted Value', hue='Dataset', ax=g.ax_marg_y, fill=True, common_norm=False, palette=palette, alpha=0.5, legend=False)
# 添加 y=x 对角线
g.ax_joint.plot([residuals_all['True Value'].min(), residuals_all['True Value'].max()],
[residuals_all['True Value'].min(), residuals_all['True Value'].max()],
color='red', linestyle='--', label='y=x')
# 设置标签和标题
g.set_axis_labels('True Values', 'Predicted Values', fontsize=14)
# 仅在主图中显示图例
g.ax_joint.legend(title='Dataset', loc='upper left')
# 在右下角添加箱线图，分别显示训练集和测试集的残差分布
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
ax_inset = inset_axes(g.ax_joint, width="40%", height="20%", loc='lower right',
bbox_to_anchor=(0.2, 0.05, 0.8, 0.8), # 调整这里的坐标来平移
bbox_transform=g.ax_joint.transAxes)
# 绘制横向箱线图，区分训练集和测试集的残差
sns.boxplot(data=residuals_all, y='Dataset', x='Residual', palette=palette, ax=ax_inset)
ax_inset.set_title('Residuals (Train vs Test)', fontsize=10)
ax_inset.set_xlabel('')
ax_inset.set_ylabel('')
ax_inset.yaxis.set_visible(False)
plt.savefig("表3.pdf", format='pdf', bbox_inches='tight')
plt.show()
这里可视化只展示了单个模型在训练集和测试集上的预测结果及其残差分布
# 分别获取训练集和测试集的残差
train_residuals = residuals_train['Residual']
test_residuals = residuals_test['Residual']
# 使用 Mann-Whitney U 检验来比较训练集和测试集的残差
stat, p_value = mannwhitneyu(train_residuals, test_residuals, alternative='two-sided')
print(f"U-statistic: {stat}")
print(f"p-value: {p_value}")
可以发现p 值非常小（远小于常用的显著性水平 0.05），这意味着训练集和测试集的残差分布之间存在显著差异，这表明模型在训练集和测试集上的表现存在明显不同，暗示模型可能存在过拟合问题，即模型在训练集上拟合得很好，但在测试集上表现较差，实际上这和前文该模型的拟合优度给出的结果一致，训练集上较高的 和较低的误差表明模型在训练集上拟合得很好，但在测试集上的较低 和较高误差表明模型无法很好地泛化，这与 Mann-Whitney U 检验中训练集和测试集残差分布的显著差异相一致，都认为模型存在过拟合的可能性
同话题下的热门内容
