The Behavioral Risk Factor Surveillance System (BRFSS) is a health-related telephone survey that is collected annually by the CDC. Each year, the survey collects responses from over 400,000 Americans on health-related risk behaviors, chronic health conditions, and the use of preventative services. It has been conducted every year since 1984. For this project, a csv of the dataset available on Kaggle for the year 2015 was used. This original dataset contains responses from 441,455 individuals and has 330 features. These features are either questions directly asked of participants, or calculated variables based on individual participant responses.

# import models
from random_forest import RandomForest
from logistic_regression import LogisticReg
from gradient_boosting import GradientBoosting

from xg_boost import ModelXGBoost
from data_cleaning import DataCleaning
from sklearn.model_selection import train_test_split
import pandas as pd
import numpy as np
import sklearn.datasets
import matplotlib.pyplot as plt
import shap
from ipywidgets import FloatProgress
import matplotlib.pyplot as plt
import seaborn as sns

Read data

data = pd.read_csv("data/diabetes.csv")

categorical_variables = ['GenHlth', 'Age', 'Education', 'Income']
data[categorical_variables] = data[categorical_variables].apply(lambda x: x.astype('category'))

display(data.head(3))

	Diabetes_binary	HighBP	HighChol	CholCheck	BMI	Smoker	Stroke	HeartDiseaseorAttack	PhysActivity	Fruits	...	AnyHealthcare	NoDocbcCost	GenHlth	MentHlth	PhysHlth	DiffWalk	Sex	Age	Education	Income
0	0.0	1.0	0.0	1.0	26.0	0.0	0.0	0.0	1.0	0.0	...	1.0	0.0	3.0	5.0	30.0	0.0	1.0	4.0	6.0	8.0
1	0.0	1.0	1.0	1.0	26.0	1.0	1.0	0.0	0.0	1.0	...	1.0	0.0	3.0	0.0	0.0	0.0	1.0	12.0	6.0	8.0
2	0.0	0.0	0.0	1.0	26.0	0.0	0.0	0.0	1.0	1.0	...	1.0	0.0	1.0	0.0	10.0	0.0	1.0	13.0	6.0	8.0

3 rows × 22 columns

Exploratory Data Analysis

