Predicting Recipe Website Traffic with Machine Learning

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import ks_2samp, pointbiserialr

df = pd.read_csv('recipe_site_traffic_2212.csv')
df.info()

df.head()

def validate_helper(df, col: str, unique_list:bool=True):
    """Prints basic info to help validate a column."""
    print(f'{col}:\nmissing values: {df[col].isna().sum()}\nunique values: {df[col].nunique()}')
    if unique_list:
        print(f'unique value list: {df[col].unique()}')
    print(df[col].describe())

# validate recipe
validate_helper(df, 'recipe', unique_list=False)
recipe_min = df['recipe'].min()
recipe_max = df['recipe'].max()
recipe_expected_sum = (recipe_max - recipe_min + 1) * (recipe_min + recipe_max) / 2
recipe_actual_sum = df['recipe'].sum()
if recipe_expected_sum == recipe_actual_sum:
    print("The 'recipe' column is a unique identifier.")
else:
    print("The 'recipe' column is NOT a unique identifier.")

# validate 'calories'
validate_helper(df, 'calories', unique_list=False)

# look at n/a values
df[df['calories'].isna()]

# validate 'carbohydrate'
validate_helper(df, 'carbohydrate', unique_list=False)

# validate 'sugar'
validate_helper(df, 'sugar', unique_list=False)

# validate 'protein'
validate_helper(df, 'protein', unique_list=False)

# validate 'category'
validate_helper(df, 'category', unique_list=True)

# query 'chicken_breast' and 'chicken'
count_before = df[df['category'].isin(['Chicken Breast', 'Chicken'])].shape[0]
print("count before modification:", count_before)

# clean 'chicken_breast'
df.loc[df['category'] == 'Chicken Breast', 'category'] = 'Chicken'
count_after = df[df['category'] == 'Chicken'].shape[0]
print("count after modification:", count_after)

# update dtype to categorical. 
df['category'] = df['category'].astype('category')
df['category'].dtype

# validate 'servings'
validate_helper(df, 'servings', unique_list=True)

# servings has object datatype and 6 unique values, suggesting one or more are formatted incorrectly
df['servings'].unique

if pd.api.types.is_string_dtype(df['servings']):
    df[df['servings'].isin(['4 as a snack', '6 as a snack'])]

    # clean 'as a snack'
    df['servings'] = df['servings'].str.replace(' as a snack', '', regex=False)

    # change dtype
    df['servings'] = df['servings'].astype(int)

# validate high traffic
validate_helper(df,'high_traffic', unique_list=True)

# binary encode, replacing 'high' as 1, and missing values as 0. 
df['high_traffic'] = df['high_traffic'].replace('High', 1).fillna(0)

# verify no more missing values
print(f"missing values remaining: {df['high_traffic'].isna().sum()}")

# change dtype to int, to prepare for machine learning
df['high_traffic'] = df['high_traffic'].astype(int)

# check unique values
print(f"new unique values: {df['high_traffic'].unique()}")

df.describe()

# drop missing values
cleaned_df = df.dropna(subset=['calories','carbohydrate','sugar','protein'])
cleaned_df.shape

# create nutrional totals
nutrient_vals = ['calories','carbohydrate','sugar','protein']
"""
nutrient_vals_totals = ['calories_total','carbohydrate_total','sugar_total','protein_total']

for nut in nutrient_vals:
    cleaned_df[nut + '_total'] = cleaned_df[nut] * cleaned_df['servings']
"""
cleaned_df.head()

## investigate outliers, by category
def calculate_z_scores_by_category(group):
    """Calculates Z-score within each category"""
    return (group - group.mean()) / group.std()

# apply scaling within each category while preserving indexes
category_scaled_df = cleaned_df.groupby('category')[['calories', 'carbohydrate', 'sugar', 'protein']].transform(calculate_z_scores_by_category)

# add categorical data
category_scaled_df['category'] = cleaned_df['category']

# transform to long format for cat plot. 
melted_df = category_scaled_df.melt(id_vars=['category'], value_vars=['calories', 'carbohydrate', 'sugar', 'protein'], 
                                    var_name='Nutrient', value_name='Z-score')

