科研绘图-python/R语言绘制基于Boruta特征选择算法的山脊图

import osimport warningsimport numpy as npimport pandas as pdimport xgboost as xgbimport matplotlib.pyplot as pltfrom scipy.stats import gaussian_kdewarnings.filterwarnings('ignore')plt.rcParams['font.sans-serif'] = ['DejaVu Sans']plt.rcParams['axes.unicode_minus'] = Falseplt.rcParams['figure.dpi'] = 300plt.rcParams['savefig.dpi'] = 300plt.rcParams['figure.figsize'] = (12, 10)

file_path = "data(山脊图).xlsx"if not os.path.exists(file_path):    clinical_dataset = simulate_clinical_data(num_patients=500, num_biomarkers=30)    clinical_dataset.to_excel(file_path, index=False)    print("已生成模拟数据并保存到文件")else:    clinical_dataset = pd.read_excel(file_path)    print(f"已从文件读取数据，数据形状: {clinical_dataset.shape}")    print(f"列名: {list(clinical_dataset.columns)}")

def auto_detect_columns(data):    target_candidates = ['outcome', 'target', 'label', 'y', 'class', 'diagnosis', 'response']    for candidate in target_candidates:        if candidate in data.columns:            target_col = candidate            feature_cols = [col for col in data.columns if col != target_col]            print(f"自动识别目标列为: '{target_col}'")            print(f"特征数量: {len(feature_cols)}")            return target_col, feature_cols    print("警告: 未找到预定义的目标列名称")    print("请从以下列中选择目标列:")    for i, col in enumerate(data.columns):        print(f"{i}: {col}")    target_col = data.columns[-1]    feature_cols = [col for col in data.columns if col != target_col]    print(f"默认使用最后一列作为目标列: '{target_col}'")    print(f"特征数量: {len(feature_cols)}")    return target_col, feature_cols

target_col, feature_cols = auto_detect_columns(clinical_dataset)

def get_xgb_feature_importance(data, target_col=target_col, feature_cols=feature_cols, n_iter=50):    X = data[feature_cols]    y = data[target_col]    print(f"使用目标列: {target_col}")    print(f"使用特征列: {feature_cols}")    imp_history = []    random_seeds = [250307 + i for i in range(n_iter)]    for seed in random_seeds:        model = xgb.XGBClassifier(objective='binary:logistic', random_state=seed, n_estimators=100, max_depth=3,                                  learning_rate=0.1, subsample=0.8, colsample_bytree=0.8, n_jobs=-1)        model.fit(X, y)        imp = model.get_booster().get_score(importance_type='gain')        imp_full = {feat: imp.get(feat, 0.0) for feat in feature_cols}        imp_history.append(imp_full)    imp_history_df = pd.DataFrame(imp_history)    for i in range(3):        shadow_col = f"shadow_{i + 1}"        imp_history_df[shadow_col] = np.random.normal(loc=imp_history_df[feature_cols].mean().mean(),                                                      scale=imp_history_df[feature_cols].std().mean(), size=n_iter)    return imp_history_df

imp_history = get_xgb_feature_importance(clinical_dataset, target_col, feature_cols, n_iter=50)

def classify_features(imp_history, feature_cols):    shadow_cols = [col for col in imp_history.columns if 'shadow' in col]    real_cols = feature_cols    shadow_95 = imp_history[shadow_cols].quantile(0.95, axis=0).mean()    feature_decisions = []    for feat in real_cols:        feat_95 = imp_history[feat].quantile(0.95)        if feat_95 > shadow_95:            decision = "Confirmed"        elif imp_history[feat].quantile(0.75) > shadow_95:            decision = "Tentative"        else:            decision = "Rejected"        feature_decisions.append(decision)    real_decision = pd.DataFrame({'feature': real_cols, 'decision': feature_decisions})    shadow_decision = pd.DataFrame({'feature': shadow_cols, 'decision': ["Shadow"] * len(shadow_cols)})    all_decision = pd.concat([real_decision, shadow_decision], ignore_index=True)    return all_decision

decision_all = classify_features(imp_history, feature_cols)print("\n特征分类统计:")print(decision_all['decision'].value_counts())

imp_history_long = imp_history.melt(var_name='feature', value_name='importance')imp_history_long = imp_history_long.merge(decision_all, on='feature', how='left')imp_history_long_clean = imp_history_long[np.isfinite(imp_history_long['importance'])].copy()feature_order = imp_history_long_clean.groupby('feature')['importance'].median().sort_values().index.tolist()

