Exploring Machine Learning Models: From Iris Classification to Diabetes Regression Analysis

Preface

In an ever-evolving field like machine learning, the pursuit of understanding and mastery is a continuous journey. This comprehensive guide is a testament to that journey, drawing inspiration from Mark Fenner’s accessible work “Machine Learning for Everyone,” yet going beyond to rework, enhance, and expand these concepts to suit a broader understanding.

Through tireless research, collaboration with intelligent tools such as ChatGPT, and leveraging my own domain knowledge, I’ve crafted a narrative that delves into the intricacies of machine learning. From classification tasks with the Iris dataset to in-depth linear regression analyses on the diabetes dataset, this guide is a thorough exploration of essential and advanced concepts.

Featuring 44 meticulously crafted images, each available for closer inspection, the visuals serve to clarify and enhance the textual content. Click to zoom in, navigate through the arrows, and immerse yourself in a detailed understanding of the subject matter.

A labor of love, this guide represents countless hours of study, writing, and editing. While I’ve done my utmost to ensure accuracy and clarity, the process of learning is an iterative one, and your engagement and feedback are highly valued. At an estimated read time of 153 minutes (approximately 2.55 hours), this guide is both a deep dive and a valuable resource.

Whether you are a seasoned professional or a curious enthusiast, I invite you to embark on this exploratory journey through the world of machine learning. Your understanding, insights, and innovations await.

Introduction: A Comprehensive Exploration of Machine Learning Models

1. Starting with a Classic: Iris Dataset: We initiate our exploration with the well-known Iris dataset, a starting point for many in the world of machine learning. By applying different classification algorithms, we set the stage for more complex analyses, building foundational knowledge, and understanding.

2. Main Focus: Diabetes Dataset: The core of our investigation lies in the exploration of the diabetes dataset. This investigation leads us into the world of regression modeling, emphasizing the following key aspects:

• Linear Regression Analysis: Simple yet powerful, linear regression serves as our benchmark, demonstrating strength in its simplicity.
• Ensemble Methods: A deeper dive into more sophisticated ensemble methods allows us to explore complexity and the potential for overfitting.
• Model Metrics: A thorough understanding of Mean Absolute Error (MAE), Mean Squared Error (MSE), Median Absolute Error, and R-squared, guides our evaluation process.
• Hyperparameter Tuning: By leveraging nested cross-validation, we optimize model performance.
• Normalization: We delve into the impact of data scaling and normalization, investigating how preprocessing affects model efficiency and interpretability.
• Visualizations and Comparative Analysis: Utilizing various visualization techniques, including residual plots and comparative error bar plots, we gain visual insights into model performance.

3. Comprehensive Approach: This investigation is more than a series of individual analyses; it’s a holistic exploration of machine learning models, where classification and regression, simple models and complex ensembles, model evaluation, and hyperparameter tuning all come together.

4. A Guided Journey: Join me as I embark on this multifaceted journey, where I start with the basics and move towards more intricate analyses. My comprehensive examination of the diabetes dataset, supplemented by the Iris dataset’s exploration, presents an opportunity to engage with the art and science of machine learning. From understanding key metrics like MAE, MSE, and R-squared to exploring model complexity and optimization, this study offers a glimpse into the dynamic and nuanced field of data analytics.

Understanding Errors in Iris Flower Classification

In predictive modeling for flower classification, particularly using the Iris dataset, understanding the various types of errors is vital. The dataset consists of 150 samples from three species of Iris flowers, and four features are recorded for each sample: the lengths and widths of the sepals and petals. From these features, we aim to classify the species of the Iris flower.

In a flower, the sepal and petal are both essential parts but they serve different functions.

1. Sepal: The sepals are the small, leaf-like parts growing at the base of the petals. They are typically green and resemble leaves. Their primary function is to protect the flower before it blooms. In the bud stage, the sepals enclose the petals and keep them safe from insects, weather, and other potential threats.

2. Petal: The petals are often the most noticeable part of a flower due to their bright and varied colors. They serve two main functions. Firstly, their vibrant colors and patterns are designed to attract pollinators (like bees, butterflies, and birds). Secondly, some petals have lines called nectar guides that guide pollinators to the nectar and pollen at the center of the flower.

In the context of the Iris dataset, these terms refer to the length and the width of the petal and sepal, which are four key features used to identify different types of Iris flowers. These features are continuous numerical values, which makes this dataset suitable for various types of machine learning algorithms, especially classification algorithms.

Sources of Errors

Below is a summary of potential sources of errors in our Iris flower classification model:

• Measurement Errors: Inaccuracies in the measurements of sepal and petal sizes can affect classification.

• Mismatch Between Model and Reality: If the assumed mathematical relationship between inputs (sepal and petal dimensions) and outputs (species) doesn’t match reality, predictions will be inaccurate.

Understanding Data Variance in Predictive Modeling

Inaccuracies in predictions, like incorrect species classification, may stem from various factors, including inherent randomness between input features and target output, over which we may have limited control.

Real-World Randomness in Iris Classification

Two Iris flowers with the same sepal and petal dimensions might not necessarily be of the same species. This inherent randomness is a crucial distinction between mathematical functions and real-world data modeling.

Unraveling Random Processes: Variance in Iris Classification

When classifying an Iris flower based on its features, the output is not a single value but a probability distribution over the possible species, reflecting the inherent variability or “variance” in the data.

Taming the Variance in Our Models

In machine learning systems, certain sources of error can be controlled. The choice of model, training/testing data split, and handling of outliers are examples of areas where we can manage variance.

Unmasking Bias in Our Models

The model’s alignment with the relationship between inputs and outputs is another source of error. High-bias models may fail to capture complex patterns, leading to systematic errors, while low-bias models can follow intricate patterns, potentially providing more accurate classification.

Assembling the Pieces: The Bias-Variance Decomposition

1. The bias of our model,
2. The variability introduced during the model training process, and
3. The inherent variability in our data.

This reordering corresponds to the components in your equation:

$ \text{Error} = \text{Bias}_{\text{Learner}} + \text{Variance}_{\text{Learner}}(\text{Training}) + \text{Variance}_{\text{Data}} $

1. $\text{Bias}_{\text{Learner}}$ corresponds to the bias of our model.
2. $\text{Variance}_{\text{Learner}}(\text{Training})$ corresponds to the variability introduced during the model training process.
3. $\text{Variance}_{\text{Data}}$ corresponds to the inherent variability in our data.

By understanding the individual contributions of bias and variance in our models, we can better diagnose potential challenges and optimize performance, particularly in complex tasks like Iris classification.

Interactions between Data and Learner

The interaction between training data and learning algorithm is essential. If the training data doesn’t represent the full diversity of Iris flowers, our model may not generalize well.

In conclusion, understanding these sources of error is essential for building robust Iris flower classification models. By focusing on these aspects, researchers can develop models that classify Iris species accurately and reliably.

Variations in Data Relationships

The accuracy of our model’s classification ability could be affected if the data has inconsistencies, such as errors in measuring the sepal and petal lengths. Moreover, the relationships between the different data features and the chosen classification algorithm’s mathematical functions play a significant role in classification accuracy.

Delving into the Complexity of Classification

The concept of classification in machine learning may seem intuitive, but there’s a great depth and complexity behind it. When working with a dataset like the Iris dataset, where the objective is to classify flowers into different species based on their features, we need to consider several aspects:

1. Establishing the Boundaries

In the context of the Iris dataset, boundaries help to distinguish between the different species based on their petal length and petal width. These boundaries could be linear or nonlinear, depending on the algorithm used. Is there a distinct demarcation between species, or do they overlap? The nature of these boundaries is often defined by the algorithm and its underlying mathematics.

2. Evaluation of the Boundaries

The effectiveness of the classification boundaries needs to be evaluated. How well do the boundaries separate the different species? Do species overlap in certain areas? If so, is this due to the nature of the data or the algorithm’s limitations?

The Mathematical Underpinnings of Classification

Each classification algorithm, such as decision trees, support vector classifiers, and logistic regression, has unique mathematical principles that guide how it establishes and evaluates boundaries. These principles determine how well the algorithm can segment the Iris species based on petal length and petal width.

Making Decisions with Established Boundaries

With the Iris dataset, classification involves questions like:

• Is this flower’s petal length and petal width indicative of Iris setosa?
• Is Iris virginica more likely than Iris versicolor, given the petal measurements?
• Among the three species, which one is most probable for this flower?

Broadening Our Classification Toolbox: A Case Study with the Iris Dataset

In the pursuit of classifying Iris species, we’ll begin with the technique of decision trees. Through this example, we understand the power and simplicity of decision trees. Such models are highly valuable in various fields and applications, turning the challenging task of classification into a series of manageable and logical questions.

Constructing a Decision Tree: An Introduction

Among many available machine learning models, the decision tree stands out for its interpretability and efficiency. Let’s explore how to classify Iris flowers based on petal length and petal width using a decision tree.

Our decision tree might start with the root question: “Is the petal length less than 2.5 cm?”. Depending on the response, the tree probes further, asking more specific questions until it reaches a prediction for the species of the Iris.

Tree-Building Algorithms

There are a number of major tree-building algorithms. ID3, C4.5, and C5.0 were made by Quinlan. CART was developed independently.

Steps in Tree-Building

In general, the tree-building algorithms use the following steps:

1. Evaluate the set of features and splits and pick a “best” feature-and-split.
2. Add a node to the tree that represents the feature-split.
3. For each descendant, work with the matching data and either:
   1. If the targets are similar enough, return a predicted target.
   2. If not, return to step 1 and repeat.

Constraints in Tree-Building

The algorithms control what splits and partitions are allowed, how feature-splits are evaluated, what makes targets in a group similar enough to form a leaf, and other limits.

The constraints commonly include an absolute limit on the depth of the tree and the minimum number of examples at a leaf to make a prediction, regardless of similarity. These constraints help prevent overfitting.

Imperfections and Overfitting in Decision Trees

While constructing a decision tree offers many advantages, we must be aware of potential pitfalls:

1. Inherent randomness: External factors not included in the model may affect the classification of the species. This randomness can introduce noise into our predictions.
2. Measurement errors: Any inaccuracies in measuring petal lengths or widths could lead to incorrect predictions. Precise measurements are crucial to building a reliable model.
3. Overfitting: Decision trees may perform exceptionally well on training data but fail to generalize to unseen data. Care must be taken to avoid overfitting by potentially limiting the tree’s depth or employing pruning techniques.

By understanding these challenges, we can take proactive measures to mitigate them, ensuring a more robust and accurate decision tree model.

Preventing Overfitting in Decision Trees

When dealing with decision trees, the model can become excessively complex, leading to overfitting, if it is left unconstrained. Overfitting results in high variance and poor generalization capabilities to new data. However, we can adopt several strategies to prevent overfitting, hence reducing the variance:

1. Limiting Tree Depth: By setting a maximum depth for the tree, we limit the number of questions the model can ask before making a prediction. This constraint simplifies the model, making it less likely to overfit the training data.

2. Setting a Minimum Number of Samples at Leaves: By requiring a certain number of examples in each leaf node, we force the model to make broader generalizations. This approach can prevent overfitting by not allowing the model to create rules for very specific or rare cases present in the training data.

3. Limiting the Number of Considered Features: By constraining the number of features the model considers when making splits, we not only simplify the model but also speed up the training process. This method can prevent the model from relying too heavily on potentially noisy or less relevant features.

These strategies introduce some bias into the model (making it less flexible and more general) in order to reduce its variance (improving its generalization capabilities). This approach reflects the trade-off between bias and variance that is inherent to all machine learning models.

How Decision Trees Work

The tree-building algorithm partitions the entire space of data (e.g., a plot with petal length on the x-axis and petal width on the y-axis) into smaller regions. It repeats this process until meeting a specific stopping criterion. This criterion might be reaching a maximum depth or a point where further splits contain too few examples to be meaningful.

Mathematically, the partitioning can be expressed as:

$$\text{target} = \sum_{R \in P} C_R I(\text{example} \in R)$$

Explaining the Equation

