<jupyter_start><jupyter_text>Machine Learning analysis- This is a Python base notebook  Kaggle's [Spotify Song Attributes](https://www.kaggle.com/geomack/spotifyclassification/home) dataset contains a number of features of songs from 2017 and a binary variable `target` that represents whether the user liked the song (encoded as 1) or not (encoded as 0). See the documentation of all the features [here](https://developer.spotify.com/documentation/web-api/reference/tracks/get-audio-features/).   Imports  Import libraries<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer, make_column_transformer

from sklearn.model_selection import (
    RandomizedSearchCV,
    cross_validate,
    train_test_split,
)
from sklearn.pipeline import make_pipeline 
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
    
from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.ensemble import RandomForestClassifier
from lightgbm.sklearn import LGBMClassifier
from xgboost import XGBClassifier
from catboost import CatBoostClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import VotingClassifier

from sklearn.metrics import (
    classification_report,
    roc_curve,
    RocCurveDisplay,
    roc_auc_score
)

import shap
from pylyrics import clean_text as ct
def mean_std_cross_val_scores(model, X_train, y_train, **kwargs):
    """
    Returns mean and std of cross validation

    Parameters
    ----------
    model :
        scikit-learn model
    X_train : numpy array or pandas DataFrame
        X in the training data
    y_train :
        y in the training data

    Returns
    ----------
        pandas Series with mean scores from cross_validation
    """

    scores = cross_validate(model, X_train, y_train, **kwargs)

    mean_scores = pd.DataFrame(scores).mean()
    std_scores = pd.DataFrame(scores).std()
    out_col = []

    for i in range(len(mean_scores)):
        out_col.append((f"%0.3f (+/- %0.3f)" % (mean_scores[i], std_scores[i])))

    return pd.Series(data=out_col, index=mean_scores.index)<jupyter_output><empty_output><jupyter_text>Reading the data CSVRead in the data CSV and store it as a pandas dataframe named `spotify_df`.<jupyter_code>spotify_df = pd.read_csv('data/spotify_df_processed.csv')#, index_col = 0 )
spotify_df.head(6)<jupyter_output><empty_output><jupyter_text>Data splittingSplit the data into train and test portions. Remove `song_title`, separate data to `X_train`, `y_train`, `X_test`, `y_test`.<jupyter_code>train_df, test_df = train_test_split(spotify_df, test_size=0.2, random_state=123)
X_train, y_train = train_df.drop(columns=["target"]), train_df["target"].astype('category')
X_test, y_test = test_df.drop(columns=["target"]), test_df["target"].astype('category')
# printing the number of observations for train and test sets
print('The number of observations for train set: ', train_df['target'].shape[0])
print('The number of observations for test set: ', test_df['target'].shape[0])<jupyter_output>The number of observations for train set:  3831
The number of observations for test set:  958<jupyter_text>I split 20% of the observations in the test data and 80% in the train data set. Overall the data set has about 4,000 observations, as it is not a very large data set, I preserved more portion for training.  Scoring metric<jupyter_code>print('The number of observations for positive targe: ', train_df["target"].sum())
print('The number of observations for negative targe: ', len(train_df["target"])-train_df["target"].sum())<jupyter_output>The number of observations for positive targe:  1993
The number of observations for negative targe:  1838<jupyter_text>Since it is a balanced data set and both positive and negative class are equally balanced, accuracy and precision-recall curves are selected as scoring metric.<jupyter_code>scoring_metric = ["accuracy", "roc_auc"]<jupyter_output><empty_output><jupyter_text>Preprocessing and transformations Here is different feature types and the transformations I will apply on each feature type.  | Transformation | Reason                                                                                                                                                                                                                                         | Features                                                                                                                                                                                                ||---------:|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|| OneHotEncoder |All of these features have fixed number of categories. `key` is circulated features like months. | `key`|| StandardScaler  | Numeric columns needs standardization | `acousticness`, `danceability`, `duration_ms`, `energy`, `instrumentalness`, `liveness`, `loudness`, `mode`, `speechiness`, `tempo`, `time_signature`, `valence`|| SimpleImputer, StandardScaler     | Analyze lyrics word using `paraphrase-distilroberta-base-v1`,  Numeric feature, processed in feature engineering  | `emb_sent_` || none     | The format are as required, represented as category features, processed in feature engineering | `genres_` || drop     | Free text column which has low correlation with target.| `song_title` || drop     | Replaced by artist genres `genres_`, do not want to be bounded by specific artist| `artist` |<jupyter_code>category_feats = ["key"]
none_feats = [col for col in X_train.columns if col.startswith('genres_')]
drop_feats = ['song_title', "artist", 'lyrics'] # Don't want to be bounded by artist


