<jupyter_start><jupyter_text>1 import required package<jupyter_code>import os
import shap
import joblib
import sklearn
import numpy as np
import pandas as pd
import xgboost as xgb
import lightgbm as lgb
from sklearn import model_selection
from sklearn.svm import SVC
from sklearn.metrics import roc_auc_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
import matplotlib.pyplot as plt
from matplotlib import cm
def check_label(data):
    n = 0
    for (i,_,_,_,j) in data:
        n = n + j
    print('the count of different state:\n1->',n,'  0->',len(data)-n,'  total->',len(data))
    return {'1':n, '0':len(data)-n, 'total':len(data)}

def index_delete(data, index = None):
    data_new = []
    for i, sub in enumerate(data):
        if not (sub[3][:,index] == 0).all():
            sub[2] = sub[2][:, index]
            sub[3] = sub[3][:, index]
            data_new.append(sub)
    return data_new<jupyter_output><empty_output><jupyter_text>2 get the index and features name<jupyter_code>index = [0, 1, 2, 3, 4, 5, 6, 8, 10, 11, 12, 13, 
         16, 19, 21, 22, 25, 27, 28, 29, 31, 32,
         34, 35, 39, 40, 41, 42, 46, 48, 49, 50, 
         51, 53, 55, 56, 57, 58, 60, 61, 63, 64, 
         70, 73, 74, 75, 77, 78, 79, 80, 82, 83, 
         84, 85, 89, 91, 92, 98, 100, 102, 105, 
         106, 107, 110, 111, 115, 116, 117, 119, 
         120, 122, 124, 131, 133, 135, 136, 138, 
         140, 141, 144, 145, 146, 147, 150, 151, 
         152, 155, 157, 159, 160, 161, 162, 163, 
         164, 165, 166, 170, 171, 173]

items = ['PLT', 'INR', 'D-Dimer', 'FIB', 'PT', 'HBA1c', 'HGB', 
         'Na', 'EO%', 'APOB', 'CK-MB', 'pO2', 'TBIL', 'Ca', 
         'TG', 'Cys-c', 'TP', 'CHOI', 'TCO2', 'gGlu', 'CO2CP', 
         'gHGB', 'PTA', 'NRBC%', 'ABE', 'MCV', 'CHE', 'Cl-', 
         'CaO2', 'PDW', 'α-HBDH', 'CRE', 'RDW-CV', 'WBC', 
         'P-LCR', 'GA%', 'K+', 'AMY', 'SOD', 'CK', 'TBA', 
         'MCH', 'UA', 'Mg', 'Lac', 'GLB', 'HDL', 'FCOHb', 
         'LDH', 'Glu', 'IG%', 'BUN', 'Cl', 'SpO2', 'PA-aDO2', 
         'FDP', 'p50', 'MetHb', 'APOA', 'DBIL', 'MCHC', 'ALB', 
         'RBC', 'eGFR', 'MPV', 'PA', 'Na+', 'pH', 'HCT', 'GGT', 
         'APOE', 'ALP', 'FO2Hb', 'AG', 'PCT', 'APTT', 'K', 'AB', 
         'LPS', 'proBNP', 'pCO2', 'ALT', 'A/G', 'LDL', 'Lp(a)', 
         'RDW-SD', 'MONO%', 'NEUT%', 'TnT', 'LYMPH%', 'BASO%', 
         'procalcitonin', 'P', 'HHb', 'AST', 'IDBIL', 'Ca2+', 
         'T', 'TT']
index_mimic = [0, 1, 2, 3, 4, 5, 6, 8, 10, 12, 13, 
               16, 19, 21, 25, 27, 28, 29, 32, 35, 
               39, 40, 42, 50, 51, 53, 57, 58, 61, 
               64, 70, 73, 74, 75, 77, 78, 79, 80, 
               83, 84, 85, 89, 91, 98, 102, 105, 
               106, 107, 116, 117, 119,120, 124, 
               133, 136, 138, 140, 144, 145, 146, 
               150, 155, 157, 159,160, 161, 163, 
               165, 166, 170, 171, 173]