• $R$ is a region of our graph
• $P$ is the group of all of the regions
• $C_R$ is the predicted class for the region $R$ (for example, one of the types of irises)
• $I$ is an indicator function which gives us a value of one if an example is in $R$ and zero otherwise

The indicator function and the predicted class together give us $C_R$ when an example is in $R$ and give us zero everywhere else.

Each region mostly contains points (flowers) of the same class, allowing for meaningful separations between classes.

In summary, tree-building algorithms work by recursively partitioning the data into smaller regions, creating a decision tree where each path from the root to a leaf node corresponds to a specific sequence of decisions that lead to a class prediction. Different decision tree methods can vary in the way they perform these partitions and in their stopping criteria, focusing on aspects like petal length and petal width to create regions with meaningful separations between classes.

Understanding Decision Tree Depth

The depth of a decision tree is a crucial factor that significantly influences its performance and complexity.

The Role of max_depth

The max_depth parameter in the DecisionTreeClassifier determines the maximum depth to which the tree can grow, thereby limiting the number of splits it can make.

Observing Decision Trees at Different Depths

By experimenting with trees of depths 1, 2, and 3, we can witness the changes in decision boundaries and complexity.

• max_depth=1: At this level, the tree is allowed to make just one split, bifurcating the feature space into two regions.

• max_depth=2: This allows three splits, one for the root node, and another split for each of its two child nodes, leading to a total of four regions in the feature space.

The Default Setting: max_depth=None

In the default setting, max_depth=None, the tree has the freedom to grow indefinitely. The depth is only constrained by the number of samples in the leaves or by meeting a particular stopping condition. This behaves as if max_depth is set to infinity, allowing the tree to make unlimited splits until it perfectly classifies the training data or reaches a stopping condition.

Understanding Impurity Measures in Decision Trees: Gini Value, Entropy, and Information Gain

The Gini index and entropy are both measures used in the construction of decision trees to evaluate the “quality” of a potential split. The goal of both measures is to quantify the uncertainty or disorder within a set of items. In the context of decision trees, these items are the data instances within a node.

The quality of the splits in a decision tree is usually evaluated using measures such as the Gini Index, Entropy, and Information Gain.

• Gini Index: This measure gives an indication of the impurity of a dataset. A Gini Index of 0 denotes perfect purity, where all items belong to a single class, while a Gini index of 0.5 in a binary classification problem indicates maximum impurity.

• Entropy: Entropy measures the disorder in the dataset. An entropy value of 0 means the set contains instances from a single class, whereas a value of 1 indicates that the instances are evenly split between two classes in a binary classification problem.

• Information Gain: Information gain is based on entropy. It measures the reduction in entropy achieved by a split. The attribute with the highest information gain is chosen as the splitting attribute at the node.

In the context of Extreme Random Forests, which we get to later, these measures help to evaluate the quality of the randomly selected split points.

Gini Value and Gini Index

Gini Value measures the impurity of a node in a decision tree. A value of 0 indicates perfect purity, while a value of 0.5 in binary classification indicates maximum impurity.

Gini Index quantifies the likelihood of misclassification of a randomly chosen element from the dataset. It’s expressed using the formula:

$$ Gini(p) = 1 - \sum (p_i)^2 $$

where $p_i$ is the probability of an object being classified to a particular class.

Entropy

Entropy measures the information disorder within a set of items. It’s 0 when all instances belong to one class and 1 when instances are evenly split between two classes in binary classification.

The formula for entropy is:

$$ Entropy(p) = - \sum p_i \log_2(p_i) $$

where $ p_i $ is the proportion of instances that are of class i.

Information Gain

This concept is based on entropy and signifies the reduction in entropy due to a split. It helps in choosing the splitting attribute at a node and measures the relative change in entropy with respect to the independent variables.

Comparative Insights

Both the Gini Index and Entropy are potent measures of impurity and are often interchangeable. Gini Index is generally faster to compute, making it preferable for large datasets. In contrast, Entropy can create a more balanced tree and may be more reliable with noisy data.

However, it’s crucial to recognize that these measures alone cannot prevent overfitting. Proper hyperparameter tuning and tree-pruning strategies must be employed to build a robust decision tree model.

The choice between Gini and Entropy may not significantly affect a decision tree’s performance, but understanding both concepts is valuable for various applications in machine learning and data analysis.

# Import necessary libraries
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import cross_val_score
from sklearn.datasets import load_iris

# Load iris dataset
iris = load_iris()

# Initialize Decision Tree Classifier with default parameters
decision_tree_classifier = DecisionTreeClassifier()

# Define the cross-validation strategy
# In this case, it's k-fold cross-validation with k=3
cv_strategy = 3

# Define the scoring metric for evaluation, in this case, accuracy
scoring_metric = 'accuracy'

print(f"Performing {cv_strategy}-Fold Cross Validation...")

# Perform cross-validation and store the scores
cv_scores = cross_val_score(decision_tree_classifier,
                                 iris.data,
                                 iris.target,
                                 cv=cv_strategy,
                                 scoring=scoring_metric)

print(f"Cross validation scores: {cv_scores}")

# Compute the mean accuracy and print
mean_accuracy = cv_scores.mean()
print(f"Mean accuracy: {mean_accuracy}")

Performing 3-Fold Cross Validation...
Cross validation scores: [0.98 0.92 0.98]
Mean accuracy: 0.96

Decision Boundaries

Decision Trees create blocky, regions of colored space from our features.

Forming the Boundaries

The partition, or the shapes of the regions, is formed from overlaid rectangles where only the topmost rectangle counts. Imagine neatly spreading out a few decks of cards on a rectangular table. The cards can overlap. When cards of the same suit are touching one another, they form a larger region for that suit.

To create these regions from points in space, we create simple yes/no answers by thresholding. For example, we test a feature Temp against a value 55. Choosing the values to split on is the trickiest aspect of implementing decision trees by hand.

plot_decision_boundary Function

The following code is responsible for visualizing the decision boundaries formed by the decision tree classifier. By plotting the boundaries in the feature space (e.g., petal length vs. petal width), we can observe how the decision tree algorithm partitions the data into different regions, each corresponding to a particular class. This graphical representation provides intuitive insights into the model’s decision-making process, allowing us to understand how the features are being utilized to distinguish between the different classes. It also aids in identifying potential areas of overfitting or underfitting, enhancing our ability to interpret and fine-tune the model.

def plot_decision_boundary(ax, data, tgt, model, dims, class_labels, colors, grid_step=.01):
    """
    This function plots the decision boundary of a classification model.

    Parameters:
    ax: matplotlib Axes object
    data: 2D numpy array
    tgt: 1D numpy array
    model: sklearn model object
    dims: list of two indices
    class_labels: dict of class numbers mapped to class names
    colors: list of colors for each class
    grid_step: float (default is 0.01)
    """
    # Select two dimensions from data
    twoD = data[:, list(dims)]

    # Calculate the minimum and maximum values for each dimension
    min_x1, min_x2 = np.min(twoD, axis=0) + 2 * grid_step
    max_x1, max_x2 = np.max(twoD, axis=0) - grid_step

    # Create a grid of points within the range of the data
    xs, ys = np.mgrid[min_x1:max_x1:grid_step, min_x2:max_x2:grid_step]
    grid_points = np.c_[xs.ravel(), ys.ravel()]

    # Fit the model to the data if it's not already fitted
    if not check_fitted(model):
        model = model.fit(twoD, tgt)

    # Make predictions at each point on the grid
    preds = model.predict(grid_points).reshape(xs.shape)

    # Create a pseudocolor plot of the predictions at the grid points
    ax.pcolormesh(xs, ys, preds, cmap=plt.cm.coolwarm)
    ax.set_xlim(min_x1, max_x1)
    ax.set_ylim(min_x2, max_x2)
    ax.set_title(f'Decision Boundary for {get_model_name(model)}')

    # Scatter plot the original data
    for i, (class_num, name) in enumerate(class_labels.items()):
        ax.scatter(*twoD[tgt==class_num].T, label=name, color=colors[i])

    # Set legend below the plot
    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=len(class_labels))

def check_fitted(model):
    """
    Returns True if model estimator has been fitted.

    Parameters:
    model: sklearn model object
    """
    return hasattr(model, "classes_")

# Define a function to get the name of a model
def get_model_name(model):
    """
    This function returns the name of a model's class as a string.

    Parameters:
    model: sklearn model object
    """
    return str(model.__class__).split('.')[-1][:-2]

The following code snippet loads the Iris dataset and splits it into training and testing sets for classification purposes. Additionally, it creates one-class variants for each species (Setosa, Versicolor, Virginica) within the dataset, transforming the problem into a one-versus-all binary classification. These steps enable both multi-class and binary classification, allowing for targeted analysis and prediction of individual Iris species.

# Import necessary libraries
from sklearn import datasets
from sklearn.model_selection import train_test_split

# Load the standard iris dataset
iris_dataset = datasets.load_iris()

# Split the full iris dataset into training and testing sets
X_train_iris_full, X_test_iris_full, y_train_iris_full, y_test_iris_full = train_test_split(
    iris_dataset.data, iris_dataset.target, test_size=0.33, random_state=21)

# Create one-class variants of the iris dataset
# Class 'setosa' vs others
X_train_setosa, X_test_setosa, y_train_setosa, y_test_setosa = train_test_split(
    iris_dataset.data, iris_dataset.target == 0, test_size=0.33, random_state=21)

# Class 'versicolor' vs others
X_train_versicolor, X_test_versicolor, y_train_versicolor, y_test_versicolor = train_test_split(
    iris_dataset.data, iris_dataset.target == 1, test_size=0.33, random_state=21)

# Class 'virginica' vs others
X_train_virginica, X_test_virginica, y_train_virginica, y_test_virginica = train_test_split(
    iris_dataset.data, iris_dataset.target == 2, test_size=0.33, random_state=21)

import matplotlib.pyplot as plt
import seaborn as sns

# Function to plot distribution
def plot_distribution(y, title):
    sns.countplot(x=y)
    plt.title(title)
    plt.show()

# Distributions of the target variables for each dataset
plot_distribution(y_train_iris_full, 'Full Iris Dataset')
plot_distribution(y_train_setosa, 'Setosa vs Others')
plot_distribution(y_train_versicolor, 'Versicolor vs Others')
plot_distribution(y_train_virginica, 'Virginica vs Others')

Applying the Concepts: Simplified Iris Problem

For demonstration purposes, we’ll use a simplified version of the iris problem. This version results in a single split, providing us with a clear understanding of how decision trees operate.

Setosa vs. Versicolor/Virginica

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# Convert the features to a dataframe
iris_setosa_data_df = pd.DataFrame(X_train_setosa[:, 2:], columns=iris_dataset.feature_names[2:])

# Add the target to the dataframe
iris_setosa_data_df['target'] = y_train_setosa

# Mapping the target values to their specific class names
target_names_dict = {True: 'setosa', False: 'versicolor & virginica'}

# Apply the mapping to the 'target' column
iris_setosa_data_df['target'] = iris_setosa_data_df['target'].map(target_names_dict)

palette_dict = {
    'setosa': (0.2980392156862745, 0.4470588235294118, 0.6901960784313725),  # Blue color
    'versicolor & virginica': (0.8666666666666667, 0.5176470588235295, 0.3215686274509804) # Orange color
}

# Plot using seaborn's pairplot function with the custom palette
sns.pairplot(iris_setosa_data_df, hue='target', palette=palette_dict)

# Add a title
plt.suptitle('Pairplot for the Setosa vs Versicolor & Virginica Classification Problem', y=1.02)

# Show the plot
plt.show()

The pairplot shows the relationship between petal length and petal width for Iris setosa compared to other species. It helps in visually identifying how well these features can distinguish Iris setosa from the other species.

# Import required libraries
import pydotplus
from sklearn import tree
from IPython.display import Image
import os

# Initialize the DecisionTreeClassifier
dtc_one_class = tree.DecisionTreeClassifier(random_state=42)