numeric_feats  = list(set(X_train.columns)
            - set(category_feats)
            - set(none_feats)
            - set(drop_feats)
        )

numeric_pipe = make_pipeline(SimpleImputer(missing_values=np.nan, strategy='mean'), StandardScaler())


preprocessor = make_column_transformer(
    ("drop", drop_feats),
    (numeric_pipe, numeric_feats),
    (OneHotEncoder(handle_unknown="ignore", sparse=False), category_feats),
    ("passthrough", none_feats)
)
preprocessor.fit(X_train, y_train);<jupyter_output><empty_output><jupyter_text>Baseline model<jupyter_code>results = {}
dummy = DummyClassifier()
baseline_pipe = make_pipeline(preprocessor, dummy)
results['dummy'] = mean_std_cross_val_scores(make_pipeline(preprocessor, dummy), X_train, y_train, 
                                             return_train_score=True, scoring=scoring_metric)
pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>Accuracy of Dummy classifier depends on class ratio. As we have a balanced data set, the accuracy of Dummy classifier is around 50%.  Linear models Model training - LogisticRegressionFirst, a linear model is used as a first real attempt. Hyperparameter tuning is also carried out for tuning to explore different values for the regularization hyperparameter. Cross-validation scores along with standard deviation and results summary is shown in below.<jupyter_code>#pipe logistic regression
pipe_logisticregression = make_pipeline(preprocessor,
                           LogisticRegression(max_iter=2000, 
                                             random_state=123))
#save in the results logistic regression score
results["LogisticReg"] = mean_std_cross_val_scores(pipe_logisticregression, 
                                                   X_train, 
                                                   y_train, 
                                                   return_train_score=True,
                                                   scoring = scoring_metric,
                                                   n_jobs=-1)
pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>hyperparameter optimization We will carry out hyperparameter optimization: C controls the regularization, and class_weight hyperparameter for tackling class imbalance.<jupyter_code>#parameters for logistic regression
param_dist_lg = {'logisticregression__C': np.linspace(2, 3, 6),
                 'logisticregression__class_weight': ['balanced', None]}

#randomized search to find the best parameters
random_search_lg = RandomizedSearchCV(
    pipe_logisticregression, 
    param_dist_lg,
    n_jobs=-1,
    return_train_score=True,
    scoring = scoring_metric,
    refit='accuracy',
    random_state=123
)
random_search_lg.fit(X_train, y_train)
print("Best parameter values are:", random_search_lg.best_params_)
print("Best cv score is:", random_search_lg.best_score_)
results['LogisticReg_opt'] = mean_std_cross_val_scores(random_search_lg, 
                                                       X_train, 
                                                       y_train, 
                                                       return_train_score=True,
                                                       scoring = scoring_metric,
                                                       n_jobs=-1) 
pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>We can see that with optimized hyperparameters, Logistic Regression is doing a bit better. However, it is obvious that we are dealing with overfitting (big gap between test and training scores and the training accuracy is almost 100%). The std is very small ranging in +- 0.01.<jupyter_code>print(
    classification_report(
        y_train,
        random_search_lg.predict_proba(X_train)[:, 1] > 0.5,
        target_names=["0", "1"],
    )
)
fpr, tpr, _ = roc_curve(
    y_train,
    random_search_lg.predict_proba(X_train)[:, 1],
    pos_label=random_search_lg.classes_[1],
)
print(
    "Area under the curve (AUC): {:.3f}".format(
        roc_auc_score(y_train, random_search_lg.predict_proba(X_train)[:, 1])
    )
)
roc_display = RocCurveDisplay(fpr=fpr, tpr=tpr).plot()<jupyter_output>Area under the curve (AUC): 0.997<jupyter_text>We have high score in AUC and classification report, which means our prediction model performance is good.  Model interpretation on Training set Most important features are listed in below.<jupyter_code>col_name_pp_all = [
    *numeric_feats,
    *random_search_lg.best_estimator_.named_steps["columntransformer"]
    .named_transformers_["onehotencoder"]
    .get_feature_names_out(),
    *none_feats
]

