期刊图片复现|Python回归分析全流程深度解析-回归拟合图、SHAP特征重要性总览图、交互作用强度气泡图、单特征依赖图、双特征交互效应图

# =========================================================================================# ====================================== 1. 环境设置 =======================================# =========================================================================================import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport osfrom statsmodels.nonparametric.smoothers_lowess import lowess

# =========================================================================================# ======================================2.颜色库=======================================# =========================================================================================COLOR_SCHEMES = {    0: {'train': 'blue', 'test': 'red', 'hist': '#4B0082', 'shap_scatter': '#00008B', 'lowess': '#9400D3',        'ci': '#D3D3D3', 'inter_low': 'blue', 'inter_low_fit': 'darkblue', 'inter_high': 'red',}

# =========================================================================================# ======================================3.边框与刻度设置辅助函数=======================================# =========================================================================================def apply_plot_styles(ax, tick_width=2, tick_length=5, spine_width=2, label_size=18):    for spine in ax.spines.values():  #遍历边框        spine.set_linewidth(spine_width)  #设置边框线宽    #设置刻度    ax.tick_params(axis='both',  #轴                   which='major',  #主刻度                   direction='in',  #朝向                   width=tick_width,  #粗细                   length=tick_length,#长                   labelsize=label_size)#字体大小

# =========================================================================================# ======================================4.结果图保存辅助函数======================================# =========================================================================================def save_fig_dual(fig, folder, filename_base, dpi=300):    os.makedirs(folder, exist_ok=True)  #如果文件夹不存在则创建    png_path = os.path.join(folder, f"{filename_base}.png")  #拼接保存路径    pdf_path = os.path.join(folder, f"{filename_base}.pdf")    fig.savefig(png_path, dpi=dpi, bbox_inches='tight')  #保存    fig.savefig(pdf_path, dpi=dpi, bbox_inches='tight')

# =========================================================================================# ======================================5.拟合线和置信区间计算函数=======================================# =========================================================================================#使用LOWESS拟合数据，通过Bootstrap生成大量拟合线，从而计算出95%的置信区间def bootstrap_lowess_ci(x, y, n_boot=200, frac=0.5, ci_level=0.95):    #少于10个，不进行拟合    if len(x) < 10:        return None, None    boot_lines = []  #用于存放每次重采样拟合出来的曲线    x_range = np.linspace(x.min(), x.max(), 100)  #在x的极值范围内均匀生成100个点，用作统一的插值横坐标    boot_lines_arr = np.array(boot_lines)  #将重采样曲线列表转换为NumPy数组，便于计算分位数    alpha = (1 - ci_level) / 2  #根据置信水平计算单侧分位数的概率值    lower_bound = np.quantile(boot_lines_arr, alpha, axis=0)  #置信下界    upper_bound = np.quantile(boot_lines_arr, 1 - alpha, axis=0)  #置信上界    return main_smoothed, (x_range, lower_bound, upper_bound)  #拟合曲线、x范围、置信上下界

# =========================================================================================# ======================================6.阈值点寻找函数=======================================# =========================================================================================#寻找曲线穿过y=0的所有交点/X坐标。通过检测y值正负号的变化，利用两点线性插值法精确计算出过零点的x坐标def find_roots(x_curve, y_curve):    roots = []  #存放根/零点的列表    sign_changes = np.where(np.diff(np.sign(y_curve)))[0]  #计算y值的符号差，找出正负号发生变化的相邻点索引    return roots  # 返回所有找到的零点列表

# =========================================================================================# ======================================7.阈值点绘制以及标注函数=======================================# =========================================================================================#算出零点，还在图表上画出垂直虚线，标上数值标签def find_and_plot_crossings(ax, x_curve, y_curve, color, x_range):    sign_changes = np.where(np.diff(np.sign(y_curve)))[0]  #找出y值正负值发生变化的位置索引    drawn_texts = []  #保存文本坐标        ax.text(x_root,  #x                y_pos,  #y                f' {x_root:.2f} ',  #文本                color='white',  #颜色                backgroundcolor=color,  #颜色                ha='center',  #水平                va='top',  #垂直                fontsize=18,  #字体大小                fontweight='bold',#加粗                bbox=dict(facecolor=color, edgecolor='none', pad=1),  #文本框                transform=ax.get_xaxis_transform())  #设置坐标变换        drawn_texts.append((x_root, y_pos))  #保存

