πŸ›£οΈ Week 09 - Lab Roadmap

Introduction to anomaly detection in Python

Author

The DS202 Team

Published

17 Mar 2026

βš™οΈ Setup

Downloading the student notebook

Click on the below button to download the student notebook.

Loading libraries

Let’s start by loading the relevant libraries!

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.cluster import KMeans, DBSCAN
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor, NearestNeighbors
from sklearn.metrics import calinski_harabasz_score
from kneed import KneeLocator
from yellowbrick.cluster import KElbowVisualizer
from prince import FAMD
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
import plotly.io as pio
import plotly.graph_objects as go
from plotly.subplots import make_subplots

Downloading the data

Download the datasets we will use for this lab.

Use the links below to download this dataset:

πŸ₯… Learning objectives

  • Understand the concept of anomaly detection
  • Deepen understanding of dimensionality reduction using FAMD
  • Learn how to implement anomaly detection using various algorithms

Part I: Introducing a new data set (10 minutes)

In this lab, we will use outlier detection to deepen our appreciation of 2000s and 2010s pop music. Using data from Spotify, we have a list of features for 919 popular singles that were released in the 1990s. Features include:

  • artist: Name of the Artist.
  • song: Name of the Track.
  • duration_ms: Duration of the track in milliseconds.
  • explicit: The lyrics or content of a song or a music video contain one or more of the criteria which could be considered offensive or unsuitable for children.
  • year: Release Year of the track.
  • popularity: The higher the value the more popular the song is.
  • danceability: Danceability describes how suitable a track is for dancing based on a combination of musical elements including tempo, rhythm stability, beat strength, and overall - regularity. A value of 0.0 is least danceable and 1.0 is most danceable.
  • energy: Energy is a measure from 0.0 to 1.0 and represents a perceptual measure of intensity and activity.
  • key: The key the track is in. Integers map to pitches using standard Pitch Class notation. E.g. 0 = C, 1 = Cβ™―/Dβ™­, 2 = D, and so on. If no key was detected, the value is -1.
  • loudness: The overall loudness of a track in decibels (dB). Loudness values are averaged across the entire track and are useful for comparing relative loudness of tracks. Loudness is the quality of a sound that is the primary psychological correlate of physical strength (amplitude). Values typically range between -60 and 0 db.
  • mode: Mode indicates the modality (major or minor) of a track, the type of scale from which its melodic content is derived. Major is represented by 1 and minor is 0.
  • speechiness: Speechiness detects the presence of spoken words in a track. The more exclusively speech-like the recording (e.g. talk show, audio book, poetry), the closer to 1.0 the attribute value. Values above 0.66 describe tracks that are probably made entirely of spoken words. Values between 0.33 and 0.66 describe tracks that may contain both music and speech, either in sections or layered, including such cases as rap music. Values below 0.33 most likely represent music and other non-speech-like tracks.
  • acousticness: A confidence measure from 0.0 to 1.0 of whether the track is acoustic. 1.0 represents high confidence the track is acoustic.
  • instrumentalness: Predicts whether a track contains no vocals. β€œOoh” and β€œaah” sounds are treated as instrumental in this context. Rap or spoken word tracks are clearly β€œvocal”. The closer the instrumentalness value is to 1.0, the greater likelihood the track contains no vocal content. Values above 0.5 are intended to represent instrumental tracks, but confidence is higher as the value approaches 1.0.
  • liveness: Detects the presence of an audience in the recording. Higher liveness values represent an increased probability that the track was performed live. A value above 0.8 provides strong likelihood that the track is live. valence: A measure from 0.0 to 1.0 describing the musical positiveness conveyed by a track. Tracks with high valence sound more positive (e.g. happy, cheerful, euphoric), while tracks with low valence sound more negative (e.g. sad, depressed, angry).
  • tempo: The overall estimated tempo of a track in beats per minute (BPM). In musical terminology, tempo is the speed or pace of a given piece and derives directly from the average beat duration.
  • genre: Genre of the track.
# Load the data set

# Print the shape attribute

πŸ—£CLASSROOM DISCUSSION