data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 70692 entries, 0 to 70691
Data columns (total 22 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   Diabetes_binary       70692 non-null  float64 
 1   HighBP                70692 non-null  float64 
 2   HighChol              70692 non-null  float64 
 3   CholCheck             70692 non-null  float64 
 4   BMI                   70692 non-null  float64 
 5   Smoker                70692 non-null  float64 
 6   Stroke                70692 non-null  float64 
 7   HeartDiseaseorAttack  70692 non-null  float64 
 8   PhysActivity          70692 non-null  float64 
 9   Fruits                70692 non-null  float64 
 10  Veggies               70692 non-null  float64 
 11  HvyAlcoholConsump     70692 non-null  float64 
 12  AnyHealthcare         70692 non-null  float64 
 13  NoDocbcCost           70692 non-null  float64 
 14  GenHlth               70692 non-null  category
 15  MentHlth              70692 non-null  float64 
 16  PhysHlth              70692 non-null  float64 
 17  DiffWalk              70692 non-null  float64 
 18  Sex                   70692 non-null  float64 
 19  Age                   70692 non-null  category
 20  Education             70692 non-null  category
 21  Income                70692 non-null  category
dtypes: category(4), float64(18)
memory usage: 10.0 MB

Definition of Categorical Variables

Diabetes_binary (scale: 0-2)
0 = no diabetes
1 = prediabetes
2 = diabetes

Age (scale: 1-13)
1 = 18-24; 2 = 25-29; 3 = 30-34; 4 = 35-39; 5 = 40-44; 6 = 45-49; 7 = 50-55; 8 = 56-59; 9 = 60-64; 10 = 65-69; 11 = 70-74; 12 = 75-79; 13 = 80 or older

GenHlth (scale: 1-5) 1 = excellent
2 = very good
3 = good
4 = fair
5 = poor

Education (scale: 1-6)
1 = Never attended school or only kindergarten
2 = Grades 1 through 8 (Elementary)
3 = Grades 9 through 11 (Some high school)
4 = Grade 12 or GED (High school graduate)
5 = College 1 year to 3 years (Some college or technical school)
6 = College 4 years or more (College graduate)

Income (scale: 1-8) 1 = less than 10K
2 = less than 25K
3 = less than 35K
4 = less than 45K
5 = less than 55K
6 = less than 65K
7 = less than 75K
8 = $75K or more

Exploratory Data Analysis

Distribution of Age, Income, Education

plt.figure(figsize=(10, 5))
plt.subplot(2,2,1)
sns.histplot(data['Age']).set(title='Distribution of Age')

plt.subplot(2,2,2)
sns.histplot(data['Income']).set(title='Distribution of Income')

plt.subplot(2,2,3)
sns.histplot(data['Education']).set(title='Distribution of Education')

plt.subplot(2,2,4)
sns.histplot(data['GenHlth']).set(title='GenHlth')

plt.subplots_adjust(left=0,
                    bottom=0, 
                    right=0.9, 
                    top=0.9,
                    wspace=0.5, 
                    hspace=0.5)
plt.show()

As shown on the histograms above, the distribution for Age and GenHlth are the closest to a normal distribution.

Boxplot of BMI

sns.boxplot(data['BMI'])

<Axes: ylabel='BMI'>

As shown on the boxplot above, there are outliers in BMI

Remove Outliers

obj = DataCleaning()
data = obj.replace_outliers_cols(df=data, cols=['BMI'])

2181 outliers were replaced for the column: BMI

Boxplot of BMI after replacing outliers with the upperbound or lowerbound values

sns.boxplot(data['BMI'])

<Axes: ylabel='BMI'>

Check for missing values

# Check for missing values in each column
print(f"The number of NaNs in each column: {data.isnull().sum()}")

The number of NaNs in each column: Diabetes_binary         0
HighBP                  0
HighChol                0
CholCheck               0
BMI                     0
Smoker                  0
Stroke                  0
HeartDiseaseorAttack    0
PhysActivity            0
Fruits                  0
Veggies                 0
HvyAlcoholConsump       0
AnyHealthcare           0
NoDocbcCost             0
GenHlth                 0
MentHlth                0
PhysHlth                0
DiffWalk                0
Sex                     0
Age                     0
Education               0
Income                  0
dtype: int64

Split data into training and testing

# split data frames into features and target
X =  data.drop(columns='Diabetes_binary')
y = data[['Diabetes_binary']]

# split to training and testing
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=100)

score_tracker = {}

Random Forest

rf = RandomForest()

# Train the model on the train set
model = rf.classifier(X=X_train, y=y_train)

# Training accuracy
prediction = rf.predict(model_object=model, X=X_train)
train_accuracy = rf.accuracy(prediction_object=prediction, y=y_train)
#print(f"Training Accuracy is: {train_accuracy}")

# Testing Accuracy
prediction = rf.predict(model_object=model, X=X_test)
test_accuracy = rf.accuracy(prediction_object=prediction, y=y_test)
#print(f"Testing Accuracy is: {test_accuracy}")


# Tune the model
tuner = rf.tune_classifier(model_object=model, X=X_train, y=y_train)
# Retrieve best parameters
best_params = rf.best_params(tuned_model_object=tuner)
# Retrieve best scores
best_score = rf.best_score(tuned_model_object=tuner)
print(f"Best params are === {best_params}")
#print(f"Best score is === {best_score}")

# Train model with the best hyper-parameters, returns model object, training prediction and testing prediction
tuned_model = rf.tuned_classifier(X=X_train, y=y_train, bestParams=tuner.best_params_)

# calculate accuracy for test set
prediction = rf.predict(model_object=tuned_model, X=X_test)
test_accuracy = rf.accuracy(prediction_object=prediction, y=y_test)
#print(f"Testing Accuracy after tuning: {test_accuracy}")

imp = rf.feature_importance(model_object=tuned_model)
sorted_imp = rf.sorted_feature_importance_indicies(model_object=tuned_model)
print(f"Index of Sorted Feature Importance === {sorted_imp}")

important_features = []
col_names = list(X_train.columns)

for x in sorted_imp:
    important_features.append(col_names[x])

#print(f"Sorted features importance === {important_features}")

most_important_feature = list(X_train.iloc[:, [sorted_imp[0]]].columns)[0]
print(f"Most important feature is === {most_important_feature}")

# cross-validation - AUC
cross_val_auc = rf.cross_validation_classifier(tuned_model, X_train, y_train)

print(f"Mean cross-validation AUC === {cross_val_auc.mean()}")

