用Python和XGBoost构建冠心病预测模型:从数据探索到生产部署的实战指南

在医疗健康领域,数据科学正以前所未有的方式重塑着我们对疾病的理解和预测能力。冠心病作为全球范围内主要的健康威胁之一,其早期识别和风险评估一直是临床实践中的关键挑战。传统的风险评估工具,如Framingham风险评分,虽然提供了基础框架,但在处理复杂、多维度的个体数据时,往往显得力不从心。这正是机器学习,尤其是像XGBoost这类强大的梯度提升框架,能够大显身手的地方。

这篇文章不是一篇泛泛而谈的技术概述,而是为数据科学家和医疗AI开发者准备的一份深度实战手册。我们将抛开那些空洞的理论,直接切入核心,手把手地带你走过构建一个高性能冠心病预测模型的完整生命周期。从如何获取和清洗一份“脏乱”的临床数据集开始,到深入特征工程的细节,再到使用XGBoost进行模型训练、调优与解释,最后探讨如何将这个模型封装成一个可供临床医生或研究人员使用的可靠工具。整个过程将伴随着可复现的Python代码和我在实际项目中踩过的“坑”与总结的经验。无论你是希望将AI技术应用于医疗研究的数据科学家,还是寻求构建可落地预测工具的开发者,这篇文章都将提供一条清晰的路径。

1. 数据获取与理解:寻找并认识你的“战场”

任何机器学习项目的成败,其基石都在于数据。对于冠心病预测,理想的数据集应包含丰富的临床指标、人口统计学信息以及明确的诊断标签。幸运的是,像UCI机器学习仓库这样的平台提供了多个公开可用的心脏病数据集,例如著名的克利夫兰数据集(Cleveland Heart Disease Dataset)。这是我们本次实战的起点。

第一步永远是理解你的数据。盲目地将数据扔进模型只会得到不可靠甚至有害的结果。我们首先使用Pandas加载并探索数据。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 假设我们使用克利夫兰数据集(需从UCI下载)
# 注意:实际数据可能没有列名,需要根据数据文档手动指定
column_names = [
    'age', 'sex', 'cp', 'trestbps', 'chol', 'fbs', 'restecg',
    'thalach', 'exang', 'oldpeak', 'slope', 'ca', 'thal', 'target'
]

# 加载数据,处理缺失值(原始数据中缺失值可能用‘?’表示)
df = pd.read_csv('processed.cleveland.data', names=column_names, na_values='?')

print(f"数据集形状: {df.shape}")
print("\n前5行数据预览:")
print(df.head())
print("\n数据基本信息:")
print(df.info())
print("\n缺失值统计:")
print(df.isnull().sum())

运行上述代码后,你会立刻对数据有一个宏观认识:样本数量、特征维度、缺失情况。接下来,我们需要深入理解每个特征的含义,这对于后续的特征工程至关重要。

特征名 描述 类型/取值范围
age 年龄 数值型
sex 性别 (1 = 男; 0 = 女) 分类(二进制)
cp 胸痛类型 (1: 典型心绞痛, 2: 非典型心绞痛, 3: 非心绞痛性疼痛, 4: 无症状) 分类(有序)
trestbps 静息血压 (mm Hg) 数值型
chol 血清胆固醇 (mg/dl) 数值型
fbs 空腹血糖 > 120 mg/dl (1 = 是; 0 = 否) 分类(二进制)
restecg 静息心电图结果 (0: 正常, 1: ST-T波异常, 2: 左心室肥厚) 分类(有序)
thalach 达到的最大心率 数值型
exang 运动诱发的心绞痛 (1 = 是; 0 = 否) 分类(二进制)
oldpeak 运动相对于休息诱发的ST段压低 数值型
slope 运动高峰ST段的斜率 (1: 上坡, 2: 平坦, 3: 下坡) 分类(有序)
ca 荧光透视着色的主要血管数 (0-3) 数值型(离散)
thal 地中海贫血 (3: 正常; 6: 固定缺陷; 7: 可逆缺陷) 分类
target 诊断结果 (0: 无冠心病; 1,2,3,4: 存在冠心病,严重程度递增) 多分类标签

