【生存分析机器学习】Python04-随机生存森林模型RSF及可视化

# pip install pandas numpy matplotlib seaborn scikit-learn scikit-survival lifelinesimport osimport pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as snsfrom sklearn.model_selection import train_test_splitfrom sklearn.preprocessing import StandardScaler, LabelEncoderfrom sklearn.metrics import roc_curve, aucfrom sklearn.calibration import calibration_curvefrom sklearn.inspection import permutation_importancefrom sksurv.ensemble import RandomSurvivalForestfrom sksurv.util import Survfrom sksurv.metrics import concordance_index_censoredfrom sksurv.metrics import cumulative_dynamic_aucfrom sksurv.metrics import brier_scoreimport scipy.stats as statsimport warningswarnings.filterwarnings('ignore')# 设置中文字体和输出目录plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = False# 创建结果文件夹desktop_path = os.path.join(os.path.expanduser("~"), "Desktop")results_dir = os.path.join(desktop_path, "Results模型-rsf-comprehensive")os.makedirs(results_dir, exist_ok=True)# 1. 数据准备和预处理def load_and_preprocess_data():"""    加载和预处理数据    """    try:# 尝试读取实际数据        data = pd.read_excel(os.path.join(desktop_path, "示例数据.xlsx"), sheet_name="示例数据")print("成功读取示例数据")    except:# 创建模拟数据print("未找到示例数据，创建模拟数据...")        np.random.seed(2025)        n_samples = 500        data = pd.DataFrame({'时间': np.random.exponential(100, n_samples),'结局': np.random.choice([0, 1], n_samples, p=[0.7, 0.3]),'指标1': np.random.normal(50, 10, n_samples),'指标2': np.random.normal(25, 5, n_samples),'指标3': np.random.normal(100, 20, n_samples),'指标4': np.random.normal(10, 2, n_samples),'指标5': np.random.normal(75, 15, n_samples),'指标6': np.random.normal(30, 6, n_samples),'指标7': np.random.normal(60, 12, n_samples),'指标8': np.random.choice(['A', 'B', 'C'], n_samples),'肥胖程度': np.random.choice(['正常', '超重', '肥胖'], n_samples),'教育水平': np.random.choice(['小学', '中学', '大学'], n_samples),'血型': np.random.choice(['A', 'B', 'O', 'AB'], n_samples)        })# 确保时间变量为正数        data['时间'] = np.abs(data['时间'])# 数据预处理    data['时间'] = pd.to_numeric(data['时间'], errors='coerce')    data['结局'] = pd.to_numeric(data['结局'], errors='coerce')    data['结局'] = data['结局'].fillna(0).astype(int)# 处理分类变量 - 使用独热编码    categorical_cols = ['指标8', '肥胖程度', '教育水平', '血型']for col in categorical_cols:if col in data.columns:# 使用独热编码            dummies = pd.get_dummies(data[col], prefix=col)            data = pd.concat([data, dummies], axis=1)            data = data.drop(columns=[col])# 删除缺失值    data = data.dropna()print(f"数据维度: {data.shape}")print(f"结局分布:\n{data['结局'].value_counts()}")return data# 2. 描述性统计分析def descriptive_analysis(data, results_dir):"""生成描述性统计分析"""print("\n=== 描述性统计分析 ===")# 数值变量描述    numeric_cols = data.select_dtypes(include=[np.number]).columns    numeric_cols = [col for col in numeric_cols if col not in ['时间', '结局']]    desc_stats = data[numeric_cols].describe().T    desc_stats['missing'] = data[numeric_cols].isnull().sum()# 按结局分组的描述if len(data['结局'].unique()) > 1:        grouped_stats = data.groupby('结局')[numeric_cols].mean()# 保存结果    desc_stats.to_csv(os.path.join(results_dir, "数值变量描述性统计.csv"))return desc_stats# 3. 划分训练集和测试集def split_data(data, test_size=0.3, random_state=2025):"""划分训练集和测试集"""    train_data, test_data = train_test_split(        data, test_size=test_size, random_state=random_state, stratify=data['结局']    )print(f"训练集样本数: {len(train_data)}")print(f"测试集样本数: {len(test_data)}")return train_data, test_data# 4. 构建随机生存森林模型def build_rsf_model(train_data, duration_col='时间', event_col='结局'):"""构建随机生存森林模型"""print("\n=== 构建随机生存森林模型 ===")# 选择数值型预测变量    numeric_cols = train_data.select_dtypes(include=[np.number]).columns    predictor_cols = [col for col in numeric_cols if col not in [duration_col, event_col]]

