<jupyter_start><jupyter_text>description 6/03/19 clinically guided aggregation modeling w/ 2 elix variablessklearn modeling using local methods of the median imputed training data using origional min/max clinically guided aggregation. note the preprocessing of data from 07.20-worst_case_model was performed in R (09.newagg2_preprocessing_med_impute.rmd). this eventually will be converted over to python, but for now works in r. preprocessing includes variable formatting (categorical to factor variables in r, train/test split, and median imputation).<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import os, sys
from pathlib import Path
import seaborn as sns
import numpy as np
import glob
from sklearn.metrics import confusion_matrix, classification_report, roc_curve, roc_auc_score, accuracy_score, auc, precision_recall_fscore_support, pairwise, f1_score, log_loss, make_scorer
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.externals.joblib import Memory
from sklearn.preprocessing import MinMaxScaler, OneHotEncoder, Imputer
from sklearn.model_selection import StratifiedKFold, GridSearchCV, RandomizedSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.utils import validation
from scipy.sparse import issparse
from scipy.spatial import distance
from sklearn import svm
from sklearn.impute import SimpleImputer

#importin xg boost and all needed otherstuff
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from xgboost import XGBClassifier #conda install -c conda-forge xgboost to install
##adding these, lets see if it helps with xgboost crash
os.environ['KMP_DUPLICATE_LIB_OK']='True'


#reducing warnings that are super common in my model
import warnings
from sklearn.exceptions import DataConversionWarning
warnings.simplefilter(action='ignore') #ignore all warnings
# warnings.simplefilter(action='ignore', category=FutureWarning)
# warnings.filterwarnings(action='ignore', category=DataConversionWarning)
# warnings.filterwarnings(action='ignore', category=DeprecationWarning)


memory = Memory(cachedir='/tmp', verbose=0)
<EMAIL> above any def fxn.

RANDOM_STATE = 15485867

%matplotlib inline
plt.style.use('ggplot')

from notebook.services.config import ConfigManager
cm = ConfigManager()
cm.update('livereveal', {
        'width': 1024,
        'height': 768,
        'scroll': True,
})


%load_ext autotime<jupyter_output>The autotime extension is already loaded. To reload it, use:
  %reload_ext autotime
time: 60.9 ms<jupyter_text>variables used to saving models<jupyter_code>date="11062019"#'03062019'
dataset="clinagg"
timewindowdays="72"
#cohort import

os.chdir('/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling') #use to change working directory
wd= os.getcwd() #'/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling'


final_pt_df2 = pd.read_csv(Path(wd + '/data/raw/csv/04042019_final_pt_df2_v.csv') , index_col=0) #only for patients with minimum vitals, 14478 icustay_id
patients= list(final_pt_df2['subject_id'].unique())
hadm_id= list(final_pt_df2['hadm_id'].unique())
icustay_id= list(final_pt_df2['icustay_id'].unique())
icustay_id= [int(x) for x in icustay_id]
#importing in all clinical_variable files
timewindowdays="72"
# date= '04042019'
date='11062019'
os.chdir(r'/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling/data/processed/')
preimp_train_df= pd.read_csv(Path(wd+'/data/processed/merged/{}_worst_df_train_preImp_{}.csv'.format(date,timewindowdays),  index_col=0))
preimp_test_df= pd.read_csv(Path(wd+'/data/processed/merged/{}_worst_df_test_preImp_{}.csv'.format(date,timewindowdays),  index_col=0))

def preprocessing(data):

    """
    1) rename columns
    2) standardize last 2 columns to be standardized
    3) convert categorical columns to proper format
    4) median impute
    """
    from sklearn.impute import SimpleImputer
        
    rename_dic={
    "('max', 'sodium')": "maxSodium" ,
    "('max', 'sodium')" : "maxSodium",
    "('min', 'sodium')" : "minSodium",
    "('max', 'calcium')" : "maxCalcium",
    "('min', 'calcium')" : "minCalcium",
    "('max', 'sodium')": "maxSodium",
    "('min', 'sodium')": "minSodium",
    "('max', 'wbc')": "maxWBC",
    "('min', 'wbc')": "minWBC",
    "bands": "ibands",
    "pco2": "ipco2"
        }
    preimp_df=preimp_df.rename(rename_dic, axis='columns')
    preimp_df.loc[preimp_df['first_admit_age']>90,"first_admit_age"]=90
    
    
    train_data=data.copy()
    weight_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).median()
    weight_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.25)#.between(train_data['col'].quantile(.25), df['col'].quantile(.75), inclusive=True)]
    weight_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.75)
    weight_iqr=weight_quant3-weight_quant1
    #print(weight_median,weight_quant3,weight_quant1, weight_iqr)

    age_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).median()
    age_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.25)
    age_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.75)
    age_iqr=age_quant3-age_quant1

    #converting to log scaled standardized data for age/weight
    train_data['weight']=train_data['weight'].apply(lambda x: (np.log(x+1)-weight_median)/weight_iqr)
    train_data['first_admit_age']=train_data['first_admit_age'].apply(lambda x: (np.log(x+1)-age_median)/age_iqr)
    
    ### onehot encoding categorical var
    cols_to_transform=['ethnicity', 'ibands', 'ipco2',
                       'any_vasoactive',"leukocyte","nitrite",
                       'pao2fio2ratio', 'vent_recieved',  "dobutamine",
                       "dopamine","epinephrine","norepinephrine",
                       "phenylephrine","rrt","vasopressin",'cancer_elix' ]
    train_data = pd.get_dummies(train_data, columns = cols_to_transform, drop_first=True)
    
    
    #binarizing and poping outcome for training data
    train_data.loc[train_data['final_bin']=="C_pos/A_full","final_bin"]=1
    train_data.loc[train_data['final_bin']=="C_neg/A_partial","final_bin"]=0
    train_data['final_bin']=pd.to_numeric(train_data['final_bin'])
    
    
    imp = SimpleImputer(missing_values=np.nan, strategy='median')
    train_data=pd.DataFrame(imp.fit_transform(train_data), columns=list(train_data))

    ## establishing training data and labels
    x_train= train_data.copy()
    z_icustay_id=x_train.pop('icustay_id')
    y_train= x_train.pop("final_bin").values
    
    return(x_train, y_train, z_icustay_id)

x_train, y_train, z_icustay_id= preprocessing(pd.merge(preimp_train_df, final_pt_df2[['icustay_id','final_bin']]))