# Fit the model to the one-class variant of the iris dataset
# Use only the first two features
dtc_one_class.fit(X_train_setosa[:, 2:], y_train_setosa)

# Generate the DOT data for the tree
dot_data_setosa = tree.export_graphviz(dtc_one_class, out_file=None,
                                       feature_names=iris_dataset.feature_names[2:],
                                       class_names=["versicolor & virginica", "setosa"],
                                       filled=True, rounded=True)

# Check if the output directory exists
if not os.path.exists('outputs'):
    # If it doesn't exist, create it
    os.makedirs('outputs')

# Use pydotplus to generate a graph from the DOT data
graph_setosa = pydotplus.graph_from_dot_data(dot_data_setosa)

# Add a title to the graph
graph_setosa.set_label('Decision Tree for the Setosa vs Versicolor/Virginica Problem')
graph_setosa.set_fontsize('20')

# Define the output path
output_path_setosa = "outputs/iris_setosa_vs_others.png"

# Write the PNG to a file
graph_setosa.write_png(output_path_setosa)

# Create an Image object using the correct file path, and display it
Image(output_path_setosa)

import numpy as np
fig, ax = plt.subplots(figsize=(6,6))
plot_decision_boundary(ax=ax, data=X_train_setosa[:, :], tgt=y_train_setosa, model=dtc_one_class,
                       class_labels={0: 'Versicolor/Virginica', 1:'Setosa'}, dims=[2,3], colors=['red','blue'])

# Add x and y labels
ax.set_xlabel('Petal Length (cm)')
ax.set_ylabel('Petal Width (cm)')

Text(0, 0.5, 'Petal Width (cm)')

Versicolor vs. Setosa/Virginica

# Import necessary libraries
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# Convert the features to a dataframe
iris_versicolor_features_df = pd.DataFrame(X_train_versicolor[:, 2:], columns=iris_dataset.feature_names[2:])

# Add the target to the dataframe
iris_versicolor_features_df['target'] = y_train_versicolor

# Define a dictionary to map the target values to their respective class names
target_names_dict = {True: 'versicolor', False: 'setosa & virginica'}

# Apply the mapping to the 'target' column
iris_versicolor_features_df['target'] = iris_versicolor_features_df['target'].map(target_names_dict)

# Define the custom palette using the RGB values
palette_dict = {
    'versicolor': (0.2980392156862745, 0.4470588235294118, 0.6901960784313725),  # Blue color
    'setosa & virginica': (0.8666666666666667, 0.5176470588235295, 0.3215686274509804) # Orange color
}

# Use seaborn's pairplot function to visualize the pairwise relationships in the dataset with the custom palette
sns.pairplot(iris_versicolor_features_df, hue='target', palette=palette_dict)

# Add a title
plt.suptitle('Pairplot for the Versicolor vs Setosa/Virginica Classification Problem', y=1.02)

# Show the plot
plt.show()

# Initialize the DecisionTreeClassifier for versicolor
dtc_versicolor = tree.DecisionTreeClassifier(random_state=42)

# Fit the model to the one-class variant of the iris dataset
dtc_versicolor.fit(X_train_versicolor[:, 2:], y_train_versicolor)

# Generate the DOT data for the tree
dot_data_versicolor = tree.export_graphviz(dtc_versicolor, out_file=None,
                                           feature_names=iris_dataset.feature_names[2:],
                                           class_names=["setosa/virginica", "versicolor"],
                                           filled=True, rounded=True)

# Use pydotplus to generate a graph from the DOT data
graph_versicolor = pydotplus.graph_from_dot_data(dot_data_versicolor)

graph_versicolor.set_label('Decision Tree for the Versicolor vs Setosa/Virginica Classification Problem')
graph_versicolor.set_fontsize('20')

# Define the output path
output_path_versicolor = "outputs/iris_versicolor_vs_others.png"

# Write the PNG to a file
graph_versicolor.write_png(output_path_versicolor)

# Display the Image
Image(output_path_versicolor)

fig, ax = plt.subplots(figsize=(6,6))
plot_decision_boundary(ax=ax, data=X_train_versicolor[:, :], tgt=y_train_versicolor, model=dtc_versicolor,
                       class_labels={0: 'Setosa/Virginica', 1:'Versicolor'}, dims=[2,3], colors=['red','blue'])

# Add x and y labels
ax.set_xlabel('Petal Length (cm)')
ax.set_ylabel('Petal Width (cm)')

Text(0, 0.5, 'Petal Width (cm)')

Virginica vs. Versicolor/Setosa

# Import necessary libraries
import pandas as pd
import seaborn as sns

# Convert the features to a dataframe
iris_virginica_train_features_df = pd.DataFrame(X_train_virginica[:, 2:], columns=iris_dataset.feature_names[2:])

# Add the target to the dataframe
iris_virginica_train_features_df['target'] = y_train_virginica

# We'll use a dictionary to map the target values to their respective class names
target_names_dict = {True: 'virginica', False: 'setosa/versicolor'}

# Apply the mapping to the 'target' column
iris_virginica_train_features_df['target'] = iris_virginica_train_features_df['target'].map(target_names_dict)

# Use seaborn's pairplot function to visualize the pairwise relationships in the dataset
sns.pairplot(iris_virginica_train_features_df, hue='target')

plt.suptitle('Pairplot for the Virginica vs Setosa/Versicolor classification problem', y=1.02)

Text(0.5, 1.02, 'Pairplot for the Virginica vs Setosa/Versicolor classification problem')

# Initialize the DecisionTreeClassifier for virginica
dtc_virginica = tree.DecisionTreeClassifier(random_state=42)

# Fit the model to the one-class variant of the iris dataset
dtc_virginica.fit(X_train_virginica[:, 2:], y_train_virginica)

# Generate the DOT data for the tree
dot_data_virginica = tree.export_graphviz(dtc_virginica, out_file=None,
                                          feature_names=iris_dataset.feature_names[2:],
                                          class_names=["setosa/versicolor", "virginica"],
                                          filled=True, rounded=True)

# Use pydotplus to generate a graph from the DOT data
graph_virginica = pydotplus.graph_from_dot_data(dot_data_virginica)

graph_virginica.set_label('Decision Tree for the Virginica vs Versicolor/Setosa Classification Problem')
graph_virginica.set_fontsize('20')
# Define the output path
output_path_virginica = "outputs/iris_virginica_vs_others.png"

# Write the PNG to a file
graph_virginica.write_png(output_path_virginica)

# Display the Image
Image(output_path_virginica)

fig, ax = plt.subplots(figsize=(6,6))
plot_decision_boundary(ax=ax, data=X_train_virginica[:, :], tgt=y_train_virginica, model=dtc_virginica,
                       class_labels={0: 'Setosa/Versicolor', 1:'Virginica'}, dims=[2,3], colors=['red','blue'])

# Add x and y labels
ax.set_xlabel('Petal Length (cm)')
ax.set_ylabel('Petal Width (cm)')

Text(0, 0.5, 'Petal Width (cm)')

# Import necessary libraries
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# Convert the features to a dataframe
iris_full_features_df = pd.DataFrame(X_train_iris_full[:, 2:], columns=iris_dataset.feature_names[2:])

# Add the target to the dataframe
iris_full_features_df['target'] = y_train_iris_full

# We'll use a dictionary to map the target values to their respective class names
target_names_dict = {0: 'setosa', 1: 'versicolor', 2: 'virginica'}

# Apply the mapping to the 'target' column
iris_full_features_df['target'] = iris_full_features_df['target'].map(target_names_dict)

# Define a color palette as a dictionary with specific RGB values
palette_dict = {
    'setosa': (0.8666666666666667, 0.5176470588235295, 0.3215686274509804),  # Orange color
    'versicolor': (0.2980392156862745, 0.4470588235294118, 0.6901960784313725), # Blue color
    'virginica': (0.3333333333333333, 0.6588235294117647, 0.40784313725490196)
}

# Use seaborn's pairplot function to visualize the pairwise relationships in the dataset with the custom palette
sns.pairplot(iris_full_features_df, hue='target', palette=palette_dict)

# Add a title
plt.suptitle('Pairplot for the Multi-class Iris Classification Problem', y=1.02)

# Show the plot
plt.show()

# Import necessary libraries
from sklearn import tree
from IPython.display import Image
import pydotplus

# Initialize the DecisionTreeClassifier for multi-class classification
dtc_iris_full = tree.DecisionTreeClassifier(random_state=42)

# Fit the model to the full iris dataset
dtc_iris_full.fit(X_train_iris_full[:, 2:], y_train_iris_full)

# Generate the DOT data for the tree
dot_data_iris_full = tree.export_graphviz(dtc_iris_full, out_file=None,
                                          feature_names=iris_dataset.feature_names[2:],
                                          class_names=iris_dataset.target_names,
                                          filled=True, rounded=True)

# Use pydotplus to generate a graph from the DOT data
graph_iris_full = pydotplus.graph_from_dot_data(dot_data_iris_full)

# Define the output path
output_path_iris_full = "outputs/iris_full.png"

# Write the PNG to a file
graph_iris_full.write_png(output_path_iris_full)

# Display the Image
Image(output_path_iris_full)

fig, ax = plt.subplots(figsize=(6,6))
plot_decision_boundary(ax=ax, data=X_train_iris_full, tgt=y_train_iris_full, model=dtc_iris_full,
                       class_labels={0: 'Setosa', 1: 'Versicolor', 2:'Virginica'}, dims=[2,3], colors=['red', 'blue', 'grey'])

# Add x and y labels
ax.set_xlabel('Petal Length (cm)')
ax.set_ylabel('Petal Width (cm)')

Text(0, 0.5, 'Petal Width (cm)')

# Import necessary libraries
from sklearn.svm import SVC
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
import matplotlib.pyplot as plt
import numpy as np

# Load the iris dataset
iris = load_iris()

# Prepare class labels and colors
class_labels = {i: name for i, name in enumerate(iris.target_names)}
colors = ['red', 'green', 'blue']

# A dictionary to hold the classifiers
# Here, only a Decision Tree Classifier (DTC) with a maximum depth of 3 is used
tree_classifiers = {'DTC': DecisionTreeClassifier(max_depth=3)}

# Create a figure and a set of subplots
fig, ax = plt.subplots(1, 1, figsize=(6, 5))

# Loop through the items in the dictionary
for name, mod in tree_classifiers.items():
    # plot_decision_boundary is a function that you have defined elsewhere in your code
    plot_decision_boundary(ax, iris.data, iris.target, mod, [0, 1], class_labels, colors)

    # Set the title of the plot to be the name of the classifier
    ax.set_title("Decision Boundary for " + name)

# Create legend handlers manually
handlers = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=c, markersize=5) for c in colors]
# Put a legend below current axis

# Add a figure text
fig.text(0.1, 0.01, "Decision boundary plot of a decision tree classifier trained on the first two features of the iris dataset.", ha='left')

# Automatically adjust subplot parameters to give specified padding
plt.tight_layout(rect=[0, 0.08, 1, 0.96])  # leaving space at the bottom for the legend

# Display the plot
plt.show()

# Import required libraries
import pydotplus
import os
from sklearn import tree
from IPython.display import Image
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split

# Load the iris dataset
iris_dataset = load_iris()

# Split the data into train and test sets
iris_train_features, iris_test_features, iris_train_target, iris_test_target = train_test_split(
    iris_dataset.data, iris_dataset.target, test_size=0.3, random_state=42)

# Initialize the DecisionTreeClassifier
dtc = tree.DecisionTreeClassifier()

# Fit the model to the iris dataset
dtc.fit(iris_train_features, iris_train_target)

# Generate the DOT data
dot_data = tree.export_graphviz(dtc, out_file=None,
                                feature_names=iris_dataset.feature_names,
                                class_names=iris_dataset.target_names,
                                filled=True, rounded=True)

# Use pydotplus to generate a graph from the DOT data
graph = pydotplus.graph_from_dot_data(dot_data)

# Define the output path
output_path = "outputs/iris_full.png"

# Check if the output directory exists
if not os.path.exists('outputs'):
    # If it doesn't exist, create it
    os.makedirs('outputs')

