<jupyter_start><jupyter_text>V. Feature Selection, Algorithm improvement, Generalization demo.In this part, we create a model, we try different algorithms and see which one delivers the best results. Then we chose the best algorithm and fine tune it. This notebook presents the following parts:    1) Model creation    2) Feature Selection (hybrid method)    3) Algorithm testing        - verify reduce dataset gives goods results    4) Improve the model (fine tune hyperparameters)        - verify improved model delivers a better metric    5) Save Model for inference    6) Show evidence the new model will generalize well<jupyter_code># import libraries
import pandas as pd
from sklearn import preprocessing
import sklearn.model_selection as ms
from sklearn import linear_model
import sklearn.metrics as sklm
import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as ss
import math
from sklearn.model_selection import train_test_split 
from sklearn.linear_model import LinearRegression
from sklearn import metrics
import xgboost as xgb
from sklearn.metrics import roc_auc_score, r2_score

from sklearn.ensemble import RandomForestRegressor

from sklearn.feature_selection import VarianceThreshold
from sklearn.feature_selection import SelectKBest, SelectPercentile

from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
from sklearn.metrics import roc_auc_score, mean_squared_error
from sklearn.feature_selection import f_classif, f_regression
from sklearn.metrics import mean_squared_log_error



%matplotlib inline
%matplotlib inline
#import data set
df=pd.read_csv('dftobeimproved.csv')
df.shape
df.head(2)<jupyter_output><empty_output><jupyter_text>**1) Model Creation**<jupyter_code># separate train and test sets
X_train, X_test, y_train, y_test = train_test_split(
    df.drop(labels=['attractiveness'], axis=1),
    df['attractiveness'],
    test_size=0.3,
    random_state=0)

X_train.shape, X_test.shape

# I keep a copy of the dataset with all the variables
# to measure the performance of machine learning models
# at the end of the notebook

X_train_original = X_train.copy()
X_test_original = X_test.copy()<jupyter_output><empty_output><jupyter_text>**2) Features Selection (hybrid method:Hybrid recursive feature addition)** From the previous notebook I know the best model is the Random Forest, so I create a hybrid method to select the best features for this algorithm.<jupyter_code># the first step of this procedure consists in building
# a machine learning algorithm using all the available features
# and then determine the importance of the features according
# to the algorithm

# set the seed for reproducibility
seed_val = 1000000000
np.random.seed(seed_val)

# build initial model using all the features
model_all_features = RandomForestRegressor()

model_all_features.fit(X_train, y_train)

# calculate the roc-auc in the test set
y_pred_test = model_all_features.predict(X_test)
r2_score_all = r2_score(y_test, y_pred_test)
print('Test all features rdf R2 = %f' % (r2_score_all))
# the second step consist of deriving the importance of 
# each feature and ranking them from the most to the least
# important

# get feature name and importance
features = pd.Series(model_all_features.feature_importances_)
features.index = X_train.columns

# sort the features by importance
features.sort_values(ascending=False, inplace=True)

# plot
features.plot.bar(figsize=(20,6))
# view the list of ordered features
features = list(features.index)
features
# next, we need to build a machine learning
# algorithm using only the most important feature

# set the seed for reproducibility
seed_val = 1000000000
np.random.seed(seed_val)

# build initial model using all the features
model_one_feature = RandomForestRegressor()

# train using only the most important feature
model_one_feature.fit(X_train[features[0]].to_frame(), y_train)

# calculate the roc-auc in the test set
y_pred_test = model_one_feature.predict(X_test[features[0]].to_frame())
r2_score_first = r2_score(y_test, y_pred_test)
print('Test one feature RandomForest R2=%f' % (r2_score_first))
# the final step consists in adding one at a time
# all the features, from the most to the least
# important, and build an xgboost at each round.

# once we build the model, we calculate the new r2
# if the new r2 is bigger than the original one
# (with one feature), then that feature that was added
# was important, and we should keep it.
# otherwise, we should remove the feature

# recursive feature addition:

# first we arbitrarily set the increase in r2
# if the increase is above this threshold,
# the feature will be kept
tol = 0.0015

print('Doing recursive feature addition')

# we initialise a list where we will collect the
# features we should keep
features_to_keep = [features[0]]

# set a counter to know how far ahead the loop is going
count = 1

