<jupyter_start><jupyter_text>You are breaking my heartThis dataset contains information on 303 patients. Several medically relevant data are available (age, sex, cholesterol, resting blood pressure...). Our task is to predict the presence of heart disease (column "target", 0 means healty, 1 means sick).This dataset is described in detail:* on Kaggle datasets: https://www.kaggle.com/ronitf/heart-disease-uci* on its original webpage: https://archive.ics.uci.edu/ml/datasets/Heart+DiseaseI've downloaded a copy of the data and made it available at the following url:<jupyter_code>DATASET_URL = 'https://raw.githubusercontent.com/ne1s0n/coding_excercises/master/data/datasets_33180_43520_heart.csv'<jupyter_output><empty_output><jupyter_text>Battleplan* load dataset* a bit of data visualization/exploration* a baseline classifier: logistic regression* improve the data, improve the classifier* a better NN classifier: let's add a layer  * doing a proper crossvalidation  Config<jupyter_code>#let's fix already our desired number of epochs
EPOCHS = 200<jupyter_output><empty_output><jupyter_text>To be sure to have the same results, we can fix the random seeds.<jupyter_code>#general random seed
from numpy.random import seed
seed(10)

#tensorflow-specific seed. If you want to know why
#we need to set two separate random seeds see
#https://blog.cmgresearch.com/2020/09/04/tensorflow-has-to-random-seeds.html
import tensorflow
tensorflow.random.set_seed(10)<jupyter_output><empty_output><jupyter_text>Data<jupyter_code>import pandas

#pandas can read a csv directly from a url
heart_data = pandas.read_csv(DATASET_URL)
print(heart_data)
#splitting features and target
features = heart_data.iloc[:,:-1]
target = heart_data.iloc[:,-1]
#take a look at what we have done
print(heart_data.columns)
print(features.shape)
print(target.shape) #beware of rank 1 arrays<jupyter_output>Index(['age', 'sex', 'cp', 'trestbps', 'chol', 'fbs', 'restecg', 'thalach',
       'exang', 'oldpeak', 'slope', 'ca', 'thal', 'target'],
      dtype='object')
(303, 13)
(303,)<jupyter_text>Train and Validation sets<jupyter_code>#we want to have the same proportion of classes in both train and validation sets
from sklearn.model_selection import StratifiedShuffleSplit

#building a StratifiedShuffleSplit object (sss among friends) with 20% data
#assigned to validation set (here called "test")
sss = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=0)

#the .split() method returns (an iterable over) two lists which can be
#used to index the samples that go into train and validation sets
for train_index, val_index in sss.split(features, target):
    features_train = features.iloc[train_index, :]
    features_val   = features.iloc[val_index, :]
    target_train   = target[train_index]
    target_val     = target[val_index]
    
#let's print some shapes to get an idea of the resulting data structure
print(features_train.shape)
print(features_val.shape)
print(target_train.shape)
print(target_val.shape)<jupyter_output>(242, 13)
(61, 13)
(242,)
(61,)<jupyter_text>Baseline predictor: logistic regression<jupyter_code>#we are building a "sequential" model, meaning that the data will 
#flow like INPUT -> ELABORATION -> OUTPUT. In other words, there 
#we will not be any loops, i.e. our output will never be recycled as
#input for the first layer. That would make the architecture
#recurrent, which we don't want since the data is tabular and
#recurrent networks are more suited for streams and sequences of data
from keras.models import Sequential

#a "dense" layer is a layer were all the data coming in are connected
#to all nodes. In our case there is only one node in the layer, and
#it receives all the features.
#Doc: https://keras.io/api/layers/core_layers/dense/
from keras.layers import Dense

#this is optional: we'll add some extra metrics which are more human
#readable than crossentropy, so to have a practical feeling about
#how well the model is doing
import keras.metrics

# 2-class logistic regression in Keras
model = Sequential()
model.add(Dense(1, activation='sigmoid', input_dim=features_train.shape[1]))

#the model is declared, but we still need to compile it to actually
#build all the internal data structures. Moreover, we still need to 
#specify optimizer, loss function and, optionally, extra metrics
model.compile(optimizer='rmsprop', loss='binary_crossentropy', 
    metrics=[
      keras.metrics.BinaryAccuracy(),
      keras.metrics.AUC(name='auc'),
    ])<jupyter_output><empty_output><jupyter_text>At this point we are ready to train. You may have noticed that so far we just described the architecture. It's now time to have the network meed train and validation data.<jupyter_code>#doc: https://keras.io/api/models/model_training_apis/#fit-method
