<jupyter_start><jupyter_text>Logistic Regression Classifier - Fitting to all stroke teams together Including shapley value analysis This notebook is based on the original notebook "logistic_regression_single_fit.ipynb", which has these aims. AimsAssess accuracy of a logistic regression classifier, using k-fold (5-fold) training/test data splits (each data point is present in one and only one of the five test sets). This notebook fits all data in a single model, with hospital ID being a one-hot encoded feature.The notebook includes:* A range of accuracy scores* Receiver operating characteristic (ROC) and Sensitivity-Specificity Curves * Identify feature weights (model coefficients)* Performing a learning rate test (relationship between training set size and accuracy) This notebook extends "logistic_regression_single_fit.ipynb" to include Shapley value analysis (comparing them to feature weights, and creating beeswarm and waterfall plots). Aims extension:1. Calculate Shapley values for the logistic regression modelShapley values are calculated per instance, per feature. To get a value per feature take the mean of all the instances. Do this for each k-fold.    1. Get average and standard deviation of Shapley values and model weights for each of the five k-fold splitsFor shapley values this becomes the mean of the 5 means (one for each of the 5 kfolds)    1. Plot mean Shapley value vs mean weight - and put an error bar of standard error (SD / root(n)) for each1. Show the top 10 (absolute value) Shapley values and model weights1. Combine all the data (so no train/test split) - refit and get Shapley values and model weights, and show a Beeswarm plot for all the data, and then show a Waterfall plot for three examples (pref in a single wide MatPlotLib plot) - one with a low probability of giving thrombolysis, one with a probability close to 0.5, and one with a high probability of giving thrombolysis.Beeswarm uses all the data, waterfall plot shows individual instances (so choose 3 instances).Using code from https://medium.com/dataman-in-ai/the-shap-with-more-elegant-charts-bc3e73fa1c0c to create the SHAP plots, and transform the log-odds to probabilities for Waterfall plot.Things I learnt: 1. SHAP plots can display either the original feature value, or the standardised value, depending on which argument you pass1. dataset needs to be a pandas DataFrame for the plot to include the feature names. The function "standardise_data()" returns a numpy array - hence why the plots weren't displaying feature names.1. Can not work out how to plot waterfall plots in subplots, in order to display them side-by-side. Extension task 11 Jan 2022Where you have the section on 'fit whole data' - just use the first train/test spilt. Can you then compare average Shap values in train and test (should eb similar!) and then pick three examples (like you did before) from the test set - so we are showing how Shap explains examples not used in training.P.S. It was interesting that time to scan was dominant on average, but was not always the most important feature.Could you compare mean Shap values with median values please.Also when calculating the Shapley values we took absolute of the mean. Can we take the mean of the absolute.  Import libraries<jupyter_code># Turn warnings off to keep notebook tidy
import warnings
warnings.filterwarnings("ignore")

import os
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import numpy as np
import pandas as pd

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import auc
from sklearn.metrics import roc_curve<jupyter_output><empty_output><jupyter_text>Install shap if necessary (e.g. running in Google Colab)<jupyter_code>try:
    import shap
except:
    !pip install shap
    import shap

# Install math for sqrt
import math<jupyter_output><empty_output><jupyter_text>Import dataData has previously been split into 5 training/test splits.<jupyter_code>data_loc = '../../data/sam_1/kfold_5fold/'
train_data, test_data = [], []

for i in range(5):
    
    train_data.append(pd.read_csv(data_loc + 'train_{0}.csv'.format(i)))
    test_data.append(pd.read_csv(data_loc + 'test_{0}.csv'.format(i)))
train_data[0]<jupyter_output><empty_output><jupyter_text>Functions   Standardise dataStandardisation subtracts the mean and divides by the standard deviation, for each feature. Here we use the sklearn built-in method for standardisation.<jupyter_code>def standardise_data(X_train, X_test):
    """
    Converts all data to a similar scale.
    Standardisation subtracts mean and divides by standard deviation
    for each feature.
    Standardised data will have a mena of 0 and standard deviation of 1.
    The training data mean and standard deviation is used to standardise both
    training and test set data.
    """
    
    # Initialise a new scaling object for normalising input data
    sc = StandardScaler() 

    # Set up the scaler just on the training set
    sc.fit(X_train)

    # Apply the scaler to the training and test sets
    train_std=sc.transform(X_train)
    test_std=sc.transform(X_test)
    
    return train_std, test_std<jupyter_output><empty_output><jupyter_text>Calculate accuracy measures<jupyter_code>def calculate_accuracy(observed, predicted):
    
    """
    Calculates a range of accuracy scores from observed and predicted classes.
    
    Takes two list or NumPy arrays (observed class values, and predicted class 
    values), and returns a dictionary of results.
    
     1) observed positive rate: proportion of observed cases that are +ve
     2) Predicted positive rate: proportion of predicted cases that are +ve
     3) observed negative rate: proportion of observed cases that are -ve
     4) Predicted negative rate: proportion of predicted cases that are -ve  
     5) accuracy: proportion of predicted results that are correct    
     6) precision: proportion of predicted +ve that are correct
     7) recall: proportion of true +ve correctly identified
     8) f1: harmonic mean of precision and recall
     9) sensitivity: Same as recall
    10) specificity: Proportion of true -ve identified:        
    11) positive likelihood: increased probability of true +ve if test +ve
    12) negative likelihood: reduced probability of true +ve if test -ve
    13) false positive rate: proportion of false +ves in true -ve patients
    14) false negative rate: proportion of false -ves in true +ve patients
    15) true positive rate: Same as recall
    16) true negative rate
    17) positive predictive value: chance of true +ve if test +ve
    18) negative predictive value: chance of true -ve if test -ve
    
    """
    
    # Converts list to NumPy arrays
    if type(observed) == list:
        observed = np.array(observed)
    if type(predicted) == list:
        predicted = np.array(predicted)
    
    # Calculate accuracy scores
    observed_positives = observed == 1
    observed_negatives = observed == 0
    predicted_positives = predicted == 1
    predicted_negatives = predicted == 0
    
    true_positives = (predicted_positives == 1) & (observed_positives == 1)
    
    false_positives = (predicted_positives == 1) & (observed_positives == 0)
    
    true_negatives = (predicted_negatives == 1) & (observed_negatives == 1)
    
    false_negatives = (predicted_negatives == 1) & (observed_negatives == 0)
    
    accuracy = np.mean(predicted == observed)
    
    precision = (np.sum(true_positives) /
                 (np.sum(true_positives) + np.sum(false_positives)))
        
    recall = np.sum(true_positives) / np.sum(observed_positives)
    
    sensitivity = recall
    
    f1 = 2 * ((precision * recall) / (precision + recall))
    
    specificity = np.sum(true_negatives) / np.sum(observed_negatives)
    
    positive_likelihood = sensitivity / (1 - specificity)
    
    negative_likelihood = (1 - sensitivity) / specificity
    
    false_positive_rate = 1 - specificity
    
    false_negative_rate = 1 - sensitivity
    
    true_positive_rate = sensitivity
    
    true_negative_rate = specificity
    
    positive_predictive_value = (np.sum(true_positives) / 
                                 np.sum(observed_positives))
    
    negative_predictive_value = (np.sum(true_negatives) / 
                                  np.sum(observed_negatives))
    
    # Create dictionary for results, and add results
    results = dict()
    
    results['observed_positive_rate'] = np.mean(observed_positives)
    results['observed_negative_rate'] = np.mean(observed_negatives)
    results['predicted_positive_rate'] = np.mean(predicted_positives)
    results['predicted_negative_rate'] = np.mean(predicted_negatives)
    results['accuracy'] = accuracy
    results['precision'] = precision
    results['recall'] = recall
    results['f1'] = f1
    results['sensitivity'] = sensitivity
    results['specificity'] = specificity
    results['positive_likelihood'] = positive_likelihood
    results['negative_likelihood'] = negative_likelihood
    results['false_positive_rate'] = false_positive_rate
    results['false_negative_rate'] = false_negative_rate
    results['true_positive_rate'] = true_positive_rate
    results['true_negative_rate'] = true_negative_rate
    results['positive_predictive_value'] = positive_predictive_value
    results['negative_predictive_value'] = negative_predictive_value
    
    return results<jupyter_output><empty_output><jupyter_text>Find model probability threshold to match predicted and actual thrombolysis use<jupyter_code>def find_threshold(probabilities, true_rate):
    
    """
    Find classification threshold to calibrate model
    """
    
    index = (1-true_rate)*len(probabilities)
    
    threshold = sorted(probabilities)[int(index)]
    
    return threshold<jupyter_output><empty_output><jupyter_text>Line intersectUsed to find point of sensitivity-specificity curve where sensitivity = specificity.<jupyter_code>def get_intersect(a1, a2, b1, b2):
    """ 
    Returns the point of intersection of the lines passing through a2,a1 and b2,b1.
    a1: [x, y] a point on the first line
    a2: [x, y] another point on the first line
    b1: [x, y] a point on the second line
    b2: [x, y] another point on the second line
    """
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)<jupyter_output><empty_output><jupyter_text>Create new folderIf a folder does not exist, create it.This is used when saving a file to a folder<jupyter_code>def create_folder_if_not_exist(output_path):

    if not os.path.exists(output_path):
        os.mkdir(output_path)
    
    return()<jupyter_output><empty_output><jupyter_text>Getting Shapley values with SHAPWe use the `explainer` and `shap_values` methods from the SHAP library to get Shapley values.A slight confusion is that we can get Shapley values alone, or with other variables.<jupyter_code>def calculate_shapley_values(model, X_train):
    # We can use the `explainer(X)` method to get objects with shap values, base value*, and feature data.
    # *The base value is “the value that would be predicted if we did not know any features for the current output.” In other words, it is the mean prediction, or mean(yhat).

    # Train explainer on Training set
    explainer = shap.Explainer(model, X_train)
    
    # To get Shapley values alone:
    shap_values = explainer.shap_values(X_train)

    # Get Shapley values along with base, and features
    shap_values_extended = explainer(X_train)
    
    return (shap_values, shap_values_extended)<jupyter_output><empty_output><jupyter_text>Part 1: Fit logistic regression model using training and test sets (use k-fold splits)<jupyter_code># Set up list to store models and calibarion threshold