# =========================================================================================# ======================================8.回归拟合图=======================================# =========================================================================================def plot_regression_results(metrics, colors, output_folder):    fig = plt.figure(figsize=(8.7, 8))  #创建画布    ax = plt.gca()  #获取当前画布坐标轴    all_y = pd.concat([metrics['train']['true'], metrics['test']['true']])  #合并训练集和测试集的真实值数据    #1:1线    plt.plot([all_y.min(),  #起点x              all_y.max()],  #终点x             [all_y.min(),  #起点y              all_y.max()],  #终点y             'k--',  #样式             lw=2,  #线宽             label='1:1 Line')  #图例名    # 设置边框和刻度    apply_plot_styles(ax)    # 保存    save_fig_dual(fig, output_folder, 'regression_fit_plot')    plt.close(fig)  #关闭

# =========================================================================================# ======================================9.特征重要性条形图与SHAP蜂巢图组合图绘制函数=======================================# =========================================================================================def plot_shap_summary(shap_values, X_test, feature_names, shap_df, base_values, colors,output_folder):    #创建画布    fig = plt.figure(figsize=(10, 10), dpi=300)    ax_sw = fig.add_axes([0.32, 0.11, 0.59, 0.77])  #定义主坐标轴在画布上的相对位置及大小    ax_bar = ax_sw.twiny()  #创建共享Y轴    y_pos = np.arange(len(feature_names))[::-1]  #生成逆序的Y轴坐标序列，使重要特征排在最上面    #绘制水平条形图    ax_bar.barh(y=y_pos,  #y                width=shap_df["mean_shap"].values,  #条形长度                height=0.6,  #高度                color=colors['hist'],  #颜色                alpha=0.5,  #透明度                edgecolor="none",  #边框线                zorder=0)  #层    ax_sw.set_xlabel("SHAP value (impact on model output)", fontsize=20)  #下x轴标题    ax_sw.set_yticks(y_pos)  #主坐标轴的Y轴刻度位置与背景条形图    ax_sw.set_yticklabels(feature_names, fontsize=12)  #Y轴刻度文本    #设置边框    ax_sw.spines['top'].set_visible(False)    ax_sw.spines['right'].set_visible(False)    ax_bar.spines['top'].set_visible(False)    ax_bar.spines['right'].set_visible(False)    apply_plot_styles(ax_sw)  #调整边框和刻度粗细    apply_plot_styles(ax_bar)  # 调整边框和刻度粗细    #保存    save_fig_dual(fig, output_folder, 'combined_shap_summary_plot')    plt.close(fig)  #关闭

# =========================================================================================# ======================================10.特征交互强度气泡热力图绘制函数=======================================# =========================================================================================def plot_bubble_heatmap(shap_interaction, feature_names, colors, output_folder):    scatter = ax.scatter(x_idx,  #x                         y_idx,  #y                         s=sizes,  #大小                         c=colors_vals,  #交互强度值                         cmap=cmap,  #颜色映射                         vmin=0,  #颜色映射的最小值                         vmax=max_val,  # 颜色映射的最大值                         alpha=1,  #透明度                         edgecolors='black',  #边缘线颜色                         linewidth=1.3,  #边缘线粗细                         zorder=2)  #层    ax.set_xticks(range(n_features))  #设置x轴刻度位置    #背景网格线    ax.grid(True,  #开启            linestyle=':',  #样式            alpha=0.6,  #透明度            zorder=1)  #层    cax = ax.inset_axes([0.05, 0.05, 0.03, 0.35])  #创建颜色条轴    cbar = plt.colorbar(scatter, cax=cax)  #生成颜色条    cbar_ticks = [0.00, max_val / 3, max_val * 2 / 3, max_val]  #颜色条刻度值    cbar.set_ticks(cbar_ticks)  #应用到颜色条上    #刻度值    cbar.ax.set_yticklabels([f'{v:.3f}' for v in cbar_ticks], fontsize=14, fontweight='bold')    y_positions = [0.4, 0.283, 0.166, 0.05]  #气泡图例位置    legend_vals = [max_val, max_val * 2 / 3, max_val / 3, 0.00]  #气泡图例显示的气泡大小    save_fig_dual(fig, output_folder, 'shap_interaction_bubble_heatmap_diag_zero')  #保存    plt.close(fig)  #关闭

