Python机器学习模型分析

模型训练与预测

一般模型训练方法

利用SkLearn库进行Logistics回归模型预测

def model_train_and_predict(y_var):
    """
    模型训练和预测函数，x_var和y_var是用于测试的数据集
    1, 需要指定机器学习模型
    2, 利用fit方法训练模型（采用训练组）
    3, 利用predict方法进行预测（采用测试组）->y_pred是预测结果
    4, 利用predict_proba方法进行预测概率->y_pred_proba是预测概率
    """
    model = LogisticRegression(random_state=0, max_iter=1000)
    model.fit(x_train, y_train)
    y_pred = model.predict(x_var)
    y_pred_proba = model.predict_proba(x_var)[:, 1]
    # 用于拟合真实结果得到TP和FP概率
    fpr, tpr, _ = metrics.roc_curve(y, y_pred)

多模型训练及预测

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# 利用sklearn库的机器学习模型
from sklearn import metrics
# 导入训练集和测试集模块
from sklearn.model_selection import train_test_split
# 导入机器学习模型
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.neighbors import KNeighborsClassifier
import xgboost, lightgbm, catboost
from sklearn.calibration import calibration_curve
# 这个用来屏蔽掉Logistics回归分析的警报信息
import warnings
warnings.filterwarnings("ignore")

# 读取数据(Excel文件)
df = pd.read_excel("C:/Users/Admin/Desktop/原始数据.xls")

# 提取特征和目标变量,直接使用DataFrame进行行索引，列索引
# features需要是列表变量
features = ["Var1", "Var2", "Var3", "Var4", "Var5"]
X = df[features].values
y = df['Result'].values

# 划分训练集和测试集
# test_size:测试集所占比例
# random_state:随机数种子(随机数种子相同，每次生成的随机数相同)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)

# 创建模型列表,默认字段 模型ID: 模型方法
models_list = {
    "Logistic": LogisticRegression(random_state=0, max_iter=1000),
    "SVM": SVC(probability=True,random_state=0),
    "GBM": GradientBoostingClassifier(random_state=0),
    "NeuralNetwork": MLPClassifier(hidden_layer_sizes=(8,), max_iter=500, random_state=0),
    "RandomForest": RandomForestClassifier(random_state=0),
    "Xgboost": xgboost.XGBClassifier(random_state=0),
    "KNN": KNeighborsClassifier(),
    "Adaboost": AdaBoostClassifier(random_state=0),
    "LightGBM": lightgbm.LGBMClassifier(n_estimators=500, random_state=0),
    "CatBoost": catboost.CatBoostClassifier(verbose=0, random_state=0) 
}
# 用于承接模型预测结果
result = {}
# 布置10*10的画布
plt.figure(figsize=(10, 10)).clf()

# 设置模型评估的使用数据
x_var = X_test
y_var = y_test

# 利用遍历模型字典进行模型训练和预测
for name, model in models_list.items():
    model.fit(X_train, y_train)
    # 这里利用测试集得到预测结概率估计
    y_proba = model.predict_proba(X_test)[:,1]
    # 计算校准曲线的真正概率和预测概率
    porb_true, porb_pred = calibration_curve(y_test, y_proba, n_bins=4, strategy='quantile')
    # 获取tpr和fpr用于绘制ROC曲线
    fpr, tpr, _ = metrics.roc_curve(y, y_proba)
    # 求得ROC曲线下面积AUC值
    auc = round(metrics.roc_auc_score(y, y_proba), 3)

    # 利用Booststrap方法计算AUC的95% CI
    bootstrapped_aucs = []
    rng = np.random.RandomState(0)
    for _ in range(1000):
        # 有放回抽样（保持测试集的标签和预测概率对应）
        indices = rng.randint(0, len(y), len(y))
        if len(np.unique(y[indices])) < 2:
            # 避免抽样后只有一类标签（无法计算AUC）
            continue
        # 计算当前抽样的AUC
        bootstrap_auc = metrics.roc_auc_score(
            y[indices], y_proba[indices]
        )
        bootstrapped_aucs.append(bootstrap_auc)
    
    # 计算95%置信区间（2.5%和97.5%分位数）
    sorted_aucs = np.array(bootstrapped_aucs)
    sorted_aucs.sort()
    auc_ci_lower = round(sorted_aucs[int(0.025 * len(sorted_aucs))], 3)
    auc_ci_upper = round(sorted_aucs[int(0.975 * len(sorted_aucs))], 3)
    
    str_auc = '(AUC=' + "%.3f"%auc + ', 95%CI:' + "%.3f"%auc_ci_lower + '-' + "%.3f"%auc_ci_upper + ')'
    plt.plot(fpr, tpr, label=name + str_auc)