single_models = []
thresholds_calibrate = []

# Set up lists for observed and predicted
observed = []
predicted_proba = []
predicted = []
shap_values_train = []
shap_values_ext_train = []
shap_values_test = []
shap_values_ext_test = []

# Set up list for feature weightsd
feature_weights = []

# Loop through k folds
for k_fold in range(5):
    
    # Get k fold split
    train = train_data[k_fold]
    test = test_data[k_fold]
    
    # Get X and y
    X_train = train.drop('S2Thrombolysis', axis=1)
    X_test = test.drop('S2Thrombolysis', axis=1)
    y_train = train['S2Thrombolysis']
    y_test = test['S2Thrombolysis']
    
    # One hot encode hospitals
    X_train_hosp = pd.get_dummies(X_train['StrokeTeam'], prefix = 'team')
    X_train = pd.concat([X_train, X_train_hosp], axis=1)
    X_train.drop('StrokeTeam', axis=1, inplace=True)
    X_test_hosp = pd.get_dummies(X_test['StrokeTeam'], prefix = 'team')
    X_test = pd.concat([X_test, X_test_hosp], axis=1)
    X_test.drop('StrokeTeam', axis=1, inplace=True)
    
    # Standardise X data
    X_train_std, X_test_std = standardise_data(X_train, X_test)

    # Define and Fit model
    model = LogisticRegression(solver='lbfgs')
    model.fit(X_train_std, y_train)
    
    # Get feature weights
    weights = model.coef_[0]
    feature_weights.append(weights)
    
    # Get predicted probabilities
    y_probs = model.predict_proba(X_test_std)[:,1]
    observed.append(y_test)
    predicted_proba.append(y_probs)
    
    # Calibrate model and get class
    true_rate = np.mean(y_test)
    threshold = find_threshold(y_probs, true_rate)
    thresholds_calibrate.append(threshold)
    y_class = y_probs >= threshold
    y_class = np.array(y_class) * 1.0
    predicted.append(y_class)
    
    # Calculate Shapley values
    ### Global influence of features in the training set
    ### Here we are looking at Shapley values for the training set, to understand the influence of features in training.
    (shapley_values, shapley_values_extended) = calculate_shapley_values(model, X_train_std)
    shap_values_train.append(shapley_values)
    shap_values_ext_train.append(shapley_values_extended)

    ### And for the Test set
    (shapley_values, shapley_values_extended) = calculate_shapley_values(model, X_test_std)
    shap_values_test.append(shapley_values)
    shap_values_ext_test.append(shapley_values_extended)
    
    # Print accuracy
    accuracy = np.mean(y_class == y_test)
    print (
        f'Run {k_fold}, accuracy: {accuracy:0.3f}, threshold {threshold:0.3f}')<jupyter_output>Run 0, accuracy: 0.833, threshold 0.474
Run 1, accuracy: 0.834, threshold 0.475
Run 2, accuracy: 0.830, threshold 0.469
Run 3, accuracy: 0.833, threshold 0.470
Run 4, accuracy: 0.833, threshold 0.470<jupyter_text>Results   Accuracy measures<jupyter_code># Set up list for results
k_fold_results = []

# Loop through k fold predictions and get accuracy measures
for i in range(5):
    results = calculate_accuracy(observed[i], predicted[i])
    k_fold_results.append(results)
    
# Put results in DataFrame
single_fit_results = pd.DataFrame(k_fold_results).T
single_fit_results
single_fit_results.T.describe()<jupyter_output><empty_output><jupyter_text>Receiver Operator Characteristic and Sensitivity-Specificity Curves Receiver Operator Characteristic Curve:<jupyter_code># Set up lists for results
k_fold_fpr = [] # false positive rate
k_fold_tpr = [] # true positive rate
k_fold_thresholds = [] # threshold applied
k_fold_auc = [] # area under curve

# Loop through k fold predictions and get ROC results 
for i in range(5):
    fpr, tpr, thresholds = roc_curve(observed[i], predicted_proba[i])
    roc_auc = auc(fpr, tpr)
    k_fold_fpr.append(fpr)
    k_fold_tpr.append(tpr)
    k_fold_thresholds.append(thresholds)
    k_fold_auc.append(roc_auc)
    print (f'Run {i} AUC {roc_auc:0.4f}')

# Show mean area under curve  
mean_auc = np.mean(k_fold_auc)
sd_auc = np.std(k_fold_auc)
print (f'\nMean AUC: {mean_auc:0.4f}')
print (f'SD AUC: {sd_auc:0.4f}')<jupyter_output>Run 0 AUC 0.9035
Run 1 AUC 0.9056
Run 2 AUC 0.9014
Run 3 AUC 0.9036
Run 4 AUC 0.9045

Mean AUC: 0.9037
SD AUC: 0.0014<jupyter_text>Sensitivity-specificity curve:<jupyter_code>k_fold_sensitivity = []
k_fold_specificity = []

for i in range(5):
    # Get classificiation probabilities for k-fold replicate
    obs = observed[i]
    proba = predicted_proba[i]
    
    # Set up list for accuracy measures
    sensitivity = []
    specificity = []
    
    # Loop through increments in probability of survival
    thresholds = np.arange(0.0, 1.01, 0.01)
    for cutoff in thresholds: #  loop 0 --> 1 on steps of 0.1
        # Get classificiation using cutoff
        predicted_class = proba >= cutoff
        predicted_class = predicted_class * 1.0
        # Call accuracy measures function
        accuracy = calculate_accuracy(obs, predicted_class)
        # Add accuracy scores to lists
        sensitivity.append(accuracy['sensitivity'])
        specificity.append(accuracy['specificity'])
    
    # Add replicate to lists
    k_fold_sensitivity.append(sensitivity)
    k_fold_specificity.append(specificity)<jupyter_output><empty_output><jupyter_text>Combined plot:<jupyter_code>fig = plt.figure(figsize=(10,5))

# Plot ROC
ax1 = fig.add_subplot(121)
for i in range(5):
    ax1.plot(k_fold_fpr[i], k_fold_tpr[i], color='orange')
ax1.plot([0, 1], [0, 1], color='darkblue', linestyle='--')
ax1.set_xlabel('False Positive Rate')
ax1.set_ylabel('True Positive Rate')
ax1.set_title('Receiver Operator Characteristic Curve')
text = f'Mean AUC: {mean_auc:.3f}'
ax1.text(0.64,0.07, text, 
         bbox=dict(facecolor='white', edgecolor='black'))
plt.grid(True)

# Plot sensitivity-specificity
ax2 = fig.add_subplot(122)
for i in range(5):
    ax2.plot(k_fold_sensitivity[i], k_fold_specificity[i])
ax2.set_xlabel('Sensitivity')
ax2.set_ylabel('Specificity')
ax2.set_title('Sensitivity-Specificity Curve')
plt.grid(True)


plt.tight_layout(pad=2)

output_path = './output/'
create_folder_if_not_exist(output_path)

#plt.savefig(output_path + 'lr_single_fit_roc_sens_spec.jpg', dpi=300)
plt.savefig(output_path + 'lr_single_fit_roc_sens_spec.png', dpi=300)

plt.show()<jupyter_output><empty_output><jupyter_text>Identify cross-over of sensitivity and specificity<jupyter_code>sens = np.array(k_fold_sensitivity).mean(axis=0)
spec = np.array(k_fold_specificity).mean(axis=0)
df = pd.DataFrame()
df['sensitivity'] = sens
df['specificity'] = spec
df['spec greater sens'] = spec > sens

# find last index for senitivity being greater than specificity 
mask = df['spec greater sens'] == False
last_id_sens_greater_spec = np.max(df[mask].index)
locs = [last_id_sens_greater_spec, last_id_sens_greater_spec + 1]
points = df.iloc[locs][['sensitivity', 'specificity']]

# Get intersetction with line of x=y
a1 = list(points.iloc[0].values)
a2 = list(points.iloc[1].values)
b1 = [0, 0]
b2 = [1, 1]

intersect = get_intersect(a1, a2, b1, b2)[0]
print(f'\nIntersect: {intersect:0.3f}')<jupyter_output>Intersect: 0.820<jupyter_text>Collate and save results<jupyter_code>hospital_results = []
kfold_result = []
threshold_results = []
observed_results = []
prob_results = []
predicted_results = []

for i in range(5):
    hospital_results.extend(list(test_data[i]['StrokeTeam']))
    kfold_result.extend(list(np.repeat(i, len(test_data[i]))))
    threshold_results.extend(list(np.repeat(thresholds_calibrate[i], len(test_data[i]))))
    observed_results.extend(list(observed[i]))
    prob_results.extend(list(predicted_proba[i]))
    predicted_results.extend(list(predicted[i]))    
    
single_model = pd.DataFrame()
single_model['hospital'] = hospital_results
single_model['observed'] = np.array(observed_results) * 1.0
single_model['prob'] = prob_results
single_model['predicted'] = predicted_results
single_model['k_fold'] = kfold_result
single_model['threshold'] = threshold_results
single_model['correct'] = single_model['observed'] == single_model['predicted']

# Save
output_path = './predictions/'
create_folder_if_not_exist(output_path)
filename = output_path + 'single_fit_rf_k_fold.csv'
single_model.to_csv(filename, index=False)<jupyter_output><empty_output><jupyter_text>Feature WeightsGet Logistic Regression weights (average across k-fold results).Sort by absolute value of weight.KP edited to calculate the std dev as well as the mean weight<jupyter_code>features = X_test.columns.values

#number of values contributing to the mean calculation
n = len(feature_weights)

# Get average feature importance from k-fold
weights = np.array(feature_weights).mean(axis = 0)
weights = pd.DataFrame(data = weights, index=features)
weights.columns = ['weight_mean']
weights['weight_std_dev'] = np.array(feature_weights).std(axis = 0)
weights['weight_std_error'] = weights['weight_std_dev']/math.sqrt(n)
weights['weight_mean_abs'] = abs(weights['weight_mean'])