# =========================================================================================# ======================================11.SHAP单因素依赖图绘制函数=======================================# =========================================================================================def plot_dependence(X_test, shap_values, feature_names, colors, save_folder):    n_features = min(9, len(feature_names))  #只选择前9个最重要特征进行绘制    #创建画布        ax2 = ax1.twinx()  #创建第二个Y轴对象        counts, bin_edges = np.histogram(x_values, bins=30)  #计算特征数据的直方图分布，返回频数和边界        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2  #计算每个区间的中心x坐标        bin_width = bin_edges[1] - bin_edges[0]  #单个区间的宽度        # 在主图层上绘制分布条形图，X坐标使用区间中心点        ax1.bar(bin_centers,  #x                counts,  #y                width=bin_width * 0.6,  #柱体宽度                align='center',  #居中对齐                color=colors['hist'],  #颜色                alpha=0.3,  #透明度                label='Distribution')  #图例标签        ax1.set_ylim(0, counts.max() * 1.1)  #左侧y轴范围        #绘制散点        ax2.scatter(x_values,  #x                    shap_vals,  #y                    alpha=0.7,  #透明度                    s=25,  #散点大小                    color=colors['shap_scatter'],  #颜色                    label='Sample',  #图例标签                    zorder=2)  #层    save_fig_dual(fig, save_folder, 'dependence_top9_grid')  #保存    plt.close(fig)  #关闭

# =========================================================================================# ======================================12.特征交互效应依赖图绘制函数=======================================# =========================================================================================def plot_interaction(X_test, shap_interaction, feature_names, colors, save_folder):    n_features = len(feature_names)  #获取特征总数    interaction_strengths = []  #用于保存任意两个特征之间的交互强度        points = ax2.scatter(x_vals, #x                             shap_vals,#y                             c=inter_vals,#根据副特征数值上色                             cmap=cmap,#颜色映射                             alpha=1,#透明度                             s=25, #大小                             zorder=2,#层                             label='Samples') #图例标签        median_val = inter_vals.median()  #计算次特征中位数，用以划分高低两个子组        #定义划分高低组信息的字典        groups = {            'low': {'mask': inter_vals <= median_val, 'color': colors['inter_low'], 'fit': colors['inter_low_fit'],'label': f'{s_name} <= {median_val:.2f}'},            'high': {'mask': inter_vals > median_val, 'color': colors['inter_high'], 'fit': colors['inter_high_fit'],'label': f'{s_name} > {median_val:.2f}'}        }        #颜色条轴        cax_auto = ax2.inset_axes([1.22, 0.0, 0.04, 1.0])        #创建颜色条        cbar = fig.colorbar(points, cax=cax_auto)        #标题        cbar.set_label(s_name, size=18, labelpad=5)        #刻度设置        cbar.ax.tick_params(labelsize=18)        cbar.ax.tick_params(axis='y', width=2, length=4, direction='in')        # 控制颜色条自身的外边框线宽        for spine in cbar.ax.spines.values():            spine.set_linewidth(1.5)        # 控制颜色条自身的外边框线宽        for spine in cbar.ax.spines.values():            spine.set_linewidth(1.5)

# =========================================================================================# ======================================13.执行部分=======================================# =========================================================================================if __name__ == '__main__':    excel_file_path = r'data.xlsx'  #原始数据    test_size = 0.3  #测试集比例    random_state = 0  #随机种子    os.makedirs(output_folder, exist_ok=True)  #是否存在，若无则自动创建一个新文件夹    df = pd.read_excel(excel_file_path)  #读取数据    y = df[target_column_name]  #目标    X = df.drop(columns=[target_column_name])  #特征    metrics = {}  #保存性能指标    #遍历数据评估    for name, X_data, y_data in [('train', X_train, y_train),                                 ('test', X_test, y_test)]:        y_pred = best_model.predict(X_data)  #预测    explainer = shap.Explainer(best_model)  #shap解析器    shap_values_obj = explainer(X_test)  # 测试集分析    #提取SHAP的基准预测值    base_values = shap_values_obj.base_values[0] if isinstance(shap_values_obj.base_values, (np.ndarray,list)) else shap_values_obj.base_values    #是否批量绘图    plot_all =True

excel_file_path = r'data.xlsx'  #原始数据

2.设置目标变量，执行部分：

target_column_name = 'FVC'  #目标变量

output_folder = r'效应图'  #结果路径

4.设置测试集的比例，执行部分：

test_size = 0.3  #测试集比例

5.设置随机种子，执行部分：

random_state = 0  #随机种子

6.设置超参数网格，执行部分：

param_grid = {'n_estimators': [10, 20, 50],'max_depth': [2, 3, 5],'learning_rate': [0.05, 0.1]}

7.设置是否要进行批量绘图，执行部分：

plot_all =True