WiDS Datathon Maastricht 2021

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import seaborn as sns

import pandas as pd
df = pd.read_csv('TrainingWiDS2021.csv', index_col=0)
df_pred = pd.read_csv('UnlabeledWiDS2021.csv', index_col=0)

print('train data:', df.shape, 'test data:', df_pred.shape)
df.head()

# a short view of missing values and unique values in the data, as well as correlation with the target
y_df = df['diabetes_mellitus']
df = df.drop(columns=['diabetes_mellitus'])
def highlight_nas(s):
    '''
    highlight values above 0.5 in a Series yellow.
    '''
    is_max = s > 0.5
    return ['background-color: yellow' if v else '' for v in is_max]

dftypes = pd.DataFrame(df.dtypes.rename('type'))
dftypes['isnull'] = df.isnull().mean()
dftypes['unique'] = df.nunique()
dftypes['corr'] = df.corrwith(y_df)
dftypes['abs_corr'] = dftypes['corr'].apply(abs)
dftypes.sort_values(by='abs_corr', ascending=False)[['type','isnull','unique','corr']][0:15]\
        .style.apply(highlight_nas, subset=['isnull']).set_caption("Tab.1")

#uniting both train and test data for feature creation and data exploration
test_id = df_pred['encounter_id']

train_size = df.shape[0]
test_size = df_pred.shape[0]
X = pd.concat([df, df_pred], axis=0)
X = X.drop(columns=['encounter_id','hospital_id'])
X.head()

%matplotlib inline
# comparing the distribution of train and test set for 3 features. try out other features as well.
import matplotlib.pyplot as plt


X['dataset']='train'
X['dataset'].iloc[train_size:]='test'
fig, ax = plt.subplots(1,3, figsize=(12,5))
X.groupby(['dataset'])['age'].plot.kde(title='Distribution of age', ax=ax[0])
X.groupby(['dataset'])['weight'].plot.kde(title='Distribution of weight', ax=ax[1])
X.groupby(['dataset'])['bmi'].plot.kde(title='Distribution of bmi', ax=ax[2])#, ind=50)
X.drop(columns=['dataset'],inplace=True)
plt.legend()
plt.suptitle( 'Fig. 1')
plt.show();
y_df.value_counts().plot.pie(title='diabetes mellitus')
plt.suptitle( 'Fig. 2')

plt.show();

import gc
del (display)

# correlation of the maximum examination values for the first 24 hours. Only those with a correlation
lab_col = [c for c in X.columns if((c.startswith("h1")) | (c.startswith("d1")))]
lab_max = [col for col in lab_col if (col.endswith('max'))&(col.startswith('d1'))]
M = df[lab_max].corr()
cond = (M>0.3) | (M<-0.3)
M = M.loc[(cond).sum(axis=1)>=2,(cond).sum(axis=0)>=2]
print(M.shape)
with pd.option_context('display.max_rows', 100, 'display.max_columns', 100):
    display(M.style.background_gradient(cmap='Blues'))

#correlation heatmap
fig, heat = plt.subplots(figsize = (12,8))

heat = sns.heatmap(M, annot=True, fmt=',.2f', center=0, cbar=False, annot_kws={'size': 8})#
plt.suptitle( 'Fig. 3 - Correlations between several maximum values of the first day');

# cleaning errors in the data, by exchanging maximum and minimum values, when min is larger than max.
min_max_feats=[f[:-4] for f in df.columns if f[-4:]=='_min']
for col in min_max_feats:
    X.loc[X[f'{col}_min'] > X[f'{col}_max'], [f'{col}_min', f'{col}_max']] = X.loc[X[f'{col}_min'] > X[f'{col}_max'], [f'{col}_max', f'{col}_min']].values

# selecting categorical values
hospital_categorical = ['hospital_admit_source','icu_admit_source','icu_stay_type','icu_type','icu_id']
diseases = ['aids', 'cirrhosis','hepatic_failure', 'immunosuppression', 'leukemia','lymphoma','solid_tumor_with_metastasis']
personal_categorical = ['gender','ethnicity']
cat_features = hospital_categorical+diseases+personal_categorical

# encoding categories, but only those which are in test and training set. the rest is set to -1.
for col in cat_features:
    cats_train = set(X[col].iloc[:train_size].unique())
    cats_test = set(X[col].iloc[train_size:].unique())
    categories = cats_train.intersection(cats_test)
    categories = [k for k in categories if k==k] #nan is never equal to itself :)
    print(col, cats_train.difference(cats_test), cats_test.difference(cats_train))

    X[col] = pd.Categorical(X[col], categories=categories).reorder_categories(categories)
    X[col] = X[col].cat.codes

