<jupyter_start><jupyter_text>[View in Colaboratory](https://colab.research.google.com/github/ckbjimmy/2018_mlw/blob/master/nb1_classification.ipynb)  Machine Learning for Clinical Predictive AnalyticsWe would like to introduce basic machine learning techniques and toolkits for clinical knowledge discovery in the workshop.The material will cover common useful algorithms for clinical prediction tasks, as well as the diagnostic workflow of applying machine learning to real-world problems. We will use [Google colab](https://colab.research.google.com/) / python jupyter notebook and two datasets:- Breast Cancer Wisconsin (Diagnostic) Database, and - pre-extracted ICU data from PhysioNet Database  to build predictive models.The learning objectives of this workshop tutorial are:- Learn how to use Google colab / jupyter notebook- Learn how to build machine learning models for clinical classification and/or clustering tasksTo accelerate the progress without obstacles, we hope that the readers fulfill the following prerequisites:- [Skillset] basic python syntax- [Requirements] Google account OR [anaconda](https://anaconda.org/anaconda/python)In part 1, we will go through the basic of machine learning for classification problems.In part 2, we will investigate more on unsupervised learning methods for clustering and visualization.In part 3, we will play with neural networks.  Part I – Classification (supervised learning)In the first part of the workshop, we would like to explore how to utilize machine learning algorithms to approach clinical predictive analytic problems, specifically, the classification problem.After going through this tutorial, we hope that you will understand how to use python and scikit-learn to design and build simple machine learning models for classification problems and how to evaluate the performance and make it be interpretable.The tutorial is modified from Dr. <NAME>'s version for 2017 workshop.In the beginning, we will install and import the packages needed for the workshop (we usually do this in python). Then we will work on a well-structured breast cancer dataset, then move on to the real world PhysioNet dataset.In the first chunk, it will take few seconds to install and import packages we need in the tutorial.Colab / jupyter notebook provides a convenient function that you can run the bash command by adding `!` before the command. e.g. `!pip install tensorflow` , `!ls`, `!clear`.<jupyter_code># try:
#     import psycopg2
# except:
#     !pip install psycopg2
  
try:
    import pydotplus
except:
    !pip install pydotplus
    
try:
    import graphviz
except:
    !apt-get install graphviz -y


from __future__ import print_function

import numpy as np
import pandas as pd
import sklearn
import sys
import datetime as dt

import matplotlib
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties # for unicode fonts

from collections import OrderedDict

from sklearn import datasets

# used to print out pretty pandas dataframes
from IPython.display import display, HTML

from sklearn.pipeline import Pipeline

# used to impute mean for data and standardize for computational stability
from sklearn.preprocessing import Imputer
from sklearn.preprocessing import StandardScaler

# logistic regression is our favourite model ever
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import LogisticRegressionCV # l2 regularized regression
from sklearn.linear_model import LassoCV

from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.ensemble import AdaBoostClassifier
from sklearn.naive_bayes import MultinomialNB
from sklearn.svm import LinearSVC
from sklearn.svm import SVC
from sklearn.calibration import CalibratedClassifierCV
from sklearn.neural_network import MLPClassifier         

from sklearn import model_selection

# used to calculate AUROC/accuracy
from sklearn import metrics
from sklearn.model_selection import cross_val_score

# used to create confusion matrix
from sklearn.metrics import confusion_matrix

%matplotlib inline<jupyter_output>Reading package lists... Done
Building dependency tree       
Reading state information... Done
graphviz is already the newest version (2.38.0-16ubuntu2).
0 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.<jupyter_text>We also defined the `plot_model_pred_2d` function for plotting. We will use this function later for visualizing the decision boundaries.<jupyter_code>def plot_model_pred_2d(mdl, X, y, feat):
    # look at the regions in a 2d plot
    # based on scikit-learn tutorial plot_iris.html

    # get minimum and maximum values
    x0_min = X[:, 0].min()
    x0_max = X[:, 0].max()
    x1_min = X[:, 1].min()
    x1_max = X[:, 1].max()

    xx, yy = np.meshgrid(np.linspace(x0_min, x0_max, 100),
                         np.linspace(x1_min, x1_max, 100))

    Z = mdl.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    # plot the contour - colouring different regions
    cs = plt.contourf(xx, yy, Z, cmap='Blues')

    # plot the individual data points - colouring by the *true* outcome
    color = y.ravel()
    plt.scatter(X[:, 0], X[:, 1], c=color, marker='o', s=40, cmap='Blues')

    plt.xlabel(feat[0])
    plt.ylabel(feat[1])
    plt.axis("tight")
    plt.colorbar()<jupyter_output><empty_output><jupyter_text>Breast cancer predictionLet's start from using the [breast cancer dataset in UCI data repository](https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.names) to have a quick look on how to do the analysis and build models using well-structured data (without missing data and other data cleaning problems). The python machine learning package `scikit-learn` has already helped us preprocess the data. We load the breast cancer dataset from `sklearn.datasets`, and show the description of this dataset.<jupyter_code># load the data
df_bc = datasets.load_breast_cancer()
print([k for k in df_bc.keys()])
print(df_bc['feature_names'])

# if you want a description of the dataset, uncomment the below line
print(df_bc['DESCR'])

# pick index of the features to use (only pick 2)
#    :Attribute Information (in order):
#        0 - CRIM     per capita crime rate by town
#        1 - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
#        2 - INDUS    proportion of non-retail business acres per town
#        3 - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
#        4 - NOX      nitric oxides concentration (parts per 10 million)
#        5 - RM       average number of rooms per dwelling
#        6 - AGE      proportion of owner-occupied units built prior to 1940
#        7 - DIS      weighted distances to five Boston employment centres
#        8 - RAD      index of accessibility to radial highways
#        9 - TAX      full-value property-tax rate per $10,000
#       10 - PTRATIO  pupil-teacher ratio by town
#       11 - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
#       12 - LSTAT    % lower status of the population
#       13 - MEDV     Median value of owner-occupied homes in $1000's<jupyter_output>['data', 'target', 'target_names', 'DESCR', 'feature_names']
['mean radius' 'mean texture' 'mean perimeter' 'mean area'
 'mean smoothness' 'mean compactness' 'mean concavity'
 'mean concave points' 'mean symmetry' 'mean fractal dimension'
 'radius error' 'texture error' 'perimeter error' 'area error'
 'smoothness error' 'compactness error' 'concavity error'
 'concave points error' 'symmetry error' 'fractal dimension error'
 'worst radius' 'worst texture' 'worst perimeter' 'worst area'
 'worst smoothness' 'worst compactness' 'worst concavity'
 'worst concave points' 'worst symmetry' 'worst fractal dimension']
Breast Cancer Wisconsin (Diagnostic) Database
=============================================

Notes
-----
Data Set Characteristics:
    :Number of Instances: 569

    :Number of Attributes: 30 numeric, predictive attributes and the class

    :Attribute Information:
        - radius (mean of distances from center to points on the perimeter)
        - texture (standard deviation of g[...]<jupyter_text>In the next chunk, we use the second and the last feature (mean texture and worst fractal dimension, `[1, 29]`) to build models to classify whether the lesion is tumor or not.We use the features and the corresponding label to build the model.For training machine learning models, we usually split our data into training and testing dataset, and only use the training subset for modeling (using `train_test_split` function). The purpose of this splitting step is to make the model less overfit and be more generalizable.For much robust models, we may further split the training set to training and validation sets. Check the figure to understand the relationship between training, validation and testing splitted sets.![split](http://www.codeproject.com/KB/AI/1146582/validation.PNG)[Source] http://magizbox.com/training/machinelearning/site/evaluation/<jupyter_code>idx = [1, 29] # choose the features for modeling, you can change the idx and see what will happen
X = df_bc['data'][:, idx] # X is your feature set
y = df_bc['target'] # y is the target, the label
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=42) # split your data to train/test sets

feat = [x for x in df_bc['feature_names'][idx]]<jupyter_output><empty_output><jupyter_text>Above, we've extracted only two features out of the breast cancer dataset as predictors.Then, we run six different machine learning algorithms to build the classifier for predicting the lesion is malignant or not.We quickly visualize the results of all the models we have presented here. In the figure, we can see the decision boundaries created by the models.These boundaries can be a source of interpretability of the clinical predictive model.We also print out their performance as measured by the **Area Under the Receiver Operator Characteristic curve (AUROC)**, which is a commonly used metrics in the machine learning world.<jupyter_code># we put all algorithms in a dictionary, you can add/remove them
# there are many algorithms in scikit-learn you can try: 
# http://scikit-learn.org/stable/supervised_learning.html#supervised-learning
# the `.fit` step is the training process
clf = dict()
clf['Logistic Regression'] = LogisticRegression(fit_intercept=True).fit(X_train,y_train)
clf['Decision Tree'] = DecisionTreeClassifier(criterion='entropy', splitter='best').fit(X_train,y_train)
clf['Gradient Boosting'] = GradientBoostingClassifier(n_estimators=10).fit(X_train, y_train)
clf['Random Forest'] = RandomForestClassifier(n_estimators=10).fit(X_train, y_train)
clf['Bagging'] = BaggingClassifier(n_estimators=10).fit(X_train, y_train)
clf['AdaBoost'] = AdaBoostClassifier(n_estimators=10).fit(X_train, y_train)

# visualization
fig = plt.figure(figsize=[16,9])

print('AUROC\tModel')
for i, curr_mdl in enumerate(clf):
    yhat = clf[curr_mdl].predict_proba(X_test)[:,1] # prediction
    score = metrics.roc_auc_score(y_test, yhat) # get AUROC
    print('{:0.3f}\t{}'.format(score, curr_mdl))
    ax = fig.add_subplot(2,3,i+1)
    plot_model_pred_2d(clf[curr_mdl], X_test, y_test, feat) # we define this function in the beginning
    plt.title(curr_mdl)
    
plt.show()<jupyter_output>AUROC	Model
0.780	Logistic Regression
0.687	Decision Tree
0.807	Gradient Boosting
0.781	Random Forest
0.804	Bagging
0.801	AdaBoost<jupyter_text>Here we can see that quantitatively, Logistic Regression, AdaBoost and Gradient Boosting have produced the highest discrimination among all the models (~0.80). The decision surfaces of these models also seem simpler, and less "noisy", which is likely the reason for the improved **generalization** on the held out test set. K-fold cross-validationNow we will include all features of the breast cancer dataset - we will no longer be able to easily visualize the models, but we will be able to better evaluate them. We'll also switch to using **5-fold cross-validation** to get a good estimate of the generalization performance of the model.For the k-fold cross-validation, you can still split the data to train/test sets to make your models much robust.The figure demonstrates the simple 4-fold cross-validation.![cv](https://upload.wikimedia.org/wikipedia/commons/1/1c/K-fold_cross_validation_EN.jpg)[Source] https://en.wikipedia.org/wiki/Cross-validation_(statistics)<jupyter_code># use all features
X = df_bc['data']
y = df_bc['target']

# use cross-validation to estimate the performance of each model
print('Acc\tAUROC\tModel')
for curr_mdl in clf:
    scores = cross_val_score(clf[curr_mdl], X, y, cv=5, scoring='accuracy')
    auc = cross_val_score(clf[curr_mdl], X, y, cv=5, scoring='roc_auc')
    print('{:0.3f}\t{:0.3f}\t{}'.format(scores.mean(), auc.mean(), curr_mdl))<jupyter_output>Acc	AUROC	Model
0.953	0.992	Logistic Regression
0.931	0.928	Decision Tree
0.947	0.981	Gradient Boosting
0.947	0.988	Random Forest
0.953	0.980	Bagging
0.946	0.988	AdaBoost<jupyter_text>We note two things here:1. by **using the entire feature set** we have dramatically improved performance of the model, indicating that there was more information contained in the other columns, and2. most of our models are performing relatively equivalently (except for the super simple model, the decision tree), and do not forget **logistic regression**!  PhysioNetWe will now practice using these models on a dataset acquired from patients admitted to intensive care units at the Beth Israel Deaconness Medical Center in Boston, MA. All patients in the cohort stayed for at least 48 hours, and the goal of the prediction task is to predict in-hospital mortality. If you're interested, you can read more about the dataset [here](http://physionet.org/challenge/2012/).The data is originally provided as hourly observations for a number of variables, and the preprocessing step involved extracting summary statistics across all these observations.The outcome is the first column `hospitalmortality` (`0` means alive while discharge).The rest of the data are features you can use to predict this binary outcome.We will again try to build the classifier using logistic regression.You may want to use other algorithms such as AdaBoost, Random Forest, Bagging, and Gradient Boosting. Pick your favourite and play with the parameters to see how well you can do.First, we start from loading the data and look at all columns in the dataset---there are 182 possible predictors and an outcome variable.This time we read the data in a slightly different way using pandas `dataframe` format. pandas `dataframe` is useful for importing the tabular data/spreadsheet.<jupyter_code>df = pd.read_csv('https://raw.githubusercontent.com/ckbjimmy/2018_mlw/master/data/PhysionetChallenge2012_data.csv')
df.head()
# check our features
print([x for x in df.columns])<jupyter_output>['hospitalmortality', 'ALPMin', 'ALPMax', 'ALPFirst', 'ALPLast', 'ALPMedian', 'ALTMin', 'ALTMax', 'ALTFirst', 'ALTLast', 'ALTMedian', 'ASTMin', 'ASTMax', 'ASTFirst', 'ASTLast', 'ASTMedian', 'Age', 'AlbuminMin', 'AlbuminMax', 'AlbuminFirst', 'AlbuminLast', 'AlbuminMedian', 'BUNMin', 'BUNMax', 'BUNFirst', 'BUNLast', 'BUNMedian', 'BilirubinMin', 'BilirubinMax', 'BilirubinFirst', 'BilirubinLast', 'BilirubinMedian', 'CholesterolFirst', 'CholesterolLast', 'CreatinineMin', 'CreatinineMax', 'CreatinineFirst', 'CreatinineLast', 'CreatinineMedian', 'DiasABPMin', 'DiasABPMax', 'DiasABPFirst', 'DiasABPLast', 'DiasABPMedian', 'FiO2Min', 'FiO2Max', 'FiO2First', 'FiO2Last', 'FiO2Median', 'GCSMin', 'GCSMax', 'GCSFirst', 'GCSLast', 'GCSMedian', 'Gender', 'GlucoseMin', 'GlucoseMax', 'GlucoseFirst', 'GlucoseLast', 'GlucoseMedian', 'HCO3Min', 'HCO3Max', 'HCO3First', 'HCO3Last', 'HCO3Median', 'HCTMin', 'HCTMax', 'HCTFirst', 'HCTLast', 'HCTMedian', 'HRMin', 'HRMax', 'HRFirst', 'HRLast', 'HRMedian', 'Height'[...]<jupyter_text>Next, we split features and outcome variables.<jupyter_code># first column [0] is target (outcome variable)
y = df['hospitalmortality'].values

# second column onward [1:] are features 
X = df.drop(['hospitalmortality'], axis=1).values

# we keep the feature name in X_header
X_header = df.drop(['hospitalmortality'],axis=1).columns
# train a simple logistic regression, and...
clf = LogisticRegression(fit_intercept=True).fit(X, y)<jupyter_output><empty_output><jupyter_text>The above error of  `ValueError: Input contains NaN, infinity or a value too large for dtype('float64').` tells us that we have missing data! (Alistair: Dramatic music plays).Good thing about Google colab---they have a button for you to search stackoverflow.We then check the `NaN` using the transpose view---there are a lot `NaN` in our dataset.<jupyter_code>df.head().T<jupyter_output><empty_output><jupyter_text>Missing data imputationIn the above transposed view, we can see a few of the missing values for all predictors. There are many methods which can be used for missing data imputation. Of those tools, many of them are available in the `imputer` module of scikit-learn.Here is a [brief tutorial of imputing missing values](http://scikit-learn.org/stable/modules/preprocessing.htmlimputation).You can [read more about the imputer module here](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.Imputer.html).We choose to take a reasonably simple approach, and replace all missing values (`'NaN'`) with the mean value of observed measurements (`strategy='mean'`).<jupyter_code># define the method of imputation
imp = Imputer(missing_values='NaN', strategy='mean', axis=0)

# learn the parameters of that imputation using our data X
# i.e., in our case, calculate the means of each column which we will use to impute data
imp = imp.fit(X)

# apply that imputation
X_transform = imp.transform(X)<jupyter_output><empty_output><jupyter_text>We can now see that our missing data has been replaced:<jupyter_code>print('Before:', end='')
print(X[0:10, 0])

print('After:', end='')
print(X_transform[0:10, 3])<jupyter_output>Before:[ nan  nan 105. 105.  nan 101.  nan  47.  nan  nan]
After:[104.30236686 104.30236686 105.         105.         104.30236686
 101.         104.30236686  47.         104.30236686 104.30236686]<jupyter_text>For linear models (e.g. logistic regression), this will be roughly equivalent to imputing the average risk for this feature for patients missing that data.We can now train a model on the data!<jupyter_code>clf = LogisticRegression(fit_intercept=True)
clf.fit(X_transform, y)<jupyter_output><empty_output><jupyter_text>We then evaluate the model using the AUROC again. An AUROC = 0.5 represents a random classifier making random predictions, while an AUROC = 1.0 represents a perfect classifier which always outputs the correct prediction.<jupyter_code>yhat = clf.predict_proba(X_transform)[:,1]
auc = metrics.roc_auc_score(y, yhat)

print('{:0.3f} - AUROC of model (training set).'.format(auc))<jupyter_output>0.877 - AUROC of model (training set).<jupyter_text>Normalizartion / Find the important featuresLet's dig into the model to determine what is driving the predictions. In healthcare applications, model interpretability is very important as it can be used to validate the model's construct validity.<jupyter_code># plot the coefficients
coef = np.row_stack([df.columns[1:], clf.coef_]).T

idxSort = np.argsort(coef[:,1])

print('Top 5 predictors negatively correlated with mortality:')
for n in range(5):
    print('{:1.2f} - {}'.format(coef[idxSort[n],1], coef[idxSort[n],0]))
    
print()
print('Top 5 predictors positively correlated with mortality:')
for n in range(5):
    print('{:1.2f} - {}'.format(coef[idxSort[-n-1],1], coef[idxSort[-n-1],0]))<jupyter_output>Top 5 predictors negatively correlated with mortality:
-0.35 - TempMedian
-0.34 - CSRU
-0.29 - BilirubinFirst
-0.24 - GCSLast
-0.21 - MechVentFirst

Top 5 predictors positively correlated with mortality:
0.42 - TroponinMeasured
0.26 - FiO2Last
0.25 - FiO2Median
0.22 - SICU
0.21 - MICU<jupyter_text>One issue with these coefficients is that they are unadjusted---that is they are based on the original scale of the data. For example, the coefficient for FiO2, shown in percentage (0-100), will likely be ~100 times higher than the coefficient for mechanical ventilation (0/1). While this is fine for modelling purposes, it also means we can't directly compare the magnitudes of the coefficients to determine the magnitude of the association between the feature and risk.Of course there is a solution! We can scale the data using `StandardScaler` function in scikit-learn.<jupyter_code># problem: the above coefficients all operate on different scales
# solution: scale the data before we use it in the model
imp = Imputer(missing_values='NaN', strategy="mean", axis=0)
imp = imp.fit(X)
X_transform = imp.transform(X)

scale = StandardScaler()
scale = scale.fit(X_transform)
X_transform = scale.transform(X_transform)
# fit the model on scaled data
clf = LogisticRegression(fit_intercept=True)
clf = clf.fit(X_transform, y)

# evaluate the model
yhat = clf.predict_proba(X_transform)[:,1]
auc = metrics.roc_auc_score(y, yhat)

print('{:0.3f} - AUROC of model (training set).'.format(auc))

# print the coefficients
coef = np.row_stack([df.columns[1:], clf.coef_]).T
idxSort = np.argsort(coef[:,1])

print('Top 5 predictors negatively correlated with mortality:')
for n in range(5):
    print('{:1.2f} - {}'.format(coef[idxSort[n],1], coef[idxSort[n],0]))
    
print()
print('Top 5 predictors positively correlated with mortality:')
for n in range(5):
    print('{:1.2f} - {}'.format(coef[idxSort[-n-1],1], coef[idxSort[-n-1],0]))<jupyter_output>0.881 - AUROC of model (training set).
Top 5 predictors negatively correlated with mortality:
-0.89 - ALPMax
-0.83 - GCSLast
-0.79 - pHMin
-0.71 - BilirubinFirst
-0.58 - ASTFirst

Top 5 predictors positively correlated with mortality:
0.73 - BilirubinMedian
0.69 - ASTMin
0.67 - MechVentDuration
0.64 - BUNLast
0.64 - pHMedian<jupyter_text>Let's look at the odds ratios for these coefficients graphically. This will give us a more intuitive comparison between the features. For those which have negative correlation with bad outcome (i.e. higher values indicate better prognosis), we reverse the coefficient on the graph so it is more comparable.<jupyter_code># plot the top 10 lowest / highest
coef = np.row_stack([df.columns[1:], clf.coef_]).T
idxSort = np.argsort(coef[:,1])

idxTop = idxSort[0:10]
idxBot = idxSort[-1:-11:-1]

colTop = [0.8906,0.1016,0.1094]
colBot = [0.2148,0.4922,0.7188]

f = plt.figure(figsize=[10,6])
ax1 = f.add_subplot(111)
ax1.plot( np.exp(-coef[idxTop,1].astype(float)), range(10), 's', markersize=10, color = colTop)

# move negatively correlated ticklabels over to the right
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position("right")
ax1.set_ylim([-1,10])

ax1.yaxis.set_ticks(range(10))
ax1.yaxis.set_ticklabels(coef[idxTop,0], color=colTop, fontsize=16)
plt.ylabel("Negatively correlated", color=colTop, fontsize=16)

ax2 = f.add_subplot(111, sharex=ax1, frameon=False)
ax2.plot( np.exp(coef[idxBot,1].astype(float)), range(10), 'o', markersize=10, color=colBot )


ax2.yaxis.set_ticks(range(10))
ax2.yaxis.set_ticklabels(coef[idxBot,0], color=colBot, fontsize=16)
ax2.set_ylim([-1,10])
plt.ylabel("Positively correlated", color=colBot, fontsize=16)

# set axes
plt.show()<jupyter_output><empty_output><jupyter_text>PipelineWe can use scikit-learn's pipeline feature to combine together all of these steps (imputation, scaling, and final model building). This is technically the same thing, but syntactically it's a lot cleaner.<jupyter_code># can use pipelines as a nicer way of doing the above
estimator = Pipeline([
        ("imputer", Imputer(missing_values='NaN',
                                  strategy="mean",
                                  axis=0)),
        ("scaler", StandardScaler()),
        ("logreg", LogisticRegression(fit_intercept=True))
    ])

mdl_pipeline = estimator.fit(X, y)

yhat = mdl_pipeline.predict_proba(X)[:,1]

auc = metrics.roc_auc_score(y, yhat)
print('{:0.3f} - AUROC of model (training set).'.format( auc ))<jupyter_output>0.881 - AUROC of model (training set).<jupyter_text>Note that we achieve the same performance as before---again, this is simply a short form of what we have done before.You'll note that in the print statement we keeps stating `"training set"`, and may be wondering what this means and why it matters. The training set corresponds to the data used to develop the model.Importantly, this means that the model has seen the data that we are using to estimate how well the model is doing. For the regression model above, the consequences of this may not be obvious. Let's switch to a more flexible model and see what happens.<jupyter_code># can use pipelines as a nicer way of doing the above
estimator = Pipeline([
        ("imputer", Imputer(missing_values='NaN',
                                  strategy="mean",
                                  axis=0)),
        ("scaler", StandardScaler()),
        ("rf", RandomForestClassifier())
    ])

mdl_pipeline = estimator.fit(X, y)

yhat = mdl_pipeline.predict_proba(X)[:,1]

auc = metrics.roc_auc_score(y, yhat)
print('{:0.3f} - AUROC of model (training set).'.format( auc ))<jupyter_output>1.000 - AUROC of model (training set).<jupyter_text>Amazing! We have an AUROC = 1.0, which is perfect, therefore we have a perfect model. Time to go home!Obviously, it's not that simple, though we wish it were :) We have used a model called a Random Forest. Omitting a few details, a random forest is a collection of decision trees. Each decision tree tries to perfectly classify a small subset of the overall data. When combined together, all of these decision trees perfectly classify the entire dataset. However, that does not mean that they work perfectly.How can we check this? Well, let's hide some data from the classifier, and see how well it does on this hidden data. This is a very common concept in machine learning: we have a training set, and we hide a test set to see how we're really doing---as we did for the breast cancer prediction task.You can either use the splitting function in `scikit-learn` (in the breast cancer example), or use the simple random indexing approach to split the data.Here we use the latter approach.We realized that the result of AUC = 1.0 comes from overfitting of training data.You may change to other algorithms in the `Pipeline` and see which algorithm has better generalizability for testing data.  Prevent from overfitting<jupyter_code># we use 80% data for training and 20% for testing
np.random.seed(42)
idxTrain = np.random.rand(X.shape[0])
idxTest = idxTrain > 0.8
idxTrain = ~idxTest

X_train = X[idxTrain, :]
y_train = y[idxTrain]
X_test = X[idxTest, :]
y_test = y[idxTest]
# can use pipelines as a nicer way of doing the above
estimator = Pipeline([
        ("imputer", Imputer(missing_values='NaN',
                                  strategy="mean",
                                  axis=0)),
        ("scaler", StandardScaler()),
        ("rf", RandomForestClassifier())
    ])

mdl_pipeline = estimator.fit(X_train, y_train)

yhat = mdl_pipeline.predict_proba(X_train)[:,1]
auc = metrics.roc_auc_score(y_train, yhat)
print('{:0.3f} - AUROC of model (training set).'.format(auc))

yhat = mdl_pipeline.predict_proba(X_test)[:,1]
auc = metrics.roc_auc_score(y_test, yhat)
print('{:0.3f} - AUROC of model (test set).'.format(auc))<jupyter_output>1.000 - AUROC of model (training set).
0.777 - AUROC of model (test set).<jupyter_text>InterpretabilityNow that we have prevented our model from memorizing the data, we can get a much better estimate of how well we are actually doing. Note that with the more flexible model, we are not necessarily doing better. However, there are many "hyperparameters" for the Random Forest (different knobs we can tweak in order for things to work better). With some optimization of these hyperparameters, we may be able to do better. Note that this improvement in performance is at the cost of some interpretability of the model. We no longer have a simple set of coefficients to interpret: we now have hundreds of decision trees simultaneously being used to determine patient outcome.All is not lost! We can still garner some meaning from the models. One common method of interpretting the model is figuring out how important each feature is in the model. The method for doing this is quite intuitive: we can mess up a single column (randomize it), and see how much worse our model does. If that column is very important, then we would expect our model to be much worse. If the column is not that important, then our model performance would be about the same.<jupyter_code>D = X_test.shape[1]
# we should repeat the data shuffling to get a decent estimate of feature importance
B = 10
feat_importance = np.zeros([B, D], dtype=float)


# let's get the baseline performance of our model
yhat = mdl_pipeline.predict_proba(X_test)[:,1]
loss_base = metrics.log_loss(y_test, yhat)


# for each feature
for d in range(D):
    for b in range(B):
        # generate a random shuffle
        idxShuffle = np.random.randint(low=0, high=D, size=X_test.shape[0])
        # note this random shuffle allows the same row to be selected more than once
        # this is intentional, and called "bootstrapping", you can read up on it online

        X_fi = X_test

        # shuffle the data
        X_fi[:, d] = X_fi[idxShuffle, d]

        # apply our model
        yhat = mdl_pipeline.predict_proba(X_fi)[:,1]

        # calculate performance on the data with a shuffled column
        loss = metrics.log_loss(y_test, yhat)

        # store this in the feature importance matrix
        # note we subtract this off the baseline AUROC
        feat_importance[b,d] = loss_base - loss
        
    if np.mod(d,10)==0:
        print('Finished {:2.1f}% of features.'.format(d*100.0/D))
        
print('Done!')<jupyter_output>Finished 0.0% of features.
Finished 5.5% of features.
Finished 11.0% of features.
Finished 16.5% of features.
Finished 22.0% of features.
Finished 27.5% of features.
Finished 33.0% of features.
Finished 38.5% of features.
Finished 44.0% of features.
Finished 49.5% of features.
Finished 54.9% of features.
Finished 60.4% of features.
Finished 65.9% of features.
Finished 71.4% of features.
Finished 76.9% of features.
Finished 82.4% of features.
Finished 87.9% of features.
Finished 93.4% of features.
Finished 98.9% of features.
Done!<jupyter_text>Now that we've calculated the importance, let's plot the most important 10 features. We can define the most important features as those which caused the largest degradation in performance.<jupyter_code># take mean degradation in performance
feat_imp_mu = np.mean(feat_importance, axis=0)
feat_imp_iqr = np.percentile(feat_importance, 75, axis=0) - np.percentile(feat_importance, 25, axis=0)

# sort to get the features which had the most impact
idxSort = np.argsort(feat_imp_mu)

plt.figure(figsize=[12,8])
# note we negate the feature importance - now higher values indicate more important
plt.barh(range(10), -feat_imp_mu[idxSort[0:10]], color='r',
         yerr = feat_imp_iqr[idxSort[0:10]],
        tick_label = X_header[idxSort[0:10]],
        align='center')
plt.show()<jupyter_output><empty_output><jupyter_text>Note that we have changed our metric from the AUROC to the logistic loss. All of the top 10 features make sense: lab values are useful in determining patient outcome.There are other methods for determining feature importance---in fact, some are included in the random forest model object. The only trick is that they are buried within the pipeline that we used earlier.<jupyter_code>rf_mdl = mdl_pipeline.named_steps['rf']

# get feature importances
feat_imp_rf = rf_mdl.feature_importances_

# same as before - let's plot it!
idxSort = np.argsort(feat_imp_rf)

plt.figure(figsize=[12,8])
plt.barh(range(10), feat_imp_rf[idxSort[0:10]],
        tick_label = X_header[idxSort[0:10]],
        align='center')
plt.show()<jupyter_output><empty_output><jupyter_text>Here we get very different feature importances. Why is this? It helps to understand how these feature importances are derived.The first feature importance (log loss approach) was done based on overall performance: potassium ended up being very useful across all patients in predicting mortality.Alternatively, the second method is based on the Gini importance (in random forest). Remember that a random forest is made up of many trees. This second method says: for each split that uses the feature, how much "purer" are the two subsequent nodes. Below, the percentage is the frequency of the outcome, where we start with 50% (equal number of both groups in that node).```50%|| split on feature here|\| \|  \|   \|    \25%  75%```Above, after we split on the feature, our nodes are "purer". That is, the left node has more negative outcomes (0s), and the right node as more positive outcomes (1s).From this, we can conclude that MAP (mean arterial pressure) is very useful for separating those who survive from those who die, which makes sens.Patients who have low MAP are very ill, and it is an effective indicator of illness. However, not all patients have low MAP. Across all patients, potassium is a more important measure in indicating illness.  K-fold cross-validation (in detail)Finally, let's review a very effective technique for estimating how well our classifier is doing---cross-validation. Recall that we split the dataset into two sets: training and test, specifically 80% training and 20% test. As we tested using a smaller dataset (only 20% of the total records), we likely will end up with a noisier estimate of performance than if we used the entire dataset. Of course we can't use the entire dataset because then our model results aren't believable. Cross-validation is a clever way to use the entire dataset and use held-out data for proper evaluation. Each observation is assigned to a group. Let's say we have $k$ groups. We then can build a model using data in the $k$-th group as the test set, and all other data as the training set. Let's take a look at some code.<jupyter_code># create the K fold indices
K = 5

idxK = np.random.permutation(X.shape[0])
idxK = np.mod(idxK,K)

print('First 10 observations:',end=' ')
print(idxK[0:10])
print('idxK shape:',end=' ')
print(idxK.shape)<jupyter_output>First 10 observations: [3 3 2 3 3 1 4 0 3 0]
idxK shape: (4000,)<jupyter_text>Above you can see we have generated an integer, up to $k-1$ (but including 0), for all observations in the training set. We will now iterate through the unique integers to "hold out" a subset of the data for each model training.<jupyter_code>xval_auc = list()

for k in range(K):
    idxTest = idxK == k
    idxTrain = idxK != k

    X_train = X[idxTrain, :]
    y_train = y[idxTrain]
    X_test = X[idxTest, :]
    y_test = y[idxTest]
    
    # train the model using all but the kth fold
    curr_mdl = estimator.fit(X_train, y_train)
    
    # evaluate using the kth fold
    yhat = curr_mdl.predict_proba(X_test)[:,1]
    curr_score = metrics.roc_auc_score(y_test, yhat)
    xval_auc.append( curr_score )
    

    print('{} - Finished fold {} of {}. AUROC {:0.3f}.'.format(dt.datetime.now(),
                                                               k+1, K, curr_score))<jupyter_output>2018-06-16 14:18:31.389958 - Finished fold 1 of 5. AUROC 0.776.
2018-06-16 14:18:31.628354 - Finished fold 2 of 5. AUROC 0.784.
2018-06-16 14:18:31.874355 - Finished fold 3 of 5. AUROC 0.770.
2018-06-16 14:18:32.120065 - Finished fold 4 of 5. AUROC 0.761.
2018-06-16 14:18:32.354571 - Finished fold 5 of 5. AUROC 0.782.<jupyter_text>Of course, you can also use the function `cross_val_score` as we used in the breast cancer prediction for cross-validation.  Hyperparameter optimizationLet's look at our estimator.<jupyter_code>estimator = Pipeline([
        ("imputer", Imputer(missing_values='NaN',
                                  strategy="mean",
                                  axis=0)),
        ("scaler", StandardScaler()),
        ("rf", RandomForestClassifier())
    ])

print(estimator.named_steps['rf'])<jupyter_output>RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)<jupyter_text>Careful observation will note a few parameters which we could play around with:* criterion - what method to evaluate how well a node is splitting the data* max_depth - a limit on how deep the tree is* max_features - how many to look at for an inidividual tree* max_leaf_nodes - maximum number of leaf or end nodes* min_impurity_split - minimum impurity to allow for a split (the criterion is impurity)* min_samples_leaf - minimum samples in a leaf node* min_samples_split - minimum samples required to split a node* min_weight_fraction_leaf - similar to min_samples_leaf, except sometimes we are using weighted data, and so this takes that into account* n_estimators - the number of trees to learnMore details of parameteres can be found in the [`scikit-learn` document about random forest classifier](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html).You'll note that all of these parameters have values - but who is to say these values are optimal? They have been set to sensible defaults - all of which perform well on a variety of problems - but there's no reason that they would perform best on our problem! What we need to do is *hyperparameter optimization*: the above are "hyperparameters" (i.e., they are not the parameters the model uses to make decisions, but they are used to configure the model itself). We can optimize these in a few ways, but the simplest to understand is grid search. Grid search is essentially a brute force method of trying a bunch of combinations and picking the one that does best.Happily, `scikit-learn` provides some tools to do grid search using cross-validation, which is exactly what we want to do.<jupyter_code>params = {'rf__max_depth': [10],
         'rf__n_estimators': [10,50]}

estimator = Pipeline([
        ("imputer", Imputer(missing_values='NaN',
                                  strategy="mean",
                                  axis=0)),
        ("scaler", StandardScaler()),
        ("rf", RandomForestClassifier())
    ])

curr_mdl = model_selection.GridSearchCV(estimator, params,
                   scoring='roc_auc',
                   verbose=2, refit=True)
curr_mdl.fit(X, y)

curr_mdl.cv_results_['mean_test_score']
for x in curr_mdl.cv_results_:
    print(x,end=': ')
    print(curr_mdl.cv_results_[x])<jupyter_output>mean_fit_time: [0.13832275 0.56821219]
std_fit_time: [0.00238568 0.00405725]
mean_score_time: [0.00850328 0.01606464]
std_score_time: [1.34339696e-05 3.31599248e-04]
param_rf__max_depth: [10 10]
param_rf__n_estimators: [10 50]
params: [{'rf__max_depth': 10, 'rf__n_estimators': 10}, {'rf__max_depth': 10, 'rf__n_estimators': 50}]
split0_test_score: [0.80967704 0.83168913]
split1_test_score: [0.77325759 0.82610966]
split2_test_score: [0.80500587 0.8269249 ]
mean_test_score: [0.79597565 0.82824189]
std_test_score: [0.01618285 0.00246108]
rank_test_score: [2 1]
split0_train_score: [0.97978629 0.99596622]
split1_train_score: [0.98850215 0.99508019]
split2_train_score: [0.98258944 0.99499918]
mean_train_score: [0.98362596 0.99534853]
std_train_score: [0.00363293 0.00043803]<jupyter_text>More models!We are now armed with the ability to train and effectively evaluate our models. Now the fun begins. Let's evaluate a few different models to see which one does best.We use all techniques and tips above in this chunk for modeling.You can even find algorithms you like in the [`scikit-learn` document](http://scikit-learn.org/stable/supervised_learning.htmlsupervised-learning), add them into `models`, and fine-tune the hyperparameters.Or you can also change the $k$ value for cross-validation.<jupyter_code># Rough timing info:
#     rf - 3 seconds per fold
# logreg - 4 seconds per fold
#  lasso - 8 seconds per fold
models = {'l1': LogisticRegression(penalty='l1', multi_class='ovr', n_jobs=-1),
          'l2': LogisticRegression(penalty='l2', multi_class='ovr', n_jobs=-1), 
          #'nb': MultinomialNB(), 
          'svc': CalibratedClassifierCV(base_estimator=LinearSVC(penalty='l2', loss='squared_hinge', C=1.0, multi_class='ovr', random_state=0, max_iter=1000), cv=5), 
          'rbf': SVC(kernel='rbf', probability=True, decision_function_shape='ovr'), 
          'rf': RandomForestClassifier(n_estimators=100, n_jobs=-1), 
          'adab': AdaBoostClassifier(n_estimators=100), 
          'gb': GradientBoostingClassifier(n_estimators=100, learning_rate=1.0, max_depth=1, random_state=0),
          'mlp':MLPClassifier(solver='adam', alpha=1e-5, hidden_layer_sizes=(128, 64, 32), random_state=42)
         }       
                
# create k-fold indices
K = 5 # number of folds
idxK = np.random.permutation(X.shape[0])
idxK = np.mod(idxK,K)

mdl_val = dict()
results_val = dict()
pred_val = dict()
tar_val = dict()

for mdl in models:
    print('=============== {} ==============='.format(mdl))
    mdl_val[mdl] = list()
    results_val[mdl] = list() # initialize list for scores
    pred_val[mdl] = list()
    tar_val[mdl] = list()

    if mdl == 'xgb':
        # no pre-processing of data necessary for xgb
        estimator = Pipeline([(mdl, models[mdl])])

    else:
        estimator = Pipeline([("imputer", Imputer(missing_values='NaN',
                                          strategy="mean",
                                          axis=0)),
                      ("scaler", StandardScaler()),
                      (mdl, models[mdl])]) 

    for k in range(K):
        # train the model using all but the kth fold
        curr_mdl = estimator.fit(X[idxK != k, :],y[idxK != k])

        # get prediction on this dataset
        if mdl == 'lasso':
            curr_prob = curr_mdl.predict(X[idxK == k, :])
        else:
            curr_prob = curr_mdl.predict_proba(X[idxK == k, :])
            curr_prob = curr_prob[:,1]
        
        pred_val[mdl].append(curr_prob)
        tar_val[mdl].append(y[idxK == k])
        
        # calculate score (AUROC)
        curr_score = metrics.roc_auc_score(y[idxK == k], curr_prob)

        # add score to list of scores
        results_val[mdl].append(curr_score)

        # save the current model
        mdl_val[mdl].append(curr_mdl)
        
        print('{} - Finished fold {} of {}. AUROC {:0.3f}.'.format(dt.datetime.now(), k+1, K, curr_score))<jupyter_output>=============== l1 ===============
{'fraction_non_alphanumeric': 0.09666011787819254, 'fraction_numerical': 0.029030394082977, 'mean_word_length': 4.299608035276825, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 93, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '458491', 'license': 'MIT', 'n_tokens_mistral': 13994, 'n_tokens_neox': 12920, 'n_words': 4875, 'path': 'nb1_classification.ipynb', 'repo_name': 'ckbjimmy/2018_mlw', 'size': '43275'}