# Write the PNG to a file
graph.write_png(output_path)

# Create an Image object using the correct file path, and display it
Image(output_path, width=600, height=600)

Understanding the Diabetes Dataset: Features and Target Variable

The diabetes dataset included in sklearn is a collection of various biometric and demographic features, utilized for regression analysis. Each feature in this dataset has been standardized, meaning that the mean has been subtracted and the result divided by the standard deviation for each column.

Standardization, also referred to as z-scoring, transforms the features so that each one has a mean of 0 and a standard deviation of 1. This process ensures that differences in the ranges of the features don’t disproportionately influence the model’s behavior. It neutralizes the effect of varying scales among features, like height measured in inches or income in dollars.

In the case of the diabetes dataset, even categorical values such as sex have been numerically recorded and then standardized. This process results in values that may seem unconventional, like negative ages, or gender coded as numerical values rather than “M” or “F.” The ultimate goal of standardization is to prepare the data for efficient use in machine learning models by creating uniformity in the feature scales.

Preprocessing the Diabetes Dataset

The diabetes dataset includes ten baseline variables such as age, sex, body mass index, average blood pressure, and six blood serum measurements. These features are used to model the progression of the disease after one year. To make the data ready for machine learning models, it’s often preprocessed with the following steps:

1. Standardization: Each feature variable is scaled to have a mean of zero and a standard deviation of one. The standardized value is calculated using the formula:

   $$ x_{\text{standardized}} = \frac{{x - \mu}}{{\sigma}} $$

• $ x_{\text{standardized}} $: This is the standardized value of the feature, representing the transformed value after applying the standardization process.
• $ x $: The original value of the feature that you want to standardize.
• $ \mu $: The mean of the feature across all instances in the dataset. It represents the average value of that particular feature.
• $ \sigma $: The standard deviation of the feature across all instances in the dataset. It represents how much the values of that particular feature deviate from the mean.

Understanding Standardization in Data Normalization

Normalization is a crucial step in preprocessing data, especially in machine learning. A common method of normalization is standardization. This process involves two main steps:

1. Centering: Adjusting the mean of the data to zero.
2. Scaling: Scaling the data to have a standard deviation of 1.

Let’s dive into an example using a one-dimensional dataset.

# Importing necessary libraries
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# Creating 1D data for standardization (evenly spaced values)
values = np.linspace(-5, 10, 20)
data_df = pd.DataFrame(values, columns=['original_values'])

# Applying standardization (centering around mean and scaling by standard deviation)
mean_original = data_df['original_values'].mean()
std_original = data_df['original_values'].std()
data_df['standardized_values'] = (data_df['original_values'] - mean_original) / std_original

# Plotting original and standardized data
fig, ax = plt.subplots(1, 1, figsize=(6, 3)) # Increased figure width
sns.stripplot(data=data_df, ax=ax)

# Rotating x-axis labels for better visibility
plt.xticks(rotation=45, ha='right')

# Optional: moving the x-axis labels slightly down for better alignment
plt.subplots_adjust(bottom=0.2)

plt.show()

display(data_df.describe().loc[['mean', 'std']])

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler

# Generate random diabetes-like data for demonstration
original = np.random.uniform(-5, 5, 100)

# Apply standardization: centering the mean to 0 and scaling the standard deviation to 1
scaler = StandardScaler()
scaled = scaler.fit_transform(original.reshape(-1,1))[:,0]

# Create bins for visualization
bins = np.floor(original).astype(np.uint8) + 5

# Create a DataFrame to store original and scaled data
df = pd.DataFrame({'original': original,
                   'scaled': scaled,
                   'bin': bins})

# Melt the DataFrame to long format for seaborn plotting
df_melted = pd.melt(df, id_vars='bin', var_name='scale', value_name='value')

# Define a color palette with distinct colors for each bin
palette = sns.color_palette("tab10", n_colors=len(df['bin'].unique()))

# Map the colors to the bins
df_melted['color'] = df_melted['bin'].map(dict(zip(df['bin'].unique(), palette)))

# Create subplots for histograms and swarmplot
fig, axes = plt.subplots(2, 1, figsize=(4, 6))

# Plot histograms for original and scaled data
sns.histplot(df['original'], bins=20, color='blue', alpha=0.5, ax=axes[0], kde=True)
axes[0].set_title('Original Data Distribution')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Frequency')

sns.histplot(df['scaled'], bins=20, color='green', alpha=0.5, ax=axes[0], kde=True)
axes[0].legend(['Original', 'Scaled'])

# Create the swarmplot with distinct colors for each bin
sns.swarmplot(x='scale', y='value', data=df_melted, hue='bin', palette=palette, ax=axes[1]).legend_.remove()
axes[1].set_title('Original and Scaled Data Swarmplot')

plt.tight_layout()
plt.show()

Standardization is a process that adjusts the scale of data without altering its shape. It centers the data by subtracting the mean, and then scales it by dividing by the standard deviation, transforming the data into units of standard deviation. Mathematically, the standardization formula is given by:

$$\text{scaled} = \frac{{\text{original} - \text{mean(original)}}}{\text{stddev(original)}}$$

This process leaves the overall shape of the data unchanged. If the original data had a uniform distribution, the standardized data will still be uniform, but with a mean of zero. The sides may appear more compressed or stretched, but the fundamental shape remains intact. Even if the original data had multiple peaks (like a two-humped camel), standardizing it would still preserve those peaks.

In practical terms, standardization is used with tools such as MinMaxScaler and StandardScaler in data preprocessing. This can be useful in machine learning applications where different features may be measured on vastly different scales, and it is necessary to bring them to a common scale for modeling.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler

# Create two normal distributions to form the two-humped camel shape
hump1 = np.random.normal(loc=-2, scale=1, size=500)
hump2 = np.random.normal(loc=2, scale=1, size=500)

# Concatenate the two distributions
two_humped_camel_data = np.concatenate((hump1, hump2)).reshape(-1, 1)

# Apply standardization using StandardScaler from scikit-learn
scaler = StandardScaler()
standardized_data = scaler.fit_transform(two_humped_camel_data)

# Plot the original and standardized data
fig, axes = plt.subplots(2, 1, figsize=(6, 10))

sns.histplot(two_humped_camel_data, bins=20, color='blue', alpha=0.5, ax=axes[0], kde=True, label='Original')
sns.histplot(standardized_data, bins=20, color='green', alpha=0.5, ax=axes[0], kde=True, label='Standardized')
axes[0].set_title('Original and Standardized Two-Humped Camel Data')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Frequency')
axes[0].legend()

sns.swarmplot(x=two_humped_camel_data[:, 0], ax=axes[1], size=3, color='blue', label='Original')
sns.swarmplot(x=standardized_data[:, 0], ax=axes[1], size=3, color='green', label='Standardized')
axes[1].set_title('Original and Standardized Data Swarmplot')
axes[1].set_xlabel('Value')
axes[1].legend()

plt.tight_layout()
plt.show()

2. Normalization: The dataset is further processed so that each instance (or feature vector) has a squared length of 1. This means that the sum of the squared values for each instance equals 1:

   $$ \sum x_i^2 = 1 $$

   • $ x_i $: The standardized value of the $i$-th feature of a particular instance.
   • The sum symbol $ \sum $: Represents the summation of the squared standardized values across all features in a particular instance.
   • The equation ensures that the sum of the squared standardized values for each instance equals 1, a process that scales the data to a specific range suitable for certain machine learning algorithms.

from sklearn.datasets import load_diabetes
import pandas as pd

# Load the diabetes dataset
diabetes_data = load_diabetes()

# Convert to DataFrame
diabetes_df = pd.DataFrame(diabetes_data['data'], columns=diabetes_data['feature_names'])

# Calculate the sum of squares for each feature across the rows
sum_of_squares = diabetes_df.apply(lambda column: (column**2).sum(), axis=0)

# Print the sum of squares for each feature
print(sum_of_squares)

age    1.0
sex    1.0
bmi    1.0
bp     1.0
s1     1.0
s2     1.0
s3     1.0
s4     1.0
s5     1.0
s6     1.0
dtype: float64

The standardization process scales the features so that they all have the same mean and standard deviation, making them directly comparable. The normalization to have a squared length of 1 further scales the data instances in a way that makes them suitable for models sensitive to the scale of the features.

These preprocessing steps are critical when applying machine learning models to the diabetes dataset. The transformed data can be effectively used for regression analysis, allowing us to predict the continuous outcome of disease progression based on the ten feature variables. Moreover, it provides a practical opportunity to delve into more advanced concepts like feature selection and optimization techniques in machine learning.

Feature Variables

The dataset contains ten baseline variables, age, sex, body mass index, average blood pressure, and six blood serum measurements:

1. Age: Age in years.
2. Sex: A binary variable indicating gender.
3. BMI (Body Mass Index): It’s a measure that quantifies an individual’s body weight in relation to their height. It’s an indicator of body fatness and often used to categorize individuals into weight classes that may lead to health problems.
4. Average Blood Pressure: This represents the patient’s average blood pressure.
5. S1: Total serum cholesterol.
6. S2: Low-density lipoproteins (LDL cholesterol).
7. S3: High-density lipoproteins (HDL cholesterol).
8. S4: Total cholesterol / HDL cholesterol ratio.
9. S5: Log of serum triglycerides level.
10. S6: Blood sugar level.

Target Variable

• Progression of the disease: The target variable in the diabetes dataset is a quantitative measure of the progression of diabetes one year after the baseline, based on the response of interest.

# Import necessary libraries
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsRegressor
import pandas as pd

# Load diabetes dataset
diabetes = load_diabetes()
diabetes_df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)

# Adding the target variable (disease progression) to the DataFrame
diabetes_df['Disease Progression'] = diabetes.target

# Renaming the columns with the provided names
column_names = {
    'age': 'Age',
    'sex': 'Sex',
    'bmi': 'BMI (Body Mass Index)',
    'bp': 'Average Blood Pressure',
    's1': 'Total Serum Cholesterol',
    's2': 'Low-Density Lipoproteins (LDL Cholesterol)',
    's3': 'High-Density Lipoproteins (HDL Cholesterol)',
    's4': 'Total Cholesterol / HDL Cholesterol Ratio',
    's5': 'Log of Serum Triglycerides Level',
    's6': 'Blood Sugar Level'
}

diabetes_df.rename(columns=column_names, inplace=True)

diabetes_df.head()

from sklearn.datasets import load_diabetes
import pandas as pd

# Load the diabetes dataset
diabetes_data = load_diabetes()

# Convert to DataFrame
diabetes_df = pd.DataFrame(diabetes_data['data'], columns=diabetes_data['feature_names'])

# Verify the mean is close to 0 for each feature column
mean_values = diabetes_df.mean()

# Set a threshold under which values will be considered zero
threshold = 1e-15

# Force values below the threshold to be exactly zero
mean_values_zeroed = mean_values.apply(lambda x: 0 if abs(x) < threshold else x)

print("Mean values (forced to zero):")
print(mean_values_zeroed)

from sklearn.preprocessing import StandardScaler

# Standardize the data
scaler = StandardScaler()
diabetes_df_standardized = pd.DataFrame(scaler.fit_transform(diabetes_df), columns=diabetes_df.columns)

# Calculate the standard deviation for the standardized data
standard_deviations_standardized = diabetes_df_standardized.std()
print("Standard deviations (standardized):")
print(standard_deviations_standardized)

# Calculate the sum of squares for each feature across the rows
sum_of_squares = diabetes_df.apply(lambda column: (column**2).sum(), axis=0)

# Print the sum of squares for each feature
print("Sum of squares:")
print(sum_of_squares)

Mean values (forced to zero):
age    0
sex    0
bmi    0
bp     0
s1     0
s2     0
s3     0
s4     0
s5     0
s6     0
dtype: int64
Standard deviations (standardized):
age    1.001133
sex    1.001133
bmi    1.001133
bp     1.001133
s1     1.001133
s2     1.001133
s3     1.001133
s4     1.001133
s5     1.001133
s6     1.001133
dtype: float64
Sum of squares:
age    1.0
sex    1.0
bmi    1.0
bp     1.0
s1     1.0
s2     1.0
s3     1.0
s4     1.0
s5     1.0
s6     1.0
dtype: float64