items_mimic = ['PLT', 'INR', 'D-Dimer', 'FIB', 'PT', 'HBA1c', 
               'HGB', 'Na', 'EO%', 'CK-MB', 'pO2', 'TBIL', 
               'Ca', 'TG', 'TP', 'CHOI', 'TCO2', 'gGlu', 
               'gHGB', 'NRBC%', 'ABE', 'MCV', 'Cl-', 'CRE', 
               'RDW-CV', 'WBC', 'K+', 'AMY', 'CK', 'MCH', 
               'UA', 'Mg', 'Lac', 'GLB', 'HDL', 'FCOHb', 
               'LDH', 'Glu', 'BUN', 'Cl', 'SpO2', 'PA-aDO2', 
               'FDP', 'MetHb', 'DBIL', 'MCHC', 'ALB', 'RBC', 
               'Na+', 'pH', 'HCT', 'GGT', 'ALP', 'AG', 'APTT', 
               'K', 'AB', 'proBNP', 'pCO2', 'ALT', 'LDL', 'MONO%', 
               'NEUT%', 'TnT', 'LYMPH%', 'BASO%', 'P', 'AST', 
               'IDBIL', 'Ca2+', 'T', 'TT']
if not os.path.exists('model/'):
    os.makedirs('model/')<jupyter_output><empty_output><jupyter_text>3 get the train, valid, test data<jupyter_code>dic_name = '../input/sample-xjtu/dic_3_1.npy'
# data = np.load(dic_name, allow_pickle=True)
data = np.load(dic_name, allow_pickle=True)
data = index_delete(data, index_mimic)
feature_num = data[0][2].shape[1]
print('features count:', feature_num)

lackitem = []
x = []
y = []
temp = []
for i in data:
    temp = np.array(i[2])
    x.append(np.reshape(temp, -1))
    y.append(np.array(i[4]))
x = np.array(x)
y = np.array(y).reshape(len(y))

items_real = []
time = [' at T1', ' at T2', ' at T3']
for i in range(3):
    for j in items_mimic:
        items_real.append(j + time[i])
data_pd = pd.DataFrame(np.hstack((x,y.reshape(-1,1))), columns = items_real + ['state'])
data_pd.to_csv('model/data.csv', index = None)
print('shape for x and y:', x.shape, y.shape, len(data))
x_train, x_test, y_train, y_test = model_selection.train_test_split(data_pd.iloc[:,:-1], data_pd.iloc[:,-1], test_size = 0.2, random_state = 42)
x_train, x_valid, y_train, y_valid = model_selection.train_test_split(x_train, y_train, test_size = 0.125, random_state = 42)
check_label(data)
result = {}<jupyter_output>features count: 72
shape for x and y: (7847, 216) (7847,) 7847
the count of different state:
1-> 3046.0   0-> 4801.0   total-> 7847<jupyter_text>4 train and save model  4.1 LR model<jupyter_code>def LR_model(x_train, y_train, x_test, y_test):
    model = LogisticRegression()  
    model.fit(x_train,y_train)
    pred_train = model.predict_proba(x_train)[:,-1]
    pred_test = model.predict_proba(x_test)[:,-1]
    result = [roc_auc_score(y_train,pred_train), roc_auc_score(y_test,pred_test)]
    return model, result, (y_test.values, pred_test)

print('running for LR:\n')
model_LR,auc_LR, result_LR = LR_model(x_train, y_train, x_test, y_test)
result['LR'] = result_LR
joblib.dump(model_LR, 'model/LR_99.model')
print(auc_LR)<jupyter_output><empty_output><jupyter_text>4.2 SVM model<jupyter_code>def SVM_model(x_train, y_train, x_test, y_test):
    model = SVC(kernel='rbf', probability=True)  
    model.fit(x_train,y_train)  
    pred_train = model.predict_proba(x_train)[:,-1]
    pred_test = model.predict_proba(x_test)[:,-1]
    result = [roc_auc_score(y_train,pred_train), roc_auc_score(y_test,pred_test)]
    return model, result, (y_test.values, pred_test)