# now we loop over all the features, in order of importance:
# remember that features is the list of ordered features
# by importance
for feature in features[1:]:
    print()
    print('testing feature: ', feature, ' which is feature ', count,
          ' out of ', len(features))
    count = count + 1

    # initialise model
    model_int = RandomForestRegressor()

    # fit model with the selected features
    # and the feature to be evaluated
    model_int.fit(
        X_train[features_to_keep + [feature] ], y_train)

    # make a prediction over the test set
    y_pred_test = model_int.predict(
        X_test[features_to_keep + [feature] ])

    # calculate the new r2
    r2_score_int = r2_score(y_test, y_pred_test)
    print('New Test R2={}'.format((r2_score_int)))

    # print the original roc-auc with all the features
    print('All features Test R2={}'.format((r2_score_first)))

    # determine the drop in the roc-auc
    diff_r2 = r2_score_int - r2_score_first

    # compare the increase in r2 with the tolerance
    # we set previously
    if diff_r2 >= tol:
        print('Increase in r2 = {}'.format(diff_r2))
        print('keep: ', feature)
        print
        # if the increase in the r2 is bigger than the threshold
        # we keep the feature and re-adjust the r2 to the new value
        # considering the added feature
        r2_score_first = r2_score_int
        
        # and we append the feature to keep to the list
        features_to_keep.append(feature)
    else:
        # we ignore the feature
        print('Increase in r2 = {}'.format(diff_r2))
        print('remove: ', feature)
        print


# now the loop is finished, we evaluated all the features
print('DONE!!')
print('total features to keep: ', len(features_to_keep))
# capture the 8 selected features
seed_val = 1000000000
np.random.seed(seed_val)

# build initial model
rdf_reduced = RandomForestRegressor()

# fit the model with the selected features
rdf_reduced.fit(X_train[features_to_keep], y_train)

# make predictions
y_pred_test = rdf_reduced.predict(X_test[features_to_keep])

# calculate roc-auc
r2_score_final = r2_score(y_test, y_pred_test)
print('Test selected features R2=%f' % (r2_score_final))
X_train[features_to_keep].columns<jupyter_output><empty_output><jupyter_text>We have reduced the features space from 18 to 8 **2) Algoirthm testing** **Random Forest Regressor** (with reduced dataset)<jupyter_code>#print RMSLE
print ('RDF reduced RMSLE is', np.sqrt(mean_squared_log_error(y_test, y_pred_test)))<jupyter_output>RDF reduced RMSLE is 0.048126906820446916<jupyter_text>**Random Forest Regressor** (with original dataset)<jupyter_code># I compare with the result of the random forest with the original datset
Original_Forest = RandomForestRegressor()
Original_Forest.fit(X_train_original,y_train)
Original_Forest_predictions = Original_Forest.predict(X_test_original)

print ('Original Random Forest Regression RMSLE is', np.sqrt(mean_squared_log_error(y_test, Original_Forest_predictions)))<jupyter_output>Original Random Forest Regression RMSLE is 0.04844630778458226<jupyter_text>Reducing the dataset by half of the features, we get the same metric.  **3) Improving the model** 3.1 Find best model hyperparameters<jupyter_code>rdf_final = RandomForestRegressor(random_state = 42)
from pprint import pprint

# Look at parameters used by our current forest
print('Parameters currently in use:\n')
pprint(rdf_final.get_params())
# Create Random Hyperparamter Grid
from sklearn.model_selection import RandomizedSearchCV
# Number of trees in random forest
n_estimators = [int(x) for x in np.linspace(start = 200, stop = 2000, num = 10)]
# Number of features to consider at every split
max_features = ['auto', 'sqrt']
# Maximum number of levels in tree
max_depth = [int(x) for x in np.linspace(10, 110, num = 11)]
max_depth.append(None)
# Minimum number of samples required to split a node
min_samples_split = [2, 5, 10]
# Minimum number of samples required at each leaf node
min_samples_leaf = [1, 2, 4]
# Method of selecting samples for training each tree
bootstrap = [True, False]
# Create the random grid
random_grid = {'n_estimators': n_estimators,
               'max_features': max_features,
               'max_depth': max_depth,
               'min_samples_split': min_samples_split,
               'min_samples_leaf': min_samples_leaf,
               'bootstrap': bootstrap}
pprint(random_grid)
# Random Search Training
# Use the random grid to search for best hyperparameters
# First create the base model to tune
rdf_final = RandomForestRegressor()
# Random search of parameters, using 3 fold cross validation, 
# search across 100 different combinations, and use all available cores
rdf_final = RandomizedSearchCV(estimator = rdf_final, param_distributions = random_grid, n_iter = 100, cv = 3, verbose=2, random_state=42, n_jobs = -1)
# Fit the random search model
rdf_final.fit(X_train[features_to_keep], y_train)
rdf_final.best_params_<jupyter_output><empty_output><jupyter_text>4.2 Updating best parameters on exisitng best algorithm<jupyter_code>from sklearn.ensemble import RandomForestRegressor
Forestimproved = RandomForestRegressor(n_estimators = 2000, min_samples_split=10, 
                                       min_samples_leaf=2, max_features='auto', max_depth=40, bootstrap=True)