注意:原始数据集中,target字段有5个值(0-4)。为了简化成一个二分类问题(有无冠心病),通常将值>0的合并为1(患病)。但需注意,这可能会丢失疾病严重程度的信息。另一种做法是将其视为多分类或有序回归问题,这取决于你的具体目标。

2. 数据预处理与特征工程:从原始数据到模型“食材”

数据预处理是模型成功的“隐形守护者”。这一步的目标是将原始数据转化为干净、规整、对模型友好的格式。我们分步进行。

2.1 处理缺失值与异常值

cathal两个特征中可能存在缺失值(用‘?’表示)。我们需要决定是删除这些样本还是填充。

# 检查缺失值
print(df.isnull().sum())

# 策略1:删除缺失值较少的行(如果缺失比例很低)
# df_clean = df.dropna()

# 策略2:根据特征含义进行填充。对于‘ca’(血管数),中位数可能是一个合理选择。
# 对于‘thal’(地中海贫血),使用众数填充,因为它是分类变量。
df['ca'].fillna(df['ca'].median(), inplace=True)
df['thal'].fillna(df['thal'].mode()[0], inplace=True) # 取第一个众数

# 检查异常值:例如,胆固醇水平为0显然不合理
print(df[df['chol'] == 0]) # 查看是否有此类异常
# 假设我们发现chol有0值,将其视为缺失或根据领域知识修正
df['chol'].replace(0, np.nan, inplace=True)
df['chol'].fillna(df['chol'].median(), inplace=True)

2.2 目标变量转换与特征类型处理

如前所述,我们将多分类目标转换为二分类。同时,需要正确处理分类特征。XGBoost本身可以处理数值输入,但将分类特征进行适当的编码(如独热编码或标签编码)有时能带来更好的效果,尤其是对于无序分类特征。

# 目标变量二值化:将>0的值视为患病(1)
df['target_binary'] = (df['target'] > 0).astype(int)

# 分离特征和目标
X = df.drop(['target', 'target_binary'], axis=1)
y = df['target_binary']

# 识别分类特征和数值特征
categorical_features = ['sex', 'cp', 'fbs', 'restecg', 'exang', 'slope', 'thal']
numerical_features = ['age', 'trestbps', 'chol', 'thalach', 'oldpeak', 'ca']

# 对分类特征进行独热编码(One-Hot Encoding)
# 注意:对于有序分类特征(如cp, restecg, slope),标签编码或保留序数关系可能更合适。
# 这里为简化,对‘cp’, ‘thal’使用独热编码,对‘slope’考虑其有序性。
from sklearn.preprocessing import OneHotEncoder, LabelEncoder
from sklearn.compose import ColumnTransformer

# 使用ColumnTransformer构建预处理管道
preprocessor = ColumnTransformer(
    transformers=[
        ('num', 'passthrough', numerical_features), # 数值特征保持不变
        ('cat', OneHotEncoder(drop='first', sparse_output=False), ['sex', 'cp', 'fbs', 'restecg', 'exang', 'thal']) # 分类特征独热编码,丢弃第一类以避免共线性
    ])

X_processed = preprocessor.fit_transform(X)
# 获取编码后的特征名称(便于后续分析)
cat_encoder = preprocessor.named_transformers_['cat']
encoded_cat_features = cat_encoder.get_feature_names_out(['sex', 'cp', 'fbs', 'restecg', 'exang', 'thal'])
all_feature_names = np.concatenate([numerical_features, encoded_cat_features])

2.3 探索性数据分析与特征衍生

在投入模型之前,可视化数据分布和特征与目标的关系能提供宝贵洞见。

# 1. 目标变量分布
sns.countplot(x=y)
plt.title('目标变量分布(0: 健康, 1: 患病)')
plt.show()

# 2. 数值特征分布与目标的关系
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.ravel()
for idx, col in enumerate(numerical_features):
    sns.histplot(data=df, x=col, hue='target_binary', kde=True, ax=axes[idx])
    axes[idx].set_title(f'{col} 分布')
plt.tight_layout()
plt.show()