print('running for SVM:\n')
model_SVM, auc_SVM, result_SVM = SVM_model(x_train, y_train, x_test, y_test)
result['SVM'] = result_SVM
joblib.dump(model_SVM, 'model/SVM_99.model')
print(auc_SVM)<jupyter_output><empty_output><jupyter_text>4.3 Xgooost model<jupyter_code>def Xgboost_model(x_train, y_train, x_test, y_test):
    xgb_train = xgb.DMatrix(x_train, label=y_train)
    xgb_valid = xgb.DMatrix(x_valid, label=y_valid)
    xgb_test = xgb.DMatrix(x_test)
    params = {'booster': 'gbtree',
            'objective': 'rank:pairwise',  # 二分类的问题
            'gamma': 0.5,  # 用于控制是否后剪枝的参数,越大越保守，一般0.1、0.2这样子。
            'max_depth': 15,  # 构建树的深度，越大越容易过拟合
            #'lambda': 0.5,  # 控制模型复杂度的权重值的L2正则化项参数，参数越大，模型越不容易过拟合。
            'subsample': 0.5,  # 随机采样训练样本
            'colsample_bytree': 0.5,  # 生成树时进行的列采样
            'min_child_weight': 1,
            # 这个参数默认是 1，是每个叶子里面 h 的和至少是多少，对正负样本不均衡时的 0-1 分类而言
            # ，假设 h 在 0.01 附近，min_child_weight 为 1 意味着叶子节点中最少需要包含 100 个样本。
            # 这个参数非常影响结果，控制叶子节点中二阶导的和的最小值，该参数值越小，越容易 overfitting。
            'scale_pos_weight': 1,
            'silent': 0,  # 设置成1则没有运行信息输出，最好是设置为0.
            'eta': 0.015,  # 如同学习率
            'nthread': 8,  # cpu 线程数
            'eval_metric': 'logloss'  # 评价方式
            }

    plst = list(params.items())
    num_rounds = 500  # 迭代次数
    watchlist = [(xgb_valid, 'valid')]
    model = xgb.train(plst, xgb_train, num_rounds, watchlist, early_stopping_rounds=200)
    # 当设置的迭代次数较大时，early_stopping_rounds 可在一定的迭代次数内准确率没有提升就停止训练
    # model = xgb.train(plst, xgb_train, num_rounds)
    pred_train = model.predict(xgb_train)
    pred_test = model.predict(xgb_test)
    result = [roc_auc_score(y_train,pred_train), roc_auc_score(y_test,pred_test)]
    return model, result, (y_test.values, pred_test)

print('running for Xgboost:\n')
model_Xgboost, auc_Xgboost, result_Xgboost = Xgboost_model(x_train, y_train, x_test, y_test)
result['Xgboost'] = result_Xgboost
joblib.dump(model_Xgboost, 'model/Xgboost_99.model')
print(auc_Xgboost)<jupyter_output>running for Xgboost:

[09:51:07] WARNING: ../src/learner.cc:573: 
Parameters: { "scale_pos_weight", "silent" } might not be used.

  This may not be accurate due to some parameters are only used in language bindings but
  passed down to XGBoost core.  Or some parameters are not used but slip through this
  verification. Please open an issue if you find above cases.