# Sort by importance (absolute weight)
weights.sort_values(by='weight_mean_abs', ascending=False, inplace=True)

# Save
weights.to_csv('output/lr_single_fit_feature_weights.csv')

# Display top 25
weights.head(25)<jupyter_output><empty_output><jupyter_text>Line chart:<jupyter_code># Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Get labels and values
labels = weights.index.values[0:25]
val = weights['weight_mean'].values[0:25]

# Plot
ax.plot(val, marker='o')
ax.set_ylabel('Feature weight (standardised features)')
ax.set_xticks(np.arange(len(labels)))
ax.set_xticklabels(labels)

# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=90, ha="right",
         rotation_mode="anchor")

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_line.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_line.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Bar chart:<jupyter_code># Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Get labels and values
labels = weights.index.values[0:25]
pos = np.arange(len(labels))
val = weights['weight_mean'].values[0:25]

# Plot
ax.bar(pos, val)
ax.set_ylabel('Feature weight (standardised features)')
ax.set_xticks(np.arange(len(labels)))
ax.set_xticklabels(labels)

# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=90, ha="right",
         rotation_mode="anchor")

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_bar.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Learning curveExamine the relationship between training data size and accuracy.<jupyter_code># Set up list to collect results
results_training_size = []
results_accuracy = []
results_all_accuracy = []

# Get maximum training size (number of training records)
max_training_size = train_data[0].shape[0]

# Construct training sizes (values closer at lower end)
train_sizes = [50, 100, 250, 500, 1000, 2500]
for i in range (5000, max_training_size, 5000):
    train_sizes.append(i)

# Loop through training sizes
for train_size in train_sizes:
    
    # Record accuracy across k-fold replicates
    replicate_accuracy = []

    for replicate in range(5):
        
        # Get training and test data (from first k-fold split)
        train = train_data[0]
        test = test_data[0]
        
        # One hot encode hospitals
        train_hosp = pd.get_dummies(train['StrokeTeam'], prefix = 'team')
        train = pd.concat([train, train_hosp], axis=1)
        train.drop('StrokeTeam', axis=1, inplace=True)
        test_hosp = pd.get_dummies(test['StrokeTeam'], prefix = 'team')
        test = pd.concat([test, test_hosp], axis=1)
        test.drop('StrokeTeam', axis=1, inplace=True) 
        
        # Sample from training data
        train = train.sample(n=train_size)

        # Get X and y
        X_train = train.drop('S2Thrombolysis', axis=1)
        X_test = test.drop('S2Thrombolysis', axis=1)
        y_train = train['S2Thrombolysis']
        y_test = test['S2Thrombolysis']        
   
        # Standardise X data
        X_train_std, X_test_std = standardise_data(X_train, X_test)

        # Define and Fit model
        model = LogisticRegression(solver='lbfgs')
        model.fit(X_train_std, y_train)

        # Predict test set
        y_pred_test = model.predict(X_test_std)

        # Get accuracy and record results
        accuracy = np.mean(y_pred_test == y_test)
        replicate_accuracy.append(accuracy)
        results_all_accuracy.append(accuracy)
    
    # Store mean accuracy across the k-fold splits
    results_accuracy.append(np.mean(replicate_accuracy))
    results_training_size.append(train_size)
    
k_fold_accuracy = np.array(results_all_accuracy).reshape(len(train_sizes), 5)<jupyter_output><empty_output><jupyter_text>Plot learning curve<jupyter_code>fig = plt.figure(figsize=(10,5))

ax1 = fig.add_subplot(121)

for i in range(5):
    ax1.plot(results_training_size, k_fold_accuracy[:, i])

ax1.set_xlabel('Training set size')
ax1.set_ylabel('Accuracy)')

# Focus on first 5000
ax2 = fig.add_subplot(122)
for i in range(5):
    ax2.plot(results_training_size, k_fold_accuracy[:, i])

ax2.set_xlabel('Training set size')
ax2.set_ylabel('Accuracy')
ax2.set_xlim(0, 5000)

plt.tight_layout()
plt.savefig('./output/lr_single_learning_curve.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Calibration and assessment of accuracy when model has high confidence<jupyter_code># Collate results in Dataframe
reliability_collated = pd.DataFrame()

# Loop through k fold predictions
for i in range(5):
    
    # Get observed class and predicted probability
    obs = observed[i]
    prob = predicted_proba[i]
    
    # Bin data with numpy digitize (this will assign a bin to each case)
    step = 0.10
    bins = np.arange(step, 1+step, step)
    digitized = np.digitize(prob, bins)
        
    # Put single fold data in DataFrame
    reliability = pd.DataFrame()
    reliability['bin'] = digitized
    reliability['probability'] = prob
    reliability['observed'] = obs
    classification = 1 * (prob > 0.5 )
    reliability['correct'] = obs == classification
    reliability['count'] = 1
    
    # Summarise data by bin in new dataframe
    reliability_summary = pd.DataFrame()

    # Add bins and k-fold to summary
    reliability_summary['bin'] = bins
    reliability_summary['k-fold'] = i

    # Calculate mean of predicted probability of thrombolysis in each bin
    reliability_summary['confidence'] = \
        reliability.groupby('bin').mean()['probability']

    # Calculate the proportion of patients who receive thrombolysis
    reliability_summary['fraction_positive'] = \
        reliability.groupby('bin').mean()['observed']
    
    # Calculate proportion correct in each bin
    reliability_summary['fraction_correct'] = \
        reliability.groupby('bin').mean()['correct']
    
    # Calculate fraction of results in each bin
    reliability_summary['fraction_results'] = \
        reliability.groupby('bin').sum()['count'] / reliability.shape[0]   
    
    # Add k-fold results to DatafRame collation
    reliability_collated = reliability_collated.append(reliability_summary)
    
# Get mean results
reliability_summary = reliability_collated.groupby('bin').mean()
reliability_summary.drop('k-fold', axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Reliability summary Plot results<jupyter_code>fig = plt.figure(figsize=(10,5))


# Plot predicted prob vs fraction psotive
ax1 = fig.add_subplot(1,2,1)

# Loop through k-fold reliability results
for i in range(5):
    mask = reliability_collated['k-fold'] == i
    k_fold_result = reliability_collated[mask]
    x = k_fold_result['confidence']
    y = k_fold_result['fraction_positive']
    ax1.plot(x,y, color='orange')
# Add 1:1 line
ax1.plot([0,1],[0,1], color='k', linestyle ='--')
# Refine plot
ax1.set_xlabel('Model probability')
ax1.set_ylabel('Fraction positive')
ax1.set_xlim(0, 1)
ax1.set_ylim(0, 1)

# Plot accuracy vs probability
ax2 = fig.add_subplot(1,2,2)
# Loop through k-fold reliability results
for i in range(5):
    mask = reliability_collated['k-fold'] == i
    k_fold_result = reliability_collated[mask]
    x = k_fold_result['confidence']
    y = k_fold_result['fraction_correct']
    ax2.plot(x,y, color='orange')
# Refine plot
ax2.set_xlabel('Model probability')
ax2.set_ylabel('Fraction correct')
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 1)

ax3 = ax2.twinx()  # instantiate a second axes that shares the same x-axis
for i in range(5):
    mask = reliability_collated['k-fold'] == i
    k_fold_result = reliability_collated[mask]
    x = k_fold_result['confidence']
    y = k_fold_result['fraction_results']
    ax3.plot(x,y, color='blue')
    
ax3.set_xlim(0, 1)
ax3.set_ylim(0, 0.5)
ax3.set_ylabel('Fraction of samples')

custom_lines = [Line2D([0], [0], color='orange', alpha=0.6, lw=2),
                Line2D([0], [0], color='blue', alpha = 0.6,lw=2)]

plt.legend(custom_lines, ['Fraction correct', 'Fraction of samples'],
          loc='upper center')

plt.tight_layout(pad=2)
plt.savefig('./output/lr_single_reliability.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Get accuracy of model when model is at least 80% confident<jupyter_code>bins = [0.1, 0.2, 0.9, 1.0]
acc = reliability_summary.loc[bins].mean()['fraction_correct']
frac = reliability_summary.loc[bins].sum()['fraction_results']

print ('For samples with at least 80% confidence:')
print (f'Proportion of all samples: {frac:0.3f}')
print (f'Accuracy: {acc:0.3f}')<jupyter_output>For samples with at least 80% confidence:
Proportion of all samples: 0.599
Accuracy: 0.896<jupyter_text>Observations (of the original notebook)* Overall accuracy = 83.2% (89.6% for those 60% samples with at least 80% confidence of model)* Using nominal threshold (50% probability), specificity (88%) is greater than sensitivity (72%)* The model can achieve 82.0% sensitivity and specificity simultaneously* ROC AUC = 0.904* Only marginal improvements are made above a training set size of 20k* Key features predicting use of thrombolysis are:    * Time from arrival to scan    * Stroke type    * Disability before stroke    * Presence of AFAnticoagulent_Yes    * Level of consciousness* The model shows good calibration of probability vs. fraction positive, without need of additional calibration  Onwards is the extension of the original notebook. Shapley valuesGet Shapley values for the Logistic Regression (average across k-fold results) for training set.Sort by absolute value. Calculate Shapley Value (one per feature) for each of the k-fold models, using the average across all instances.Do this in two ways:1) take mean of the raw Shapley Values1) take mean of the absolute of the Shapley ValuesStore the Shapley Values (one per feature for each k-fold) in a list. Also store the standard deviation.<jupyter_code>shap_value_mean = [] #calculates the mean of the raw shapley values
shap_value_std = [] #calculates the standard deviation of the raw shapley values
shap_value_abs_mean = [] #calculates the mean of the absolute shapley values
shap_value_abs_std = [] #calculates the standard deviation of the absolute shapley values

# calculate mean of all the instances, to get a shapley value for each feature for each kfold