# plot
g = sns.catplot(x='category', y='Z-score', hue='Nutrient', kind='box', data=melted_df, height=5, aspect=2)
g.set_xticklabels(rotation=0)
g.set_axis_labels("Category", "Z-score")
plt.title('Nutritional Value Z-Scores for Each Category')
plt.show()

## remove outliers based on calculated bounds per category
## updated this to keep outlier removal contained to filtered_df and not updating cleaned_df

# empty dataframe for filtered data
filtered_df = pd.DataFrame()

# loop through each category
for cat in cleaned_df['category'].unique():
    category_data = cleaned_df[cleaned_df['category'] == cat]

    for col in nutrient_vals:
        # calculate Q1, Q3, and IQR
        Q1 = category_data[col].quantile(0.25)
        Q3 = category_data[col].quantile(0.75)
        IQR = Q3 - Q1

        # filter outliers based on bounds
        lower_bound = Q1 - 1.5 * IQR
        upper_bound = Q3 + 1.5 * IQR
        category_data = category_data[(category_data[col] >= lower_bound) & (category_data[col] <= upper_bound)]
        
    # append current category to filter_df
    filtered_df = pd.concat([filtered_df, category_data], ignore_index=True)
    
# replace cleaned_df with filtered_df
#cleaned_df = filtered_df

# show remaining boxplots
#cleaned_df.shape

# print proportion of data set exceeds 1.5x IQR bounds
outlier_proportion = (cleaned_df.shape[0] - filtered_df.shape[0]) / cleaned_df.shape[0]
print(outlier_proportion)

## updated this so bottom graph uses filtered_df instead of cleaned_df, since outliers are not being removed. 

# apply z-scores, per category
original_grouped_z_scores = df.groupby('category')[['calories', 'carbohydrate', 'sugar', 'protein']].transform(calculate_z_scores_by_category)
cleaned_grouped_z_scores = filtered_df.groupby('category')[['calories', 'carbohydrate', 'sugar', 'protein']].transform(calculate_z_scores_by_category)

# add 'category' column
original_grouped_z_scores['category'] = df['category']
cleaned_grouped_z_scores['category'] = filtered_df['category']

# melt the dataframes for plotting
original_melted = original_grouped_z_scores.melt(id_vars=['category'], value_vars=['calories', 'carbohydrate', 'sugar', 'protein'],
                                                 var_name='Nutrient', value_name='Z-score')
cleaned_melted = cleaned_grouped_z_scores.melt(id_vars=['category'], value_vars=['calories', 'carbohydrate', 'sugar', 'protein'],
                                               var_name='Nutrient', value_name='Z-score')
# concat data
data = pd.concat([original_melted.assign(dataset='Original'), cleaned_melted.assign(dataset='Outliers Filtered')], ignore_index=True)

# map boxplots to data
g = sns.FacetGrid(data, row="dataset", height=5, aspect=3, sharex=False) 
g.map_dataframe(sns.boxplot, x='category', y='Z-score', hue='Nutrient', palette='colorblind')
g.set_xticklabels(rotation=0)
g.set_axis_labels("Category", "Z-score")
g.add_legend()

g.fig.suptitle('Nutritional Value Z-Scores for Each Category, Before/After Identifying Outliers', fontsize=16, y=1.05)

plt.show()

# determine if target variable distribution differs signficantly due to removed outliers
ks_stat, ks_p = ks_2samp(df['high_traffic'], cleaned_df['high_traffic'])
if ks_p < 0.05:
    sig = 'The distributions are signficantly different'
else:
    sig = 'No statistical difference in the distributions. '
print(f'KS Statistic (high_traffic): {ks_stat}, P-value: {ks_p}. {sig}')

# overal count plot
sns.countplot(x='high_traffic', hue='high_traffic', data=cleaned_df, palette='colorblind', dodge=False)
plt.title('High Traffic Distribution, Overall')
plt.xticks(ticks=[0,1],labels=['Normal', 'High'], rotation=45)
plt.xlabel('High Traffic')
plt.ylabel('Count')
plt.legend().remove() 
plt.show()

plt.figure(figsize=(12, 6))

# prepare data
category_counts = pd.crosstab(cleaned_df['category'], cleaned_df['high_traffic'])

# high traffic vs category
ax1 = plt.subplot(1, 2, 1)
category_counts.plot(kind='bar', color=sns.color_palette("colorblind"), stacked=False, ax=ax1)
plt.title('High Traffic Distribution, Per Category')
plt.xlabel('Category')
plt.ylabel('Count')
plt.xticks(rotation=45)
plt.legend(title='Traffic Status', labels=['Normal', 'High'])

