Python绘制GAM模型非线性交互效应等高线热力图(带有峰值区间和谷值区间)

    # 创建画布    fig, axes = plt.subplots(num_rows,                             num_cols,                             figsize=(fig_width, fig_height))    plt.subplots_adjust(hspace=0.25, wspace=0.25)  #调整子图之间的水平和垂直间距    if num_rows == 1 and num_cols == 1: # 如果只有1行1列        axes = np.array([[axes]]) # 将axes包装成二维数组以便统一索引    elif num_rows == 1: # 如果只有1行多列        axes = axes.reshape(1, -1) # 将axes重塑为1行N列的二维数组    elif num_cols == 1: # 如果只有多行1列        axes = axes.reshape(-1, 1) # 将axes重塑为N行1列的二维数组

for i, row_val in enumerate(row_unique_vals): # 遍历每一个行分组值        for j, col_val in enumerate(col_unique_vals): # 遍历每一个列分组值            ax = axes[i, j] # 获取当前行列对应的子图对象            base_sub_df = df[                (df[GROUP_ROW_COL_NAME] == row_val) &  # 匹配当前行分组值                (df[GROUP_COL_COL_NAME] == col_val)  # 匹配当前列分组值                ]  # 根据当前行列分组条件提取基础数据子集            #计算X轴特征的阈值            x_threshold_low = np.percentile(base_sub_df[ANALYSIS_X_COL_NAME], 5)  #计算第5百分位数            x_threshold_high = np.percentile(base_sub_df[ANALYSIS_X_COL_NAME], 95)  #计算第95百分位数            #计算Y轴特征的阈值            y_threshold_low = np.percentile(base_sub_df[ANALYSIS_Y_COL_NAME], 5)  #计算第5百分位数            y_threshold_high = np.percentile(base_sub_df[ANALYSIS_Y_COL_NAME], 95)  #计算第95百分位数            #样本量检查            if len(sub_df) < 10:  #检查筛选后的样本量是否过少                ax.axis('off')  #如果样本不足，关闭当前子图的坐标轴显示                continue  # 跳过

            X_matrix = sub_df[feature_cols].values # 提取特征矩阵数据            Y_target = sub_df[TARGET_COL_NAME].values # 提取目标变量数据            # 坐标轴范围            bounds = {                'x': (sub_df[ANALYSIS_X_COL_NAME].min(), sub_df[ANALYSIS_X_COL_NAME].max()),                'y': (sub_df[ANALYSIS_Y_COL_NAME].min(), sub_df[ANALYSIS_Y_COL_NAME].max())            }            # 划分训练集和测试集            X_train, X_test, y_train, y_test = train_test_split(X_matrix, Y_target, test_size=0.2, random_state=42)            # GAM模型拟合            gam_terms = te(analysis_x_idx, analysis_y_idx, n_splines=10)            for k in range(len(feature_cols)): # 遍历所有特征列索引                if k != analysis_x_idx and k != analysis_y_idx: # 如果不是X轴或Y轴对应的特征                    gam_terms += s(k) # 为其他特征添加平滑样条项            # 构建并拟合GAM模型            gam = LinearGAM(gam_terms).fit(X_train, y_train)            y_pred_test = gam.predict(X_test) # 在测试集上进行预测            r2 = r2_score(y_test, y_pred_test) # 计算R2            rmse = np.sqrt(mean_squared_error(y_test, y_pred_test)) # RMSE            p_values = gam.statistics_['p_values'] # 获取模型的P值统计信息            interaction_p_value = p_values[0] if len(p_values) > 0 else 1.0 #获取交互项的P值，如果没有则设为1.0

            x_min, x_max = bounds['x'] #X轴边界            y_min, y_max = bounds['y'] #Y轴边界            xi = np.linspace(x_min, x_max, 100) # 在X轴范围内生成100个网格点            yi = np.linspace(y_min, y_max, 100) # 在Y轴范围内生成100个网格点            Xi, Yi = np.meshgrid(xi, yi) # 生成二维网格坐标矩阵            Xi_flat = Xi.ravel() # 将X网格矩阵展平为一维数组            Yi_flat = Yi.ravel() # 将Y网格矩阵展平为一维数组            n_grid = len(Xi_flat) # 获取网格点的总数量            grid_matrix[:, analysis_x_idx] = Xi_flat # 将分析的X特征列替换为网格值            grid_matrix[:, analysis_y_idx] = Yi_flat # 将分析的Y特征列替换为网格值            Zi = gam.predict(grid_matrix).reshape(Xi.shape) # 对网格数据进行预测并重塑为网格形状

            z_flat = Zi.flatten() # 将预测结果Z展平            p05 = np.percentile(z_flat, 5) # 计算第5百分位数            p10 = np.percentile(z_flat, 10) # 计算第10百分位数            p90 = np.percentile(z_flat, 90) # 计算第90百分位数            p95 = np.percentile(z_flat, 95) # 计算第95百分位数            current_levels = sorted(list(set(current_levels))) # 去重并排序分级阈值            if len(current_levels) < 6: # 如果分级数量不足6个（由于数值重复等原因）                current_levels = np.linspace(z_flat.min(), z_flat.max(), 6) # 则线性生成6个分级            is_bottom_row = (i == num_rows - 1) # 判断是否为最后一行            is_left_col = (j == 0) # 判断是否为第一列            title_str = f"{row_val}-{col_val}" # 生成子图标题            # 调用绘图函数绘制子图            draw_single_plot_content(ax, Xi, Yi, Zi, r2, rmse, title_str, bounds, levels=current_levels)            if is_bottom_row: # 如果是最后一行                ax.set_xlabel(ANALYSIS_X_COL_NAME, fontsize=18, fontweight='bold') # 设置X轴标签            if is_left_col: # 如果是第一列                ax.set_ylabel("BH (m)", fontsize=18, fontweight='bold') # 设置Y轴标签

   # 在图形上添加一个自定义位置的坐标轴用于放置图例    legend_ax = fig.add_axes([0.12,  #最左侧                              0.90,  #下边界                              0.75,  #该坐标轴的宽度占据画布总宽度                              0.03])  #该坐标轴的高度占据画布总高度    legend_ax.axis('off') # 关闭图例坐标轴的显示    # 定义图例信息列表    legend_info = [        (CURRENT_COLORS[4], "95~100%"),        (CURRENT_COLORS[3], "90~95%"),        (CURRENT_COLORS[0], "0~5%"),        (CURRENT_COLORS[1], "5~10%")    ]    start_x = 0.015# 图例起始X坐标    for idx, (color, label) in enumerate(legend_info): # 遍历图例信息        # 添加图例文本        legend_ax.text(start_x + idx * 0.24 + 0.12,  #X轴坐标                       0.02,  #Y坐标                       label,  #标签文本                       transform=legend_ax.transAxes,  #坐标系                       fontsize=24,  #字体大小                       va='bottom',  #垂直对齐方式                       fontweight='bold')  #加粗

SELECTED_SCHEME = 20

2.设置原始数据文件路径：

df = pd.read_excel(r"Data.xlsx")  #读取数据

SAVE_DIR = r"峰值区间"  #结果保存路径

4.定义目标变量：

TARGET_COL_NAME = 'Target_CE'  # 目标变量

GROUP_ROW_COL_NAME = 'Income_Group'  # 行分组GROUP_COL_COL_NAME = 'Climate_Zone'  # 列分组

6.定义要分析的特征，就是x、y轴：

ANALYSIS_X_COL_NAME = 'BD'  # X轴对应的特征名称ANALYSIS_Y_COL_NAME = 'BH'  # Y轴对应的特征名称

plt.savefig(fr"Analysis_Result_Scheme{SELECTED_SCHEME}.png",dpi=300, bbox_inches='tight')plt.savefig(fr"Analysis_Result_Scheme{SELECTED_SCHEME}.pdf",bbox_inches='tight')