Forestimproved.fit(X_train[features_to_keep], y_train)
Forestimproved_predictions = Forestimproved.predict(X_test[features_to_keep])
print ('Forest Improved Regression RMSLE is', np.sqrt(mean_squared_log_error(y_test, Forestimproved_predictions)))<jupyter_output>Forest Improved Regression RMSLE is 0.046995274291810746<jupyter_text>4.3 Check for **improvement on RMSLE**<jupyter_code>print ('Original Random Forest Regression RMSLE is', np.sqrt(mean_squared_log_error(y_test, Original_Forest_predictions)))
print ('RDF reduced RMSLE is', np.sqrt(mean_squared_log_error(y_test, y_pred_test)))
print ('Forest Improved Regression RMSLE is', np.sqrt(mean_squared_log_error(y_test, Forestimproved_predictions)))
print ('RMSLE improvement from Original:',(np.sqrt(mean_squared_log_error(y_test, Original_Forest_predictions))) - (np.sqrt(mean_squared_log_error(y_test, Forestimproved_predictions))))<jupyter_output>Original Random Forest Regression RMSLE is 0.04844630778458226
RDF reduced RMSLE is 0.048126906820446916
Forest Improved Regression RMSLE is 0.046995274291810746
RMSLE improvement from Original: 0.0014510334927715154<jupyter_text>**5. Save model improved**<jupyter_code># Saving LogRec model
from sklearn.externals import joblib
from joblib import dump, load
dump(Forestimproved,'Forestimproved.joblib')<jupyter_output><empty_output><jupyter_text>**6. Verification that the model created will generalize well**<jupyter_code>def print_metrics(y_train, y_test, n_parameters):
    ## First compute R^2 and the adjusted R^2
    r2 = sklm.r2_score(y_train, y_test)
    r2_adj = r2 - (n_parameters - 1)/(y_train.shape[0] - n_parameters) * (1 - r2)
    
    ## Print the usual metrics and the R^2 values
    print('Mean Square Error      = ' + str(sklm.mean_squared_error(y_train, y_test)))
    print('Root Mean Square Error = ' + str(math.sqrt(sklm.mean_squared_error(y_train, y_test))))
    print('Mean Absolute Error    = ' + str(sklm.mean_absolute_error(y_train, y_test)))
    print('Median Absolute Error  = ' + str(sklm.median_absolute_error(y_train, y_test)))
    print('R^2                    = ' + str(r2))
    print('Adjusted R^2           = ' + str(r2_adj))
   
y_score = Forestimproved.predict(X_test[features_to_keep]) 
print_metrics(y_test, y_score, 28)<jupyter_output>Mean Square Error      = 0.005740182869225286
Root Mean Square Error = 0.07576399454374938
Mean Absolute Error    = 0.05952382991064678
Median Absolute Error  = 0.04926644898238958
R^2                    = 0.6490595812272444
Adjusted R^2           = 0.6248876646281005<jupyter_text>Given the difference between R^2 and Adjusted R2 is so small, we can safely infer that the model will generalize well.<jupyter_code>#Lets look at the residuals plot 

## Calculate residuals (y-yhat)
df['residuals']=y_train - y_test
def resid_qq(y_train, y_test):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test, y_score)
    ## now make the residual plots
    ss.probplot(resids, plot = plt)
    plt.title('Residuals vs. predicted values')
    plt.xlabel('Predicted values')
    plt.ylabel('Residual')
    
resid_qq(y_train, y_test)<jupyter_output><empty_output><jupyter_text>As with the histogram, the Q-Q Normal plot indicates the residuals are close to Normally distributed, show some skew (deviation from the straight line). This is particularly for large residuals.<jupyter_code>def hist_resids(y_test, y_score):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test, y_score)
    ## now make the residual plots
    sns.distplot(resids)
    plt.title('Histogram of residuals')
    plt.xlabel('Residual value')
    plt.ylabel('count')
    
hist_resids(y_test, y_score)<jupyter_output><empty_output><jupyter_text>This plot looks reasonable. The residual values are distributed near both sides of 0.<jupyter_code>def resid_plot(y_test, y_score):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test, y_score)
    ## now make the residual plots
    sns.regplot(y_score, resids, fit_reg=False)
    plt.title('Residuals vs. predicted values')
    plt.xlabel('Predicted values')
    plt.ylabel('Residual')

resid_plot(y_test, y_score)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07969009407858328, 'fraction_numerical': 0.02863862755949087, 'mean_word_length': 4.922572716099959, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '79704', 'license': 'MIT', 'n_tokens_mistral': 4823, 'n_tokens_neox': 4430, 'n_words': 1429, 'path': '6. Feature Selection and Algorithm improvement.ipynb', 'repo_name': 'titoniubo/x_biamow', 'size': '14456'}