apache_cols = [col for col in X.columns if col.endswith('apache')]
d1_cols = [col for col in X.columns if col.startswith('d1_')]
base = set(col[3:-4] for col in d1_cols)

# the quality of kNN imputation. Additionally to the missing values, we impute known unseen values, to measure the imputation quality by the prediction on these.
def kNN_score(knn_model, Xfit, yfit):
    """ checks how well a kNN model can predict the target in yfit given the features in Xfit
    returns: the average 3-fold cross-validation score (r^2)"""

    kfold = KFold(n_splits=3)
    avg_score=0
    for train_index, test_index in kfold.split(Xfit, yfit):
        X_train, X_test = Xfit.iloc[train_index], Xfit.iloc[test_index]
        y_train, y_test = yfit.iloc[train_index], yfit.iloc[test_index]
        knn_model.fit(X_train, y_train)
        score = knn_model.score(X_test, y_test, sample_weight=None)
        print(score,end=', ')
        avg_score+=score/3
    print('score = ', avg_score)
    return avg_score

# decision whether or not a value is imputed depends on 2 things, the quality of imputation, and how many values can be imputed with the model
def try_to_impute_target(X, target, features, min_score=0.85):
    X_k = X[features]
    y_k = X[target]
    cond = (~X_k.isnull().any(axis=1))
    cond_test = cond&(y_k.isnull())
    cond_train = cond&(~y_k.isnull())
    Xfit = X_k[cond_train]
    yfit = y_k[cond_train]
    imputable_rows = cond_test.sum()
    training_rows = cond_train.sum()
    print(target, imputable_rows, training_rows )
    if imputable_rows>100:
        knn_model = KNeighborsRegressor(n_neighbors=50)
        if kNN_score(knn_model, Xfit, yfit)>min_score:
            print('imputing', target)
            knn_model.fit(Xfit, yfit)
            X['imp_'+target] = y_k#.isnull().astype(int)
            X['imp_'+target][cond_test] = knn_model.predict(X_k[cond_test])

# we set the features on which the values we want to impute depend. These are the features that kNN uses to find the closest values.
imputers = {

'h1_bun_max':['d1_bun_max','d1_bun_min'],#,'d1_bun_mean'],#0.973
'h1_bun_min':['d1_bun_max','d1_bun_min','bun_apache','h1_bun_max'],#'d1_bun_mean',#0.993
'd1_hemaglobin_min':['d1_wbc_max','d1_wbc_min','d1_hematocrit_max','d1_hematocrit_min'],#,'d1_wbc_mean',,'d1_hematocrit_mean'0.932
'd1_hemaglobin_max':[ 'd1_wbc_max','d1_wbc_min','d1_hematocrit_max','d1_hematocrit_min','d1_hemaglobin_min'],#,'d1_wbc_mean',,'d1_hematocrit_mean'#0.942
'd1_hematocrit_min':['d1_hemaglobin_max','d1_hemaglobin_min'],#'d1_hemaglobin_mean'],#0.941
'd1_hematocrit_max':['d1_hemaglobin_max','d1_hemaglobin_min', 'd1_hematocrit_min'],#,'d1_hemaglobin_mean'0.968

'glucose_apache':['d1_glucose_max','d1_glucose_min','d1_creatinine_max','d1_albumin_max','d1_bun_max' ],#,'d1_glucose_mean',,'d1_glucose_range',0.899
'h1_hco3_min':['d1_hco3_max','d1_creatinine_max', 'd1_hco3_min' ],#'d1_hco3_mean',0.855
'h1_hco3_max':['d1_hco3_max','d1_creatinine_max', 'd1_hco3_min','h1_hco3_min' ],#'d1_hco3_mean',
'h1_creatinine_min':[ 'd1_creatinine_max', 'd1_creatinine_min', 'creatinine_apache',
                      'd1_hco3_max', 'd1_bun_max','h1_hco3_max' ],#0.848
'h1_creatinine_max':[ 'd1_creatinine_max', 'd1_creatinine_min', 'creatinine_apache',
                      'd1_hco3_max', 'd1_bun_max','h1_hco3_max' ,'h1_creatinine_min'],#0.874
'h1_arterial_ph_max':['d1_arterial_ph_max','d1_arterial_ph_min','ph_apache']#,'d1_hco3_min','d1_lactate_max''d1_arterial_pco2_max'],'paco2_for_ph_apache'
}
           #

