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
%%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)
%%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)
%%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)
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);
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)
%%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)
%%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)
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)
%%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);
%%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);
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();
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)
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))
#TODO: add prediction error vs. complexity plot