# append score to dictionary
score_tracker['Random Forest'] = cross_val_auc.mean()

Best params are === {'criterion': 'entropy', 'max_depth': 10, 'max_features': 6, 'min_samples_leaf': 2, 'min_samples_split': 2, 'n_estimators': 510}
Index of Sorted Feature Importance === [13  0  3 18  1 16 20 15  6 14 19 17 10  2  7  4  8  9  5 12 11]
Most important feature is === GenHlth
Mean cross-validation AUC === 0.8287197394531376

GBM

gb = GradientBoosting()

# Train the model on the train set
model = gb.classifier(X=X_train, y=y_train)
# Training accuracy
prediction_train = gb.predict(model_object=model, X=X_train)
train_accuracy = gb.accuracy(prediction_object=prediction_train, y=y_train)
#print(f"Training Accuracy is: {train_accuracy}")

# Testing accuracy
prediction_test = gb.predict(model_object=model, X=X_test)
test_accuracy = gb.accuracy(prediction_object=prediction_test, y=y_test)
#print(f"Testing Accuracy is: {test_accuracy}")

#tune model
tuner = gb.tune_classifier(model_object=model, X=X_train, y=y_train)
# Retrieve best parameters
best_params = gb.best_params(tuned_model_object=tuner)
# Retrieve best scores
best_score = gb.best_score(tuned_model_object=tuner)

print(f"Best params are === {best_params}")
#print(f"Best score is === {best_score}")

tuned_model = gb.tuned_classifier(X=X_train, y=y_train, bestParams=best_params)

# Testing accuracy
prediction_test = gb.predict(model_object=tuned_model, X=X_test)
test_accuracy = gb.accuracy(prediction_object=prediction_test, y=y_test)
#print(f"Testing Accuracy after tuning is: {test_accuracy}")

# feature importance
imp = gb.feature_importance(model_object=tuned_model)
sorted_imp = gb.sorted_feature_importance_indicies(model_object=tuned_model)
#print(f"Index of Sorted Feature Importance === {sorted_imp}")

important_features = []
col_names = list(X_train.columns)

for x in sorted_imp:
    important_features.append(col_names[x])

print(f"Sorted features importance === {important_features}")


most_important_feature = list(X_train.iloc[:, [sorted_imp[0]]].columns)[0]
print(f"Most important feature is === {most_important_feature}")

# cross-validation - AUC
cross_val_auc = gb.cross_validation_classifier(tuned_model, X_train, y_train)

print(f"Mean cross-validation AUC === {cross_val_auc.mean()}")

# append score to dictionary
score_tracker['GBM'] = cross_val_auc.mean()

Best params are === {'criterion': 'squared_error', 'learning_rate': 0.001, 'max_depth': 90, 'max_features': 2, 'min_samples_leaf': 3, 'min_samples_split': 2, 'n_estimators': 360}
Sorted features importance === ['BMI', 'GenHlth', 'Age', 'HighBP', 'HighChol', 'Income', 'PhysHlth', 'DiffWalk', 'Education', 'MentHlth', 'HeartDiseaseorAttack', 'PhysActivity', 'Sex', 'Smoker', 'Fruits', 'Veggies', 'HvyAlcoholConsump', 'CholCheck', 'Stroke', 'NoDocbcCost', 'AnyHealthcare']
Most important feature is === BMI
Mean cross-validation AUC === 0.8263080767962554

Logistic Regression

lr = LogisticReg()

#train the model on the training set
lr_model = lr.classifier(X=X_train, y=y_train)

# predict binary
prediction = lr.predict(model_object=lr_model, X=X_train)
# Accuracy training set
accuracy_train = lr.accuracy(model_object=lr_model, X=X_train, y=y_train)
#print(f"Training Accuracy is: {accuracy_train}")
# Accuracy testing set
accuracy_test = lr.accuracy(model_object=lr_model, X=X_test, y=y_test)
#print(f"Testing Accuracy is: {accuracy_test}")

# tune the model
tuner = lr.tune(model_object=lr_model, X=X_train, y=y_train)
# Retrieve best parameters
best_params = lr.best_params(tuned_model_object=tuner)
# Retrieve best scores
best_score = lr.best_score(tuned_model_object=tuner)
print(f"Best params are === {best_params}")
#print(f"Best score is === {best_score}")