from sklearn.model_selection import KFold
from sklearn.neighbors import KNeighborsRegressor

for target, features in imputers.items():

    try_to_impute_target(X, target, features)

# second round of imputations, with features that are not hand selected.
imputers={}

for col in base:
    if col+'_apache' in apache_cols:
        imputers[col+'_apache']=['d1_'+col+'_max','d1_'+col+'_min']
        imputers['h1_'+col+'_min']=['d1_'+col+'_max','d1_'+col+'_min',col+'_apache']
        imputers['h1_'+col+'_max']=['d1_'+col+'_max','d1_'+col+'_min',col+'_apache','h1_'+col+'_min']
    else:
        imputers['h1_'+col+'_min']=['d1_'+col+'_max','d1_'+col+'_min']
        imputers['h1_'+col+'_max']=['d1_'+col+'_max','d1_'+col+'_min','h1_'+col+'_min']

for target, features in imputers.items():
    try_to_impute_target(X, target, features)
#problems with: h1_sysbp_min,

# new features will be the average and difference between the max and min values of each measure
for col in lab_col:
   # X['glucomaxdiv'+col] = X['d1_glucose_max']/X[col]
  #  X['glucomindiv'+col] = X['d1_glucose_min']/X[col]
    if col.endswith('max'):
        X[col[:-4]+'_range'] = X[col]-X[col[:-4]+'_min']
        X[col[:-4]+'_mean'] = (X[col]+X[col[:-4]+'_min'])/2

# relating examination values to the average value of patients in the same weight group
def relate_to_meanweight(X_df, col, colnew):

    albucol = X_df.groupby('weightgroup')[col].mean()
    X_df[colnew] = np.nan
    cond = ~X_df[col].isnull()& ~X_df.weightgroup.isnull()
    X_df[colnew][cond] = X_df[cond][[col,'weightgroup']].apply(lambda x: x[col]/albucol.loc[x['weightgroup']], axis=1)
    return X_df

# relating examination values to the median value of patients in the same age group
def relate_to_medianage(X_df, col, colnew):

    albucol = X_df.groupby('age')[col].median()
    X_df[colnew] = np.nan
    cond = ~X_df[col].isnull()& ~X_df.age.isnull()
    X_df[colnew][cond] = X_df[cond][[col,'age']].apply(lambda x: x[col]/albucol.loc[x['age']], axis=1)
    return X_df

# creating new features from relations of various examination values
def get_extra_cols1(X_df):
    X_df = relate_to_medianage(X_df, 'd1_albumin_mean','albuage')
    X_df = relate_to_medianage(X_df, 'd1_creatinine_mean','creage')
    X_df = relate_to_medianage(X_df, 'd1_glucose_mean','glucoage')
    #X_df['albuweight'] = X_df['d1_albumin_mean']/X_df['weight']
    X_df['weightgroup']=pd.cut(X_df.weight,[0,20,30,40,50,60,70,80,140])
    X_df = relate_to_meanweight(X_df, 'd1_albumin_mean','albuweight')
    X_df = relate_to_meanweight(X_df, 'd1_creatinine_mean','creaweight')
    X_df = relate_to_meanweight(X_df, 'd1_glucose_mean','glucoweight')
    return X_df

