Predicting abalone age and dealing with multicollinearity.

!pip install pingouin &> /dev/null
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sp
from sklearn import preprocessing
from sklearn.linear_model import LinearRegression
from sklearn import metrics
import pingouin as pg
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split

abalone = pd.read_csv('./data/abalone.csv')
display(abalone.head())
sns.pairplot(abalone, hue = 'sex');

# Violin Plot of sex distributed by age
sns.violinplot(x='sex', y='age', data=abalone)
plt.show()

sexlista = ['Males','Females','Infants']

for i,k in zip(sexlista,['M','F','I']):
    dfinn = abalone[abalone.sex==k]
    print(f'Distribution of weight by age for {i}')
    df_age = dfinn[['sex','age','whole_wt']].groupby('age').mean()
    display(df_age.sort_values(by='whole_wt', ascending = False).head(3))
    df_age.plot(kind = 'bar',title = f'Distribution of weight by age for {i}')
    plt.show()

# Print regression result

def regression_metrics(y_true, y_pred):
    # Regression metrics
    ev = metrics.explained_variance_score(y_true, y_pred)
    #msle = metrics.mean_squared_log_error(y_true, y_pred) 
    #penalizes an under-predicted estimate greater than an over-predicted estimate.
    mae = metrics.mean_absolute_error(y_true, y_pred)
    mdae = metrics.median_absolute_error(y_true, y_pred)
    mse = metrics.mean_squared_error(y_true, y_pred) 
    r2 = metrics.r2_score(y_true, y_pred)
    rmse = np.sqrt(mse)
    acc = (rmse/mae-1) #accuracy: For 0 (better) to sqrt(n)-1
    #acd = results.rsquared_adj # adjusted coefficient of determination
    #cod = results.rsquared  # coefficient of determination
    return ev, mae, mdae, mse, r2, rmse, acc

def regression_results(y_true, y_pred):
    ev, mae, mdae, mse, r2, rmse, acc = regression_metrics(y_true, y_pred)
    print('Explained Variance: ', round(ev, 6))    
    #print('Mean Square Log Error: ', round(msle, 6))
    # print('Mean Absolute Error: ', round(mae, 6))
    print('Median Absolute Error: ', round(mdae, 6))
    # print('Mean Squared Error: ', round(mse, 8))
    # print('Root Mean Squared Error: ', round(rmse, 6))
    print('Accuracy Acc(RMSE/MAE - 1): ', round(acc, 6))
    #print('Adjusted Coefficient of Determination (adjusted R2): ', round(acd, 6))
    #print('Coefficient of Determination: ', round(cod, 6))
        
    
# # Create dataframe with comparative results (actual vs predicted)    
def creates_df_compare(y_test,y_pred):
    df_compare = pd.DataFrame({'actual': y_test, 'pred': y_pred})
    df_compare['error'] = df_compare['actual'] - df_compare['pred']
    df_compare['residue'] = np.abs(df_compare['actual'] - df_compare['pred'])
    df_compare['residue_sqr'] = df_compare['residue']**2
    avg_actual = df_compare['actual'].mean()
    df_compare['var'] = (df_compare['actual'] - avg_actual)**2
    return df_compare


# Average Relative Error Chart - Real vs Predicted
def real_vs_predicted(yhat,y):
    ev, mae, mdae, mse, r2, rmse, acc = regression_metrics(y, yhat)
    x = list(range(len(y)))
    print("Mean Absolut Error:",mae)
    print("Mean Squared Error:", mse)
    print("Root Mean Squared Error:", rmse)
    print("R-Squared:", r2)
    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(16,8))
    fig.suptitle('Original vs Predicted')

    ax1.plot(x, y, color="blue", label="original")
    ax1.set_ylabel('Original vs Predicted')
    ax1.plot(x, yhat, color="red", label="original")

    ax2.plot(x, y, color="blue", label="original")
    ax2.set_ylabel('Original')

    ax3.plot(x, yhat, color="red", label="predicted",)
    ax3.set_ylabel('Predicted')
    plt.show()
    

cols = ['sex', 'length', 'diameter', 'height', 'whole_wt', 'shucked_wt','viscera_wt', 'shell_wt']
numcols = cols[1:]

df_norm = abalone[cols].copy()

# To eliminate possible problems of scale difference between the variables, we will apply the standardization.
scaler = preprocessing.StandardScaler()
# Apply standardization to numeric columns
df_norm[numcols] = scaler.fit_transform(df_norm[numcols])
# drop_first=True to prevent the linear regression matrix from being invertible (non-singular)
# we must have one column less than the categories
x_scaled = pd.get_dummies(df_norm, drop_first=True)

y = abalone['age']

regressor = LinearRegression()
results = regressor.fit(x_scaled,y)
y_pred = results.predict(x_scaled)

regression_results(y, y_pred)

display(creates_df_compare(y, y_pred).round(4))

real_vs_predicted(y_pred, y)

fig, ax = plt.subplots(figsize=(10,5))
dfcorr = abalone.corr()
sns.heatmap(dfcorr, annot = True, cmap = 'PuBu');

fig, ax = plt.subplots(figsize=(10,5))
sns.heatmap(abalone.pcorr(), annot = True, cmap = 'PuBu');

coeff_df = pd.DataFrame(results.coef_, x_scaled.columns, columns=['coefficient'])

coeff_df.reset_index(inplace=True)

coeff_df = coeff_df.rename(columns = {'index':'feature'})

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

ax = sns.barplot(x = 'coefficient', y =  'feature', 
                 data = coeff_df.sort_values(by='coefficient', 
                                             ascending = False))

coeff_df['coeff_abs'] = np.abs(coeff_df['coefficient'])

coeff_df['coeff_perc'] = coeff_df['coeff_abs'] / np.sum(coeff_df['coeff_abs'])

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

ax = sns.barplot(x = 'coeff_perc', y =  'feature', 
                 data = coeff_df.sort_values(by='coeff_abs', 
                                             ascending = False))

# split train and test
X_train, X_test, y_train, y_test = train_test_split( x_scaled, y, test_size=0.25, random_state=42)

pca = PCA(n_components = x_scaled.shape[1])
pca_data = pca.fit_transform(X_train)

percent_var_explained = pca.explained_variance_/(np.sum(pca.explained_variance_))
cumm_var_explained = np.cumsum(percent_var_explained)

plt.plot(cumm_var_explained)
plt.grid()
plt.xlabel("n_components")
plt.ylabel("% variance explained")
plt.show()