for kfold in range(5):
    shap_value_mean.append(np.array(shap_values_train[kfold]).mean(axis = 0))
    shap_value_std.append(np.array(shap_values_train[kfold]).std(axis = 0))
    shap_value_abs_mean.append(np.absolute(np.array(shap_values_train[kfold])).mean(axis = 0))
    shap_value_abs_std.append(np.absolute(np.array(shap_values_train[kfold])).std(axis = 0))

# shap_value_mean has length 5. 
# shap_value_mean[0] has the shapley values for the individual features for k-fold 0 (taking a mean across all of the instances)
# shap_value_mean[1] has the shapley values for the individual features for k-fold 1 (taking a mean across all of the instances)<jupyter_output><empty_output><jupyter_text>Now calculate a Shapley Value for each feature (taking average across the k folds).Do this two ways.1) take the mean across the 5 k-folds (so a mean of the 5 means that were calculated from the raw values). Then take the absolute value "shap_values_mean_mean_abs"2) take the mean across the 5 k-folds (so a mean of the 5 means that were calculated from absolute values). "shap_values_abs_mean_mean"Calculate the standard deviation and standard error for each way (this represents the range of the means across the 5 k-folds)<jupyter_code>features = X_test.columns.values

#number of values contributing to the mean calculation
n = len(feature_weights)

# Get average feature importance across the 5 k-folds

# Calculated from raw shapley values
shap_value_mean_mean = np.array(shap_value_mean).mean(axis = 0)
shap_values_train_df = pd.DataFrame(data = shap_value_mean_mean, index=features)
shap_values_train_df.columns = ['shap_values_mean_mean']
shap_values_train_df['shap_values_mean_mean_std_dev'] = np.array(shap_value_mean).std(axis = 0)
shap_values_train_df['shap_values_mean_mean_std_error'] = shap_values_train_df['shap_values_mean_mean_std_dev']/math.sqrt(n)
shap_values_train_df['shap_values_mean_mean_abs'] = abs(shap_values_train_df['shap_values_mean_mean']) #takes the absolute of the mean

# Calculated from absolute shapley values
shap_values_train_df['shap_values_abs_mean_mean'] = np.array(shap_value_abs_mean).mean(axis = 0)# takes the mean of the mean of the absolute values
shap_values_train_df['shap_values_abs_mean_mean_std_dev'] = np.array(shap_value_abs_mean).std(axis = 0)# takes the std dev of the means of the absolute values
shap_values_train_df['shap_values_abs_mean_mean_std_error'] = shap_values_train_df['shap_values_abs_mean_mean_std_dev']/math.sqrt(n)

# Save
shap_values_train_df.to_csv('output/lr_single_fit_shapley_values_train.csv')<jupyter_output><empty_output><jupyter_text>Show top 5 features based on the Shapley Value as calculated from raw values<jupyter_code># Sort by importance, based on shapley value calculated by taking the absolute value of the mean of the mean.
shap_values_train_df.sort_values(by='shap_values_mean_mean_abs', ascending=False, inplace=True)

# Display top 25
shap_values_train_df.head(5)<jupyter_output><empty_output><jupyter_text>Show top 5 features based on the Shapley Value as calculated from absolute values<jupyter_code># Sort by importance, based on shapley value calculated by taking the mean of the mean of the absolute.
shap_values_train_df.sort_values(by='shap_values_abs_mean_mean', ascending=False, inplace=True)

# Display top 25
shap_values_train_df.head(5)<jupyter_output><empty_output><jupyter_text>Check the effect of either calculating the absolute(mean(shapley_value)) vs mean(absolute(shapley_value))<jupyter_code># Set up figure
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)

# Plot points
x = shap_values_train_df['shap_values_abs_mean_mean']
y = shap_values_train_df['shap_values_mean_mean_abs']
ax.scatter(x, y)

ax.set_title('Mean shapley value: mean(absolute) vs. absolute(mean)')
ax.set_xlabel('Shapley value (mean(absolute))')  
ax.set_ylabel('Shapley value (absolute(mean))')

# Add x=y line to the plot
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)

# and set the axis the same range and dimension
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>The feature S2BrainImagingTime_min is an outlier for both feature weights and shapley values.S2BrainImagingTime_min has a Shapley Value of -1.151812, when the rest of the features are in the range of -0.03 to 0.289.  Plot mean Shapley Value vs mean Feature WeightInclude an error bar showing the standard error (SD / root(n))<jupyter_code># join the two sets of data (shapley value, and feature weight) in a single dataframe based on the feature

shap_and_weights_df = pd.concat([shap_values_train_df['shap_values_mean_mean_abs'], 
                                 shap_values_train_df['shap_values_mean_mean_std_error'],
                                 shap_values_train_df['shap_values_abs_mean_mean'], 
                                 shap_values_train_df['shap_values_abs_mean_mean_std_error'],
                                 weights['weight_mean'], weights['weight_std_error']],
                                 axis=1)
def plot_scatter_with_error_bars(df, x_column, y_column, title):
    # plot data

    # Set up figure
    fig = plt.figure(figsize=(8,8))
    ax = fig.add_subplot(111)

    # Plot points
    x = df[f"{x_column}_mean"]#'shap_values_mean']
    y = df[y_column]#'weight_mean']
    ax.scatter(x, y)

    # Plot error bars
    x_err = df['shap_values_std_error']
    y_err = df['weight_std_error']
    ## you can use color ="r" for red or skip to default as blue
    plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

    ax.set_title(title)#'Shapley value vs. Feature weight for each feature')
    ax.set_xlabel('Shapley value')
    ax.set_ylabel('Feature weight')  

    plt.tight_layout()
    #plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
    plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
    plt.show()
    
    return()<jupyter_output><empty_output><jupyter_text>Plot the Shapley values as calculated from taking the absolute of the mean.<jupyter_code># Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df['shap_values_mean_mean_abs']