x_train, y_train, z_icustay_id= preprocessing(pd.merge(preimp_train_df, final_pt_df2[['icustay_id','final_bin']]))
x_train
test_df
imp = SimpleImputer(missing_values=np.nan, strategy='median')
imp.fit(x_test_df)
pd.DataFrame(imp.fit_transform(x_test_df), columns=list(x_test_df)).head()
pd.DataFrame(test_df)
pd.DataFrame(imp.fit_transform(x_test_df), columns=list(x_test_df)).head()
imp = SimpleImputer(missing_values=np.nan, strategy='median')
imp.fit(preimp_train_df)<jupyter_output><empty_output><jupyter_text>r code to translate<jupyter_code># cols_to_transform=['ethnicity', 'ibands', 'ipco2',
#                    'any_vasoactive',"leukocyte","nitrite",
#                    'pao2fio2ratio', 'vent_recieved',  "dobutamine",
#                    "dopamine","epinephrine","norepinephrine",
#                    "phenylephrine","rrt","vasopressin",'cancer_elix' ]<jupyter_output><empty_output><jupyter_text>functions related to gower distance, preprocessing, finding train samples, etc...<jupyter_code>def continuous_filter(df):  
    global all_xy
    categorical=['gender',
     'ethnicity_black',
     'ethnicity_hispanic',
     'ethnicity_unknown/other',
     'ethnicity_white/nonhispanic',
     'ibands_absent',
     'any_vasoactive_True',
     'leukocyte_1',
     'nitrite_1',
     'pao2fio2ratio(200, 333]',
     'pao2fio2ratio_(333, 475]',
     'pao2fio2ratio_(475, 3000]',
     'vent_recieved_1',
     'vent_recieved_2',
     'dobutamine_True',
     'dopamine_True',
     'epinephrine_True',
     'norepinephrine_True',
     'phenylephrine_True',
     'rrt_True',
     'vasopressin_True',
     'ipco2_absent',
     "cancer_elix_True"]
    all_xy_label= list(all_xy)
    in_both= list(set(categorical)& set(all_xy))

    ##restricting all_xy to only continuous variables
    all_xy_cont = df.drop(in_both, axis=1)
    return(all_xy_cont)
#folder = '/Users/xuzhenxing/Documents/mimic_AKI_data/real_time_prediction/features/all/dropped/xy'
# folder = './xy'

def preprocessing(folder, time_interval, isnormalized=True):
    """Data preprocessing, Preprocessing  missing data with mean imputation; Normalize continous feature with MinMaxScaler;
    Normalize categorical feature with OneHotEncoder.

    Args:
        folder: dir path of source data;
        time_interval: interval of time, can be 24,48,72,96,120,144.
    Returns:
        x: features
        y: lables

    """

    all_xy = pd.read_csv(os.path.join(folder, 'all_{}hours_test_individualization_1thousand.csv'.format(time_interval)), index_col=0)
    # print (all_xy.shape)
    # print (all_xy.columns)

    medi = ['diuretics', 'nsaid', 'radio', 'angiotensin']
    pat = ['gender', 'age', 'ethnicity']
    # Total 9 comorbidity
    comm = ['congestive_heart_failure', 'peripheral_vascular', 'hypertension',
            'diabetes', 'liver_disease', 'mi', 'cad', 'cirrhosis', 'jaundice']

    # Total 8 chartevents
    chart = ['DiasBP_min', 'DiasBP_max', 'DiasBP_first', 'DiasBP_last', 'DiasBP_slope', 'DiasBP_avg',
             'Glucose_min', 'Glucose_max', 'Glucose_first', 'Glucose_last', 'Glucose_slope', 'Glucose_avg',
             'HeartRate_min', 'HeartRate_max', 'HeartRate_first', 'HeartRate_last', 'HeartRate_slope', 'HeartRate_avg',
             'MeanBP_min', 'MeanBP_max', 'MeanBP_first', 'MeanBP_last', 'MeanBP_slope', 'MeanBP_avg',
             'RespRate_min', 'RespRate_max', 'RespRate_first', 'RespRate_last', 'RespRate_slope', 'RespRate_avg',
             'SpO2_min', 'SpO2_max', 'SpO2_first', 'SpO2_last', 'SpO2_slope', 'SpO2_avg',
             'SysBP_min', 'SysBP_max', 'SysBP_first', 'SysBP_last', 'SysBP_slope', 'SysBP_avg',
             'Temp_min', 'Temp_max', 'Temp_first', 'Temp_last', 'Temp_slope', 'Temp_avg']

    # Total 12 labvents
    lab = ['BICARBONATE_first', 'BICARBONATE_last', 'BICARBONATE_min', 'BICARBONATE_max', 'BICARBONATE_avg',
           'BICARBONATE_slope', 'BICARBONATE_count',
           'BUN_first', 'BUN_last', 'BUN_min', 'BUN_max', 'BUN_avg', 'BUN_slope', 'BUN_count',
           'CHLORIDE_first', 'CHLORIDE_last', 'CHLORIDE_min', 'CHLORIDE_max', 'CHLORIDE_avg', 'CHLORIDE_slope',
           'CHLORIDE_count',
           'CREATININE_first', 'CREATININE_last', 'CREATININE_min', 'CREATININE_max', 'CREATININE_avg',
           'CREATININE_slope', 'CREATININE_count',
           'HEMOGLOBIN_first', 'HEMOGLOBIN_last', 'HEMOGLOBIN_min', 'HEMOGLOBIN_max', 'HEMOGLOBIN_avg',
           'HEMOGLOBIN_slope', 'HEMOGLOBIN_count',
           'INR_first', 'INR_last', 'INR_min', 'INR_max', 'INR_avg', 'INR_count',
           'PLATELET_first', 'PLATELET_last', 'PLATELET_min', 'PLATELET_max', 'PLATELET_avg', 'PLATELET_slope',
           'PLATELET_count',
           'POTASSIUM_first', 'POTASSIUM_last', 'POTASSIUM_min', 'POTASSIUM_max', 'POTASSIUM_avg', 'POTASSIUM_slope',
           'POTASSIUM_count',
           'PT_first', 'PT_last', 'PT_min', 'PT_max', 'PT_avg', 'PT_count',
           'PTT_first', 'PTT_last', 'PTT_min', 'PTT_max', 'PTT_avg', 'PTT_count',
           'WBC_first', 'WBC_last', 'WBC_min', 'WBC_max', 'WBC_avg', 'WBC_slope', 'WBC_count',
           'CALCIUM_first', 'CALCIUM_last', 'CALCIUM_min', 'CALCIUM_max', 'CALCIUM_avg', 'CALCIUM_count'
           ]

    if time_interval != 24:  # The 24h data lack of the feature 'CALCIUM_slope'
        lab.append('CALCIUM_slope')
    subset = medi + pat + comm + ['avg_urine'] + ['egfr_min'] + ['label'] # note that ['avg_urine'] + ['egfr_min'] is important, ignoring if they are empty.

    all_xy = all_xy.dropna(subset=subset)

    # print ('after dropping nan in the catergorical variables, the shape is {}'.format(all_xy.shape))

    all_conti_x = all_xy[chart + lab + ['avg_urine'] + ['egfr_min'] + ['age']]
    # print (all_conti_x.shape)
    # print (all_conti_x)
    all_categ_x = all_xy[['gender'] + ['ethnicity'] + medi + comm]
    # print (all_categ_x.shape)
    # print (all_categ_x)

    # Using mean imputer after drop the nan data in medication, patient demographic data, avg_ureine, egfr_min and label
    imp = Imputer(strategy='mean', axis=0)
    all_conti_x_fitted = imp.fit_transform(all_conti_x)

    def normalize(all_conti_x_fitted, all_categ_x):
        # using the MinMaxScaler to normalization the all_x
        min_max_scaler = MinMaxScaler()
        all_conti_x_fitted = min_max_scaler.fit_transform(all_conti_x_fitted)
        # print (all_conti_x_fitted.shape, all_conti_x_fitted)
        # all_conti_x = DataFrame(all_conti_x_fitted, columns=all_conti_x.columns)
        # print (all_conti_x.shape)

        onehot_enc = OneHotEncoder(sparse=False)  # dense format
        all_categ_x_fitted = onehot_enc.fit_transform(all_categ_x)
        # print (all_categ_x_fitted.shape, all_categ_x_fitted)
        return all_conti_x_fitted, all_categ_x_fitted

    if isnormalized:
        all_conti_x_fitted, all_categ_x_fitted = normalize(all_conti_x_fitted, all_categ_x)

    x = np.hstack((all_conti_x_fitted, all_categ_x_fitted))
    # y = all_xy['label']
    # x = np.array(x)
    # y = np.array(y)
    # print (x.shape, y.shape)
    # return x, y
    y = all_xy['label']
    z_icustay_id = y.index
    x = np.array(x)
    y = np.array(y)
    z_icustay_id = np.array(z_icustay_id)

    print (x.shape, y.shape)
    return x, y, z_icustay_id, all_xy