# 3. 特征相关性热图(处理后的数值部分)
X_processed_df = pd.DataFrame(X_processed, columns=all_feature_names)
# 计算相关性时,只选取数值型部分(包括编码后的0/1变量)
corr_matrix = X_processed_df.corr()
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, annot=False, cmap='coolwarm', center=0)
plt.title('特征相关性热图')
plt.show()

通过EDA,你可能会发现一些有趣的现象,比如thalach(最大心率)在患病组中平均值可能较低,oldpeak(ST段压低)在患病组中更高。这些观察可以指导后续的特征选择,甚至启发你创建新的特征,例如创建风险比值特征agechol的比值,或者trestbpsthalach的比值,这些复合特征有时能比单一特征提供更强的预测信号。

3. 模型构建:让XGBoost发挥威力

数据准备就绪后,我们进入核心环节——构建XGBoost分类器。XGBoost因其卓越的性能、速度和防止过拟合的能力而广受欢迎。

3.1 数据分割与基线模型

首先,我们需要将数据划分为训练集、验证集和测试集。验证集用于在训练过程中调整超参数,测试集用于最终评估模型的泛化能力。

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, confusion_matrix, classification_report

# 划分训练+验证集 和 测试集 (80% / 20%)
X_train_val, X_test, y_train_val, y_test = train_test_split(
    X_processed, y, test_size=0.2, random_state=42, stratify=y
)

# 进一步将训练+验证集划分为训练集和验证集 (75% / 25% of 80% = 60% / 20% of total)
X_train, X_val, y_train, y_val = train_test_split(
    X_train_val, y_train_val, test_size=0.25, random_state=42, stratify=y_train_val
)

print(f"训练集大小: {X_train.shape}")
print(f"验证集大小: {X_val.shape}")
print(f"测试集大小: {X_test.shape}")

现在,我们建立一个简单的XGBoost基线模型,使用默认参数。

import xgboost as xgb

# 创建DMatrix是XGBoost的高效数据格式,但不是必须的
# dtrain = xgb.DMatrix(X_train, label=y_train)
# dval = xgb.DMatrix(X_val, label=y_val)

# 使用sklearn API,更符合常规工作流
baseline_model = xgb.XGBClassifier(
    objective='binary:logistic', # 二分类逻辑回归
    eval_metric='logloss', # 评估指标
    use_label_encoder=False,
    random_state=42
)

baseline_model.fit(X_train, y_train)

# 在验证集上评估
y_val_pred = baseline_model.predict(X_val)
y_val_pred_proba = baseline_model.predict_proba(X_val)[:, 1]

print("=== 基线模型在验证集上的表现 ===")
print(f"准确率: {accuracy_score(y_val, y_val_pred):.4f}")
print(f"精确率: {precision_score(y_val, y_val_pred):.4f}")
print(f"召回率: {recall_score(y_val, y_val_pred):.4f}")
print(f"F1分数: {f1_score(y_val, y_val_pred):.4f}")
print(f"AUC分数: {roc_auc_score(y_val, y_val_pred_proba):.4f}")
print("\n分类报告:")
print(classification_report(y_val, y_val_pred))
print("\n混淆矩阵:")
print(confusion_matrix(y_val, y_val_pred))

基线模型的结果给了我们一个起点。但通常,默认参数远非最优。接下来我们需要进行超参数调优。

3.2 超参数调优:寻找最佳配置

XGBoost有许多超参数可以调整。手动网格搜索(GridSearchCV)或随机搜索(RandomizedSearchCV)是常用方法。这里我们使用RandomizedSearchCV,因为它能在更短的时间内探索更大的参数空间。

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import uniform, randint

# 定义参数分布
param_distributions = {
    'n_estimators': randint(100, 500),
    'max_depth': randint(3, 10),
    'learning_rate': uniform(0.01, 0.3), # 从0.01到0.31的均匀分布
    'subsample': uniform(0.6, 0.4), # 从0.6到1.0
    'colsample_bytree': uniform(0.6, 0.4),
    'gamma': uniform(0, 0.5), # 最小损失减少阈值
    'reg_alpha': uniform(0, 1.0), # L1正则化
    'reg_lambda': uniform(1, 2.0), # L2正则化
    'min_child_weight': randint(1, 10)
}

# 创建模型
xgb_clf = xgb.XGBClassifier(objective='binary:logistic', use_label_encoder=False, random_state=42)

