A Practical Approach to Building and Evaluating Machine Learning Models

# Eg. importing libraries (For Machine Learning, you will need additional libraries)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import sys
color_pal = sns.color_palette()


from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, make_scorer
from sklearn.model_selection import train_test_split
from xgboost import XGBRegressor
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler, OrdinalEncoder
from sklearn.kernel_ridge import KernelRidge
from sklearn.svm import SVR
from sklearn import model_selection
from sklearn.model_selection import GridSearchCV


print(f"Python Version: {sys.version}")
print(f"Numpy Version: {np.__version__}")
print(f"Pandas Version: {pd.__version__}")
print(f"Seaborn Version: {sns.__version__}")
color_pal

Python Version: 3.11.9 (main, Apr 19 2024, 16:48:06) [GCC 11.2.0]
Numpy Version: 1.26.4
Pandas Version: 2.2.2
Seaborn Version: 0.13.2

# Load data
df = pd.read_csv('./insurance.csv')
df.head()

# Check shape
# We have 1338 rows and 7 columns
df.shape

(1338, 7)

# Missing values. Count null values
df.isna().sum()

age         0
sex         0
bmi         0
children    0
smoker      0
region      0
charges     0
dtype: int64

# Count duplicates
df.duplicated().sum()

1

# View duplicates. Keep =False will include the first occurrence.
# I will not remove this as it's only 1 and could be a legitimate value.
df[df.duplicated(keep=False)]

# Check the datatypes
df.dtypes

age           int64
sex          object
bmi         float64
children      int64
smoker       object
region       object
charges     float64
dtype: object