y = shap_and_weights_df['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df['shap_values_mean_mean_std_error']
y_err = shap_and_weights_df['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (taking absolute of the mean) vs. Feature weight for each feature')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Remove any features that are an outlier. For this case, remove 'S2BrainImagingTime_min'<jupyter_code>shap_and_weights_df_mod = shap_and_weights_df.drop('S2BrainImagingTime_min')
# Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df_mod['shap_values_mean_mean_abs']
y = shap_and_weights_df_mod['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df_mod['shap_values_mean_mean_std_error']
y_err = shap_and_weights_df_mod['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (taking absolute of the mean) vs. Feature weight for each feature [removed outlier]')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Remove any more features that are an outlier. In this case remove any features that have a feature weight outside -1 to 1, or above 1 for shapley value<jupyter_code>mask = ((shap_and_weights_df_mod["weight_mean"] > -1) & 
        (shap_and_weights_df_mod["weight_mean"] < 1))
shap_and_weights_df_mod = shap_and_weights_df_mod[mask]

mask = shap_and_weights_df_mod["shap_values_mean_mean_abs"] < 1
shap_and_weights_df_mod = shap_and_weights_df_mod[mask]
# Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df_mod['shap_values_mean_mean_abs']
y = shap_and_weights_df_mod['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df_mod['shap_values_mean_mean_std_error']
y_err = shap_and_weights_df_mod['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (taking absolute of the mean) vs. Feature weight for each feature [removed outliers]')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Plot the Shapley values as calculated from the absolute values.<jupyter_code># Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df['shap_values_abs_mean_mean']
y = shap_and_weights_df['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df['shap_values_abs_mean_mean_std_error']
y_err = shap_and_weights_df['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (calculating mean from the absolute values) vs. Feature weight for each feature')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Remove any features that are an outlier. For this case, remove 'S2BrainImagingTime_min'<jupyter_code>shap_and_weights_df_mod = shap_and_weights_df.drop('S2BrainImagingTime_min')
# Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df_mod['shap_values_abs_mean_mean']
y = shap_and_weights_df_mod['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df_mod['shap_values_abs_mean_mean_std_error']
y_err = shap_and_weights_df_mod['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (calculating mean from the absolute values) vs. Feature weight for each feature [removed outliers]')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Remove any more features that are an outlier. In this case remove any features that have a feature weight outside -1 to 1, or above 1 for shapley value<jupyter_code>mask = ((shap_and_weights_df["weight_mean"] > -1) & 
        (shap_and_weights_df["weight_mean"] < 1))
shap_and_weights_df_mod = shap_and_weights_df[mask]

mask = shap_and_weights_df_mod["shap_values_abs_mean_mean"] < 1
shap_and_weights_df_mod = shap_and_weights_df_mod[mask]
# Set up figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

# Plot points
x = shap_and_weights_df_mod['shap_values_abs_mean_mean']
y = shap_and_weights_df_mod['weight_mean']
ax.scatter(x, y)

# Plot error bars
x_err = shap_and_weights_df_mod['shap_values_abs_mean_mean_std_error']
y_err = shap_and_weights_df_mod['weight_std_error']
## you can use color ="r" for red or skip to default as blue
plt.errorbar(x, y, xerr=x_err, yerr=y_err, fmt="o")#, color="r")

ax.set_title('Shapley value (calculating mean from the absolute values) vs. Feature weight for each feature [removed outliers]')
ax.set_xlabel('Shapley value')
ax.set_ylabel('Feature weight')  

plt.tight_layout()
#plt.savefig('output/lr_single_fit_feature_weights_bar.jpg', dpi=300)
plt.savefig('output/lr_single_fit_feature_weights_vs_shapley_values.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>1. Show the top 10 (absolute value) Shapley values and Feature weights<jupyter_code>shap_and_weights_df = shap_and_weights_df.sort_values(by=['shap_values_mean_mean_abs'], ascending=False)
shap_and_weights_10_df = pd.DataFrame(data = shap_and_weights_df.index[0:10], 
                                      columns = ['absolute of mean shapley values ranked'])

shap_and_weights_df = shap_and_weights_df.sort_values(by=['shap_values_abs_mean_mean'], ascending=False)
shap_and_weights_10_df['mean of absolute shapley values ranked'] = shap_and_weights_df.index[0:10]

shap_and_weights_df['weight_mean_abs'] = shap_and_weights_df['weight_mean'].abs() 
shap_and_weights_df = shap_and_weights_df.sort_values(by=['weight_mean_abs'], ascending=False)
shap_and_weights_10_df['absolute feature weights ranked'] = shap_and_weights_df.index[0:10]

shap_and_weights_10_df<jupyter_output><empty_output><jupyter_text>Part 2: Fit a logistic regresion model on just the first k-fold split (create SHAP plots)Only use the first k fold train/test split to fit a logistic regression model, and recalculate the Shapley values and model weightsCompare average Shap values in train and test (should be similar!) Pick three examples from the test set to show how Shap explains examples not used in training.<jupyter_code>chosen_k_fold = 3

# Get k fold split
train = train_data[chosen_k_fold]
test = test_data[chosen_k_fold]

# Get X and y
X_train = train.drop('S2Thrombolysis', axis=1)
X_test = test.drop('S2Thrombolysis', axis=1)
y_train = train['S2Thrombolysis']
y_test = test['S2Thrombolysis']

# One hot encode hospitals
X_train_hosp = pd.get_dummies(X_train['StrokeTeam'], prefix = 'team')
X_train = pd.concat([X_train, X_train_hosp], axis=1)
X_train.drop('StrokeTeam', axis=1, inplace=True)
X_test_hosp = pd.get_dummies(X_test['StrokeTeam'], prefix = 'team')
X_test = pd.concat([X_test, X_test_hosp], axis=1)
X_test.drop('StrokeTeam', axis=1, inplace=True)

# Standardise X data
X_train_std, X_test_std = standardise_data(X_train, X_test)

# The above function (standardise_data) returns a numpy array. 
# Using X_train_std in the SHAP plots means no feature names are used on the plot
# To get feature names to appear on the plots, save the numpy array as a pandas Dataframe (with feature names as columns)
X_train_std_df = pd.DataFrame(X_train_std, columns=X_train.columns)
X_test_std_df = pd.DataFrame(X_test_std, columns=X_test.columns)

# Define and Fit model
model = LogisticRegression(solver='lbfgs')
model.fit(X_train_std, y_train)

# Get feature weights
weights = model.coef_[0]

# Get predicted probabilities
y_probs = model.predict_proba(X_test_std)[:,1]

# Calibrate model and get class
true_rate = np.mean(y_test)
threshold = find_threshold(y_probs, true_rate)
y_class = y_probs >= threshold
y_class = np.array(y_class) * 1.0

# Calculate Shapley values
### Global influence of features in the training set
### Here we are looking at Shapley values for the training set, to understand the influence of features in training.
(shapley_values_train, shapley_values_train_extended) = calculate_shapley_values(model, X_train_std_df)

### And for the Test set
(shapley_values_test, shapley_values_test_extended) = calculate_shapley_values(model, X_test_std_df)

# Print accuracy
accuracy = np.mean(y_class == y_test)
print (
    f'Accuracy: {accuracy:0.3f}, threshold {threshold:0.3f}')<jupyter_output>Accuracy: 0.833, threshold 0.470<jupyter_text>Compare average Shap values in train and test (should be similar - but guess what!)  Create dataframe for shap values for the train and test set<jupyter_code># Calculate shap value for train and test set (take average across all instances)
shap_value_train_mean = np.array(shapley_values_train).mean(axis = 0)
shap_value_train_median = np.median(np.array(shapley_values_train), axis = 0)
shap_value_train_std = (np.array(shapley_values_train).std(axis = 0))
shap_value_test_mean = (np.array(shapley_values_test).mean(axis = 0))
shap_value_test_median = np.median(np.array(shapley_values_test), axis = 0)
shap_value_test_std = (np.array(shapley_values_test).std(axis = 0))

# store in dataframe
shap_values_train_test_df = pd.DataFrame(data = shap_value_train_mean,
                                        columns = ["shap_value_train_mean"])
shap_values_train_test_df["shap_value_train_median"] = shap_value_train_median
shap_values_train_test_df["shap_value_train_std"] = shap_value_train_std
shap_values_train_test_df["shap_value_test_mean"] = shap_value_test_mean
shap_values_train_test_df["shap_value_test_median"] = shap_value_test_median
shap_values_train_test_df["shap_value_test_std"] = shap_value_test_std
shap_values_train_test_df.set_index(X_train.columns.values, inplace=True)<jupyter_output><empty_output><jupyter_text>To look at the values for a certain feature, use this code:<jupyter_code>shap_values_train_test_df.loc["S2BrainImagingTime_min"]<jupyter_output><empty_output><jupyter_text>Plot scatter plot of Shapley Values from train vs test (repeat for both mean and median shapley value)<jupyter_code># Set up figure
fig = plt.figure(figsize=(12,12))

# Mean plot
ax = fig.add_subplot(121)

# Plot points
x = shap_values_train_test_df["shap_value_train_mean"]
y = shap_values_train_test_df["shap_value_test_mean"]
ax.scatter(x, y)

ax.set_title('Mean shapley value: train vs. test')
ax.set_xlabel('Shapley value Train (mean)')
ax.set_ylabel('Shapley value Test (mean)')  

# Add x=y line to the plot
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)

# and set the axis the same range and dimension
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)

# Median plot
ax = fig.add_subplot(122)

# Plot points
x = shap_values_train_test_df["shap_value_train_median"]
y = shap_values_train_test_df["shap_value_test_median"]
ax.scatter(x, y)

ax.set_title('Median shapley value: train vs. test')
ax.set_xlabel('Shapley value Train (median)')
ax.set_ylabel('Shapley value Test (median)')  

# Add x=y line to the plot
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)

# and set the axis the same range and dimension
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)


plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Expected a nice linear x=y plot. There appears to be three outliers.Let's investigate.Rank the features by their absolute difference in Shapley Value between train and test<jupyter_code>difference = (shap_values_train_test_df["shap_value_train_median"]-shap_values_train_test_df["shap_value_test_median"])
shap_values_train_test_df["difference"] = difference.abs()
shap_values_train_test_df.sort_values(by="difference", ascending=False, inplace=True)<jupyter_output><empty_output><jupyter_text>Show top 10 features with greatest difference in Shapley Value between train and test (to identify the outliers)<jupyter_code>shap_values_train_test_df.head(10)<jupyter_output><empty_output><jupyter_text>Remove the 3 features with the greatest difference between train and test (the three outliers)<jupyter_code>shap_values_train_test_df_mod = shap_values_train_test_df.drop(shap_values_train_test_df.index[0:3])
# Set up figure
fig = plt.figure(figsize=(15,15))

ax = fig.add_subplot(121)

# Plot points
x = shap_values_train_test_df_mod["shap_value_train_mean"]
y = shap_values_train_test_df_mod["shap_value_test_mean"]
ax.scatter(x, y)

# Add x=y line to the plot
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)

# and set the axis the same range and dimension
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)

ax.set_title('Mean shapley value: train vs. test')
ax.set_xlabel('Shapley value Train (mean)')
ax.set_ylabel('Shapley value Test (mean)')  


ax = fig.add_subplot(122)

# Plot points
x = shap_values_train_test_df_mod["shap_value_train_median"]
y = shap_values_train_test_df_mod["shap_value_test_median"]
ax.scatter(x, y)

# Add x=y line to the plot
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)

# and set the axis the same range and dimension
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)

ax.set_title('Median shapley value: train vs. test')
ax.set_xlabel('Shapley value Train (median)')
ax.set_ylabel('Shapley value Test (median)')  

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>What's the top 10 features for Shapley value and feature weight for this model? For both the train and test set?(I look at this because for the k-fold data, feature "S2BrainImagingTime_min" was an outlier for the global importance, but when looked at individual instances, this feature was often third in importance.<jupyter_code>features = X_test.columns.values

# Create a dataframe with shapley values ordered by absolute value (for training set)
shapley_values_train_ordered = pd.DataFrame(data = shapley_values_train.mean(axis = 0), index=features)
shapley_values_train_ordered.columns = ['shap_values_mean']
shapley_values_train_ordered['shap_values_mean_abs'] = abs(shapley_values_train_ordered['shap_values_mean'])
# Sort by importance (absolute)
shapley_values_train_ordered.sort_values(by='shap_values_mean_abs', ascending=False, inplace=True)

# Create a dataframe with shapley values ordered by absolute value (for test set)
shapley_values_test_ordered = pd.DataFrame(data = shapley_values_test.mean(axis = 0), index=features)
shapley_values_test_ordered.columns = ['shap_values_mean']
shapley_values_test_ordered['shap_values_mean_abs'] = abs(shapley_values_test_ordered['shap_values_mean'])
# Sort by importance (absolute)
shapley_values_test_ordered.sort_values(by='shap_values_mean_abs', ascending=False, inplace=True)

# Create a dataframe with feature weights ordered by absolute value
feature_weights_ordered = pd.DataFrame(data = weights, index=features)
feature_weights_ordered.columns = ['weight_mean']
feature_weights_ordered['weight_mean_abs'] = abs(feature_weights_ordered['weight_mean'])
# Sort by importance (absolute weight)
feature_weights_ordered.sort_values(by='weight_mean_abs', ascending=False, inplace=True)

# Join the results fo teh top 10 shapley value features and top 10 feature weights in a dataframe
shap_and_weights_10_df = pd.DataFrame(data = shapley_values_train_ordered.index[0:10], 
                                      columns = ['absolute shapley values train ranked'])
shap_and_weights_10_df['absolute shapley values test ranked'] = shapley_values_test_ordered.index[0:10]
shap_and_weights_10_df['absolute feature weights ranked'] = feature_weights_ordered.index[0:10]

shap_and_weights_10_df.head(10)<jupyter_output><empty_output><jupyter_text>Show a Beeswarm plot for all the dataUsing code from https://samuel-book.github.io/samuel-1/random_forest/characterise_contentious_patients.htmlobservations  Shap bee swarm plotKP: still working on this subplot, as they are overlapping. See below for separate plots.<jupyter_code># Set up figure
fig = plt.figure(figsize=(60,5))

print ("LHS: Training dataset. RHS Test dataset")

# Left plot of the training set
ax = fig.add_subplot(121)