tuned_model = lr.tuned_model(X=X_train, y=y_train, bestParams=tuner.best_params_)
# calculate accuracy for test set
prediction = lr.predict(model_object=tuned_model, X=X_test)
test_accuracy = lr.accuracy(model_object=tuned_model, X=X_test, y=y_test)
#print(f"Testing Accuracy after tuning: {test_accuracy}")

feature_importance, sorted_idx, pos = lr.feature_importance(model_object=tuned_model)

# cross-validation - AUC
cross_val_auc = lr.cross_validation(tuned_model, X_train, y_train)

print(f"Mean cross-validation AUC === {cross_val_auc.mean()}")

# append score to dictionary
score_tracker['Logistic Regression'] = cross_val_auc.mean()

Best params are === {'solver': 'newton-cg', 'max_iter': 100, 'C': 100}
Mean cross-validation AUC === 0.8271197031081078

XGBoost

xg = ModelXGBoost()

# Train the model on the train set
model = xg.classifier(X=X_train, y=y_train)
# Training accuracy
prediction_train = xg.predict(model_object=model, X=X_train)
train_accuracy = xg.accuracy(prediction_object=prediction_train, y=y_train)
#print(f"Training Accuracy is: {train_accuracy}")

# Testing accuracy
prediction_test = xg.predict(model_object=model, X=X_test)
test_accuracy = xg.accuracy(prediction_object=prediction_test, y=y_test)
#print(f"Testing Accuracy is: {test_accuracy}")

#tune model
tuner = xg.tune_classifier(model_object=model, X=X_train, y=y_train)
# Retrieve best parameters
best_params = xg.best_params(tuned_model_object=tuner)
# Retrieve best scores
best_score = xg.best_score(tuned_model_object=tuner)

print(f"Best params are === {best_params}")
#print(f"Best score is === {best_score}")

tuned_model = xg.tuned_classifier(X=X_train, y=y_train, bestParams=best_params)

# Testing accuracy
prediction_test = xg.predict(model_object=tuned_model, X=X_test)
test_accuracy = xg.accuracy(prediction_object=prediction_test, y=y_test)
#print(f"Testing Accuracy after tuning is: {test_accuracy}")

# feature importance
imp = xg.feature_importance(model_object=tuned_model)
sorted_imp = xg.sorted_feature_importance_indicies(model_object=tuned_model)
#print(f"Index of Sorted Feature Importance === {sorted_imp}")

important_features = []
col_names = list(X_train.columns)

for x in sorted_imp:
    important_features.append(col_names[x])

#print(f"Sorted features importance === {important_features}")


most_important_feature = list(X_train.iloc[:, [sorted_imp[0]]].columns)[0]
print(f"Most important feature is === {most_important_feature}")

# cross-validation - AUC
cross_val_auc = xg.cross_validation_classifier(tuned_model, X_train, y_train)

print(f"Mean cross-validation AUC === {cross_val_auc.mean()}")

# append score to dictionary
score_tracker['XGBoost'] = cross_val_auc.mean()

Best params are === {'subsample': 0.8, 'reg_alpha': 0, 'n_estimators': 160, 'min_child_weight': 4, 'max_depth': 2, 'learning_rate': 0.05, 'gamma': 0.4, 'colsample_bytree': 0.7}
Most important feature is === HighBP
Mean cross-validation AUC === 0.8293038596929987

Best Performing model

best_model = max(score_tracker.items(), key=lambda k: k[1])
print(f"The best model is: {best_model}")

The best model is: ('XGBoost', 0.8293038596929987)

The best performing model was XGBoost with a AUC of 0.8293

SHAP Values from XGBoost

# shap values
shap_values, X = xg.shap_values(model_object=tuned_model, X=X_train)
shap.summary_plot(shap_values, X, plot_type="bar")

The plot above shows the plot of the SHAP values. This shows the most important features that the model used to make predictions. The most important features in determining a person has diabetes is GenHlth(General Health), BMI, HighBP (High blood pressure), Age, and HighChol (High Cholestrol)

Confusion Matrix for XGBoost model

from sklearn.metrics import confusion_matrix
from sklearn import metrics

y_pred = prediction_test
y_true = y_test
confusion_matrix = confusion_matrix(y_true, y_pred)

cm_display = metrics.ConfusionMatrixDisplay(confusion_matrix = confusion_matrix, display_labels = [False, True])

cm_display.plot()
plt.show()