build_shap.py

import pandas as pd
import numpy as np
import lightgbm as lgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, roc_auc_score, confusion_matrix
import matplotlib.pyplot as plt
import shap
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("SHAP 分析：全局和局部解释 LightGBM 模型")
print("=" * 60)

# 1. 读取数据
print("\n[1] 读取数据...")
df_base = pd.read_csv('customer/customer_base.csv')
df_behavior = pd.read_csv('customer/customer_behavior_assets.csv')

# 排序
df_behavior = df_behavior.sort_values(['customer_id', 'stat_month']).reset_index(drop=True)

print(f"    客户基础数据：{len(df_base)} 条")
print(f"    客户行为数据：{len(df_behavior)} 条")

# 2. 构造时序特征函数
def create_time_series_features(group, base_month, lookback=6):
    """
    为每个客户在特定基线月份构造时序特征
    """
    months_list = sorted(group['stat_month'].unique())
    if base_month not in months_list:
        return None

    base_idx = months_list.index(base_month)
    if base_idx < lookback - 1:
        return None

    feature_months = months_list[base_idx - lookback + 1: base_idx + 1]
    feature_data = group[group['stat_month'].isin(feature_months)].copy()

    if len(feature_data) < lookback:
        return None

    future_months = []
    base_year, base_m = int(base_month.split('-')[0]), int(base_month.split('-')[1])
    for i in range(1, 4):
        m = base_m + i
        y = base_year
        while m > 12:
            m -= 12
            y += 1
        future_months.append(f"{y}-{m:02d}")

    future_data = group[group['stat_month'].isin(future_months)]

    features = {}

    base_data = feature_data[feature_data['stat_month'] == base_month].iloc[0]
    features['total_assets'] = base_data['total_assets']
    features['deposit_balance'] = base_data['deposit_balance']
    features['financial_balance'] = base_data['financial_balance']
    features['fund_balance'] = base_data['fund_balance']
    features['insurance_balance'] = base_data['insurance_balance']
    features['product_count'] = base_data['product_count']
    features['app_login_count'] = base_data['app_login_count']
    features['app_financial_view_time'] = base_data['app_financial_view_time']
    features['credit_card_monthly_expense'] = base_data['credit_card_monthly_expense']
    features['investment_monthly_count'] = base_data['investment_monthly_count']
    features['financial_repurchase_count'] = base_data['financial_repurchase_count']

    assets = feature_data.sort_values('stat_month')['total_assets'].values

    if len(assets) >= 2:
        features['asset_growth_rate_1m'] = (assets[-1] - assets[-2]) / (assets[-2] + 1)
    if len(assets) >= 4:
        features['asset_growth_rate_3m'] = (assets[-1] - assets[-4]) / (assets[-4] + 1)
    if len(assets) >= 6:
        features['asset_growth_rate_6m'] = (assets[-1] - assets[-6]) / (assets[-6] + 1)

    features['asset_mean_6m'] = np.mean(assets)
    features['asset_std_6m'] = np.std(assets)
    features['asset_cv_6m'] = np.std(assets) / (np.mean(assets) + 1)
    features['current_vs_avg'] = (assets[-1] - np.mean(assets)) / (np.mean(assets) + 1)

    if len(assets) >= 2:
        x = np.arange(len(assets))
        slope, intercept = np.polyfit(x, assets, 1)
        features['trend_slope'] = slope
        y_pred = slope * x + intercept
        ss_res = np.sum((assets - y_pred) ** 2)
        ss_tot = np.sum((assets - np.mean(assets)) ** 2)
        features['trend_r_squared'] = 1 - (ss_res / (ss_tot + 1))

    if len(assets) >= 6:
        recent_3_mean = np.mean(assets[-3:])
        earlier_3_mean = np.mean(assets[-6:-3])
        features['asset_momentum'] = (recent_3_mean - earlier_3_mean) / (earlier_3_mean + 1)

    if len(future_data) > 0:
        max_assets_future = future_data['total_assets'].max()
        features['target'] = 1 if max_assets_future >= 1000000 else 0
    else:
        return None

    features['base_month'] = base_month
    features['customer_id'] = group['customer_id'].iloc[0]

    return features

# 3. 构造训练数据
print("\n[2] 构造时序特征和目标变量...")

base_months = ['2024-07', '2024-08', '2024-09', '2024-10', '2024-11',
               '2024-12', '2025-01', '2025-02', '2025-03']

all_samples = []

