model_notebook.py

# -*- coding: utf-8 -*-
"""explainable-risk-scoring-for-chronic-disease.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1lzMGz8RCUrmCLxPmcW9wzLpXXDOAkKbO

## Imports
"""

# Commented out IPython magic to ensure Python compatibility.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import gdown
# %matplotlib inline
# import warnings

# sklearn preprocessing imports
from sklearn.preprocessing import OneHotEncoder, StandardScaler, OrdinalEncoder
from sklearn.impute import KNNImputer
from sklearn.model_selection import train_test_split, GridSearchCV, RandomizedSearchCV , cross_val_score
from sklearn import metrics

from statsmodels.stats.outliers_influence import variance_inflation_factor

#sklearn model imports
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import AdaBoostClassifier, RandomForestClassifier,RandomForestRegressor,GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score , confusion_matrix , classification_report ,mean_absolute_error ,  r2_score , root_mean_squared_error



from sklearn.naive_bayes import GaussianNB
from xgboost import XGBClassifier

# warnings.filterwarnings("ignore")
sns.set_style('darkgrid')

# Constants

RANDOM_STATE = 32

"""## Load Dataset"""

url = 'https://drive.google.com/drive/folders/1o-9kSIeL52_OFId4OmjyicYTBg9eW9O_'
gdown.download_folder(url, quiet=True, use_cookies=False)

df = pd.read_csv('dataset/heart_disease_uci.csv')

"""## EDA

## Feature description

- **id**: (Unique id for each patient)
- **age**: (Age of the patient in years)
- **origin**: (place of study)
- **sex**: (Male/Female)
- **cp**: chest pain type ([typical angina, atypical angina, - non-anginal, asymptomatic])
- **trestbps**: resting blood pressure (resting blood - pressure (in mm Hg on admission to the hospital))
- **chol**: (serum cholesterol in mg/dl)
- **fbs**: (if fasting blood sugar > 120 mg/dl)
- **restecg**: (resting electrocardiographic results) -- Values: [normal, stt abnormality, lv hypertrophy]
- **thalach**: maximum heart rate achieved
- **exang**: exercise-induced angina (True/ False)
- **oldpeak**: ST depression induced by exercise relative to - rest
- **slope**: the slope of the peak exercise ST segment
- **ca**: number of major vessels (0-3) colored by - fluoroscopy
- **thal**: [normal; fixed defect; reversible defect]
- **num**: the predicted attribute [0=no heart disease; 1,2,3,4 = stages of heart disease ]
"""

df

print('Females:')
total_females = len(df.loc[(df['sex'] == 'Female')])
print(f"Total Females: {total_females}")

undiseased_females = len(df.loc[(df['sex'] == 'Female') & (df['num'] != 0)])
print(f"Undiseased Females: {undiseased_females}")


print('\n\nMales:')
total_male = len(df.loc[(df['sex'] == 'Male')])
print(f"Total Male: {total_male}")

undiseased_male = len(df.loc[(df['sex'] == 'Male') & (df['num'] != 0)])
print(f"Undiseased Male: {undiseased_male}")

df.isna().sum()

"""**Above we have lot's of missing values:**
- Either we need to use KNN to impute these values

- Once the model is created we may consider to remove (if it doesn't impact the model negatively) the columns or impute the values with KNN
"""

# df.boxplot(column=['age', 'trestbps', 'chol', 'thalch', 'oldpeak', 'ca', 'num'])
df.boxplot(column=['age', 'trestbps', 'chol', 'thalch', 'oldpeak', 'ca', 'num'])

plt.figure(figsize=(10, 6))
sns.boxplot(
    data=df,
    x='cp',
    y='num'
)

# display(df.info())

cat_columns = ['sex', 'dataset', 'cp', 'fbs', 'restecg', 'exang', 'slope', 'thal']

print('Unique items in column: ')

for column in cat_columns:
  print(f'{column} : {df[column].nunique()}')

sns.pairplot(
    df[[
        "ca",
        "oldpeak",
        "thalch",
        "chol",
        "age",
        "trestbps",
        "num"
    ]],
    diag_kind="kde"
)
plt.show()

plt.figure(figsize=(8,6))
plt.title("Age vs Disease Stage")
sns.barplot(x = 'num', y = 'age', data = df)
plt.ylabel('Age')
plt.xlabel('Stage')