history = model.fit(
    x = features_train, 
    y = target_train, 
    epochs = 10, 
    validation_data = (features_val, target_val), 
    verbose=1)
#Taking a look to what's inside the returned history dictionary
print(history.history.keys())
#function to plot losses and metrics evolution
import matplotlib.pyplot as plt
def plot_loss_history(h, title):
  for metric in h.history.keys():
    #ignoring metrics on validation set, which are implied when
    #plotting on training set
    if metric.startswith('val_'):
      continue
    
    #if we get here we found a metric on the training set,
    #let's plot it
    plt.plot(h.history[metric], label = "Train set")
    plt.plot(h.history["val_" + metric], label = "Validation set")
    plt.xlabel('Epochs')
    plt.title(title + ' - ' + metric)
    plt.legend()
    plt.show()
plot_loss_history(history, 'Logistic (10 epochs)')
#putting verbose to 0 to avoid filling the screen
history2 = model.fit(features_train, target_train, epochs=(EPOCHS - 10), 
                     validation_data=(features_val, target_val), verbose=0)
#putting together the whole history
for k in history.history.keys():
  history.history[k] += history2.history[k]

#and plotting again
plot_loss_history(history, 'Logistic (' + str(EPOCHS) + ' epochs)')<jupyter_output><empty_output><jupyter_text>Taking a look to the [confusion matrix](https://www.analyticssteps.com/blogs/what-confusion-matrix)<jupyter_code>#sklearn provides a function to compute the confusion
#matrix given two numpy arrays: the true labels and
#the predicted labels
#doc: https://scikit-learn.org/stable/modules/generated/sklearn.metrics.confusion_matrix.html
from sklearn.metrics import confusion_matrix

#first: let's get the predictions for our validation set
predictions = model.predict(features_val)

#taking a look at the first few ones
print(predictions[1:5])<jupyter_output>[[0.89097434]
 [0.7931739 ]
 [0.9552659 ]
 [0.98716545]]<jupyter_text>At this stage predictions are numbers in the [0,1] range. That's expected, since it's the normal output of a sigmoid (a.k.a. a logistic function).To change them in labels (healty/sick) we need to discretyze. Everything above 0.5 goes to one, everything below goes to zero. In this way it will be comparable to true labels, which are naturally in the {zero, one} set.<jupyter_code>import numpy

#using numpy "where" function to do substitutions
#doc: https://numpy.org/doc/stable/reference/generated/numpy.where.html
predicted_labels = numpy.where(predictions > 0.5, 1, 0)

#computing and printing the matrix
con_mat_df = confusion_matrix(target_val, predicted_labels)
print(con_mat_df)<jupyter_output>[[21  7]
 [ 2 31]]<jupyter_text>Improvement: data normalization<jupyter_code>#getting an idea about features averages, sd
avg = features_train.mean()
std = features_train.std()
print('Feature means')
print(avg)
print('\nFeature standard deviations')
print(std)
#IMPORTANT: normalizing features using the same weights for both
#train and validation test (which are computed ON THE TRAIN SET)
features_train = (features_train - avg)/std
features_val = (features_val - avg)/std<jupyter_output><empty_output><jupyter_text>Improvement: class balancing<jupyter_code>#comparing the number of samples for each class
N_TOT = target_train.shape[0]
N_DISEASE = target_train.sum()
N_NORMAL = N_TOT - N_DISEASE

print('Total samples: ' + str(N_TOT) + ' (normal:' + str(N_NORMAL) + ' diseases:' + str(N_DISEASE) + ')')<jupyter_output>Total samples: 242 (normal:110 diseases:132)<jupyter_text>We propose a "in-house" class balancing function. There's a more complicated (and effective) one in sklearn called [compute_class_weight()](https://scikit-learn.org/stable/modules/generated/sklearn.utils.class_weight.compute_class_weight.html)<jupyter_code>#https://www.tensorflow.org/api_docs/python/tf/keras/Model#fit and
#whatch for argument class_weight

#handy function, to be reused later
def get_weights(target):
  #counting the number of instances for each class
  class_1_cnt = target_train.sum()
  class_0_cnt = target.shape[0] - target_train.sum()

  #we get weights as inverse of class ratios
  weight_for_0 = (class_0_cnt + class_1_cnt) / class_0_cnt
  weight_for_1 = (class_0_cnt + class_1_cnt) / class_1_cnt

  #and we are done
  return({0: weight_for_0, 1: weight_for_1})


#let's we invoke the function right away
class_weight = get_weights(target_train)
print('Computed weights:' + str(class_weight))<jupyter_output>Computed weights:{0: 2.2, 1: 1.8333333333333333}<jupyter_text>Improvement: better model<jupyter_code>#let's keep it simple: adding a Dense ReLU layer

# 2-class logistic regression in Keras
model2 = Sequential()
model2.add(Dense(10, activation='relu', input_dim=features_train.shape[1]))
model2.add(Dense(1, activation='sigmoid'))

#the model is declared, but we still need to compile it to actually
#build all the data structures
model2.compile(optimizer='rmsprop', loss='binary_crossentropy', 
    metrics=[
      keras.metrics.BinaryAccuracy(),
      keras.metrics.AUC(),
    ])<jupyter_output><empty_output><jupyter_text>Train again!<jupyter_code>#train with normalized data, class weights, improved model, same epochs
history_m2 = model2.fit(
    features_train, target_train, 
    epochs=EPOCHS, 
    validation_data=(features_val, target_val), 
    class_weight = class_weight,
    verbose=0)
#take a look at our results
plot_loss_history(history_m2, 'Improved NN')<jupyter_output><empty_output><jupyter_text>Confusion matrix for second model<jupyter_code>predictions = model2.predict(features_val)
predicted_labels = numpy.where(predictions > 0.5, 1, 0)
con_mat_df = confusion_matrix(target_val, predicted_labels)
print(con_mat_df)
#a direct comparison of loss functions
plt.subplot(1, 2, 1)
plt.plot(history.history['loss'], label='Model 1')
plt.plot(history_m2.history['loss'], label='Model 2')
plt.title('Train set loss')
plt.legend()
plt.xlabel('Epochs')

plt.subplot(1, 2, 2)
plt.plot(history.history['val_loss'], label='Model 1')
plt.plot(history_m2.history['val_loss'], label='Model 2')
plt.title('Validation set loss')
plt.legend()
plt.xlabel('Epochs')

plt.show()
#taking a look at final val accuracy

print('Binary accuracy on validation set')
print('Model 1) ' + str(history.history['val_binary_accuracy'][-1]))
print('Model 2) ' + str(history_m2.history['val_binary_accuracy'][-1]))

print('\nLoss on validation set')
print('Model 1) ' + str(history.history['val_loss'][-1]))
print('Model 2) ' + str(history_m2.history['val_loss'][-1]))<jupyter_output>Binary accuracy on validation set
Model 1) 0.8524590134620667
Model 2) 0.868852436542511