for customer_id, group in df_behavior.groupby('customer_id'):
    for base_month in base_months:
        result = create_time_series_features(group, base_month, lookback=6)
        if result is not None:
            all_samples.append(result)

df_samples = pd.DataFrame(all_samples)
print(f"    构造样本数：{len(df_samples)}")

# 4. 合并静态特征
print("\n[3] 合并静态特征...")

static_features = ['age', 'monthly_income', 'gender', 'occupation_type',
                   'lifecycle_stage', 'marriage_status', 'city_level']

df_static = df_base[['customer_id'] + static_features].copy()

df_train = df_samples.merge(df_static, on='customer_id', how='left')
print(f"    合并后样本数：{len(df_train)}")

# 5. 数据预处理
print("\n[4] 数据预处理...")

target_col = 'target'
feature_cols = [col for col in df_train.columns if col not in ['target', 'customer_id', 'base_month']]

X = df_train[feature_cols].copy()
y = df_train[target_col].copy()

from sklearn.preprocessing import LabelEncoder
categorical_cols = ['gender', 'occupation_type', 'lifecycle_stage', 'marriage_status', 'city_level']

for col in categorical_cols:
    le = LabelEncoder()
    X[col] = le.fit_transform(X[col].astype(str))

X = X.fillna(0)

print(f"    特征数：{len(feature_cols)}")
print(f"    正样本数：{y.sum()} ({y.mean()*100:.2f}%)")
print(f"    负样本数：{len(y) - y.sum()} ({(1-y.mean())*100:.2f}%)")

# 6. 划分训练集和测试集
print("\n[5] 划分数据集...")

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)

print(f"    训练集：{len(X_train)} 条")
print(f"    测试集：{len(X_test)} 条")

# 7. 训练 LightGBM 模型
print("\n[6] 训练 LightGBM 模型...")

model = lgb.LGBMClassifier(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    random_state=42,
    class_weight='balanced',
    verbosity=-1
)

model.fit(X_train, y_train)

# 8. 模型评估
print("\n[7] 模型评估...")

y_pred = model.predict(X_test)
y_pred_proba = model.predict_proba(X_test)[:, 1]

print("\n分类报告:")
print(classification_report(y_test, y_pred, target_names=['未达 100 万', '达到 100 万']))

auc_score = roc_auc_score(y_test, y_pred_proba)
print(f"AUC 分数：{auc_score:.4f}")

cm = confusion_matrix(y_test, y_pred)
print("\n混淆矩阵:")
print(cm)

# 9. SHAP 全局解释
print("\n[8] SHAP 全局解释...")

# 创建 SHAP 解释器
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# 获取正类的 SHAP 值
if isinstance(shap_values, list):
    shap_values_positive = shap_values[1]
else:
    shap_values_positive = shap_values

print(f"    SHAP 值形状：{shap_values_positive.shape}")

# 保存 SHAP 值到文件
shap_values_df = pd.DataFrame(shap_values_positive, columns=feature_cols)
shap_values_df['customer_id'] = df_train.iloc[X_test.index]['customer_id'].values
shap_values_df['base_month'] = df_train.iloc[X_test.index]['base_month'].values
shap_values_df['true_label'] = y_test.values
shap_values_df['pred_prob'] = y_pred_proba
shap_values_df.to_csv('shap_values.csv', index=False)
print("    SHAP 值已保存到 shap_values.csv")

# SHAP 摘要图（蜂群图）
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

print("\n[9] 生成 SHAP 全局解释可视化...")

plt.figure(figsize=(14, 10))
shap.summary_plot(
    shap_values_positive, 
    X_test, 
    feature_names=feature_cols,
    plot_type='dot',
    show=False
)
plt.title('SHAP 摘要图：特征对预测的全局影响\n（SHAP 值越大，对预测影响越大）', fontsize=14, pad=20)
plt.tight_layout()
plt.savefig('shap_summary_plot.png', dpi=150, bbox_inches='tight')
plt.close()
print("    SHAP 摘要图已保存到 shap_summary_plot.png")

# SHAP 重要性图（条形图）
plt.figure(figsize=(14, 10))
shap.summary_plot(
    shap_values_positive, 
    X_test, 
    feature_names=feature_cols,
    plot_type='bar',
    show=False
)
plt.title('SHAP 特征重要性：基于 SHAP 值绝对值', fontsize=14, pad=20)
plt.tight_layout()
plt.savefig('shap_importance_plot.png', dpi=150, bbox_inches='tight')
plt.close()
print("    SHAP 重要性图已保存到 shap_importance_plot.png")