# More detailed information
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1338 entries, 0 to 1337
Data columns (total 7 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   age       1338 non-null   int64  
 1   sex       1338 non-null   object 
 2   bmi       1338 non-null   float64
 3   children  1338 non-null   int64  
 4   smoker    1338 non-null   object 
 5   region    1338 non-null   object 
 6   charges   1338 non-null   float64
dtypes: float64(2), int64(2), object(3)
memory usage: 73.3+ KB

# View Statistics for numerical data. 
df.describe()

# view Box plots for numeric data. Charges range is too large, so we need to plot each box plot separately.
sns.boxplot(df)
plt.show()

# Create Boxplots

color_pal = sns.color_palette()
# Set figure size
sns.set(rc={'figure.figsize': (12.0, 8.0)})

columns = ['age', 'bmi', 'children', 'charges']

# Number of plots and setup grid layout
n_cols = 2
n_rows = (len(columns) + n_cols - 1) // n_cols  # Calculate rows needed = (4 + 2 - 1) // 2 = 2

fig, axes = plt.subplots(n_rows, n_cols, sharex=False, figsize=(12, 8))

#print(axes)
# Flatten axes array for easy iteration
axes = axes.flatten()

# Create boxplots for each column
for i, column in enumerate(columns):
    sns.boxplot(y=df[column], ax=axes[i], color=color_pal[i])
    axes[i].set_title(column)

plt.suptitle('Box Plot for Insurance Data', fontsize=16)

# Delete any unused subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

# We can use the same idea for the plots above and create histogram with kde plot to see distribution
# This is some repeating code, so we can put this in a function to be able to create each. I will leave as is for now for simplicity.
color_pal = sns.color_palette()
# Set figure size
sns.set(rc={'figure.figsize': (12.0, 8.0)})

columns = ['age', 'bmi', 'children', 'charges']

# Number of plots and setup grid layout
n_cols = 2
n_rows = (len(columns) + n_cols - 1) // n_cols  # Calculate rows needed = (4 + 2 - 1) // 2 = 2

fig, axes = plt.subplots(n_rows, n_cols, sharex=False, figsize=(12, 8))

#print(axes)
# Flatten axes array for easy iteration
axes = axes.flatten()

# Create boxplots for each column
for i, column in enumerate(columns):
    sns.histplot(df[column], kde=True, bins=10, ax=axes[i], color=color_pal[i])
    axes[i].set_title(column)

plt.suptitle('Box Plot for Insurance Data', fontsize=16)

# Delete any unused subplots if any
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

# Let's do a pairplot and color by sex
sns.pairplot(df, hue='sex')
plt.show()

# View correlation with heatmap
# Calculate Correlation first
df_corr = df[columns].corr()
df_corr

# Use heatmap to visualize correlation
sns.heatmap(df_corr, annot=True)
plt.show()

from scipy import stats
# Apply Pearson Correlation coefficient to see if there is any statistical significance between bmi and charges and age and charges
# HO: There is no correlation (correlation coefficient is 0)
# H1: There is a non zero correlation

pearson_coef, p_value = stats.pearsonr(df['bmi'], df['charges'])
print("The Pearson Correlation Coefficient is", pearson_coef, " with a P-value of P =", p_value)

# Make decision
alpha = 0.005

if p_value <= alpha:
    print("Reject HO: There is a non zero correlation")
else:
    print("Fail to reject HO: There is no correlation")

The Pearson Correlation Coefficient is 0.19834096883362887  with a P-value of P = 2.459085535116766e-13
Reject HO: There is a non zero correlation

# Check if you need new features that you can get externally or combining existing features
# For simplicity we will add an extra feature age_range, that adds categorical data for specified age range
# 0-18: Bin 0 (Under 18)
# 18-30: Bin 1
# 31-50: Bin 2
# 51-70: Bin 3
# Above 70
# Define age bins and labels
bins = [0, 18, 30, 50, 70, float('inf')]  # Age ranges. Use float('inf') so that we don't have to manually set upperbound
labels = [0, 1, 2, 3, 4]  # Corresponding labels
df['age_range'] = pd.cut(df['age'], bins=bins, labels=labels, right=False ) # right=False to get distinct ranges, no overlaps
df.head()

# Create variables for our features and Target
FEATURES = ['age', 'sex', 'bmi', 'children', 'smoker', 'region', 'age_range']
TARGET = 'charges'

X = df[FEATURES]
y = df[TARGET]

# Split into Features and Target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# View the size of splits
print(f'Data Shape: {df.shape}')
print(f'Train Data Shape: {X_train.shape}')
print(f'Test Data Shape: {X_test.shape}')

Data Shape: (1338, 8)
Train Data Shape: (1070, 7)
Test Data Shape: (268, 7)

# Normalize and Encode the data
preprocessor = ColumnTransformer(
    transformers = [
        ('num', StandardScaler(), ['age', 'bmi', 'children']),
        ('cat_nom', OneHotEncoder(), ['sex', 'smoker', 'region']),
        ('cat_ord', OrdinalEncoder(), ['age_range'])
    ],
    remainder='passthrough'
)

# Add models to test
models = []
models.append(('LR', LinearRegression()))
models.append(('KR', KernelRidge()))
models.append(('SVR', SVR()))

# Using KFold to evaluate models
def RMSE(y_true, y_pred):
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    return rmse

def my_scorer():
    return make_scorer(RMSE, greater_is_better=False) # If using R2, you would need to change this to True

results = {}

for name, model in models:
    # Create a pipeline that first applies preprocessing and then fits the model
    pipeline = Pipeline(steps=[
        ('preprocessor', preprocessor),
        ('model', model)
    ])

    kfold = model_selection.KFold(n_splits=10, random_state=7, shuffle=True)
    score_results = model_selection.cross_val_score(pipeline, X_train, y_train, cv=kfold, scoring=my_scorer())
    results[name] = abs(score_results.mean())

for key, value in results.items():
    print(f"{key} => {value}")

# Sort and print the models by their MSE performance (lower MSE is better)
sorted_results = dict(sorted(results.items(), key=lambda item: item[1]))
print("\nSorted Models (Lower RMSE is Better):")
for key, value in sorted_results.items():
    print(f"{key} => MSE: {value:.6f}")

LR => 6135.912678882409
KR => 6136.46623117466
SVR => 12558.969016331775

Sorted Models (Lower RMSE is Better):
LR => MSE: 6135.912679
KR => MSE: 6136.466231
SVR => MSE: 12558.969016

# Ignore warnings from here. This will take a while
import warnings
warnings.filterwarnings('ignore')

# LR and KR are close. We don't have hyperparameters to finetune for LR
# Use GridSearch with KR and see if we can finetune the parameters and then compare both
pipeline = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('model', KernelRidge())
])