category_counts = pd.crosstab(cleaned_df['servings'], cleaned_df['high_traffic'])

# high traffic vs servings
ax2 = plt.subplot(1, 2, 2)
category_counts.plot(kind='bar', color=sns.color_palette("colorblind"), stacked=False, ax=ax2)
plt.title('High Traffic Distribution, Per Serving')
plt.xlabel('Servings')
plt.ylabel('Count')
plt.xticks(rotation=0)
plt.legend(title='Traffic Status', labels=['Normal', 'High'])


plt.tight_layout()
plt.show()

plt.figure(figsize=(12, 18))

ax1 = plt.subplot(4, 2, 1)
sns.histplot(data=cleaned_df, x='calories', kde=True, ax=ax1)
plt.title('Calories, Histplot')

ax2 = plt.subplot(4, 2, 2)
sns.boxplot(data=cleaned_df, y='calories', ax=ax2)
plt.title('Calories, Boxplot')

ax3 = plt.subplot(4, 2, 3)
sns.histplot(data=cleaned_df, x='carbohydrate', kde=True, ax=ax3)
plt.title('Carbohydrate, Histplot')

ax4 = plt.subplot(4, 2, 4)
sns.boxplot(data=cleaned_df, y='carbohydrate', ax=ax4)
plt.title('Carbohydrate, Boxplot')

ax5 = plt.subplot(4, 2, 5)
sns.histplot(data=cleaned_df, x='sugar', kde=True, ax=ax5)
plt.title('Sugar, Histplot')

ax6 = plt.subplot(4, 2, 6)
sns.boxplot(data=cleaned_df, y='sugar', ax=ax6)
plt.title('Sugar, Boxplot')

ax7 = plt.subplot(4, 2, 7)
sns.histplot(data=cleaned_df, x='protein', kde=True, ax=ax7)
plt.title('Protein, Histplot')

ax8 = plt.subplot(4, 2, 8)
sns.boxplot(data=cleaned_df, y='protein', ax=ax8)
plt.title('Protein, Boxplot')

plt.tight_layout()
plt.show()

# numeric heatmap
heatmap_vals = ['high_traffic'] + nutrient_vals
numeric_vals = nutrient_vals
numeric = cleaned_df[heatmap_vals]
sns.heatmap(numeric.corr(),annot=True).set(title='Correlation Heatmap Between Numeric Variables and High Traffic')
plt.show()   

# compile correlation data into dataframe
correlation_data = []

# linear correlations with 'high_traffic'
for col in numeric_vals:
    correlation, p_value = pointbiserialr(cleaned_df['high_traffic'], cleaned_df[col])
    correlation_data.append({'Variable': col, 'Type': 'Original', 'Correlation': correlation, 'P-Value': p_value})
    print(f"{col} to high_traffic correlation: {correlation:.4f} (P-value: {p_value:.4f})")

# linear correlations with 'high_traffic', per category
for cat in cleaned_df['category'].unique():
    category_df = cleaned_df[cleaned_df['category'] == cat]
    for col in numeric_vals:
        correlation, p_value = pointbiserialr(category_df['high_traffic'], category_df[col])
        print(f"{cat},{col} to high_traffic correlation: {correlation:.4f} (P-value: {p_value:.4f})")

# apply log transformation
cleaned_df_log = cleaned_df.copy()
for col in numeric_vals:
    cleaned_df_log[col] = np.log(cleaned_df_log[col] + 1)

# log correlation with  'high_traffic'
for col in numeric_vals:
    correlation, p_value = pointbiserialr(cleaned_df_log['high_traffic'], cleaned_df_log[col])
    correlation_data.append({'Variable': col, 'Type': 'Log-Transformed', 'Correlation': correlation, 'P-Value': p_value})
    print(f"Log-transformed {col} to high_traffic correlation: {correlation:.4f} (P-value: {p_value:.4f})")
# log correlation with 'high_traffic', per category
for cat in cleaned_df_log['category'].unique():
    category_df = cleaned_df_log[cleaned_df_log['category'] == cat]
    for col in numeric_vals:
        correlation, p_value = pointbiserialr(category_df['high_traffic'], category_df[col])
        print(f"{cat}, Log-transformed {col} to high_traffic correlation: {correlation:.4f} (P-value: {p_value:.4f})")

