<jupyter_start><jupyter_text>Regression - predicting house prices> A notebook exploring a housing dataset and experimenting with various regression models.- toc: true - badges: false- comments: true- categories: [regression]- image: images/regression.png- author: <NAME>  SummaryThe aim of this project is to experiment with data manipulation, adequately prepare a dataset for use and apply various regression models to the manipulated data in order to build an accurate regression model. I also employed hyperparameter tuning to the best model to further increase performance.All in all, I developed a regression model with an R2 score of ~0.87  Full code can be found here.  Library imports<jupyter_code>#collapse-hide
import numpy as np
import pandas as pd
import sklearn as sk
from sklearn import preprocessing
import matplotlib.pyplot as plt
from sklearn import linear_model
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
import seaborn as sns
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
import warnings
from sklearn.svm import LinearSVR
from sklearn.preprocessing import StandardScaler
from sklearn import svm, datasets
from sklearn.model_selection import GridSearchCV
from sklearn.inspection import permutation_importance

warnings.filterwarnings(action='ignore')
pd.set_option("max_columns", 25)<jupyter_output><empty_output><jupyter_text>Data input and data cleaningThe first stage is to use pandas to read in the dataset into a dataframe, investigate the data, and replace any null values with 0's.<jupyter_code>df = pd.read_csv('../projects/regression/kc_house_data.csv')
df.fillna(0, inplace=True)
df<jupyter_output><empty_output><jupyter_text>Next, I wanted simplify the 'date' feature to just capture the year the house was recorded, and then use this in conjunction with the 'year_built' feature to create one new feature; 'age'.<jupyter_code>d =[]
for i in df['date'].values:
    d.append(i[:4])
    
df['date'] = d
df['date']=df['date'].astype(float)

df['age'] = df['date'] - df['yr_built']<jupyter_output><empty_output><jupyter_text>The next stage in data cleaning was to remove any features which were obviously redundant. Here, I removed the 'date' and 'yr_built' features which I had already simplified, alongside the 'id' feature.<jupyter_code>df = df.drop(["date", "id", 'yr_built'],  axis=1)<jupyter_output><empty_output><jupyter_text>Removing correlated featuresNext, I wrote some code to remove pairs of features which showed a high degree of correlation (>80%). I did this as my understanding is that for regression problems, multicollinearity is bad. After doing some more research, it seems that multicollinearity does not fundamentally worsen a model's peformance, but rather makes interpretation of coefficients/analysis of the features much harder. I actually found a slight decrease in peformance when removing the correlated features, so for optimal peformance this code can be commented out. I have kept the code in case I ever want to dig deeper into interpreting the influence of certain features on the model.<jupyter_code># corr_features =[]

# for i , r in df.corr().iterrows():
#     k=0
#     for j in range(len(r)):
#         if i!= r.index[k]:
#             if r.values[k] >=0.5:
#                 corr_features.append([i, r.index[k], r.values[k]])
#         k += 1

# feat =[]
# for i in corr_features:
#     if i[2] >= 0.8:
#         feat.append(i[0])
#         feat.append(i[1])

# df.drop(list(set(feat)), axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Visualise the distributions of each featureAt this stage, I thought it would be useful to have a look at what each feature's distribution looks like.<jupyter_code>df.hist(figsize=(30,20))
plt.show()<jupyter_output><empty_output><jupyter_text>Removing outliersI was conflicted about including this code. On one hand, removing outliers definitely improves the model's peformance. However, I am also of the opinion that almost all data holds value, so throwing some away doesn't sit great with me.In the end I became a sucker for peformance, and decided to remove the outliers from the datasets. Beyond just peformance gains, there is rationale for doing this. Regressions are designed to summarise the data, and having outliers greatly swing the summary of the vast majority of the dataset isn't always logical. Sometimes new, outlying data points won't be described well by a model that has omitted outliers, but it should boost peformance in the general case.In summary, this issue is something that needs to be considered in a case by case basis. As long as there is explained reason behind the decision to remove outliers (or not), then there is good justification to go either way based on the circumstances.In terms of implementation, as the price follows a somewhat gaussian distribution, I used a standard deviation cut-off approach to removing outliers. In short, each price was normalized, and if this value was greater than or less than 3 then I removed that row from the dataset. I've also visualised the outliers, where blue is low outliers and red is high outliers.<jupyter_code>zscore = []
outlier =[]
threshold = 3