(Your class teacher will mediate this discussion)

  • How would you explore this data if dimensionality reduction was not an option?

  • How would you answer: what are the most common types of songs one can find on this data set?

Part II: Factor Analysis of Mixed Data (FAMD) (20 minutes)

Let’s create a list of different musical attributes, and filter the data frame to only include said attributes:

# Define which columns are numeric and which are categorical

# Create a list of musical attributes

# Create a new filtered data frame

# Ensure categorical columns are treated as categorical

🎯 ACTION POINTS

  1. To try to make sense of the sheer number of combinations of attributes, we will run FAMD and apply it to our data set:
# Option 1: Using FAMD (recommended for mixed data)
# FAMD handles both numeric and categorical variables automatically
pipe = Pipeline([
    ("famd", FAMD(
        n_components = 10,
        n_iter=3,
        random_state=42,
        engine='sklearn'
    ))
])

# Call the fit_transform method on 'pipe' (create a new object called 'components')
  1. How much β€˜information’ are we keeping after compressing the data with FAMD?

If we take the first 5 principle components, we can explain around 70% of the variance among our features

# Get the explained variance from FAMD
# In prince, you need to calculate the proportion from eigenvalues
eigenvalues = pipe.named_steps["famd"].eigenvalues_
total_inertia = eigenvalues.sum()
explained_variance_ratio = eigenvalues / total_inertia

# Create a data frame with cumulative variance
variance_explained = pd.DataFrame({
    "component": range(1, 11),  # Changed to 11 to match 10 components
    "variance": np.cumsum(explained_variance_ratio[:11])
})

# Create the plot
  1. Let’s focus on the first two components, as they are common plotting practices.

Think of the plot below as a ~45% compressed 2D representation of our 10 musical attributes.

# The FAMD transformation already returns a DataFrame of component scores.
# We therefore simply copy it and rename the columns to more informative labels
# (FAMD1, FAMD2, …) rather than recreating the DataFrame, which could introduce
# alignment issues.

components_df = components.copy()
components_df.columns = [f"FAMD{i+1}" for i in range(components_df.shape[1])]

# Add information on artist and track to the data frame
output_df = pd.concat([components_df, hits[["song","artist"]].reset_index(drop=True)], axis=1)

# Set renderer to browser
pio.renderers.default = "browser"

# Your existing code...
fig = px.scatter(output_df, 
                 x="FAMD1", 
                 y="FAMD2",
                 hover_data={"artist": True, 
                            "song": True},
                 opacity=0.5)

fig.update_layout(
    xaxis_title="FAMD1",
    yaxis_title="FAMD2",
    width=900,
    height=700
)

fig.show()  # This will now open in your default browser

πŸ§‘β€πŸ« TEACHING MOMENT

Your class teacher will now guide the conversation and explain the plot below. If needed, they will recap how FAMD works.

# Get the FAMD model
famd_model = pipe.named_steps["famd"]

# Calculate correlations between original variables and components
list_dfs = []

for i in range(4):  # First 4 components
    component_values = components_df[f"FAMD{i+1}"]
    
    loadings = []
    for col in hits_attrs.columns:
        if hits_attrs[col].dtype in ['float64', 'int64']:
            # Numeric variables - calculate correlation
            corr = hits_attrs[col].corr(component_values)
        else:
            # Categorical variables - calculate correlation with numeric encoding
            corr = hits_attrs[col].astype('category').cat.codes.corr(component_values)
        loadings.append(corr)
    
    df = pd.DataFrame({
        "features": hits_attrs.columns,
        "loadings": loadings,
        "component": i + 1
    })
    list_dfs.append(df)

# Concatenate all dataframes
loadings_df = pd.concat(list_dfs, ignore_index=True)

# Create absolute loadings column
loadings_df["abs_loadings"] = np.abs(loadings_df["loadings"])

# Create the plot
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