# Import necessary libraries
import seaborn as sns
import pandas as pd
from sklearn.datasets import load_diabetes

# Renaming the target column to represent Disease Progression
column_names['Disease Progression'] = 'Disease Progression (One Year After Baseline)'

# Apply the renaming
diabetes_df.rename(columns=column_names, inplace=True)

unique_values = diabetes_df['Sex'].unique()

palette = sns.color_palette('husl', len(unique_values))

custom_palette = dict(zip(unique_values, palette))

g = sns.pairplot(diabetes_df, hue='Sex', palette=custom_palette, height=2.5, aspect=1)

# Adjust x-axis labels
for ax in g.axes[-1, :]:
    label = ax.get_xlabel()
    ax.set_xlabel(label, rotation=45, ha='center')  # 'ha' sets the horizontal alignment
    ax.xaxis.labelpad = 10  # Adds padding to move the label away from the axis

# Adjust y-axis labels
for ax in g.axes[:, 0]:
    label = ax.get_ylabel()
    ax.set_ylabel(label, rotation=0, va='center')  # 'va' sets the vertical alignment
    ax.yaxis.labelpad = 125  # Adds padding to move the label away from the axis

# Show the plot
plt.show()

Linear Regression

Linear Regression (LR) is a fundamental statistical technique used to model the relationship between a dependent variable and one or more independent variables by fitting a linear equation to the observed data. Despite its simplicity, LR has a broad application across various fields, from mathematics to science.

Evaluating Model Performance with Scikit-Learn’s Metrics

When building machine learning models, the choice of evaluation metrics is vital. It guides the optimization and gives insights into how well the model is performing. Scikit-learn, a popular machine learning library, provides a comprehensive set of metrics for different machine learning tasks, such as classification, regression, and clustering.

Classification Metrics

For classification tasks, the following metrics are often employed:

• Accuracy: Measures the proportion of correctly classified instances.
• F1 Score: Harmonic mean of precision and recall, a balanced measure.
• Precision: Ratio of true positive predictions to the total predicted positives.
• Recall: Ratio of true positive predictions to the total actual positives.
• ROC AUC: Area under the Receiver Operating Characteristic Curve, a comprehensive measure of model’s ability to distinguish between classes.

Regression Metrics

In regression, the following metrics provide different perspectives on error and fit:

• R-squared: Explains how well the independent variables explain the variability of the dependent variable.
• Mean Squared Error (negated): Average squared difference between predicted and actual values.
• Mean Absolute Error (negated): Average absolute difference between predicted and actual values.

Note: The negation in some metrics aligns them with scikit-learn’s convention where higher scores are better.

Clustering Metrics

Clustering tasks may utilize metrics such as:

• Adjusted Rand Index: Measures the similarity of two assignments, ignoring permutations.
• Homogeneity Score: Checks whether the clustering algorithm assigns only those data points that are members of a single class to the same cluster.

Accessing Available Scorers

You can access a complete list of available scorers in scikit-learn using the following code snippet:

from sklearn import metrics

# List of available scorers
scorers = metrics.get_scorer_names()

# Print the available scorers
for scorer in scorers:
    print(scorer)

accuracy
adjusted_mutual_info_score
adjusted_rand_score
average_precision
balanced_accuracy
completeness_score
explained_variance
f1
f1_macro
f1_micro
f1_samples
f1_weighted
fowlkes_mallows_score
homogeneity_score
jaccard
jaccard_macro
jaccard_micro
jaccard_samples
jaccard_weighted
matthews_corrcoef
max_error
mutual_info_score
neg_brier_score
neg_log_loss
neg_mean_absolute_error
neg_mean_absolute_percentage_error
neg_mean_gamma_deviance
neg_mean_poisson_deviance
neg_mean_squared_error
neg_mean_squared_log_error
neg_median_absolute_error
neg_negative_likelihood_ratio
neg_root_mean_squared_error
normalized_mutual_info_score
positive_likelihood_ratio
precision
precision_macro
precision_micro
precision_samples
precision_weighted
r2
rand_score
recall
recall_macro
recall_micro
recall_samples
recall_weighted
roc_auc
roc_auc_ovo
roc_auc_ovo_weighted
roc_auc_ovr
roc_auc_ovr_weighted
top_k_accuracy
v_measure_score

Scikit-learn’s rich set of evaluation metrics offers flexibility in assessing models across various tasks. By understanding and wisely choosing the right metrics, you can align the model’s optimization with your specific objectives and gain deeper insights into its performance.

def print_default_scorer(model):
    """Prints the default scoring metric for a given scikit-learn model."""
    try:
        default_scorer_description = model.score.__doc__.splitlines()[0]
        print(f"Default scorer for {type(model).__name__}: {default_scorer_description}")
    except AttributeError:
        print(f"{type(model).__name__} does not have a 'score' method.")

# Example usage
from sklearn.linear_model import LinearRegression

linear_regression = LinearRegression()
print_default_scorer(linear_regression)

Default scorer for LinearRegression: Return the coefficient of determination of the prediction.

Understanding Regression Metrics: A Closer Look at R², MAE, MSE, and Median Absolute Error

When working with regression models, understanding the different evaluation metrics is essential for interpreting the performance of the model. Let’s delve into four key metrics and examine how they impact our model assessments.

1. R² (R-squared)

$R^2$ is the default metric for all regressors in many machine learning libraries, such as scikit-learn. Representing the coefficient of determination, R² measures how well the predictions of the model fit the actual observed values. It provides a score between 0 and 1, where higher values indicate a more accurate model. Let’s peel back the layers and understand what $ R^2 $ represents.

Two Models, One Metric

Imagine two models:

1. Our Model: This is the predictive model we’ve built. We’ll measure its performance by calculating the sum of squared differences between our predictions and the actual values:

   $$ \text{ErrorsModel} = \sum_i (\text{Predictions}_i - \text{TrueValues}_i)^2 $$

2. Baseline Model: A simple model that only predicts the mean of the target variable. The sum of squared differences for this model is:

3. $$ \text{ErrorsBaseline} = \sum_i (\text{MeanPredictions}_i - \text{TrueValues}_i)^2 $$

# Sample predictions from our model
predictions = np.array([1, 2, 3])
print("Predictions from our model:", predictions)

# Baseline predictions using the mean value
mean_predictions = np.array([2, 2, 2])
print("Baseline predictions (using mean):", mean_predictions)

# Actual observed values
true_values = np.array([2, 3, 4])
print("Actual observed values:", true_values)

# Calculating the errors for our model
errors_our_model_array = (predictions - true_values)**2
print("Squared errors for our model:", errors_our_model_array)

# Summing the squared errors for our model
errors_our_model = np.sum(errors_our_model_array)
print("Sum of squared errors for our model:", errors_our_model)

# Calculating the errors for the baseline model
errors_baseline_array = (mean_predictions - true_values)**2
print("Squared errors for the baseline model:", errors_baseline_array)

# Summing the squared errors for the baseline model
errors_baseline = np.sum(errors_baseline_array)
print("Sum of squared errors for the baseline model:", errors_baseline)

# You can further compare these two values to see how well your model is performing compared to a simple baseline.

Predictions from our model: [1 2 3]
Baseline predictions (using mean): [2 2 2]
Actual observed values: [2 3 4]
Squared errors for our model: [1 1 1]
Sum of squared errors for our model: 3
Squared errors for the baseline model: [0 1 4]
Sum of squared errors for the baseline model: 5

What Does $ R^2 $ Tell Us?

The $ R^2 $ value is computed using:

$ R^2 = 1 - \frac{\text{ErrorsOurModel}}{\text{ErrorsBaseline}} $

This metric illustrates how much better our predictive model is compared to a simplistic model that only predicts the mean. A higher $ R^2 $ value (closer to 1) suggests a more effective model in capturing the underlying data patterns.

In our journey through machine learning, the $ R^2 $ metric stands as a valuable tool in interpreting the performance of regression models. It puts our model’s errors into perspective by comparing them with a basic baseline, providing a normalized measure of how well our model fits the data. This understanding helps us make informed decisions in selecting and fine-tuning our models, striving for a balance between simplicity and predictive power.

Understanding $ R^2 $ - Coefficient of Determination

The $ R^2 $ value, also known as the coefficient of determination, is a popular metric used to evaluate the performance of a regression model. In simple terms, $ R^2 $ tells us how well our model is doing compared to a naive baseline model that always predicts the mean target value. Here’s a breakdown of how it’s calculated:

1. Compute the Sum of Squared Errors for Our Model: This is the sum of the squares of the differences between the predicted values and the actual values.
2. Compute the Sum of Squared Errors for the Baseline Model: This baseline model always predicts the mean target value, and we calculate its sum of squared errors.
3. Calculate the Ratio of Errors: Divide the sum of squared errors of our model by the sum of squared errors of the baseline model.
4. Compute $ R^2 $: The final $ R^2 $ value is given by $ 1 $ minus the above ratio.

The $ R^2 $ value ranges from 0 to 1. A value of 1 indicates perfect predictions, while a value of 0 suggests that the model is no better than the baseline model.

Let’s illustrate this with a code example:

import numpy as np

# Sample predictions and true values
predictions = np.array([1, 2, 3])
print("Predictions from our model:", predictions)

true_values = np.array([2, 3, 4])
print("Actual observed values:", true_values)

mean_predictions = np.array([2, 2, 2]) # Assuming mean as constant prediction
print("Baseline predictions (using mean):", mean_predictions)

# Errors for our model
errors_our_model_array = (predictions - true_values)**2
print("Squared errors for our model:", errors_our_model_array)

# Sum of squared errors for our model
errors_our_model = np.sum(errors_our_model_array)
print("Sum of squared errors for our model:", errors_our_model)

# Errors for the baseline (mean) model
errors_baseline_array = (mean_predictions - true_values)**2
print("Squared errors for the baseline model:", errors_baseline_array)

# Sum of squared errors for the baseline (mean) model
errors_baseline = np.sum(errors_baseline_array)
print("Sum of squared errors for the baseline model:", errors_baseline)

# Calculating R^2
r_2 = 1 - (errors_our_model / errors_baseline)
print("Calculated R^2 value: {:.2f}".format(r_2))

Predictions from our model: [1 2 3]
Actual observed values: [2 3 4]
Baseline predictions (using mean): [2 2 2]
Squared errors for our model: [1 1 1]
Sum of squared errors for our model: 3
Squared errors for the baseline model: [0 1 4]
Sum of squared errors for the baseline model: 5
Calculated R^2 value: 0.40

In this example, the calculated $ R^2 $ value is $ 0.40 $.

Delving into $ R^2 $: Beyond Linear Regression

The $ R^2 $ metric is generally well understood within the context of linear regression. However, it’s essential to explore how it behaves with models other than linear regression, especially when it leads to negative values.

Properties in Linear Regression

1. Range of Values: In linear regression, the $\frac{\text{SSE}_{\text{OurModel}}}{{\text{SSE}_{\text{mean}}}} $ term will always be between 0 and 1.

2. Interpretation: When the model is identical to the mean (constant term), the value is 1. When the model makes no errors, the value is 0. Linear regression, when fit on training data, can’t perform worse than the mean model.

When Using Non-linear Models

1. Negative $ R^2 $ Values: If a non-linear model performs worse than the mean, the $ R^2 $ value could become negative. This is counterintuitive, as squared numbers are typically positive.
2. Normalization by Baseline Model: The ratio between the SSEs provides a normalized performance comparison against a standard baseline model. This can be seen as a comparison between the mean squared errors (MSE) of the models.
3. Generalized Understanding: $ R^2 $ can be viewed as a measure of the MSE that is normalized by the MSE obtained from a simple baseline model predicting the mean. It provides a way to compare different models by analyzing the ratios of their MSEs or SSEs.