ax = shap.summary_plot(shap_values = shapley_values_train, 
                       features = X_train_std_df.values,
                       feature_names = X_train_std_df.columns.values, 
                       cmap = plt.get_cmap('nipy_spectral'), 
                       show=False)

# Right plot of the test set
ax = fig.add_subplot(122)

ax = shap.summary_plot(shap_values = shapley_values_test, 
                       features = X_test_std_df.values,
                       feature_names = X_test_std_df.columns.values, 
                       cmap = plt.get_cmap('nipy_spectral'), 
                       show=False)
plt.tight_layout()

# Trying to get the subplots to not overlay
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)
#left  = 0.125  # the left side of the subplots of the figure
#right = 0.9    # the right side of the subplots of the figure
#bottom = 0.1   # the bottom of the subplots of the figure
#top = 0.9      # the top of the subplots of the figure
#wspace = 0.2   # the amount of width reserved for blank space between subplots
#hspace = 0.2   # the amount of height reserved for white space between subplots

plt.savefig('output/shap_beeswarm.png', dpi=300)
fig = plt.figure(figsize=(6,6))

print ("Training dataset")

shap.summary_plot(shap_values = shapley_values_train, 
                  features = X_train_std_df.values,
                  feature_names = X_train_std_df.columns.values, 
                  cmap = plt.get_cmap('nipy_spectral'), 
                  show=False)
plt.tight_layout()

plt.savefig('output/shap_beeswarm.png', dpi=300)
fig = plt.figure(figsize=(6,6))

print ("Test dataset")

shap.summary_plot(shap_values = shapley_values_test, 
                  features = X_test_std_df.values,
                  feature_names = X_test_std_df.columns.values, 
                  cmap = plt.get_cmap('nipy_spectral'), 
                  show=False)
plt.tight_layout()

plt.savefig('output/shap_beeswarm.png', dpi=300)<jupyter_output>Test dataset<jupyter_text>Remove the outlier S2BrainImagingTime_min, and replot the beeswarm<jupyter_code>shapley_values_train_mod = np.delete(shapley_values_train, X_test.columns.get_loc("S2BrainImagingTime_min"), 1)
X_train_std_df_mod = X_train_std_df.drop("S2BrainImagingTime_min", axis=1)

shapley_values_test_mod = np.delete(shapley_values_test, X_test.columns.get_loc("S2BrainImagingTime_min"), 1)
X_test_std_df_mod = X_test_std_df.drop("S2BrainImagingTime_min", axis=1)
fig = plt.figure(figsize=(6,6))

print ("Training dataset")

shap.summary_plot(shap_values = shapley_values_train_mod, 
                  features = X_train_std_df_mod.values,
                  feature_names = X_train_std_df_mod.columns.values, 
                  cmap = plt.get_cmap('nipy_spectral'), show=False)
plt.tight_layout()

plt.savefig('output/shap_beeswarm.png', dpi=300)
fig = plt.figure(figsize=(6,6))

print ("Training dataset")

shap.summary_plot(shap_values = shapley_values_test_mod, 
                  features = X_test_std_df_mod.values,
                  feature_names = X_test_std_df_mod.columns.values, 
                  cmap = plt.get_cmap('nipy_spectral'), show=False)
plt.tight_layout()

plt.savefig('output/shap_beeswarm.png', dpi=300)<jupyter_output>Training dataset<jupyter_text>Plot Waterfall plots for three instances: Taken from the test set, an example with a low probability of giving thrombolysis, one with a probability close to 0.5, and one with a high probability of giving thrombolysis. Use code for waterfall from 05_keras_tabular-waterfall.ipynb(If possible, plot in a single wide MatPlotLib plot)  Identify three instances from the test set, to cover the range for y_probs: 0, 0.5, 1.Take the first occurance that satisfies the criteria (by using "[0][0]")<jupyter_code>location_low_probability = np.where(y_probs < 0.01)[0][0]
location_middle_probability = np.where((y_probs > 0.49) & (y_probs < 0.51))[0][0]
location_high_probability = np.where(y_probs > 0.99)[0][0]<jupyter_output><empty_output><jupyter_text>Waterfall plot for low probability of thrombolysis (showing log-odds)<jupyter_code>shap.plots.waterfall(shapley_values_test_extended[location_low_probability], 
                     max_display=15)<jupyter_output><empty_output><jupyter_text>Waterfall plots (showing log-odds), one for each instance (low, middle, high change for thrombolysis)<jupyter_code>fig = plt.figure(figsize=(10,5))

ax1 = fig.add_subplot(131)
shap.plots.waterfall(shapley_values_test_extended[location_low_probability], max_display=15)#, show=False)
ax1.title.set_text('Low thrombolysis probability')

ax2 = fig.add_subplot(132)
shap.plots.waterfall(shapley_values_test_extended[location_middle_probability], max_display=15)#, show=False)
ax2.title.set_text('Medium thrombolysis probability')

ax3 = fig.add_subplot(133)
shap.plots.waterfall(shapley_values_test_extended[location_high_probability], max_display=15)#, show=False)
ax3.title.set_text('High thrombolysis probability')

fig.set_figheight(10)
fig.set_figwidth(50)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>(using content from "https://medium.com/dataman-in-ai/the-shap-with-more-elegant-charts-bc3e73fa1c0c")By default, SHAP explains models in terms of their margin output before the logistic link function [https://shap.readthedocs.io/en/latest/example_notebooks/api_examples/plots/waterfall.html], hence the x axis units in the waterfall plot are log-odds units, and not probability.You can convert the log-odd to a probability of [0,1] by using the logistic sigmoid function, which is expit(x) = 1/(1+exp(-x)), the inverse of the logit function. In other words, 1(1+exp(-4.894) ) = 0.992. Here's a utility function that does the conversion from the log-odds to probabilities (with thanks for Dataman in AI).<jupyter_code>def shap_transform_scale(original_shap_values, Y_pred, which):
    
    #Compute the transformed base value, which consists in applying the logit function to the base value
    from scipy.special import expit #Importing the logit function for the base value transformation
    untransformed_base_value = original_shap_values.base_values[-1]
   
    #Computing the original_explanation_distance to construct the distance_coefficient later on
    original_explanation_distance = np.sum(original_shap_values.values, axis=1)[which]
    
    base_value = expit(untransformed_base_value ) # = 1 / (1+ np.exp(-untransformed_base_value))

    #Computing the distance between the model_prediction and the transformed base_value
    distance_to_explain = Y_pred[which] - base_value

    #The distance_coefficient is the ratio between both distances which will be used later on
    distance_coefficient = original_explanation_distance / distance_to_explain

    #Transforming the original shapley values to the new scale
    shap_values_transformed = original_shap_values / distance_coefficient

    #Finally resetting the base_value as it does not need to be transformed
    shap_values_transformed.base_values = base_value
    shap_values_transformed.data = original_shap_values.data
    
    #Now returning the transformed array
    return shap_values_transformed<jupyter_output><empty_output><jupyter_text>Let’s apply the conversion, then do the waterfall plots  Waterfall plot for low probability of thrombolysis (showing probability)<jupyter_code>print("The prediction is ", y_probs[location_low_probability])
shap_values_transformed = shap_transform_scale(shapley_values_test_extended, y_probs, location_low_probability)
shap.plots.waterfall(shap_values_transformed[location_low_probability])<jupyter_output><empty_output><jupyter_text>Waterfall plots (showing probability), one for each instance (low, middle, high chance for thrombolysis)<jupyter_code>fig = plt.figure(figsize=(10,5))

expected_value = shapley_values_test_extended.base_values[0]

ax1 = fig.add_subplot(131)
shap_values_transformed = shap_transform_scale(shapley_values_test_extended, y_probs, location_low_probability)
shap.plots.waterfall(shap_values_transformed[location_low_probability], max_display=15)#, show=False)
ax1.title.set_text('Low thrombolysis probability')

ax2 = fig.add_subplot(132)
shap_values_transformed = shap_transform_scale(shapley_values_test_extended, y_probs, location_middle_probability)
shap.plots.waterfall(shap_values_transformed[location_middle_probability], max_display=15)#, show=False)
ax2.title.set_text('Medium thrombolysis probability')

ax3 = fig.add_subplot(133)
shap_values_transformed = shap_transform_scale(shapley_values_test_extended, y_probs, location_high_probability)
shap.plots.waterfall(shap_values_transformed[location_high_probability], max_display=15)#, show=False)
ax3.title.set_text('High thrombolysis probability')

#fig.set_figheight(15)
fig.set_figwidth(25)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Decision plots showing log-odds  For low probability thrombolysis<jupyter_code>instance = location_low_probability

chosen_instance = X_test.iloc[instance] #used for the feature value (in brackets on graph)
shap_values_instance = shapley_values_test[instance] #used for charted value
expected_value = shapley_values_test_extended.base_values[0]#explainer.expected_value #The starting point
print(f"The expected value is {expected_value}")
print(f"Chart is dislaying in brackets the raw (not standardised) feature values (model was trained on standardised values)")

shap.decision_plot(expected_value, shap_values_instance, chosen_instance)
instance = location_low_probability

chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
expected_value = shapley_values_test_extended.base_values[0]#explainer.expected_value
print(f"The expected value is {expected_value}")
print(f"Chart is dislaying in brackets the standardised feature values (model was trained on standardised values)")

shap.decision_plot(expected_value, shap_values_instance, chosen_instance)<jupyter_output><empty_output><jupyter_text>Decision plots showing three instances on a subplot (log-odds)<jupyter_code>fig = plt.figure(figsize=(10,5))

expected_value = shapley_values_test_extended.base_values[0]#explainer.expected_value

ax1 = fig.add_subplot(131)
instance = location_low_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, show=False)
ax1.title.set_text('Low thrombolysis probability (displaying log-odds)')

ax2 = fig.add_subplot(132)
instance = location_middle_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, show=False)
ax2.title.set_text('Medium thrombolysis probability (displaying log-odds)')

ax3 = fig.add_subplot(133)
instance = location_high_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, show=False)
ax3.title.set_text('High thrombolysis probability (displaying log-odds)')