plt.figure(figsize=(8,6))
plt.title("Cholestrol vs Disease Stage")
sns.barplot(x = 'num', y = 'chol', data = df)
plt.ylabel('Cholestrol')
plt.xlabel('Stage')

plt.figure(figsize=(8,6))
plt.title("Max Heart Rate vs Disease Stage")
sns.barplot(x = 'num', y = 'thalch', data = df)
plt.ylabel('Max Heart Rate')
plt.xlabel('Stage')

plt.figure(figsize=(8,6))
plt.title("Chest Pain vs Disease Stage")
sns.barplot(x = 'cp', y = 'num', data = df)
plt.xlabel('Chest Pain')
plt.ylabel('Stage')

plt.figure(figsize=(8,6))
plt.title("ca: number of major vessels vs Disease Stage")
sns.barplot(x = 'num', y = 'ca', data = df)
plt.ylabel('ca: number of major vessels')
plt.xlabel('Stage')

plt.figure(figsize=(8,6))
plt.title("Blood Pressure vs Disease Stage")
sns.barplot(x = 'num', y = 'trestbps', data = df)
plt.ylabel('Blood Pressure')
plt.xlabel('Stage')

plt.figure(figsize=(10, 8))
sns.heatmap(df.corr(numeric_only=True), annot=True, cmap='coolwarm', fmt=".2f")
plt.show()

df[df.duplicated(keep=False)]

"""No duplicates in the dataframe

## Data transformation
"""

# df.drop(columns=['thal'], inplace=True, axis=1)

df[ df['chol'] == 0]

"""### Encoding data (One Hot Encoding)"""

encoded_df = pd.get_dummies(df, columns=cat_columns, drop_first=True)

print(f'Origin columns length: {len(df.columns)}, encoded length: {len(encoded_df.columns)}')
encoded_df

encoded_df.drop(columns='id', inplace=True)

X = encoded_df
Y = X.pop('num')

# from statsmodels.tools.tools import add_constant
# X_with_const = add_constant(X)

# X_with_const.columns



# vif_data = pd.DataFrame()
# vif_data["feature"] = X_with_const.columns
# vif_data["VIF"] = [variance_inflation_factor(X_with_const.values, i)
#                    for i in range(len(X_with_const.columns))]

# display(X)
# display(Y)

# temp_df = df.copy()

# if('trestbps' in df.columns):
#   df['bp_category'] = pd.cut(df['trestbps'],
#                                    bins=[0, 120, 130, 300],
#                                    labels=[0, 1, 2]).astype(float)
#   df.drop('trestbps', inplace=True, axis=1)

# df

"""Removed binning, since it did not improved the model any further, and we may miss some important information as well

### Outliers Treatment
"""

df.boxplot(column=['trestbps', 'chol', 'thalch'])

"""It is impossible for a living person to have 0 cholestrol or blood pressure, so I am replacing those values with nan and later on it will be imputed with KNN"""

# df['chol'][df['chol'] == 0]
# df['trestbps'][df['trestbps'] == 0]

df['chol'] = df['chol'].replace(0, np.nan)
df['trestbps'] = df['trestbps'].replace(0, np.nan)

df.boxplot(column=['trestbps', 'chol', 'thalch'])

"""- Cholestrol: Now, since a normal human being can have 400 or even higher cholestrol levels, we are not capping them.


"""

df.boxplot(column=['oldpeak', 'ca'])

upper_limit = df['oldpeak'].quantile(0.99)
df['oldpeak'] = np.where(df['oldpeak'] > upper_limit, upper_limit, df['oldpeak'])

"""Capping oldpeak (ST Depression) to 99 percentile of the oldpeak data to avoid outliers"""

df.boxplot(column=['oldpeak', 'ca'])

"""### Train Test splitting"""

# Data split for train and test:
# Train: 80%, Test: 20%
x_train, x_test, y_train, y_test = train_test_split(X, Y, test_size=0.2, random_state=RANDOM_STATE)

"""### Scaling"""

scaler = StandardScaler()

x_train_scaled = scaler.fit_transform(x_train)

x_test_scaled = scaler.transform(x_test)

np.sum(np.isnan(x_train_scaled))

"""### Missing value using KNN"""

# Imputed missing values with KNN imputer

knn_imputer = KNNImputer(n_neighbors=2)

x_train_imputed = knn_imputer.fit_transform(x_train_scaled)

x_test_imputed = knn_imputer.transform(x_test_scaled)