# 随机搜索
random_search = RandomizedSearchCV(
    estimator=xgb_clf,
    param_distributions=param_distributions,
    n_iter=50, # 迭代次数
    scoring='roc_auc', # 以AUC作为优化目标
    cv=5, # 5折交叉验证
    verbose=1,
    random_state=42,
    n_jobs=-1 # 使用所有CPU核心
)

random_search.fit(X_train_val, y_train_val) # 使用完整的训练+验证集进行交叉验证搜索

print("最佳参数:", random_search.best_params_)
print("最佳交叉验证AUC分数: {:.4f}".format(random_search.best_score_))

# 获取最佳模型
best_xgb_model = random_search.best_estimator_

提示:对于医疗预测模型,我们通常更关注召回率(尽可能找出所有患者)和AUC(模型整体排序能力),而不仅仅是准确率。在RandomizedSearchCV中,将scoring设置为‘recall’‘roc_auc’可以引导搜索向这些目标优化。

3.3 处理类别不平衡

医疗数据中,患病与健康的样本数量常常不平衡。XGBoost提供了scale_pos_weight参数来调整正负样本的权重。一个简单的设置方法是scale_pos_weight = (负样本数) / (正样本数)

# 计算类别权重
neg_count = np.sum(y_train == 0)
pos_count = np.sum(y_train == 1)
scale_pos_weight = neg_count / pos_count
print(f"负/正样本比例: {scale_pos_weight:.2f}")

# 将scale_pos_weight加入参数调优范围,或直接用于最终模型
balanced_model = xgb.XGBClassifier(
    **random_search.best_params_, # 使用之前找到的最佳参数
    scale_pos_weight=scale_pos_weight,
    objective='binary:logistic',
    use_label_encoder=False,
    random_state=42
)
balanced_model.fit(X_train, y_train)

4. 模型评估、解释与部署

一个模型不仅在验证集上表现好,更需要我们理解它为何做出决策,并且能够稳定地运行在生产环境中。

4.1 在测试集上进行最终评估

使用从未参与过训练和调优的测试集,对最佳模型进行最终、公正的评估。

# 使用最佳模型(考虑类别平衡的)在测试集上预测
y_test_pred = balanced_model.predict(X_test)
y_test_pred_proba = balanced_model.predict_proba(X_test)[:, 1]

print("=== 最佳模型在测试集上的最终表现 ===")
print(f"准确率: {accuracy_score(y_test, y_test_pred):.4f}")
print(f"精确率: {precision_score(y_test, y_test_pred):.4f}")
print(f"召回率: {recall_score(y_test, y_test_pred):.4f}")
print(f"F1分数: {f1_score(y_test, y_test_pred):.4f}")
print(f"AUC分数: {roc_auc_score(y_test, y_test_pred_proba):.4f}")

# 绘制ROC曲线
from sklearn.metrics import roc_curve, auc
fpr, tpr, thresholds = roc_curve(y_test, y_test_pred_proba)
roc_auc = auc(fpr, tpr)

