机器学习——13种模型+统计显著性检验+贝叶斯优化+SHAP高级可视化（回归问题）
·
感谢关注，一起来学习干货。
本期摘要
1.EDA分析
2.绘制13种机器学习算法的性能对比
3.统计显著性检验
4.针对最优机器学习模型做SHAP解释分析
第一部分：探索性数据分析（EDA）
糖尿病数据集探索性数据分析（EDA）的Python代码。整个过程将被分解为十个阶段，每个阶段都会有详细的介绍和逐行代码注释，确保您能理解每一行代码的作用，并能将其应用到其他数据集中。
阶段 0：用户配置
作用与目的： 这个区域是整个脚本的“控制中心”。用户唯一需要修改的地方就在这里。通过将这些易变的参数（如文件路径、目标列名）集中放置在顶部，代码变得高度可复用。用户无需阅读和修改后续的复杂逻辑，只需在这里提供自己数据集的信息，即可运行整个分析流程。
excel_file_path = r'F:\2025-8-21-公众号Python机器学习ml.csv'# 定义数据文件的完整路径。r'' 表示原始字符串，可以防止路径中的反斜杠 \ 被误解为转义字符。
target_column_name = 'risk_score'# 定义您数据中作为目标变量（Y值，即需要预测的量）的列的名称。
sheet_name = None# 如果您的Excel文件有多个工作表(sheet)，请在这里指定要读取的表名（字符串格式，如 'Sheet1'）；如果文件只有一个工作表或您想读取第一个，则保持为 None。
阶段 1：准备工作 - 导入库
作用与目的： 这个阶段导入了所有后续分析和建模所需的Python库。一个好的编程习惯是把所有的导入语句都放在脚本的开头，这样可以清晰地看到代码依赖了哪些工具。
import pandas as pd # 导入pandas库，通常简写为pd，是数据分析的核心库，用于处理表格数据（DataFrame）。
import numpy as np # 导入numpy库，通常简写为np，提供强大的多维数组对象和数学函数，是科学计算的基础。
import matplotlib.pyplot as plt # 导入matplotlib的pyplot模块，简写为plt，是Python最常用的绘图库。
import matplotlib # 导入matplotlib主库，用于进行一些全局配置。
import os # 导入os库，用于与操作系统交互，如此处用来创建文件夹（os.makedirs）。
import warnings # 导入warnings库，用于控制程序运行时产生的警告信息的显示方式。
import re # 导入re库，用于正则表达式操作，如此处用于清理文件名中的非法字符。
from statsmodels.nonparametric.smoothers_lowess import lowess # 从statsmodels库导入lowess平滑函数，用于在散点图上绘制平滑的趋势线。
import shap # 导入shap库，用于模型解释，计算SHAP值来量化每个特征对模型预测的贡献。
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score # 从scikit-learn库导入数据划分、网格搜索和交叉验证工具。
import xgboost as xgb # 导入xgboost库，一个非常流行且高效的梯度提升框架。
import lightgbm as lgb # 导入lightgbm库，另一个高效、快速的梯度提升框架。
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error # 从scikit-learn库导入用于评估回归模型性能的指标函数。
import matplotlib.colors as mcolors # 导入matplotlib的颜色模块，用于处理颜色。
import seaborn as sns # 导入seaborn库，它基于matplotlib，提供了更美观、更高级的统计图形。
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor, AdaBoostRegressor, ExtraTreesRegressor # 导入多种集成学习回归模型。
from sklearn.svm import SVR # 导入支持向量回归模型。
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet # 导入线性回归模型及其变体。
from sklearn.tree import DecisionTreeRegressor # 导入决策树回归模型。
from sklearn.neighbors import KNeighborsRegressor # 导入K近邻回归模型。
import optuna # 导入optuna库，一个用于超参数自动优化的框架。
from optuna.samplers import TPESampler # 从optuna导入TPE采样器，一种先进的优化算法。
import time # 导入时间库，可以用来计算代码执行时间。
from sklearn.preprocessing import StandardScaler # 从scikit-learn导入StandardScaler，用于特征标准化。
from sklearn.pipeline import Pipeline # 从scikit-learn导入Pipeline，用于将多个处理步骤串联成一个工作流。
阶段 2：环境配置与辅助函数定义
作用与目的： 这个阶段设置了代码的运行环境。主要包括：1）抑制不必要的警告信息，保持输出界面的整洁。2）配置matplotlib绘图库，确保图表中的中文字符和负号能够正常显示。3）定义一个辅助函数sanitize_filename，用于处理文件名，确保在保存文件时不会因为非法字符而出错。
python
# 忽略警告
warnings.filterwarnings("ignore") # 设置警告过滤器，忽略所有类型的警告信息，使输出更简洁。
#matplotlib.use('TkAgg') # 这行被注释掉了。它的作用是设置matplotlib的后端，'TkAgg'是一个常用的图形界面后端，在某些环境下需要手动指定。
plt.rcParams['font.serif'] = ['Times New Roman'] # 设置matplotlib绘图时使用的衬线字体为'Times New Roman'。
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置matplotlib绘-图时使用的无衬线字体为'SimHei'（黑体），这是为了正确显示中文标签。
plt.rcParams['axes.unicode_minus'] = False# 设置matplotlib在有中文的环境下能够正常显示坐标轴的负号。
warnings.filterwarnings("ignore", category=RuntimeWarning) # 专门忽略运行时（RuntimeWarning）类型的警告。
warnings.filterwarnings("ignore", category=FutureWarning) # 专门忽略未来（FutureWarning）类型的警告，这些通常是关于库未来版本变化的提示。
#显示中文
plt.rcParams['font.family'] = 'SimHei'# 再次确认matplotlib绘图时使用的字体族为'SimHei'，确保中文显示。
#显示负号
plt.rcParams['axes.unicode_minus'] = False# 再次确认matplotlib正常显示负号。
defsanitize_filename(name): # 定义一个名为 sanitize_filename 的函数，接收一个字符串name作为参数。
return re.sub(r'[\\/*?:"<>|]', '_', name) # 使用正则表达式(re.sub)查找name中所有非法文件名字符（\/*?:"<>|），并用下划线'_'替换它们。
阶段 3：数据加载与预处理
作用与目的： 这是数据分析的起点。该阶段负责：1）根据用户配置的路径加载数据文件。2）进行关键的验证，确保目标列存在于数据中，如果不存在则终止程序并给出提示。3）将数据分割为特征集（X）和目标变量（y）。4）移除已知的、对模型无用的列（如ID、日期等）。
python
print("--- 开始执行任务 ---") # 在控制台打印一条消息，表示程序开始运行。
print(f"正在从 '{excel_file_path}' 加载数据...") # 打印提示信息，告知用户正在加载哪个文件。
if sheet_name: # 检查用户是否在配置中指定了工作表名称（sheet_name是否不为None）。
df = pd.read_csv(excel_file_path, sheet_name=sheet_name) # 如果指定了，则使用pandas的read_csv函数读取指定工作表的数据。
else: # 如果没有指定工作表名称。
df = pd.read_csv(excel_file_path) # 则默认读取CSV文件的全部内容。
print("数据加载成功！") # 如果数据成功加载，打印成功信息。
if target_column_name notin df.columns: # 检查用户指定的目标列名是否存在于数据框(df)的列名中。
print(f"错误：目标列 '{target_column_name}' 不存在于数据中！") # 如果不存在，打印错误信息。
print(f"可用的列为: {df.columns.tolist()}") # 并列出所有可用的列名，方便用户检查和修正配置。
exit() # 退出程序，因为缺少目标列无法继续分析。
y = df[target_column_name] # 从数据框(df)中取出目标列，赋值给变量y。
X = df.drop(columns=[target_column_name]) # 从数据框(df)中删除目标列，剩下的所有列作为特征集，赋值给X。
# 删除无用的特征列'user_id', 'date'
X.drop(columns=['user_id', 'date'], inplace=True, errors='ignore') # 从特征集X中删除'user_id'和'date'这两列。inplace=True表示直接在X上修改；errors='ignore'表示如果这些列不存在，也不会报错。
print("已删除无用特征列 'user_id' 和 'date'。") # 打印确认信息。
feature_names = X.columns.tolist() # 获取特征集X中所有列的名称，并转换为一个列表，方便后续使用。
阶段 4：探索性数据分析 (EDA) - 准备工作
作用与目的： 这个阶段为正式的EDA分析做准备。它会创建一个专门的文件夹来存放所有生成的图表和报告，这样可以保持工作目录的整洁。同时，它会设置全局的绘图风格和分辨率，确保所有输出的图表都具有统一、高质量的外观。
print("\n" + "=" * 60) # 打印一个换行符和一行等号，作为视觉分隔符。
print("🔮 开始探索性数据分析 (EDA) 🔮🔮") # 打印EDA部分的开始标志。
print("=" * 60) # 再次打印一行等号。
# 创建EDA输出文件夹
eda_output_folder = "EDA_Analysis_Results"# 定义一个字符串变量，作为存放EDA结果的文件夹名称。
os.makedirs(eda_output_folder, exist_ok=True) # 使用os库创建这个文件夹。exist_ok=True表示如果文件夹已经存在，则不执行任何操作，也不会报错。
print(f"📁 EDA分析结果将保存到: '{eda_output_folder}' 文件夹") # 打印提示信息，告知用户结果将保存在哪里。
# 设置绘图样式
plt.style.use('default') # 设置matplotlib使用默认的绘图风格。
sns.set_palette("husl") # 设置seaborn的调色板为"husl"，这是一套色彩丰富且均衡的颜色方案。
plt.rcParams['figure.dpi'] = 100# 设置在屏幕上显示图形的分辨率为100 DPI (每英寸点数)。
plt.rcParams['savefig.dpi'] = 300# 设置保存图形时的分辨率为300 DPI，以获得更高质量的图片文件。
阶段 5：EDA - 数据概览与统计
作用与目的： 这是EDA的第一步，目的是对数据有一个宏观的了解。它会分析：
- 基本信息
数据集有多大（行和列的数量）。
- 数据类型
数据集中包含哪些类型的数据（如整数、浮点数、对象等）。
- 缺失值
数据中是否存在缺失值，哪些列缺失最严重。这是数据清洗的关键步骤。
- 统计摘要
对于数值型特征，计算它们的核心统计量（如均值、标准差、最大/最小值等）。
所有这些分析结果都会被打印出来，并保存到CSV文件中，以备后续查阅。
python
# ===== 1.1.1 数据概览与统计 =====
print(f"\n⚡ 1.1.1 数据概览与统计信息分析") # 打印该子阶段的标题。
print("-" * 40) # 打印短分隔线。
# 基本信息
print(f"📊 数据集基本信息:") # 打印子标题。
print(f" - 数据集形状: {df.shape}") # 打印数据集的形状，即 (行数, 列数)。
print(f" - 总样本数: {len(df):,}") # 打印总样本数（行数），:, 格式化数字，添加千位分隔符。
print(f" - 总特征数: {len(df.columns)}") # 打印总特征数（列数）。
print(f" - 目标变量: {target_column_name}") # 打印目标变量的名称。
# 数据类型分析
print(f"\n📈 数据类型分布:") # 打印子标题。
dtype_counts = df.dtypes.value_counts() # 统计数据框中每种数据类型的列的数量。
for dtype, count in dtype_counts.items(): # 遍历统计结果。
print(f" - {dtype}: {count} 个特征") # 打印每种数据类型及其对应的列数。
# 缺失值分析
print(f"\n🔍 缺失值分析:") # 打印子标题。
missing_values = df.isnull().sum() # 计算每一列的缺失值（null）总数。
missing_percentage = (missing_values / len(df)) * 100# 计算每一列缺失值占总行数的百分比。
missing_summary = pd.DataFrame({ # 将缺失值统计结果整合成一个pandas DataFrame。
'Missing_Count': missing_values, # 第一列为缺失值数量。
'Missing_Percentage': missing_percentage # 第二列为缺失值百分比。
}).sort_values('Missing_Count', ascending=False) # 按缺失值数量从高到低排序。
print(f" - 总缺失值数: {missing_values.sum():,}") # 打印整个数据集的总缺失值数。
print(f" - 包含缺失值的特征数: {(missing_values > 0).sum()}") # 打印存在缺失值的列的总数。
if (missing_values > 0).sum() > 0: # 检查是否存在任何缺失值。
print(f"\n 前10个缺失值最多的特征:") # 如果有，则打印标题。
for idx, (feature, row) inenumerate(missing_summary.head(10).iterrows()): # 遍历缺失值最多的前10列。
if row['Missing_Count'] > 0: # 确保这一列确实有缺失值。
print(f" {idx + 1:2d}. {feature:<20}: {row['Missing_Count']:>6,} ({row['Missing_Percentage']:5.2f}%)") # 格式化打印排名、列名、缺失数和百分比。
else: # 如果没有任何缺失值。
print(" ✅ 数据集中无缺失值") # 打印一条确认信息。
# 保存缺失值分析结果
missing_summary.to_csv(os.path.join(eda_output_folder, 'missing_values_analysis.csv'), # 将缺失值分析的DataFrame保存为CSV文件。
index=True, encoding='utf-8-sig') # index=True保留索引列，encoding='utf-8-sig'确保中文在Excel中打开时不会乱码。
# 基本统计量
print(f"\n📊 数值特征基本统计量:") # 打印子标题。
numeric_columns = df.select_dtypes(include=[np.number]).columns # 筛选出所有数值类型的列。
desc_stats = df[numeric_columns].describe() # 对所有数值列计算描述性统计量（均值、标准差、分位数等）。
print(desc_stats.round(4)) # 将统计结果四舍五入到4位小数后打印。
# 保存统计量
desc_stats.to_csv(os.path.join(eda_output_folder, 'descriptive_statistics.csv'), # 将描述性统计结果保存为CSV文件。
index=True, encoding='utf-8-sig') # 同样，保留索引并使用正确的编码。
阶段 6：EDA - 目标变量分布分析
作用与目的： 深入分析要预测的目标变量 y。理解目标变量的分布至关重要，因为它可能影响模型的选择和评估。此阶段通过：
- 统计信息
计算均值、标准差、偏度（skewness）和峰度（kurtosis）等，从数值上描述其分布形态。
- 可视化
绘制四种图表（直方图、核密度图、Q-Q图、箱线图）从不同角度直观展示
y的分布，并判断其是否接近正态分布、是否存在异常值。
python
# ===== 1.1.2 目标变量分布分析 =====
print(f"\n🌟 1.1.2 目标变量 '{target_column_name}' 分布分析") # 打印该子阶段的标题。
print("-" * 40) # 打印分隔线。
# 目标变量基本统计
target_stats = y.describe() # 计算目标变量y的描述性统计量。
print(f"📈 目标变量统计信息:") # 打印子标题。
print(f" - 均值: {target_stats['mean']:.4f}") # 打印均值，保留4位小数。
print(f" - 标准差: {target_stats['std']:.4f}") # 打印标准差，保留4位小数。
print(f" - 最小值: {target_stats['min']:.4f}") # 打印最小值，保留4位小数。
print(f" - 最大值: {target_stats['max']:.4f}") # 打印最大值，保留4位小数。
print(f" - 中位数: {target_stats['50%']:.4f}") # 打印中位数（50%分位数），保留4位小数。
print(f" - 偏度: {y.skew():.4f}") # 打印偏度，衡量分布的不对称性。
print(f" - 峰度: {y.kurtosis():.4f}") # 打印峰度，衡量分布的尖锐程度。
# 绘制目标变量分布图
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12)) # 创建一个2x2的子图网格，整个图形大小为15x12英寸。
# 直方图
ax1.hist(y, bins=50, alpha=0.7, color='skyblue', edgecolor='black', density=True) # 在第一个子图(ax1)上绘制y的直方图，分为50个区间，设置透明度和颜色，density=True表示纵轴为概率密度。
ax1.axvline(y.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {y.mean():.3f}') # 在图上画一条红色的虚线表示均值。
ax1.axvline(y.median(), color='orange', linestyle='--', linewidth=2, label=f'Median: {y.median():.3f}') # 在图上画一条橙色的虚线表示中位数。
ax1.set_xlabel(f'{target_column_name}') # 设置ax1的x轴标签。
ax1.set_ylabel('Density') # 设置ax1的y轴标签。
ax1.set_title(f'Distribution of {target_column_name}', fontweight='bold') # 设置ax1的标题，并加粗。
ax1.legend() # 显示图例（均值和中位数的标签）。
ax1.grid(True, alpha=0.3) # 显示网格线，并设置透明度。
# 核密度估计图
y.plot.kde(ax=ax2, color='purple', linewidth=2) # 在第二个子图(ax2)上绘制y的核密度估计(KDE)图，这是一种平滑的直方图。
ax2.axvline(y.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {y.mean():.3f}') # 同样，在KDE图上画出均值线。
ax2.axvline(y.median(), color='orange', linestyle='--', linewidth=2, label=f'Median: {y.median():.3f}') # 同样，画出中位数线。
ax2.set_xlabel(f'{target_column_name}') # 设置ax2的x轴标签。
ax2.set_ylabel('Density') # 设置ax2的y轴标签。
ax2.set_title(f'Kernel Density Estimation of {target_column_name}', fontweight='bold') # 设置ax2的标题。
ax2.legend() # 显示图例。
ax2.grid(True, alpha=0.3) # 显示网格线。
# Q-Q图检验正态性
from scipy import stats # 从scipy库导入stats模块，用于统计检验。
stats.probplot(y, dist="norm", plot=ax3) # 在第三个子图(ax3)上绘制Q-Q图，用于直观地比较y的分布与正态分布的差异。
ax3.set_title(f'Q-Q Plot: {target_column_name} vs Normal Distribution', fontweight='bold') # 设置ax3的标题。
ax3.grid(True, alpha=0.3) # 显示网格线。
# 箱线图
box_data = ax4.boxplot(y, patch_artist=True, labels=[target_column_name]) # 在第四个子图(ax4)上绘制y的箱线图，用于展示数据的分位数和异常值。
box_data['boxes'][0].set_facecolor('lightgreen') # 设置箱线图的箱体颜色为浅绿色。
ax4.set_ylabel(f'{target_column_name}') # 设置ax4的y轴标签。
ax4.set_title(f'Box Plot of {target_column_name}', fontweight='bold') # 设置ax4的标题。
ax4.grid(True, alpha=0.3) # 显示网格线。
# 添加统计信息到箱线图
stats_text = f'Mean: {y.mean():.3f}\nStd: {y.std():.3f}\nSkew: {y.skew():.3f}\nKurt: {y.kurtosis():.3f}'# 创建一个包含均值、标准差、偏度和峰度信息的字符串。
ax4.text(0.02, 0.98, stats_text, transform=ax4.transAxes, fontsize=10, # 在ax4的左上角添加这个文本框。
verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8)) # 设置文本框的样式。
plt.suptitle(f'Target Variable Analysis: {target_column_name}', fontsize=16, fontweight='bold') # 为整个2x2图形设置一个总标题。
plt.tight_layout() # 自动调整子图的布局，防止标签重叠。
plt.savefig(os.path.join(eda_output_folder, 'target_variable_analysis.png'), dpi=300, bbox_inches='tight') # 将生成的图形保存到之前创建的文件夹中。
plt.show() # 在屏幕上显示图形。
阶段 7：EDA - 正态性检验
作用与目的： 这是对目标变量分析的补充。阶段6的可视化提供了直观判断，而此阶段则使用严格的统计检验方法（Shapiro-Wilk, D'Agostino, Jarque-Bera）来从数学上判断目标变量是否服从正态分布。检验结果（特别是p-value）可以为是否需要对目标变量进行变换（如对数变换）提供依据。
python
# 正态性检验
print(f"\n🔬 正态性检验结果:") # 打印子标题。
try: # 使用try-except块来捕获可能发生的错误（例如，库未安装）。
from scipy.stats import shapiro, normaltest, jarque_bera # 导入三种不同的正态性检验函数。
# Shapiro-Wilk检验（适合小样本）
iflen(y) <= 5000: # Shapiro-Wilk检验通常在样本量小于5000时效果较好。
shapiro_stat, shapiro_p = shapiro(y.sample(min(5000, len(y)))) # 对y（或其子样本）进行Shapiro-Wilk检验。
print(f" - Shapiro-Wilk检验: 统计量={shapiro_stat:.4f}, p-value={shapiro_p:.4e}") # 打印检验的统计量和p值。
# D'Agostino正态性检验
dagostino_stat, dagostino_p = normaltest(y) # 对y进行D'Agostino's K-squared检验。
print(f" - D'Agostino正态性检验: 统计量={dagostino_stat:.4f}, p-value={dagostino_p:.4e}") # 打印检验结果。
# Jarque-Bera检验
jb_stat, jb_p = jarque_bera(y) # 对y进行Jarque-Bera检验。
print(f" - Jarque-Bera检验: 统计量={jb_stat:.4f}, p-value={jb_p:.4e}") # 打印检验结果。
if dagostino_p < 0.05: # 通常以p值是否小于0.05作为判断标准。
print(f" ⚠️ 目标变量可能不符合正态分布 (p < 0.05)") # 如果p值小于0.05，则拒绝原假设，认为数据不服从正态分布。
print(f" 💡 建议考虑对数转换、Box-Cox变换等方法") # 给出相应的建议。
else: # 如果p值大于或等于0.05。
print(f" ✅ 目标变量基本符合正态分布 (p >= 0.05)") # 无法拒绝原假设，可以认为数据服从正态分布。
except Exception as e: # 如果try块中的代码出现任何异常。
print(f" ❌ 正态性检验失败: {e}") # 打印错误信息。
阶段 8：EDA - 特征与目标变量关系分析
作用与目的： 这个阶段的核心是找出哪些特征与我们的目标变量 y 关系最密切。这有助于我们理解业务、进行特征选择和建立更有效的模型。主要方法是：
- 计算相关性
计算所有数值特征与目标变量之间的皮尔逊相关系数。
- 排序和展示
列出与目标变量相关性最强和最弱的特征。
- 热力图
通过热力图可视化所有特征之间的相关性，帮助发现多重共线性问题。
- 散点图
对相关性最强的几个特征，绘制它们与目标变量的散点图，并拟合一条回归线，直观地观察它们之间的关系形态。
python
# ===== 1.1.3 特征与目标变量关系分析 =====
print(f"\n✨ 1.1.3 特征与目标变量关系分析") # 打印该子阶段的标题。
print("-" * 40) # 打印分隔线。
# 计算相关性矩阵
print(f"📊 正在计算特征相关性...") # 打印提示信息。
numeric_features = X.select_dtypes(include=[np.number]).columns # 再次获取所有数值特征的列名。
correlation_data = pd.concat([X[numeric_features], y], axis=1) # 将数值特征X和目标变量y合并成一个新的DataFrame。
correlation_matrix = correlation_data.corr() # 计算这个新DataFrame中所有列之间的相关性矩阵。
# 与目标变量的相关性
target_correlations = correlation_matrix[target_column_name].drop(target_column_name).sort_values(key=abs, ascending=False) # 从相关性矩阵中提取目标变量那一列，去掉它自身的相关性（值为1），并按相关系数的绝对值进行降序排序。
print(f"\n🔗 与目标变量相关性最强的前10个特征:") # 打印标题。
for i, (feature, corr) inenumerate(target_correlations.head(10).items()): # 遍历排序后最前面的10个特征。
print(f" {i + 1:2d}. {feature:<20}: {corr:>8.4f}") # 格式化打印排名、特征名和相关系数。
print(f"\n🔗 与目标变量相关性最弱的前5个特征:") # 打印标题。
for i, (feature, corr) inenumerate(target_correlations.tail(5).items()): # 遍历排序后最后面的5个特征。
print(f" {i + 1:2d}. {feature:<20}: {corr:>8.4f}") # 格式化打印排名、特征名和相关系数。
# 绘制相关性热力图
plt.figure(figsize=(14, 12)) # 创建一个新的图形，设置尺寸为14x12英寸。
mask = np.triu(np.ones_like(correlation_matrix.corr(), dtype=bool)) # 创建一个上三角矩阵的掩码(mask)，用于在热力图中只显示下三角部分，避免信息冗余。
sns.heatmap(correlation_matrix.corr(), mask=mask, annot=True, cmap='coolwarm', center=0, # 使用seaborn绘制热力图。mask应用掩码，annot=True在格子上显示数值，cmap设置颜色方案。
square=True, fmt='.2f', cbar_kws={"shrink": .8}, # square=True使格子为正方形，fmt='.2f'设置数值格式，cbar_kws调整颜色条。
xticklabels=True, yticklabels=True) # 确保x和y轴的标签都显示。
plt.title('Feature Correlation Heatmap', fontsize=16, fontweight='bold', pad=20) # 设置图表标题。
plt.xticks(rotation=45, ha='right') # 将x轴的标签旋转45度，防止重叠。
plt.yticks(rotation=0) # 保持y轴标签水平。
plt.tight_layout() # 自动调整布局。
plt.savefig(os.path.join(eda_output_folder, 'correlation_heatmap.png'), dpi=300, bbox_inches='tight') # 保存热力图。
plt.show() # 显示热力图。
# 绘制与目标变量相关性最强的特征散点图
top_features = target_correlations.head(6).index.tolist() # 获取与目标变量相关性最强的前6个特征的名称。
fig, axes = plt.subplots(2, 3, figsize=(18, 12)) # 创建一个2x3的子图网格。
axes = axes.flatten() # 将2x3的子图数组展平为一维数组，方便遍历。
for i, feature inenumerate(top_features): # 遍历这6个特征。
ax = axes[i] # 获取当前要绘制的子图。
# 散点图
ax.scatter(X[feature], y, alpha=0.6, s=30, color='steelblue', edgecolors='white', linewidth=0.5) # 在子图上绘制当前特征与目标变量y的散点图。
# 添加回归线
z = np.polyfit(X[feature], y, 1) # 使用numpy计算一元线性回归的系数。
p = np.poly1d(z) # 根据系数创建一个多项式函数。
ax.plot(X[feature], p(X[feature]), "r--", alpha=0.8, linewidth=2) # 将回归线绘制在散点图上。
# 计算R²
y_pred_simple = p(X[feature]) # 使用拟合的回归线进行预测。
r2_simple = r2_score(y, y_pred_simple) # 计算这个简单线性回归的R²分数。
ax.set_xlabel(feature, fontsize=11) # 设置子图的x轴标签。
ax.set_ylabel(target_column_name, fontsize=11) # 设置子图的y轴标签。
ax.set_title(f'{feature} vs {target_column_name}\nCorr: {target_correlations[feature]:.3f}, R²: {r2_simple:.3f}', # 设置子图的标题，包含相关系数和R²值。
fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3) # 显示网格线。
plt.suptitle('Top 6 Features vs Target Variable', fontsize=16, fontweight='bold') # 为整个图形设置总标题。
plt.tight_layout() # 自动调整布局。
plt.savefig(os.path.join(eda_output_folder, 'top_features_scatter_plots.png'), dpi=300, bbox_inches='tight') # 保存散点图。
plt.show() # 显示图形。
# 保存相关性分析结果
target_correlations.to_csv(os.path.join(eda_output_folder, 'target_correlations.csv'), # 将目标变量与各特征的相关性排序结果保存到CSV文件。
header=['Correlation'], encoding='utf-8-sig') # 设置表头和编码。
阶段 9：EDA - 特征分布分析
作用与目的： 这个阶段关注于每个数值特征自身的分布情况。通过为每个特征绘制直方图，我们可以了解：
- 集中趋势
数据的中心位置在哪里（均值、中位数）。
- 离散程度
数据的分布范围有多广。
- 分布形态
数据是单峰、双峰还是其他形态？是否对称（通过偏度skewness判断）？ 这有助于发现数据中可能存在的问题，例如数据倾斜严重，可能需要进行变换。
python
# ===== 1.1.4 特征分布分析 =====
print(f"\n🔮 1.1.4 特征分布分析") # 打印该子阶段的标题。
print("-" * 40) # 打印分隔线。
# 数值特征分布分析
print(f"📊 正在分析 {len(numeric_features)} 个数值特征的分布...") # 打印提示信息。
# 绘制所有数值特征的直方图
n_features = len(numeric_features) # 获取数值特征的数量。
n_cols = 4# 设定每行显示4个子图。
n_rows = (n_features + n_cols - 1) // n_cols # 计算需要的行数，这是一种向上取整的技巧。
if n_features > 0: # 检查是否存在数值特征。
fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, 5 * n_rows)) # 根据计算的行列数创建子图网格。
axes = axes.flatten() if n_rows > 1else [axes] if n_rows == 1else axes # 将子图数组展平，以方便使用索引进行遍历。
for i, feature inenumerate(numeric_features): # 遍历所有数值特征。
if i < len(axes): # 确保索引在子图数组范围内。
ax = axes[i] # 获取当前子图。
# 直方图
X[feature].hist(bins=30, alpha=0.7, ax=ax, color='lightblue', edgecolor='black') # 在当前子图上绘制该特征的直方图。
ax.axvline(X[feature].mean(), color='red', linestyle='--', linewidth=2, # 在图上画出均值线。
label=f'Mean: {X[feature].mean():.2f}')
ax.axvline(X[feature].median(), color='orange', linestyle='--', linewidth=2, # 在图上画出中位数线。
label=f'Median: {X[feature].median():.2f}')
ax.set_title(f'{feature}\nSkew: {X[feature].skew():.2f}', fontsize=10, fontweight='bold') # 设置子图标题，并显示该特征的偏度。
ax.set_xlabel(feature, fontsize=9) # 设置x轴标签。
ax.set_ylabel('Frequency', fontsize=9) # 设置y轴标签。
ax.legend(fontsize=8) # 显示图例。
ax.grid(True, alpha=0.3) # 显示网格线。
# 隐藏多余的子图
for i inrange(len(numeric_features), len(axes)): # 遍历那些没有被用来绘图的多余子图。
axes[i].set_visible(False) # 将它们隐藏。
plt.suptitle('Distribution of Numeric Features', fontsize=16, fontweight='bold') # 设置总标题。
plt.tight_layout() # 自动调整布局。
plt.savefig(os.path.join(eda_output_folder, 'feature_distributions.png'), dpi=300, bbox_inches='tight') # 保存图形。
plt.show() # 显示图形。
阶段 10：EDA - 离群值检测
作用与目的： 离群值（Outliers）是数据中与其他观测值显著不同的数据点，它们可能会严重影响某些模型的性能。此阶段使用箱线图和IQR（四分位距）方法来系统地检测每个数值特征中的离群值。
- 可视化
箱线图能直观地展示离群值。
- 量化分析
计算每个特征的离群值数量和比例，并按比例排序，帮助我们快速定位离群值问题最严重的特征。
python
# 绘制箱线图检测离群值
iflen(numeric_features) > 0: # 检查是否存在数值特征。
fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, 5 * n_rows)) # 再次创建子图网格。
axes = axes.flatten() if n_rows > 1else [axes] if n_rows == 1else axes # 展平子图数组。
outlier_summary = {} # 创建一个空字典，用于存储每个特征的离群值分析结果。
for i, feature inenumerate(numeric_features): # 遍历所有数值特征。
if i < len(axes): # 确保索引在子图数组范围内。
ax = axes[i] # 获取当前子图。
# 箱线图
box_data = ax.boxplot(X[feature], patch_artist=True, labels=[feature]) # 在子图上绘制该特征的箱线图。
box_data['boxes'][0].set_facecolor('lightgreen') # 设置箱体颜色。
# 计算离群值
Q1 = X[feature].quantile(0.25) # 计算第一四分位数（Q1）。
Q3 = X[feature].quantile(0.75) # 计算第三四分位数（Q3）。
IQR = Q3 - Q1 # 计算四分位距（IQR）。
lower_bound = Q1 - 1.5 * IQR # 定义离群值的下界。
upper_bound = Q3 + 1.5 * IQR # 定义离群值的上界。
outliers = X[feature][(X[feature] < lower_bound) | (X[feature] > upper_bound)] # 筛选出所有低于下界或高于上界的数据点，即离群值。
outlier_percentage = len(outliers) / len(X[feature]) * 100# 计算离群值占该特征总数据点的百分比。
outlier_summary[feature] = { # 将该特征的离群值信息存入字典。
'count': len(outliers), # 离群值数量。
'percentage': outlier_percentage, # 离群值百分比。
'lower_bound': lower_bound, # 下界。
'upper_bound': upper_bound # 上界。
}
ax.set_title(f'{feature}\nOutliers: {len(outliers)} ({outlier_percentage:.1f}%)', # 设置子图标题，显示离群值数量和百分比。
fontsize=10, fontweight='bold')
ax.set_ylabel(feature, fontsize=9) # 设置y轴标签。
ax.grid(True, alpha=0.3) # 显示网格线。
# 隐藏多余的子图
for i inrange(len(numeric_features), len(axes)): # 遍历多余的子图。
axes[i].set_visible(False) # 将它们隐藏。
plt.suptitle('Box Plots for Outlier Detection', fontsize=16, fontweight='bold') # 设置总标题。
plt.tight_layout() # 自动调整布局。
plt.savefig(os.path.join(eda_output_folder, 'feature_boxplots.png'), dpi=300, bbox_inches='tight') # 保存图形。
plt.show() # 显示图形。
# 输出离群值分析结果
print(f"\n🚨 离群值检测结果 (使用IQR方法):") # 打印子标题。
outlier_df = pd.DataFrame(outlier_summary).T # 将存储离群值信息的字典转换为DataFrame，并转置（.T）使特征名为行索引。
outlier_df = outlier_df.sort_values('percentage', ascending=False) # 按离群值百分比降序排列。
print(f" {'特征':<20}{'离群值数量':<10}{'离群值比例':<10}{'下界':<12}{'上界':<12}") # 打印格式化的表头。
print("-" * 70) # 打印分隔线。
for feature, row in outlier_df.head(10).iterrows(): # 遍历离群值比例最高的前10个特征。
print(
f" {feature:<20}{row['count']:<10.0f}{row['percentage']:<9.1f}% {row['lower_bound']:<12.2f}{row['upper_bound']:<12.2f}") # 格式化打印每一行的信息。
# 保存离群值分析结果
outlier_df.to_csv(os.path.join(eda_output_folder, 'outlier_analysis.csv'), # 将离群值分析的DataFrame保存为CSV文件。
index=True, encoding='utf-8-sig') # 保留索引并使用正确编码。
阶段 11：EDA - 特征重要性初步分析
作用与目的： 在建模之前，初步评估哪些特征可能更重要，有助于后续的特征选择。这里用了两种简单而有效的方法：
- 方差分析
方差大的特征通常包含更多的信息。如果一个特征的方差接近于0，说明它在所有样本上都差不多，对区分样本可能没什么用。
- 互信息
互信息可以捕捉特征与目标变量之间的任何关系（包括非线性关系），比只衡量线性关系的相关系数更强大。得分高的特征与目标变量的关联性更强。
python
# ===== 1.1.5 特征重要性初步分析 =====
print(f"\n⚡⚡ 1.1.5 特征重要性初步分析") # 打印该子阶段的标题。
print("-" * 40) # 打印分隔线。
# 使用方差分析特征重要性
feature_variance = X[numeric_features].var().sort_values(ascending=False) # 计算所有数值特征的方差，并按方差大小降序排列。
print(f"📊 特征方差排序 (前10个):") # 打印标题。
for i, (feature, variance) inenumerate(feature_variance.head(10).items()): # 遍历方差最大的前10个特征。
print(f" {i + 1:2d}. {feature:<20}: {variance:>12.4f}") # 格式化打印排名、特征名和方差值。
# 使用互信息评估特征重要性
try: # 使用try-except块处理可能的错误。
from sklearn.feature_selection import mutual_info_regression # 从scikit-learn导入用于回归问题的互信息计算函数。
print(f"\n🔍 正在计算互信息...") # 打印提示信息。
mi_scores = mutual_info_regression(X[numeric_features], y, random_state=0) # 计算所有数值特征与目标变量y之间的互信息分数。
mi_scores_df = pd.DataFrame({ # 将结果整理成DataFrame。
'feature': numeric_features, # 第一列为特征名。
'mutual_info_score': mi_scores # 第二列为对应的互信息分数。
}).sort_values('mutual_info_score', ascending=False) # 按互信息分数降序排列。
print(f"🔗 互信息得分排序 (前10个):") # 打印标题。
for i, row in mi_scores_df.head(10).iterrows(): # 遍历分数最高的前10个特征。
print(f" {i + 1:2d}. {row['feature']:<20}: {row['mutual_info_score']:>12.6f}") # 格式化打印排名、特征名和互信息分数。
# 保存互信息结果
mi_scores_df.to_csv(os.path.join(eda_output_folder, 'mutual_information_scores.csv'), # 将互信息分数保存为CSV文件。
index=False, encoding='utf-8-sig') # index=False表示不保存行索引。
except Exception as e: # 如果计算失败。
print(f" ❌ 互信息计算失败: {e}") # 打印错误信息。
阶段 12：EDA - 生成数据质量报告
作用与目的： 这是整个EDA流程的收尾工作。它将前面所有分析的关键结论汇总到一个文本文件中，形成一份简洁的数据质量报告。这份报告可以快速地让其他人（或未来的自己）了解这份数据的基本情况、主要特点和潜在问题，非常便于沟通和存档。
python
# ===== 1.1.6 生成数据质量报告 =====
print(f"\n🌟 1.1.6 生成数据质量报告") # 打印该子阶段的标题。
print("-" * 40) # 打印分隔线。
# 创建综合数据质量报告
quality_report_path = os.path.join(eda_output_folder, 'data_quality_report.txt') # 定义报告文件的完整路径。
withopen(quality_report_path, 'w', encoding='utf-8') as f: # 使用'with open'语法打开文件进行写入（'w'），可以确保文件使用后被自动关闭。encoding='utf-8'支持中文。
f.write("=" * 60 + "\n") # 写入一行等号作为报告头。
f.write("📊 数据质量与探索性分析报告\n") # 写入报告主标题。
f.write("=" * 60 + "\n\n") # 写入分隔线和空行。
f.write("1. 数据集基本信息:\n") # 写入第一部分标题。
f.write(f" - 数据集形状: {df.shape}\n") # 写入数据集形状。
f.write(f" - 总样本数: {len(df):,}\n") # 写入总样本数。
f.write(f" - 总特征数: {len(df.columns)}\n") # 写入总特征数。
f.write(f" - 数值特征数: {len(numeric_features)}\n") # 写入数值特征数。
f.write(f" - 目标变量: {target_column_name}\n\n") # 写入目标变量名。
f.write("2. 数据质量评估:\n") # 写入第二部分标题。
f.write(f" - 总缺失值: {missing_values.sum():,}\n") # 写入总缺失值数。
f.write(f" - 缺失值比例: {(missing_values.sum() / (df.size) * 100):.2f}%\n") # 写入总缺失值比例（df.size是总单元格数）。
f.write(f" - 包含缺失值的特征数: {(missing_values > 0).sum()}\n\n") # 写入包含缺失值的特征数。
第二部分：13种模型性能对比
在完成初步的数据探索性分析之后，接下来的核心目标是进行一次广泛的模型性能横向评测。代码中预定义了13种不同原理的回归模型，涵盖了线性模型、支持向量机、近邻算法、决策树以及多种先进的集成学习算法（如随机森林、梯度提升树等）。进行这种“海选”式评估的目的，并非为了立即确定一个最优模型，而是为了快速筛选出在此特定数据集上表现出众、具有进一步优化潜力的模型类别。这种方法可以有效避免在不适合该数据特性的模型上投入过多的调优精力。
阶段一：数据集划分与多模型性能对比
为确保评估的全面性和可靠性，每个模型都将从多个维度进行考量：包括在训练集和测试集上的性能指标（用以判断过拟合程度）、交叉验证得分（用以评估模型的稳定性和鲁棒性）以及训练所需时间（作为衡量模型效率的实用性指标）。整个训练、预测、评估的过程被设计为自动化循环，最终将所有模型的性能数据整合到一个结构化的表格中，为后续的可视化分析和模型选择奠定坚实基础。
print("\n正在划分训练集和验证集...") # 在控制台输出提示信息，表明数据集划分步骤开始。
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3,
random_state=0) # 调用train_test_split函数划分数据。X为特征集，y为目标变量。test_size=0.3指定30%的数据作为测试集，其余70%为训练集。random_state=0设定随机种子，确保每次代码运行时划分结果保持一致，便于实验复现。
print(f"数据集划分完成。训练集: {X_train.shape[0]}样本, 验证集: {X_test.shape[0]}样本。") # 输出划分完成后训练集与测试集的样本数量。
# 定义13种机器学习模型
models = {
'XGBoost': xgb.XGBRegressor(random_state=0, eval_metric='rmse'), # 定义XGBoost回归模型，设定随机种子和评估指标。
'LightGBM': lgb.LGBMRegressor(random_state=0, verbose=-1), # 定义LightGBM回归模型，设定随机种子并抑制训练过程中的冗余输出。
'Random Forest': RandomForestRegressor(random_state=0), # 定义随机森林回归模型，设定随机种子。
'Gradient Boosting': GradientBoostingRegressor(random_state=0), # 定义梯度提升回归模型，设定随机种子。
'Extra Trees': ExtraTreesRegressor(random_state=0), # 定义极端随机树回归模型，设定随机种子。
'AdaBoost': AdaBoostRegressor(random_state=0), # 定义AdaBoost回归模型，设定随机种子。
'Support Vector Machine': SVR(), # 定义支持向量回归（SVR）模型。
'Linear Regression': LinearRegression(), # 定义线性回归模型。
'Ridge Regression': Ridge(random_state=0), # 定义岭回归模型，它是线性回归的L2正则化版本。
'Lasso Regression': Lasso(random_state=0), # 定义Lasso回归模型，它是线性回归的L1正则化版本。
'ElasticNet': ElasticNet(random_state=0), # 定义弹性网络回归模型，它结合了L1和L2正则化。
'Decision Tree': DecisionTreeRegressor(random_state=0), # 定义决策树回归模型，设定随机种子。
'K-Nearest Neighbors': KNeighborsRegressor() # 定义K近邻回归模型。
}
print(f"\n开始训练和评估 {len(models)} 种机器学习模型...") # 输出提示信息，告知将要开始对所有定义的模型进行训练和评估。
# 存储模型性能结果
model_results = [] # 初始化一个空列表，用于收集每个模型运行后的性能指标字典。
trained_models = {} # 初始化一个空字典，用于存储每个训练完成的模型对象，以备后续调用。
for name, model in models.items(): # 启动一个循环，遍历上方定义的models字典中的每一个模型。
print(f"\n正在训练 {name}...") # 输出当前正在训练的模型名称。
start_time = time.time() # 记录当前时间，作为模型训练的起始时间。
# 对于需要缩放的模型，创建Pipeline
if name in ['Support Vector Machine', 'Linear Regression', 'Ridge Regression',
'Lasso Regression', 'ElasticNet', 'K-Nearest Neighbors']: # 判断当前模型是否对特征尺度敏感（如SVM、线性模型、KNN等）。
# 创建一个包含缩放和模型的管道
pipeline = Pipeline([
('scaler', StandardScaler()), # 定义管道的第一步：使用StandardScaler对数据进行特征标准化处理。
('model', model) # 定义管道的第二步：将标准化后的数据传递给当前模型。
])
trainable_model = pipeline # 对于这些模型，将整个管道作为训练对象。
else:
# 树模型不需要缩放
trainable_model = model # 对于像决策树和其集成模型这类对特征尺度不敏感的模型，直接使用原始模型对象。
try: # 使用try-except代码块来捕获并处理模型训练中可能出现的异常，避免因单个模型失败而中断整个流程。
# 训练模型
trainable_model.fit(X_train, y_train) # 调用fit方法，使用训练集数据对模型（或管道）进行训练。
# 预测
y_pred_train = trainable_model.predict(X_train) # 使用训练好的模型对训练集自身进行预测，用于评估拟合程度。
y_pred_test = trainable_model.predict(X_test) # 使用训练好的模型对测试集进行预测，用于评估泛化能力。
# 计算评估指标
train_r2 = r2_score(y_train, y_pred_train) # 计算训练集上的决定系数R²。
test_r2 = r2_score(y_test, y_pred_test) # 计算测试集上的决定系数R²。
train_rmse = np.sqrt(mean_squared_error(y_train, y_pred_train)) # 计算训练集上的均方根误差RMSE。
test_rmse = np.sqrt(mean_squared_error(y_test, y_pred_test)) # 计算测试集上的均方根误差RMSE。
train_mae = mean_absolute_error(y_train, y_pred_train) # 计算训练集上的平均绝对误差MAE。
test_mae = mean_absolute_error(y_test, y_pred_test) # 计算测试集上的平均绝对误差MAE。
# 计算交叉验证分数
cv_scores = cross_val_score(trainable_model, X_train, y_train, cv=5, scoring='r2') # 在训练集上执行5折交叉验证，评估指标为R²。
cv_mean = cv_scores.mean() # 计算5次交叉验证R²分数的平均值。
cv_std = cv_scores.std() # 计算5次交叉验证R²分数的标准差。
training_time = time.time() - start_time # 计算当前时间与起始时间的差值，得到模型训练耗时。
# 存储结果
model_results.append({ # 将当前模型的所有性能指标打包成一个字典，并添加到结果列表中。
'Model': name,
'Train_R2': train_r2,
'Test_R2': test_r2,
'Train_RMSE': train_rmse,
'Test_RMSE': test_rmse,
'Train_MAE': train_mae,
'Test_MAE': test_mae,
'CV_R2_Mean': cv_mean,
'CV_R2_Std': cv_std,
'Training_Time': training_time
})
# 保存训练好的模型
trained_models[name] = model # 将训练好的原始模型对象存入字典，键为模型名称。
print(f" {name} 训练完成:") # 输出训练完成的确认信息。
print(f" 训练集 R2: {train_r2:.4f}, 测试集 R2: {test_r2:.4f}") # 输出模型在训练集和测试集上的R²分数。
print(f" 测试集 RMSE: {test_rmse:.4f}, 测试集 MAE: {test_mae:.4f}") # 输出模型在测试集上的RMSE和MAE。
print(f" 交叉验证 R2: {cv_mean:.4f} (±{cv_std:.4f})") # 输出交叉验证的R²均值和标准差。
print(f" 训练时间: {training_time:.2f}秒") # 输出模型的训练耗时。
except Exception as e: # 如果在try代码块中发生任何异常。
print(f" {name} 训练失败: {str(e)}") # 输出失败的模型名称及具体的错误信息。
continue# 跳过当前循环的剩余部分，继续训练下一个模型。
# 创建结果DataFrame
results_df = pd.DataFrame(model_results) # 将收集了所有模型结果的列表转换为一个pandas DataFrame。
results_df = results_df.sort_values('Test_R2', ascending=False) # 根据测试集R²（Test_R2）这一列进行降序排序，以评估模型性能。
print("\n所有模型训练完成！模型性能排名（按测试集R2排序）:") # 输出所有模型训练结束的总结性信息。
print("=" * 80) # 输出一行分隔符以增强可读性。
for idx, row in results_df.iterrows(): # 遍历排序后的结果DataFrame的每一行。
print(
f"{row.name + 1:2d}. {row['Model']:20s} | 测试R2: {row['Test_R2']:.4f} | 测试RMSE: {row['Test_RMSE']:.4f} | CV_R2: {row['CV_R2_Mean']:.4f}") # 格式化输出一个简洁的性能排名列表。
阶段二：模型性能可视化对比分析
此阶段通过绘制一组条形图，从四个关键维度——模型准确性（R²分数）、预测误差（RMSE）、模型稳定性（交叉验证得分）和计算效率（训练时间）——对所有模型进行了并列比较。这种多维度的可视化不仅能清晰地展示出哪个模型的预测最准确，还能揭示模型是否存在过拟合风险（通过比较训练集与测试集的R²差距）、性能是否稳健（通过观察交叉验证得分的均值和误差棒）以及其在实际部署中的可行性（通过考量训练时间）。在此基础上，进一步为排名前五的优胜模型生成了一张综合性能雷达图。雷达图是一种高级的可视化技术，它能够将多个不同量纲的性能指标（经过归一化处理）呈现在同一张图中，形成每个模型的独特“性能指纹”。通过对比不同模型在雷达图上围成的区域面积和形状，可以极其直观地对模型的综合实力做出判断，为最终的模型选择提供强有力的决策支持。
# ===== 绘制模型性能对比图 =====
print("\n正在绘制模型性能对比图...") # 输出提示信息，表明开始进入可视化阶段。
# 创建输出目录
output_dir = "model_comparison_results"# 定义一个用于存放输出图片的文件目录名称。
os.makedirs(output_dir, exist_ok=True) # 创建该目录，如果目录已存在则不执行任何操作。
# 设置图表样式
plt.style.use('default') # 将matplotlib的绘图风格设置为默认。
sns.set_palette("husl") # 使用seaborn的"husl"调色板，以获得美观且区分度高的颜色。
# 图1: R2分数对比
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12)) # 创建一个2x2的子图网格，用于容纳四个对比图，并设定整个画布的大小。
# R2分数对比
models_sorted = results_df['Model'].tolist() # 从已排序的DataFrame中获取模型名称列表。
train_r2_sorted = results_df['Train_R2'].tolist() # 从已排序的DataFrame中获取训练集R²分数列表。
test_r2_sorted = results_df['Test_R2'].tolist() # 从已排序的DataFrame中获取测试集R²分数列表。
x_pos = np.arange(len(models_sorted)) # 为每个模型生成一个在x轴上的位置索引。
width = 0.35# 定义分组柱状图中每个柱子的宽度。
bars1 = ax1.bar(x_pos - width / 2, train_r2_sorted, width, label='Train R²', alpha=0.8) # 在第一个子图(ax1)上，绘制训练集R²的柱状图，位置向左微调。
bars2 = ax1.bar(x_pos + width / 2, test_r2_sorted, width, label='Test R²', alpha=0.8) # 在第一个子图(ax1)上，绘制测试集R²的柱状图，位置向右微调，形成对比。
ax1.set_xlabel('Models', fontsize=12) # 设置ax1的x轴标签。
ax1.set_ylabel('R² Score', fontsize=12) # 设置ax1的y轴标签。
ax1.set_title('Model R² Score Comparison', fontsize=14, fontweight='bold') # 设置ax1的标题。
ax1.set_xticks(x_pos) # 指定x轴刻度的位置。
ax1.set_xticklabels(models_sorted, rotation=45, ha='right') # 将x轴的刻度标签设置为模型名称，并旋转45度以防重叠。
ax1.legend() # 显示图例（'Train R²' 和 'Test R²'）。
ax1.grid(True, alpha=0.3) # 显示背景网格线。
# 添加数值标签
for bar in bars1: # 遍历代表训练集R²的每个柱子。
height = bar.get_height() # 获取柱子的高度。
ax1.text(bar.get_x() + bar.get_width() / 2., height + 0.01, # 在柱子顶端中心位置添加文本。
f'{height:.3f}', ha='center', va='bottom', fontsize=8) # 文本内容为R²分数，格式化为三位小数。
for bar in bars2: # 遍历代表测试集R²的每个柱子。
height = bar.get_height() # 获取柱子的高度。
ax1.text(bar.get_x() + bar.get_width() / 2., height + 0.01, # 同样在柱子顶端添加数值标签。
f'{height:.3f}', ha='center', va='bottom', fontsize=8)
# RMSE对比
test_rmse_sorted = results_df['Test_RMSE'].tolist() # 从已排序的DataFrame中获取测试集RMSE列表。
bars3 = ax2.bar(models_sorted, test_rmse_sorted, alpha=0.8, color='coral') # 在第二个子图(ax2)上绘制RMSE的柱状图。
ax2.set_xlabel('Models', fontsize=12) # 设置ax2的x轴标签。
ax2.set_ylabel('RMSE', fontsize=12) # 设置ax2的y轴标签。
ax2.set_title('Model RMSE Comparison (Test Set)', fontsize=14, fontweight='bold') # 设置ax2的标题。
ax2.tick_params(axis='x', rotation=45) # 旋转x轴的刻度标签。
ax2.grid(True, alpha=0.3) # 显示背景网格线。
for bar in bars3: # 遍历代表RMSE的每个柱子。
height = bar.get_height() # 获取柱子的高度。
ax2.text(bar.get_x() + bar.get_width() / 2., height + height * 0.01, # 在柱子顶端添加数值标签。
f'{height:.3f}', ha='center', va='bottom', fontsize=8)
# 交叉验证分数对比
cv_means = results_df['CV_R2_Mean'].tolist() # 获取交叉验证R²均值列表。
cv_stds = results_df['CV_R2_Std'].tolist() # 获取交叉验证R²标准差列表。
bars4 = ax3.bar(models_sorted, cv_means, yerr=cv_stds, alpha=0.8, color='lightgreen', capsize=5) # 在第三个子图(ax3)上绘制柱状图，并使用yerr参数添加表示标准差的误差棒。
ax3.set_xlabel('Models', fontsize=12) # 设置ax3的x轴标签。
ax3.set_ylabel('Cross-Validation R² Score', fontsize=12) # 设置ax3的y轴标签。
ax3.set_title('Model Cross-Validation Performance Comparison', fontsize=14, fontweight='bold') # 设置ax3的标题。
ax3.tick_params(axis='x', rotation=45) # 旋转x轴的刻度标签。
ax3.grid(True, alpha=0.3) # 显示背景网格线。
# 训练时间对比
training_times = results_df['Training_Time'].tolist() # 获取训练时间列表。
bars5 = ax4.bar(models_sorted, training_times, alpha=0.8, color='gold') # 在第四个子图(ax4)上绘制训练时间的柱状图。
ax4.set_xlabel('Models', fontsize=12) # 设置ax4的x轴标签。
ax4.set_ylabel('Training Time (seconds)', fontsize=12) # 设置ax4的y轴标签。
ax4.set_title('Model Training Time Comparison', fontsize=14, fontweight='bold') # 设置ax4的标题。
ax4.tick_params(axis='x', rotation=45) # 旋转x轴的刻度标签。
ax4.grid(True, alpha=0.3) # 显示背景网格线。
plt.tight_layout() # 自动调整所有子图的布局，以避免元素重叠。
plt.savefig(os.path.join(output_dir, 'model_performance_comparison.png'), dpi=300, bbox_inches='tight') # 将整个2x2的图形以300 DPI的高分辨率保存到指定目录。
plt.show() # 在屏幕上显示该图形。
# 综合性能雷达图
print("正在绘制综合性能雷达图...") # 输出提示信息。
top_5_models = results_df.head(5) # 从排序后的结果中选取表现最好的前5个模型。
# 准备雷达图数据
categories = ['Test_R2', 'CV_R2_Mean', 'Training_Speed'] # 定义雷达图的三个性能评估维度。
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection='polar')) # 创建一个使用极坐标系(projection='polar')的子图，这是绘制雷达图的基础。
# 归一化数据（训练时间需要反向归一化，时间越短越好）
normalized_data = {} # 初始化一个字典，用于存储归一化后的数据。
for model_name in top_5_models['Model']: # 遍历前5个模型的名称。
model_row = top_5_models[top_5_models['Model'] == model_name].iloc[0] # 提取该模型对应的性能数据行。
# 归一化到0-1之间
test_r2_norm = (model_row['Test_R2'] - results_df['Test_R2'].min()) / ( # 对Test_R2进行最小-最大归一化，使其值落在0-1之间。
results_df['Test_R2'].max() - results_df['Test_R2'].min())
cv_r2_norm = (model_row['CV_R2_Mean'] - results_df['CV_R2_Mean'].min()) / ( # 对交叉验证R²均值进行最小-最大归一化。
results_df['CV_R2_Mean'].max() - results_df['CV_R2_Mean'].min())
# 训练时间反向归一化（时间越短得分越高）
time_norm = 1 - (model_row['Training_Time'] - results_df['Training_Time'].min()) / ( # 对训练时间进行反向归一化，时间最短的模型得分接近1，最长的接近0。
results_df['Training_Time'].max() - results_df['Training_Time'].min())
normalized_data[model_name] = [test_r2_norm, cv_r2_norm, time_norm] # 将该模型归一化后的三个性能值存入字典。
# 绘制雷达图
angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist() # 计算每个维度在雷达图上对应的角度。
angles += angles[:1] # 复制第一个角度值并追加到列表末尾，以使雷达图的路径闭合。
colors = plt.cm.Set3(np.linspace(0, 1, len(top_5_models))) # 从Set3调色板中为5个模型生成一组不同的颜色。
for i, model_name inenumerate(top_5_models['Model']): # 遍历前5个模型。
values = normalized_data[model_name] # 获取该模型归一化后的性能值。
values += values[:1] # 同样地，复制第一个值并追加到末尾以闭合图形。
ax.plot(angles, values, 'o-', linewidth=2, label=model_name, color=colors[i]) # 绘制代表该模型性能的多边形轮廓线。
ax.fill(angles, values, alpha=0.25, color=colors[i]) # 填充该多边形区域，并设置透明度，以增强视觉效果。
ax.set_xticks(angles[:-1]) # 在雷达图上设置每个维度轴线的位置。
ax.set_xticklabels(['Test R²', 'Cross-Validation R²', 'Training Speed']) # 为每个维度轴线设置标签。
ax.set_ylim(0, 1) # 设置雷达图的数值范围为0到1。
ax.set_title('Top 5 Models Comprehensive Performance Radar Chart', size=16, fontweight='bold', pad=20) # 设置雷达图的标题。
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0)) # 显示图例，并将其放置在图表的右侧外部。
ax.grid(True) # 显示雷达图的网格线。
plt.savefig(os.path.join(output_dir, 'top5_models_radar_chart.png'), dpi=300, bbox_inches='tight') # 将雷达图以高分辨率保存为图片。
plt.show() # 在屏幕上显示雷达图。
第三部分：统计显著性检验
在完成初步的模型性能排序之后，一个关键的问题随之而来：排名第一的模型比排名第二的模型好，这种“好”在统计学上是显著的吗？或者说，这种性能上的微小差异是否可能仅仅是由于数据划分的随机性造成的？为了回答这个问题，我们需要引入统计显著性检验。
此阶段的目的是对排名前列的模型进行成对的、严谨的统计比较。它不再仅仅关注性能指标的均值，而是深入分析性能指标的分布。通过一系列的统计计算，我们可以量化地评估两个模型性能差异的置信度。
阶段 1: 初始化与设置
这个阶段主要是为接下来的统计分析做准备工作。它包括打印一个清晰的标题、创建用于存储结果的文件夹，以及导入所有必需的统计分析库。
作用解释: 此阶段的代码确保了分析环境的整洁和所需工具的可用性。通过创建一个专门的文件夹，所有生成的报告和图表都会被整齐地组织在一起，方便后续查找。导入 scipy 和 statsmodels 库则为后续进行 t检验、Wilcoxon检验和多重比较校正等高级统计操作提供了函数支持。
# ===== 统计显著性检验 (Statistical Significance Testing) =====
print("\n" + "=" * 70)
print("📊 开始统计显著性检验分析 📊")
print("=" * 70)
# 创建统计检验输出文件夹
stats_test_folder = os.path.join(output_dir, "statistical_significance_tests")
os.makedirs(stats_test_folder, exist_ok=True)
print(f"📁 统计检验结果将保存到: '{stats_test_folder}' 文件夹")
# 导入必要的统计库
from scipy import stats
from scipy.stats import ttest_rel, wilcoxon, normaltest, levene
import itertools
from statsmodels.stats.multitest import multipletests
阶段 2: 收集详细的交叉验证分数
在这个阶段，代码将对性能排名前五的模型进行更详细、更严谨的交叉验证。之前的模型评估可能使用了较少的折数或不同的随机种子，而这里我们使用统一的10折交叉验证来为每个模型收集10个性能得分。这些得分是后续所有统计检验的基础。
作用解释: 为了进行可靠的统计比较，我们需要每个模型在一系列相同的数据子集上的性能表现，而不是单一的平均分。通过对每个模型运行10折交叉验证，我们得到了10个R²分数。这个分数列表（或分布）比单一的平均分包含了更多的信息，如性能的稳定性和波动范围，这对于判断模型间的差异是否真实存在至关重要。
# ===== 收集详细的交叉验证分数 =====
print("\n🔬 2.1.1 收集详细交叉验证分数进行统计检验")
print("-" * 50)
# 选择Top 5模型进行详细分析
top_5_models_names = results_df.head(5)['Model'].tolist()
print(f"📈 将对以下Top 5模型进行统计显著性检验:")
for i, model_name inenumerate(top_5_models_names):
print(f" {i + 1}. {model_name}")
# 设置交叉验证参数
from sklearn.model_selection import StratifiedKFold, KFold
cv_folds = 10# 使用10折交叉验证获得更稳定的结果
random_state = 42
print(f"\n🔄 正在进行 {cv_folds} 折交叉验证收集详细分数...")
# 存储每个模型的详细CV分数
detailed_cv_scores = {}
cv_statistics = {}
for model_name in top_5_models_names:
print(f" -> 正在评估: {model_name}")
# 获取对应的模型
if model_name in trained_models:
model = trained_models[model_name]
else:
# 如果模型不在已训练模型中，重新创建
if model_name == 'XGBoost':
model = xgb.XGBRegressor(random_state=random_state, eval_metric='rmse')
elif model_name == 'LightGBM':
model = lgb.LGBMRegressor(random_state=random_state, verbose=-1)
elif model_name == 'Random Forest':
model = RandomForestRegressor(random_state=random_state)
elif model_name == 'Gradient Boosting':
model = GradientBoostingRegressor(random_state=random_state)
elif model_name == 'Extra Trees':
model = ExtraTreesRegressor(random_state=random_state)
elif model_name == 'AdaBoost':
model = AdaBoostRegressor(random_state=random_state)
elif model_name == 'Support Vector Machine':
model = SVR()
elif model_name == 'Linear Regression':
model = LinearRegression()
elif model_name == 'Ridge Regression':
model = Ridge(random_state=random_state)
elif model_name == 'Lasso Regression':
model = Lasso(random_state=random_state)
elif model_name == 'ElasticNet':
model = ElasticNet(random_state=random_state)
elif model_name == 'Decision Tree':
model = DecisionTreeRegressor(random_state=random_state)
elif model_name == 'K-Nearest Neighbors':
model = KNeighborsRegressor()
try:
# 执行详细的交叉验证
kf = KFold(n_splits=cv_folds, shuffle=True, random_state=random_state)
cv_scores = cross_val_score(model, X_train, y_train, cv=kf, scoring='r2', n_jobs=-1)
# 存储结果
detailed_cv_scores[model_name] = cv_scores
# 计算统计量
cv_statistics[model_name] = {
'mean': cv_scores.mean(),
'std': cv_scores.std(),
'min': cv_scores.min(),
'max': cv_scores.max(),
'median': np.median(cv_scores),
'q25': np.percentile(cv_scores, 25),
'q75': np.percentile(cv_scores, 75),
'cv_coefficient': cv_scores.std() / cv_scores.mean() if cv_scores.mean() != 0else np.inf
}
print(f" ✅ 完成 - 平均R²: {cv_scores.mean():.4f} (±{cv_scores.std():.4f})")
except Exception as e:
print(f" ❌ 失败: {str(e)}")
continue
# 显示详细统计信息
print(f"\n📊 详细交叉验证统计结果:")
print(f"{'模型':<20}{'均值':<8}{'标准差':<8}{'中位数':<8}{'最小值':<8}{'最大值':<8}{'变异系数':<10}")
print("-" * 85)
for model_name, stats_dict in cv_statistics.items():
print(f"{model_name:<20}{stats_dict['mean']:<8.4f}{stats_dict['std']:<8.4f} "
f"{stats_dict['median']:<8.4f}{stats_dict['min']:<8.4f}{stats_dict['max']:<8.4f} "
f"{stats_dict['cv_coefficient']:<10.4f}")
阶段 3: 正态性检验 (Normality Test)
在选择统计检验方法之前，我们需要了解数据的分布特性。这个阶段使用 Shapiro-Wilk 检验来判断每个模型收集到的10个交叉验证分数是否符合正态分布。
作用解释: 这是统计分析中的一个关键步骤。像“配对t检验”这样的参数检验方法，其有效性的前提是数据（或其差异）来自正态分布。如果数据不符合正态分布，那么t检验的结果可能不可靠，我们就需要选择非参数检验方法（如 Wilcoxon 符号秩检验）。此阶段为后续选择正确的检验方法提供了依据。
# ===== 2.1.2 正态性检验 =====
print(f"\n🔬 2.1.2 正态性检验")
print("-" * 30)
normality_results = {}
print(f"{'模型':<20}{'Shapiro-Wilk':<15}{'p-value':<12}{'正态分布':<10}")
print("-" * 60)
for model_name, scores in detailed_cv_scores.items():
# Shapiro-Wilk检验
stat, p_value = stats.shapiro(scores)
is_normal = p_value > 0.05
normality_results[model_name] = {'statistic': stat, 'p_value': p_value, 'is_normal': is_normal}
status = "✅ 是"if is_normal else"❌ 否"
print(f"{model_name:<20}{stat:<15.4f}{p_value:<12.4e}{status:<10}")
阶段 4: 配对t检验 (Paired t-test)
这个阶段执行核心的统计比较。它将性能最好的模型（基准模型）与其余四个模型进行逐一的配对t检验。
作用解释: 配对t检验的目的是判断两个模型的性能均值差异是否具有统计显著性。换句话说，最佳模型的性能优势是真实存在的，还是仅仅由数据抽样的随机性造成的？检验会给出一个p值，如果p值小于一个阈值（通常是0.05），我们就可以认为两个模型之间存在显著差异。此外，代码还计算了科恩 d (Cohen's d) 效应量，它衡量了差异的大小，这比单纯的“是/否”显著性判断提供了更丰富的信息。
# ===== 2.1.3 配对t检验 =====
print(f"\n📈 2.1.3 配对t检验 (Paired t-test)")
print("-" * 40)
# 获取最佳模型（排名第一）
best_model = top_5_models_names[0]
print(f"🏆 基准模型（最佳）: {best_model}")
# 进行配对t检验
ttest_results = {}
print(f"\n配对t检验结果:")
print(f"{'模型对比':<30}{'t统计量':<12}{'p-value':<12}{'显著性':<10}{'效应量(d)':<12}")
print("-" * 80)
for i, model_name inenumerate(top_5_models_names[1:], 1):
if model_name in detailed_cv_scores:
# 配对t检验
t_stat, p_value = ttest_rel(detailed_cv_scores[best_model], detailed_cv_scores[model_name])
# 计算Cohen's d效应量
mean_diff = detailed_cv_scores[best_model].mean() - detailed_cv_scores[model_name].mean()
pooled_std = np.sqrt((detailed_cv_scores[best_model].var() + detailed_cv_scores[model_name].var()) / 2)
cohens_d = mean_diff / pooled_std if pooled_std != 0else0
# 判断显著性
is_significant = p_value < 0.05
significance = "🔴 显著"if is_significant else"⚪ 不显著"
ttest_results[f"{best_model} vs {model_name}"] = {
't_statistic': t_stat,
'p_value': p_value,
'is_significant': is_significant,
'cohens_d': cohens_d,
'mean_difference': mean_diff
}
comparison_name = f"{best_model} vs {model_name}"
print(f"{comparison_name:<30}{t_stat:<12.4f}{p_value:<12.4e}{significance:<10}{cohens_d:<12.4f}")
阶段 5: Wilcoxon符号秩检验
此阶段进行 Wilcoxon 符号秩检验，它是一种非参数检验方法。
作用解释: 这是配对t检验的备用方案。如果之前的正态性检验显示数据不符合正态分布，那么 Wilcoxon 检验的结果会比 t检验更可靠。它不依赖于数据的分布假设，而是通过比较得分对的秩次来判断是否存在差异。同时进行两种检验并比较它们的结果，可以使我们的结论更加稳健。
# ===== 2.1.4 Wilcoxon符号秩检验 =====
print(f"\n📊 2.1.4 Wilcoxon符号秩检验 (非参数检验)")
print("-" * 45)
wilcoxon_results = {}
print(f"Wilcoxon符号秩检验结果:")
print(f"{'模型对比':<30}{'W统计量':<12}{'p-value':<12}{'显著性':<10}")
print("-" * 70)
for i, model_name inenumerate(top_5_models_names[1:], 1):
if model_name in detailed_cv_scores:
# Wilcoxon符号秩检验
try:
w_stat, p_value = wilcoxon(detailed_cv_scores[best_model], detailed_cv_scores[model_name])
is_significant = p_value < 0.05
significance = "🔴 显著"if is_significant else"⚪ 不显著"
wilcoxon_results[f"{best_model} vs {model_name}"] = {
'w_statistic': w_stat,
'p_value': p_value,
'is_significant': is_significant
}
comparison_name = f"{best_model} vs {model_name}"
print(f"{comparison_name:<30}{w_stat:<12.4f}{p_value:<12.4e}{significance:<10}")
except Exception as e:
print(f"{model_name:<30}{'Error':<12}{'---':<12}{'❌ 失败':<10}")
阶段 6: 多重比较校正 (Multiple Comparison Correction)
由于我们进行了多次统计检验（将最佳模型与其余4个模型比较），这个阶段对得到的p值进行校正。这里使用的是 Bonferroni 校正方法。
作用解释: 当同时进行多次统计检验时，偶然得到“显著”结果（即假阳性）的概率会增加。例如，在0.05的显著性水平下，进行20次检验，你可能期望有1次是由于纯粹的偶然性而显示为“显著”。多重比较校正（如 Bonferroni）通过提高判断显著性的门槛（实际上是调整p值）来控制这种错误的累积。这使得我们的结论更加保守和可信。
# ===== 2.1.5 多重比较校正 =====
print(f"\n🔧 2.1.5 多重比较校正 (Bonferroni校正)")
print("-" * 40)
# 收集所有p值进行多重比较校正
t_test_p_values = [result['p_value'] for result in ttest_results.values()]
wilcoxon_p_values = [result['p_value'] for result in wilcoxon_results.values()]
if t_test_p_values:
# Bonferroni校正
rejected_t, corrected_p_t, _, _ = multipletests(t_test_p_values, method='bonferroni')
print(f"配对t检验 - Bonferroni校正结果:")
print(f"{'模型对比':<30}{'原始p值':<12}{'校正p值':<12}{'校正后显著性':<15}")
print("-" * 75)
for i, (comparison, result) inenumerate(ttest_results.items()):
significance = "🔴 显著"if rejected_t[i] else"⚪ 不显著"
print(f"{comparison:<30}{result['p_value']:<12.4e}{corrected_p_t[i]:<12.4e}{significance:<15}")
if wilcoxon_p_values:
rejected_w, corrected_p_w, _, _ = multipletests(wilcoxon_p_values, method='bonferroni')
print(f"\nWilcoxon检验 - Bonferroni校正结果:")
print(f"{'模型对比':<30}{'原始p值':<12}{'校正p值':<12}{'校正后显著性':<15}")
print("-" * 75)
for i, (comparison, result) inenumerate(wilcoxon_results.items()):
significance = "🔴 显著"if rejected_w[i] else"⚪ 不显著"
print(f"{comparison:<30}{result['p_value']:<12.4e}{corrected_p_w[i]:<12.4e}{significance:<15}")
阶段 7: 可视化统计检验结果
这个阶段将前面分析得到的数值结果转换成直观的图表，包括箱线图、效应量条形图和p值热力图。
作用解释: 图表是沟通分析结果的最有效方式。"一图胜千言"，这些图表帮助我们：
- 箱线图 (Box Plot):
直观地比较不同模型性能得分的分布情况，包括中位数、波动范围和异常值。图上还标记了显著性水平（
*,**,***），将统计检验结果与数据分布结合起来。 - 效应量图 (Effect Size Plot):
展示了最佳模型与其他模型性能差异的“大小”，而不仅仅是“有无差异”。这有助于判断一个统计上显著的差异在实际应用中是否足够重要。
- p值热力图 (P-value Heatmap):
以矩阵的形式清晰地展示所有模型两两比较的p值，颜色深浅代表p值大小，使得快速识别哪些模型对之间存在显著差异变得非常容易。
# ===== 2.1.6 可视化统计检验结果 =====
print(f"\n📊 2.1.6 生成统计检验可视化")
print("-" * 35)
# 绘制交叉验证分数箱线图
plt.figure(figsize=(14, 8))
cv_data = [detailed_cv_scores[model] for model in top_5_models_names if model in detailed_cv_scores]
cv_labels = [model for model in top_5_models_names if model in detailed_cv_scores]
bp = plt.boxplot(cv_data, labels=cv_labels, patch_artist=True, showmeans=True)
# 美化箱线图
colors = ['lightblue', 'lightgreen', 'lightcoral', 'lightyellow', 'lightpink']
for patch, color inzip(bp['boxes'], colors[:len(bp['boxes'])]):
patch.set_facecolor(color)
plt.xlabel('Models', fontsize=12, fontweight='bold')
plt.ylabel('Cross-Validation R² Score', fontsize=12, fontweight='bold')
plt.title('Cross-Validation Score Distribution for Top 5 Models', fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(True, alpha=0.3)
# 添加统计显著性标记
for i, model inenumerate(cv_labels[1:], 1):
comparison_key = f"{best_model} vs {model}"
if comparison_key in ttest_results:
p_val = ttest_results[comparison_key]['p_value']
if p_val < 0.001:
significance_mark = "***"
elif p_val < 0.01:
significance_mark = "**"
elif p_val < 0.05:
significance_mark = "*"
else:
significance_mark = "ns"
# 在箱线图上方添加显著性标记
y_pos = max([max(scores) for scores in cv_data]) * (1.02 + i * 0.03)
plt.text(i + 1, y_pos, significance_mark, ha='center', va='bottom', fontweight='bold', fontsize=12)
plt.tight_layout()
plt.savefig(os.path.join(stats_test_folder, 'cv_scores_boxplot_with_significance.png'),
dpi=300, bbox_inches='tight')
plt.show()
# 绘制效应量图
if ttest_results:
plt.figure(figsize=(12, 6))
comparisons = list(ttest_results.keys())
effect_sizes = [ttest_results[comp]['cohens_d'] for comp in comparisons]
p_values = [ttest_results[comp]['p_value'] for comp in comparisons]
# 创建颜色映射
colors = ['red'if p < 0.05else'gray'for p in p_values]
bars = plt.bar(range(len(comparisons)), effect_sizes, color=colors, alpha=0.7)
# 添加效应量解释线
plt.axhline(y=0.2, color='green', linestyle='--', alpha=0.5, label='Small effect (d=0.2)')
plt.axhline(y=0.5, color='orange', linestyle='--', alpha=0.5, label='Medium effect (d=0.5)')
plt.axhline(y=0.8, color='red', linestyle='--', alpha=0.5, label='Large effect (d=0.8)')
plt.xlabel('Model Comparisons', fontsize=12, fontweight='bold')
plt.ylabel("Cohen's d (Effect Size)", fontsize=12, fontweight='bold')
plt.title("Effect Size Analysis for Model Comparisons", fontsize=14, fontweight='bold')
plt.xticks(range(len(comparisons)), [comp.replace(f"{best_model} vs ", "") for comp in comparisons],
rotation=45, ha='right')
plt.legend()
plt.grid(True, alpha=0.3)
# 添加数值标签
for i, (bar, effect_size) inenumerate(zip(bars, effect_sizes)):
plt.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.01,
f'{effect_size:.3f}', ha='center', va='bottom', fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(stats_test_folder, 'effect_size_analysis.png'),
dpi=300, bbox_inches='tight')
plt.show()
# 绘制p值热力图
iflen(top_5_models_names) > 1:
plt.figure(figsize=(10, 8))
# 创建p值矩阵
n_models = len(top_5_models_names)
p_value_matrix = np.ones((n_models, n_models))
for i, model1 inenumerate(top_5_models_names):
for j, model2 inenumerate(top_5_models_names):
if i != j and model1 in detailed_cv_scores and model2 in detailed_cv_scores:
iff"{model1} vs {model2}"in ttest_results:
p_value_matrix[i, j] = ttest_results[f"{model1} vs {model2}"]['p_value']
eliff"{model2} vs {model1}"in ttest_results:
p_value_matrix[i, j] = ttest_results[f"{model2} vs {model1}"]['p_value']
# 绘制热力图
mask = p_value_matrix == 1# 遮盖对角线
sns.heatmap(p_value_matrix, mask=mask, annot=True, fmt='.3e', cmap='RdYlBu_r',
xticklabels=top_5_models_names, yticklabels=top_5_models_names,
cbar_kws={'label': 'p-value'})
plt.title('P-value Matrix for Model Comparisons (t-test)', fontsize=14, fontweight='bold')
plt.xlabel('Models', fontsize=12)
plt.ylabel('Models', fontsize=12)
plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)
plt.tight_layout()
plt.savefig(os.path.join(stats_test_folder, 'pvalue_heatmap.png'), dpi=300, bbox_inches='tight')
plt.show()
第四部分：最优模型的贝叶斯优化
阶段一：选择最优模型
在进行精细化调优之前，我们首先需要确定要优化哪个模型。这个阶段的代码就是从之前模型竞赛的结果中，找出冠军模型。
# 选择最优模型并进行贝叶斯优化 =====
best_model_name = results_df.iloc[0]['Model']
print(f"\n最优模型为: {best_model_name}")
print(f"测试集 R² 分数: {results_df.iloc[0]['Test_R2']:.4f}")
# 获取最优模型
best_model = trained_models[best_model_name]
print(f"\n开始对最优模型 {best_model_name} 进行贝叶斯优化...")
代码作用解释：
-
从已排序的性能结果表
results_df中，获取第一行（最优）的模型名称。 -
根据名称，从预训练模型字典
trained_models中获取对应的模型对象。 -
打印出选定的最优模型及其原始性能，作为优化的起点。
阶段二：定义目标函数 (Objective Function)
这是贝叶斯优化的核心。我们需要定义一个“评估标准”，告诉优化工具（Optuna）如何判断一组超参数的好坏。这个函数接收一套由优化器“建议”的参数，运行模型，然后返回一个代表性能的分数。
# 定义贝叶斯优化的目标函数
defobjective(trial):
"""贝叶斯优化的目标函数"""
if best_model_name == 'XGBoost':
params = {
'n_estimators': trial.suggest_int('n_estimators', 100, 1000),
'max_depth': trial.suggest_int('max_depth', 3, 15),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
'subsample': trial.suggest_float('subsample', 0.6, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
'reg_alpha': trial.suggest_float('reg_alpha', 0, 10),
'reg_lambda': trial.suggest_float('reg_lambda', 0, 10),
'random_state': 0
}
model = xgb.XGBRegressor(**params)
elif best_model_name == 'LightGBM':
params = {
'n_estimators': trial.suggest_int('n_estimators', 100, 1000),
'max_depth': trial.suggest_int('max_depth', 3, 15),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
'subsample': trial.suggest_float('subsample', 0.6, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
'reg_alpha': trial.suggest_float('reg_alpha', 0, 10),
'reg_lambda': trial.suggest_float('reg_lambda', 0, 10),
'random_state': 0,
'verbose': -1
}
model = lgb.LGBMRegressor(**params)
elif best_model_name == 'Random Forest':
params = {
'n_estimators': trial.suggest_int('n_estimators', 50, 500),
'max_depth': trial.suggest_int('max_depth', 3, 20),
'min_samples_split': trial.suggest_int('min_samples_split', 2, 20),
'min_samples_leaf': trial.suggest_int('min_samples_leaf', 1, 10),
'max_features': trial.suggest_categorical('max_features', ['sqrt', 'log2', None]),
'random_state': 0
}
model = RandomForestRegressor(**params)
elif best_model_name == 'Support Vector Machine':
params = {
'C': trial.suggest_float('C', 0.1, 100, log=True),
'gamma': trial.suggest_categorical('gamma', ['scale', 'auto']),
'kernel': trial.suggest_categorical('kernel', ['rbf', 'linear', 'poly']),
'epsilon': trial.suggest_float('epsilon', 0.01, 1.0)
}
model = SVR(**params)
elif best_model_name == 'Gradient Boosting':
params = {
'n_estimators': trial.suggest_int('n_estimators', 50, 500),
'max_depth': trial.suggest_int('max_depth', 3, 15),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
'subsample': trial.suggest_float('subsample', 0.6, 1.0),
'random_state': 0
}
model = GradientBoostingRegressor(**params)
else:
# 对于其他模型，返回原始性能
return results_df[results_df['Model'] == best_model_name]['CV_R2_Mean'].values[0]
# 使用交叉验证评估模型性能
cv_scores = cross_val_score(model, X_train, y_train, cv=5, scoring='r2')
return cv_scores.mean()
代码作用解释：
-
定义一个
objective函数，供 Optuna 调用以评估每组超参数。 -
函数内部，使用
trial.suggest_*方法为不同模型的超参数定义搜索范围。 -
使用5折交叉验证（
cross_val_score）来计算当前参数组合的平均R²分数，以获得稳定的性能评估。 -
返回这个平均分数，Optuna将尝试最大化该值。
阶段三：执行优化过程
定义好蓝图后，现在开始“施工”。我们创建一个Optuna的“研究”（Study），并让它执行指定次数的试验。
# 运行贝叶斯优化
print("开始贝叶斯优化超参数搜索...")
study = optuna.create_study(direction='maximize', sampler=TPESampler(seed=0))
study.optimize(objective, n_trials=100, show_progress_bar=True)
print("\n贝叶斯优化完成！")
print(f"最佳参数: {study.best_params}")
print(f"最佳交叉验证 R² 分数: {study.best_value:.4f}")
代码作用解释：
-
创建一个Optuna的
study对象，设定优化方向为最大化 (maximize)。 -
调用
study.optimize()方法，让它执行objective函数100次，自动搜索最佳参数。 -
优化完成后，从
study对象中提取并打印找到的最佳参数 (best_params) 和最高分 (best_value)。
阶段四：最终模型的训练与评估
优化过程找到了最佳“配方”，但还没有做出最终的“蛋糕”。这个阶段，我们用这个最佳配方，在所有训练数据上训练一个最终模型，并用独立的测试集来检验其真实水平。
# 使用最佳参数重新训练模型
if best_model_name == 'XGBoost':
optimized_model = xgb.XGBRegressor(**study.best_params)
elif best_model_name == 'LightGBM':
optimized_model = lgb.LGBMRegressor(**study.best_params)
elif best_model_name == 'Random Forest':
optimized_model = RandomForestRegressor(**study.best_params)
elif best_model_name == 'Support Vector Machine':
optimized_model = SVR(**study.best_params)
elif best_model_name == 'Gradient Boosting':
optimized_model = GradientBoostingRegressor(**study.best_params)
else:
optimized_model = best_model
# 训练优化后的模型
optimized_model.fit(X_train, y_train)
# 评估优化后的模型
y_pred_train_opt = optimized_model.predict(X_train)
y_pred_test_opt = optimized_model.predict(X_test)
opt_train_r2 = r2_score(y_train, y_pred_train_opt)
opt_test_r2 = r2_score(y_test, y_pred_test_opt)
opt_test_rmse = np.sqrt(mean_squared_error(y_test, y_pred_test_opt))
opt_test_mae = mean_absolute_error(y_test, y_pred_test_opt)
代码作用解释：
-
使用找到的最佳参数
study.best_params来创建一个新的、最终的模型实例。 -
在 全部 训练数据上训练这个优化后的模型 (
optimized_model.fit)。 -
使用训练好的模型对训练集和测试集进行预测。
-
计算模型在未见过的测试集上的最终性能指标（R², RMSE, MAE）。
阶段五：展示优化成果——性能对比
为了证明我们的努力没有白费，这个阶段将优化前后的性能指标进行并排比较，效果一目了然。
print(f"\n优化前后性能对比 ({best_model_name}):")
print("=" * 50)
print(f"{'指标':<15}{'优化前':<12}{'优化后':<12}{'提升':<10}")
print("=" * 50)
original_test_r2 = results_df[results_df['Model'] == best_model_name]['Test_R2'].values[0]
original_test_rmse = results_df[results_df['Model'] == best_model_name]['Test_RMSE'].values[0]
original_test_mae = results_df[results_df['Model'] == best_model_name]['Test_MAE'].values[0]
print(
f"{'测试集 R²':<15}{original_test_r2:<12.4f}{opt_test_r2:<12.4f}{((opt_test_r2 - original_test_r2) / original_test_r2 * 100):>+7.2f}%")
print(
f"{'测试集 RMSE':<15}{original_test_rmse:<12.4f}{opt_test_rmse:<12.4f}{((opt_test_rmse - original_test_rmse) / original_test_rmse * 100):>+7.2f}%")
print(f"\n=== 任务完成 ===")
print(f"最优模型: {best_model_name}")
print(f"最终测试集 R² 分数: {opt_test_r2:.4f}")
print(f"最终测试集 RMSE: {opt_test_rmse:.4f}")
代码作用解释
-
从原始结果表
results_df中获取优化前的性能分数。 -
以表格形式清晰地打印优化前、优化后的性能指标，并计算提升百分比。
-
最后，总结并报告最终模型的性能。
阶段六：结果的可视化与保存
最后一步，我们将模型的预测效果用图形化方式展示出来，并将这张图和最终模型保存下来，以便报告和未来使用。
print("\n正在绘制训练集与验证集的回归拟合图...") # 打印开始绘制拟合图的提示
y_train_pred = y_pred_train_opt # 使用模型对训练集进行预测
y_test_pred = y_pred_test_opt # 使用模型对测试集进行预测
# 计算训练集的指标
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))
mae_train = mean_absolute_error(y_train, y_train_pred)
# 计算验证集的指标
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
mae_test = mean_absolute_error(y_test, y_test_pred)
print(f"\n训练集评估指标: R2={r2_train:.4f}, RMSE={rmse_train:.4f}, MAE={mae_train:.4f}")
print(f"验证集评估指标: R2={r2_test:.4f}, RMSE={rmse_test:.4f}, MAE={mae_test:.4f}")
plt.figure(figsize=(8, 8), dpi=150) # 创建一个新的图形窗口，并设置大小和分辨率
plt.scatter(y_train, y_train_pred, alpha=0.5, label='Train', color='blue') # 绘制训练集的真实值与预测值的散点图
plt.scatter(y_test, y_test_pred, alpha=0.7, label='Validation', color='red', marker='^') # 绘制测试集的真实值与预测值的散点图
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=2, label='1：1 (y=x)') # 绘制一条y=x的参考线
plt.xlabel('Actual Values', fontsize=12) # 设置X轴标签
plt.ylabel('Predicted Values', fontsize=12) # 设置Y轴标签
#最优模型拟合
plt.title(f'{best_model_name} Regression Fit Plot', fontsize=14, fontweight='bold') # 设置图表标题
plt.legend(loc='upper left') # 显示图例，并放在左上角
plt.grid(True) # 显示网格线
# 创建一个将在图上显示的文本框
metrics_text = (
f'test dataset:\n'
f'R^2 = {r2_test:.4f}\n'
f'RMSE = {rmse_test:.4f}\n'
f'MAE = {mae_test:.4f}\n\n'
f'train dataset:\n'
f'R^2 = {r2_train:.4f}\n'
f'RMSE = {rmse_train:.4f}\n'
f'MAE = {mae_train:.4f}')
# 将文本框放置在图的右下角
plt.text(0.95, 0.05, metrics_text, transform=plt.gca().transAxes, fontsize=10,
verticalalignment='bottom', horizontalalignment='right',
bbox=dict(boxstyle='round,pad=0.5', fc='wheat', alpha=0.5))
# 定义当前路径是主输出文件夹的路径
current_path = os.getcwd() # 获取当前脚本所在的目录路径
output_main_folder = os.path.join(current_path, 'output') # 定义主输出文件
os.makedirs(output_main_folder, exist_ok=True) # 创建主输出文件夹，如果已存在则不报错
fit_plot_path = os.path.join(output_main_folder, 'regression_fit_plot.png') # 拼接拟合图的完整保存路径
plt.savefig(fit_plot_path, bbox_inches='tight') # 保存图形，bbox_inches='tight'可以裁剪掉空白边缘
print(f"回归拟合图已保存至: '{fit_plot_path}'") # 打印保存路径
plt.show() # 显示图形
model= optimized_model # 将优化后的模型赋值给变量model，方便后续使用
代码作用解释
-
绘制一张“真实值 vs. 预测值”的散点图，以可视化模型性能。
-
图中包含一条
y=x的完美预测参考线，以及训练集（蓝点）和测试集（红三角）的预测点。 -
在图的右下角添加一个文本框，显示详细的性能指标。
-
将生成的图像保存到本地的
output文件夹中，并在屏幕上显示出来。 -
将最终的优化模型赋值给变量
model，以便后续直接使用。
第五部分：最优模型的shap高级可视化
阶段 1: 初始化与设置
这个阶段是所有增强版SHAP可视化的准备工作。它首先打印一个清晰的标题，然后创建一个新的子文件夹 enhanced_shap_visualizations，专门用于存放接下来生成的所有高级图表和报告。
作用解释: 此阶段的代码确保了分析环境的整洁。通过创建一个独立的文件夹，所有增强版的可视化结果都会被集中管理，与之前生成的基础图表和报告分开，使得文件结构更加清晰，方便用户查找特定类型的分析结果。
print("\n=== 开始增强版SHAP可视化分析 ===")
# 创建增强版可视化的输出文件夹
enhanced_shap_folder = os.path.join(output_main_folder, 'enhanced_shap_visualizations')
os.makedirs(enhanced_shap_folder, exist_ok=True)
print(f"增强版SHAP可视化将保存到: '{enhanced_shap_folder}'")
阶段 2: 生成SHAP决策图 (Decision Plot)
此阶段生成SHAP决策图，这是一种可以同时可视化多个样本预测路径的强大工具。
作用解释: 决策图从底部（模型的基准预测值）开始，向上逐一展示每个特征是如何将预测值“推”向最终结果的。每一条线代表一个样本的预测过程。这种图的优势在于能够清晰地展示样本群体的总体趋势，以及个体样本之间的差异。通过观察线条的走向和聚集情况，我们可以理解模型在做出相似或不同预测时，是否依赖于相同的特征模式。
# 7.1 SHAP决策图（Decision Plot）
print("\n7.1 正在生成SHAP决策图...")
try:
# 选择前10个最重要的特征和前50个样本进行决策图可视化
n_samples_decision = min(50, len(X_test))
n_features_decision = min(10, len(sorted_features))
# 创建决策图
plt.figure(figsize=(12, 8))
shap.decision_plot(
base_value=shap_values_obj.base_values[0],
shap_values=shap_values_sorted[:n_samples_decision, :n_features_decision],
features=X_test_sorted.iloc[:n_samples_decision, :n_features_decision],
feature_names=sorted_features[:n_features_decision],
show=False
)
plt.title('SHAP Decision Plot - Top 10 Features for 50 Samples', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_decision_plot.png'), dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP决策图已生成")
except Exception as e:
print(f" ✗ SHAP决策图生成失败: {e}")
阶段 3: 生成SHAP瀑布图 (Waterfall Plot)
这个阶段为单个样本生成SHAP瀑布图，它详细地分解了单次预测的构成。
作用解释: 瀑布图是解释单个预测最直观的方式之一。它从底部的模型基准值 E[f(x)] 开始，然后像搭积木一样，将每个特征的贡献（正向或负向）逐一累加，最终得到顶部的最终预测值 f(x)。红色箭头表示该特征“推高”了预测值，蓝色箭头表示“拉低”。通过这个图，我们可以精确地看到对于某一个特定的样本，哪些特征是其预测结果的主要驱动因素。
# 7.2 SHAP瀑布图（Waterfall Plot）
print("\n7.2 正在生成SHAP瀑布图...")
try:
# 为前3个样本分别生成瀑布图
for sample_idx inrange(min(3, len(X_test))):
plt.figure(figsize=(10, 6))
shap.plots.waterfall(
shap.Explanation(
values=shap_values_sorted[sample_idx],
base_values=shap_values_obj.base_values[0],
data=X_test_sorted.iloc[sample_idx],
feature_names=sorted_features
),
show=False
)
plt.title(f'SHAP Waterfall Plot - Sample {sample_idx + 1}', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, f'shap_waterfall_plot_sample_{sample_idx + 1}.png'),
dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP瀑布图已生成")
except Exception as e:
print(f" ✗ SHAP瀑布图生成失败: {e}")
阶段 4: 生成SHAP力图 (Force Plot)
此阶段为单个样本生成力图，这是另一种解释单次预测的强大可视化。
作用解释: 力图将所有特征的贡献想象成作用在一个平衡点（基准值）上的“力”。推高预测值的特征（红色）向右推，拉低预测值的特征（蓝色）向左推。所有力作用的结果，就是最终的预测值。这个图的优点是它非常紧凑，能在一行内清晰地展示出正负贡献力量的博弈过程，并直观地显示出哪些特征的“力量”最大。
# 7.3 SHAP力图（Force Plot）
print("\n7.3 正在生成SHAP力图...")
try:
# 为单个样本生成力图
for sample_idx inrange(min(3, len(X_test))):
# 四舍五入处理，保留两位小数
expected_value_rounded = np.round(explainer.expected_value, 2)
shap_values_rounded = np.round(shap_values[sample_idx], 2)
feature_values_rounded = X_test_sorted.iloc[sample_idx].round(2)
shap_force_plot = shap.force_plot(
expected_value_rounded,
shap_values_rounded,
feature_values_rounded,
feature_names=sorted_features,
matplotlib=True,
show=False
)
plt.title(f'SHAP Force Plot - Sample {sample_idx + 1}', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, f'shap_force_plot_sample_{sample_idx + 1}.png'),
dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP力图已生成")
except Exception as e:
print(f" ✗ SHAP力图生成失败: {e}")
阶段 5: 生成SHAP热力图 (Heatmap)
这个阶段生成SHAP热力图，用于展示特征重要性在不同样本间的分布和模式。
作用解释: 热力图的x轴是特征，y轴是不同的样本。图中每个单元格的颜色代表了该特征对该样本的SHAP值大小（通常红色表示正贡献，蓝色表示负贡献）。通过观察热力图，我们可以发现一些有趣的模式，例如：
- 某些特征是否对所有样本都有相似的影响？
（整列颜色一致）
- 是否存在样本集群，它们对特征的反应模式相似？
（几行颜色模式相似）
- 特征重要性是否集中在少数几个样本上？
这种图为我们提供了一个从样本和特征两个维度同时审视模型行为的视角。
# 7.4 SHAP热力图
print("\n7.4 正在生成SHAP热力图...")
try:
# 选择前20个样本和前15个特征
n_samples_heatmap = min(20, len(X_test))
n_features_heatmap = min(15, len(sorted_features))
plt.figure(figsize=(12, 8))
shap.plots.heatmap(
shap.Explanation(
values=shap_values_sorted[:n_samples_heatmap, :n_features_heatmap],
data=X_test_sorted.iloc[:n_samples_heatmap, :n_features_heatmap].values,
feature_names=sorted_features[:n_features_heatmap]
),
show=False
)
plt.title('SHAP Heatmap - Feature Importance Across Samples', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_heatmap.png'), dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP热力图已生成")
except Exception as e:
print(f" ✗ SHAP热力图生成失败: {e}")
阶段 6: 进行SHAP特征聚类分析
这个阶段利用SHAP值对特征进行聚类，以发现哪些特征对模型的影响模式相似。
作用解释: 传统的特征聚类通常基于特征本身的数值相关性，而这里的聚类是基于特征对模型预测的贡献模式。代码使用层次聚类算法，将平均SHAP重要性相似的特征聚集在一起，并用树状图（Dendrogram）进行可视化。通过这个图，我们可以识别出冗余的特征（即对模型有相似作用的特征），这对于特征选择和简化模型非常有帮助。
# 7.6 特征聚类分析
print("\n7.6 正在进行SHAP特征聚类分析...")
try:
from sklearn.cluster import KMeans
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.metrics import silhouette_score
# 使用SHAP值进行特征聚类
shap_importance = np.abs(shap_values_sorted).mean(axis=0)
# 层次聚类
plt.figure(figsize=(12, 8))
linkage_matrix = linkage(shap_importance.reshape(-1, 1), method='ward')
dendrogram(linkage_matrix, labels=sorted_features, orientation='left')
plt.title('Feature Clustering Based on SHAP Importance', fontsize=14, fontweight='bold')
plt.xlabel('Distance')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_feature_clustering.png'), dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP特征聚类分析已完成")
except Exception as e:
print(f" ✗ SHAP特征聚类分析失败: {e}")
阶段 7: 生成SHAP值分布分析
此阶段为最重要的几个特征绘制其SHAP值的分布直方图。
作用解释: 与依赖图不同，这个分析不关心SHAP值与特征值本身的关系，而是纯粹关注SHAP值的分布形态。通过直方图，我们可以看出：
-
一个特征的贡献是单峰的（影响模式单一）还是多峰的（可能存在多种影响模式）？
-
分布是对称的，还是偏斜的（例如，该特征主要起正向作用，很少起负向作用）？
-
分布的离散程度如何（影响的强度是否稳定）？
这为深入理解单个特征的行为提供了另一个维度的信息。
# 7.7 SHAP值分布分析
print("\n7.7 正在生成SHAP值分布分析...")
try:
# 为前6个最重要的特征绘制SHAP值分布
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
for i, feature_name inenumerate(sorted_features[:6]):
ax = axes[i]
shap_vals = shap_values_sorted[:, i]
# 绘制直方图
ax.hist(shap_vals, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
ax.axvline(shap_vals.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {shap_vals.mean():.3f}')
ax.axvline(0, color='black', linestyle='-', linewidth=1, alpha=0.5)
ax.set_title(f'{feature_name}', fontsize=12, fontweight='bold')
ax.set_xlabel('SHAP Value')
ax.set_ylabel('Frequency')
ax.legend()
ax.grid(True, alpha=0.3)
plt.suptitle('SHAP Value Distributions for Top 6 Features', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_value_distributions.png'), dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP值分布分析已完成")
except Exception as e:
print(f" ✗ SHAP值分布分析失败: {e}")
阶段 8: 生成SHAP特征相关性矩阵
这个阶段计算并可视化特征SHAP值之间的相关性矩阵。
作用解释: 这与传统的特征相关性矩阵有本质区别。传统方法衡量的是特征原始数值的线性关系，而这里衡量的是特征对模型预测的贡献之间的线性关系。如果两个特征的SHAP值高度正相关，意味着当一个特征对预测产生正向贡献时，另一个特征也倾向于产生正向贡献（反之亦然）。这可以帮助我们发现模型在内部是如何“组合”使用特征的，揭示出比原始数据相关性更深层次的关联模式。
# 7.8 SHAP特征相关性矩阵
print("\n7.8 正在生成SHAP特征相关性矩阵...")
try:
# 计算SHAP值的相关性矩阵
shap_corr = np.corrcoef(shap_values_sorted.T)
plt.figure(figsize=(12, 10))
mask = np.triu(np.ones_like(shap_corr, dtype=bool))
sns.heatmap(shap_corr, mask=mask, annot=True, cmap='coolwarm', center=0,
square=True, fmt='.2f', cbar_kws={"shrink": .5},
xticklabels=sorted_features, yticklabels=sorted_features)
plt.title('SHAP Values Correlation Matrix', fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_correlation_matrix.png'), dpi=300, bbox_inches='tight')
plt.show()
print(" ✓ SHAP特征相关性矩阵已生成")
except Exception as e:
print(f" ✗ SHAP特征相关性矩阵生成失败: {e}")
阶段 9: 进行SHAP稳定性分析
此阶段通过Bootstrap（自助法）重采样技术来评估特征重要性（平均绝对SHAP值）的稳定性。
作用解释: 模型的训练结果会受到训练数据随机性的影响。为了检验我们的特征重要性排名是否稳健，此阶段的代码多次从测试集中有放回地抽取80%的样本，并重新计算SHAP重要性。最后，它会计算每个特征重要性的均值和标准差。
- 误差棒图
直观地展示了每个特征重要性的波动范围。误差棒越短，说明该特征的重要性越稳定。
- 变异系数 (CV)
是一个标准化的稳定性度量（标准差/均值），CV值越小，特征越稳定。
这个分析帮助我们识别那些在不同数据子集上始终保持重要的核心特征，以及那些重要性波动较大的、可能不太可靠的特征。
# 7.9 SHAP稳定性分析
print("\n7.9 正在进行SHAP稳定性分析...")
try:
# 使用bootstrap方法分析SHAP值的稳定性
n_bootstrap = 50
bootstrap_means = []
for _ inrange(n_bootstrap):
# 随机抽样
sample_indices = np.random.choice(len(X_test), size=int(0.8 * len(X_test)), replace=True)
bootstrap_shap = shap_values_sorted[sample_indices]
bootstrap_mean = np.abs(bootstrap_shap).mean(axis=0)
bootstrap_means.append(bootstrap_mean)
bootstrap_means = np.array(bootstrap_means)
# 计算均值和标准差
mean_importance = bootstrap_means.mean(axis=0)
std_importance = bootstrap_means.std(axis=0)
cv_importance = std_importance / mean_importance # 变异系数
# 绘制稳定性图
plt.figure(figsize=(12, 8))
x_pos = np.arange(len(sorted_features))
plt.errorbar(x_pos, mean_importance, yerr=std_importance, fmt='o', capsize=5, capthick=2)
plt.xticks(x_pos, sorted_features, rotation=45, ha='right')
plt.xlabel('Features')
plt.ylabel('Mean Absolute SHAP Value')
plt.title('SHAP Value Stability Analysis (Bootstrap)', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(enhanced_shap_folder, 'shap_stability_analysis.png'), dpi=300, bbox_inches='tight')
plt.show()
# 保存稳定性数据
stability_df = pd.DataFrame({
'Feature': sorted_features,
'Mean_SHAP': mean_importance,
'Std_SHAP': std_importance,
'CV_SHAP': cv_importance
}).sort_values('CV_SHAP')
stability_df.to_csv(os.path.join(enhanced_shap_folder, 'shap_stability_analysis.csv'),
index=False, encoding='utf-8-sig')
print(" ✓ SHAP稳定性分析已完成")
print(f" 前5个最稳定的特征（按变异系数排序）:")
for i, row in stability_df.head().iterrows():
print(f" {row['Feature']}: CV = {row['CV_SHAP']:.4f}")
except Exception as e:
print(f" ✗ SHAP稳定性分析失败: {e}")
阶段 10: 生成综合SHAP分析报告
这是最后一个阶段，负责将所有数值化的SHAP分析结果汇总到一个易于阅读的文本报告中。
作用解释: 图表非常直观，但一个结构化的文本报告同样重要。此阶段生成的 comprehensive_shap_report.txt 文件，系统地总结了：
-
模型的基本信息。
-
最终的特征重要性排名及其数值。
-
SHAP值的关键统计数据（最大、最小、标准差等）。
-
每个特征是正向贡献还是负向贡献的比例。
-
基于以上分析得出的初步结论和建议。
这个报告为整个SHAP分析提供了一个可存档、可分享的书面总结，方便他人快速理解分析的核心发现。
# 7.10 综合SHAP报告
print("\n7.10 正在生成综合SHAP分析报告...")
try:
# 创建综合报告
report_path = os.path.join(enhanced_shap_folder, 'comprehensive_shap_report.txt')
withopen(report_path, 'w', encoding='utf-8') as f:
f.write("=== 综合SHAP分析报告 ===\n")
f.write("=" * 50 + "\n\n")
f.write("1. 模型基本信息:\n")
f.write(f" - 模型类型: {best_model_name}\n")
f.write(f" - 测试样本数: {len(X_test)}\n")
f.write(f" - 特征数量: {len(sorted_features)}\n")
f.write(f" - 模型基准值: {shap_values_obj.base_values[0]:.4f}\n\n")
f.write("2. 特征重要性排名（按平均绝对SHAP值）:\n")
for i, (feature, importance) inenumerate(zip(sorted_features, shap_df["mean_shap"].values)):
f.write(f" {i + 1:2d}. {feature:20s}: {importance:.6f}\n")
f.write("\n")
f.write("3. SHAP值统计信息:\n")
f.write(f" - 最大SHAP值: {shap_values_sorted.max():.6f}\n")
f.write(f" - 最小SHAP值: {shap_values_sorted.min():.6f}\n")
f.write(f" - SHAP值总体标准差: {shap_values_sorted.std():.6f}\n")
f.write(f" - 平均绝对SHAP值: {np.abs(shap_values_sorted).mean():.6f}\n\n")
f.write("4. 特征贡献度分析:\n")
positive_contrib = (shap_values_sorted > 0).sum(axis=0) / len(shap_values_sorted)
for i, (feature, pos_ratio) inenumerate(zip(sorted_features, positive_contrib)):
f.write(f" {feature:20s}: {pos_ratio * 100:5.1f}% 样本为正贡献\n")
f.write("\n")
f.write("5. 分析建议:\n")
f.write(" a) 关注前5个最重要的特征进行模型优化\n")
f.write(" b) 对于贡献度方向不一致的特征，建议进一步分析其条件依赖关系\n")
f.write(" c) 可以考虑特征工程来提升模型稳定性\n")
f.write(" d) 建议定期重新评估特征重要性以确保模型的时效性\n")
print(f" ✓ 综合SHAP分析报告已保存至: {report_path}")
except Exception as e:
print(f" ✗ 综合SHAP分析报告生成失败: {e}")
该文章案例
数据请加微信免费获取。
注：本代码全程Python语言实现，拿到代码后，先用示例数据复现跑通，确认环境没问题后，再上自己的数据。
【数据，请加微信免费获取】
如果你对类似于这样的文章感兴趣。
欢迎关注、点赞、转发
更多推荐
17
0- 0
已为社区贡献9条内容
所有评论(0)