In essence, while the $ R^2 $ value has a clear interpretation within linear regression, it may be more complex with other types of models. The key insight is that it serves as a normalized measure of error, allowing comparison between models based on their performance relative to a baseline.

The $ R^2 $ value’s flexibility and ability to capture complex relationships make it a powerful metric for evaluating regression models, even when the underlying relationships are not linear. It’s a reminder that the understanding and application of statistical measures often require considering the specific context and nature of the data and models being used.

Manually Computing $ R^2 $ Values: A Comparison with Sklearn

In this section, we explore how to compute the $ R^2 $ value both manually and using the Scikit-learn library. By doing this, we gain a deeper understanding of the calculations involved and ensure that our manual computation aligns with the function provided by the library.

The Experiment

The experiment being conducted here is comparing the performance of a baseline model, which predicts the mean of the targets, with the performance of a more complex model trained on the training data.

1. Using Sklearn’s DummyRegressor

The DummyRegressor provides a simple baseline model that always predicts the mean. This can be used to calculate the $ R^2 $ value for our test set.

2. Manual Computations

We can also manually compute the sum of squared errors (SSE) using the mean values of the training and test targets, and compare these values with the SSE calculated by the sklearn model.

The Calculation

The following calculations are performed:

1. Calculate $ R^2 $ using Sklearn: This value is stored in base_r2_sklearn.
2. Compute $ R^2 $ manually using SSE: The manual calculation uses the Sum of Squared Errors (SSE) of the baseline predictions (sse_base_preds) and the SSE of the mean test (sse_mean_test). The ratio of these SSEs is subtracted from 1 to get the manual $ R^2 $.

Code Explanation

• sse_base_preds/sse_mean_test: The ratio of the SSE of the baseline predictions to the SSE of the mean of the test targets.
• 1-(sse_base_preds/sse_mean_test): The manual calculation of $ R^2 $, subtracting the ratio from 1.

The Counterintuitive Result

The result is that the manual calculation of $ R^2 $ is very close to the value computed by Sklearn, which might be surprising.

Why is this Happening?

Sklearn’s $ R^2 $ is calculating its base model using the mean of the true values we are testing against. This means that the comparison is not between my_model.fit(train) and mean_model.fit(train), but rather between my_model.fit(train) and mean_model.fit(test). This might be considered counterintuitive since you might expect that the mean model would be fit on the training data.

from sklearn.dummy import DummyRegressor
from sklearn.metrics import r2_score
from sklearn.datasets import load_diabetes
import numpy as np

diabetes = load_diabetes()

# Split into training and testing sets
X_train_diabetes, X_test_diabetes, y_train_diabetes, y_test_diabetes = train_test_split(
    diabetes.data, diabetes.target, test_size=0.2, random_state=42
)

# Fit the DummyRegressor to the training data
baseline_model = DummyRegressor(strategy='mean')
baseline_model.fit(X_train_diabetes, y_train_diabetes)
baseline_predictions = baseline_model.predict(X_test_diabetes)

# Calculate R2 using Sklearn
base_r2_sklearn = r2_score(y_test_diabetes, baseline_predictions)

# Manual computations
base_errors = baseline_predictions - y_test_diabetes
sse_baseline_predictions = np.dot(base_errors, base_errors)

train_mean_errors = np.mean(y_train_diabetes) - y_test_diabetes
sse_mean_train = np.dot(train_mean_errors, train_mean_errors)

test_mean_errors = np.mean(y_test_diabetes) - y_test_diabetes
sse_mean_test = np.dot(test_mean_errors, test_mean_errors)

print("Sklearn train-mean model SSE (on test):", sse_baseline_predictions)
print("Manual train-mean model SSE (on test):", sse_mean_train)
print("Manual test-mean model SSE (on test):", sse_mean_test)

ratio_sse = sse_baseline_predictions / sse_mean_test
manual_r2 = 1 - ratio_sse

print("Ratio of SSE (baseline predictions to mean test):", ratio_sse)
print("Manual R2:", manual_r2)
print("R2 score using Sklearn:", base_r2_sklearn)

Sklearn train-mean model SSE (on test): 477176.4776942276
Manual train-mean model SSE (on test): 477176.4776942276
Manual test-mean model SSE (on test): 471535.5056179775
Ratio of SSE (baseline predictions to mean test): 1.0119629847785423
Manual R2: -0.011962984778542296
R2 score using Sklearn: -0.011962984778542296

This exercise demonstrates how the manual computation of the $ R^2 $ value can give us insight into the behavior of the metric, allowing us to understand what’s happening behind the scenes when we call a library function. It also illustrates a nuanced point about how $ R^2 $ is computed in Sklearn, specifically how the baseline model is determined using the mean of the test targets, not the training targets. This detail leads to the counterintuitive result observed.

from sklearn.dummy import DummyRegressor
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
import pandas as pd
from collections import defaultdict

# WARNING! Don't try this at home, boys and girls!
# We are fitting on the *test* set... to mimic the behavior of sklearn R^2.

# Fit DummyRegressor on the test set
test_base = DummyRegressor(strategy='mean')
test_base.fit(X_test_diabetes, y_test_diabetes)
test_base_preds = test_base.predict(X_test_diabetes)
test_base_mse = mean_squared_error(test_base_preds, y_test_diabetes)

# Models to evaluate
models = [KNeighborsRegressor(n_neighbors=3), LinearRegression()]

results = defaultdict(dict)

# Train and test each model
for m in models:
    preds = (m.fit(X_train_diabetes, y_train_diabetes).predict(X_test_diabetes))
    mse = mean_squared_error(preds, y_test_diabetes)
    r2 = r2_score(y_test_diabetes, preds)
    model_name = type(m).__name__
    results[model_name]['R^2'] = r2
    results[model_name]['MSE'] = mse

# Create a DataFrame with the results
df = pd.DataFrame(results).T
df['Norm_MSE'] = df['MSE'] / test_base_mse
df['1-R^2'] = 1 - df['R^2']

display(df)

1. Fitting on Test Set: A DummyRegressor is fit on the test set using the mean strategy. This is done to mimic the behavior of Sklearn’s $ R^2 $ computation, and it’s not a common practice (hence the warning comment).
2. Model Evaluation: Two models (KNeighborsRegressor and LinearRegression) are trained on the training set and evaluated on the test set. For each model, the Mean Squared Error (MSE) and $ R^2 $ are computed.
3. Normalized MSE: The MSE for each model is normalized by the MSE of the DummyRegressor that was fit on the test set.
4. Results Display: A Pandas DataFrame is used to display the results, including the normalized MSE and $ 1 - R^2 $.

The code illustrates that the $ R^2 $ computed by Sklearn is equivalent to the MSE of the model normalized by the fit-to-the-test-sample mean model. If we knew the mean of our test targets, the value tells us how well we would perform in comparison to predicting that known mean.

This exercise showcases an uncommon and counterintuitive practice of fitting a model on the test set to mimic the behavior of Sklearn’s $ R^2 $. It is valuable as an illustration but is not a typical approach in practice.

Recommendations on $ R^2 $

The coefficient of determination, $ R^2 $, is a commonly used metric in regression analysis, but it comes with some nuances and complexities that may make it less desirable for certain applications. Below, we summarize the key points:

1. Scientific and Statistical Baggage:

   • $ R^2 $ has a lot of scientific and statistical interpretations that may not be directly applicable.
   • Terms like “percent”, “linear”, or “explained” should be viewed with skepticism when applied to Sklearn’s implementation of $ R^2 $.
   • Some interpretations are true under specific circumstances but not universally.

2. Multiple Formulas for $ R^2 $:

   • Sklearn uses one particular formula, but there are multiple formulas for $ R^2 $.
   • These formulas are equivalent for linear models with an intercept, but not always in other scenarios.
   • This inconsistency can lead to confusion.

3. Relationship to a Normalized MSE:

   • $ R^2 $ has a simple relationship to a very unusual thing: a normalized MSE computed on a test-sample-trained mean model.
   • This relationship can create misunderstandings about what $ R^2 $ actually represents.

Considering the complexities associated with $ R^2 $, the recommendation is to use Mean Squared Error (MSE) or Root Mean Squared Error (RMSE) instead. If normalization is desired, you can compare the regression model with a training-set-trained mean model. This approach avoids the ambiguity and potential misconceptions associated with $ R^2 $.

2. Mean Absolute Error (MAE)

MAE calculates the average absolute difference between predicted and actual values. It treats errors equally, so an error of 10 has a penalty of 10. This metric is ideal when you want to give equal weight to all errors and are less concerned with outliers.

3. Mean Squared Error (MSE)

MSE takes the differences between predicted and actual values, squares them, and then averages the results. This leads to larger errors being penalized more severely—an error of 10 results in a penalty of 100. If your analysis requires putting significant emphasis on larger errors, MSE may be the right choice.

4. Median Absolute Error

The Median Absolute Error offers a robust approach, focusing on the median of the errors rather than the mean. It protects against extreme outliers, making it suitable for scenarios where you are willing to tolerate some large errors, provided the overall trend is consistent.

Choosing the appropriate regression metric depends on the specific requirements of your project. While MAE provides a balanced view, treating all errors equally, MSE is more sensitive to larger discrepancies. Median Absolute Error offers protection against extreme outliers.

These metrics offer different perspectives on the model’s performance, and understanding their nuances helps in selecting the right one for your needs. Whether you are interested in penalizing larger errors more harshly or need a more balanced view, these metrics provide valuable insights into the quality and reliability of your regression models.

Residual Plots

Residual plots are graphical techniques for evaluating regressors, serving as a regression analog of confusion matrices. Here, we will look into two types of plots:

1. Error Plots: This consists of plotting actual (true) target values against predicted values. A perfect prediction would result in points along the line (y = x). Any deviation from this line represents an error.

2. Regression Errors Plot: A plot showing the difference between predicted and actual values, highlighting where the model is predicting too high (positive error) or too low (negative error). Connecting bars from data points to the perfect line are used to visualize the errors.

Let’s proceed to create these plots using a simple example:

import pandas as pd
import matplotlib.pyplot as plt
import itertools as it

# Sample data
ape_df = pd.DataFrame({'predicted': [4, 2, 9],
                       'actual': [3, 5, 7]})
ape_df['error'] = ape_df['predicted'] - ape_df['actual']
ape_df.index.name = 'example'
display(ape_df)

# Function to plot regression errors
def regression_errors(figsize, predicted, actual, errors='all'):
    fig, axes = plt.subplots(1, 2, figsize=figsize, sharex=True, sharey=True)
    df = pd.DataFrame({'actual': actual, 'predicted': predicted})

    for ax, (x, y) in zip(axes, it.permutations(['actual', 'predicted'])):
        ax.plot(df[x], df[y], '.', label='data') # Plotting data points
        ax.plot(df['actual'], df['actual'], '-', label='perfection') # Line of perfection
        ax.legend()
        ax.set_xlabel('{} Value'.format(x.capitalize()))
        ax.set_ylabel('{} Value'.format(y.capitalize()))
        ax.set_aspect('equal')

    axes[1].yaxis.tick_right()
    axes[1].yaxis.set_label_position("right")

    # Show connecting bars from data to perfect line
    if errors == 'all':
        errors = range(len(df))
    if errors:
        acts = df.actual.iloc[errors]
        preds = df.predicted.iloc[errors]
        axes[0].vlines(acts, preds, acts, 'r') # Vertical lines for errors
        axes[1].hlines(acts, preds, acts, 'r') # Horizontal lines for errors

    plt.show()

# Usage
regression_errors((6, 3), ape_df.predicted, ape_df.actual)

Interpretation of Residual Plots

In the context of regression analysis, residual plots help provide insights into the performance of a model by visualizing the discrepancies between predicted values and actual values.

Left-hand Graph: In this plot, with actual values on the x-axis and predicted values on the y-axis, the differences between prediction and reality are represented vertically. The orange line (y = x) symbolizes the ideal case where the predicted value equals the true actual value, meaning an error of zero. Any deviation from this line is an error.