def perf_model(pipe, param_grid, name, X_train, X_test,
               y_train, y_test, scoring, verbose=0):
    gs = GridSearchCV(estimator=pipe, param_grid=param_grid, scoring=scoring, cv=5, n_jobs=-1, verbose=verbose)
    gs.fit(X_train, y_train)

    y_train_pred = gs.predict(X_train)
    y_test_pred = gs.predict(X_test)

    acc_train = accuracy_score(y_true=y_train, y_pred=y_train_pred)
    acc_test = accuracy_score(y_true=y_test, y_pred=y_test_pred)

    fpr, tpr, _ = roc_curve(y_train, gs.predict_proba(X_train)[:, 1])
    auc_train = auc(fpr, tpr)

    fpr, tpr, _ = roc_curve(y_test, gs.predict_proba(X_test)[:, 1])
    auc_test = auc(fpr, tpr)

    confmat_train = confusion_matrix(y_true=y_train, y_pred=y_train_pred)
    confmat_test = confusion_matrix(y_true=y_test, y_pred=y_test_pred)

    print (' best parameter: ', gs.best_params_)
    print (' training acc:%.2f auc:%.2f ' % (acc_train, auc_train))
    print (' testing acc:%.2f auc:%.2f ' % (acc_test, auc_test))

    print (' train confusion matrix:\n', confmat_train)
    print (' testing confusion matrix:\n', confmat_test)
    print (' classification report:\n', classification_report(y_test, y_test_pred))

    train_report = np.array(precision_recall_fscore_support(y_train, y_train_pred))
    train_class1_report = train_report[:, 1]
    train_metrics = list(train_class1_report[:-1])
    train_metrics.extend([acc_train, auc_train])
    print ('training metrics: precision, recall, f1-score, acc, auc')
    print (train_metrics)

    test_report = np.array(precision_recall_fscore_support(y_test, y_test_pred))
    test_class1_report = test_report[:, 1]
    test_metrics = list(test_class1_report[:-1])
    test_metrics.extend([acc_test, auc_test])
    print ('test metrics: precision, recall, f1-score, acc, auc')
    print (test_metrics)

    return train_metrics, test_metrics
    """
    plt.xlabel("False Positive Rate")
    plt.ylabel("True Positive Rate (recall)")

    plt.plot(fpr, tpr, label="acc:%f auc:%f" % (acc_test, auc_test))
    plt.legend(loc="best")
    plt.show()
    plt.close()

    precision, recall, _ = precision_recall_curve(y_train, gs.predict_proba(X_train)[:,1])
    average_precision = average_precision_score(y_test, gs.predict_proba(X_test)[:,1])
    plt.xlabel("precision")
    plt.ylabel("recall")
    plt.step(precision, recall, where='post', label='AP={0:0.2f}'.format(average_precision))
    plt.legend(loc="best")
    plt.show()
    plt.close()
    """

def try_dbdt(X_train, X_test, y_train, y_test, scoring):
    gbm = GradientBoostingClassifier(learning_rate=0.05, n_estimators=120, min_samples_leaf=60,
                                     max_features=9, subsample=0.7, random_state=10)

    param_grid = {'max_depth': list(range(3, 14, 2)), 'min_samples_split': list(range(100, 801, 200))}
    train_metrics, test_metrics = perf_model(gbm, param_grid, 'GBDT', X_train, X_test, y_train, y_test, scoring, 0)
    return train_metrics, test_metrics
#issue im having is that 

def try_models_cross(X_train, X_test, y_train, y_test, scoring):#  select data cross 5 Fold
    # X_train, X_test, y_train, y_test = train_test_split(X, Y, train_size=0.7, stratify=Y, random_state=RANDOM_STATE)
    # """
    # print ('\n\nLinear Logistic Regression with L1 Penalty')
    # lgr_l1_train_metrics, lgr_l1_test_metrics = try_lgr_l1(X_train, X_test, y_train, y_test, scoring)
    #
    # print ('\n\nLinear Logistic Regression with L2 Penalty')
    # lgr_l2_train_metrics, lgr_l2_test_metrics = try_lgr_l2(X_train, X_test, y_train, y_test, scoring)
    #
    # print ('\n\nStochastic Gradient Descent')
    # Elastic_train_metrics, Elastic_test_metrics = try_sgd(X_train, X_test, y_train, y_test, scoring)
    #
    # print ('\n\nRandom Forest')
    # rf_train_metrics, rf_test_metrics = try_rf(X_train, X_test, y_train, y_test, scoring)
    # #
    print ('\n\nGradient Boosting Decision tree')
    xgboost_train_metrics, xgboost_test_metrics = try_dbdt(X_train, X_test, y_train, y_test, scoring)<jupyter_output><empty_output><jupyter_text>importing and formatting dataset<jupyter_code>os.chdir('/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling') #use to change working directory
wd= os.getcwd() #'/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling'

final_pt_df2 = pd.read_csv(Path(wd + '/data/raw/csv/04042019_final_pt_df2_v.csv') , index_col=0) #only for patients with minimum vitals
patients= list(final_pt_df2['subject_id'].unique())
hadm_id= list(final_pt_df2['hadm_id'].unique())
icustay_id= list(final_pt_df2['icustay_id'].unique())
icustay_id= [int(x) for x in icustay_id]
train_data= pd.read_csv("/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling/models/imputation/{}_median_imputed_train_{}.csv".format(date,timewindowdays)) #two class training data
test_data= pd.read_csv("/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling/models/imputation/{}_median_imputed_test_{}.csv".format(date,timewindowdays)) #two class training data

#/Users/geickelb1/Documents/GitHub/mimiciii-antibiotics-modeling/models/imputation/04042019_newagg2_median_imputed_test.csv
final_pt_df2['icustay_id'].nunique()<jupyter_output><empty_output><jupyter_text>light data reformatting for model most data are already converted to median type zscores, however weight and admit age still need to be converted.<jupyter_code>def preprocessing(data):
    train_data=data.copy()
    weight_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).median()
    weight_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.25)#.between(train_data['col'].quantile(.25), df['col'].quantile(.75), inclusive=True)]
    weight_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.75)
    weight_iqr=weight_quant3-weight_quant1
    #print(weight_median,weight_quant3,weight_quant1, weight_iqr)

    age_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).median()
    age_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.25)
    age_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.75)
    age_iqr=age_quant3-age_quant1

    #converting to log scaled standardized data for age/weight
    train_data['weight']=train_data['weight'].apply(lambda x: (np.log(x+1)-weight_median)/weight_iqr)
    train_data['first_admit_age']=train_data['first_admit_age'].apply(lambda x: (np.log(x+1)-age_median)/age_iqr)
    
    ### onehot encoding categorical var
    cols_to_transform=['ethnicity', 'ibands', 'ipco2',
                       'any_vasoactive',"leukocyte","nitrite",
                       'pao2fio2ratio', 'vent_recieved',  "dobutamine",
                       "dopamine","epinephrine","norepinephrine",
                       "phenylephrine","rrt","vasopressin",'cancer_elix' ]
    train_data = pd.get_dummies(train_data, columns = cols_to_transform, drop_first=True)
    
    
    #binarizing and poping outcome for training data
    train_data.loc[train_data['final_bin']=="C_pos/A_full","final_bin"]=1
    train_data.loc[train_data['final_bin']=="C_neg/A_partial","final_bin"]=0
    train_data['final_bin']=pd.to_numeric(train_data['final_bin'])
    
    ## establishing training data and labels
    x_train= train_data.copy()
    z_icustay_id=x_train.pop('icustay_id')
    y_train= x_train.pop("final_bin").values
    
    return(x_train, y_train, z_icustay_id)