price_mean = np.mean(df['price'])
price_std = np.std(df['price'])

for i in df['price']:
        z = (i-price_mean)/price_std
        zscore.append(z)
        if np.abs(z) > threshold:
            outlier.append(i)


plt.figure(figsize = (10,6))
sns.distplot(zscore, kde=False)
plt.axvspan(xmin = -3 ,xmax= min(zscore),alpha=0.2, color='blue', label='Lower Outliers')
plt.axvspan(xmin = 3 ,xmax= max(zscore),alpha=0.2, color='red', label='Upper Outliers')
plt.show()

dj=[]
for i in df.price:
    if i in set(outlier):
        dj.append(0.0)
    else:
        dj.append(i)
        
df['P'] = dj

df = df.drop(df[df['P'] == 0.0].index)<jupyter_output><empty_output><jupyter_text>After the outliers were removed, I split the dataset into a feature set and a price set, and removed the redundant features appropriatley.<jupyter_code>X = df.drop(['price','P'], axis=1)
Y = df['price']<jupyter_output><empty_output><jupyter_text>Preparing the training and test datasetsThe first step was to easily split the dataframes into train and test sets, on a 70:30 ratio.<jupyter_code>x_train, x_test, y_train, y_test = train_test_split(X, Y, test_size = 0.3)<jupyter_output><empty_output><jupyter_text>Next, to deal with categorical values (only zipcode in this case), I employed one-hot encoding to categorical features in the feature sets.<jupyter_code>train_zip = pd.get_dummies(x_train.zipcode, prefix="zipcode")
x_train = x_train.drop('zipcode',axis=1)
x_train = x_train.join(train_zip)

test_zip = pd.get_dummies(x_test.zipcode, prefix="zipcode")
x_test = x_test.drop('zipcode',axis=1)
x_test = x_test.join(test_zip)<jupyter_output><empty_output><jupyter_text>The final stage before model selection was to scale the data. I went for the StandardScaler, which implements a basic scaling based off the mean and standard deviation of the datasets. It is important to note that the scaling fit to the training set is the same scaling applied to the test set. This is to avoid the model learning anything about the test set before we let it attack it.<jupyter_code>#collapse-output
scaler = StandardScaler()
scaler.fit(x_train)
scaler.transform(x_train)
scaler.transform(x_test)<jupyter_output><empty_output><jupyter_text>Model selectionTo start I defined a set of models which I wanted to include in my model evaluation. These included the ones I was most familiar with, or had heard good things about.<jupyter_code>lr = linear_model.LinearRegression()
ridge = linear_model.Ridge()
lasso = linear_model.Lasso()
e_net = linear_model.ElasticNet()
r_forest = RandomForestRegressor(n_jobs=-1)
k_neighbours = KNeighborsRegressor()
d_tree = DecisionTreeRegressor()


models = [lr, ridge, lasso, e_net, r_forest, k_neighbours, d_tree]
model_names = ['Linear', 'Ridge', 'Lasso', 'Elastic Net', 'Random Forest', 'KNN', 'Decision Tree']<jupyter_output><empty_output><jupyter_text>Next, I created a set of cross-validation sets via the sklearn KFold validator class. This allowed me to generate a preliminary R2 score for each model via cross-validation, without having to fully test each model (which would have taken an age!).<jupyter_code>for i in range(len(models)):
    kfold = KFold(n_splits=7)
    result = cross_val_score(models[i], x_train, y_train, cv=kfold, scoring='r2')
    print(model_names[i] + ': ' + str(result.mean()) )<jupyter_output>Linear: 0.834841451524125
