TDM 30100: Project 5 - Linear Regression Project

import pandas as pd

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

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"]

spotify_popularity_data = spotify_popularity_data.drop(columns=drop_cols)

# For YOU to do: List columns still in spotify_popularity_data after removing drop_cols

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

# Print shape of X and y
print(_____) # For YOU to do
print(____) #For YOU to do

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)

# For YOU to do: print X_train shape

# For YOU to do: print X_test shape

# For YOU to do: print y_train shape

# For YOU to do: print y_test shape

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("_____") # For YOU to fill
plt.ylabel("______") # For YOU to fill
plt.title("_____") # For YOU to fill
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("______") # For YOU to fill in
plt.ylabel("______")  # For YOU to fill in
plt.title("_______") # For YOU to fill in
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=__): # For YOU to fill in
    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]

# For YOU to do: print 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).")


# STUDENT TODO
X_train = X_train.drop(index=________) # For YOU to fill in
y_train = y_train.drop(index=________) # For YOU to fill in

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 model

print(_______) # For YOU to do

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("________") # For YOU to fill in
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()

import matplotlib.pyplot as plt

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