Loss on validation set
Model 1) 0.42718806862831116
Model 2) 0.42888498306274414<jupyter_text>Improvement: fine tuning number of units We want to add more units to our hidden layer. How many? This is our first true HYPERPARAMETER. To tune it we'll implement a five-folds crossvalidation scheme. BUT we'll need a new, pristine piece of data to actually evaluate the performances. It's time to introduce...  Train, validation and TEST sets!<jupyter_code>#support function to have a leaner code down below:
#input: train set, test set, number of units in the hidden layer
#output: train history object

def train_NN(feat_tr, feat_val, tar_tr, tar_val, n_units):
    #same 1 hidden layer model as above
    m = Sequential()
    m.add(Dense(units = n_units, activation='relu', input_dim=feat_tr.shape[1]))
    m.add(Dense(1, activation='sigmoid'))

    #the model is declared, but we still need to compile it to actually
    #build all the data structures
    m.compile(optimizer='rmsprop', loss='binary_crossentropy', 
        metrics=[
          keras.metrics.BinaryAccuracy(),
          keras.metrics.AUC(name='auc'),
        ])
    
    #compute class weights for this specific split
    cw = get_weights(tar_tr)

    #ready to train!
    h = m.fit(
      feat_tr, tar_tr, 
      epochs=EPOCHS, 
      validation_data=(feat_val, tar_val), 
      class_weight = cw,
      verbose=0)
    
    #and we are done
    return(h)