plt.figure()
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (area = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve')
plt.legend(loc="lower right")
plt.show()

4.2 模型可解释性:SHAP值分析

对于医疗应用,模型的“黑箱”特性是不可接受的。医生需要知道模型是基于哪些因素做出预测的。SHAP (SHapley Additive exPlanations) 是目前最流行的模型解释工具之一。

import shap

# 创建一个SHAP解释器
explainer = shap.TreeExplainer(balanced_model)
# 计算测试集样本的SHAP值
shap_values = explainer.shap_values(X_test)

# 1. 特征重要性摘要图
shap.summary_plot(shap_values, X_test, feature_names=all_feature_names, plot_type="bar")

# 2. 详细SHAP摘要图(显示特征值与SHAP值的关系)
shap.summary_plot(shap_values, X_test, feature_names=all_feature_names)

# 3. 对单个预测进行解释(例如,测试集中第一个样本)
sample_idx = 0
shap.force_plot(explainer.expected_value, shap_values[sample_idx,:], X_test[sample_idx,:], feature_names=all_feature_names)

SHAP图能清晰地告诉你,对于整个模型(摘要图)或单个预测(力图),每个特征是如何影响最终结果的。例如,你可能会发现thal(地中海贫血类型)、ca(血管数)和cp(胸痛类型)是驱动预测的最重要特征。这不仅能验证模型的医学合理性,也能为临床医生提供决策支持。

4.3 模型部署:从Jupyter Notebook到生产API

模型训练和评估完成后,我们需要将其封装起来,以便集成到更大的应用系统中。这里介绍两种简单实用的方法:使用pickle序列化模型,以及构建一个简单的Flask API。

第一步:保存模型和预处理管道

import joblib # 比pickle更适合存储大型numpy数组

# 保存预处理管道和模型
model_assets = {
    'preprocessor': preprocessor,
    'model': balanced_model,
    'feature_names': all_feature_names.tolist()
}

joblib.dump(model_assets, 'coronary_heart_disease_predictor.pkl')
print("模型资产已保存至 'coronary_heart_disease_predictor.pkl'")

第二步:创建预测函数和Flask API

创建一个名为predict_service.py的文件:

# predict_service.py
import joblib
import numpy as np
import pandas as pd
from flask import Flask, request, jsonify

app = Flask(__name__)

# 加载模型资产
model_assets = joblib.load('coronary_heart_disease_predictor.pkl')
preprocessor = model_assets['preprocessor']
model = model_assets['model']
feature_names = model_assets['feature_names']

def preprocess_input(input_dict):
    """将传入的字典转换为模型可处理的DataFrame格式"""
    # 确保顺序与训练时一致
    df_input = pd.DataFrame([input_dict])
    # 注意:这里假设输入字典的键与原始特征名一致
    # 实际应用中需要更健壮的校验和错误处理
    return df_input

@app.route('/predict', methods=['POST'])
def predict():
    try:
        # 获取JSON数据
        data = request.get_json(force=True)
        # 预处理
        df_input = preprocess_input(data)
        # 应用相同的预处理转换
        X_input_processed = preprocessor.transform(df_input)
        # 预测
        prediction = model.predict(X_input_processed)[0]
        prediction_proba = model.predict_proba(X_input_processed)[0][1] # 患病概率

        # 构建返回结果
        result = {
            'prediction': int(prediction),
            'probability': float(prediction_proba),
            'interpretation': '高风险' if prediction == 1 else '低风险'
        }
        return jsonify(result), 200

    except Exception as e:
        return jsonify({'error': str(e)}), 400

if __name__ == '__main__':
    # 在生产环境中,应使用WSGI服务器如Gunicorn
    app.run(host='0.0.0.0', port=5000, debug=False)

现在,你可以通过发送一个HTTP POST请求到http://localhost:5000/predict来使用模型。请求体应是一个JSON对象,包含所有必要的特征。

# 示例使用curl进行测试
curl -X POST http://localhost:5000/predict \
  -H "Content-Type: application/json" \
  -d '{"age": 55, "sex": 1, "cp": 3, "trestbps": 130, "chol": 250, "fbs": 0, "restecg": 0, "thalach": 150, "exang": 0, "oldpeak": 1.2, "slope": 2, "ca": 0, "thal": 3}'

4.4 持续集成与监控

模型部署上线并非终点。你需要建立监控机制来跟踪模型在生产环境中的表现,防止因数据分布变化(数据漂移)而导致性能下降。可以定期(如每月)用新收集的带标签数据评估模型性能,并设置性能下降的警报阈值。

此外,考虑将整个流程(数据预处理、训练、评估、部署)脚本化,并使用如GitLab CI/CD、Jenkins或Airflow等工具进行自动化。这确保了模型的可复现性和可维护性。

构建一个真正可靠、可解释且可部署的冠心病预测模型,远不止是调出一个高AUC分数。它涉及对数据的深刻理解、严谨的工程实践以及对应用场景的持续关注。希望这份指南能为你提供一个坚实的起点。在实际项目中,你还需要与领域专家(心血管医生)紧密合作,确保每一个特征的处理、每一个预测结果都经得起医学逻辑的推敲。毕竟,我们构建的不仅是一个模型,更是一个可能影响他人健康决策的工具。

Logo

Agent 垂直技术社区,欢迎活跃、内容共建。

更多推荐