for idx, component in enumerate([1, 2, 3, 4]):
    comp_data = loadings_df[loadings_df["component"] == component].sort_values("loadings")
    
    # Create horizontal bar plot with color based on absolute value
    bars = axes[idx].barh(comp_data["features"], comp_data["loadings"])
    
    # Color bars using viridis colormap
    norm = plt.Normalize(vmin=0, vmax=comp_data["abs_loadings"].max())
    colors = plt.cm.viridis(norm(comp_data["abs_loadings"]))
    
    for bar, color in zip(bars, colors):
        bar.set_color(color)
    
    # Formatting
    axes[idx].set_xlabel("Loadings")
    axes[idx].set_ylabel("Features")
    axes[idx].set_title(f"Component {component}")
    axes[idx].axvline(x=0, color='black', linestyle='-', linewidth=0.8)
    axes[idx].grid(axis='x', alpha=0.3)

plt.tight_layout()
plt.show()

πŸ—£οΈ Discussion:

  • How does Figure 3 help you interpret Figure 2?

  • How does the above help you think about the attributes of the most common type of songs?

Part III: Anomaly detection techniques (1 hour)

πŸ‘₯ IN PAIRS, go through the action points below and discuss your impressions and key takeaways.

🎯 ACTION POINTS

  1. Take a look at the clusters identified by DBSCAN! We will employ a method that can help you determine values for the epsilon neighbourhood and minimum samples hyperparameters. We adapted this code from here.
# Set min_samples equal to 2 times the number of dimensions
min_samples = 4

# Instantiate nearest neighbours model, setting n_neighbors equal to min_sample
nearest_neighbors = NearestNeighbors(n_neighbors=min_samples)

# Fit the model to the first two principle components
neighbors = nearest_neighbors.fit(output_df[["FAMD1","FAMD2"]])

# Extract the distances and indices from the nearest neighbours model
distances, indices = neighbors.kneighbors(output_df[["FAMD1","FAMD2"]])

# Sort the distances from the 4th dimension of the distances matrix
distances = np.sort(distances[:,min_samples-1], axis=0)
# Identify the knee point
i = np.arange(len(distances))
knee = KneeLocator(
    i,
    distances,
    S=1,
    curve="convex",
    direction="increasing",
    interp_method="polynomial"
)

eps = distances[knee.knee]

# Plot the k-distance graph
plt.figure(figsize=(8,6))
plt.plot(i, distances, color="steelblue", linewidth=2)

# Horizontal line showing the recommended epsilon
plt.axhline(eps, color="red", linestyle="--", label=f"Suggested eps β‰ˆ {eps:.3f}")

# Optional: vertical knee location
plt.axvline(knee.knee, color="grey", linestyle=":", alpha=0.7)

plt.xlabel("Points (sorted by distance)")
plt.ylabel(f"{min_samples}-NN distance")
plt.title("k-distance graph for choosing DBSCAN epsilon")
plt.legend()
plt.tight_layout()
plt.show()

print(f"We should set the epsilon neighbourhood value to ~ {np.round(eps,4)}!")
# Instantiate a DBSCAN model
dbscan = DBSCAN(eps = eps, min_samples = min_samples)

# Fit the model to the first two principle component features
_ = dbscan.fit(output_df[["FAMD1","FAMD2"]])

Alternatively try out your own values for eps and min_samples

# Code here
# Prepare the data
to_plot = output_df.copy()
to_plot["dbscan"] = [str(lab) for lab in dbscan.labels_]
to_plot["dbscan_outlier"] = np.where(to_plot["dbscan"] == "-1", "Yes", "No")

# Create the interactive plot
fig = px.scatter(to_plot, 
                 x="FAMD1",  # Changed from PC1 to FAMD1
                 y="FAMD2",  # Changed from PC2 to FAMD2
                 color="dbscan_outlier",
                 hover_data={"song": True,
                            "artist": True,
                            "dbscan_outlier": False,  # Hide from hover
                            "FAMD1": ":.3f",
                            "FAMD2": ":.3f"},
                 labels={"dbscan_outlier": "Outlier"},
                 color_discrete_map={"Yes": "#F8766D", "No": "#00BFC4"})  # ggplot2 default colors

# Customize layout
fig.update_layout(
    xaxis_title="FAMD1",
    yaxis_title="FAMD2",
    width=900,
    height=700,
    legend_title_text="Outlier"
)

