←Back to the bookChapter 27 · Case Study: Spotify Track Features

Chapter 27 · Case Study

Spotify Track Features 🎧

One messy file of 915 songs, taken from raw to analysis-ready. We audit, describe, visualize, then clean every problem in turn, explaining the reason for each decision.

Statistics, Data Science and AI: A Visual Handbook · John Fisher · 2026

⚙️ Setup & load

The data loads from the book's data/ folder (local) with a GitHub raw-URL fallback (Colab).

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
CYAN="#0891b2"; AMBER="#d97706"; PURPLE="#7c3aed"; GREEN="#059669"; PINK="#db2777"; INK="#1a2138"; GRID="#e6e9f2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
pd.set_option("display.max_columns", 30)

BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try:    df = pd.read_csv("../../data/spotify_tracks.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"spotify_tracks.csv")
raw_rows = len(df)
print("loaded:", df.shape)

loaded: (915, 15)

BEAT 1 · FIRST CONTACT

Audit before anything else

Look at the raw file: shape, dtypes, duplicates, missingness, and how the categories were actually typed in. We fix nothing yet, we just find out what we are dealing with.

print("shape:", df.shape)
print("\nduplicate rows:", int(df.duplicated().sum()))
miss = df.isna().sum(); print("\nmissing:\n", miss[miss>0].to_string())
print("\ngenre typed", df["genre"].nunique(), "ways:", sorted(df["genre"].dropna().unique()))
print("\nmode typed", df["mode"].nunique(), "ways:", sorted(df["mode"].dropna().unique()))
print("\ntempo == 0 (impossible):", int((df["tempo"]==0).sum()))
df.head(3)

shape: (915, 15)

duplicate rows: 8

missing:
 tempo         18
key           12
popularity    73

genre typed 19 ways: ['EDM', 'Electronic', 'Hip-Hop', 'POP', 'Pop', 'R&B', 'classical', 'country', 'electronic', 'hip hop', 'hip-hop', 'hiphop', 'jazz', 'pop', 'pop ', 'r and b', 'r&b', 'rnb', 'rock']

mode typed 8 ways: ['M', 'Maj', 'Major', 'Min', 'Minor', 'm', 'major', 'minor']

tempo == 0 (impossible): 23

	track_id	track_name	artist	genre	danceability	energy	loudness	tempo	valence	duration_ms	key	mode	popularity	audio_quality	mood
0	TRK00063	Silent Summer	Foxglove	pop	0.890	0.707	-21.31	143.448	0.436	215091	0.0	Major	10.0	standard	happy
1	TRK00677	Velvet Gravity	The Paper Kites	rock	0.532	0.808	-16.58	117.839	0.274	280893	4.0	Major	20.0	low	sad
2	TRK00140	Broken Waves	Northbound	r&b	0.722	0.368	-38.31	136.771	0.476	210682	7.0	Major	38.0	high	chill

What the audit reveals: eight duplicate rows; missing tempo, key, and popularity; a genre column typed nineteen different ways for eight real genres; mode typed eight ways for two; and a tempo of 0 BPM, which is impossible and therefore missing data in disguise. That is the to-do list.

BEAT 2 · DESCRIBE

Numbers first

Descriptive statistics on the numeric columns. The headline is the scales: these features live in completely different number ranges, which will matter the moment we model.

num = ["danceability","energy","loudness","tempo","valence","duration_ms","popularity"]
desc = df[num].describe().loc[["mean","std","min","max"]].round(2)
print(desc.to_string())
print("\nNotice the ranges: danceability 0-1, loudness -40-0 dB, tempo 0-300,")
print("duration_ms in the hundreds of thousands. Four different worlds.")

      danceability  energy  loudness   tempo  valence  duration_ms  popularity
mean          0.56    0.56    -27.63  115.30     0.52    239581.64       28.67
std           0.19    0.22     10.32   32.68     0.21    595349.75       16.69
min           0.08    0.02    -40.00    0.00     0.01     97419.00        0.00
max           0.95    0.98     -0.25  299.00     1.00  18000000.00       81.00

Notice the ranges: danceability 0-1, loudness -40-0 dB, tempo 0-300,
duration_ms in the hundreds of thousands. Four different worlds.

BEAT 3 · VISUALIZE

Let the charts talk

Three pictures: the distribution of popularity, a correlation heatmap of the audio features, and the one strong relationship it reveals. (We plot the raw data first, then re-plot after cleaning.)

fig, ax = plt.subplots(1, 2, figsize=(11,3.8))
ax[0].hist(df["popularity"].dropna(), bins=24, color=CYAN, edgecolor="white")
ax[0].set_title(f"popularity (skew {df['popularity'].skew():+.2f})"); ax[0].set_yticks([])
# correlation heatmap
corr = df[num].corr()
im = ax[1].imshow(corr, cmap="coolwarm", vmin=-1, vmax=1)
ax[1].set_xticks(range(len(num))); ax[1].set_xticklabels(num, rotation=45, ha="right", fontsize=8)
ax[1].set_yticks(range(len(num))); ax[1].set_yticklabels(num, fontsize=8)
for i in range(len(num)):
    for j in range(len(num)):
        ax[1].text(j,i,f"{corr.iloc[i,j]:.2f}",ha="center",va="center",fontsize=7,
                   color="white" if abs(corr.iloc[i,j])>0.6 else "#1a2138")
ax[1].set_title("audio-feature correlations"); ax[1].grid(False)
plt.colorbar(im, ax=ax[1], fraction=0.046); plt.tight_layout(); plt.show()

r = df["energy"].corr(df["loudness"])
fig, ax = plt.subplots(figsize=(6,4))
ax.scatter(df["energy"], df["loudness"], s=14, color=PURPLE, alpha=0.5, edgecolor="white", linewidth=0.3)
ax.set_xlabel("energy"); ax.set_ylabel("loudness (dB)")
ax.set_title(f"energy vs loudness:  r = {r:.2f}"); plt.tight_layout(); plt.show()
print(f"energy and loudness move together almost perfectly (r = {r:.2f}).")

energy and loudness move together almost perfectly (r = 0.94).

What the visuals say: popularity is mildly right-skewed; almost every audio feature is uncorrelated with every other, with one loud exception, energy and loudness correlate at r ≈ 0.94. That near-duplication is multicollinearity: the two columns carry nearly the same information. Worth knowing now, because we may later drop or combine one of them.

BEAT 4 · CLEAN

Duplicates, categories, missing values

Now we fix what the audit found, in order. Deduplicate (Ch 19), standardize the categories (Ch 19), then handle the missing values (Ch 20), each with its reason.

clean = df.drop_duplicates().copy()
print(f"deduplicate: {raw_rows} -> {len(clean)} rows  ({raw_rows-len(clean)} exact duplicates removed)")

# standardize genre (19 spellings -> 8) and mode (8 -> 2)
gmap = {"pop":"pop","Pop":"pop","POP":"pop","pop ":"pop",
        "hip-hop":"hip-hop","hip hop":"hip-hop","Hip-Hop":"hip-hop","hiphop":"hip-hop",
        "r&b":"r&b","R&B":"r&b","rnb":"r&b","r and b":"r&b",
        "electronic":"electronic","Electronic":"electronic","EDM":"electronic"}
clean["genre"] = clean["genre"].str.strip().map(lambda g: gmap.get(g, g))
mmap = {"Major":"Major","major":"Major","M":"Major","Maj":"Major",
        "Minor":"Minor","minor":"Minor","m":"Minor","Min":"Minor"}
clean["mode"] = clean["mode"].map(lambda m: mmap.get(m, m))
# two more messy categoricals: an ORDINAL (audio_quality) and a NOMINAL (mood)
clean["audio_quality"] = clean["audio_quality"].str.strip().str.lower().replace({"std":"standard"})
clean["mood"] = clean["mood"].str.strip().str.lower()
print("genres now:", clean["genre"].nunique(), "| modes now:", sorted(clean["mode"].dropna().unique()))
print("audio_quality:", sorted(clean["audio_quality"].unique()), "| mood:", sorted(clean["mood"].unique()))

deduplicate: 915 -> 907 rows  (8 exact duplicates removed)
genres now: 8 | modes now: ['Major', 'Minor']
audio_quality: ['high', 'low', 'standard'] | mood: ['chill', 'happy', 'sad']

# missing values, each treated on its merits
# 1) tempo: 0 is impossible -> recode to NaN, then median-impute (a measurement gap)
clean["tempo"] = clean["tempo"].replace(0, np.nan)
tempo_med = clean["tempo"].median()
n_tempo = int(clean["tempo"].isna().sum())
clean["tempo"] = clean["tempo"].fillna(tempo_med)
# 2) key: nominal (musical key) -> impute the mode
clean["key"] = clean["key"].fillna(clean["key"].mode()[0])
# 3) popularity is the OUTCOME we care about -> do not invent it; drop those rows
n_pop = int(clean["popularity"].isna().sum())
clean = clean.dropna(subset=["popularity"])
print(f"tempo: recoded 0->NaN and median-imputed {n_tempo} values (median {tempo_med:.1f})")
print(f"key:   mode-imputed the missing musical keys")
print(f"popularity: dropped {n_pop} rows missing the outcome (never impute a target)")
print("rows now:", len(clean))

tempo: recoded 0->NaN and median-imputed 41 values (median 117.6)
key:   mode-imputed the missing musical keys
popularity: dropped 73 rows missing the outcome (never impute a target)
rows now: 834

Why these three differ: a missing tempo is a measurement gap, so we estimate it (median, which resists outliers). A missing key is a category, so we fill the most common one. But popularity is the thing we ultimately want to explain, and inventing values for your outcome biases everything downstream, so those rows are dropped instead. Same problem, three different right answers (Chapter 20).

BEAT 5 · OUTLIERS

Find the impossible durations

A box plot of duration_ms screams: a handful of tracks are minutes, even hours, too long. We separate genuine long songs from data errors and remove only the errors (Ch 21).

dur_min = clean["duration_ms"]/60000
q1,q3 = clean["duration_ms"].quantile([.25,.75]); iqr=q3-q1; fence=q3+1.5*iqr
print(f"IQR upper fence: {fence/60000:.1f} min;  rows above it: {int((clean['duration_ms']>fence).sum())}")
errors = clean[clean["duration_ms"]>600000]   # > 10 minutes = almost certainly corrupt
print("impossible durations (minutes):", (errors["duration_ms"]/60000).round(1).to_list())

fig, ax = plt.subplots(figsize=(7,3.4))
ax.boxplot(dur_min, orientation="horizontal", widths=0.6)
ax.scatter(errors["duration_ms"]/60000, [1]*len(errors), color=PINK, zorder=5, s=40, label="data errors")
ax.set_xlabel("duration (minutes)"); ax.set_title("duration: 3 impossible values"); ax.legend()
ax.set_yticks([]); plt.tight_layout(); plt.show()

before = len(clean); clean = clean[clean["duration_ms"]<=600000].copy()
print(f"removed {before-len(clean)} corrupt rows; {len(clean)} remain")

IQR upper fence: 6.3 min;  rows above it: 19
impossible durations (minutes): [300.0, 25.0, 32.0]

removed 3 corrupt rows; 831 remain

Remove vs keep: the IQR rule flags a couple dozen long tracks, but most are legitimate (extended mixes, classical pieces). Only three are impossible, up to five hours, so those are the data errors we delete. We keep the genuine long songs: an outlier is not automatically an error, and judgment, not a blind rule, decides (Chapter 21).

BEAT 6 · TRANSFORM, SCALE & ENCODE

Make it model-ready

The shaping steps. A log transform straightens the skewed duration (Ch 22); scaling puts every feature on the same footing; and we encode by type (Ch 24): the ordinal audio_quality as 1/2/3, the nominal genre and mood as one-hot columns.

from sklearn.preprocessing import StandardScaler

# log transform the skewed duration
before_sk = clean["duration_ms"].skew()
clean["log_duration"] = np.log(clean["duration_ms"])
after_sk = clean["log_duration"].skew()
fig, ax = plt.subplots(1,2,figsize=(10,3.4))
ax[0].hist(clean["duration_ms"]/60000, bins=30, color=AMBER, edgecolor="white")
ax[0].set_title(f"duration, raw (skew {before_sk:+.2f})"); ax[0].set_xlabel("minutes"); ax[0].set_yticks([])
ax[1].hist(clean["log_duration"], bins=30, color=GREEN, edgecolor="white")
ax[1].set_title(f"log(duration) (skew {after_sk:+.2f})"); ax[1].set_yticks([])
plt.tight_layout(); plt.show()
print(f"log transform pulled skew from {before_sk:+.2f} to {after_sk:+.2f}: near-symmetric.")

log transform pulled skew from +0.93 to +0.04: near-symmetric.

# scale the audio features so none dominate by sheer magnitude
feat = ["danceability","energy","loudness","tempo","valence","log_duration"]
scaled = pd.DataFrame(StandardScaler().fit_transform(clean[feat]), columns=feat, index=clean.index)
print("BEFORE scaling (std of each feature):")
print(clean[feat].std().round(2).to_string())
print("\nAFTER standardizing (every std = 1, every mean = 0):")
print(scaled.std().round(2).to_string())

# encode the categoricals by TYPE (Ch 24)
# nominal genre -> one-hot (no order)
genre_oh = pd.get_dummies(clean["genre"], prefix="genre")
# ORDINAL audio_quality -> integer codes that PRESERVE the order low < standard < high
qorder = {"low":1, "standard":2, "high":3}
clean["audio_quality_code"] = clean["audio_quality"].map(qorder)
# NOMINAL mood -> one-hot (happy/sad/chill have no order)
mood_oh = pd.get_dummies(clean["mood"], prefix="mood")
print("ordinal map :", qorder)
print("mood one-hot:", list(mood_oh.columns))
print("\nmean popularity by audio_quality (the ordinal tracks it):")
print(clean.groupby("audio_quality", observed=True)["popularity"].mean()
      .reindex(["low","standard","high"]).round(1).to_string())

BEFORE scaling (std of each feature):
danceability     0.19
energy           0.22
loudness        10.35
tempo           26.73
valence          0.21
log_duration     0.29

AFTER standardizing (every std = 1, every mean = 0):

danceability    1.0
energy          1.0
loudness        1.0
tempo           1.0
valence         1.0
log_duration    1.0
ordinal map : {'low': 1, 'standard': 2, 'high': 3}
mood one-hot: ['mood_chill', 'mood_happy', 'mood_sad']

mean popularity by audio_quality (the ordinal tracks it):
audio_quality
low         19.3
standard    26.8
high        40.1

# assemble the analysis-ready table
ready = pd.concat([scaled,
                   clean[["popularity","audio_quality_code"]].set_index(scaled.index),
                   genre_oh, mood_oh], axis=1)
print("ANALYSIS-READY:", ready.shape)
print(f"\njourney: {raw_rows} raw rows  ->  {len(ready)} clean rows, {ready.shape[1]} model-ready columns")
ready.head(3)

ANALYSIS-READY: (831, 19)

journey: 915 raw rows  ->  831 clean rows, 19 model-ready columns

	danceability	energy	loudness	tempo	valence	log_duration	popularity	audio_quality_code	genre_classical	genre_country	genre_electronic	genre_hip-hop	genre_jazz	genre_pop	genre_r&b	genre_rock	mood_chill	mood_happy	mood_sad
0	1.728965	0.647654	0.601988	0.938509	-0.396648	0.128207	10.0	2	False	False	False	False	False	True	False	False	False	True	False
1	-0.135705	1.105685	1.059423	-0.020261	-1.165114	1.064286	20.0	1	False	False	False	False	False	False	False	True	False	False	True
2	0.853924	-0.889699	-1.042067	0.688530	-0.206904	0.055571	38.0	3	False	False	False	False	False	False	True	False	True	False	False

BEAT 7 · EXPLORE THE CLEAN DATA

Now we can trust what we see

With the file clean, the describing (Part II) and visualizing (Part III) toolkits finally tell the truth. Measurements first, then the charts that were not trustworthy on the raw data.

# descriptive measurements (Ch 8-11) on the cleaned columns
def measure(x):
    return pd.Series({
        "mean":   x.mean(),   "median": x.median(), "mode": x.round().mode().iloc[0],
        "SD":     x.std(),    "IQR":    x.quantile(.75)-x.quantile(.25),
        "CV %":   x.std()/x.mean()*100,
        "Q1":     x.quantile(.25), "Q3": x.quantile(.75),
        "skew":   x.skew(),   "kurtosis": x.kurtosis()})
M = pd.DataFrame({"tempo": measure(clean["tempo"]),
                  "popularity": measure(clean["popularity"])}).round(2)
print(M.to_string())

           tempo  popularity
mean      118.38       28.68
median    117.56       26.00
mode      118.00       22.00
SD         26.73       16.66
IQR        32.69       24.00
CV %       22.58       58.09
Q1        101.43       16.00
Q3        134.12       40.00
skew        0.65        0.55
kurtosis    3.77       -0.37

Reading the measurements: popularity is far more variable than tempo relative to its center (CV 58% vs 23%, Ch 9), so songs differ much more in popularity than in tempo. Both are mildly right-skewed (mean > median, Ch 11). And tempo is leptokurtic (kurtosis ≈ +3.8): a sharp peak around 118 BPM with heavier tails, while popularity is flatter. The five-number summary (min, Q1, median, Q3, max) frames where the bulk of each variable sits (Ch 10).

# Ch 14 (categorical) + Ch 16 (grouped): genre counts, and the folklore check
gc = clean["genre"].value_counts()
vm = clean.groupby("mode")["valence"].mean()
fig, ax = plt.subplots(1, 2, figsize=(11,3.8))
ax[0].bar(gc.index, gc.values, color=CYAN, edgecolor="white")
ax[0].set_title("genre distribution (clean)"); ax[0].set_ylabel("songs")
ax[0].tick_params(axis="x", rotation=35)
ax[1].bar(vm.index, vm.values, color=[AMBER, PURPLE], edgecolor="white", width=0.55)
ax[1].set_ylim(0,0.7); ax[1].set_title("mean valence: Major vs Minor")
for i,v in enumerate(vm.values): ax[1].text(i, v+0.02, f"{v:.3f}", ha="center", fontweight="bold")
plt.tight_layout(); plt.show()
print("Folklore says major keys sound happier. Here the gap is tiny:", round(vm.max()-vm.min(),3))

Folklore says major keys sound happier. Here the gap is tiny: 0.007

# Ch 15/16: box plots of energy by genre reveal a pattern the raw data hid
order = clean.groupby("genre")["energy"].median().sort_values(ascending=False).index
data  = [clean.loc[clean["genre"]==g, "energy"].values for g in order]
hi = {"pop","rock","electronic"}
fig, ax = plt.subplots(figsize=(9,4))
bp = ax.boxplot(data, patch_artist=True, tick_labels=list(order))
for patch, g in zip(bp["boxes"], order):
    patch.set_facecolor(CYAN if g in hi else PURPLE); patch.set_alpha(0.55)
ax.set_ylabel("energy"); ax.set_title("energy by genre: two natural clusters")
ax.tick_params(axis="x", rotation=25); plt.tight_layout(); plt.show()
print("High-energy (cyan): pop / rock / electronic ~0.71  |  mellow (purple): the rest ~0.40")

High-energy (cyan): pop / rock / electronic ~0.71  |  mellow (purple): the rest ~0.40

What the clean charts reveal: genres are reasonably balanced (Ch 14); the major-vs-minor "happier key" folklore is not supported here (the valence gap is about 0.006, Ch 16); and energy splits the genres into two clear clusters, an upbeat group (pop, rock, electronic) and a mellow group, a real structure that was not safe to read until the duplicates, missing values, and errors were gone (Ch 15/16). Exploration on clean data is where understanding actually happens.

---

🎯 Recap: what we did and why

• Deduplicated 8 exact-copy rows (Ch 19), they would double-count songs.
• Standardized genre (19 spellings → 8) and mode (8 → 2) so groups are real (Ch 19).
• Recoded tempo 0 → NaN and median-imputed; mode-imputed key; dropped rows missing the outcome popularity (Ch 20).
• Removed 3 impossible durations while keeping genuine long songs (Ch 21).
• Log-transformed the skewed duration, scaled all features, and one-hot encoded genre (Ch 22 & 24).

Along the way the data told us something: energy and loudness are near-duplicates (r ≈ 0.94), and audio features barely predict popularity here. That is a finding, not a cleaning step, and it is exactly what exploratory analysis is for.

Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher

←Back to the bookChapter 27 · Case Study: Spotify Track Features