## heatmaps for correlation and p-values
# convert correlation data to dataframe
correlation_df = pd.DataFrame(correlation_data)

# plot
fig, ax = plt.subplots(1, 2, figsize=(16, 6))

# correlation
correlation_pivot = correlation_df.pivot("Variable", "Type", "Correlation")
sns.heatmap(correlation_pivot, annot=True, fmt=".3f", cmap='vlag', center=0, vmin=-1, vmax=1, ax=ax[0])
ax[0].set_title('Correlation of Numeric Values with High Traffic')

# p-value
correlation_pivot = correlation_df.pivot("Variable", "Type", "P-Value")
sns.heatmap(correlation_pivot, annot=True, fmt=".4f", cmap='vlag_r', center=0.5, vmin=0, vmax=1, ax=ax[1])
ax[1].set_title('P-Values of Numeric Values with High Traffic')

plt.show()

def plot_logistic_fit(data, x, y, axn, transformed=False, ):
    sns.regplot(x=x, y=y, data=data, logistic=True, ci=None, scatter_kws={'alpha':0.5}, ax=ax[axn])
    title = 'Logistic Regression Fit of ' + ('Log-Transformed ' if transformed else '') + x
    ax[axn].set_title(title)
    plt.xlabel('Log-Transformed ' + x if transformed else x)
    plt.ylabel(y)

fig, ax = plt.subplots(1, 2, figsize=(16, 6))
plot_logistic_fit(cleaned_df, 'protein', 'high_traffic', axn=0)
plot_logistic_fit(cleaned_df_log, 'protein', 'high_traffic', axn=1,transformed=True)
plt.show()

fig, ax = plt.subplots(1, 2, figsize=(14, 7))

# pie chart for 'category'
cleaned_df['category'].value_counts().plot.pie(autopct='%1.1f%%', ax=ax[0])
ax[0].set_title('Category Proportions')
ax[0].set_ylabel('')


# pie chart for 'serving'
cleaned_df['servings'].value_counts().plot.pie(autopct='%1.1f%%',ax=ax[1])
ax[1].set_title('Servings Proportions')
ax[1].set_ylabel('')

plt.tight_layout()
plt.show()

# dummy categories
category_dummies = pd.get_dummies(cleaned_df['category'])

categories = []
correlations = []
p_values = []

# calc corr and p-value
for column in category_dummies:
    corr, p_val = pointbiserialr(category_dummies[column], cleaned_df['high_traffic'])
    categories.append(column)
    correlations.append(corr)
    p_values.append(p_val)

# store results in idf
results_df_cleaned = pd.DataFrame({
    'Category': categories,
    'Correlation': correlations,
    'P-Value': p_values
})

# correlation bar chart with p-values
fig, ax = plt.subplots(figsize=(12, 6))
bars = sns.barplot(x='Category', y='Correlation', data=results_df_cleaned, ax=ax)
ax.set_title('Correlation of Categories with High Traffic')
ax.set_ylabel('Correlation Coefficient')
ax.set_ylim([-1, 1])
ax.set_xticklabels(results_df_cleaned['Category'], rotation=45, ha="right")

"""
# p-value anotations
for bar, p_value in zip(bars.patches, results_df_cleaned['P-Value']):
    y = bar.get_height()
    x = bar.get_x() + bar.get_width() / 2
    ax.text(x, y, f'p={p_value:.9f}', ha='center', va='bottom' if y >= 0 else 'top', color='black', fontsize=9)
"""
plt.show()

# dummy categories
category_dummies = pd.get_dummies(cleaned_df['category'])
cat_corr_df = pd.concat([cleaned_df, category_dummies], axis=1)
categories = []
correlations = []
p_values = []

# calculate correlation and p-values for each category
for column in category_dummies.columns:
    correlation, p_value = pointbiserialr(cat_corr_df[column], cat_corr_df['high_traffic'])
    categories.append(column)
    correlations.append(correlation)
    p_values.append(p_value)

# create dataframe
correlation_df = pd.DataFrame({
    'Category': categories,
    'Correlation': correlations,
    'P-Value': p_values
})