"""## Model Selection
For this dataset I would select multiple models and from them I would pick the best performing one
- Logistic Regression (OVR)
- Random Forest
- Gradient Boosting
- XGBoost
- KNN
- Naive Baise

### Logistic Regression
"""

LR_MODEL_TEST = LogisticRegression(random_state=RANDOM_STATE, solver='liblinear', max_iter=1000, multi_class='ovr')

C = [0.001, 0.01, 0.1, 1, 10]

cv_folds = 5

param_grid = {
  'C': C,
  'l1_ratio': [0, 0.1, 0.2, 0.3, 0.8, 0.9, 1]
}

gs_model = GridSearchCV(estimator=LR_MODEL_TEST, param_grid=param_grid, cv=5, scoring='accuracy', refit=True, return_train_score = True)

lr_model = gs_model.fit(x_train_imputed, y_train)

cv_results = pd.DataFrame(lr_model.cv_results_)
cv_results.columns

# Getting the important columns out
cv_results_important = cv_results[['param_C', 'param_l1_ratio',
                                   'mean_train_score','std_train_score',
                                   'mean_test_score','std_test_score',
                                   'rank_test_score']]
cv_results_important.sort_values('rank_test_score').head()

print("Best Param: ", lr_model.best_params_)
print("Best Score: ", lr_model.best_score_)

lr_best_C = lr_model.best_params_.get('C')
lr_best_l1_ratio = lr_model.best_params_.get('l1_ratio')

# Final Model with the best possible params

lr_model = LogisticRegression(random_state=RANDOM_STATE, solver='liblinear', max_iter=1000, multi_class='ovr', C=lr_best_C, l1_ratio=lr_best_l1_ratio)

lr_model.fit(x_train_imputed, y_train)

# Train and Test Accuracy

lr_predict_prob_train = lr_model.predict_proba(x_train_imputed)

lr_predict_prob_test = lr_model.predict_proba(x_test_imputed)

print("Train Data Accurancy = ",round(100*metrics.roc_auc_score(y_train, lr_predict_prob_train, multi_class='ovr'),4))
print("Test Data Accurancy = ",round(100*metrics.roc_auc_score(y_test, lr_predict_prob_test, multi_class='ovr'),4))

"""### Random Forest"""

rf_model_test = RandomForestClassifier(n_estimators=100, random_state=RANDOM_STATE)


cv_folds = 5

param_grid = {
  'n_estimators': [100, 140],
  'min_samples_leaf': [2, 3, 4],
  'max_depth': [4, 6, 10]
}

gs_model = GridSearchCV(estimator=rf_model_test, param_grid=param_grid, cv=5, scoring='accuracy', refit=True, return_train_score = True)

rf_model_final = gs_model.fit(x_train_imputed, y_train)

cv_results = pd.DataFrame(rf_model_final.cv_results_)

# display(cv_results.columns)

# Getting the important columns out
cv_results_important = cv_results[['param_min_samples_leaf', 'param_max_depth','param_n_estimators', 'mean_train_score','std_train_score', 'mean_test_score','std_test_score', 'rank_test_score']]

cv_results_important.sort_values('rank_test_score').head(10)

print("Best Param: ", rf_model_final.best_params_)
print("Best Score: ", rf_model_final.best_score_)

rf_model = RandomForestClassifier(n_estimators=140, min_samples_leaf=3, max_depth=6, random_state=RANDOM_STATE)

rf_model.fit(x_train_imputed, y_train)

rf_predicted_train = rf_model.predict_proba(x_train_imputed)
rf_predicted_test = rf_model.predict_proba(x_test_imputed)

print("Train Data Accurancy = ",round(100*metrics.roc_auc_score(y_train, rf_predicted_train, multi_class='ovr'),4))
print("Test Data Accurancy = ",round(100*metrics.roc_auc_score(y_test, rf_predicted_test, multi_class='ovr'),4))

"""We got above accuracy with feature encoding and it's performing very badly, but since for trees algorithms categorical feature encoding is not mandatory, we will try again without feature encoding."""

df

""" **Since the Random forest model was not giving optimal result because of overfitting, let's try to optimize it for RF by imputing missing values and keeping the categorical columns**"""

sample_df = df.copy()