def xy_preprocessing(data):
    train_data=data.copy()
    weight_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).median()
    weight_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.25)#.between(train_data['col'].quantile(.25), df['col'].quantile(.75), inclusive=True)]
    weight_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","weight"]+1).quantile(0.75)
    weight_iqr=weight_quant3-weight_quant1
    #print(weight_median,weight_quant3,weight_quant1, weight_iqr)

    age_median=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).median()
    age_quant1=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.25)
    age_quant3=np.log(train_data.loc[train_data['final_bin']=="C_neg/A_partial","first_admit_age"]+1).quantile(0.75)
    age_iqr=age_quant3-age_quant1

    #converting to log scaled standardized data for age/weight
    train_data['weight']=train_data['weight'].apply(lambda x: (np.log(x+1)-weight_median)/weight_iqr)
    train_data['first_admit_age']=train_data['first_admit_age'].apply(lambda x: (np.log(x+1)-age_median)/age_iqr)
    
    ### onehot encoding categorical var
    cols_to_transform=['ethnicity', 'ibands', 'ipco2', 'any_vasoactive',"leukocyte","nitrite", 'pao2fio2ratio', 'vent_recieved',  "dobutamine","dopamine","epinephrine","norepinephrine","phenylephrine","rrt","vasopressin",'cancer_elix' ]
    train_data = pd.get_dummies(train_data, columns = cols_to_transform, drop_first=True)
    
    
    #binarizing and poping outcome for training data
    train_data.loc[train_data['final_bin']=="C_pos/A_full","final_bin"]=1
    train_data.loc[train_data['final_bin']=="C_neg/A_partial","final_bin"]=0
    train_data['final_bin']=pd.to_numeric(train_data['final_bin'])
    
    ## establishing training data and labels

    all_xy= train_data.copy()
    all_xy=all_xy.set_index('icustay_id').rename(columns={'final_bin':"label"})
    
    return(all_xy)
x_train, y_train, z_icustay_id = preprocessing(train_data)
x_test, y_test, z_icustay_id_test= preprocessing(test_data)
#for local modeling
all_xy=xy_preprocessing(train_data)
all_xy_test=xy_preprocessing(test_data)
list(all_xy)<jupyter_output><empty_output><jupyter_text>optional qc<jupyter_code>#x_train.iloc[1:5, 25:45]
#x_train.iloc[1:5, 35:65]
#x_train.iloc[1:5, 10:30]<jupyter_output><empty_output><jupyter_text>looking at correlation of all variables<jupyter_code>corr = x_train.corr().abs()

plt.figure(figsize=(25, 20))
ax = sns.heatmap(
    corr, 
    vmin=-1, vmax=1, center=0,
    cmap=sns.diverging_palette(20, 220, n=200),
    square=True
)
ax.set_xticklabels(
    ax.get_xticklabels(),
    rotation=45,
    horizontalalignment='right'
);
sol = (corr.where(np.triu(np.ones(corr.shape), k=1).astype(np.bool)).stack().sort_values(ascending=False))
cor_df=pd.DataFrame(sol)#.sort_values(kind="quicksort") #[-10:0])
cor_df=cor_df.reset_index()
cor_df=cor_df.rename(columns={'level_0': 'corx', 'level_1': 'cory', 0:'corr'})
cor_df2=cor_df[(cor_df['corx']!=cor_df['cory']) & (cor_df['corr']>0.7)].sort_values('corr', ascending=False)
cor_df2.head()<jupyter_output><empty_output><jupyter_text>DROPING one of the 2 columns with correlation >0.7 	corx	cory	corr0	ipco2_absent	pao2fio2Ratio_(475, 3000]	0.8724181	maxWBC	minWBC	0.8023732	bun	creatinine	0.7208613	maxSodium	minSodium	0.704233<jupyter_code>x_train.drop(columns=list(cor_df2['corx']), inplace=True, errors='raise')
x_test.drop(columns=list(cor_df2['corx']), inplace=True, errors='raise')
all_xy.drop(columns=list(cor_df2['corx']), inplace=True, errors='raise')
all_xy_test.drop(columns=list(cor_df2['corx']), inplace=True, errors='raise')<jupyter_output><empty_output><jupyter_text>formatting x and y for modleing<jupyter_code>#x=np.array(x_train.iloc[:,[1,2,3,4,5,6,7,8,9,38,39,40,41]]).copy() #copy of x_train
#x=np.array(x_train.iloc[:,[1,2,3,4,5,6,7,8,9]]).copy() #copy of x_train
#x=np.array(x_train.iloc[:,38:])
x=np.array(x_train.copy())

#x=np.array(train_data.iloc[:,[1,2,3,4]]).copy() #copy of x_train
#train_data.iloc[:,[1,2,3,4,5]] ###drastically reducing my dataframe size to test algorithm
y=y_train.copy() #copy of y_train

##all_xy: train data with finalbin:label and index=icustay_id
#all_xy=train_data.copy().set_index("icustay_id").rename(columns={'final_bin':"label"}) #

time_interval=4
print(len(x_train),len(x_test))<jupyter_output><empty_output><jupyter_text>Modelbuilding* step1) hypertune xgb on 5fold cv.* step2) test entire trainset and predict trainset.* step3) run hypertuned model on 5fold cv with lr and get overall metrics.* step4) local model testing  step1) XGB hypertuning<jupyter_code>def evaluate(model, test_features, test_labels):
    from sklearn.metrics import log_loss
    
    y_hat = model.predict(test_features)
    errors = abs(y_hat - test_labels)
    mape = 100 * np.mean(errors / test_labels)
    accuracy = 100 - mape
    auc=roc_auc_score(test_labels, y_hat)
    loss= log_loss(test_labels, y_hat)
    
    print ('the AUC is: {:0.2f}'.format(auc))
    print ('the logloss is: {:0.2f}'.format(loss))
    print(confusion_matrix(test_labels, y_hat))
    
    return loss


# base_model = RandomForestClassifier(n_estimators = 10, random_state = 42)
# base_model.fit(x, y)
# base_auc = evaluate(base_model, x, y)

# best_random = rf_random.best_estimator_
# random_auc = evaluate(best_random, x, y)