def get_extra_cols2(X_df, suff='mean'):
    X_df['bilialbu'+suff] = X_df['d1_bilirubin_'+suff]/X_df['d1_albumin_'+suff]
    X_df['albucalc'+suff] = X_df['d1_albumin_'+suff]/X_df['d1_calcium_'+suff]
    X_df['creatigluc'+suff] = X_df['d1_creatinine_'+suff]/X_df['d1_glucose_'+suff]
    X_df['buncrea'+suff] = X_df['d1_bun_'+suff]/X_df['d1_creatinine_'+suff]
    X_df['glucohem'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_hemaglobin_'+suff]
    X_df['glucohemo2'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_hemaglobin_'+suff]*X_df['d1_hco3_'+suff]
    X_df['hemahema'+suff] = X_df['d1_hemaglobin_'+suff]/X_df['d1_hematocrit_'+suff]
    X_df['hemoglucocrit'+suff] = X_df['d1_hemaglobin_'+suff]/X_df['d1_hematocrit_'+suff]/X_df['d1_glucose_'+suff]
    X_df['lactohco3'+suff] = X_df['d1_lactate_'+suff]*X_df['d1_hco3_'+suff]
    X_df['gluspo2'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_spo2_'+suff]
    X_df['lactoglu'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_spo2_'+suff]+X_df['d1_lactate_'+suff]/X_df['d1_hco3_'+suff]
    X_df['gluco_inr'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_inr_'+suff]
    X_df['hco3pco2'+suff] = X_df['d1_hco3_'+suff]/X_df['d1_arterial_pco2_'+suff]
    X_df['lactogludiv'+suff] = X_df['d1_lactate_'+suff]/X_df['d1_glucose_'+suff]
    X_df['glucobun'+suff] = X_df['d1_glucose_'+suff]/X_df['d1_bun_'+suff]
    return X_df

def get_extra_cols3(X_df, suff='mean'):
    X_df['weight/height'] = X_df['weight']/X_df['height']
    X_df['bilibmi'] = X_df['d1_bilirubin_mean']/X_df['bmi']
    X_df['creatibmi'] = X_df['d1_creatinine_mean']/X_df['bmi']
    X_df['wbcbmi']=X_df['d1_wbc_mean']/X_df['bmi']
    X_df = X_df.drop(columns=['weightgroup'])

    return X_df
X = get_extra_cols1(X)
X = get_extra_cols2(X, suff='mean')
X = get_extra_cols2(X, suff='max')
X = get_extra_cols3(X)

# splitting the concatenated train and test set again into train and test
train_df = X.iloc[:train_size]
test_df = X.iloc[train_size:]

# splitting the training data into train and test
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(train_df, y_df, test_size=0.1,  random_state=42, shuffle=True, stratify=y_df)

import lightgbm as lgb
from sklearn import metrics
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

from scipy import stats, special

pds = {'num_leaves': (10, 500),
       'feature_fraction': (0.1, 0.5),
       'feature_fraction_bynode':(0.1, 1),
       'learning_rate':(0.008, 0.05),
       'max_depth': (6, 1000 ),
       'min_data_in_leaf':(10, 250),
       'lambda_l2':(0,5),
       'lambda_l1':(0,5)
         }

# objective function for hyper parameter optimization
def hyp_lgbm(num_leaves,  feature_fraction, feature_fraction_bynode,  learning_rate, max_depth, min_data_in_leaf, lambda_l2, lambda_l1):

        params = {'objective':'binary',
                  'early_stopping_round':100,
                  'num_iterations': 500,
                  'metric':'auc',
                  'verbose':-1,
                  'is_unbalance': True,
                 'min_gain_to_split':0.01,
                 'feature_pre_filter':False}
        params["num_leaves"] = int(round(num_leaves))
        params['feature_fraction'] = max(min(feature_fraction, 1), 0)
        params['feature_fraction_bynode'] = max(min(feature_fraction_bynode, 1), 0)
        params['learning_rate']=learning_rate
        params['max_depth'] = int(round(max_depth))
        params['min_data_in_leaf']  = int(round(min_data_in_leaf))
        params['lambda_l2'] = lambda_l2
        params['lambda_l1'] = lambda_l1

        cv_results = lgb.cv(params, train_data, nfold=5, seed=17, categorical_feature=cat_features,
                            verbose_eval=None, stratified=True)
        iters = len(cv_results['auc-mean'])
        print(iters, np.mean(cv_results['auc-mean'][-100:]),np.std(cv_results['auc-mean'][-100:]))
        params['num_iterations']=iters
        best_iter.append(iters)
        return np.max(cv_results['auc-mean'])

# hyperparameter optimization
from bayes_opt import BayesianOptimization
best_iter=[]
train_data = lgb.Dataset(X_train, label=y_train, categorical_feature=cat_features)
optimizer = BayesianOptimization(hyp_lgbm, pds, random_state=7)
optimizer.maximize(init_points=5, n_iter=10)

# using the parameters of the best round of hyperparameter optimization for future models
param = {'objective':'binary',
                  'metric':'auc',
                  'verbose':-1,
                 # 'is_unbalance': True,
                 'min_gain_to_split':0.01,
                 'feature_pre_filter':False
        }

best_params=optimizer.max['params']

param.update(best_params)

best_params_rounded = {
        'num_leaves': round(best_params['num_leaves']),
         'feature_fraction':max(min(best_params['feature_fraction'], 1), 0),
         'feature_fraction_bynode':max(min(best_params['feature_fraction_bynode'], 1), 0),
         'max_depth':int(round(best_params['max_depth'])),
        'min_data_in_leaf':int(round(best_params['min_data_in_leaf']))
         }
param.update(best_params_rounded)
param['metric'] = ['auc', 'binary_logloss']
param['num_iterations'] = best_iter[np.argmax([k['target'] for k in optimizer.res])]

# training model with set of optimized hyperparameters on our training set and testing on our test set
train_data = lgb.Dataset(X_train, label=y_train, categorical_feature=cat_features)#, weight=w)
bst = lgb.train(param, train_data, categorical_feature=cat_features)#, valid_sets=valid_data, early_stopping_rounds=10)#feature_name=list(X_df.columns),

ypred = bst.predict(X_test)
fpr, tpr, thresholds = metrics.roc_curve(y_test, ypred)
roc_auc = metrics.auc(fpr, tpr)
print(roc_auc) #0.87233

# area under the receiver operating characteristic (ROC) curve as quantifier of model quality
rocplot = metrics.RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=roc_auc, estimator_name='diabetes estimator')

rocplot.plot()
plt.suptitle( 'Fig. 4 - ROC curve for diabetes estimator');

import shap

shap.initjs()

# shapley values to measure feature importance
shap_values = shap.TreeExplainer(bst).shap_values(X_test)

shap.summary_plot(shap_values, X_test)
bst.feature_importance('gain')
plt.suptitle( 'Fig. 5 - Feature importance from shapley values');

# using shapley values to chose the best features for the model
con_cols = zip(X_train.columns,abs(shap_values[0]).mean(axis=0) )
#sorted(con_cols, key=lambda x: x[1], reverse=True)[0:400]
concols = [col for col,val in con_cols if val>0.001]
pro_cols = zip(X_train.columns,abs(shap_values[1]).mean(axis=0) )
#procols = sorted(pro_cols, key=lambda x: x[1], reverse=True)[0:400]
procols = [col for col,val in pro_cols if val>0.00]
bestcon = X_train.columns[list(set(abs(shap_values[0]).argmax(axis=1)))]
bestpro = X_train.columns[list(set(abs(shap_values[1]).argmax(axis=1)))]
take_best_features2 = set(concols).union(set(procols)).union(set(bestcon)).union(set(bestpro))
len(take_best_features2)

# using feature importance of the lgbm trees by looking at features that created the highest gain.
# comparing this to the best features of a model with randomized target and using those features which are
# better with the real target than with the randomized target
y_rand = np.array(y_train).copy()
np.random.shuffle(y_rand)
train_data = lgb.Dataset(X_train, label=y_rand, categorical_feature=cat_features)#, weight=w)
bst_rand = lgb.train(param, train_data, categorical_feature=cat_features)
rand_imp1 = bst_rand.feature_importance( importance_type='gain')
np.random.shuffle(y_rand)
bst_rand = lgb.train(param, train_data, categorical_feature=cat_features)
rand_imp2 = bst_rand.feature_importance( importance_type='gain')

best_feats = zip(bst.feature_importance( importance_type='gain') - (rand_imp1+rand_imp2)/2, X_train.columns)
sorted_best = sorted(best_feats, key=lambda x: x[0], reverse=True)
take_best_features = [a[1] for a in sorted_best if a[0]>0]
best_values =  [a[0] for a in sorted_best if a[0]>0]

print(X_train.shape[1], len(take_best_features))

# plotting the lgbm feature importance for the target, the randomized target and the difference between both.
fig, axes = plt.subplots(1,3, figsize=(35,10))
lgb.plot_importance(bst, importance_type='gain', max_num_features=20,  ax=axes[0]);
lgb.plot_importance(bst_rand, importance_type='gain', max_num_features=20,  ax=axes[1], title='feature importance for randomized target');
axes[2].barh(y=np.flipud(take_best_features[:20]), width=np.flipud(best_values[:20]), height=0.2)
axes[2].set_title('new feature importance')
axes[2].grid(True)
plt.suptitle( 'Fig. 6 - Feature importance by gain')
plt.show;


#print(len(take_best_features))

#
take_best_features = list(take_best_features)
cat_best_features = list(set(cat_features).intersection(set(take_best_features)))

# training and validating the model again with subset of selected features
train_data = lgb.Dataset(X_train[take_best_features], label=y_train, categorical_feature=cat_best_features)
bst = lgb.train(param, train_data, categorical_feature=cat_best_features)

ypred = bst.predict(X_test[take_best_features], num_iteration=bst.best_iteration)
fpr, tpr, thresholds = metrics.roc_curve(y_test, ypred)
metrics.auc(fpr, tpr) #0.86917

explainer = shap.TreeExplainer(bst)
shap_values = explainer.shap_values(X_test[take_best_features])

shap.summary_plot(shap_values, X_test[take_best_features])

expected_value = explainer.expected_value
if isinstance(expected_value, list):
    expected_value = expected_value[1]
print(f"Explainer expected value: {expected_value}")

select = range(50)
features = X_test[take_best_features].iloc[select]
shap_values_sub = explainer.shap_values(features)[1]

y_pred = (shap_values_sub.sum(1) + expected_value) > 0
misclassified = y_pred != y_test.iloc[select]

#shapley decision plot
shap.decision_plot(expected_value, shap_values_sub, features,  highlight=misclassified);
plt.suptitle( 'Fig. 7 - Decision  plot');
plt.show();

# shapley force plot
k = np.array(select)[misclassified][3]
shap.force_plot(expected_value, shap_values_sub[k,:], X_test[take_best_features].iloc[k,:], link='logit')

# shapley dependence plots
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(23,5))
shap.dependence_plot("d1_glucose_max", shap_values[1], X_test[take_best_features], ax=axes[0], show=False, title='(a)')
shap.dependence_plot('age', shap_values[1], X_test[take_best_features], ax=axes[1], show=False, title='(b)')
shap.dependence_plot("bmi", shap_values[1], X_test[take_best_features], ax=axes[2], show=False, title='(c)')
fig.suptitle( 'Fig. 8')
plt.show()

# defining hyper parameters for the final model training (chosen from previous hyper parameter testing rounds)
evals = {}
evals['lgb'] = {}
lparams = {}

seeds = [2, 3, 4]
nfold = [5, 5, 5]
lr = [0.03, 0.05, 0.01]
leaves = [77, 50, 280]
max_depth = [-1, 300, 300]
ff = [0.35, 0.22, 0.3]
ff_bynode = [0.65, 0.6, 0.5]
stratified=[False, True, True]
shuffle=[True, True, True]
n_estimators = [600, 450, 1600]
min_data_in_leaf = [20,20,240]
min_child_weight=[45, 45, 1e-3]
early_stopping_rounds = 100
verbose_eval = 0
reg_alpha = [0.7, 1, 0]
reg_lambda = [0.8, 1, 19]
verbose = -1
n_jobs = 4

for k in [0,1,2]:
    lparams[k] = dict(boosting_type='gbdt',
               objective='binary',
               metric='auc',
               learning_rate=lr[k],
               num_leaves=leaves[k],
               max_depth=max_depth[k],
               feature_fraction=ff[k],
               feature_fraction_bynode=ff_bynode[k],
               n_estimators=n_estimators[k],
               min_child_weight=min_child_weight[k],
               min_data_in_leaf=min_data_in_leaf[k],
               lambda_l1=reg_alpha[k],
               lambda_l2=reg_lambda[k],
               verbose=verbose,
               seed=seeds[0],
               n_jobs=n_jobs,)

test_preds = np.zeros(len(test_df))
dtrain = lgb.Dataset(train_df[take_best_features], y_df)
dtest = test_df[take_best_features].copy()
testlen = test_df.shape[0]

# predicting and ensembling several lgbm models, averaging over predictions for each fold and averaging over ranked prediction for each model
for i, seed in enumerate(seeds):
    print(f'Training Model with SEED : {seed}')
    evals['lgb'][i] = lgb.cv(lparams[i],
                             dtrain,
                             nfold=nfold[i],
                             stratified=stratified[i],
                             shuffle=shuffle[i],
                             num_boost_round=n_estimators,
                             early_stopping_rounds=early_stopping_rounds,
                             verbose_eval=verbose_eval,
                             return_cvbooster=True,
                             seed = seed,
                             show_stdv=True)
    print(f'SEED {i} Average fold  AUC {np.round(max(evals["lgb"][i]["auc-mean"]),5)}')
    test_preds += stats.rankdata(np.mean(evals['lgb'][i]['cvbooster'].predict(dtest, num_iteration=evals['lgb'][i]['cvbooster'].best_iteration), axis=0)) / (testlen * len(seeds))

submission = pd.DataFrame.from_dict({
    'encounter_id':test_id,
    'diabetes_mellitus':test_preds,
})
submission.to_csv('submission.csv', index=False)