sample_df['ca'].fillna(-1, inplace=True)
sample_df['oldpeak'].fillna(-1, inplace=True)
# sample_df['trestbps'].fillna(-1, inplace=True)
sample_df['chol'].fillna(-1, inplace=True)
sample_df['fbs'].fillna(-1, inplace=True)
sample_df['thalch'].fillna(-1, inplace=True)
sample_df['slope'].fillna('missing', inplace=True)
sample_df['thal'].fillna('missing', inplace=True)
# sample_df['restecg'].fillna('missing', inplace=True)

sample_df = sample_df.dropna(subset=['exang', 'restecg'])


X_fr = sample_df
Y_fr = X_fr.pop('num')

X_fr.drop(columns='id', inplace=True)

# display(X_fr)

# display(X_fr.isna().sum())

x_train_fr, x_test_fr, y_train_fr, y_test_fr = train_test_split(
    X_fr, Y_fr, test_size=0.2, random_state=RANDOM_STATE
)

# categorical coummns
categorical_cols = ['sex', 'dataset', 'cp', 'fbs', 'restecg', 'exang', 'slope', 'thal']

encoder = OrdinalEncoder(handle_unknown='use_encoded_value', unknown_value=-1)

x_train_fr[categorical_cols] = encoder.fit_transform(x_train_fr[categorical_cols])
x_test_fr[categorical_cols] = encoder.transform(x_test_fr[categorical_cols])

# reconfirm to check if all cols are int
x_train_fr = x_train_fr.apply(pd.to_numeric, errors='coerce')
x_test_fr = x_test_fr.apply(pd.to_numeric, errors='coerce')


rf_model = RandomForestClassifier(n_estimators=120, min_samples_leaf=2, max_depth=6, random_state=RANDOM_STATE)

rf_model.fit(x_train_fr, y_train_fr)

rf_predicted_train = rf_model.predict_proba(x_train_fr)
rf_predicted_test = rf_model.predict_proba(x_test_fr)

print("Train Data Accurancy = ",round(100*metrics.roc_auc_score(y_train_fr, rf_predicted_train, multi_class='ovr'),4))
print("Test Data Accurancy = ",round(100*metrics.roc_auc_score(y_test_fr, rf_predicted_test, multi_class='ovr'),4))

"""We see a improvement of around 5-6% after training the model without onehot encoding, but still need improvement

### XGBoost
"""

xg_model_test = XGBClassifier(
    objective='multi:softprob', # for multi-class probability prediction
    n_estimators=100,           # number of boosting rounds (trees)
    learning_rate=0.1,          # step size shrinkage
    max_depth=3,                # maximum depth of a tree
    seed=27
  )

C = [0.001, 0.01, 0.1, 1, 10]

cv_folds = 5

param_grid = {
  'n_estimators': [100, 140],
  # 'min_samples_leaf': [2, 3, 4],
  # 'max_depth': [4, 5, 6, 8, 10, 12]
}

xg_model = GridSearchCV(estimator=xg_model_test, param_grid=param_grid, cv=5, scoring='accuracy', refit=True, return_train_score = True)

xg_model_final = xg_model.fit(x_train_imputed, y_train)

# display(cv_results = pd.DataFrame(xg_model_final.cv_results_))
cv_results.columns

# Getting the important columns out
cv_results_important = cv_results[['param_max_depth', 'param_min_samples_leaf',
                                   'mean_train_score','std_train_score',
                                   'mean_test_score','std_test_score',
                                   'rank_test_score']]
cv_results_important.sort_values('rank_test_score').head()

print("Best Param: ", rf_model_final.best_params_)
print("Best Score: ", rf_model_final.best_score_)

# Final Model with the best possible params

xg_model = XGBClassifier(
    objective='multi:softprob', # for multi-class probability prediction
    n_estimators=100,           # number of boosting rounds (trees)
    learning_rate=0.1,          # step size shrinkage
    max_depth=3,                # maximum depth of a tree
    seed=27
  )

xg_model.fit(x_train_imputed, y_train)

# Train and Test Accuracy

xg_predict_prob_train = xg_model.predict_proba(x_train_imputed)

xg_predict_prob_test = xg_model.predict_proba(x_test_imputed)

print("Train Data Accurancy = ",round(100*metrics.roc_auc_score(y_train, xg_predict_prob_train, multi_class='ovr'),4))
print("Test Data Accurancy = ",round(100*metrics.roc_auc_score(y_test, xg_predict_prob_test, multi_class='ovr'),4))

"""This model is also not performing well, there is a huge gap between train and test ROC accuracy

### Combined Models

Let's create combined report for all the classification models, and then select highest performing one.
"""