# print('Improvement of {:0.2f}%.'.format( 100 * (random_auc - base_auc) / base_auc))
def hypertuning_fxn(X, y, nfolds, model , param_grid, base_model, scoring="neg_log_loss", gridsearch=True, n_iter=20): 
    if gridsearch==True:
        grid_search = GridSearchCV(estimator= model,
                                         param_grid=param_grid,
                                         cv=nfolds,
                                         scoring=scoring,
                                         return_train_score=True,
                                         n_jobs = -1)
    else:
        grid_search = RandomizedSearchCV(estimator= model,
                                         param_distributions= param_grid,
                                         n_iter=n_iter,
                                         cv=nfolds,
                                         scoring=scoring,
                                         return_train_score=True,
                                         n_jobs = -1)
    grid_search.fit(X, y)    
    
    print("Grid scores on development set:")
    means = grid_search.cv_results_['mean_test_score']
    stds = grid_search.cv_results_['std_test_score']
    
    for mean, std, params in zip(means, stds, grid_search.cv_results_['params']):
        print("%0.3f (+/-%0.03f) for %r"
              % (mean, std * 2, params))
        
    #grid_search.best_params_
    print(grid_search.best_score_)
    print("\n")
    print(grid_search.best_params_)
    
    print('\n base model:')
    base_model = base_model#(random_state = 42)
    base_model.fit(x, y)
    base_auc = evaluate(base_model, x, y)
    
    print('\n hypertuned model:')
    best_random = grid_search.best_estimator_
    random_auc = evaluate(best_random, x, y)

    print('logloss change of {:0.2f}%. after hypertuning on training set (may be overfit)'.format( 100 * (random_auc - base_auc) / base_auc))
    
    print(grid_search.best_estimator_)
    
    return(grid_search)
###xgboost
model= XGBClassifier(n_estimators=100, min_child_weight=2, #changed: GridSearchCV ->RandomizedSearchCV
                                              gamma=0, subsample=0.8, colsample_bytree=0.8,
                                              objective='binary:logistic', n_jobs=-1, seed=27)
scale_pos_weight = [0.1, 1, 5, 10]
max_depth = [1, 2, 3, 4, 5]
learning_rate=[0.01, 0.1, 0.5, 1]
param_grid = {'scale_pos_weight': scale_pos_weight, 'max_depth' : max_depth, "learning_rate":learning_rate}

base_model=XGBClassifier(random_state = 42)
xgboost_hyper=hypertuning_fxn(x, y, nfolds=5, model=model , param_grid=param_grid, base_model= base_model, scoring="neg_log_loss", n_iter=20, gridsearch=True)
###rf
# Number of trees in random forest
#n_estimators = [100, 1000]#[int(x) for x in np.linspace(start = 10, stop = 1000, num = 5)]
# Number of features to consider at every split
max_features = [3,'auto', 10 ]
# Maximum number of levels in tree
max_depth = [5,10, 25]#[int(x) for x in np.linspace(5, 110, num = 5)]
#max_depth.append(None)
# Minimum number of samples required to split a node
min_samples_split = [2, 5, 10]
# Minimum number of samples required at each leaf node
min_samples_leaf = [2, 5, 10]
# Method of selecting samples for training each tree
#bootstrap = [True, False]

#class_weight is either a dictionary of each class to a uniform weight for that class (e.g., {1:.9, 2:.5, 3:.01}), or is a string telling sklearn how to automatically determine this dictionary.
class_weight= [None,{0:1, 1:4}, {0:(1/np.bincount(y))[0], 1:(1/np.bincount(y))[1]}]


param_grid = {#'n_estimators': n_estimators,
               'max_features': max_features,
               'max_depth': max_depth,
               'min_samples_split': min_samples_split,
               'min_samples_leaf': min_samples_leaf,
               'class_weight': class_weight}


model= RandomForestClassifier(criterion='entropy')
base_model=RandomForestClassifier(random_state = 42, criterion='entropy')

rf_hyper=hypertuning_fxn(x, y, nfolds=10, model=model , param_grid=param_grid, base_model= base_model, scoring="neg_log_loss",n_iter = 30, gridsearch=True)

class_weight<jupyter_output><empty_output><jupyter_text>Hypertune SVC<jupyter_code>model= svm.SVC(probability=True)
kernel = ['linear', 'rbf', 'poly']
gamma = [0.1, 1, 'auto'] #Kernel coefficient for ‘rbf’, ‘poly’ and ‘sigmoid’. default=’auto’ uses 1 / n_features
C = [0.1, 1, 10, 100] #Penalty parameter C of the error term.
degree = [0, 1, 2]
class_weight=['balanced', None]

param_grid = {'kernel': kernel,
              'gamma': gamma,
              'C': C,
              'degree': degree,
              'class_weight':class_weight}

base_model=svm.SVC(probability=True)

svc_hyper=hypertuning_fxn(x, y, nfolds=4, model=model , param_grid=param_grid, base_model= base_model, scoring="neg_log_loss", n_iter=10, gridsearch=False)
###knn
from sklearn.neighbors import KNeighborsClassifier
model= KNeighborsClassifier()

n_neighbors = [3,4,5, 8, 10, 25]
weights=['uniform']
p=[1,2] #1= mmanhattan, 2= euclidian


param_grid = {'n_neighbors': n_neighbors,
              'weights': weights,
              'p': p}

base_model=KNeighborsClassifier()

knn_hyper=hypertuning_fxn(x, y, nfolds=10, model=model , param_grid=param_grid, base_model= base_model, scoring="neg_log_loss", n_iter=40, gridsearch=True)
knn_hyper.best_params_
evaluate(knn_hyper, x, y)
evaluate(base_model, x, y)
base_model.get_params
#KNeighborsClassifier(RANDOM_STATE=42)<jupyter_output><empty_output><jupyter_text>Hypertuned Model Initialization<jupyter_code># xgboost = XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
#        colsample_bytree=0.8, gamma=0, learning_rate=0.1, max_delta_step=0,
#        max_depth=4, min_child_weight=2, missing=None, n_estimators=100,
#        n_jobs=-1, nthread=None, objective='binary:logistic',
#        random_state=0, reg_alpha=0, reg_lambda=1, scale_pos_weight=1,
#        seed=27, silent=True, subsample=0.8)
# #xgboost.fit(x, y)

# logreg = LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=10, fit_intercept=True,
#                         intercept_scaling=1, class_weight='balanced', random_state=None)

# #logreg.fit(x, y)

# rf= RandomForestClassifier(bootstrap=False, class_weight={0: 1, 1: 4},
#             criterion='entropy', max_depth=10, max_features='auto',
#             max_leaf_nodes=None, min_impurity_decrease=0.0,
#             min_impurity_split=None, min_samples_leaf=2,
#             min_samples_split=2, min_weight_fraction_leaf=0.0,
#             n_estimators=600, n_jobs=None, oob_score=False,
#             random_state=None, verbose=0, warm_start=False)
# #rf.fit(x,y)

# # from sklearn.naive_bayes import GaussianNB
# # gnb =GaussianNB()
# # nb_y_pred = gnb.fit(x, y)

# svc= svm.SVC(C=100, cache_size=200, class_weight='balanced', coef0=0.0,
#           decision_function_shape='ovr', degree=0, gamma=1, kernel='linear',
#           max_iter=-1, probability=True, random_state=None, shrinking=True,
#           tol=0.001, verbose=False)
# #svc.fit(x, y)

xgboost= xgboost_hyper.best_estimator_

logreg = LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=10, fit_intercept=True,
                        intercept_scaling=1, class_weight='balanced', random_state=None)

rf= rf_hyper.best_estimator_

svc=svc_hypter.best_estimator_