# 10. SHAP 局部解释
print("\n[10] SHAP 局部解释...")

# 选择一个正样本（预测为高价值客户）
positive_indices = np.where((y_test == 1) & (y_pred_proba > 0.9))[0]
if len(positive_indices) > 0:
    sample_idx_pos = positive_indices[0]
    sample_pos = X_test.iloc[[sample_idx_pos]]
    
    print(f"\n选择正样本（高价值客户）索引：{sample_idx_pos}")
    print(f"    真实标签：{y_test.iloc[sample_idx_pos]}")
    print(f"    预测概率：{y_pred_proba[sample_idx_pos]:.4f}")
    
    # 保存样本数据
    sample_pos_df = sample_pos.copy()
    sample_pos_df['customer_id'] = df_train.iloc[X_test.index[sample_idx_pos]]['customer_id']
    sample_pos_df['base_month'] = df_train.iloc[X_test.index[sample_idx_pos]]['base_month']
    sample_pos_df.to_csv('sample_positive.csv', index=False)
    print("    正样本数据已保存到 sample_positive.csv")
    
    # Force Plot
    plt.figure(figsize=(16, 6))
    shap.force_plot(
        explainer.expected_value[1] if isinstance(explainer.expected_value, list) else explainer.expected_value,
        shap_values_positive[sample_idx_pos],
        sample_pos,
        feature_names=feature_cols,
        matplotlib=True,
        show=False
    )
    plt.title(f'SHAP 力导向图：正样本（预测概率 = {y_pred_proba[sample_idx_pos]:.4f}）', fontsize=14)
    plt.tight_layout()
    plt.savefig('shap_force_plot_positive.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("    正样本 SHAP 力导向图已保存到 shap_force_plot_positive.png")
    
    # Waterfall Plot
    plt.figure(figsize=(14, 10))
    shap.plots.waterfall(
        shap.Explanation(
            values=shap_values_positive[sample_idx_pos],
            base_values=explainer.expected_value[1] if isinstance(explainer.expected_value, list) else explainer.expected_value,
            data=sample_pos.values[0],
            feature_names=feature_cols
        ),
        max_display=15,
        show=False
    )
    plt.title(f'SHAP 瀑布图：正样本（预测概率 = {y_pred_proba[sample_idx_pos]:.4f}）', fontsize=14, pad=20)
    plt.tight_layout()
    plt.savefig('shap_waterfall_plot_positive.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("    正样本 SHAP 瀑布图已保存到 shap_waterfall_plot_positive.png")
    
    # 打印特征贡献
    print("\n    正样本特征贡献（Top 10）:")
    shap_contrib_pos = pd.DataFrame({
        'feature': feature_cols,
        'shap_value': shap_values_positive[sample_idx_pos],
        'feature_value': sample_pos.values[0]
    })
    shap_contrib_pos['abs_shap'] = np.abs(shap_contrib_pos['shap_value'])
    shap_contrib_pos = shap_contrib_pos.sort_values('abs_shap', ascending=False)
    
    print(shap_contrib_pos[['feature', 'feature_value', 'shap_value']].head(10).to_string(index=False))

# 选择一个负样本
negative_indices = np.where((y_test == 0) & (y_pred_proba < 0.1))[0]
if len(negative_indices) > 0:
    sample_idx_neg = negative_indices[0]
    sample_neg = X_test.iloc[[sample_idx_neg]]
    
    print(f"\n选择负样本（普通客户）索引：{sample_idx_neg}")
    print(f"    真实标签：{y_test.iloc[sample_idx_neg]}")
    print(f"    预测概率：{y_pred_proba[sample_idx_neg]:.4f}")
    
    # 保存样本数据
    sample_neg_df = sample_neg.copy()
    sample_neg_df['customer_id'] = df_train.iloc[X_test.index[sample_idx_neg]]['customer_id']
    sample_neg_df['base_month'] = df_train.iloc[X_test.index[sample_idx_neg]]['base_month']
    sample_neg_df.to_csv('sample_negative.csv', index=False)
    print("    负样本数据已保存到 sample_negative.csv")
    
    # Force Plot
    plt.figure(figsize=(16, 6))
    shap.force_plot(
        explainer.expected_value[1] if isinstance(explainer.expected_value, list) else explainer.expected_value,
        shap_values_positive[sample_idx_neg],
        sample_neg,
        feature_names=feature_cols,
        matplotlib=True,
        show=False
    )
    plt.title(f'SHAP 力导向图：负样本（预测概率 = {y_pred_proba[sample_idx_neg]:.4f}）', fontsize=14)
    plt.tight_layout()
    plt.savefig('shap_force_plot_negative.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("    负样本 SHAP 力导向图已保存到 shap_force_plot_negative.png")
    
    # Waterfall Plot
    plt.figure(figsize=(14, 10))
    shap.plots.waterfall(
        shap.Explanation(
            values=shap_values_positive[sample_idx_neg],
            base_values=explainer.expected_value[1] if isinstance(explainer.expected_value, list) else explainer.expected_value,
            data=sample_neg.values[0],
            feature_names=feature_cols
        ),
        max_display=15,
        show=False
    )
    plt.title(f'SHAP 瀑布图：负样本（预测概率 = {y_pred_proba[sample_idx_neg]:.4f}）', fontsize=14, pad=20)
    plt.tight_layout()
    plt.savefig('shap_waterfall_plot_negative.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("    负样本 SHAP 瀑布图已保存到 shap_waterfall_plot_negative.png")
    
    # 打印特征贡献
    print("\n    负样本特征贡献（Top 10）:")
    shap_contrib_neg = pd.DataFrame({
        'feature': feature_cols,
        'shap_value': shap_values_positive[sample_idx_neg],
        'feature_value': sample_neg.values[0]
    })
    shap_contrib_neg['abs_shap'] = np.abs(shap_contrib_neg['shap_value'])
    shap_contrib_neg = shap_contrib_neg.sort_values('abs_shap', ascending=False)
    
    print(shap_contrib_neg[['feature', 'feature_value', 'shap_value']].head(10).to_string(index=False))

# 11. 全局特征重要性对比
print("\n[11] 全局特征重要性对比...")

# LightGBM 特征重要性
feat_importance_lgbm = pd.DataFrame({
    'feature': feature_cols,
    'importance': model.feature_importances_
}).sort_values('importance', ascending=False)

# SHAP 特征重要性（基于绝对值均值）
feat_importance_shap = pd.DataFrame({
    'feature': feature_cols,
    'importance': np.mean(np.abs(shap_values_positive), axis=0)
}).sort_values('importance', ascending=False)

# 保存对比表
comparison_df = feat_importance_lgbm[['feature', 'importance']].merge(
    feat_importance_shap[['feature', 'importance']],
    on='feature',
    suffixes=('_lgbm', '_shap')
)
comparison_df.to_csv('feature_importance_comparison.csv', index=False)
print("    特征重要性对比已保存到 feature_importance_comparison.csv")

print("\nTop 10 特征重要性对比:")
print("-" * 80)
print(f"{'排名':<6} {'特征 (LightGBM)':<30} {'LightGBM重要性':<15} {'特征 (SHAP)':<30} {'SHAP重要性':<15}")
print("-" * 80)
for i in range(min(10, len(feat_importance_lgbm))):
    lgbm_feat = feat_importance_lgbm.iloc[i]['feature']
    lgbm_imp = feat_importance_lgbm.iloc[i]['importance']
    shap_feat = feat_importance_shap.iloc[i]['feature']
    shap_imp = feat_importance_shap.iloc[i]['importance']
    print(f"{i+1:<6} {lgbm_feat:<30} {lgbm_imp:<15.1f} {shap_feat:<30} {shap_imp:<15.4f}")

print("\n" + "=" * 60)
print("SHAP 分析完成！")
print("=" * 60)
print("\n生成的文件:")
print("  1. shap_values.csv - SHAP 值数据")
print("  2. shap_summary_plot.png - SHAP 摘要图（全局解释）")
print("  3. shap_importance_plot.png - SHAP 重要性图（全局解释）")
print("  4. sample_positive.csv - 正样本数据")
print("  5. sample_negative.csv - 负样本数据")
print("  6. shap_force_plot_positive.png - 正样本力导向图（局部解释）")
print("  7. shap_force_plot_negative.png - 负样本力导向图（局部解释）")
print("  8. shap_waterfall_plot_positive.png - 正样本瀑布图（局部解释）")
print("  9. shap_waterfall_plot_negative.png - 负样本瀑布图（局部解释）")
print(" 10. feature_importance_comparison.csv - LightGBM vs SHAP 特征重要性对比")