[0]	valid-logloss:0.67674
[1]	valid-logloss:0.66205
[2]	valid-logloss:0.64942
[3]	valid-logloss:0.63667
[4]	valid-logloss:0.62485
[5]	valid-logloss:0.61254
[6]	valid-logloss:0.60057
[7]	valid-logloss:0.58898
[8]	valid-logloss:0.57773
[9]	valid-logloss:0.56977
[10]	valid-logloss:0.55867
[11]	valid-logloss:0.54614
[12]	valid-logloss:0.53736
[13]	valid-logloss:0.52533
[14]	valid-logloss:0.51764
[15]	valid-logloss:0.50876
[16]	valid-logloss:0.49947
[17]	valid-logloss:0.48968
[18]	valid-logloss:0.48233
[19]	valid-logloss:0.47468
[20]	valid-logloss:0.46612
[21]	valid-logloss:0.45644
[22]	valid-logloss:0.44852
[23]	valid-logloss:[...]<jupyter_text>4.4 Lightgbm model<jupyter_code>def Lightgbm_model(x_train, y_train, x_test, y_test):
    lgb_train = lgb.Dataset(x_train, label=y_train)
    lgb_valid = lgb.Dataset(x_valid, label=y_valid)
    params = {
        'boosting_type': 'gbdt',
        'objective': 'binary',
        # 'metric': {'auc','binary_logloss'},
        'metric': {'binary_logloss'},
        'num_leaves': 30,
        'max_depth': 15,
        'max_bin': 50,
        'min_data_in_leaf': 15,
        'learning_rate': 0.015,
        'feature_fraction': 0.5,
        'bagging_fraction': 0.5,
        'bagging_freq': 5,
        'scale_pos_weight': 50
    }

    MAX_ROUNDS = 500
    # model = lgb.train(params, lgb_train, num_boost_round=MAX_ROUNDS)
    model = lgb.train(params, lgb_train, valid_sets=lgb_valid, 
                      num_boost_round=MAX_ROUNDS,early_stopping_rounds=200)
    # pred_train = model.predict(x_train, num_iteration = model.best_iteration)
    # pred_test = model.predict(x_test, num_iteration = model.best_iteration)
    pred_train = model.predict(x_train)
    pred_test = model.predict(x_test)
    result = [roc_auc_score(y_train,pred_train), roc_auc_score(y_test,pred_test)]
    return model, result, (y_test.values, pred_test)
print('running for Lightgbm:\n')
model_Lightgbm, auc_Lightgbm, result_Lightgbm = Lightgbm_model(x_train, y_train, x_test, y_test)
joblib.dump(model_Lightgbm, 'model/Lightgbm_99.model')
model_Lightgbm.save_model('model/Lightgbm_99.txt')
result['Lightgbm'] = result_Lightgbm
print(auc_Lightgbm)<jupyter_output><empty_output><jupyter_text>4.5 save results<jupyter_code>np.save('model/result.npy', result)
model = ['LR', 'SVM', 'Xgboost', 'Lightgbm']
auc = np.array([auc_LR, auc_SVM, auc_Xgboost, auc_Lightgbm])
auc_pd = {'model': model, 'auc_train':auc[:,0], 'test_auc':auc[:,1]}
auc_pd = pd.DataFrame(auc_pd)
auc_pd.to_csv('model/auc.csv')<jupyter_output><empty_output><jupyter_text>5 plot features importance for Xgboost model<jupyter_code>importance_value = pd.DataFrame(model_Xgboost.get_fscore().items(), columns= ['features_names', 'F_score'])
importance_value.to_csv('model/features_importance.csv', index = None)
fig, ax = plt.subplots(1,1,figsize=(5,10))
xgb.plot_importance(model_Xgboost, ax=ax, max_num_features=20, 
                    grid=False, xlabel='Features contributiton score')
plt.savefig('model/features_importance.svg', bbox_inches = 'tight')
def draw_bar(key_name,key_values,size=(5,8), N = 20):
    norm = plt.Normalize(min(key_values),max(key_values))
    norm_values = norm(key_values[-N:])
    map_vir = cm.get_cmap(name='coolwarm')
    colors = map_vir(norm_values)
    fig,ax = plt.subplots(1,1, figsize=size) #调用figure创建一个绘图对象
    plt.barh(key_name[-N:], key_values[-N:],height=0.5,color=colors,edgecolor='black') # edgecolor边框颜色  
    # ax.spines['left'].set_visible(False)
    ax.spines['right'].set_visible(False) 
    ax.spines['top'].set_visible(False) 
    # ax.spines['bottom'].set_visible(False) 
    sm = cm.ScalarMappable(cmap=map_vir, norm = norm)  # norm设置最大最小值
    plt.xlabel('Feature contributiton (%)')
    plt.ylabel('Features')
    plt.title('Feature importance')
    
    position=fig.add_axes([0.9, 0.15, 0.035, 0.695])#位置[左,下,右,上]
    cbar=plt.colorbar(sm,cax=position)
    # cbar.ax.set_yticklabels(['1','2'])
    plt.savefig('model/features_importance_colorbar.svg', bbox_inches = 'tight')