def reset_model(model_name, hardcode=True):
    if hardcode==True:
        if model_name== 'xgboost':
            model = XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
               colsample_bytree=0.8, gamma=0, learning_rate=0.1, max_delta_step=0,
               max_depth=4, min_child_weight=2, missing=None, n_estimators=100,
               n_jobs=-1, nthread=None, objective='binary:logistic',
               random_state=0, reg_alpha=0, reg_lambda=1, scale_pos_weight=1,
               seed=27, silent=True, subsample=0.8)
            

        elif model_name== 'logreg':
            model = LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=10, fit_intercept=True,
                                    intercept_scaling=1, class_weight='balanced', random_state=None)

        elif model_name== 'rf':
            model = RandomForestClassifier(bootstrap=False, class_weight={0: 1, 1: 4},
                criterion='entropy', max_depth=10, max_features='auto',
                max_leaf_nodes=None, min_impurity_decrease=0.0,
                min_impurity_split=None, min_samples_leaf=2,
                min_samples_split=2, min_weight_fraction_leaf=0.0,
                n_estimators=600, n_jobs=None, oob_score=False,
                random_state=None, verbose=0, warm_start=False)

        elif model_name== 'svc':
            model = svm.SVC(C=100, cache_size=200, class_weight='balanced', coef0=0.0,
                  decision_function_shape='ovr', degree=0, gamma=1, kernel='linear',
                  max_iter=-1, probability=True, random_state=None, shrinking=True,
                  tol=0.001, verbose=False)
            
        elif model_name== 'knn':
            model = KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
                   metric_params=None, n_jobs=None, n_neighbors=25, p=1,
                   weights='uniform')
    
    else:
            if model_name== 'xgboost':
                model = xgboost_hyper.best_estimator_

            elif model_name== 'logreg':
                model = LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=10, fit_intercept=True,
                                        intercept_scaling=1, class_weight='balanced', random_state=None)

            elif model_name== 'rf':
                model = rf_hyper.best_estimator_

            elif model_name== 'svc':
                model = svc_hyper.best_estimator_
                
            elif model_name== 'knn':
                model = knn_hyper.best_estimator_
    return(model)
def get_auc_score(model,train_index, x=x,y=y):
    y_pred_proba = model.predict_proba(x[train_index])[:, 1] 
    roc_score=roc_auc_score(y[train_index], y_pred_proba)
    return(roc_score)<jupyter_output><empty_output><jupyter_text>youden index and plotting functions<jupyter_code>def optimal_youden_index(fpr, tpr, thresholds, tp90=True):
    """
    inputs fpr, tpr, thresholds from metrics.roc(),
    outputs the clasification threshold, roc dataframe, and the index of roc dataframe for optimal youden index
    """
    #making dataframe out of the thresholds
    roc_df= pd.DataFrame({"thresholds": thresholds,"fpr":fpr, "tpr": tpr})
    roc_df.iloc[0,0] =1
    roc_df['yuden']= roc_df['tpr']-roc_df['fpr']
    
    if tp90==True:
        idx= roc_df[roc_df['tpr']>=0.9]['yuden'].idxmax() #changed this so now finds optimial yuden threshold but tp>=90%
    else:
        idx=roc_df['yuden'].idxmax() #MAX INDEX
    
    youden_threshold=roc_df.iloc[idx,0] #threshold for max youden
    return(youden_threshold, roc_df, idx)
    
def plot_roc(fpr, tpr, roc_auc, roc_df, idx, save=False,model_name=None, folder_name=None, file_name=None):
    plt.title('ROC with optimal Youden Index')
    plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
    plt.legend(loc = 'lower right')
    plt.plot([0, 1], [0, 1],'r--')
    
    #finding the point on the line given threshold 0.5 (finding closest row in roc_df)
    og_idx=roc_df.iloc[(roc_df['thresholds']-0.5).abs().argsort()[:1]].index[0]
    plt.plot(roc_df.iloc[og_idx,1], roc_df.iloc[og_idx,2],marker='x', markersize=5, color="g")
    plt.annotate(s="P(>=0.5)",xy=(roc_df.iloc[og_idx,1]+0.02, roc_df.iloc[og_idx,2]-0.04),color='g') #textcoords
    
    
    plt.plot(roc_df.iloc[idx,1], roc_df.iloc[idx,2],marker='o', markersize=5, color="r") ##
    plt.annotate(s="M_Youden",xy=(roc_df.iloc[idx,1]+0.02, roc_df.iloc[idx,2]-0.04),color='r' ) #textcoords
    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.ylabel('True Positive Rate')
    plt.xlabel('False Positive Rate')
    
    if save==True:
        if folder_name != None:
            address = 'figures/{}/'.format(folder_name)
        else:
            address = 'figures/'
        if not os.path.exists(address):
            os.makedirs(address)
        plt.savefig(address+"/{}_{}.png".format(model_name,file_name),bbox_inches='tight')
    else: pass
    
    plt.show()
    
def plot_table_as_fig(table_in, col_labels, row_labels, save=False,model_name=None,folder_name=None, file_name=None,figsize=(6,1)):
    
    fig = plt.figure(figsize=figsize)
    table = plt.table(cellText = table_in, 
                  colLabels = col_labels,
                  rowLabels = row_labels,
                  loc='best')
    plt.axis("tight")
    plt.axis('off')
    if save==True:
        if folder_name != None:
            address = 'figures/{}/'.format(folder_name)
        else:
            address = 'figures/'
        if not os.path.exists(address):
            os.makedirs(address)
        plt.savefig(address+"/{}_{}.png".format(model_name,file_name),bbox_inches='tight')
    else: pass
    
    plt.show()
    