# Hyperparameter grid for the KernelRidge model.
# Look up the documentation to decide on parameters.
param_grid = {
    'model__alpha': [1.0, 2.0, 2.5, 1.5],
    'model__kernel': ['linear', 'sigmoid', 'rbf', 'poly'],
    'model__gamma': [None, 1, 0.1],
    'model__degree': [1, 2, 3],
    'model__coef0': [0.0, 0.1, 0.5]
}

gridSearch = GridSearchCV(pipeline, param_grid, cv=10, scoring=my_scorer())
gridSearch.fit(X_train, y_train)
# Look for the best score
best_score = abs(gridSearch.best_score_) # Since my_scorer() uses greater_is_better=False. The negative version is being presented here
print(f"The best score is: {best_score}") 
best_params = gridSearch.best_params_
print(f"The best params are: {best_params}") 
# Select the best estimator with best params for Kernel Ridge. This will be used below for model evaluation
kr = gridSearch.best_estimator_

The best score is: 4928.332175217688
The best params are: {'model__alpha': 2.5, 'model__coef0': 0.5, 'model__degree': 2, 'model__gamma': 1, 'model__kernel': 'poly'}

# We have to fit the LR model first as we didn't do that above. The Kernel Ridgle kr was already created above
pipeline_lr = Pipeline(steps=[
     ('preprocessor', preprocessor),
     ('model', LinearRegression())
])
pipeline_lr.fit(X_train, y_train)

# Make predictions with LR and KR
y_pred_lr = pipeline_lr.predict(X_test)
y_pred_kr = kr.predict(X_test)

print(f"Score for LR: {RMSE(y_test, y_pred_lr)}")
print(f"Score for KR: {RMSE(y_test, y_pred_kr)}")

Score for LR: 5801.634782167118
Score for KR: 4531.34843081154

# Let us visualize the Actual Values vs Predicted for LR and use a line of equality to see how close the predictions are to this line.
# Define a function so we don't end up duplicating code
def plot_scatter(y_test, y_pred, model_name):
    plt.figure(figsize=(10, 8))
    plt.scatter(y_test, y_pred)
    plt.xlabel('Actual Values')
    plt.ylabel('Predicted Values')
    plt.title(f'Actual vs Predicted for {model_name} ')
    plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], color='red', linestyle='--')
    plt.show()

# Plot for Linear Regression
plot_scatter(y_test, y_pred_lr, "Linear Regression")

# Plot for Kernel Ridge
plot_scatter(y_test, y_pred_kr, "Kernel Ridge")

# The plots can be also lined up beside each other using subplots.

fig, axes = plt.subplots(1, 2, figsize=(12, 8))
axes[0].scatter(y_test, y_pred_lr)
axes[0].set_xlabel('Actual Values')
axes[0].set_ylabel('Predicted Values')
axes[0].set_title(f'Actual vs Predicted for Linear Regression')
axes[0].plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], color='red', linestyle='--')
axes[1].scatter(y_test, y_pred_kr)
axes[1].set_xlabel('Actual Values')
axes[1].set_ylabel('Predicted Values')
axes[1].set_title(f'Actual vs Predicted for Kernel Ridge')
axes[1].plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], color='red', linestyle='--')
plt.show()

# Save model
import pickle

with open('kr_model.pkl', 'wb') as file:
    pickle.dump(kr, file)