Right-hand Graph: By swapping the axes, the differences between prediction and reality become horizontal. This graph answers the question: “For a given prediction, how do we perform?”

Example Usage with Diabetes Dataset & Linear Regression

Let’s apply this visualization method to a diabetes dataset using Linear Regression:

from sklearn.linear_model import LinearRegression

lr = LinearRegression()
preds = (lr.fit(X_train_diabetes, y_train_diabetes).predict(X_test_diabetes))
regression_errors((8, 4), preds, y_test_diabetes, errors=[-20])

Insights:

• When the actual value ranges from 200 to 250, the model seems to consistently predict low.
• When the predicted value is near 200, real values vary widely, ranging from 50 to 300.

This observation illustrates the ability to identify patterns and biases in the predictions, allowing for improvements in model tuning.

Introduction to Residual Plots

Residual plots are used to visualize the residuals, which are the differences between the actual values and the predicted values. They are different from error plots, where the error is calculated as the difference between the predicted values and the actual values. In residual plots, the calculation is reversed: $ \text{residuals} = \text{actual} - \text{predicted} $

Understanding Residuals in Predictive Modeling

In predictive modeling, two terms often arise: error and residual. Both are related but have distinct interpretations:

• Error: The difference between the predicted value and the actual value. An error of 2 means our prediction was over by 2. It reflects what happened in our prediction, whether we overestimated or underestimated the value.

• Residual: The adjustment needed to correct our prediction to match the actual value. A residual of 2 means we need to subtract 2 to arrive at the correct value. It signifies what adjustment is needed to align with reality.

This can be made concrete by considering the following example:

import pandas as pd

# Create a DataFrame with actual and predicted values
prediction_data = pd.DataFrame({
    'Actual Values': [3, 5, 7],
    'Predicted Values': [4, 2, 9]
})

# Calculate error as (Predicted - Actual)
prediction_data['Error (Predicted - Actual)'] = prediction_data['Predicted Values'] - prediction_data['Actual Values']

# Calculate residuals as (Actual - Predicted)
prediction_data['Residuals (Actual - Predicted)'] = prediction_data['Actual Values'] - prediction_data['Predicted Values']

# Set the index name
prediction_data.index.name = 'Example'

# Display the DataFrame
display(prediction_data)

Understanding the difference between residuals and errors is key when interpreting residual plots. While errors highlight the discrepancies between predictions and reality, residuals provide the mirror image, showing how reality deviates from the model’s predictions. Both perspectives can offer unique insights into the model’s behavior.

Residual Plots

Residual plots graph the predicted value against the residual for that prediction. These plots are helpful for understanding what adjustments are needed to improve predictions. In essence, they highlight what’s left over (the residue) after a prediction.

The two primary graphs in this context are:

1. Actual against Predicted: This shows the relationship between the actual values and the predicted values by the model.
2. Predicted against Residual: This plots the predicted values against the residuals, allowing you to see where the model is overestimating or underestimating.

By comparing these graphs, you can identify systematic patterns in your residuals, which can lead to insights into potential biases or inconsistencies in your model. Such analysis plays a vital role in refining the model for better performance.

def regression_residuals(ax, predicted, actual, show_errors=None, right=False):
    '''Plot regression residuals.

    Parameters:
    ax : matplotlib.axes.AxesSubplot
        The axes to plot on.
    predicted : array-like
        The predicted values.
    actual : array-like
        The actual values.
    show_errors : str or sequence of indices, optional
        If 'all', it shows vertical lines for all residuals. If a sequence, only for those indices.
    right : bool, optional
        If True, positions the y-axis tick marks to the right.
    '''

    # Create DataFrame to hold actual and predicted values
    df = pd.DataFrame({'actual': actual, 'predicted': predicted})
    # Calculate the error for each prediction
    df['error'] = df.actual - df.predicted

    # Plot predicted vs. error as dots
    ax.plot(df.predicted, df.error, '.')
    # Plot zero line for reference
    ax.plot(df.predicted, np.zeros_like(predicted), '-')

    # Optional setting for y-axis on the right
    if right:
        ax.yaxis.tick_right()
        ax.yaxis.set_label_position("right")

    ax.set_xlabel('Predicted Value')
    ax.set_ylabel('Residual')

    # If specified, show vertical lines to visualize the errors
    if show_errors == 'all':
        show_errors = range(len(df))
    if show_errors:
        preds = df.predicted.iloc[show_errors]
        errors = df.error.iloc[show_errors]
        ax.vlines(preds, 0, errors, 'r')

    # Debugging print statements
    # print("DataFrame with actual, predicted, and error values:")
    # print(df.head())

# Example usage
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
ax1.plot(ape_df.predicted, ape_df.actual, 'r.', # pred vs actual
         [0, 10], [0, 10], 'b-') # perfect line
ax1.set_xlabel('Predicted')
ax1.set_ylabel('Actual')
regression_residuals(ax2, ape_df.predicted, ape_df.actual, 'all', right=True)

The following code segment demonstrates how to compare the residual plots of two different regression models: a Linear Regression model and a k-Nearest Neighbors (k-NN) Regressor. The comparison is done using the same training and testing data for both models. Residual plots are a powerful way to visualize the discrepancies between predicted and actual values, helping to understand how the models might be over- or under-predicting.

from sklearn import linear_model, neighbors
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# Load the diabetes dataset
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target

# Split the data into training and testing sets
X_train_diabetes, X_test_diabetes, y_train_diabetes, y_test_diabetes = train_test_split(X, y, test_size=0.2, random_state=42)

# Define models: Linear Regression and k-NN Regressor
models = [
    (linear_model.LinearRegression(), 'Linear Regression Residuals'),
    (neighbors.KNeighborsRegressor(), 'k-NN Regressor Residuals')
]

# Create subplots to plot the residual plots side by side
fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharex=True, sharey=True)
fig.tight_layout()

# Loop through the models and plot the residual plots
for (model, title), ax in zip(models, axes):
    # Fit the model and predict the test data
    predicted_values = model.fit(X_train_diabetes, y_train_diabetes).predict(X_test_diabetes)

    # Call the regression_residuals function (defined earlier) to plot
    regression_residuals(ax, predicted_values, y_test_diabetes, show_errors=range(1, 21), right=True)

    # Set the title for each subplot
    ax.set_title(title)

plt.show()

Here, both models are trained on a diabetes dataset and then used to predict the target values for the test data. Residual plots are drawn for each model, with the residuals referred to as “predictive residuals.” This indicates that they are computed from the test set, aligning with the common practice of evaluating models based on unseen data. The section underscores the significance of residual plots not just for conventional statistical analyses but also for practical insights into model performance, failures, and potential improvements. Through standardization, we make our data more compatible for various machine learning algorithms, enabling better learning and prediction.

Understanding Linear Combinations

Linear combinations are a mathematical way of describing how different elements combine together, each multiplied by a certain amount, and then summed. This concept is present in many daily situations.

Recipe Example

Imagine you are baking a cake, and the recipe calls for specific quantities of different ingredients. You can think of the recipe as a linear combination of these ingredients:

• Flour: 2 cups at $0.50 per cup
• Sugar: 1 cup at $0.75 per cup
• Eggs: 3 eggs at $0.25 per egg
• Butter: 1 stick at $1.00 per stick

The total cost of the ingredients can be calculated as a linear combination:

$ (2 \times 0.50) + (1 \times 0.75) + (3 \times 0.25) + (1 \times 1.00) = 3.50 $

In this example, the quantities of the ingredients are like the “weights,” and the cost per unit is what gets multiplied by each weight. The total cost is the sum of these products.

quantity = np.array([2, 1, 3, 1])
costs = np.array([.50, .75, .25, 1.00])

np.sum(quantity * costs)

3.5

You can use np.dot to calculate the same result.

print(quantity.dot(costs))
print(np.dot(quantity, costs))
print(quantity @ costs)

3.5
3.5
3.5

Grocery Receipt Example

Consider a simple grocery receipt as an alternative example containing three items: wine, oranges, and muffins. The total cost can be calculated as a linear combination of the quantities and unit costs:

• Wine: 2 bottles at $12.50 per bottle
• Oranges: 12 oranges at $0.50 per orange
• Muffins: 3 muffins at $1.75 per muffin

The total cost is calculated as follows:

$ (2 \times 12.50) + (12 \times 0.50) + (3 \times 1.75) = 36.25 $

In this scenario, the quantities of each item act like the “weights,” and the unit costs are what gets multiplied by each weight. The final cost is the sum of these individual products, demonstrating a real-world example of a linear combination.

quantity = np.array([2, 12, 3])
costs = np.array([12.50, 0.50, 1.75])

np.sum(quantity * costs)

36.25

print(quantity.dot(costs), np.dot(quantity, costs), quantity @ costs)

36.25 36.25 36.25

Dot Product in Linear Combination

The concept of a linear combination can be expressed using a dot product, a fundamental operation in linear algebra. It consists of two main parts:

1. Pairwise Multiplication: Each corresponding pair of elements from two sequences is multiplied together. In a shopping context, this can be thought of as multiplying the quantity of each item by its cost to get a series of subtotals.
2. Summation: The subtotals obtained from the pairwise multiplication are added up to get the final result.

In mathematical notation, if you have quantities $ q $ and costs $ c $, the dot product is expressed as: $ \sum_{i} q_i c_i $

• The symbol $ \Sigma $ (a capital Greek sigma) indicates the summation.
• The expression $ q_i c_i $ represents the multiplication of corresponding elements.
• The index $ i $ ties the pieces together, iterating through each element in sequence.

In simpler terms, the dot product is the “sum product” of the quantities and costs. In programming, libraries like NumPy provide functions to easily calculate the dot product, allowing for efficient computation of linear combinations.

for q_i, c_i in zip(quantity, costs):
  print(f"{q_i} {c_i:5.2f} --> {q_i * c_i:5.2f}")

print("Total:", sum(q*c for q,c in zip(quantity,costs)))

2 12.50 --> 25.00
12  0.50 -->  6.00
3  1.75 -->  5.25
Total: 36.25

Normally we use NumPy, via np.dot, to do the work for us.

Understanding Simple Averages and Weighted Averages

When we talk about a simple average (often called the mean), we’re referring to the central value derived from summing up all the numbers in a set and then dividing by the count of those numbers. For instance, to find the average of the numbers $10$, $20$, and $30$, we’d sum them up and divide by $3$, which is the count.

Mathematically, for values $10$, $20$, and $30$: $ \text{Simple Average} = \frac{10 + 20 + 30}{3} $

This method inherently assumes that each value contributes equally, or has equal “weight”, to the average.

In contrast, a weighted average provides the flexibility to assign different importance or “weights” to each value. Instead of treating each value equally, we multiply each value by its respective weight, sum up the products, and then divide by the sum of the weights.

import numpy as np

# Define the values and calculate the number of elements
values = np.array([10.0, 20.0, 30.0])
num_elements = len(values)

# Create a weights array filled with equal weights (1/number of elements)
weights = np.full_like(values, 1/num_elements)

# Display values and weights
print("Values:", values)
print("Weights:", weights)

# Compute the mean of the values using the numpy mean function
print("Mean via numpy mean function:", np.mean(values))

# Compute the weighted average using the dot product of weights and values
print("Weighted average via dot product:", np.dot(weights, values))

Values: [10. 20. 30.]
Weights: [0.33333333 0.33333333 0.33333333]
Mean via numpy mean function: 20.0
Weighted average via dot product: 20.0

Weighted Sums and Weighted Averages

The mean of a set of numbers can be interpreted as a special case of a weighted sum—where each value is multiplied by an equal weight. When we assign different weights to these values, we delve into the realm of weighted averages.

In a basic average, each value has an equal contribution, which can be thought of as each value having an equal “weight”. However, in real-world scenarios, we often encounter situations where certain values have more significance or importance than others. This is where weighted averages come into play.

For instance, consider three options with values: 10, 20, and 30. If we decide to assign weights of 0.5, 0.25, and 0.25 to these options respectively, it indicates:

• The first option holds twice the significance as the other two.
• The latter two options are of equal significance.