# 5. 模型性能评估def evaluate_model(rsf, train_data, test_data, predictor_cols, X_train, y_train, results_dir):"""评估模型性能"""print("\n=== 模型性能评估 ===")# 准备测试数据    X_test = test_data[predictor_cols].values    y_test = np.array([(bool(row['结局']), row['时间'])for _, row in test_data.iterrows()],                      dtype=[('event', '?'), ('time', '<f8')])# 训练集C-index    train_pred = rsf.predict(X_train)    train_cindex = concordance_index_censored(y_train['event'], y_train['time'], train_pred)[0]# 测试集C-index    test_pred = rsf.predict(X_test)    test_cindex = concordance_index_censored(y_test['event'], y_test['time'], test_pred)[0]print(f"训练集C-index: {train_cindex:.3f}")print(f"测试集C-index: {test_cindex:.3f}")# 时间依赖性ROC曲线    time_points = [12, 24, 36]  # 1年、2年、3年    auc_values = calculate_time_dependent_auc(rsf, X_test, y_test, time_points, results_dir)# Brier分数    brier_scores = calculate_brier_score(rsf, test_data, X_test, time_points, results_dir)return {'train_cindex': train_cindex,'test_cindex': test_cindex,'auc_values': auc_values,'brier_scores': brier_scores,'test_pred': test_pred,'X_test': X_test,'y_test': y_test    }def calculate_time_dependent_auc(rsf, X_test, y_test, time_points, results_dir):"""计算时间依赖性AUC"""print("\n=== 时间依赖性AUC计算 ===")# 使用scikit-survival的cumulative_dynamic_auc函数    try:# 获取风险分数        risk_scores = rsf.predict(X_test)# 计算时间依赖性AUC        auc_values = {}for time_point in time_points:# 使用cumulative_dynamic_auc计算特定时间点的AUC# 注意：这里我们使用一个时间范围，然后取最接近我们目标时间点的值times = np.percentile(y_test['time'], np.linspace(5, 95, 20))            iauc, mean_auc = cumulative_dynamic_auc(y_test, y_test, risk_scores, times)# 找到最接近目标时间点的索引            time_idx = np.argmin(np.abs(times - time_point))            auc_value = iauc[time_idx]            auc_values[time_point] = auc_valueprint(f"{time_point}个月AUC: {auc_value:.3f}")    except Exception as e:print(f"使用cumulative_dynamic_auc计算AUC时出错: {e}")print("使用替代方法计算AUC...")# 替代方法：使用生存函数计算        auc_values = calculate_auc_alternative(rsf, X_test, y_test, time_points)return auc_valuesdef calculate_auc_alternative(rsf, X_test, y_test, time_points):"""替代方法计算AUC"""    auc_values = {}# 获取生存函数    surv_funcs = rsf.predict_survival_function(X_test)for time_point in time_points:# 计算在特定时间点的风险分数（1-生存概率）        risk_scores = []for sf in surv_funcs:# 找到最接近的时间点times = sf.x            time_idx = np.argmin(np.abs(times - time_point))            risk_scores.append(1 - sf.y[time_idx])        risk_scores = np.array(risk_scores)# 创建二元标签        y_true_binary = (y_test['time'] <= time_point) & y_test['event']if y_true_binary.sum() > 0:  # 确保有阳性样本            fpr, tpr, _ = roc_curve(y_true_binary, risk_scores)            roc_auc = auc(fpr, tpr)            auc_values[time_point] = roc_aucprint(f"{time_point}个月AUC: {roc_auc:.3f}")else:print(f"{time_point}个月: 无事件发生，无法计算AUC")            auc_values[time_point] = np.nanreturn auc_valuesdef calculate_brier_score(rsf, test_data, X_test, time_points, results_dir):"""计算Brier分数"""print("\n=== Brier分数计算 ===")    brier_scores = {}# 获取生存函数    surv_funcs = rsf.predict_survival_function(X_test)for time_point in time_points:# 找到最接近的时间点times = surv_funcs[0].x        time_idx = np.argmin(np.abs(times - time_point))# 计算在该时间点的生存概率        surv_probs = []for sf in surv_funcs:            surv_probs.append(sf.y[time_idx])        surv_probs = np.array(surv_probs)# 计算观察到的结果        observed = ((test_data['时间'] <= time_point) & (test_data['结局'] == 1)).astype(int)# 计算Brier分数        brier_score = np.mean((observed - (1 - surv_probs)) ** 2)        brier_scores[time_point] = brier_scoreprint(f"{time_point}个月Brier分数: {brier_score:.3f}")# 保存Brier分数    brier_df = pd.DataFrame(list(brier_scores.items()), columns=['Time', 'Brier_Score'])    brier_df.to_csv(os.path.join(results_dir, "Brier分数.csv"), index=False)return brier_scores# 6. 可视化分析def create_visualizations(rsf, train_data, test_data, predictor_cols, performance_metrics, results_dir):"""创建所有可视化图表"""print("\n=== 创建可视化图表 ===")# 准备测试数据    X_test = performance_metrics['X_test']    y_test = performance_metrics['y_test']

# 6.1 变量重要性图    plot_variable_importance(rsf, X_test, y_test, predictor_cols, results_dir)# 6.2 生存曲线    plot_survival_curves(rsf, X_test, results_dir)# 6.3 校准曲线    plot_calibration_curve(rsf, test_data, X_test, results_dir)