<jupyter_start><jupyter_text>Used vehicles price prediction - A Machine Learning Project by <NAME>*first full ML project in Python*  Abstract The purpose of this project is to explore and understand the process of predicting the price of used vehicles in the U.S.A, the most important characteristics that influence the sale price, and finally to build a supervised machine learning model to generate predictions. The best-performing model resulting from this analysis is a **feed-forward neural network** supported by Tensorflow / Keras.The data for this project was obtained from Kaggle datasets, which was scraped / queried from Craiglist.This project originated from a personal inquire: As I started my MBA journey in the U.S.A, how much would it cost to buy a car for my personal use? Can I budget for it? This question could trivially be answered for any-one specific manufacturer or model, but I wanted to build on my Python and ML skills and finally, I chose to transform this curiosity into a full-fledged project.As such, the **goals** for these project were defined as:* First, and most important goal is: To create a regression model capable for predicting a used vehicle price.* Second, and depending on the the results of the models, I  will attempt to obtain information what are the most important attributes (model weights).Exploratory data analysis revealed early the complexity of the task: There is a substantial amount of manufacturers and even more models with multiple layers of variation and complexity. More importantly, analysis of the model runs and the residuals showed the instability of some of the models as all results were heavily influenced by the presence of substantially large outlier errors. Considering all that was said, the *Mean Average Error* is my chosen KPI for evaluating residual, or in other words **how successful is this model in predicting the price of a used car**. However to CHOOSE a model the preferred metric was the *Root Mean Squared Error* as it is much better at capturing extreme deviations and outliers, therefore signaling a far more stable model.The models used in this analysis were:* **Scikit-Learn supported models**:    * Lasso Regression    * Standard Linear Regression    * Random Forest Regressor    * Gradient Boosted Regressor* **Tensorflow-Keras supported model**:    * Feed-forward Neural Network with 3 hidden layers**SCOPE**: * For time reasons I will **exclude** any kind of **image analysis**. Vehicle images (could be extracted from each post) will be ignored.* Manufacturers with less than 100 samples in the dataset will be excluded* Car models with less than 30 samples in the datasets will be excluded* In-depth NLP tasks are not in scope of this analysis, but some simple tasks (bag-of-words, etc) may be used for feature extraction.* Unique data points such as VIN, which could be used in history analysis, will also be excluded.* These exercise will attempt to create a single model to predict ALL used vehicles regardless of manufacturer or model. Under some circumstances (Tesla's for example) vehicles exhibit substantially different behaviours on price based on their brand and a specific Manufacturer-Regression-Model could be a much better approach. However, for simplicity's sake, **I'll do a "one for all" model**.**Finally** using Root Mean Squared Error (**RMSE**) as our main KPI, the best performing model was a **Tensorflow/Keras Feed-Forward Neural Network** achieving a RMSE of 2,894 USD and a Mean Average Error (**MAE**) of 1,775 USD.NOTE: Irregularities and inflationary pressures affected the price of used vehicles throughout the year 2021. The effects of this is neither reflected nor considered in this exercise.This analysis was carried out in May-June 2021  Table of Contents1. [Initialization](Initialization)    1. [Imports](Imports)    1. [Functions](Functions)1. [Data_exploration](Data_exploration)    1. [Handling_missing_data](Handling_missing_data)    2. [Exploratory_data_analysis](Exploratory_data_analysis)1. [Model_training](Model_training)    1. [Linear_models](Linear_models)    2. [Ensemble_regressors](Ensemble_regressors)    3. [Keras_model](Keras_model)1. [Model_evaluation](Model_evaluation)  Initialization  Imports<jupyter_code># standard imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import seaborn as sns

# scikit-learn imports
## scikit tools
from sklearn.feature_selection import mutual_info_regression, f_regression
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error

## joblib imports for saving / loading objects
from joblib import dump, load

## scikit models 
from sklearn.linear_model import LassoCV, LinearRegression
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

# tensorflow imports
from tensorflow.keras.models import load_model, save_model, Sequential
from tensorflow.keras.layers import Dense, Dropout, Input
from tensorflow.keras.callbacks import EarlyStopping

# setting style
plt.style.use("seaborn-colorblind")<jupyter_output><empty_output><jupyter_text>The first thing I'll do after loading the data will be to analyze and discover columns that only have 1 value OR where each value is unique. Excluding geographical coordinates, features such as "ID" or otherwise unique identifiers do not assist for analysis.<jupyter_code># reading file and info
df = pd.read_csv("data/vehicles.csv")
df.nunique()<jupyter_output><empty_output><jupyter_text>We can observe that **ID** and **URL** are unique for each vehicle posting, so we can go ahead and discard them since they do not provide any useful insights.Additionally, image analysis is not in the scope of this review, so I will also ignore **image_url**.Finally, as mentioned in scope, I will delete VIN, since for the scope for of this project I will not be querying the car's history.<jupyter_code># deleting useless columns.
df = df.drop(columns = ["id", "url", "image_url", "VIN"])<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code># Defining function to encode a dataframe.

def gen_array(dataframe : pd.DataFrame, cat_enc = None, num_enc = None) -> np.array:
    '''
    Generate a fit-ready array based on a dataframe.
    
    This function takes a pandas dataframe, separates it into "object" and "not object" slices.
    The slices are fit to their respect encoders, and finally both slices are concatenated.
    
    If encoders are not provided, they are instantiated inside the function, otherwise, use the provided encoders.
    
    Parameters
    ----------
    dataframe : Pandas Dataframe object
        Data structure to slice and encode.
    cat_enc : Scikit-learn categorical encoder (OneHot, Ordinal, or other.)
        Encoder for categorical data.
    num_enc : Scikit-learn numerical encoder (StandardScaler, MinMaxScaler, or other.)
        Encoder for numerical data.
    Returns
    -------
    a numpy array with the encoded attributes
    generated encoders
    
    '''
    if cat_enc is None and num_enc is None:
        
        # creating encoders for the data
        cat_encoder = OneHotEncoder(sparse = False, handle_unknown = "ignore") # setting Sparse to false because i'll need to concatenate the numerical variables later
        num_encoder = StandardScaler()
        
        cat_encoder.fit(dataframe.select_dtypes(include = "object"))
        num_encoder.fit(dataframe.select_dtypes(exclude = "object"))
        
    elif cat_enc is None:
        # scenario where only one encoder is provided
        num_encoder = num_enc
        cat_encoder = OneHotEncoder(sparse = False, handle_unknown = "ignore")
        cat_encoder.fit(dataframe.select_dtypes(include = "object"))
        
    elif num_enc is None:
        # scenario where only one encoder is provided
        cat_encoder = cat_enc
        num_encoder = StandardScaler()
        num_encoder.fit(dataframe.select_dtypes(exclude = "object"))
        
    else:
        cat_encoder = cat_enc
        num_encoder = num_enc

    cat_slice = cat_encoder.transform(dataframe.select_dtypes(include = "object"))
    num_slice = num_encoder.transform(dataframe.select_dtypes(exclude = "object"))
    
    array = np.concatenate([cat_slice, num_slice], axis=1)
    
    return array, cat_encoder, num_encoder

def inverse_trans(array : np.array, cat_enc, num_enc):
    '''
    Inverse an encoded array or arrays into it's parent un-encoded parent object or a Pandas Dataframe (more than one sample).
    
    This function takes an array as an input and performs inverse transforms to return the array to its parent components.
    The design of this function is limited to the transformations made in the function gen_array, thus this function is EXTREMELY LIMITED IN SCOPE.
    
    Parameters
    ----------
    array : numpy array object
        Data structure to slice and inverse_encode.
    cat_enc : Scikit-learn categorical encoder (OneHot, Ordinal, or other.)
        Encoder for categorical data.
    num_enc : Scikit-learn numerical encoder (StandardScaler, MinMaxScaler, or other.)
        Encoder for numerical data.   
    Returns
    -------
    a numpy array with the inversed_encoded attributes
    '''
    cat_features = array[:,:-2]
    num_features = array[:,-2:]
    
    if cat_features.shape[0] == 1:
        
        cat_features = cat_enc.inverse_transform(cat_features.reshape(1, -1))
        num_features = num_enc.inverse_transform(num_features.reshape(1, -1))
        
        fin_array = np.concatenate([cat_features, num_features], axis=1)
        
        return fin_array
        
    else:
        
        cat_features = cat_enc.inverse_transform(cat_features)
        num_features = num_enc.inverse_transform(num_features)
        
        fin_array = np.concatenate([cat_features, num_features], axis=1)
        
        return fin_array<jupyter_output><empty_output><jupyter_text>Data_exploration  Handling_missing_data<jupyter_code>df.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 426880 entries, 0 to 426879
Data columns (total 22 columns):
 #   Column        Non-Null Count   Dtype  
---  ------        --------------   -----  
 0   region        426880 non-null  object 
 1   region_url    426880 non-null  object 
 2   price         426880 non-null  int64  
 3   year          425675 non-null  float64
 4   manufacturer  409234 non-null  object 
 5   model         421603 non-null  object 
 6   condition     252776 non-null  object 
 7   cylinders     249202 non-null  object 
 8   fuel          423867 non-null  object 
 9   odometer      422480 non-null  float64
 10  title_status  418638 non-null  object 
 11  transmission  424324 non-null  object 
 12  drive         296313 non-null  object 
 13  size          120519 non-null  object 
 14  type          334022 non-null  object 
 15  paint_color   296677 non-null  object 
 16  description   426810 non-null  object 
 17  county        0 non-null       float64
 18  stat[...]<jupyter_text>Df.info() provides a good insight into the missing information, but let's try a different way to see how much missing data do we really have<jupyter_code># calculating all missing data points and normalizing to df size
(df.isna().sum() * 100 / df.shape[0]).to_frame().T # setting the count as an horizontally-oriented frame for better visualization<jupyter_output><empty_output><jupyter_text>As we can see, columns such as:* County -> 100% missing data* size -> 71.76% missing datahave large amounts of missing data that we can outright consider impossible to fix, so we will just delete them.<jupyter_code># I'll delete the columns with too much missing data to be recoverable
df = df.drop(columns = ["county", "size"])
# let's start exploring some of the remaining information to understand what kind of information do we have at hand
## To begin with, let's start exploring REGION & REGION_URL
df.info()
# Top 10 records in region
df["region"].unique()[:10]
# top 10 records in region URL
df["region_url"].unique()[:10]<jupyter_output><empty_output><jupyter_text>In general, region and region_url seems to be related to the seller's location and both hold pretty much the same information.So, for the time being, I will delete only one of them.Additionally, since region already holds geographical information, we will discard lat and long.<jupyter_code>df = df.drop(columns=["region_url", "lat", "long"])<jupyter_output><empty_output><jupyter_text>Now, before we continue exploring the rest of columns with missing data or with possibly relevant features, let's verify if the description column holds any easily extractable information.Remember we won't any in-depth NLP tasks for this project<jupyter_code>df["description"].unique()[:2]<jupyter_output><empty_output><jupyter_text>Unfortunately it seems the description field holds a lot of information, but in-depth NLP is not in the scope of this project.As such, I will delete description.<jupyter_code>df = df.drop(columns = "description")<jupyter_output><empty_output><jupyter_text>Now, **price** is our target so we will not perform any analysis in it.Let's explore some more features!<jupyter_code># Analyzing condition
df["condition"].value_counts(dropna = False)
# It seems condition is a simply a string describing the state of the vehicle, possibly from a dropdown menu. 
# In this case, instead of getting rid of the non-provided data we can fill with a placeholder string
df["condition"] = df["condition"].fillna(value = "Not relevant")
df["condition"].value_counts(dropna = False)
df["cylinders"].value_counts(dropna = False)<jupyter_output><empty_output><jupyter_text>For simplicity's sake, I'll provide cylinders the same treatment as "condition". But there seems to be a trend here, so I'll check some of the other variables.<jupyter_code>df["cylinders"] = df["cylinders"].fillna(value = "Not relevant")
df["cylinders"].value_counts(dropna = False)
df["drive"].value_counts(dropna = False)
df["type"].value_counts(dropna = False)
df["paint_color"].value_counts(dropna = False)
df["fuel"].value_counts(dropna = False)
df["title_status"].value_counts(dropna = False)
df["transmission"].value_counts(dropna = False)<jupyter_output><empty_output><jupyter_text>Drive, type, fuel, paint_color and the like have similarly high amounts of missing data, however we won't want to discard them as they might hold relevance for many sellers, so we will give them the same treatment as both condition and cylinders.In summary, the "Not Relevant" label is a placeholder for when someone does not wish to think about what they're choosing.<jupyter_code>df["drive"] = df["drive"].fillna(value = "Not relevant")
df["type"] = df["type"].fillna(value = "Not relevant")
df["paint_color"] = df["paint_color"].fillna(value = "Not relevant")
df["fuel"] = df["fuel"].fillna(value = "Not relevant")
df["title_status"] = df["title_status"].fillna(value = "Not relevant")
df["transmission"] = df["transmission"].fillna(value = "Not relevant")
# Let's recalculate the missing to see how much we have done
(df.isna().sum() * 100 / df.shape[0]).to_frame().T<jupyter_output><empty_output><jupyter_text>We have gotten rid of remarkable amount of missing data, but we still have some missing data. At this point I evaluated my options:The largest amount of remaining missing data comes from "MANUFACTURER". Now, in general terms this is considered a CRITICAL FEATURE: Regardless of statistical signals, due to SME everyone knows you can't possibly buy a car without knowing its manufacturer.Additionally, it seems we have 6055 different "models" to possibly infer the manufacturer. Although technically possible, I will consider this **unfeasible** for a one person, exploratory project. As such, **I have decided to delete all the rows where we have no manufacturer information**.Since there are no further actionables on the manufacturer front, I decided to **delete all the remaining records**.<jupyter_code>df = df.dropna()<jupyter_output><empty_output><jupyter_text>Exploratory_data_analysis To start the exploratory data analysis, let's explore how does our target variable (price) is distributed!<jupyter_code>plt.figure(figsize=(20,6))
sns.histplot(data = df["price"], bins = 20)
plt.ticklabel_format(style='plain', axis='x')<jupyter_output><empty_output><jupyter_text>Strangely enough we can see that something's with how our data is being represented. Why is our X axis going so far away? At first glance it seems we have some massive **OUTLIERS** that are preventing an accurate representation of the data.At this point, and based on personal criteria, I decided to visualize an outlier as **any vehicle worth more than 50,000**.<jupyter_code># Let's find out how many outliers we have, as a percentage
df["price"][df["price"] > 50000].shape[0] / df.shape[0] * 100<jupyter_output><empty_output><jupyter_text>Based on these results we can say that 2.78% of our data consists of outliers. In this respect, and considering my personal definition of outliers in regards to price, **I will delete all records for vehicles whose prices is above 50,000**. These prices are either typos, placeholders, or vehicles that cannot be trivially valued (classics, restorations, etc).<jupyter_code>df = df[df["price"] < 50000]<jupyter_output><empty_output><jupyter_text>I'll also do an additional sanity check and see how many vehicles are worth less than 1000 USD. I'll consider vehicles below 1k outliers if they are less than 5% of the remaining info.<jupyter_code>df[df["price"] < 1000].shape[0] / df.shape[0] * 100<jupyter_output><empty_output><jupyter_text>It's a difficult decision since **10%** of ALL samples are below 1,000 USD.After thoughtful consideration I decided to eliminate these samples because:* Vehicles worth / being sold at less than 1,000 USD are mostly broken cars or vehicles with major faults* There's a high possibility these are typos here (put 900 instead of 9,000) but this cannot be corrected if true.Among other things. As such, I will **delete ALL VEHICLES WITH A LISTED SALES PRICE OF LESS THAN 1,000 USD**<jupyter_code>df = df[df["price"] > 1000]
# Let's try visualizing the price distribution now that are rid of outliers 
plt.figure(figsize=(20,6))
ax = sns.histplot(data = df["price"], bins = 35, kde = True, stat = "probability")
# adding ticks every 5,000 USD
ax.xaxis.set_major_locator(MultipleLocator(2500))<jupyter_output><empty_output><jupyter_text>As it happens with most price related distributions, we can see that the majority of sales / prices are on the left, cheaper side. In particular, ~35% of sales are below 10,000 USD. A  mental sanity check says this is reasonable as these are used vehicles.Another interesting insight is that, despite the right-skew of the distribution, the slope of the descent or in other words, how far the data goes to the other side, is not very steep. I interpret this as good because it seems we have a captured and understood the underlying behaviour of the data.Now, my next step would to dive deeper into some of the features and see what else we can find and learn!<jupyter_code>df.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 346855 entries, 27 to 426879
Data columns (total 16 columns):
 #   Column        Non-Null Count   Dtype  
---  ------        --------------   -----  
 0   region        346855 non-null  object 
 1   price         346855 non-null  int64  
 2   year          346855 non-null  float64
 3   manufacturer  346855 non-null  object 
 4   model         346855 non-null  object 
 5   condition     346855 non-null  object 
 6   cylinders     346855 non-null  object 
 7   fuel          346855 non-null  object 
 8   odometer      346855 non-null  float64
 9   title_status  346855 non-null  object 
 10  transmission  346855 non-null  object 
 11  drive         346855 non-null  object 
 12  type          346855 non-null  object 
 13  paint_color   346855 non-null  object 
 14  state         346855 non-null  object 
 15  posting_date  346855 non-null  object 
dtypes: float64(2), int64(1), object(13)
memory usage: 45.0+ MB<jupyter_text>Although most of our data is in text strings, let's try running a correlation analysis to see any / importance of the few numerical features we have<jupyter_code>plt.figure(figsize=(20,6))
sns.heatmap(df.corr(), annot = True, cmap="PuOr")<jupyter_output><empty_output><jupyter_text>These preliminary observations are interesting and the results for odometer align with one would expect! That is, the more a car has been on the road, the less it is worth.Nevertheless, there is a problem with year here. Year's correlation is OK, but the direction is counter to logic. The problem might stem that, since these records do not belong a specific year, the car's manufacturing year isn't the most accurate metric for the cars age. In this case we need to calculate a proxy for the vehicles age.<jupyter_code># first, we'll convert the posting date to an actual timestamp
df["posting_date"] = pd.to_datetime(df["posting_date"], utc = True)
df["posting_date"]
# In this case what we really care is neither the cars year of manufacture NOR the posting date (this would always be *now* for a new prediction)
# The real variable is the cars AGE which would be -> Current_year - cars_mfg_year
df["age"] = df["posting_date"].dt.year - df["year"] # Extract the year from posting date and substract the cars MFG age to get its "age" at sale time
df["age"].value_counts()
# Let's run correlations again to see if there's any difference
plt.figure(figsize=(20,6))
sns.heatmap(df.corr(), annot = True, cmap="PuOr")<jupyter_output><empty_output><jupyter_text>As desired, the feature "age" captures the correct direction of the relationship as one would expect.In this situation, as the designer of this exercise, I have chosen to preserve **age** as the real predicting variable to maintain as it:* contains both MFG Year and year of inquire as part of a single variable* reduces the complexity (even if only by 1 feature) of the modelAs such, I will delete both "year" and "posting date"<jupyter_code>df = df.drop(columns = ["year", "posting_date"])<jupyter_output><empty_output><jupyter_text>Now that we have created the correct proxy variable, let's try visualizing how price behaves in regards to our numerical variables!<jupyter_code>sns.jointplot(data = df, x = "odometer", y = "price", height = 6);
sns.jointplot(data = df, x = "age", y = "price", height = 6);<jupyter_output><empty_output><jupyter_text>The relationships are difficult to visualize with the jointplots, but with some effort we can observe hints at the relationships captured by the correlations. After analyzing the numerical features, I will now start exploring the categorical data!In this case there's no correlations, but I will start by visualizing distributions based on countplots.<jupyter_code># separating column names for easier plotting
cat_cols = list(df.select_dtypes(include="object").columns)
num_cols = list(df.select_dtypes(include="float64"))
# This will be used only to plot the categorical columns so 6 x 2
fig, ax = plt.subplots(2,6, figsize=(25,10))

# separating rows for easier processing
row_1 = ax[0]
row_2 = ax[1]

# creating plots for the first row
for index, column in enumerate(cat_cols[:6]):
    # creating ax figure in order of value
    sns.countplot(ax=row_1[index], data = df, x=column, order = df[column].value_counts().index)
    # setting plot title
    row_1[index].set_title(column)
    # hiding xaxis since they become unreadable
    row_1[index].get_xaxis().set_visible(False)

# creating plots for the second row
for index, column in enumerate(cat_cols[6:]):
    # creating ax figure in order of value
    sns.countplot(ax=row_2[index], data = df, x=column, order = df[column].value_counts().index)
    # setting plot title
    row_2[index].set_title(column)
    # hiding xaxis since they become unreadable
    row_2[index].get_xaxis().set_visible(False)
    
fig.suptitle("Count of unique instances per feature");<jupyter_output><empty_output><jupyter_text>At first glance there are some interesting findings:* **TITLE_STATUS, FUEL & TRANSMISSION**: are HIGHLY unbalanced features. Although there is a large number of samples from which to draw conclusions, these imbalances make me think these variables will have negligible predicting power.* **MODEL**: This feature does not display large counts. Let's do a deep dive and see what we find.* **REGION & STATE**: Hold very similar (if not the same) information, yet they have different distributions. Let's see if we can choose to keep only one of them, but for now we'll hold it.<jupyter_code>df["model"].nunique()
df["model"].value_counts()<jupyter_output><empty_output><jupyter_text>It seems we have an important problem with model.There are too many unique models! And it's not a topic of unique vehicle models, but, the model description is far too heterogenous. Let's see what we can do about it!**A POSTERIORI CORRECTION** In this case I will be looking at how many models have less than 30 samples<jupyter_code>df["model"].value_counts()[df["model"].value_counts() < 30]<jupyter_output><empty_output><jupyter_text>Wow, there are over 19,000 unique instances of models with less than 30 appearances. This kind of information will definitely pollute our models' learning.Let's remove all these features from the base.<jupyter_code># To be sure, let's find out exactly how many records we are going to delete
df["model"].value_counts()[df["model"].value_counts() < 30].sum()
# extracting the string representation of all the models which less than 30 appearances
bad_models = list(df["model"].value_counts()[df["model"].value_counts() < 30].index)

# converting into an array for faster processing
bad_models = np.array(bad_models)

# creating a boolean mask. Before further processing we'll verify if this number matches what we expect
mask = df["model"].isin(bad_models)
mask.sum()
# Let's also precalculate how many rows we'll have at the hend
df.shape[0] - mask.sum()
# Now that i'm sure the mask is ok, I'll reserve it.
mask = ~mask

# Finally, let's do the filtering out the polluting models from the data
df = df[mask]

# verifying correct shape
df.shape[0]
df["model"].nunique()<jupyter_output><empty_output><jupyter_text>Now it's time to do the same but for manufacturers, let's see how many (if any) have less than 150 samples<jupyter_code>df["manufacturer"].value_counts()<jupyter_output><empty_output><jupyter_text>Fortunately it seems only 3 manufacturers match our criteria for elimination!As such, I will exclude these manufacturers from the predictions models.<jupyter_code># let's extract only those with less than 150 samples
bad_mfg = list(df["manufacturer"].value_counts()[df["manufacturer"].value_counts() < 150].index)
#once again, convert to array
bad_mfg = np.array(bad_mfg)

# Let's create the mask
mask_mfg = df["manufacturer"].isin(bad_mfg)

# Let's invert it
mask_mfg = ~mask_mfg

# finally let's pass it

df = df[mask_mfg]

# Let's see our final shape
df.shape
# Finding unique models per manufacturer
df.groupby("manufacturer")["model"].nunique().sort_values(ascending = False)<jupyter_output><empty_output><jupyter_text>I have greatly reduced the complexity of the information contained in model! Hopefully this is enough to move forward! After our findings and explorationsFortunately, Scikit-learn provides some convenient tools for this:* f_regression -> which is basically a one-way Anova.* mutual_info_regression -> Analogue to its class counterpart, but for continuous variablesbut first, from this point onwards, I will create the Train, Validation and Test splits! This is necessary as we **want to avoid inferences that we are not sure would apply to a ML / DL model that has not seen new data**. As such, the relevant statistics will be obtained using training sets only!<jupyter_code># Separating target and features
y = df["price"]
X_parent = df.drop(columns = "price")
# creating first split. The result here is the untouched test splits. The "splits" are a preliminary split that will be separated again for train and validation
X_split, X_test, y_split, y_test = train_test_split(X_parent, y, train_size = 0.8, random_state = 20) # adding random state for reproducibility
# creating the 2nd split. This one creates the train and validation splits
X_train, X_val, y_train, y_val = train_test_split(X_split, y_split, train_size = 0.75, random_state = 20) # adding random state for reproducibility
# first we need to encode the categorical and features since we can't process them as it is
## instantiating scalers
ord_enc = OrdinalEncoder()
norm_sca = StandardScaler()

## separating categorical and numerical features
cat_feats = X_train.select_dtypes(include="object")
num_feats = X_train.select_dtypes(include="float64")

## fitting scalers
ord_enc.fit(cat_feats)
norm_sca.fit(num_feats)

## re-assigning variables
cat_feats = ord_enc.transform(cat_feats)
num_feats = norm_sca.transform(num_feats)

## concatenating arrays
X_stats = np.concatenate([cat_feats, num_feats], axis=1)
f_arrray = f_regression(X_stats, y_train)
mut_info = mutual_info_regression(X_stats, y_train)
# creating a results table
cat_cols = list(X_train.select_dtypes(include="object").columns) # extreacting categorical columns
num_cols = list(X_train.select_dtypes(include="float64")) # extracting numerical columns
cols = cat_cols + num_cols # setting columns in correct order and then merging lists
pd.DataFrame([f_arrray[1], mut_info], index = ["f_regression", "mutual_regre"], columns = cols) #I'm only keeping the p-values of f_regression<jupyter_output><empty_output><jupyter_text>Before explaining the results found from this analysis I will provide some contextual information:* f_regression from sklearn is an univariate test that mimics an ANOVA F-TEST but for a continuous problem (like price). More importantly, and as outlined in the functions page itself **IT SHOULDNT BE USED AS A STAND ALONE SELECTION METHOD**. In essence, this trains to calculate how much of the variance of target (price) is explained by the variance of a feature. The value presented in the results table is the p-value associated to the test.* mutual information for regression is a nonparametrics (does not assume any distribution for the data) test that attempts to find the features that contribute to most to a change in the target feature (price)* Interestingly, while many of the features had large amounts of missing values which were replaced with placeholder strings, they reflect different behaviours, providing some confidence that the transformation's did not distort the results too much.That said, the results of both tests provide some interesting, but inconclusive insights:1. Both methods lack consistency in determining which features might be relevant / contribute more to the variation in price, HOWEVER they do have some agreement into **which features DO NOT contribute much to the variation of the price**. In this sense, the results can be used in reverse to remove features that add little prediction power.2. Based on the previous point the following features could be considered of LOW relevance:    1. **TITLE_STATUS**: Has a high p-value and low mutual_information score; by both metrics this feature does not contribute much to price. Unfortunately I'm unfamiliar with these intricancies of this feature and what it represents. **On principle, this feature will be deleted from all splits**     2. **Region / State**: Has a high p-value and low mutual_information score. by both metrics this feature does not contribute much to price. Interestingly both features capture very similar (if not the same information), as such, one can consider it to be slightly redundant. **On principle, one of these features will be deleted from all splits.   After these analysis, we have some strong signals about which features to get rid of. These inputs can be combined with the observed in the EDA to finally get rid of some features!<jupyter_code># deleting unnecessary vars and cleaning namespace
del ord_enc 
del norm_sca
del cat_feats
del num_feats
del X_stats
del cat_cols
del num_cols
del cols
del X_parent
del X_train
del X_val
del X_test
del y
del y_val
del y_train
del y_test<jupyter_output><empty_output><jupyter_text>All that said, the following features will not be considered any further for the purposes of this analysis:* Region* Title Status<jupyter_code>df = df.drop(columns = ["region", "title_status"])<jupyter_output><empty_output><jupyter_text>At this point I consider the EDA and feature selection process finished and I can start training models!  Model_training A few considerations will be explained before entering the training, fitting and evaluating models phase.Although I'm not aiming for a specific model or solution, one the sub-goals of this project is to build a Neural Network with Keras to create the prediction model. Nevertheless, training a Neural Network is computationally expensive, thus I chose the following approach:* I would first training sklearn Linear models (with default parameters) and use them to:    * Find out baseline performance    * Estimate feature weights and possible run a 2nd trial of feature selection* Then I would run ensemble models (with default parameters) to:    * Obtain a secondary baseline performance based on averaged results    * Obtain feature weights (often included in the attributes of tree-based models* Later, after performing some final adjustments, I would create the neural network.* Finally, I would compare the results of all models and choose what works best and fine-tune hyperparameters if necessary.<jupyter_code>## Separating target and features
y = df["price"]
X_parent = df.drop(columns = "price")

## creating TRUE SPLITS. The result here is the untouched test splits. The "splits" are a preliminary split that will be separated again for train and validation
X_split, X_test, y_split, y_test = train_test_split(X_parent, y, train_size = 0.8, random_state = 40) # adding random state for reproducibility

## creating the 2nd split. This one creates the train and validation splits
X_train, X_val, y_train, y_val = train_test_split(X_split, y_split, train_size = 0.75, random_state = 40) # adding random state for reproducibility

## creating a copy of train set for future purposes
X_save = X_train
# encoding training arrays 
X_train, cat_encoder, num_encoder = gen_array(X_train)
# encoding valdition array with pre-trained encoders
X_val, cat_encoder, num_encoder = gen_array(X_val, cat_enc = cat_encoder, num_enc = num_encoder)
# recreating standarized indexes to be able to explore faulty values
y_val2 = y_val.reset_index().drop(columns = "index").squeeze()
y_test2 = y_test.reset_index().drop(columns = "index").squeeze()
# creating a list of column names in the order of concatenation for easier analysis
cat_cols = list(X_parent.select_dtypes(include = "object").columns)
num_cols = list(X_parent.select_dtypes(exclude = "object").columns)
cols = cat_cols + num_cols<jupyter_output><empty_output><jupyter_text>Linear_models These are the steps I will follow with the linear models:* Instantiate every model in a corresponding variable* Fit all models, one by one* Use the fitted models to generate predictions* Calculate and store residuals for visualization purposes* Calculate the following errors for each model:    * Mean Absolute Error (MAE)    * Mean Squared Error (MSE)    * Root Mean Squared Error (RMSE)<jupyter_code>lasso = LassoCV(verbose = 1, n_jobs = -1)
linear = LinearRegression(n_jobs = -1)
lasso.fit(X_train, y_train)
linear.fit(X_train, y_train)
# obtaining predictions
lasso_pred = lasso.predict(X_val)
linear_pred = linear.predict(X_val)
# calculating residuals
lasso_res = lasso_pred - y_val2
linear_res = linear_pred - y_val2
# creating a summary table to visualize results
MAE = [mean_absolute_error(y_val, lasso_pred), mean_absolute_error(y_val, linear_pred)]
MSE = [mean_squared_error(y_val, lasso_pred), mean_squared_error(y_val, linear_pred)]
RMSE = [np.sqrt(mean_squared_error(y_val, lasso_pred)), np.sqrt(mean_squared_error(y_val, linear_pred))]

pd.DataFrame({"MAE": MAE, "MSE": MSE, "RMSE": RMSE}, index = ["Lasso", "Linear"])<jupyter_output><empty_output><jupyter_text>Interestingly it seems a simple linear model with no regularization is better at finding relationships within the data! Let's try plotting the residuals to see how this looks!<jupyter_code>fig, ax = plt.subplots(1,2, figsize=(24,8))
sns.lineplot(ax = ax[0], data = lasso_res, color = "green")
sns.lineplot(ax = ax[1], data = linear_res, color = "purple");<jupyter_output><empty_output><jupyter_text>Observing our results, and more importantly the residuals we can observe there are some MASSIVE outliers. This tells me two things:* The model's are currently unstable, but it seems to be related to some very specific values that may be fixable or that can be omitted from the training sets..* Unstability aside, linear models seem to describe the data decently well. Strangely, though, the linear model with NO regularization seems to be more effective at predictions.That said, let's explore the causes of this instability and if we can fix our process a-posteriori.<jupyter_code># Let's find the data points that are causing massive instability in our models.
## The outlier detection will be generated bi-directionally although the plots indicate the problem is that rediction understimates value.
## We only want the indeces since i'm only interested in retrieving them from the parent array
unstable_vals_idx = lasso_res[(lasso_res < -50000) | (lasso_res > 50000)].index

## retriving the actual unstable values
unstable_features = X_val[unstable_vals_idx,:]
unstable = inverse_trans(unstable_features, cat_enc = cat_encoder, num_enc = num_encoder)
pd.DataFrame(unstable, columns = cols)<jupyter_output><empty_output><jupyter_text>Unfortunately analyzing the outliers yielded no results of significance. From observing this points I can't detect any considerable trends, so I will move forward with what i have.  Ensemble_regressors These are the steps I will follow ensemble regressors:* Instantiate every model in a corresponding variable* Fit all models, one by one* Use the fitted models to generate predictions* Calculate and store residuals for visualization purposes* Calculate the following errors for each model:    * Mean Absolute Error (MAE)    * Mean Squared Error (MSE)    * Root Mean Squared Error (RMSE)<jupyter_code>forest = RandomForestRegressor(verbose = 1, n_jobs = -1, random_state = 20) #Random state for reproducibility
boost = GradientBoostingRegressor(verbose = 1, random_state = 20) #Random state for reproducibility
forest.fit(X_train, y_train)
boost.fit(X_train, y_train)
# obtaining predictions
forest_pred = forest.predict(X_val)
boost_pred = boost.predict(X_val)
# calculating residuals
forest_res = forest_pred - y_val2.values
boost_res = boost_pred - y_val2.values
# creating a summary table to visualize results
MAE = [mean_absolute_error(y_val, forest_pred), mean_absolute_error(y_val, boost_pred)]
MSE = [mean_squared_error(y_val, forest_pred), mean_squared_error(y_val, boost_pred)]
RMSE = [np.sqrt(mean_squared_error(y_val, forest_pred)), np.sqrt(mean_squared_error(y_val, boost_pred))]

pd.DataFrame({"MAE": MAE, "MSE": MSE, "RMSE": RMSE}, index = ["Random_Forest", "Boosted_Trees"])
fig, ax = plt.subplots(1,2, figsize=(24,8))
sns.lineplot(ax = ax[0], data = forest_res, color = "green")
sns.lineplot(ax = ax[1], data = boost_res, color = "purple");<jupyter_output><empty_output><jupyter_text>Interestingly, the non-parametric seem to be much better at detecting trends in the underlying data. This leads me to two conclusions* the data is NOT linearly-related, which means linear models are bad estimators.* The exact nature of the relation is undefined, so we need something that can approximate anything...a Neural Network!  Keras_model<jupyter_code>model = Sequential()

# Creating input layer
model.add(Input(shape = (1802,), name = "input_lay"))

# creating first layer
model.add(Dense(900, activation = "relu"))
# regularization on first one to dropout all the possibly useless attributes
model.add(Dropout(0.4))

# creating a second layer with a substantially smaller size
model.add(Dense(100, activation = "relu"))
# regularizing
model.add(Dropout(0.1))

# third layer with a substantially smaller size
model.add(Dense(11, activation = "relu"))
# no regularization in the third layer

# final layer
model.add(Dense(1))

# compiling
model.compile(optimizer = "adam", loss = "mse", metrics = [tf.keras.metrics.MeanAbsoluteError(), tf.keras.metrics.RootMeanSquaredError()])
stop = EarlyStopping(monitor = "val_loss", patience = 3, mode = "min")
model.fit(x = X_train, y = y_train.values, validation_data = (X_val, y_val.values), batch_size = 128, epochs = 100, callbacks = [stop])<jupyter_output>Epoch 1/100
1287/1287 [==============================] - 20s 16ms/step - loss: 68727584.0000 - mean_absolute_error: 5289.6025 - root_mean_squared_error: 8290.2100 - val_loss: 19639276.0000 - val_mean_absolute_error: 2947.0520 - val_root_mean_squared_error: 4431.6221
Epoch 2/100
1287/1287 [==============================] - 20s 15ms/step - loss: 18088816.0000 - mean_absolute_error: 2840.9019 - root_mean_squared_error: 4253.0947 - val_loss: 15109528.0000 - val_mean_absolute_error: 2508.0359 - val_root_mean_squared_error: 3887.0977
Epoch 3/100
1287/1287 [==============================] - 20s 15ms/step - loss: 15470683.0000 - mean_absolute_error: 2594.5283 - root_mean_squared_error: 3933.2788 - val_loss: 13837034.0000 - val_mean_absolute_error: 2396.9192 - val_root_mean_squared_error: 3719.8164
Epoch 4/100
1287/1287 [==============================] - 20s 15ms/step - loss: 14489422.0000 - mean_absolute_error: 2500.2266 - root_mean_squared_error: 3806.4973 - val_loss: 12964396.0000 - val_mean[...]<jupyter_text>Model_evaluation Now that we have trained a couple of models we can go ahead test all of them versus the test set to determine which is / are the best for our task.<jupyter_code># Encoding the test set for processing.

X_test, cat_encoder, num_encoder = gen_array(X_test, cat_enc = cat_encoder, num_enc = num_encoder)
# creating predictions
forest_test_pred = forest.predict(X_test)
boost_test_pred = boost.predict(X_test)
lasso_test_pred = lasso.predict(X_test)
Keras_test_pred = model.predict(X_test)
# list of results 
MAE = [mean_absolute_error(y_test, forest_test_pred), mean_absolute_error(y_test, boost_test_pred), 
       mean_absolute_error(y_test, lasso_test_pred), mean_absolute_error(y_test, Keras_test_pred)]
MSE = [mean_squared_error(y_test, forest_test_pred), mean_squared_error(y_test, boost_test_pred), 
       mean_squared_error(y_test, lasso_test_pred), mean_squared_error(y_test, Keras_test_pred)]
RMSE = [np.sqrt(mean_squared_error(y_test, forest_test_pred)), np.sqrt(mean_squared_error(y_test, boost_test_pred)), 
       np.sqrt(mean_squared_error(y_test, lasso_test_pred)), np.sqrt(mean_squared_error(y_test, Keras_test_pred))]
pd.DataFrame({"MAE": MAE, "MSE": MSE, "RMSE": RMSE}, index = ["Forest", "Gradient Trees", "Lasso", "Keras"])
# Interestingly our models produce the expected results as per the validation sets behaviour. As such, I will attempt to find the best model via residual visualization
forest_test_res = y_test.values - forest_test_pred
boost_test_res = y_test.values - boost_test_pred
lasso_test_res = y_test.values - lasso_test_pred
Keras_test_res = y_test.values - Keras_test_pred.reshape(-1,)
fig, ax = plt.subplots(1,4, figsize = (24,6))

sns.lineplot(ax = ax[0], data = forest_test_res, color = "green")
ax[0].set_title("Forest Estimator")

sns.lineplot(ax = ax[1], data = boost_test_res, color = "red")
ax[1].set_title("Gradient Boosted Trees Estimator")

sns.lineplot(ax = ax[2], data = lasso_test_res, color = "purple")
ax[2].set_title("Lasso Regression Estimator")

sns.lineplot(ax = ax[3], data = Keras_test_res, color = "blue")
ax[3].set_title("Keras Feed-Forward NN Estimator")
# saving keras model
save_model(model = model, filepath = "Keras Model")
# saving categorical and numerical encoders
cat_encoder_file = "categorical_encoder.pkl"
num_encoder_file = "numerical_encoder.pkl"
categories_file = "categories_file.pkl"

with open(cat_encoder_file, "wb") as file:
    dump(cat_encoder, file)

with open(num_encoder_file, "wb") as file:
    dump(num_encoder, file)
    
with open(categories_file, "wb") as file:
    dump(X_save, file)<jupyter_output><empty_output><jupyter_text>Feature_Importance<jupyter_code>features_imp = pd.DataFrame(data = forest.feature_importances_[:-2].reshape(1,-1), columns = cat_encoder.get_feature_names())
# For every column we will consider the feature importance to be the sum of the feature importance of all its members
manuf_imp = features_imp.iloc[:,:36].T.sum() 
model_imp = features_imp.iloc[:,36:1692].T.sum() 
condition_imp = features_imp.iloc[:,1692:1699].T.sum() 
cylinders_imp = features_imp.iloc[:,1699:1708].T.sum()
fuel_imp = features_imp.iloc[:,1708:1714].T.sum()
transmission_imp = features_imp.iloc[:,1714:1718].T.sum()
drive_imp = features_imp.iloc[:,1718:1722].T.sum()
type_imp = features_imp.iloc[:,1722:1736].T.sum()
color_imp = features_imp.iloc[:,1736:1749].T.sum()
state_imp = features_imp.iloc[:,1749:1800].T.sum()

values = {"manufacturer": manuf_imp, "model":model_imp, "condition": condition_imp, "cylinders":cylinders_imp,
         "fuel": fuel_imp, "transmission":transmission_imp, "drive":drive_imp, "type":type_imp, "color":color_imp, "state":state_imp}
# creating the final visualization dataframe
feature_importance = pd.DataFrame(values)
numerical_importance = pd.DataFrame(forest.feature_importances_[-2:].reshape(1,-1),columns = ["odometer", "age"])
feature_importances = pd.concat([feature_importance, numerical_importance], axis=1)
feature_importances<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08562452605351532, 'fraction_numerical': 0.02151446213844654, 'mean_word_length': 5.042152114151067, 'pattern_counts': {'":': 19, '<': 172, '<?xml version=': 0, '>': 174, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '41', 'dir': 'jupyter-structured-clean-dedup', 'id': '589074', 'license': 'MIT', 'n_tokens_mistral': 14227, 'n_tokens_neox': 13021, 'n_words': 5215, 'path': 'notebook/Used_cars_Dev.ipynb', 'repo_name': 'alexkingsing/used_cars', 'size': '46181'}