Such weightage can arise from various factors:

1. Value Importance: The first option might have twice the impact or value as the other two in a given context.
2. Probability: The weights could signify that the first option is twice as likely to occur in a random event compared to the other two.

Interestingly, these interpretations—value importance and probability—are closely related. When applied to costs or quantities, the distinctions often blur, and we can view them as two sides of the same coin.

import numpy as np

# Define the values
values = np.array([10, 20, 30])

# Assign weights to each value
weights = np.array([.5, .25, .25])

# Compute the weighted average using dot product
weighted_average = np.dot(weights, values)

# Display the result
print(f"Given values: {values}")
print(f"Assigned weights: {weights}")
print(f"Weighted average: {weighted_average}")

Given values: [10 20 30]
Assigned weights: [0.5  0.25 0.25]
Weighted average: 17.5

Expected Value in a Game of Chance

A powerful application of the weighted average concept is in the realm of probability, particularly when calculating the expected value. In probabilistic terms, the expected value gives us the long-run average value of repetitions of the same experiment or scenario.

Imagine a simple dice game:

• You roll a standard six-sided die.
• If it lands on an odd number, you earn $1.00.
• If it lands on an even number, you lose $0.50.

The expected value is computed by multiplying each outcome (earn or lose) by its probability, then summing these products. Mathematically, the expected value $ E(game) $ is represented as: $ E(game) = \sum_{i} p_i v_i $ Where:

• $ p_i $ are the probabilities of the respective outcomes.
• $ v_i $ are the corresponding payouts (or values) of those outcomes.

If you were to play this dice game 100 times, the amount you’d expect to earn on average would be the expected value multiplied by 100. Over a small number of games, your earnings might vary, but as you play more games, your average earnings will gravitate towards the expected value.

import numpy as np

# Define the payoffs and probabilities
payoffs = np.array([1.0, -0.5])  # [odd (win), even (loss)]
probs = np.array([0.5, 0.5])     # Equal probability for both outcomes

# Compute expected value
expected_value = np.dot(payoffs, probs)

print(f"Expected value per game: ${expected_value:.2f}")

# Simulate the game for 10,000 rounds
def is_even(num):
    """Check if a number is even."""
    return int(num % 2 == 0)  # Convert to 1 for even and 0 for odd

total_winnings = 0.0
num_games = 10000

for _ in range(num_games):
    die_roll = np.random.randint(1, 7)
    total_winnings += payoffs[is_even(die_roll)]  # Add if even, subtract if odd

print(f"Total winnings after {num_games} games: ${total_winnings:.2f}")

Expected value per game: $0.25
Total winnings after 10000 games: $2491.00

Sum of Squares

The “sum of squares” refers to the total sum obtained by squaring each element of a set of numbers individually and then adding those squared values together. Mathematically, if you have a set of numbers $ a_1, a_2, a_3, \dots, a_n $, the sum of squares is:

$ \sum_{i=1}^{n} a_i^2 $

It’s called “sum of squares” because each element is multiplied by itself, resulting in its square.

To compute the sum of squares in Python using NumPy, we can simply square the array element-wise and then sum the squared values. This is efficient as it utilizes vectorized operations. Here’s an example:

import numpy as np

# Define an array of values
values = np.array([5, -3, 2, 1])

# Compute the squares of each value
squared_values = values ** 2

# Sum the squared values to obtain the sum of squares
sum_of_squares = np.sum(squared_values)

# Alternatively, the sum of squares can be computed using the dot product
dot_product_result = np.dot(values, values)

print("Squared values:", squared_values)
print("Sum of squares:", sum_of_squares)
print("Using dot product:", dot_product_result)

Squared values: [25  9  4  1]
Sum of squares: 39
Using dot product: 39

In the above code, the “sum of squares” is computed in two ways:

1. By squaring each element and then summing the result.
2. Using the dot product of the array with itself.

Both approaches yield the same result, as the dot product of a vector with itself essentially computes the sum of the squares of its elements.

Sum of Squared Errors (SSE)

When we’re making predictions using models or algorithms, there will inevitably be some disparity between our predicted values and the actual values. This difference is called an “error”. Errors can be positive (overestimation) or negative (underestimation). To quantify these errors without having them cancel each other out, we square each error. The aggregate of these squared values is known as the “sum of squared errors”.

Given actual values $a$ and predicted values $p$, the error $e$ for an individual data point is: $ e = p - a $

The squared error $e^2$ is: $ e^2 = (p - a)^2 $

The sum of squared errors (SSE) for all data points is: $ SSE = \sum (p_i - a_i)^2 $ where the summation is over all data points.

import numpy as np
import pandas as pd

# Given errors
errors = np.array([5, -5, 3.2, -1.1])

# Squared errors
squared_errors = errors ** 2

# Displaying errors and squared errors
df = pd.DataFrame({'Errors': errors, 'Squared Errors': squared_errors})

# Adding a total row for squared errors
total_squared_errors = squared_errors.sum()
total_row = pd.DataFrame({'Errors': ['Total'], 'Squared Errors': [total_squared_errors]})
df = pd.concat([df, total_row], ignore_index=True)

display(df)

# Sum of squared errors using dot product
sse = np.dot(errors, errors)
print("\nSum of Squared Errors using np.dot:", sse)

Sum of Squared Errors using np.dot: 61.45

Relation to Distances

You might find this reminiscent of the Pythagorean theorem: $c^2 = a^2 + b^2$. Here, $a$ and $b$ are like our errors - they’re the “distances” between actual and predicted values. In essence, when we talk about errors in predictive modeling, we’re referring to the “distance” or “discrepancy” between our predictions and reality.

In sum, understanding the concept of sum of squared errors is pivotal, as it’s a foundational stone in numerous machine learning algorithms. It offers a systematic way to measure and optimize the performance of models.

Finding the ‘Best’ Line: The Essence of Linear Regression

The task of linear regression is analogous to drawing a straight line through a set of dots on a piece of paper. While this concept might seem simple, multiple lines can be drawn through more than two dots, leading to the question: which line is the “best”?

The term “best” here needs definition. In the context of linear regression, the best line is the one that minimizes the vertical distance between the dots (data points) and the line itself. This concept can be visualized by imagining the line as an Olympic high-jump bar that can be raised or lowered while staying parallel to the ground.

Finding the best line becomes a matter of seeking a “happy medium” where the line is neither too high nor too low, balancing the vertical distance to all the points. Mathematically, this can be represented as lines of the form $ y = c $, where $ c $ is a constant that defines the height of our “bar.”

This intuitive understanding leads to a systematic and calculable way to compare different lines, determining which one best fits the data. The idea of finding the best-fitting line is foundational to linear regression and forms the basis for understanding more complex models in machine learning.

import numpy as np
from pprint import pprint
import pandas as pd
import matplotlib.pyplot as plt

def axis_helper(ax, lims):
    ax.set_xlim(lims)
    ax.set_xticks([])
    ax.set_ylim(lims)
    ax.set_yticks([])
    ax.set_aspect('equal')

def visualize_errors():
    # Sample data: two (x, y) points
    D = np.array([[3, 5], [4, 2]])
    x, y = D[:, 0], D[:, 1]

    #print("\nSample Data (x, y):")
    #pprint(D) # Using pprint to print aligned data

    # Define horizontal lines to represent predictions
    horizontal_lines = np.array([1, 2, 3, 3.5, 4, 5])
    results = []
    fig, axes = plt.subplots(1, 6, figsize=(10, 5))

    for h_line, ax in zip(horizontal_lines, axes.flat):
        #print("\n" + "-"*40) # Divider for easier viewing
        #print(f"Predicting with horizontal line (y = {h_line})") # Debugging print statement

        # Apply styling
        axis_helper(ax, (0, 6))
        ax.set_title(f"Prediction: y = {h_line}")

        # Plot the data
        red_dots = ax.plot(x, y, 'ro')[0]
        ax.set_xlabel('X Values') # X-axis label
        ax.set_ylabel('Y Values') # Y-axis label

        # Plot the prediction line
        yellow_line = ax.axhline(h_line, color='y')

        # Plot the errors
        predictions = h_line
        error_lines = ax.vlines(x, predictions, y, color='g') # Store error lines

        # Calculate the errors and their sum of squares
        errors = y - predictions
        #print("Errors (y - prediction):", errors) # Debugging print statement
        sse = np.dot(errors, errors)
        #print("Sum of Squared Errors (SSE):", sse) # Debugging print statement
        distance = np.sqrt(sse)
        #print("Distance (Square Root of SSE):", distance) # Debugging print statement for distance

        # Store results, including original points
        results.append((f"y = {h_line}", D.tolist(), errors.tolist(), errors.sum(), sse, distance))

    # Add a single legend for all subplots
    fig.legend([red_dots, yellow_line, error_lines], ['Data Points', 'Prediction Line', 'Errors'], loc='center left', bbox_to_anchor=(1, 0.5))

    col_labels = "Prediction Line", "Original Points", "Errors", "Sum", "SSE", "Distance"
    plt.tight_layout() # Adjust layout to prevent overlap
    plt.show()
    display(pd.DataFrame.from_records(results, columns=col_labels, index="Prediction Line"))

# Call the function to visualize
visualize_errors()

Expanding Linear Regression to Sloped Lines

In the exploration of linear regression, the restriction of drawing only horizontal lines is lifted, allowing the drawing of sloped lines. This development recognizes that a sloped line can represent a cluster of data points better if there is an overall growth or descent, akin to a plane taking off or landing.

The equation for the sloped line takes the classic form $ y = mx + b $, where the variables $ m $ (slope) and $ b $ (intercept) are adjusted to minimize the distance between the line and the data points. The distance from a point to the line can be represented as $ \text{distance}(m \cdot x + b, y) $, and the total distance across all data is the sum of these individual distances.

Even though a sloped line is considered, there may still be scenarios where a horizontal line is the best fit for the data. In that case, the slope $ m $ would be set to zero.

def process(D, model, ax):
    # Extracting x (input) and y (actual output) from the data
    x, y = D[:, 0], D[:, 1]
    m, b = model # Slope and intercept of the model

    # Model equation as a string
    model_str = f"y = {m} * x + {b}"
    ax.set_title(model_str) # Setting the model equation as the title

    # Other styling and plotting settings
    axis_helper(ax, (0, 8))
    red_dots = ax.plot(x, y, 'ro')[0] # Plotting the data points

    # Plotting the prediction line
    helper_xs = np.array([0, 8])
    helper_line = m * helper_xs + b
    yellow_line = ax.plot(helper_xs, helper_line, color='y')[0]

    # Calculating and plotting the errors
    predictions = m * x + b
    errors = y - predictions # Calculating errors
    error_lines = ax.vlines(x, predictions, y, color='g')

    # Calculating sum of squared errors (SSE) and distance
    sse = np.dot(errors, errors)
    distance = np.sqrt(sse)

    # Returning the results including the original data points and model equation
    return (D, model_str, errors, errors.sum(), sse, distance), (red_dots, yellow_line, error_lines)

# Data and models definition
D = np.array([[3, 5], [4, 2]])
lines_mb = np.array([[1, 0], [1, 1], [1, 2], [-1, 8], [-3, 14]])

# Column labels for the result DataFrame including original data points and model equation
col_labels = ("Original Points", "Model", "Raw Errors", "Sum", "SSE", "Total Distance")

# Results collection and handle collection
results = []
handles = []
fig, axes = plt.subplots(1, 5, figsize=(12, 6))
records = [process(D, mod, ax) for mod, ax in zip(lines_mb, axes.flat)]
results, handle_list = zip(*records)
handles = handle_list[0]  # Assuming all handles are the same for all subplots

# Add a single legend for all subplots with all three handles
fig.legend(handles, ['Data Points', 'Prediction Line', 'Errors'], loc='center left', bbox_to_anchor=(0.85, 0.5))

# Adjust layout to prevent overlap
plt.tight_layout(rect=[0, 0, 0.85, 1])

# Displaying the result DataFrame
df = pd.DataFrame.from_records(results, columns=col_labels)
display(df)