from sklearn.pipeline import Pipeline

import numpy as np, random
from sklearn.model_selection import StratifiedKFold

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

model_configs = [
    {
        "name": "LogisticRegression",
        "pipeline": Pipeline([
            ('scaler', StandardScaler()),
            ('imputer', KNNImputer(n_neighbors=5)),
            ('model', LogisticRegression(multi_class='ovr', random_state=RANDOM_STATE, max_iter=1000))
        ]),
        "params": {
            "model__C": [0.01, 0.1, 1, 10],
            "model__solver": ["liblinear", "lbfgs"]
        }
    },
    {
        "name": "RandomForest",
        "pipeline": Pipeline([
            ('scaler', StandardScaler()),
            ('imputer', KNNImputer(n_neighbors=5)),
            ('model', RandomForestClassifier(random_state=RANDOM_STATE))
        ]),
        "params": {
            "model__n_estimators": [100, 140, 200],
            "model__max_depth": [4, 6, 8],
            "model__min_samples_leaf": [3, 5]
        }
    },
    {
        "name": "GradientBoosting",
        "pipeline": Pipeline([
            ('scaler', StandardScaler()),
            ('imputer', KNNImputer(n_neighbors=5)),
            ('model', GradientBoostingClassifier(random_state=RANDOM_STATE))
        ]),
        "params": {
              "model__n_estimators": [100, 200],
              "model__learning_rate": [0.05, 0.1],
              "model__max_depth": [3, 5]
        }
    },
    {
        "name": "AdaBoost",
        "pipeline": Pipeline([
            ('scaler', StandardScaler()),
            ('imputer', KNNImputer(n_neighbors=5)),
            ('model', AdaBoostClassifier(random_state=RANDOM_STATE))
        ]),
        "params": {
              "model__n_estimators": [50,100,200],
              "model__learning_rate": [0.05,0.1,1]
        }
    },
    {
        "name": "KNeighbors",
        "pipeline": Pipeline([
            ('scaler', StandardScaler()),
            ('imputer', KNNImputer(n_neighbors=5)),
            ('model', KNeighborsClassifier())
        ]),
        "params": {
        "model__n_neighbors": [3,5,7,9],
        "model__weights": ["uniform", "distance"]
    }
    }
]


def train_and_evaluate(X_train, y_train, X_test, y_test, cv_strategy):
    # model_configs = model_configs
    results = []
    best_overall_model = None
    best_roc_auc = 0.0
    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

    for config in model_configs:
        print(f"Processing: {config['name']}...")

        # Initialize GridSearchCV
        grid_search = GridSearchCV(
            estimator=config['pipeline'],
            param_grid=config['params'],
            cv=cv_strategy,
            scoring='accuracy',
            n_jobs=-1
        )

        # Fit using raw data (Pipeline handles the rest)
        grid_search.fit(X_train, y_train)

        # Get the best estimator (the full pipeline)
        best_pipe = grid_search.best_estimator_

        # Calculate ROC-AUC scores
        prob_train = best_pipe.predict_proba(X_train)
        prob_test = best_pipe.predict_proba(X_test)

        roc_train = metrics.roc_auc_score(y_train, prob_train, multi_class='ovr')
        roc_test = metrics.roc_auc_score(y_test, prob_test, multi_class='ovr')

        results.append({
            "Model": config['name'],
            "Best_Params": grid_search.best_params_,
            "CV_Accuracy": grid_search.best_score_,
            "ROC_Train": roc_train * 100,
            "ROC_Test": roc_test * 100,
            "Pipeline_Object": best_pipe
        })

        # Track the best model based on Test ROC-AUC
        if (roc_test * 100) > best_roc_auc:
            best_roc_auc = roc_test * 100
            best_overall_model = best_pipe

    # Sort results for easier reading
    results_df = pd.DataFrame(results).sort_values(by="ROC_Test", ascending=False)
    return results_df, best_overall_model

# Make sure to use the RAW X_train/X_test from your split, not the imputed ones!
summary_report, best_heart_model = train_and_evaluate(x_train, y_train, x_test, y_test, cv)

"""### Models accuracy report"""

# display(summary_report[['Model', 'CV_Accuracy', 'ROC_Train', 'ROC_Test']])

plt.figure(figsize=(10, 8))
sns.barplot(x='Model', y='ROC_Train', data=summary_report)
sns.barplot(x='Model', y='ROC_Test', data=summary_report)