def plot_ridgeplot(data, feature_order, color_scheme, save_path):    fig, ax = plt.subplots(figsize=(12, 10))    scale = 0.8    y_offsets = {feat: i * scale for i, feat in enumerate(feature_order)}    for feat in feature_order:        feat_data = data[data['feature'] == feat]['importance']        kde = gaussian_kde(feat_data)        x_vals = np.linspace(feat_data.min(), feat_data.max(), 500)        y_vals = kde(x_vals)        y_vals = y_vals / y_vals.max() * scale * 0.7        y_base = y_offsets[feat]        decision = data[data['feature'] == feat]['decision'].iloc[0]        color = color_scheme.get(decision, (0.8, 0.8, 0.8, 0.6))        ax.fill_between(x_vals, y_base, y_base + y_vals, color=color[:3], alpha=color[3], edgecolor=None, linewidth=0)        ax.plot(x_vals, y_base + y_vals, color='black', linewidth=0.8, linestyle='-')    ax.grid(True, color='gray', alpha=0.3, linestyle='-', linewidth=0.5)    ax.set_yticks([y_offsets[feat] for feat in feature_order])    ax.set_yticklabels(feature_order, fontsize=10)    ax.set_xlabel('Importance (Gain)', fontsize=12)    ax.set_ylabel('Features', fontsize=12)    ax.set_title('Feature Importance Ridge Plot Based on XGBoost', fontsize=14, pad=20)    legend_elements = [plt.Rectangle((0, 0), 1, 1, facecolor=color[:3], alpha=color[3], edgecolor='black', label=dec)                       for dec, color in color_scheme.items()]    ax.legend(handles=legend_elements, loc='right', frameon=True, fontsize=10)    plt.tight_layout()    plt.savefig(save_path, bbox_inches='tight', dpi=300)    plt.close()