# Customize hover template for cleaner tooltips
fig.update_traces(
    hovertemplate="<b>%{customdata[0]}</b><br>" +
                  "Artist: %{customdata[1]}<br>" +
                  "FAMD1: %{x:.3f}<br>" +
                  "FAMD2: %{y:.3f}<br>" +
                  "<extra></extra>"
)

fig.write_html("dbscan_outliers.html", auto_open=True)

πŸ—£ Discussion: How well do you think DBSCAN performs at anomaly detection on the two principle components?

🎯 ACTION POINTS

  1. Take a look at the clusters identified by k-means:

In Week 08, we used the elbow method to guide the choice of the number of clusters. Here we introduce an alternative metric, the Calinski–Harabasz (CH) score, which evaluates clustering quality by comparing between-cluster separation with within-cluster compactness. Higher scores indicate better-defined clusters so the number of clusters that maximises the CH score is typically considered the most appropriate.

# Instantiate a k-means model
model = KMeans(random_state=42)

# Instantiate a visualizer from the yellowbrick library
visualizer = KElbowVisualizer(model, k=(2,10), metric = "calinski_harabasz", timings=False)

# Fit the visualizer to the first two principle components
visualizer.fit(output_df[["FAMD1","FAMD2"]])

# Finalize and render the figure
visualizer.show() 
# Instantiate a model
kmeans = KMeans(n_clusters=3, random_state=42)

# Fit the model to the data
_ = kmeans.fit(output_df[["FAMD1","FAMD2"]])

# Create a data frame based on customers by adding the cluster labels
to_plot["kmeans"] = [str(i) for i in kmeans.labels_]

# Create the plot

πŸ—£ Discussion: How well do you think k-means performs at anomaly detection on the two principal components?

🎯 ACTION POINTS

  1. Take a look at the clusters identified by the isolation forest:
# Instantiate a model
pipe_iso = Pipeline([("scaler", StandardScaler()), ("isoforest", IsolationForest(random_state=123))])

# Fit model to training data
pipe_iso.fit(output_df[["FAMD1","FAMD2"]])

# Calculate the anomaly scores for the same data frame
to_plot["isoforest"] = pipe_iso.score_samples(output_df[["FAMD1","FAMD2"]])

# Thresholds to try out
iso_ths = [-0.7, -0.65, -0.6, -0.55]

# Create variables that exceed each threshold
for th in iso_ths:
    to_plot[f"iso_th_{th}"] = to_plot["isoforest"] <= th

# Keep only the track, artist, first two components and various isoforest threshold variables
feats_to_keep = to_plot.columns[to_plot.columns.str.contains("song|artist|FAMD[1-2]|iso_th")]

# Create a melted data frame to plot
to_plot_melted = (
    to_plot
    .filter(items = feats_to_keep)
    .melt(id_vars = feats_to_keep[feats_to_keep.str.contains("song|artist|FAMD")])
    .rename(columns = {"variable":"th", "value": "isoforest_outlier"})
    .assign(th = lambda x: x["th"].str.replace("iso_th_","Threshold = "))
)

# Convert boolean to Yes/No for better display
to_plot_melted["isoforest_outlier"] = to_plot_melted["isoforest_outlier"].map({True: "Yes", False: "No"})

# Create the plot

πŸ—£ Discussion: What is the relationship between the anomaly score and the number of outliers in the data?

🎯 ACTION POINTS

  1. Let’s see if the Local Outlier Factor (LOF) performs better than DBSCAN/Isolation Forest.

We use the LocalOutlierFactor() function to calculate local outlier factors:

# Instantiate a Local Outlier Factor model, setting nearest neighbors to 10
lof = LocalOutlierFactor(n_neighbors = 10)

# Fit the model to the first two principle components
lof.fit(output_df[["FAMD1","FAMD2"]])

# Append the negative outlier factor score to the data frame, using absolute values
to_plot["lof"] = np.abs(lof.negative_outlier_factor_)
# Create the plot

πŸ—£ Discussion: Does LOF perform better than DBSCAN or isolation forests to detect β€˜anomalous’ samples?