Linear Regression with Spotify Data

import pandas as pd

spotify_popularity_data = pd.read_csv("/anvil/projects/tdm/data/spotify/linear_regression_popularity.csv")

spotify_popularity_data.head()

drop_cols = [
    "Unnamed: 0", "Unnamed: 0.1", "track_id", "track_name", "available_markets", "href",
    "album_id", "album_name", "album_release_date", "album_type",
    "artists_names", "artists_ids", "principal_artist_id",
    "principal_artist_name", "artist_genres", "analysis_url", "duration_min"]

# Drop unnecessary columns
spotify_popularity_data = spotify_popularity_data.drop(columns=drop_cols)

for col in spotify_popularity_data:
    print(col)

popularity
disc_number
duration_ms
explicit
track_number
album_total_tracks
principal_artist_followers
acousticness
danceability
energy
instrumentalness
key
liveness
loudness
mode
speechiness
tempo
time_signature
valence
year

# Target and features
y = spotify_popularity_data["popularity"].copy()
X = spotify_popularity_data.drop(columns=["popularity"]).copy()

# Quick check
print(X.shape)
print(y.shape)

(1500, 19)
(1500,)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

print("X_train shape:", X_train.shape)
print("X_test  shape:", X_test.shape)
print("y_train shape:", y_train.shape)
print("y_test  shape:", y_test.shape)

X_train shape: (1200, 19)
X_test  shape: (300, 19)
y_train shape: (1200,)
y_test  shape: (300,)

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(8,5))
sns.histplot(y, bins=30, kde=True, color="skyblue")
plt.xlabel("Popularity")
plt.ylabel("Count")
plt.title("Distribution of Track Popularity")
plt.show()

import matplotlib.pyplot as plt
import seaborn as sns

# Convert duration to minutes for training data
duration_min_train = X_train["duration_ms"] / 60000

plt.figure(figsize=(8,5))
sns.scatterplot(x=duration_min_train, y=y_train, alpha=0.6)
sns.regplot(x=duration_min_train, y=y_train, scatter=False, color="red", ci=None)

plt.xlabel("Duration (minutes)")
plt.ylabel("Popularity") )
plt.title("Popularity vs. Duration (train set)") )
plt.show()

import pandas as pd
import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor

# Convert booleans to ints
bool_cols = X_train.select_dtypes(include=["bool"]).columns
if len(bool_cols):
    X_train[bool_cols] = X_train[bool_cols].astype(int)


def calculate_vif_iterative(X, thresh=10.0):
    X_ = X.astype(float).copy()
    while True:
        vif_df = pd.DataFrame({
            "variable": X_.columns,
            "VIF": [variance_inflation_factor(X_.values, i) for i in range(X_.shape[1])]
        }).sort_values("VIF", ascending=False).reset_index(drop=True)

        max_vif = vif_df["VIF"].iloc[0]
        worst = vif_df["variable"].iloc[0]

        if (max_vif <= thresh) or (X_.shape[1] <= 1):
            return X_, vif_df.sort_values("VIF")

        print(f"Dropping '{worst}' with VIF={max_vif:.2f}")
        X_ = X_.drop(columns=[worst])

# Run iterative VIF filtering
result_vif = calculate_vif_iterative(X_train, thresh=10.0)

# Split into the filtered dataset and the VIF summary
X_train = result_vif[0]
vif_summary = result_vif[1]

Dropping 'year' with VIF=320.29
Dropping 'time_signature' with VIF=92.09
Dropping 'disc_number' with VIF=46.33
Dropping 'energy' with VIF=19.81
Dropping 'danceability' with VIF=19.10
Dropping 'tempo' with VIF=12.17

vif_summary

import numpy as np
import pandas as pd
import statsmodels.api as sm

# Align and clean training data
X_train_cook = X_train.loc[X_train.index.intersection(y_train.index)]
y_train_cook = y_train.loc[X_train_cook.index]

# Keep only rows without missing/infinite values
mask = np.isfinite(X_train_cook).all(1) & np.isfinite(y_train_cook.to_numpy())
X_train_cook, y_train_cook = X_train_cook.loc[mask], y_train_cook.loc[mask]

# Fit model on the cleaned data
ols_model_cook = sm.OLS(y_train_cook, sm.add_constant(X_train_cook, has_constant="add")).fit()

# Cook's Distance values for each observation
cooks_distance = ols_model_cook.get_influence().cooks_distance[0]
cooks_threshold = 4 / len(X_train_cook)

# Identify outlier indices
cooks_outliers = X_train_cook.index[cooks_distance > cooks_threshold]

print(f"Flagged {len(cooks_outliers)} outliers (Cook's D > 4/n).")


X_train = X_train.drop(index=cooks_outliers)
y_train = y_train.drop(index=cooks_outliers)