def classifier_eval(model, x=x, y=y, proba_input=False,pos_label=1, print_default=True, save=False,model_name=None, folder_name=None):
    import sklearn.metrics as metrics
    from sklearn.metrics import precision_score, roc_auc_score, f1_score, recall_score

    """
    catchall classification evaluation function. will print/save the following:
    
    print/save the following:
        ROC curve marked with threshold for optimal youden (maximizing tpr+fpr with constraint that tpr>0.9)

        using 0.5 threshold:
            confusion matrix
            classification report
            npv
            accuracy


        using optimal youden (maximizing tpr+fpr with constraint that tpr>0.9):
            confusion matrix
            classification report
            npv
            accuracy
    
    output: 
        outputs modelname, auc, precision, recall, f1, and npv to a dictionary. 
    
    notes:
    youden's J statistic:
    J= sensitivity + specificity -1
    (truepos/ truepos+falseneg) + (true neg/ trueneg + falsepos) -1
    
    """
    if save==True: #making folder if one doesn't exist
        if folder_name != None:
            address = '../figures/{}/'.format(folder_name)
        else:
            address = 'train/'
        if not os.path.exists(address):
            os.makedirs(address)
    
    if proba_input==True:  #incorporating classifier_eval2() functionality into this (ie allowing user to input a y_proba instead of a model)
        y_proba= model
        y_pred=[1 if y >= 0.5 else 0 for y in y_proba]
    
    else:
        model_name=type(model).__name__

        y_pred = model.predict(x)
        y_proba = model.predict_proba(x)[:,1]
        
    fpr, tpr, thresholds = metrics.roc_curve(y, y_proba, pos_label=pos_label)
    roc_auc = metrics.auc(fpr, tpr)
    
    #gathering the optimal youden_index and df of tpr/fpr for auc and index of that optimal youden. idx is needed in the roc
    youden_threshold, roc_df, idx= optimal_youden_index(fpr, tpr, thresholds,tp90=True)

    #plotting roc
    plot_roc(fpr, tpr, roc_auc, roc_df, idx, save=save, model_name=model_name,folder_name=folder_name, file_name='roc')
    plt.show(), plt.close()
    
    #printing npv, recall, precision, accuracy
    npv=confusion_matrix(y, y_pred)[0,0]/sum(np.array(y_pred)==0)*100
    prec= precision_score(y_true=y, y_pred= y_pred, pos_label=pos_label)
    recall= recall_score(y_true=y, y_pred= y_pred, pos_label=pos_label)
    f1= f1_score(y_true=y, y_pred= y_pred, pos_label=pos_label)
    
    if print_default==True: ###can opt to not print the 0.5 classification threshold classification report/conf matrix
        #plotting confusion matrixs
        print("\n******* Using 0.5 Classification Threshold *******\n")
        print(confusion_matrix(y, y_pred))
        print ('the Accuracy is: {:01.2f}'.format(accuracy_score(y, y_pred)))
        print ("npv: {:01.2f}".format(npv))
        print ('the classification_report:\n', classification_report(y,y_pred))
    else:
        pass
    
    #### YOUDEN ADJUSTMENT #####

    print("\n******* Using Optimal Youden Classification Threshold *******\n")
    print("\nthe Youden optimal index is : {:01.2f}".format(youden_threshold))

    y_pred_youden = [1 if y >= youden_threshold else 0 for y in y_proba]

    npv_y=confusion_matrix(y, y_pred_youden)[0,0]/sum(np.array(y_pred_youden)==0)*100
    prec_y= precision_score(y_true=y, y_pred= y_pred_youden, pos_label=pos_label)*100
    recall_y= recall_score(y_true=y, y_pred= y_pred_youden, pos_label=pos_label)*100
    f1_y= f1_score(y_true=y, y_pred= y_pred_youden, pos_label=pos_label)*100
    auc_y=roc_auc_score(y_true=y, y_score= y_proba)*100
    
    ##plotting and saving confusion matrix
    confusion_youden=confusion_matrix(y, y_pred_youden)
    
    plot_table_as_fig(confusion_youden,
                  col_labels=['predicted_neg','predicted_pos'],
                  row_labels=['true_neg',"true_pos"],
                  save=save,
                  figsize=(6,1),
                  model_name=model_name,
                  folder_name=folder_name,
                  file_name='y_confusion')
    plt.show(), plt.close()
   
    #print(confusion_matrix(y, y_pred_youden))
    print("the Accuracy is: {:01.2f}".format(accuracy_score(y, y_pred_youden)))
    print("npv: {:01.2f}".format(npv_y))
    
    ###formatting classification report to be compatable with matplotlib table
    report_youden=classification_report(y,y_pred_youden,output_dict=True)  
    report_youden = pd.DataFrame.from_dict(report_youden).transpose()[['precision','recall','f1-score','support']]
    report_youden = np.round(report_youden,2)
    
    ##plotting and saving classification report
    plot_table_as_fig(table_in=np.array(report_youden),#classification_report(y, xgboost.predict(x))),
                      col_labels=['precision','recall','f1-score','support'],
                      row_labels=['neg',"pos","micro_avg","macro_avg",'weighted_avg'],
                      figsize=(15,5),
                      save=save,
                      model_name=model_name,
                      folder_name=folder_name,
                      file_name='y_report')
    plt.show(), plt.close()
    
    youden_dic= {'model':model_name, 'auc':auc_y, 'precision':prec_y, 'recall':recall_y, 'f1':f1_y, 'npv':npv_y}
    return(youden_dic)<jupyter_output><empty_output><jupyter_text>testing global model  test entire trainset and predict trainset. * step1) hypertune xgb on 5fold cv.     * step2) run hypertuned model on 5fold cv with lr and get overall metrics.* step3) test entire train set and predict testset.* step4) local model testing thresholds: * Decreasing thresholds on the decision function used to compute    fpr and tpr. `thresholds[0]` represents no instances being predicted    and is arbitrarily set to `max(y_score) + 1`.<jupyter_code>#setting up test table

#print(pd.DataFrame({'xg':xg_cv_auc, 'lr':lr_cv_auc, 'rf':rf_cv_auc, 'svc':svc_cv_auc}))
#pd.DataFrame({'model':[],'auc':[], 'precision':[], 'recall':[], 'f1':[]})
test_summary_df= pd.DataFrame({'model':[],'auc':[], 'precision':[], 'recall':[], 'f1':[], 'npv':[]})
test_summary_df<jupyter_output><empty_output><jupyter_text>model fitting<jupyter_code>xgboost = reset_model('xgboost', hardcode=False)
xgboost.fit(x, y)

logreg = reset_model('logreg', hardcode=False)
logreg.fit(x, y)

rf= reset_model('rf', hardcode=False)
rf.fit(x,y)

# from sklearn.naive_bayes import GaussianNB
# gnb =GaussianNB()
# nb_y_pred = gnb.fit(x, y)

from sklearn import svm
svc= reset_model('svc', hardcode=False)
svc.fit(x, y)

knn= reset_model('knn', hardcode=False)
knn.fit(x,y)<jupyter_output><empty_output><jupyter_text>global model test set evaluation<jupyter_code>#old
#svc_eval= classifier_eval(svc, x=np.array(x_test), y=y_test, save=True, model_name='svc', folder_name='clinical_agg')
svc_eval= classifier_eval(svc, x=np.array(x_test), y=y_test, save=False, model_name='svc', folder_name='clinical_agg_elix')
svc_eval= classifier_eval(svc, x=np.array(x_test), y=y_test, save=True, model_name='svc', folder_name='clinical_agg_elix')
xgboost_eval= classifier_eval(xgboost, x=np.array(x_test), y=y_test, save=False, model_name='xgboost', folder_name='clinical_agg_elix')
xgboost2=XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
       colsample_bytree=0.8, gamma=0, learning_rate=0.1, max_delta_step=0,
       max_depth=2, min_child_weight=2, missing=None, n_estimators=100,
       n_jobs=-1, nthread=None, objective='binary:logistic',
       random_state=0, reg_alpha=0, reg_lambda=1, scale_pos_weight=1,
       seed=27, silent=True, subsample=0.8)

xgboost2.fit(x, y)

xgboost2_eval= classifier_eval(xgboost2, x=np.array(x_test), y=y_test, save=False, model_name='xgboost2', folder_name='clinical_agg_24')

xgboost_eval= classifier_eval(xgboost, x=np.array(x_test), y=y_test, save=False, model_name='xgboost', folder_name='clinical_agg_elix')
xgboost_eval= classifier_eval(xgboost, x=np.array(x_test), y=y_test, save=False, model_name='xgboost', folder_name='clinical_agg_elix')