#fig.set_figheight(15)
fig.set_figwidth(25)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Decision plots showing three instances on a subplot (probability)<jupyter_code>fig = plt.figure(figsize=(10,5))

expected_value = shapley_values_test_extended.base_values[0]#explainer.expected_value

ax1 = fig.add_subplot(131)
instance = location_low_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, link='logit',show=False)
ax1.title.set_text('Low thrombolysis probability (displaying probability)')

ax2 = fig.add_subplot(132)
instance = location_middle_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, link='logit', show=False)
ax2.title.set_text('Medium thrombolysis probability (displaying probability)')

ax3 = fig.add_subplot(133)
instance = location_high_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shapley_values_test[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, link='logit', show=False)
ax3.title.set_text('High thrombolysis probability (displaying probability)')

#fig.set_figheight(15)
fig.set_figwidth(25)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Included this for completeness, using the same method for the Waterfall plots to transform the log-odds to probability<jupyter_code>instance_location = location_low_probability

print("The prediction is ", y_probs[instance_location])

shap_values_transformed = shap_transform_scale(shapley_values_test_extended, y_probs, instance_location)
shap_values_transformed_instance = shap_values_transformed[instance_location]

chosen_instance = X_test_std_df.iloc[instance_location]

expected_value = shap_values_transformed.base_values
print(f"The expected value is {expected_value}")

shap.decision_plot(expected_value, shap_values_transformed_instance.values, chosen_instance)<jupyter_output><empty_output><jupyter_text>(1) Global Interpretability (1.1) Bar plot for feature importance<jupyter_code>shap.plots.bar(shapley_values_test_extended, max_display=10) # default is max_display=12<jupyter_output><empty_output><jupyter_text>(1.2) Cohort plotGet an error when I run this code. I believe it is an error within the ".cohorts()" procedure as I can replicate the error just running "shapley_values_extended_df.cohorts(2)".Change the cell below from Markdown to Code to work on it. shap.plots.bar(shapley_values_extended_df.cohorts(2).abs.mean(0))  (1.3) Heatmap plot<jupyter_code>shap.plots.heatmap(shapley_values_test_extended[0:500])
shap.plots.heatmap(shapley_values_test_extended[500:1000])<jupyter_output><empty_output><jupyter_text>(2) Local Interpretability  (2.1) Waterfall plots (see earlier section)  (2.2) The bar plot for an individual case<jupyter_code>shap.plots.bar(shapley_values_test_extended[location_low_probability]) # For the low probability of thrombolysis instance<jupyter_output><empty_output><jupyter_text>(2.3) Decision plots (see earlier section)  This next section is now redundant. It was coded Waterfall plot, as I couldn't get the waterfall plot to include the Feature Names on the y axis. Thought that this was a issue for logistic regression shapley values. In fact it was because the standardise fuction returned a numpy array for X_train_std. The function needed a dataframe in order to use the column values for the chart. So I saved the X_train_std as a dataframe and the feature names now exist on the charts.****************************************** 2. Plotting waterfall using own code Code to create waterfall plot for shapley values as calculated by DeepExplainer (from the SHAP library)We will use the data as calculated in the previous code cell. This is data for a single instance:1) chosen_instance = X.iloc[instance]2) shap_values_instance = shapley_values[instance]3) explainer.expected_valueIn addition to the variables set up above (choosen_instance, shap_values_instance, expected_value), we also want the predicted output value from the model. Use this to compare with the culmulative sum of the shapley value. Import libraries<jupyter_code>import matplotlib.patches as mpatches<jupyter_output><empty_output><jupyter_text>Define function to create waterfall plot<jupyter_code>def waterfall_plot(chosen_instance, shap_values_instance, expected_value, 
                   predicted_value, max_display=-1, style_arrows=False, set_show=True):
    """
    NEED TO IMPROVE THIS FUNCTION SO THAT LABEL PLACEMENT ON THE PLOT, AND PLOT SIZE ADJUSTMENT WORKS FOR ALL PLOT SIZES
    AT THE MOMENT, I AM CUSTOMISING THE ADJUSTMENTS FOR EACH GRAPH.
    THE CODE THAT NEEDS ADJUSTING IS MARKED WITH #***
    
    chosen_instance: [pandas dataframe] the index is feature name, value is 
                     feature value for instance
    shap_values_instance: [numpy array] the shapley values for each feature for 
                          the instance
    expected_value: [float] the output value that is expected for the dataset 
                    (e.g. median)
    predicted_value: [float] the predicted output for the instance
    max_display: [integer] the number of features to include in the waterfall
                  plot. If nothing is passed, all features are plotted
    style_arrows: [boolean] whether to display as a horizontal bar chart 
                  (False), or using arrows (True)
    set_show: [boolean] set to false when creating a subplot (default is True)
    
    """
    
    # set up pandas dataframe with the data required for the waterfall chart
    data = pd.DataFrame(chosen_instance)
    data.columns = ["instance"]
    data["shap"] = shap_values_instance
    data["shap_absolute"] = np.abs(shap_values_instance)
    
    # limit the number of features to display to the number available
    max_display = min(max_display, data.shape[0])

    # if max_display not passed, set to total number of features available
    if max_display == -1:
        max_display = data.shape[0]
        
    # sort the rows on the absolute shapley values (smallest value at top)
    data = data.sort_values(by=["shap_absolute"], 
                                      ascending=True)

    # if want to display fewer features than the full set
    if max_display < data.shape[0]:

        # select features ranked by most influence
        waterfall_plot = data.tail(max_display)    

        # represent the remaining features as one row
        waterfall_plot_remaining = data.iloc[0:-max_display]
        summary_remaining = pd.DataFrame(
                                    columns=waterfall_plot_remaining.columns)
        summary_remaining["instance"] = [np.nan]
        summary_remaining["shap"] = [waterfall_plot_remaining.shap.sum()]
        summary_remaining["shap_absolute"] = (
                                [np.abs(waterfall_plot_remaining.shap.sum())])
        index_title = f"remaining_{waterfall_plot_remaining.shape[0]}_features"
        summary_remaining.index = [index_title]
        waterfall_plot = summary_remaining.append(waterfall_plot)

    else:
        waterfall_plot = data.copy(deep=True)
        
    # create two new columns containing the x coordinates for the lhs and rhs of
    # each bar    
    # one side of the first bar is the expected value plus the bar's shapley 
    # value. The remaining bars take as one side the culmulative sum of the
    # shapley values and the expected value.
    shap_series = waterfall_plot.shap.copy(deep=True)
    shap_series[0] = shap_series[0] + expected_value#explainer.expected_value
    shap_cumsum = shap_series.cumsum()

    # the other side of each bar is the value from the previous bar (for the 
    # first bar, this is the expected value)
    shap_cumsum_prev = shap_cumsum.copy(deep=True)
    shap_cumsum_prev[1:] = shap_cumsum_prev[0:-1]
    shap_cumsum_prev[0] = expected_value#explainer.expected_value

    # the shapley values can be either sign, so the LHS is the min of both
    # series, and the RHS is the max of both series
    # add these to the waterfall_plot DataFrame
    lhs_bar = pd.DataFrame([shap_cumsum_prev, shap_cumsum]).min()
    rhs_bar = pd.DataFrame([shap_cumsum_prev, shap_cumsum]).max()
    waterfall_plot["lhs_bar"] = lhs_bar
    waterfall_plot["rhs_bar"] = rhs_bar
    
    model_predict = shap_cumsum[-1]

    # add the expected and predicted value for this instance to the DataFrame

    # 220117 removed for now as they can be in different units: log-odds and probability
#    abs_diff = np.abs(model_predict-predicted_value)
#    if abs_diff > 0.000001:
#        print(f"There is a difference between the predicted value from the model"
#              f" and the culmulative sum of the shapley values: {abs_diff}")

    model_predict = np.repeat(model_predict, waterfall_plot.shape[0])
    waterfall_plot["model_predict"] = model_predict    
#    exp_value = np.repeat(explainer.expected_value, waterfall_plot.shape[0])
    exp_value = np.repeat(expected_value, waterfall_plot.shape[0])
    waterfall_plot["expected_value"] = exp_value
        
    # the plot will show blue for positive shapley value, and red for negative
    c = waterfall_plot.shap.copy(deep = True)
    mask_lt = c < 0
    mask_gt = c >= 0
    c[mask_lt] = "r"
    c[mask_gt] = "b"
    waterfall_plot["bar_colour"] = c
    
    # set up other parameters for the waterfall plot:
    # number of features
    n_features = waterfall_plot.shape[0]
    # positions for the y axis items
    pos = np.arange(0.5, n_features * 0.5 + 0.5, 0.5)
    # labels for the y axis (feature name and feature value)
    labels = waterfall_plot.index.astype(str) + " (" + waterfall_plot.instance.astype(str) + ")"
    
    # CONSTANTS TO ADJUST AXIS RANGE AND LABEL POSITIONS
    # find the range of the x axis, so can use this to determine the constants 
    # to use in order to place the labels on the graph
    x_axis_min = lhs_bar.min()
    x_axis_max = rhs_bar.max()
    x_axis_range = abs(x_axis_min) + abs(x_axis_max)
    
    # constant to use for the lhs label
#    lhs_label_adjustment = 3/x_axis_range
    lhs_label_adjustment = x_axis_range/45
    # constant to use for the rhs label
