<jupyter_start><jupyter_text>Background to osteoarthritis case studyThis is day 4 from the [5-day JADS NHS PROMs data science case study](https://github.com/jads-nl/execute-nhs-proms/blob/master/README.md). To recap the previous lectures:- In **lecture 1** we focused on data understanding and explored various possible outcome parameters.- In **lecture 2** we constructed a combined outcome parameter using cut-off points for pain and physical functioning.- In **lecture 3** we performed regression and linear modeling on the outcome paramter `t1_eq_vas`.Today we are going to focus on classifcation.  Modeling: clustering & classfication  Learning objectives Modeling: classificationRecall we defined different outcome Y suitable for classification in lecture 2:- `y_mcid`: outcome is good (True/1) if `t1_oks_score` is above the threshold or `delta_oks_score` is larger than MCID- `y_t1_pain_good` and `y_t1_functioning_good`: combination of two binary outcomes, yielding a total of 4 classes.We will use the second outcome since it is more challenging and versatile (see [this presentation (in Dutch)](https://kapitan.net/wp-content/uploads/2018/11/181108_data_driven_healthcare_congres.pdf) for more details):- perform a binary classification, defining good outcome where both pain and functioning is good at `t1`- perform multiclass classification.We will compare the performance of single estimators (SVM, Decision tree) with ensemble learning (Random Forest, Gradient Boosting). Grid search and cross validation is done for model assessment. Some examples are given to visualize scoring functions of models to aid in model selection. Python: Hands-on Machine Learning (2nd edition)- [Classification (chapter 3)](https://github.com/ageron/handson-ml2/blob/master/03_classification.ipynb)- [Support-vector machines (chapter 5)](https://github.com/ageron/handson-ml2/blob/master/05_support_vector_machines.ipynb)- [Decision trees (chapter 6)](https://github.com/ageron/handson-ml2/blob/master/06_decision_trees.ipynb)- [Ensemble learning and random forests (chapter 7)](https://github.com/ageron/handson-ml2/blob/master/07_ensemble_learning_and_random_forests.ipynb) Python: scikit-learn and matplotlib- [Demonstration of multi-metric evaluation on cross_val_score and GridSearchCV](https://scikit-learn.org/stable/auto_examples/model_selection/plot_multi_metric_evaluation.htmlsphx-glr-auto-examples-model-selection-plot-multi-metric-evaluation-py)- [Visualizations](https://scikit-learn.org/stable/visualizations.htmlvisualizations)- [Real Python Matplotlib guide](https://realpython.com/python-matplotlib-guide/)- [Matplotlib anatomy](https://matplotlib.org/3.2.2/gallery/showcase/anatomy.html)- [Yellowbrick](https://www.slideshare.net/RebeccaBilbro/learning-machine-learning-with-yellowbrick)Let's start by re-running the relevant code from previous lectures, after which we can continue extending our pipelines to experiment with different classifiers.<jupyter_code>import warnings
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.impute import SimpleImputer
from sklearn.model_selection import StratifiedShuffleSplit

#supressing warnings for readability
warnings.filterwarnings("ignore")

# To plot pretty figures directly within Jupyter
%matplotlib inline

# choose your own style: https://matplotlib.org/3.1.0/gallery/style_sheets/style_sheets_reference.html
plt.style.use('seaborn-whitegrid')

# Go to town with https://matplotlib.org/tutorials/introductory/customizing.html
# plt.rcParams.keys()
mpl.rc('axes', labelsize=14, titlesize=14)
mpl.rc('figure', titlesize=20)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

# contants for figsize
S = (8,8)
M = (12,12)
L = (14,14)

# pandas options
pd.set_option("display.max.columns", None)
pd.set_option("display.max.rows", None)
pd.set_option("display.precision", 2)

# read data
df = pd.read_parquet('https://github.com/jads-nl/execute-nhs-proms/blob/master/data/interim/knee-provider.parquet?raw=true')
# handy function to select oks columns
def oks_questions(t='t0'):
  return [
    col for col in df.columns if col.startswith(f"oks_{t}") and not col.endswith("_score")
]

# replace sentinel values in oks columns
# note we are doing imputation on original dataframe (rather than in pipeline later on)
# so we can perform it prior to StratefiedShuffleSplit
oks_no9 = oks_questions('t0') + oks_questions('t1')
impute_oks = SimpleImputer(missing_values=9, strategy="most_frequent")
df.loc[:, oks_no9] = impute_oks.fit_transform(df[oks_no9])

# group columns t0
age_band = ["age_band"]
gender = ["gender"]
age_band_categories = sorted([x for x in df.age_band.unique() if isinstance(x, str)])
comorb = [
    "heart_disease",
    "high_bp",
    "stroke",
    "circulation",
    "lung_disease",
    "diabetes",
    "kidney_disease",
    "nervous_system",
    "liver_disease",
    "cancer",
    "depression",
    "arthritis",
]
boolean = ["t0_assisted", "t0_previous_surgery", "t0_disability"]
eq5d = ["t0_mobility", "t0_self_care", "t0_activity", "t0_discomfort", "t0_anxiety"]
eq_vas = ["t0_eq_vas"]
categorical = ["t0_symptom_period", "t0_previous_surgery", "t0_living_arrangements"]
oks_score = ["oks_t0_score"]

# add number of comorbidities as extra feature
impute_comorb = SimpleImputer(missing_values=9, strategy="constant", fill_value=0)
df.loc[:, comorb] = impute_comorb.fit_transform(df[comorb])
df["n_comorb"] = df.loc[:, comorb].sum()


# define outcome Y
CUT_OFF_PAIN = 4
CUT_OFF_FUNCTIONING = 26

for t in ("t0", "t1"):
    df[f"oks_{t}_pain_total"] = df[f"oks_{t}_pain"] + df[f"oks_{t}_night_pain"]
    df[f"oks_{t}_functioning_total"] = (
        df.loc[:, [col for col in oks_questions(t) if "pain" not in col]]
        .sum(axis=1)
    )
    df[f"y_{t}_pain_good"] = df[f"oks_{t}_pain_total"].apply(
        lambda s: True if s >= CUT_OFF_PAIN else False
    )
    df[f"y_{t}_functioning_good"] = df[f"oks_{t}_functioning_total"].apply(
        lambda s: True if s >= CUT_OFF_FUNCTIONING else False
    )

# define binary outcome parameter
df["y_binary"] = np.logical_and(df.y_t1_pain_good, df.y_t1_functioning_good)

# Only using 1 split for stratefied sampling, more folds are used later on in cross-validation
split = StratifiedShuffleSplit(n_splits=1, test_size=0.3, random_state=42)
for train_index, test_index in split.split(df, df["y_binary"]):
    df_train = df.loc[train_index]
    df_test = df.loc[test_index]
    
y_train_pain_good = df_train.y_t1_pain_good
y_train_pain_good = df_train.y_t1_functioning_good
y_train_binary = df_train.y_binary

y_test_pain_good = df_test.y_t1_pain_good
y_test_pain_good = df_test.y_t1_functioning_good
y_test_binary = df_test.y_binary
y0 = pd.crosstab(df_train.y_t0_pain_good, df_train.y_t0_functioning_good, normalize=True)
y0
y1 = pd.crosstab(df_train.y_t1_pain_good, df_train.y_t1_functioning_good, normalize=True)
y1
# same pipeline as lecture 3
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OrdinalEncoder, OneHotEncoder
from sklearn.compose import ColumnTransformer


# preprocessing pipelines for specific columns
age_band_pipe = Pipeline(
    steps=[
        ("impute", SimpleImputer(missing_values=None, strategy="most_frequent")),
        ("ordinal", OrdinalEncoder(categories=[age_band_categories])),
    ]
)
gender_pipe = Pipeline(
    steps=[
        ("impute", SimpleImputer(missing_values=np.nan, strategy="most_frequent")),
        ("onehot", OneHotEncoder()),
    ]
)

# ColumnTransformer on all included columns.
# Note columns that are not specified are dropped by default
transformers = {
    "age": ("age", age_band_pipe, age_band),
    "gender": ("gender", gender_pipe, gender),
    "comorb": (
        "comorb",
        'passthrough',
        comorb,
    ),
    "categorical": (
        "categorical",
        SimpleImputer(missing_values=9, strategy="most_frequent"),
        boolean + eq5d + categorical,
    ),
    "oks": (
        "oks",
        'passthrough',
        oks_questions('t0'),
    ),
    "eq_vas": ("eqvas", SimpleImputer(missing_values=999, strategy="median"), eq_vas),
}
prep = ColumnTransformer(
    transformers=[v for _, v in transformers.items()])
# list of columns for convenience
# https://stackoverflow.com/questions/54646709/sklearn-pipeline-get-feature-name-after-onehotencode-in-columntransformer
X_train = prep.fit_transform(df_train)
X_test = prep.fit_transform(df_test)
X_columns = pd.Series(
    age_band
    + prep.named_transformers_["gender"]["onehot"].get_feature_names().tolist()
    + comorb
    + boolean
    + eq5d
    + categorical
    + oks_questions()
    + oks_score
    + eq_vas
)
X_columns<jupyter_output><empty_output><jupyter_text>Classification Data mining goalsRecall from our earlier discussions (phase 1: Business Understanding) that we want to predict the outcome of surgery. To make this more explicit, for now we choose to predict the outcome being in the bottom-right quadrant of `y1` i.e. with good functioning and less pain. The base-rate `P / (P + N) = 71%`.  Single estimators: Stochastic Gradient Decent, Decision Tree (CART) and SVM Stochastic Gradient Descent<jupyter_code>%%time
from sklearn.linear_model import SGDClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import cross_val_score


sgd = Pipeline(
    steps=[("prep", prep), ("sgd", SGDClassifier(loss="hinge", penalty="l2", tol=1e-3))]
)

# https://scikit-learn.org/stable/tutorial/statistical_inference/putting_together.html
sgd_parameters = {"sgd__max_iter": [10]} #, 100, 1000]}

sgd_search = GridSearchCV(sgd, sgd_parameters, cv=5)
_ = sgd_search.fit(df_train, y_train_binary)  # assign output of fit to dummy parameter to prevent printing in notebook
confusion_matrix(y_train_binary, sgd_search.best_estimator_.predict(df_train), normalize='all').round(3)
cross_val_score(sgd_search.best_estimator_, df_train, y_train_binary, cv=3, scoring="roc_auc").round(3)<jupyter_output><empty_output><jupyter_text>Decision Tree (CART)<jupyter_code>%%time
from sklearn.tree import DecisionTreeClassifier

cart = Pipeline(
    steps=[
        ("prep", prep),
        ("cart", DecisionTreeClassifier()),
    ]
)

cart_parameters = {
'cart__max_depth': [10, 30],
'cart__min_samples_leaf': [0.02, 0.05, 0.1, 0.2],
}

cart_search = GridSearchCV(cart, cart_parameters, cv=5)
_ = cart_search.fit(df_train, y_train_binary);
confusion_matrix(y_train_binary, cart_search.best_estimator_.predict(df_train), normalize='all').round(3)
cross_val_score(cart_search.best_estimator_, df_train, y_train_binary, cv=3, scoring="roc_auc").round(3)<jupyter_output><empty_output><jupyter_text>LinearSVC<jupyter_code>%%time
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC

svm = Pipeline(
    steps=[
        ("prep", prep),
        ("scaler", StandardScaler()),
        ("clf", LinearSVC(C=1, loss="hinge")),
    ]
)
_ = svm.fit(df_train, y_train_binary);
confusion_matrix(y_train_binary, svm.predict(df_train), normalize='all').round(3)
cross_val_score(svm, df_train, y_train_binary, cv=3, scoring="roc_auc").round(3)<jupyter_output><empty_output><jupyter_text>ROC curveIt is more intuitive to compare different estimators visually, i.e. visualizing the performance metrics. scikit-learn 0.22 comes with a [new plotting API](https://scikit-learn.org/stable/auto_examples/release_highlights/plot_release_highlights_0_22_0.htmlsphx-glr-auto-examples-release-highlights-plot-release-highlights-0-22-0-py), so it's much easier than before to quickly produce charts. We start with the good ol' ROC curve (although later on I will stress not to use this!).<jupyter_code>from sklearn.metrics import plot_roc_curve, plot_precision_recall_curve

clfs = {
    "sgd": sgd_search.best_estimator_,
    "cart": cart_search.best_estimator_,
    "svm": svm,
}

fig_roc, ax_roc = plt.subplots(figsize=S)
plt.plot([0, 1], [0, 1], 'k--')
for k,v in clfs.items():
    plot_roc_curve(v, df_test, y_test_binary, name=k, ax=ax_roc);<jupyter_output><empty_output><jupyter_text>PR curveThe results from the ROC curve are not very exciting. Let's look at the precision recall curve.<jupyter_code>fig_prc, ax_prc = plt.subplots(figsize=S)
plt.plot([0,1], [ y1.iloc[1,1],  y1.iloc[1,1]], 'k--')
for k,v in clfs.items():
    plot_precision_recall_curve(v, df_test, y_test_binary, name=k, ax=ax_prc)<jupyter_output><empty_output><jupyter_text>Ensemble methods: Random Forest & GradientBoosted TreesWe are going to do the same - train classifiers with grid search and cross-validation - but using an ensemble of learners this time.<jupyter_code>%%time
from sklearn.ensemble import RandomForestClassifier


rf = Pipeline(
    steps=[("prep", prep), ("rf", RandomForestClassifier(random_state=42, n_jobs=-1))]
)
rf_parameters = {
    "rf__max_samples": [1000],  # don't use the whole dataset: takes too long!
    "rf__n_estimators": [100],
    "rf__criterion": ["gini"],
    "rf__max_features": ["auto"],
    "rf__min_samples_leaf": [1, 2, 5],
    "rf__max_depth": [10, 15],
    "rf__oob_score": [True, False],
}
rf_search = GridSearchCV(rf, rf_parameters, cv=5)
_ = rf_search.fit(df_train, y_train_binary)
# take very long, so not run

# %%time
# from sklearn.ensemble import GradientBoostingClassifier

# gb = Pipeline(
#     steps=[("prep", prep), ("gb", GradientBoostingClassifier())]
# )
# gb_parameters = {
#     "gb__learning_rate": [1.0],
#     "gb__n_estimators": [100, 200],
#     "gb__min_samples_leaf": [1, 2, 5],
#     "gb__max_features": ["auto"],
#     "gb__max_depth": [2, 3, 5],
# }
# gb_search = GridSearchCV(gb, gb_parameters, cv=5)
# _ = gb_search.fit(df_train, y_train_binary)<jupyter_output><empty_output><jupyter_text>XGBoostThe GradientBoostingClassifier in sklearn is quite slow. [XGboost](https://xgboost.readthedocs.io/en/latest/index.html) is a better and faster implementation. After installing it separately, you can use the sklearn API to use it in our existing pipeline.<jupyter_code>%%time
import xgboost

xgb = Pipeline(steps=[("prep", prep), ("xgb", xgboost.XGBClassifier(random_state=42))])
_ = xgb.fit(df_train, y_train_binary)
confusion_matrix(y_train_binary, xgb.predict(df_train), normalize='all').round(3)
ensembles = {
    "rf": rf_search.best_estimator_,
#     "gb": gb_search.best_estimator_,
    "xgb": xgb,
}

fig_roc2, ax_roc2 = plt.subplots(figsize=S)
plt.plot([0, 1], [0, 1], 'k--')
for k,v in ensembles.items():
    plot_roc_curve(v, df_train, y_train_binary, name=k, ax=ax_roc2)
fig_prc2, ax_prc2 = plt.subplots(figsize=S)
plt.plot([0, 1], [y1.iloc[1,1], y1.iloc[1,1]], 'k--')
for k,v in ensembles.items():
    plot_precision_recall_curve(v, df_train, y_train_binary, name=k, ax=ax_prc2)<jupyter_output><empty_output><jupyter_text>But hang on, don't get too excited yet: aren't we overfitting?<jupyter_code>fig_prc2a, ax_prc2a = plt.subplots(figsize=S)
plt.plot([0, 1], [y1.iloc[1,1], y1.iloc[1,1]], 'k--')
for k,v in ensembles.items():
    plot_precision_recall_curve(v, df_test, y_test_binary, name=k, ax=ax_prc2a)<jupyter_output><empty_output><jupyter_text>Ensemble methods: balancedThe results are better than single estimators, but can we make it better still? We will try balanced Random Forest, using the `class_weight` option in scikit-learn and comparing that next to the `BalancedRandomForestClassifier` from the [imbalanced-learn](https://imbalanced-learn.readthedocs.io/en/stable/) package.<jupyter_code>%%time
from sklearn.ensemble import RandomForestClassifier


brf = Pipeline(
    steps=[("prep", prep), ("rf", RandomForestClassifier(random_state=42, n_jobs=-1))]
)
brf_parameters = {
    "rf__max_samples": [1000],  # don't use the whole dataset: takes too long!
    "rf__n_estimators": [100],
    "rf__criterion": ["gini"],
    "rf__max_features": ["auto"],
    "rf__min_samples_leaf": [1, 2, 5],
    "rf__max_depth": [10, 15],
    "rf__oob_score": [True, False],
    "rf__class_weight": ["balanced"],
}
brf_search = GridSearchCV(brf, brf_parameters, cv=5)
_ = brf_search.fit(df_train, y_train_binary)
%%time
from imblearn.ensemble import BalancedRandomForestClassifier


ibrf = Pipeline(
    steps=[("prep", prep), ("rf", BalancedRandomForestClassifier(random_state=42, n_jobs=-1))]
)
ibrf_parameters = {
    "rf__max_samples": [1000],  # don't use the whole dataset: takes too long!
    "rf__n_estimators": [100],
    "rf__criterion": ["gini"],
    "rf__max_features": ["auto"],
    "rf__min_samples_leaf": [1, 2, 5],
    "rf__max_depth": [10, 15],
    "rf__oob_score": [True, False],
}
ibrf_search = GridSearchCV(ibrf, ibrf_parameters, cv=5)
_ = ibrf_search.fit(df_train, y_train_binary)
ensembles["brf"] = brf_search.best_estimator_
ensembles["ibrf"] = ibrf_search.best_estimator_

fig_prc3, ax_prc3 = plt.subplots(figsize=S)
plt.plot([0, 1], [y1.iloc[1,1], y1.iloc[1,1]], 'k--')
for k,v in ensembles.items():
    plot_precision_recall_curve(v, df_train, y_train_binary, name=k, ax=ax_prc3);<jupyter_output><empty_output><jupyter_text>Try the inverse: predicting patients with not in good-good quadrant<jupyter_code>%%time
# see what performance is for outlier detection
y_train_binary_inverted = np.invert(y_train_binary)

xgb_inv = Pipeline(steps=[("prep", prep), ("xgb", xgboost.XGBClassifier())])
_ = xgb_inv.fit(df_train, y_train_binary_inverted)
fig_prc4, ax_prc4 = plt.subplots(figsize=S)
plt.plot([0, 1], [1 - y1.iloc[1,1], 1 - y1.iloc[1,1]], 'k--')
plot_precision_recall_curve(xgb_inv, df_train, y_train_binary_inverted, name="xgb_inv", ax=ax_prc4);<jupyter_output><empty_output><jupyter_text>Clearly more difficult to detect 'outliers' i.e. patients with less than optimal outcome.  Bonus: YellowbrickTo further optimze your workflow, you can use [Yellowbrick](https://www.scikit-yb.org/en/latest/index.html) for visualizing results.<jupyter_code>from yellowbrick.classifier import ROCAUC, PrecisionRecallCurve

roc = ROCAUC(rf_search,)
roc.fit(df_train, y_train_binary)
roc.score(df_train, y_train_binary)
roc.show();
prc = PrecisionRecallCurve(xgb_inv,)
prc.fit(df_train, y_train_binary_inverted)
prc.score(df_train, y_train_binary_inverted)
prc.show();<jupyter_output><empty_output><jupyter_text>Discussion: model assessment and evaluation Greediness of algorithmsEven if you have a good result, i.e. a good score with the test set, you should beware of the greediness of machine learning algorithms: the solution is _at best_ a local optimum, and sometimes not even that. Depending on the random state at initiazation, the results will vary. To illustrate this, let's do 5 runs of XBGBoost (being the best model we have looked at) with different random states.<jupyter_code>runs = [{"random_state": n} for n in range(2, 52, 10)]
for run in runs:
    xgb = Pipeline(steps=[("prep", prep), ("xgb", xgboost.XGBClassifier(random_state=run['random_state'], subsample=0.7))])
    run["model"] = xgb.fit(df_train, y_train_binary)
fig_random, ax_random = plt.subplots(figsize=S)
ax_random.set_ylim(0.6, 1.0)
plt.plot([0, 1], [y1.iloc[1,1], y1.iloc[1,1]], 'k--')
for run in runs:
    plot_precision_recall_curve(run['model'], df_test, y_test_binary, name=k, ax=ax_random, alpha=0.5)<jupyter_output><empty_output><jupyter_text>Feature importanceMany tree-based methods can compute feature importance. Also here, note that this is a relative importance and that due to the greediness you can not draw very hard conclusions (features by swap places depending on the parameter settings). It is useful to get an intuitive feel which features contain the most information for your classification task.<jupyter_code>fig_importance, ax_importance = plt.subplots(figsize=S)
xgboost.plot_importance(xgb['xgb'], ax=ax_importance);
# To DO: fix feature names in plot
list(zip(xgb['xgb'].get_booster().feature_names, X_columns))<jupyter_output><empty_output><jupyter_text>Bias-variance trade-offOne of the nice features of XGBoost is that it has [explicitly defined the complexity function](https://xgboost.readthedocs.io/en/latest/tutorials/model.html):$$\Omega(f) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2$$That means that we can plot the results from the gridsearch into a prediction error vs. complexity graph to visualize the bias-variance trade-off.<jupyter_code>#TODO: add prediction error vs. complexity plot<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12434265493684572, 'fraction_numerical': 0.01725070034894579, 'mean_word_length': 5.7489220563847425, 'pattern_counts': {'":': 41, '<': 61, '<?xml version=': 0, '>': 63, 'https://': 21, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '317827', 'license': 'MIT', 'n_tokens_mistral': 7428, 'n_tokens_neox': 7059, 'n_words': 1672, 'path': 'notebooks/4.0-modeling-classification.ipynb', 'repo_name': 'jads-nl/nhs-proms', 'size': '20347'}