xgboost_hyper
xgboost_hyper.best_params_
##without elix
#xgboost_eval= classifier_eval(xgboost, x=np.array(x_test), y=y_test, save=True, model_name='xgboost', folder_name='clinical_agg')
rf_eval= classifier_eval(rf, x=np.array(x_test), y=y_test, save=False, model_name='rf', folder_name='clinical_agg_elix')
#without elix
#rf_eval= classifier_eval(rf, x=np.array(x_test), y=y_test, save=True, model_name='rf', folder_name='clinical_agg')
# logreg_eval= classifier_eval(logreg, x=np.array(x_test), y=y_test)
logreg_eval= classifier_eval(logreg, x=np.array(x_test), y=y_test, save=False, model_name='logreg', folder_name='clinical_agg_elix')
#knn_eval= classifier_eval(knn, x=np.array(x_test), y=y_test, save=False, model_name='knn', folder_name='clinical_agg')
knn_eval= classifier_eval(knn, x=np.array(x_test), y=y_test, save=False, model_name='knn', folder_name='clinical_agg_elix')
knn_eval<jupyter_output><empty_output><jupyter_text>Ensemble of all models<jupyter_code>from sklearn.ensemble import VotingClassifier
#create a dictionary of our models
estimators=[("xgboost", xgboost), ('rf', rf), ('log_reg', logreg), ('svc',svc)]
#create our voting classifier, inputting our models
ensemble = VotingClassifier(estimators, voting='soft', n_jobs=-1)
# If ‘hard’, uses predicted class labels for majority rule voting.
# Else if ‘soft’, predicts the class label based on the argmax of the sums of the predicted probabilities,
# which is recommended for an ensemble of well-calibrated classifiers.

#weights: array-like, shape (n_classifiers,), optional (default=`None`)
#Sequence of weights (float or int) to weight the occurrences of predicted class labels (hard voting) or class probabilities before averaging (soft voting).
#Uses uniform weights if None.
ensemble.fit(x, y)#, sample_weight=np.array([0.67289604, 1.94595562]))
ensemble2 = VotingClassifier(estimators, voting='hard', n_jobs=-1)
ensemble2.fit(x, y)
y_pred = ensemble2.predict(np.array(x_test))
#y_proba = ensemble2.predict_proba(np.array(x_test))[:,1]
from sklearn import metrics
from sklearn.metrics import precision_score, roc_auc_score, f1_score, recall_score
pos_label=1
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred, pos_label=pos_label)
roc_auc = metrics.auc(fpr, tpr)

#gathering the optimal youden_index and df of tpr/fpr for auc and index of that optimal youden. idx is needed in the roc
youden_threshold, roc_df, idx= optimal_youden_index(fpr, tpr, thresholds,tp90=True)

#plotting roc
plot_roc(fpr, tpr, roc_auc, roc_df, idx, save=False, model_name=ensemble2,folder_name=None, file_name='roc')
plt.show(), plt.close()

#printing npv, recall, precision, accuracy
npv=confusion_matrix(y_test, y_pred)[0,0]/sum(np.array(y_pred)==0)*100
prec= precision_score(y_true=y_test, y_pred= y_pred, pos_label=pos_label)*100
recall= recall_score(y_true=y_test, y_pred= y_pred, pos_label=pos_label)*100
f1= f1_score(y_true=y_test, y_pred= y_pred, pos_label=pos_label)*100
acc=accuracy_score(y_test, y_pred)*100

print("npv:", npv,'\n')
print("prec:", prec,'\n')
print("recall:", recall,'\n')
print("f1:", f1,'\n')
print("acc:", acc,'\n')

hard_vote_summary={'model':"hard_voting_classifier",'auc':acc, 'precision':prec, 'recall':recall, 'f1':f1, 'npv':npv}
# test_summary_df= pd.DataFrame({'model':[],'auc':[], 'precision':[], 'recall':[], 'f1':[], 'npv':[]})
# test_summary_df
ensemble_eval= classifier_eval(ensemble, x=np.array(x_test), y=y_test, save=True, model_name='model_ensemble', folder_name='clinical_agg_elix')
np.unique(y)
from sklearn.utils.class_weight import compute_class_weight
compute_class_weight('balanced', np.unique(y), y)
#classifier_eval(gnb, x=x_test, y=y_test)
test_summary_df= pd.DataFrame([rf_eval,
                             logreg_eval,
                             xgboost_eval,
                             svc_eval,
                            ensemble_eval,
                              hard_vote_summary])
test_summary_df.set_index('model').round(decimals=2)<jupyter_output><empty_output><jupyter_text>variable importance<jupyter_code>def saveplot(figure_name,folder_name=None):
    """
    simple function for saving plots
    """
    if folder_name != None:
        address = 'figures/{}/'.format(folder_name)
    else:
        address = 'figures/'
    if not os.path.exists(address):
        os.makedirs(address)
    plt.savefig(address+"/{}.png".format(figure_name),bbox_inches='tight')
def var_imp(model,folder_name,model_name, n_var=4, save=True):
    model_name=type(model).__name__
    plot_title= "Top {} {} {} Variable Importance".format(n_var, folder_name,model_name)
    feat_importances = pd.Series(model.feature_importances_, index=x_train.columns)
    topn=feat_importances.nlargest(n_var).sort_values()
    ax=topn.plot(kind='barh', x='doop', title=plot_title)#.xlabel("xlab")
    ax.set_xlabel("Variable Importance")
    if save==True:
        saveplot(figure_name=plot_title, folder_name=folder_name)
    return(topn)

#     ###
#     imp= model.feature_importances_
#     var_index=[ x for x in range(0,len(rf.feature_importances_))]
#     variables=list(x_train)
#     return(pd.DataFrame({"imp":imp, 'index':var_index, 'variable': variables}).sort_values('imp', ascending=False))

#var_imp(rf,plot_title='RF_' n_var=4)
#var_imp(rf,"clinical_agg","RF", n_var=6, save=False)
var_imp(rf,"clinical_agg_elix","RF", n_var=20, save=True)
var_imp(xgboost2,"clinical_agg","xgboost", n_var=20, save=False)
var_imp(xgboost,"clinical_agg_elix_72","xgboost", n_var=10, save=True)
test_summary_df['model']
print(len(x_train),len(x_test))
len(x)<jupyter_output><empty_output><jupyter_text>model pickling<jupyter_code>import pickle
from sklearn import model_selection

# save the model to disk
model_list=[rf,logreg, xgboost,svc ]
for element in model_list:#test_summary_df['model']:
    filename = 'models/{}_{}_{}.sav'.format(date,dataset,type(element).__name__)
    #os.makedirs(os.path.dirname(filename), exist_ok=True)
    pickle.dump(element, open(filename, 'wb'))
    
    
# xgboost = reset_model('xgboost')
# xgboost.fit(x, y)

# logreg = reset_model('logreg')
# logreg.fit(x, y)

# rf= reset_model('rf')
# rf.fit(x,y)

# # from sklearn.naive_bayes import GaussianNB
# # gnb =GaussianNB()
# # nb_y_pred = gnb.fit(x, y)

# from sklearn import svm
# svc= reset_model('svc')
# svc.fit(x, y)
#$dataset= "clinagg_elix"
filename
os.path.exists(filename)
import pickle
from sklearn import model_selection

# save the model to disk
model_list=[rf,logreg, xgboost,svc ]
for element in model_list:#test_summary_df['model']:  
    filename = 'models/{}_{}_{}.sav'.format(date,dataset,type(element).__name__)
    pickle.dump(element, open(filename, 'wb'))

# # load the model from disk
# loaded_model = pickle.load(open(filename, 'rb'))
# result = loaded_model.score(X_test, Y_test)
# print(result)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13721300789454624, 'fraction_numerical': 0.01951908452234374, 'mean_word_length': 4.427378144654088, 'pattern_counts': {'":': 12, '<': 81, '<?xml version=': 0, '>': 92, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '423098', 'license': 'MIT', 'n_tokens_mistral': 21079, 'n_tokens_neox': 20035, 'n_words': 3482, 'path': 'notebooks/outdated_notebooks/08.20-imputation.ipynb', 'repo_name': 'geickelb/mimiciii-antibiotics-modeling', 'size': '55235'}