# plot
fig, ax = plt.subplots(1, 2, figsize=(16, 6))

# corr
correlation_matrix = correlation_df.pivot_table(index='Category', values='Correlation', aggfunc='sum')
sns.heatmap(correlation_matrix, annot=True, fmt=".3f", cmap='vlag', ax=ax[0])
ax[0].set_title('Correlation of Categories with High Traffic')

# p-value
p_values_matrix = correlation_df.pivot_table(index='Category', values='P-Value', aggfunc='sum')
sns.heatmap(p_values_matrix, annot=True, cmap='vlag_r', center=0.5, vmin=0, vmax=1, ax=ax[1])
ax[1].set_title('P-Values of Categories with High Traffic')

plt.tight_layout()
plt.show()

from sklearn.model_selection import train_test_split, GridSearchCV, StratifiedKFold
from sklearn.preprocessing import OneHotEncoder, FunctionTransformer, MinMaxScaler
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, precision_score

feature_cols = X.columns
resultdict = {}

model = logreg2

if len(model.coef_.shape) > 1:
    coefs = model.coef_[0]
else:
    coefs = model.coef_    

for i, col_name in enumerate(feature_cols):
    resultdict[col_name] = coefs[i]
    
plt.figure(figsize=(10, 8))
plt.bar(resultdict.keys(), resultdict.values())
plt.xticks(rotation=90)
plt.title('Feature Importance in Logistic Regression Model')
plt.ylabel('Coefficient Value')
plt.show()

feature_cols = X.columns
resultdict = {}

model = rfc2

importances = model.feature_importances_

for i, col_name in enumerate(feature_cols):
    resultdict[col_name] = importances[i]
    
plt.figure(figsize=(10, 8))
plt.bar(resultdict.keys(), resultdict.values())
plt.xticks(rotation=90)
plt.title('Feature Importance in Random Forest Model')
plt.ylabel('Coefficient Value')
plt.show()

# removed to focus on 2 models. 
# LinearSVC
"""
# define hyperparameter grid
grid = {
    'C': [0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 1],
    'penalty': ["l1", "l2"],
    'loss': ['hinge','squared_hinge']
}

# grid cross validate and fit
svm = LinearSVC(random_state=RAND_STATE)
svm_cv = GridSearchCV(svm, grid, cv=cv, scoring=SCORING_METHOD, verbose=1)
svm_cv.fit(X_train, y_train)

# display results
print(f'Best Score: {svm_cv.best_score_}')
print(f'Best Hyperparameters: {svm_cv.best_params_}')
print(f'Std deviation of CV scores for the best hyperparameters: {svm_cv.cv_results_["std_test_score"][svm_cv.best_index_]}')

# unpack best_params to create model
svm2 = LinearSVC(**svm_cv.best_params_)
svm2.fit(X_train, y_train)

# evaluate
y_pred_svm = svm2.predict(X_test)

print('Classification report:\n', classification_report(y_test, y_pred_svm))
print('Confusion matrix:\n', confusion_matrix(y_test, y_pred_svm))
"""

models = list(results.keys())
precision_scores = [results[model]['Precision'] for model in models]
accuracy_scores = [results[model]['Accuracy'] for model in models]


fig, axs = plt.subplots(2, 1, figsize=(10, 8), sharex=False)

# accuracy
sns.barplot(x=models, y=accuracy_scores, palette='colorblind', ax=axs[0])
axs[0].axhline(y=0.598, color='r', linestyle='--', label='Current Business Practice (60%)')
axs[0].set_title('Goal 1 KPI: Accuracy Scores')
axs[0].set_ylabel('Accuracy')
axs[0].legend()
for i, score in enumerate(accuracy_scores):
    axs[0].text(i, score - 0.100, f'{score:.2f}', color='white', ha='center', fontsize= 14)

# precision
sns.barplot(x=models, y=precision_scores, palette='colorblind', ax=axs[1])
axs[1].axhline(y=0.8, color='r', linestyle='--', label='KPI Goal (80%)')
axs[1].set_title('Goal 2 KPI: Precision Scores')
axs[1].set_ylabel('Precision')
axs[1].legend()
for i, score in enumerate(precision_scores):
    axs[1].text(i, score - 0.100, f'{score:.2f}', color='white', ha='center', fontsize= 14)

plt.tight_layout()
plt.show()