library(rstudioapi)setwd(dirname(rstudioapi::getSourceEditorContext()$path))library(readxl)library(xgboost)library(dplyr)library(tidyr)library(ggplot2)library(ggridges)library(scales)library(parallel)library(tibble)  # 加载tibble包，提供rownames_to_column函数options(warn=-1)file_path <- "data(山脊图).xlsx"if (!file.exists(file_path)) {  stop("文件不存在，请先生成数据")} else {  clinical_dataset <- read_excel(file_path)  cat(sprintf("已从文件读取数据，数据形状: %d x %d\n", nrow(clinical_dataset), ncol(clinical_dataset)))  cat(sprintf("列名: %s\n", paste(colnames(clinical_dataset), collapse=", ")))}auto_detect_columns <- function(data) {  target_candidates <- c('outcome', 'target', 'label', 'y', 'class', 'diagnosis', 'response')  for (candidate in target_candidates) {    if (candidate %in% colnames(data)) {      target_col <- candidate      feature_cols <- colnames(data)[colnames(data) != target_col]      cat(sprintf("自动识别目标列为: '%s'\n", target_col))      cat(sprintf("特征数量: %d\n", length(feature_cols)))      return(list(target_col = target_col, feature_cols = feature_cols))    }  }  cat("警告: 未找到预定义的目标列名称\n")  cat("请从以下列中选择目标列:\n")  for (i in seq_along(colnames(data))) {    cat(sprintf("%d: %s\n", i-1, colnames(data)[i]))  }  target_col <- colnames(data)[length(colnames(data))]  feature_cols <- colnames(data)[colnames(data) != target_col]  cat(sprintf("默认使用最后一列作为目标列: '%s'\n", target_col))  cat(sprintf("特征数量: %d\n", length(feature_cols)))  return(list(target_col = target_col, feature_cols = feature_cols))}detect_result <- auto_detect_columns(clinical_dataset)target_col <- detect_result$target_colfeature_cols <- detect_result$feature_colsget_xgb_feature_importance <- function(data, target_col, feature_cols, n_iter=50) {  X <- as.matrix(data[, feature_cols])  y <- as.numeric(data[, target_col][[1]])  cat(sprintf("使用目标列: %s\n", target_col))  cat(sprintf("使用特征列: %s\n", paste(feature_cols, collapse=", ")))  imp_history <- list()  random_seeds <- 250307 + 0:(n_iter-1)  for (seed in random_seeds) {    params <- list(      objective = "binary:logistic",      seed = seed,      max_depth = 3,      eta = 0.1,      subsample = 0.8,      colsample_bytree = 0.8,      nthread = detectCores()    )    dtrain <- xgb.DMatrix(data = X, label = y)    model <- xgb.train(      params = params,      data = dtrain,      nrounds = 100,      verbose = 0    )    # 核心修正：移除R中不存在的importance_type参数    imp <- xgb.importance(feature_names = feature_cols, model = model)    imp_df <- as.data.frame(imp)    imp_full <- setNames(rep(0, length(feature_cols)), feature_cols)    imp_full[imp_df$Feature] <- imp_df$Gain    imp_history[[as.character(seed)]] <- imp_full  }  imp_history_df <- do.call(rbind, imp_history) %>% as.data.frame()  for (i in 1:3) {    shadow_col <- sprintf("shadow_%d", i)    mean_val <- mean(unlist(imp_history_df[, feature_cols]))    std_val <- sd(unlist(imp_history_df[, feature_cols]))    imp_history_df[[shadow_col]] <- rnorm(n_iter, mean = mean_val, sd = std_val)  }  return(imp_history_df)}imp_history <- get_xgb_feature_importance(clinical_dataset, target_col, feature_cols, n_iter=50)classify_features <- function(imp_history, feature_cols) {  shadow_cols <- colnames(imp_history)[grepl("shadow", colnames(imp_history))]  real_cols <- feature_cols  shadow_95 <- mean(apply(imp_history[, shadow_cols], 2, quantile, 0.95))  feature_decisions <- sapply(real_cols, function(feat) {    feat_95 <- quantile(imp_history[[feat]], 0.95)    if (feat_95 > shadow_95) {      "Confirmed"    } else if (quantile(imp_history[[feat]], 0.75) > shadow_95) {      "Tentative"    } else {      "Rejected"    }  })  real_decision <- data.frame(    feature = real_cols,    decision = feature_decisions,    stringsAsFactors = FALSE  )  shadow_decision <- data.frame(    feature = shadow_cols,    decision = rep("Shadow", length(shadow_cols)),    stringsAsFactors = FALSE  )  all_decision <- rbind(real_decision, shadow_decision)  return(all_decision)}decision_all <- classify_features(imp_history, feature_cols)cat("\n特征分类统计:\n")print(table(decision_all$decision))# 修正后：加载tibble包后可正常使用rownames_to_columnimp_history_long <- imp_history %>%  rownames_to_column("seed") %>%  pivot_longer(cols = -seed, names_to = "feature", values_to = "importance") %>%  left_join(decision_all, by = "feature") %>%  filter(is.finite(importance))feature_order <- imp_history_long %>%  group_by(feature) %>%  summarise(median_importance = median(importance), .groups = "drop") %>%  arrange(median_importance) %>%  pull(feature)imp_history_long$feature <- factor(imp_history_long$feature, levels = feature_order)original_colors <- c(  "Confirmed" = rgb(0.0078, 0.7373, 0.7569, 0.6),  "Tentative" = rgb(1.0, 0.7569, 0.0275, 0.6),  "Rejected" = rgb(0.8980, 0.2235, 0.2078, 0.6),  "Shadow" = rgb(0.4588, 0.4588, 0.4588, 0.6))plot_ridgeplot <- function(data, feature_order, color_scheme, save_path) {  ggplot(data, aes(x = importance, y = feature, fill = decision)) +    geom_density_ridges(      scale = 0.8,      alpha = 0.6,      size = 0.8,      color = "black",      rel_min_height = 0.01,      from = min(imp_history_long$importance)    ) +    scale_fill_manual(values = color_scheme, drop = FALSE) +    theme_bw() +    theme(      text = element_text(family = "DejaVu Sans"),      axis.text = element_text(size = 10),      axis.title = element_text(size = 12),      plot.title = element_text(size = 14, margin = margin(b = 20), hjust = 0.5),      legend.position = "right",      legend.text = element_text(size = 10),      panel.grid = element_line(color = "gray", alpha = 0.1, linewidth = 0.3)    ) +    labs(      x = "Importance (Gain)",      y = "Features",      title = "Feature Importance Ridge Plot Based on XGBoost",      fill = "Decision"    ) +    coord_cartesian(clip = "off")  ggsave(    filename = save_path,    width = 12,    height = 10,    dpi = 300,    bg = "white",    device = "png",    limitsize = FALSE  )}plot_ridgeplot(  data = imp_history_long,  feature_order = feature_order,  color_scheme = original_colors,  save_path = "xgboost_ridgeplot_original.png")cat(sprintf("\n山脊图已保存为: xgboost_ridgeplot_original.png\n"))