Flagged 61 outliers (Cook's D > 4/n).

import statsmodels.api as sm

def stepwise_aic(X_train, y_train, max_vars=None, verbose=True):
    X = X_train.select_dtypes(include=[np.number]).astype(float).copy()
    y = pd.Series(y_train).astype(float)
    mask = np.isfinite(X).all(1) & np.isfinite(y)
    X, y = X.loc[mask], y.loc[mask]

    remaining, selected = set(X.columns), []
    current_aic = np.inf

    while remaining and (max_vars is None or len(selected) < max_vars):
        aics = [(sm.OLS(y, sm.add_constant(X[selected + [c]], has_constant='add')).fit().aic, c)
                for c in remaining]
        best_aic, best_var = min(aics, key=lambda t: t[0])

        if best_aic + 1e-6 < current_aic:
            selected.append(best_var); remaining.remove(best_var); current_aic = best_aic
            if verbose: print(f"+ {best_var} (AIC {best_aic:.2f})")
            # backward
            improved = True
            while improved and len(selected) > 1:
                backs = [(sm.OLS(y, sm.add_constant(X[[v for v in selected if v != d]], has_constant='add')).fit().aic, d)
                         for d in selected]
                back_aic, drop_var = min(backs, key=lambda t: t[0])
                if back_aic + 1e-6 < current_aic:
                    selected.remove(drop_var); remaining.add(drop_var); current_aic = back_aic
                    if verbose: print(f"- {drop_var} (AIC {back_aic:.2f})")
                else:
                    improved = False
        else:
            break

    model = sm.OLS(y, sm.add_constant(X[selected], has_constant='add')).fit()
    return selected, model

results_feature_selection = stepwise_aic(X_train, y_train, max_vars=None, verbose=True)

selected_cols = results_feature_selection[0]   # list of features
model = results_feature_selection[1]           # fitted OLS model

print(selected_cols)

+ principal_artist_followers (AIC 8054.83)
+ duration_ms (AIC 8001.12)
+ loudness (AIC 7961.86)
+ album_total_tracks (AIC 7932.94)
+ acousticness (AIC 7909.56)
+ explicit (AIC 7887.99)
+ mode (AIC 7873.37)
+ valence (AIC 7861.57)
+ track_number (AIC 7858.87)
+ key (AIC 7857.61)
+ speechiness (AIC 7857.16)

print(model.summary())

import matplotlib.pyplot as plt
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np

# Residuals
resid = model.resid
z = (resid - np.mean(resid)) / np.std(resid, ddof=1)

# Histogram with normal curve
plt.hist(z, bins=30, density=True, alpha=0.6, color="skyblue", edgecolor="black")
x = np.linspace(z.min(), z.max(), 100)
plt.plot(x, stats.norm.pdf(x, 0, 1), "r", linewidth=2)
plt.title("Histogram of Standardized Residuals with Normal Curve")
plt.show()

# Q–Q Plot
sm.qqplot(z, line="45", fit=True)
plt.title("Q–Q Plot of Residuals")
plt.show()

# Durbin–Watson independence test
from statsmodels.stats.stattools import durbin_watson

dw_stat = durbin_watson(resid)
print("[Independence: Durbin–Watson]")
print(f"Statistic={dw_stat:.4f}")
print("H0: Residuals are independent (no autocorrelation).")

if 1.5 < dw_stat < 2.5:
    print("-> Pass (approx. independent)")
else:
    print("-> FAIL: possible autocorrelation")
print()

[Independence: Durbin–Watson]
Statistic=1.9603
H0: Residuals are independent (no autocorrelation).
-> Pass (approx. independent)

import matplotlib.pyplot as plt

plt.scatter(model.fittedvalues, resid, alpha=0.6)
plt.axhline(0, color="red", linestyle="--")
plt.title("Residuals vs. Fitted Values")
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.show()

from sklearn.metrics import r2_score
import statsmodels.api as sm

# Making sure test data has the same selected columns
X_test = X_test[selected_cols].copy()
mask = np.isfinite(X_test).all(1)
X_test = X_test.loc[mask]
y_test = y_test.loc[X_test.index]

# Predict with statsmodels OLS
y_pred = model.predict(sm.add_constant(X_test, has_constant='add'))

# R² on the test set
r2 = r2_score(y_test, y_pred)
print(f"Test R²: {r2:.3f}")

Test R²: 0.250

import pandas as pd

alpha = 0.05  # significance threshold

# Collect pieces from the statsmodels OLS result
coef = model.params.rename("coef")
pval = model.pvalues.rename("p_value")
ci = model.conf_int()

# Assemble table
coef_tbl = (
    pd.concat([coef, pval, ci], axis=1)
      .drop(index="const", errors="ignore")             # drop intercept row
      .assign(
          effect=lambda d: d["coef"].apply(lambda x: "Positive" if x > 0 else "Negative"),
          significant=lambda d: d["p_value"] < alpha
      )
      [["coef", "effect", "p_value", "significant"]]
      .sort_values("p_value")                            # most significant first
)

# rounding
coef_tbl_rounded = coef_tbl.round({"coef": 4, "p_value": 4})
print(coef_tbl_rounded)

 coef    effect  p_value  significant
principal_artist_followers  0.0000  Positive   0.0000         True
loudness                    0.6164  Positive   0.0000         True
duration_ms                -0.0000  Negative   0.0000         True
acousticness                5.4814  Positive   0.0000         True
album_total_tracks         -0.1788  Negative   0.0000         True
explicit                    2.7251  Positive   0.0000         True
mode                       -1.9647  Negative   0.0001         True
valence                    -3.5808  Negative   0.0002         True
track_number               -0.1283  Negative   0.0300         True
key                         0.1134  Positive   0.0777        False
speechiness                -4.6218  Negative   0.1195        False

import pandas as pd
import statsmodels.api as sm

# Students will fill these in with the new song's known features
new_song = pd.DataFrame([{
    "principal_artist_followers": 5_000_000,
    "duration_ms": 210_000,
    "loudness": -6.0,
    "album_total_tracks": 10,
    "acousticness": 0.20,
    "explicit": 0,
    "mode": 1,
    "valence": 0.40,
    "track_number": 3,
    "key": 5,
    "speechiness": 0.05
}])

# Ensure the model sees the same columns it was trained on
new_song = new_song.reindex(columns=selected_cols)
new_song = new_song.fillna(X_train[selected_cols].median(numeric_only=True))
pred_pop = model.predict(sm.add_constant(new_song, has_constant="add"))[0]
print(f"Predicted popularity: {pred_pop:.1f}")

Predicted popularity: 68.9