#    rhs_label_adjustment = 0.5/x_axis_range
    rhs_label_adjustment = x_axis_range/80
    # constant to use for the lhs x_lim for axis
    #x_left*1.05
    x_axis_lhs_adjustment = 0.95
    if x_axis_min < 0:
        x_axis_lhs_adjustment = 1.05
    
    # constant to use for the rhs x_lim for axis
    #x_right*0.95
    x_axis_rhs_adjustment = 1.05
    if x_axis_max < 0:
        x_axis_rhs_adjustment = 0.95

    ### ### ### create the plot ### ### ###
    
    # set size of plot (width is constant, height is based on number of features)
    fig = plt.figure(figsize=(20, 0.4 * (max_display + 1)))
    ax1 = fig.add_subplot(111)

    # for each feature of the model create a horizontal bar
    for i in range(n_features):
        
        # get the position for the bar
        lhs_bar = waterfall_plot.lhs_bar.iloc[i]
        rhs_bar = waterfall_plot.rhs_bar.iloc[i]

        # allocate colour (based on sign of shapley value)
        colour = waterfall_plot.bar_colour.iloc[i]

        if style_arrows:
            # create plot with arrows
            
            # direction of the arrow depends on the sign of the shapley value
            if waterfall_plot.shap.iloc[i] < 0:
                x_tail = rhs_bar
                x_head = lhs_bar
            else:
                x_tail = lhs_bar
                x_head = rhs_bar

            # height of arrow the same regardless of sign of shapley value
            y_tail = (i * 0.5) + 0.5
            y_head = (i * 0.5) + 0.5
            
            # create arrow
            arrow = mpatches.FancyArrowPatch((x_tail, y_tail), 
                                             (x_head, y_head),
                                              mutation_scale=10,
                                              color = colour,
                                              ec = colour)
            # plot arrow
            ax1.add_patch(arrow)

        else: 
            # create plot with bars
            ax1.barh((i * 0.5) + 0.5, 
                     rhs_bar - lhs_bar, 
                     left = lhs_bar, 
                     height = 0.3, 
                     align = 'center', 
                     edgecolor = colour, 
                     color = colour, 
                     alpha = 0.8)


        # add text (the shapley value) to each bar
        # location depends on sign of shapley value
        
            #*** THE ADJUSTMENTS HERE (CONSTANTS FOR text_xaxis and test_yaxis)
            #*** NEED TO BE DONE A DIFFERENT WAY. AT THE MOMENT NEED TO CHANGE
            #*** AND CUSTOMISE FOR EACH GRAPH DEPENDING OF GRAPH SIZE
        if waterfall_plot.shap.iloc[i] < 0:
            text_xaxis = lhs_bar - lhs_label_adjustment
        else:
            text_xaxis = rhs_bar + rhs_label_adjustment
        text_yaxis = (i * 0.5) + 0.4

        ax1.text(text_xaxis, 
                 text_yaxis,
                 str(round(waterfall_plot.shap.iloc[i], 2)), 
                 fontsize=8,
                 color=colour)

    # expand xlim and ylim to make labels easier to read

        #*** THE ADJUSTMENTS HERE (CONSTANTS FOR set_ylim and set_xlim)
        #*** NEED TO BE DONE A DIFFERENT WAY. AT THE MOMENT NEED TO CHANGE
        #*** AND CUSTOMISE FOR EACH GRAPH DEPENDING OF GRAPH SIZE, AND ALSO 
        #*** FOR WHETHER THE X AXIS IS + or - AT IT'S EXTREMES
        
    ax1.set_ylim(0,np.max(pos)*1.02)
    x_left, x_right = ax1.get_xlim()
#    ax1.set_xlim(x_left*1.05, x_right*0.95)
    ax1.set_xlim(x_left*x_axis_lhs_adjustment, 
                 x_right*x_axis_rhs_adjustment)
 
    x_left, x_right = ax1.get_xlim()
    
    #convert values into axis units so the arrow points to the right place
    expected_value_xaxis = ((waterfall_plot.expected_value[0] -x_left) / 
                            (x_right - x_left))
    model_predict_xaxis = ((waterfall_plot.model_predict[0] -x_left) / 
                           (x_right - x_left))

    # set properties for the arrow
    arrow_args = dict(arrowstyle="->")

    # annotate the expected value on the bottom of the x axis
    ax1.annotate(f'E[f(x)] = {waterfall_plot.expected_value[0]}', 
                 xy=(expected_value_xaxis, 0), 
                 xycoords='axes fraction',
                 xytext=(-70, -40), 
                 textcoords='offset points',
                 ha="left", va="bottom",
                 arrowprops=arrow_args)

    # annotate the predicted value on the top of the x axis
    ax1.annotate(f'f(x)= {waterfall_plot.model_predict[0]}', 
                 xy=(model_predict_xaxis, 1), 
                 xycoords='axes fraction',
                 xytext=(-80, 20), 
                 textcoords='offset points',
                 ha="left", va="bottom",
                 arrowprops=arrow_args)

    ax1.set_yticks(pos)
    ax1.set_yticklabels(labels)
    ax1.set_title('Shapley values for each feature for a single instance')

    plt.show()<jupyter_output><empty_output><jupyter_text>Example waterfall plot for a patient with low probability thrombolysis<jupyter_code>instance = location_low_probability

chosen_instance = X_test.iloc[instance]
shap_values_instance = shapley_values_test[instance]
expected_value = shapley_values_test_extended.base_values[0]#np.repeat(explainer.expected_value, len(shapley_values))
chosen_predict = y_probs[instance]

waterfall_plot(chosen_instance, 
               shap_values_instance, 
               expected_value,
               chosen_predict)<jupyter_output><empty_output><jupyter_text>Remove the features that do not contribute to the predicted value (have a shapley value = 0)<jupyter_code>mask = shap_values_instance == 0
mask = np.invert(mask)
chosen_instance_filtered = chosen_instance[mask]
shap_values_instance_filtered = shap_values_instance[mask]
expected_value_filtered =  shapley_values_test_extended.base_values[0]#expected_value[0]
expected_value_filtered

waterfall_plot(chosen_instance_filtered, 
               shap_values_instance_filtered, 
               expected_value_filtered,
               chosen_predict)<jupyter_output><empty_output><jupyter_text>Or can specify the top n features to display<jupyter_code>waterfall_plot(chosen_instance, 
               shap_values_instance, 
               expected_value,
               chosen_predict,
               max_display=15)
### Example waterfall plot for a patient with middle probability thrombolysis
instance = location_middle_probability

chosen_instance = X_test.iloc[instance]
shap_values_instance = shapley_values_test[instance]
expected_value = shapley_values_test_extended.base_values[0]#np.repeat(explainer.expected_value, len(shapley_values))
chosen_predict = y_probs[instance]

waterfall_plot(chosen_instance, 
               shap_values_instance, 
               expected_value,
               chosen_predict, 
               max_display=15)
### Example waterfall plot for a patient with high probability thrombolysis
instance = location_high_probability

chosen_instance = X_test.iloc[instance]
shap_values_instance = shapley_values_test[instance]
expected_value = shapley_values_test_extended.base_values[0]#np.repeat(explainer.expected_value, len(shapley_values))
chosen_predict = y_probs[instance]

waterfall_plot(chosen_instance, 
               shap_values_instance, 
               expected_value,
               chosen_predict, 
               max_display=15)<jupyter_output><empty_output><jupyter_text>KP still working on this: Creating a subplot of waterfallsChange the cell below from markdown to code to continue working on it fig = plt.figure(figsize=(10,5))expected_value = shapley_values_test_extended.base_values[0]explainer.expected_valueax1 = fig.add_subplot(131)instance = location_low_probabilitychosen_instance = X_test_std_df.iloc[instance]shap_values_instance = shapley_values_test[instance]chosen_predict = y_probs[instance]waterfall_plot(chosen_instance,                shap_values_instance,                expected_value,               chosen_predict,                max_display=15,               show=False)ax1.title.set_text('Low thrombolysis probability (displaying log-odds)')ax2 = fig.add_subplot(132)instance = location_middle_probabilitychosen_instance = X_test_std_df.iloc[instance]shap_values_instance = shapley_values_test[instance]chosen_predict = y_probs[instance]waterfall_plot(chosen_instance,                shap_values_instance,                expected_value,               chosen_predict,                max_display=15,               show=False)ax2.title.set_text('Medium thrombolysis probability (displaying log-odds)')ax3 = fig.add_subplot(133)instance = location_high_probabilitychosen_instance = X_test_std_df.iloc[instance]shap_values_instance = shapley_values_test[instance]chosen_predict = y_probs[instance]waterfall_plot(chosen_instance,                shap_values_instance,                expected_value,               chosen_predict,                max_display=15,               show=False)ax3.title.set_text('High thrombolysis probability (displaying log-odds)')fig.set_figheight(15)fig.set_figwidth(25)plt.tight_layout()plt.show()  Decision plots showing three instances on a subplot (probability) KP (need to check why the top value in the low probability chart is way off the the -ve, when should settle on 0). DIDN'T NEED TO WORK ON THIS AS FOUND CAN PASS AN ARGUMENT TO SHAP.DECISION_PLOT "Like the force plot, the decision plot supports link='logit' to transform log odds to probabilities."Kept the code here to mark what I was working on.<jupyter_code>fig = plt.figure(figsize=(10,5))

expected_value = shap_values_transformed.base_values

ax1 = fig.add_subplot(131)
instance = location_low_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shap_values_transformed[instance]
shap.decision_plot(expected_value, shap_values_instance.values, chosen_instance, show=False)
ax1.title.set_text('Low thrombolysis probability')

ax2 = fig.add_subplot(132)
instance = location_middle_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shap_values_transformed.values[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, show=False)
ax2.title.set_text('Medium thrombolysis probability')

ax3 = fig.add_subplot(133)
instance = location_high_probability
chosen_instance = X_test_std_df.iloc[instance]
shap_values_instance = shap_values_transformed.values[instance]
shap.decision_plot(expected_value, shap_values_instance, chosen_instance, show=False)
ax3.title.set_text('High thrombolysis probability')

#fig.set_figheight(15)
fig.set_figwidth(25)

plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0937531926486781, 'fraction_numerical': 0.014643615270225785, 'mean_word_length': 4.612512742099898, 'pattern_counts': {'":': 0, '<': 270, '<?xml version=': 0, '>': 271, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '66', 'dir': 'jupyter-structured-clean-dedup', 'id': '542793', 'license': 'MIT', 'n_tokens_mistral': 30265, 'n_tokens_neox': 28672, 'n_words': 7257, 'path': 'samuel_1/log_reg_shap/wip_logistic_regression_single_fit_shap-kfold3.ipynb', 'repo_name': 'samuel-book/samuel-2', 'size': '88093'}