data = {
    "Importance": random_search_lg.best_estimator_.named_steps[
        "logisticregression"
    ].coef_[0],
}
feat_importance = pd.DataFrame(data=data, index=col_name_pp_all).sort_values(
    by="Importance", ascending=False
)

feat_importance_200 = pd.concat([feat_importance[0:100], feat_importance[-100:]])
feat_importance_200["rank"] = feat_importance_200["Importance"].rank(ascending=False)
feat_importance_200["side"] = np.where(
    feat_importance_200["Importance"] > 0, "pos", "neg"
)
feat_importance_200<jupyter_output><empty_output><jupyter_text>Most of the important features are determined by `genres`. Excluding all the features from artist, the other important features are `key`, `danceability` and `duration`. However, the ranking of importance are pretty low, starting from 423.<jupyter_code>feat_importance_reset = feat_importance.reset_index()
feat_importance_reset[feat_importance_reset['index'].str.match('^(?![genres])')]
import altair as alt

base = (
    alt.Chart(feat_importance_200.reset_index())
    .mark_bar()
    .encode(
        x=alt.X(
            "index",
        ),
        y="Importance:Q",
        color="side",
    )
    .properties(height=200, width=800, title="Top 200 important features")
)

brush = alt.selection_interval(encodings=["x"])
lower = (
    base.encode(
        x=alt.X(
            "index",
            axis=alt.Axis(labels=False, title="Features"),
            sort=alt.SortField(field="rank", order="ascending"),
        )
    )
    .properties(height=60, width=800, title="Drag the plot in below to zoom")
    .add_selection(brush)
)

upper = base.encode(
    alt.X(
        "index",
        scale=alt.Scale(domain=brush),
        axis=alt.Axis(title=""),
        sort=alt.SortField(field="rank", order="ascending"),
    )
)

upper & lower<jupyter_output><empty_output><jupyter_text>Top 200 most important features are based on genres.   Non-linear models Model training - RandomForestClassifier, XGBClassifier, LGBMClassifier, CatBoostClassifierSecond, four non-linear model are trained aside from the linear model above.  After that, feature selection and hyperparameter tuning will carried out in later stage. Cross-validation scores along with standard deviation and results summary is shown in below.<jupyter_code># Random Forest pipe
pipe_rf = make_pipeline(preprocessor, RandomForestClassifier(random_state=123))

# XGBoost pipe
pipe_xgb = make_pipeline(
    preprocessor,
    XGBClassifier(
        random_state=123, eval_metric="logloss", verbosity=0, use_label_encoder=False
    ),
)

# LGBM Classifier pipe
pipe_lgbm = make_pipeline(preprocessor, LGBMClassifier(random_state=123))

# Catboost pipe
pipe_catb = make_pipeline(preprocessor, CatBoostClassifier(verbose=0, random_state=123))

models = {
    "RandomForest": pipe_rf,
    "XGBoost": pipe_xgb,
    "LGBM": pipe_lgbm,
    "Cat_Boost": pipe_catb,
}

# summarize mean cv scores in result_non_linear
for (name, model) in models.items():
    results[name] = mean_std_cross_val_scores(
        model, X_train, y_train, return_train_score=True, scoring=scoring_metric
    )
pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>All the non-linear models are overfitting as all the training scores are close to 1. Compared to most of the models, `LGBM` and `Cat_Boost` are two of the most balanced model in performance, it has the highest accuracy and comparatively less overfitting (the gap between train and validation scores is smaller).Stability of scores is more or less stable, with standard deviation ranging in around 0.06. The fit time of `LGBM` is much shorter than `Cat_Boost`, which is important for the model application, therefore, `LGBM` is by far the most suitable model.  Model interpretation on Training setSince `feature_importance_` is not supported on the models, `shap` is used to examine the most important features of our best non-linear models.<jupyter_code>pipe_lgbm.fit(X_train, y_train);
X_train_enc = pd.DataFrame(
    data=preprocessor.transform(X_train),
    columns=col_name_pp_all,
    index=X_train.index,
)
lgbm_explainer = shap.TreeExplainer(pipe_lgbm.named_steps["lgbmclassifier"])
train_lgbm_shap_values = lgbm_explainer.shap_values(X_train_enc)
shap.summary_plot(train_lgbm_shap_values, X_train_enc, plot_type="bar")<jupyter_output>LightGBM binary classifier with TreeExplainer shap values output has changed to a list of ndarray<jupyter_text>We examined the most important features of `LGBMClassifier` with SHAP methods. Apart from genre, the results suggests that the most important features are `instrumentalness`, `danceability` and `loudness`. Comparing to the linear model, which suggested all `genre_` is the most important features. non-linear model has more diversity.  Model AveragingAll the models are overfitted. To ease the fundamental trade off, `Averaging` is attempted to see if a higher validation score and shorter fit time can be achieved.<jupyter_code>avg_classifiers = {
    "random forest": pipe_rf,
    "XGBoost": pipe_xgb,
    "LightGBM": pipe_lgbm
}

averaging_model = VotingClassifier(
    list(avg_classifiers.items()), voting="soft"
) 

results['averaging'] = mean_std_cross_val_scores(
        averaging_model, X_train, y_train, return_train_score=True, scoring=scoring_metric#, cv=2
    )
pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>According to the result, averaging the models neither help the overfitting issue nor scoring. However, the fit time is slower than `LGBMClassifier`. `LGBMClassifier` is still the best model.  Results on the test set The best performing model `LGBMClassifier` is used on the test data and report test scores. Summary is shown as below. Furthermore, some test predictions and corresponding explanation with SHAP force plots are drawn for further study.<jupyter_code>#predictions = pipe_catb.predict(X_test)
print("Test Set accuracy: ", round(pipe_lgbm.score(X_test, y_test), 3))<jupyter_output>Test Set accuracy:  0.954<jupyter_text>The test set score was similar to the validation score. Therefore, the prediction is promising.   Interpretation and feature importances on Test set<jupyter_code># Encoding X_test for SHAP force plot
X_test_enc = pd.DataFrame(
    data=preprocessor.transform(X_test),
    columns=col_name_pp_all,
    index=X_test.index,
)
X_test_enc.shape
# Create an explainer for X_test_enc
test_lgbm_shap_values = lgbm_explainer.shap_values(X_test_enc)
shap.summary_plot(test_lgbm_shap_values, X_test_enc, plot_type="bar")
y_test.iloc[8]
# Force Plot for prediction of 4th row (no_default)
shap.force_plot(
    lgbm_explainer.expected_value[1],
    test_lgbm_shap_values[1][8,:],
    X_test_enc.iloc[8,:],
    matplotlib=True,
)<jupyter_output><empty_output><jupyter_text>- As seen from the plot, the raw score is much smaller than the base value, which predict accurately as **Like** (class of '1').  - `danceability` and `genres_dance_pop` are pushing the prediction towards a higher score, and `tempo` is pushing the prediction towards a lower score.<jupyter_code>y_test.iloc[10]
# Force Plot for prediction of 2nd row (no_default)
shap.force_plot(
    lgbm_explainer.expected_value[1],
    test_lgbm_shap_values[1][4,:],
    X_test_enc.iloc[4,:],
    matplotlib=True,
)<jupyter_output><empty_output><jupyter_text>- As seen from the plot, the raw score is much smaller than the base value, which predict accurately as **Unlike** (class of '0').  - `genres_dance_pop` is pushing the prediction towards a higher score, and `genres_modern rock` and `instrumentalness` are pushing the prediction towards a lower score.  Summary of machine learning results<jupyter_code>pd.DataFrame(results)
pd.DataFrame(results).to_csv('data/model_results.csv',index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11004526837680535, 'fraction_numerical': 0.009107566285837466, 'mean_word_length': 3.4490529848957086, 'pattern_counts': {'":': 8, '<': 72, '<?xml version=': 0, '>': 74, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '179587', 'license': 'MIT', 'n_tokens_mistral': 5430, 'n_tokens_neox': 5207, 'n_words': 1516, 'path': 'spotify_user_behaviour_predictor/_build/jupyter_execute/6_Machine_Learning.ipynb', 'repo_name': 'MacyChan/spotify-user-behaviour-predictor', 'size': '18556'}