Ridge: 0.834788741056295
Lasso: 0.8348479388140282
Elastic Net: 0.6154086881998547
Random Forest: 0.868900272979691
KNN: 0.4468420436933553
Decision Tree: 0.742939379397493<jupyter_text>Hyperparameter tuningRandom Forest Regression came out top, with the highest R2 score. As such, in the intertest of time, I made the assumption that once tuned, this would still be the highest peforming model.In this step, I chose the three most prudent parameters which I thought would affect the peformance of the model; number of trees, criterion (function which measures the quality of a split) and number of features considered when looking for optimal splits.<jupyter_code>parameters = {
    'n_estimators': [10, 50, 100],
    'criterion': ['mse', 'mae'],
    'max_features': ['auto', 'sqrt', 'log2'],
            }<jupyter_output><empty_output><jupyter_text>Next, I used GridSearchCV to find the optimal set of my selected parameters; 100, mse and auto respectivley. Admittedly the first of these, number of estimators, was always going to be the largest value, so this could have been omitted from the search. If I wanted to optimise efficiency I could have plotted the peformance as a function of the number of estimators and found the point at which peformance gains become minimal.<jupyter_code>grid = GridSearchCV(r_forest, param_grid = parameters,n_jobs=-1, scoring='r2', verbose=2)
grid.fit(x_train, y_train)
print(grid.best_score_)
print(grid.best_params_)<jupyter_output>Fitting 5 folds for each of 18 candidates, totalling 90 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 8 concurrent workers.
[Parallel(n_jobs=-1)]: Done  25 tasks      | elapsed:   13.6s
[Parallel(n_jobs=-1)]: Done  90 out of  90 | elapsed: 23.0min finished
0.8686275134918929
{'criterion': 'mse', 'max_features': 'auto', 'n_estimators': 100}<jupyter_text>Final modelNow that I have selected the optimal model and tuned its parameters, its a simple process to train and use the model on the dataset. I've printed the R2 score here for completeness.Note that I have included an extra parameter 'n_jobs' which defines how many concurrent processes to run; -1 uses all available CPU cores so *shoud* speed it up!<jupyter_code>r_forest = RandomForestRegressor(criterion='mse', max_features='auto', n_estimators=100,n_jobs=-1)
r_forest.fit(x_train,y_train)
predictions = r_forest.predict(x_test)
score = r2_score(y_test,predictions)
print(score)<jupyter_output>0.8701934934737463<jupyter_text>Finally, I can use the model to make a new price prediction of a new dummy data point (for ease I have made a new dummy point that is the average values across all the testing set data points).<jupyter_code>new_pred = r_forest.predict([np.array(x_test).mean(axis=0)])
print(new_pred)<jupyter_output>[568417.02]<jupyter_text>I did consider visualising some of the decision trees in the regressor at this point, however the depth of the trees, in addition to the number of trees makes this fairly difficult to interpret, so I have left this out.  Feature importanceAs an added bonus, I have taken a look at how significant each feature was in the model. Firstly, I plotted all feature importances, and saw that the zipcode features (which were numerous, at the end of the feature list) had low significance.<jupyter_code>feature_importances = r_forest.feature_importances_
plt.bar([x for x in range(len(feature_importances))], feature_importances)<jupyter_output><empty_output><jupyter_text>Removing these from the chart allows us a better view of feature significance.<jupyter_code>importances_no_zips = feature_importances[0:17]

relevant_features = x_train.columns.tolist()[0:17]

print(relevant_features)

fig = plt.figure(figsize=[30,10])
plt.bar(relevant_features, importances_no_zips)<jupyter_output>['bedrooms', 'bathrooms', 'sqft_living', 'sqft_lot', 'floors', 'waterfront', 'view', 'condition', 'grade', 'sqft_above', 'sqft_basement', 'yr_renovated', 'lat', 'long', 'sqft_living15', 'sqft_lot15', 'age']
{'fraction_non_alphanumeric': 0.08678919258561106, 'fraction_numerical': 0.021442035815268615, 'mean_word_length': 5.539804827940421, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 78, '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': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '282708', 'license': 'Apache-2.0', 'n_tokens_mistral': 4076, 'n_tokens_neox': 3762, 'n_words': 1467, 'path': '_notebooks/2021-06-09-Regression.ipynb', 'repo_name': 'ah161652/data_science_projects', 'size': '12736'}