# 绘制k=1的虚线作为辅助线
plt.plot([0,1], [0,1], "--")
plt.legend(loc="lower right", fontsize=12)
plt.grid(alpha=0.3)
plt.show()

预测评估与指标

		真实	结果
		阳性	阴性
预测	阳性	真阳性（True Positive）	假阳性（False Positive）
结果	阴性	假阴性（False Negative）	真阴性（True Negative）

• TPR（True Positive Rate，真阳性率）= TP/(TP+FN），所有阳性结果中被正确预测成阳性的比例
• FNR（False Negative Rate，假阴性率）=FN/(TP+FN)，和真阳性合为1
• FPR（False Positive Rate，假阳性率）=FP/(FP+TN)，所有阴性结果中被错误预测为阳性的比例
• TNR（True Negative Rate，真阴性率）=TN/(FP+TN)，和假阳性合为1
• 灵敏度（Sensitivity，敏感性），即TPR，在阳性患者可正确预测出阳性病例的概率
• 特异性（Specificity，特异性），即TNR，在阴性患者中可正确预测出阴性病例的概率
• 精确度（Precision），TP/(TP+FP)，预测阳性病例结果中预测正确的概率（即可以预测正确）
• 召回率（Recall），即灵敏度/TPR（即可以预测完全）
• F-Score（F1-Measure，F值），2*(精确度*召回率)/(精确度+召回率)
• 准确率（Accuracy），(TP+TN)/n，所有预测结果中正确的概率

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

models_list = {
    "LR": LogisticRegression(random_state=0, max_iter=1000),
    "SVM": SVC(probability=True,random_state=0),
    "GBM": GradientBoostingClassifier(random_state=0),
    "NN": MLPClassifier(hidden_layer_sizes=(64, 32), max_iter=500, random_state=0),
    "RF": RandomForestClassifier(random_state=0),
    "XGB": xgboost.XGBClassifier(random_state=0),
    "KNN": KNeighborsClassifier(),
    "ADB": AdaBoostClassifier(random_state=0),
    "LGBM": lightgbm.LGBMClassifier(n_estimators=500, random_state=0),
    "CB": catboost.CatBoostClassifier(verbose=0, random_state=0) 
}

x_var = X_train
y_var = y_train

for name, model in models_list.items():
    v = {}
    model.fit(X_train, y_train)
    y_proba = model.predict_proba(x_var)[:,1]
    y_pred = model.predict(x_var)

    # 计算指标（与表4完全对应）
    v["AUC"] = roc_auc_score(y_var, y_proba)
    v["Brier 分数"] = round(brier_score_loss(y_var, y_proba), 3)
    v["准确率"] = round(accuracy_score(y_var, y_pred), 3)
    v["敏感性"] = round(recall_score(y_var, y_pred), 3)
    tn, fp, fn, tp = confusion_matrix(y_var, y_pred).ravel()
    v["特异性"] = round(tn / (tn + fp), 3)
    v["精确度"] = round(precision_score(y_var, y_pred), 3)
    v["F1分数"] = round(f1_score(y_var, y_pred), 3)

    result[name] = v
    print(name + "模型计算完毕")

print(result)

临床曲线与预后

• 受试者工作特征曲线（Receiver Operating Characteristic Curve，ROCC），横轴为假阳性率 FPR纵轴为真阳性率 TPR，可判断不同诊断阈值α下预测精度。(0,0) 代表预测结果全部为阴；(0,1) 代表可完全判断正确；(1,0) 代表完全判断错误；(1,1) 代表预测结果全部为阳。即曲线越偏左上预测精度越高，反之右下精度最低，斜率为1的直线代表预测精度趋近于随机猜测

• AUC值（Area Under Curve），即ROC的曲线下面积，可以衡量该模型的优劣，愈趋近于1代表优势越大，反之优势越小

• 准确召回率曲线（Precision Recall Curve， PRC），横轴为精确率Precision纵轴为召回率Recall，两个变量均为愈趋近于1优势越大，故曲线越偏右上预测精度越高，反之左下精度越低，斜率为-1的直线代表预测精度趋近于随机猜测

  plt.figure(figsize=(10, 10)).clf()
  
  x_var = X_test
  y_var = y_test
  # 
  for name, model in models_list.items():
      model.fit(X, y)
      # 模型预测的正类概率
      y_pred = model.predict_proba(X)[:,1]
      # 计算PR曲线和AUC
      precision, recall, _ = precision_recall_curve(y, y_pred)
      pr_auc = auc(recall, precision)
  
      # 利用Booststrap方法计算AUC的95% CI
      bootstrapped_aucs = []
      rng = np.random.RandomState(0)  # 固定随机种子，保证可复现
  
      for _ in range(1000):
          # 有放回抽样（保持测试集的标签和预测概率对应）
          indices = rng.randint(0, len(y), len(y))
          if len(np.unique(y[indices])) < 2:
              # 避免抽样后只有一类标签（无法计算AUC）
              continue
          # 计算当前抽样的AUC
           # 计算抽样数据的PR曲线和AUC
          bs_precision, bs_recall, _ = precision_recall_curve(y[indices], y_pred[indices])
          bs_auc = auc(bs_recall, bs_precision)
          bootstrapped_aucs.append(bs_auc)
      
      # 计算95%置信区间（2.5%和97.5%分位数）
      sorted_aucs = np.array(bootstrapped_aucs)
      sorted_aucs.sort()
      auc_ci_lower = round(sorted_aucs[int(0.025 * len(sorted_aucs))], 3)
      auc_ci_upper = round(sorted_aucs[int(0.975 * len(sorted_aucs))], 3)
      # auc_ci = (auc_ci_lower, auc_ci_upper)  # 置信区间元组
      
      str_auc = '(AUC=' + "%.3f"%pr_auc + ', 95%CI:' + "%.3f"%auc_ci_lower + '-' + "%.3f"%auc_ci_upper + ')'
      plt.plot(recall, precision, label=name + str_auc)
  
  plt.plot([0,1],[1,0], "--")
  plt.legend(loc="lower left", fontsize=12)
  plt.grid(alpha=0.3)
  plt.show()

• 决策曲线分析（Decision Curve Analysis，DCA），横轴为阈值概率纵轴为净获益，代表某种因素达到阈值（x）时患者出现阳性症状的概率记为Pi，当Pi超过阈值Pt（判断临床阳性）时采取某种干预措施，干预后临床获益减去弊端即为净获益（y）。此外DCA还有两条虚线：横线代表所有样本均为阴性即无需处理净获益为0，斜线表示所有样本均为阳性且接受干预。同一变量阈值下，纵轴的净获益越高，代表临床决策获益越明显，即曲线愈偏向右上优势愈大，能够评估不同风险阈值下当前模型带来的净收益

• 临床影响曲线（Clinical Impact Curve，CIC）

SHAP可视化分析

计算单个模型的SHAP值，并设置解释器

import shap
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 设置模型
model = catboost.CatBoostClassifier(iterations=100,depth=6,random_seed=0)
model.fit(X,y)

# 设置SHAP解释器
explainer = shap.KernelExplainer(model.predict,X)
shap_values = explainer.shap_values(X,nsamples=100)

# 测试模块
# exp = explainer.explain(X_test)
# print(type(shap_values[1][0]))
# print(type(exp))

# 创建Explanation对象
shap_explanation = shap.Explanation(
    values=shap_values[1],  # 类别1的SHAP值
    base_values=explainer.expected_value,  # 类别1的基础值
    data=X_test,  # 原始特征值
    feature_names=features  # 特征名称
)

# 绘制SHAP的热点图和柱状图
shap.summary_plot(shap_values, X, feature_names=features, max_display=30)
shap.summary_plot(shap_values,X,feature_names=features, plot_type="bar", max_display=30)
# 绘制SHAP的力图
shap.initjs()
shap.force_plot(explainer.expected_value, shap_values[0,:],features)

---

参考资料

灵敏度和特异度 - 维基百科，自由的百科全书

决策曲线（DCA）和临床影响曲线（CIC）如何用于指导临床实践？ - 知乎

决策曲线分析法（Decision Curve Analysis，DCA）曲线 | Public Library of Bioinformatics

评价指标整理：Precision, Recall, F-score, TPR, FPR, TNR, FNR, AUC, Accuracy - 山竹小果 - 博客园

[Day 15] SHAP實作：實戰演練SHAP解釋方法 - iT 邦幫忙::一起幫忙解決難題，拯救 IT 人的一天