#we now further split the train set to do a proper crossvalidation using
#again scikit-learn, but this time we want the indexes for each split
#so we are going to instantiate a StratifiedKFold object
#https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.StratifiedKFold.html#sklearn.model_selection.StratifiedKFold
from sklearn.model_selection import StratifiedKFold

skf = StratifiedKFold(n_splits = 5)

#as done above, the .split() method returns (an iterable over) two lists which 
#can be used to index the samples that go into train and test sets

#a loop tracker, useful for indexing and printing messages
fold = 0

#room to store all the training histories
all_histories = {}

#let's explore these possible number of units
layer_units = [2, 4, 8, 16, 32]

#loop over folds
for train_index_cv, test_index_cv in skf.split(features_train, target_train):
    features_train_cv = features_train.iloc[train_index_cv, :]
    features_test_cv  = features_train.iloc[test_index_cv, :]
    target_train_cv   = target_train.iloc[train_index_cv]
    target_test_cv    = target_train.iloc[test_index_cv]

    #user interface
    fold += 1
    print('Doing fold ' + str(fold))

    #room to store all the training histories
    all_histories[fold] = {}

    #loop over considered number of units
    for lu in layer_units:
      #user interface
      print(' - training with lu=' + str(lu))

      #training the network, storing the training history
      all_histories[fold][lu] = train_NN(
          features_train_cv, features_test_cv, 
          target_train_cv, target_test_cv, lu)

#all_histories is a bit messy, order of indexing is:
#fold -> number of units -> actual history object
#let's take a look at all these indexes:
print(all_histories.keys())
print(all_histories[1].keys())
print(all_histories[1][16].history.keys())
#for each considered number of units we have five different executions
#let's take a look at one
lu = 16
plt.plot(all_histories[1][lu].history['val_loss'], label = 'Fold 1')
plt.plot(all_histories[2][lu].history['val_loss'], label = 'Fold 2')
plt.plot(all_histories[3][lu].history['val_loss'], label = 'Fold 3')
plt.plot(all_histories[4][lu].history['val_loss'], label = 'Fold 4')
plt.plot(all_histories[5][lu].history['val_loss'], label = 'Fold 5')
plt.xlabel('Epochs')
plt.title('Test set Loss for ' + str(lu) + ' units')
plt.legend()
plt.show()
#let's forget history and focus on metrics for last epoch
#also: let's average over folds

rows = []
for lu in layer_units:
  for fold in all_histories.keys():
    row = []
    row.append(lu)
    row.append(all_histories[fold][lu].history['val_loss'][-1])
    row.append(all_histories[fold][lu].history['val_binary_accuracy'][-1])
    row.append(all_histories[fold][lu].history['val_auc'][-1])
    rows.append(row)

#converting to pandas
df = pandas.DataFrame(rows, columns=["LU", "loss", "accuracy", "AUC"])

#average over folds
df = df.groupby(df['LU']).aggregate('mean')

#and take a look to the numbers
print(df)
#selecting LU with best loss
LU_selected = df.loc[:, 'loss'].idxmin()
print('Best performance with ' + str(LU_selected) + ' units in the hidden layer')
#train again, this time using the full train set, and evaluating in the validation set
history_m3 = train_NN(
    features_train, features_val,
    target_train, target_val,
    LU_selected)

#taking a look at the final validation loss, even if the comparison it's a 
#bit unfair (can you tell why?)
print('Model 1) ' + str(history.history['val_loss'][-1]))
print('Model 1) ' + str(history_m2.history['val_loss'][-1]))
print('Model 1) ' + str(history_m3.history['val_loss'][-1]))<jupyter_output>0.42718806862831116
0.42888498306274414
0.3832351863384247
{'fraction_non_alphanumeric': 0.10778804906047973, 'fraction_numerical': 0.023019757114373755, 'mean_word_length': 5.178424785367675, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 63, 'https://': 12, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '232039', 'license': 'MIT', 'n_tokens_mistral': 5669, 'n_tokens_neox': 5276, 'n_words': 1635, 'path': 'lab_day3/day3_code01 heart disease crossv.ipynb', 'repo_name': 'ne1s0n/coding_excercises', 'size': '16551'}