{"cells":[{"metadata":{},"cell_type":"markdown","source":"\n<h1><center>APTOS 2019 Blindness Detection</center></h1>\n<h2><center>Diabetic retinopathy - SHAP model explainability</center></h2>\n\n\nIn this work, I'll train a baseline ResNet50, evaluate the model, and use SHAP model explainability technique to help us better understand our model's predictions, and how we could further improve its performance.\n\n#### About [SHAP](https://github.com/slundberg/shap) from its source:\n\n<img src=\"https://raw.githubusercontent.com/slundberg/shap/master/docs/artwork/shap_diagram.png\" width=\"400\">\n\n##### SHAP (SHapley Additive exPlanations) is a unified approach to explain the output of any machine learning model. SHAP connects game theory with local explanations, uniting several previous methods [1-7] and representing the only possible consistent and locally accurate additive feature attribution method based on expectations (see our [papers](https://github.com/slundberg/shap#citations) for details)."},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"## Dependencies"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import os\nimport shap\nimport random\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom sklearn.utils import class_weight\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import confusion_matrix, cohen_kappa_score\nfrom keras.models import Model\nfrom keras import optimizers, applications\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.callbacks import EarlyStopping, ReduceLROnPlateau\nfrom keras.layers import Dense, Dropout, GlobalAveragePooling2D, Input\n\n# Set seeds to make the experiment more reproducible.\nfrom tensorflow import set_random_seed\ndef seed_everything(seed=0):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    set_random_seed(seed)\n\nseed = 0\nseed_everything(seed)\n\n%matplotlib inline\nsns.set(style=\"whitegrid\")\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true},"cell_type":"markdown","source":"## Load data"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_kg_hide-input":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/aptos2019-blindness-detection/train.csv')\ntest = pd.read_csv('../input/aptos2019-blindness-detection/test.csv')\nprint('Number of train samples: ', train.shape[0])\nprint('Number of test samples: ', test.shape[0])\n\n# Preprocecss data\ntrain[\"id_code\"] = train[\"id_code\"].apply(lambda x: x + \".png\")\ntest[\"id_code\"] = test[\"id_code\"].apply(lambda x: x + \".png\")\ntrain['diagnosis'] = train['diagnosis'].astype('str')\ndisplay(train.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Model parameters\nBATCH_SIZE = 8\nEPOCHS = 40\nWARMUP_EPOCHS = 2\nLEARNING_RATE = 1e-4\nWARMUP_LEARNING_RATE = 1e-3\nHEIGHT = 320\nWIDTH = 320\nCANAL = 3\nN_CLASSES = train['diagnosis'].nunique()\nES_PATIENCE = 5\nRLROP_PATIENCE = 3\nDECAY_DROP = 0.5","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train test split"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_val = train_test_split(train, test_size=0.2, random_state=seed)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data generator"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"train_datagen=ImageDataGenerator(rescale=1./255, \n                                 rotation_range=360,\n                                 horizontal_flip=True,\n                                 vertical_flip=True)\n\ntrain_generator=train_datagen.flow_from_dataframe(\n    dataframe=X_train,\n    directory=\"../input/aptos2019-blindness-detection/train_images/\",\n    x_col=\"id_code\",\n    y_col=\"diagnosis\",\n    class_mode=\"categorical\",\n    batch_size=BATCH_SIZE,\n    target_size=(HEIGHT, WIDTH),\n    seed=0)\n\nvalidation_datagen = ImageDataGenerator(rescale=1./255)\n\nvalid_generator=validation_datagen.flow_from_dataframe(\n    dataframe=X_val,\n    directory=\"../input/aptos2019-blindness-detection/train_images/\",\n    x_col=\"id_code\",\n    y_col=\"diagnosis\",\n    class_mode=\"categorical\", \n    batch_size=BATCH_SIZE,   \n    target_size=(HEIGHT, WIDTH),\n    seed=0)\n\ntest_datagen = ImageDataGenerator(rescale=1./255)\n\ntest_generator = test_datagen.flow_from_dataframe(  \n        dataframe=test,\n        directory = \"../input/aptos2019-blindness-detection/test_images/\",\n        x_col=\"id_code\",\n        batch_size=1,\n        class_mode=None,\n        shuffle=False,\n        target_size=(HEIGHT, WIDTH),\n        seed=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_model(input_shape, n_out):\n    input_tensor = Input(shape=input_shape)\n    base_model = applications.ResNet50(weights=None, \n                                       include_top=False,\n                                       input_tensor=input_tensor)\n    base_model.load_weights('../input/resnet50/resnet50_weights_tf_dim_ordering_tf_kernels_notop.h5')\n\n    x = GlobalAveragePooling2D()(base_model.output)\n    x = Dropout(0.5)(x)\n    x = Dense(2048, activation='relu')(x)\n    x = Dropout(0.5)(x)\n    final_output = Dense(n_out, activation='softmax', name='final_output')(x)\n    model = Model(input_tensor, final_output)\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train top layers"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"model = create_model(input_shape=(HEIGHT, WIDTH, CANAL), n_out=N_CLASSES)\n\nfor layer in model.layers:\n    layer.trainable = False\n\nfor i in range(-5, 0):\n    model.layers[i].trainable = True\n    \nclass_weights = class_weight.compute_class_weight('balanced', np.unique(train['diagnosis'].astype('int').values), train['diagnosis'].astype('int').values)\n\nmetric_list = [\"accuracy\"]\noptimizer = optimizers.Adam(lr=WARMUP_LEARNING_RATE)\nmodel.compile(optimizer=optimizer, loss='categorical_crossentropy',  metrics=metric_list)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"STEP_SIZE_TRAIN = train_generator.n//train_generator.batch_size\nSTEP_SIZE_VALID = valid_generator.n//valid_generator.batch_size\n\nhistory_warmup = model.fit_generator(generator=train_generator,\n                                     steps_per_epoch=STEP_SIZE_TRAIN,\n                                     validation_data=valid_generator,\n                                     validation_steps=STEP_SIZE_VALID,\n                                     epochs=WARMUP_EPOCHS,\n                                     class_weight=class_weights,\n                                     verbose=1).history","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fine-tune the complete model"},{"metadata":{"_kg_hide-input":false,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"for layer in model.layers:\n    layer.trainable = True\n\nes = EarlyStopping(monitor='val_loss', mode='min', patience=ES_PATIENCE, restore_best_weights=True, verbose=1)\nrlrop = ReduceLROnPlateau(monitor='val_loss', mode='min', patience=RLROP_PATIENCE, factor=DECAY_DROP, min_lr=1e-6, verbose=1)\n\ncallback_list = [es, rlrop]\noptimizer = optimizers.Adam(lr=LEARNING_RATE)\nmodel.compile(optimizer=optimizer, loss='categorical_crossentropy',  metrics=metric_list)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"history_finetunning = model.fit_generator(generator=train_generator,\n                                          steps_per_epoch=STEP_SIZE_TRAIN,\n                                          validation_data=valid_generator,\n                                          validation_steps=STEP_SIZE_VALID,\n                                          epochs=EPOCHS,\n                                          callbacks=callback_list,\n                                          class_weight=class_weights,\n                                          verbose=1).history","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model loss graph "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"history = {'loss': history_warmup['loss'] + history_finetunning['loss'], \n           'val_loss': history_warmup['val_loss'] + history_finetunning['val_loss'], \n           'acc': history_warmup['acc'] + history_finetunning['acc'], \n           'val_acc': history_warmup['val_acc'] + history_finetunning['val_acc']}\n\nsns.set_style(\"whitegrid\")\nfig, (ax1, ax2) = plt.subplots(2, 1, sharex='col', figsize=(20, 14))\n\nax1.plot(history['loss'], label='Train loss')\nax1.plot(history['val_loss'], label='Validation loss')\nax1.legend(loc='best')\nax1.set_title('Loss')\n\nax2.plot(history['acc'], label='Train accuracy')\nax2.plot(history['val_acc'], label='Validation accuracy')\nax2.legend(loc='best')\nax2.set_title('Accuracy')\n\nplt.xlabel('Epochs')\nsns.despine()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model Evaluation\n\n## Confusion Matrix"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# Create empty arays to keep the predictions and labels\nlastFullTrainPred = np.empty((0, N_CLASSES))\nlastFullTrainLabels = np.empty((0, N_CLASSES))\nlastFullValPred = np.empty((0, N_CLASSES))\nlastFullValLabels = np.empty((0, N_CLASSES))\n\n# Add train predictions and labels\nfor i in range(STEP_SIZE_TRAIN+1):\n    im, lbl = next(train_generator)\n    scores = model.predict(im, batch_size=train_generator.batch_size)\n    lastFullTrainPred = np.append(lastFullTrainPred, scores, axis=0)\n    lastFullTrainLabels = np.append(lastFullTrainLabels, lbl, axis=0)\n\n# Add validation predictions and labels\nfor i in range(STEP_SIZE_VALID+1):\n    im, lbl = next(valid_generator)\n    scores = model.predict(im, batch_size=valid_generator.batch_size)\n    lastFullValPred = np.append(lastFullValPred, scores, axis=0)\n    lastFullValLabels = np.append(lastFullValLabels, lbl, axis=0)\n    \n    \nlastFullComPred = np.concatenate((lastFullTrainPred, lastFullValPred))\nlastFullComLabels = np.concatenate((lastFullTrainLabels, lastFullValLabels))\ncomplete_labels = [np.argmax(label) for label in lastFullComLabels]\n\ntrain_preds = [np.argmax(pred) for pred in lastFullTrainPred]\ntrain_labels = [np.argmax(label) for label in lastFullTrainLabels]\nvalidation_preds = [np.argmax(pred) for pred in lastFullValPred]\nvalidation_labels = [np.argmax(label) for label in lastFullValLabels]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, sharex='col', figsize=(24, 7))\nlabels = ['0 - No DR', '1 - Mild', '2 - Moderate', '3 - Severe', '4 - Proliferative DR']\ntrain_cnf_matrix = confusion_matrix(train_labels, train_preds)\nvalidation_cnf_matrix = confusion_matrix(validation_labels, validation_preds)\n\ntrain_cnf_matrix_norm = train_cnf_matrix.astype('float') / train_cnf_matrix.sum(axis=1)[:, np.newaxis]\nvalidation_cnf_matrix_norm = validation_cnf_matrix.astype('float') / validation_cnf_matrix.sum(axis=1)[:, np.newaxis]\n\ntrain_df_cm = pd.DataFrame(train_cnf_matrix_norm, index=labels, columns=labels)\nvalidation_df_cm = pd.DataFrame(validation_cnf_matrix_norm, index=labels, columns=labels)\n\nsns.heatmap(train_df_cm, annot=True, fmt='.2f', cmap=\"Blues\", ax=ax1).set_title('Train')\nsns.heatmap(validation_df_cm, annot=True, fmt='.2f', cmap=sns.cubehelix_palette(8), ax=ax2).set_title('Validation')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Quadratic Weighted Kappa"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Train Cohen Kappa score: %.3f\" % cohen_kappa_score(train_preds,train_labels, weights='quadratic'))\nprint(\"Validation Cohen Kappa score: %.3f\" % cohen_kappa_score(validation_preds, validation_labels, weights='quadratic'))\nprint(\"Complete set Cohen Kappa score: %.3f\" % cohen_kappa_score(train_preds+validation_preds, train_labels+validation_labels, weights='quadratic'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# SHAP Model explainability\n\n#### About SHAP's DeepExplainer from the [source repository](https://github.com/slundberg/shap#deep-learning-example-with-deepexplainer-tensorflowkeras-models): \n- Deep SHAP is a high-speed approximation algorithm for SHAP values in deep learning models that builds on a connection with [DeepLIFT](https://arxiv.org/abs/1704.02685) described in the SHAP NIPS paper. The implementation here differs from the original DeepLIFT by using a distribution of background samples instead of a single reference value, and using Shapley equations to linearize components such as max, softmax, products, divisions, etc.\n\n### First let's see the images that we will explain"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"n_explain = 2\nvalid_generator.batch_size = 10 # background dataset\nbackground, lbls = next(valid_generator)\n\nsns.set_style(\"white\")\nplt.figure(figsize=[8, 8])\nfor index, image in enumerate(background[:n_explain]):\n    plt.subplot(n_explain, 1, index+1)\n    plt.imshow(image)\n    plt.title(\"Image %s, Label: %s\" % (index, np.argmax(lbls[index])))\n    \nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Now the SHAP explanation"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# explain predictions of the model on \"n_explain\" images\ne = shap.DeepExplainer(model, background)\nshap_values = e.shap_values(background)\n\n# plot the feature attributions\nshap.image_plot(shap_values, -background[:n_explain], labels=lbls, hspace=0.1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- The plot above explains five outputs (our five levels of diabetic retinopathy 0-5) for three different images. Red pixels increase the model's output while blue pixels decrease the output. The input images are shown on the left (they are black because most of the pixels are greater than 0), and as nearly transparent grayscale backings behind each of the explanations. The sum of the SHAP values equals the difference between the expected model output (averaged over the background dataset, here I'm using 10 images) and the current model output. \n- Note that for the images that the label is \"1.0\" (the correct one), we a greater pink area.\n- Labels that have as much pink area as the correct one are labels that our model probably doesn't have a high confidence prediction."},{"metadata":{},"cell_type":"markdown","source":"## Let's try on a few more images"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"n_explain = 3\nbackground, lbls = next(valid_generator)\n\nsns.set_style(\"white\")\nplt.figure(figsize=[12, 12])\nfor index, image in enumerate(background[:n_explain]):\n    plt.subplot(n_explain, 1, index+1)\n    plt.imshow(image)\n    plt.title(\"Image %s, Label: %s\" % (index, np.argmax(lbls[index])))\n    \nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# explain predictions of the model on \"n_explain\" images\ne = shap.DeepExplainer(model, background)\nshap_values = e.shap_values(background)\n\n# plot the feature attributions\nshap.image_plot(shap_values, -background[:n_explain], labels=lbls, hspace=0.1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Apply model to test set and output predictions"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"test_generator.reset()\nSTEP_SIZE_TEST = test_generator.n//test_generator.batch_size\npreds = model.predict_generator(test_generator, steps=STEP_SIZE_TEST)\npredictions = [np.argmax(pred) for pred in preds]\n\nfilenames = test_generator.filenames\nresults = pd.DataFrame({'id_code':filenames, 'diagnosis':predictions})\nresults['id_code'] = results['id_code'].map(lambda x: str(x)[:-4])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predictions class distribution"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.subplots(1, 1, sharex='col', figsize=(24, 8.7))\nsns.countplot(x=\"diagnosis\", data=results, palette=\"GnBu_d\")\nsns.despine()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"results.to_csv('submission.csv', index=False)\nresults.head(10)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}