importance = pd.read_csv('model/features_importance.csv').values
index_sorted = np.argsort(importance[:,1])
importance_relative = importance[index_sorted,1]/sum(importance[index_sorted,1])*100
draw_bar(importance[index_sorted,0], importance_relative.tolist(), N = 20)<jupyter_output><empty_output><jupyter_text>6 shap explanation research for Xgboost model<jupyter_code>shap.initjs()
explainer = shap.TreeExplainer(model_Xgboost)
shap_values = explainer.shap_values(x_train)
print(shap_values.shape)
if not os.path.exists('shap/'):
    os.makedirs('shap/')
y_base = explainer.expected_value
print(y_base)
xgb_train = xgb.DMatrix(x_train, label=y_train)
y_pred = model_Xgboost.predict(xgb_train)
print(y_pred.mean())
j = 9

player_explainer = {}
player_explainer['feature'] = x_train.columns
player_explainer['feature_value'] = x_train.values[j]
player_explainer['shap_value'] = shap_values[j]
player_explainer = pd.DataFrame(player_explainer)
player_explainer.to_csv('shap/negative_sample.csv', index = None)

f = plt.gcf()
features = [col + ' = ' + str(round(x_train.values[j, ind], 2)) 
            for ind, col in enumerate(x_train.columns)]
t = shap.force_plot(explainer.expected_value, shap_values[j], 
                    feature_names=features,
                    out_names='output value', 
                    text_rotation = 0, show = False, matplotlib=True)
t.savefig('shap/base_value_negative.svg', bbox_inches = 'tight')
j = 4996

player_explainer = {}
player_explainer['feature'] = x_train.columns
player_explainer['feature_value'] = x_train.values[j]
player_explainer['shap_value'] = shap_values[j]
player_explainer = pd.DataFrame(player_explainer)
player_explainer.to_csv('shap/postive_sample.csv', index = None)

f = plt.gcf()
features = [col + ' = ' + str(round(x_train.values[j, ind], 2)) 
            for ind, col in enumerate(x_train.columns)]
t = shap.force_plot(explainer.expected_value, shap_values[j], 
                    feature_names=features,
                    out_names='output value', 
                    text_rotation = 0, show = False, matplotlib=True)
t.savefig('shap/base_value_postive.svg', bbox_inches = 'tight')
f = plt.gcf()
shap.summary_plot(shap_values, x_train, show = False)
# shap.summary_plot(shap_values, x_train, max_display=20, show = False)
plt.savefig('shap/shap_value.svg',bbox_inches = 'tight')
shap.dependence_plot('PLT at T3', shap_values, x_train, interaction_index=None, show=False)
plt.savefig('shap/shap_value_for_plt.svg', bbox_inches='tight')
shap.dependence_plot('PLT at T2', shap_values, x_train, interaction_index='D-Dimer at T2', show=False)
plt.savefig('shap/shap_value_for_plt_and_D_Dimer.svg', bbox_inches='tight')
shap.dependence_plot('PLT at T2', shap_values, x_train, interaction_index='PCT at T2', show=False)
plt.savefig('shap/shap_value_for_plt_and_pct.svg', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16496521230249042, 'fraction_numerical': 0.05344950907081084, 'mean_word_length': 3.893413935730194, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '106041', 'license': 'MIT', 'n_tokens_mistral': 7041, 'n_tokens_neox': 6384, 'n_words': 1198, 'path': 'sepsis_machine_learning.ipynb', 'repo_name': 'huawenbo/ERSW-project', 'size': '15379'}