<jupyter_start><jupyter_text>Lab 6 - Feature Selection Analysis   <NAME>   10.16.2020   OverviewThis lab has the purpose of comparing four different models in regards to feature selection. The models used are: Logistic regression including all features, used as the base model; Logistic regression with stepwise selection based on p-values; Random forest based on principal compnent analysis feature selection; Random Forest based on feature importance.Logistic regression is a model that uses features to predict the output of a categorical variable.Random forest is an ensemble learning technique that consists of a collection of a number of decision trees. It outputs the class that is the mode of the classes for a categorical target variable.For this lab, Random forest will be used in combination with two different methods: feature importance and principal component analysis (PCA).Gini importance is a feature selection based on the random forest classifier that gives a score based on importance to each feature.PCA is a technique used to reduce the dimension of datasets, creating new variables that is a combination of the orignal variables while still being able to explain a great portion of the variance within the dataset.The performance measure used to compare the models will be the [accuracy score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html) imported from sklearn.metrics and a comparison of the ROC/AUC curve between all models.  DataThe caravan data set used in this lab contains information on customers of an insurance company. It contains 86 features (variables). The data was collected to attempt to answer if it's possible to predict who would be interested in buying a carvan insurance policy with the target variable being CARAVAN.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
from sklearn import linear_model
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
from sklearn.metrics import accuracy_score
from sklearn.feature_selection import SelectFromModel
%matplotlib inline
# Load the caravan dataset
url = 'https://raw.githubusercontent.com/WillKoehrsen/feature-selector/master/data/caravan-insurance-challenge.csv'

df = pd.read_csv(url)
df.head()
all_vars = df.drop(['ORIGIN'], axis=1)
#Group only the input variables 
x_vars = all_vars.drop(['CARAVAN'], axis=1)
#Get the target variable
y = df['CARAVAN']

print(y[0:5])
x_vars.head()<jupyter_output>0    0
1    0
2    0
3    0
4    0
Name: CARAVAN, dtype: int64<jupyter_text>Exploratory Data Analysis<jupyter_code># Check for missing values
missing = np.isnan(all_vars.values).any()
if ( missing == False):
    print('No missing values')
else:
    print('Oh uh. There are missing values!')
# Summary of the input variables
x_vars.describe()
# Histogram of the target variable
y.hist()

# Set the title and labels
plt.xlabel("Policy Status")
plt.ylabel("Number of People")
plt.title("Number of People with Caravan's Insured vs Not Insured")

# show the plot
plt.show()<jupyter_output><empty_output><jupyter_text>The histogram above shows imbalance in the dataset with respect to the target variable. Should the data be balanced before attempting to create any model?<jupyter_code>#Attempt to capture only the highly correlated features ( > 0.6)
#https://stackoverflow.com/questions/26463714/pandas-get-combination-of-columns-where-correlation-is-high
corr_matrix = x_vars.corr().abs()

indices = np.where(corr_matrix > 0.6)
indices = [(corr_matrix.index[x], corr_matrix.columns[y]) for x, y in zip(*indices)
                                        if x != y and x < y]
len(indices)
print(indices)<jupyter_output>[('MOSTYPE', 'MOSHOOFD'), ('MGEMOMV', 'MFALLEEN'), ('MGEMOMV', 'MFWEKIND'), ('MGODPR', 'MGODGE'), ('MRELGE', 'MRELOV'), ('MRELGE', 'MFALLEEN'), ('MRELOV', 'MFALLEEN'), ('MRELOV', 'MAUT0'), ('MFALLEEN', 'MFWEKIND'), ('MOPLHOOG', 'MOPLLAAG'), ('MOPLHOOG', 'MSKA'), ('MOPLMIDD', 'MOPLLAAG'), ('MOPLLAAG', 'MSKA'), ('MOPLLAAG', 'MSKC'), ('MBERHOOG', 'MSKA'), ('MBERHOOG', 'MZFONDS'), ('MBERHOOG', 'MZPART'), ('MBERARBG', 'MSKC'), ('MHHUUR', 'MHKOOP'), ('MAUT1', 'MAUT0'), ('MZFONDS', 'MZPART'), ('MINKM30', 'MINKGEM'), ('MINK7512', 'MINKGEM'), ('PWAPART', 'AWAPART'), ('PWABEDR', 'AWABEDR'), ('PWALAND', 'AWALAND'), ('PPERSAUT', 'APERSAUT'), ('PBESAUT', 'ABESAUT'), ('PMOTSCO', 'AMOTSCO'), ('PVRAAUT', 'AVRAAUT'), ('PAANHANG', 'AAANHANG'), ('PTRACTOR', 'ATRACTOR'), ('PWERKT', 'AWERKT'), ('PBROM', 'ABROM'), ('PLEVEN', 'ALEVEN'), ('PPERSONG', 'APERSONG'), ('PGEZONG', 'AGEZONG'), ('PWAOREG', 'AWAOREG'), ('PBRAND', 'ABRAND'), ('PZEILPL', 'AZEILPL'), ('PPLEZIER', 'APLEZIER'), ('PFIETS', 'AFIETS'), ('PINB[...]<jupyter_text>The code above shows that there are 44 pairs of features with a correlation greater than |0.6|. This is an indication of the presence of collinearity in the dataset.  Models<jupyter_code>#Firt split the data 70% train and 30% test

X_train, X_test, y_train, y_test = train_test_split(x_vars, y, test_size=0.3, random_state=0)
## BALANCED DATA SET using SMOTE (COMMENT this whole cell to used imbalanced data)
from imblearn.over_sampling import SMOTE
#
smt = SMOTE()
X_train, y_train = smt.fit_sample(X_train, y_train)<jupyter_output>/usr/local/lib/python3.6/dist-packages/sklearn/externals/six.py:31: FutureWarning: The module is deprecated in version 0.21 and will be removed in version 0.23 since we've dropped support for Python 2.7. Please rely on the official version of six (https://pypi.org/project/six/).
  "(https://pypi.org/project/six/).", FutureWarning)
/usr/local/lib/python3.6/dist-packages/sklearn/utils/deprecation.py:144: FutureWarning: The sklearn.neighbors.base module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.neighbors. Anything that cannot be imported from sklearn.neighbors is now part of the private API.
  warnings.warn(message, FutureWarning)
/usr/local/lib/python3.6/dist-packages/sklearn/utils/deprecation.py:87: FutureWarning: Function safe_indexing is deprecated; safe_indexing is deprecated in version 0.22 and will be removed in version 0.24.
  warnings.warn(msg, category=FutureWarning)<jupyter_text>Logistic regression using all features<jupyter_code># Create logist regression object
logreg = linear_model.LogisticRegression(max_iter=10000)

# Train the model using the training sets
logreg.fit(X_train, y_train)

# Make predictions using the testing set
y_pred = logreg.predict(X_test)

#Print the accuracy score of the logist model
print('Logistic Regression Accuracy: {:.2f}'.format(logreg.score(X_test, y_test)))

# ROC Curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
logit_roc_auc = roc_auc_score(y_test, logreg.predict(X_test))
fpr, tpr, thresholds = roc_curve(y_test, logreg.predict_proba(X_test)[:,1])
plt.figure()
plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()<jupyter_output><empty_output><jupyter_text>Logistic Regression with stepwise based on p-values<jupyter_code>#https://datascience.stackexchange.com/questions/24405/how-to-do-stepwise-regression-using-sklearn/24447#24447

import statsmodels.api as sm

X = x_vars


def stepwise_selection(X, y, 
                       initial_list=[], 
                       threshold_in=0.01, 
                       threshold_out = 0.05, 
                       verbose=True):
    """ Perform a forward-backward feature selection 
    based on p-value from statsmodels.api.OLS
    Arguments:
        X - pandas.DataFrame with candidate features
        y - list-like with the target
        initial_list - list of features to start with (column names of X)
        threshold_in - include a feature if its p-value < threshold_in
        threshold_out - exclude a feature if its p-value > threshold_out
        verbose - whether to print the sequence of inclusions and exclusions
    Returns: list of selected features 
    Always set threshold_in < threshold_out to avoid infinite looping.
    See https://en.wikipedia.org/wiki/Stepwise_regression for the details
    """
    included = list(initial_list)
    while True:
        changed=False
        # forward step
        excluded = list(set(X.columns)-set(included))
        new_pval = pd.Series(index=excluded)
        for new_column in excluded:
            model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included+[new_column]]))).fit()
            new_pval[new_column] = model.pvalues[new_column]
        best_pval = new_pval.min()
        if best_pval < threshold_in:
#            best_feature = new_pval.argmin()
            best_feature = new_pval.idxmin()
            included.append(best_feature)
            changed=True
            if verbose:
                print('Add  {:30} with p-value {:.6}'.format(best_feature, best_pval))

        # backward step
        model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
        # use all coefs except intercept
        pvalues = model.pvalues.iloc[1:]
        worst_pval = pvalues.max() # null if pvalues is empty
        if worst_pval > threshold_out:
            changed=True
            worst_feature = pvalues.argmax()
            included.remove(worst_feature)
            if verbose:
                print('Drop {:30} with p-value {:.6}'.format(worst_feature, worst_pval))
        if not changed:
            break
    return included

result = stepwise_selection(X, y)

print('resulting features:')
print(result)
new_xvars = x_vars[['PPERSAUT', 'MKOOPKLA', 'PWAPART', 'APLEZIER', 'MOPLHOOG', 'PBRAND', 'MBERBOER', 'MRELGE', 'PWALAND', 'ABRAND', 'AZEILPL', 'MINK123M', 'PBYSTAND', 'PGEZONG', 'AGEZONG', 'MHHUUR']]
new_xvars.head()
# Split the data 70% train and 30% test

#X_train, X_test, y_train, y_test = train_test_split(new_xvars, y, test_size=0.3, random_state=0)
# Create logist regression object
logreg_pval = linear_model.LogisticRegression(max_iter=10000)

# Train the model using the training sets
logreg_pval.fit(X_train, y_train)

# Make predictions using the testing set
y_pred = logreg_pval.predict(X_test)

#Print the accuracy score of the logist model
print('Logistic Regression Accuracy: {:.2f}'.format(logreg_pval.score(X_test, y_test)))
# ROC Curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
logit_roc_auc = roc_auc_score(y_test, logreg_pval.predict(X_test))
fpr, tpr, thresholds = roc_curve(y_test, logreg_pval.predict_proba(X_test)[:,1])
plt.figure()
plt.plot(fpr, tpr, label='Logistic Regression (p-val based) (area = %0.2f)' % logit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()<jupyter_output><empty_output><jupyter_text>Random Forest based on Principal component analysis<jupyter_code>from sklearn.preprocessing import StandardScaler

# standardize the data 
X_std = StandardScaler().fit_transform(x_vars)
mean_vec = np.mean(X_std, axis=0)
cov_mat = (X_std - mean_vec).T.dot((X_std - mean_vec)) / (X_std.shape[0]-1)
print('Covariance matrix \n%s' %cov_mat)
#Perform eigendecomposition on covariance matrix
cov_mat = np.cov(X_std.T)
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
print('Eigenvectors \n%s' %eig_vecs)
print('\nEigenvalues \n%s' %eig_vals)
pca = PCA().fit(all_vars)
plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
pca = PCA(0.90).fit(all_vars)
pca.n_components_<jupyter_output><empty_output><jupyter_text>From the results above, using PCA we were able to reduce the predicting variables from 85 to 13 and still be able to explain about 90% of the variance.<jupyter_code>#https://towardsdatascience.com/machine-learning-step-by-step-6fbde95c455a
# Split the data 70% train and 30% test

X_train_std, X_test_std, y_train, y_test = train_test_split(X_std, y, test_size=0.3, random_state=0)

## Balance Standardized data (COMMENT the next line to used imbalanced data)
X_train_std, y_train = smt.fit_sample(X_train_std, y_train)

# Reduce the number of features from 85 to 13 in the X_train_std and X_test_std data
pca = PCA(n_components=13)
pca.fit(X_train_std)
X_train_pca = pca.transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
# Fit train data to Random Forest Model

from sklearn.ensemble import RandomForestClassifier

rfc_pca = RandomForestClassifier(n_estimators=1000, random_state=0)
rfc_pca.fit(X_train_pca, y_train)

y_pred = rfc_pca.predict(X_test_pca)
accuracy_score(y_test, y_pred)
# ROC Curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
rfc_pca_roc_auc = roc_auc_score(y_test, rfc_pca.predict(X_test_pca))
fpr, tpr, thresholds = roc_curve(y_test, rfc_pca.predict_proba(X_test_pca)[:,1])
plt.figure()
plt.plot(fpr, tpr, label='Random Forest with PCA (area = %0.2f)' % rfc_pca_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()<jupyter_output><empty_output><jupyter_text>Random Forest based on feature importance feature selection<jupyter_code>#Firt split the data 70% train and 30% test

X_train, X_test, y_train, y_test = train_test_split(x_vars, y, test_size=0.3, random_state=0)

## Balance data (COMMENT the next line to used imbalanced data)
X_train, y_train = smt.fit_sample(X_train, y_train)

# Create a list of feature names
feat_labels = x_vars.columns
# Create a random forest classifier
clf = RandomForestClassifier(n_estimators=1000, random_state=0)

# Train the classifier
clf.fit(X_train, y_train)


# Print the name and gini importance of each feature
for feature in zip(feat_labels, clf.feature_importances_):
    print(feature)
# Create a selector object that will use the random forest classifier to identify
# features that have an importance of more than 0.02
sfm = SelectFromModel(clf, threshold=0.02)

# Train the selector
sfm.fit(X_train, y_train)
# Print the names of the most important features
for feature_list_index in sfm.get_support(indices=True):
    print(feat_labels[feature_list_index])
# Transform the data to create a new dataset containing only the most important features
# Note: We have to apply the transform to both the training X and test X data.
X_important_train = sfm.transform(X_train)
X_important_test = sfm.transform(X_test)
# Create a new random forest classifier for the most important features
clf_important = RandomForestClassifier(random_state=0, n_jobs=-1)

# Train the new classifier on the new dataset containing the most important features
clf_important.fit(X_important_train, y_train)
# Apply The Full Featured Classifier To The Test Data
y_important_pred = clf_important.predict(X_important_test)

# View The Accuracy Of Our Limited Feature (2 Features) Model
accuracy_score(y_test, y_important_pred)
# ROC Curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
rfc_feat_roc_auc = roc_auc_score(y_test, clf_important.predict(X_important_test))
fpr, tpr, thresholds = roc_curve(y_test, clf_important.predict_proba(X_important_test)[:,1])
plt.figure()
plt.plot(fpr, tpr, label='Random Forest with PCA (area = %0.2f)' % rfc_feat_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11270147488953887, 'fraction_numerical': 0.016366917667558652, 'mean_word_length': 5.032282282282282, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '642408', 'license': 'Apache-2.0', 'n_tokens_mistral': 5661, 'n_tokens_neox': 5316, 'n_words': 1538, 'path': 'feature_selection.ipynb', 'repo_name': 'kiehhub/notebooks', 'size': '16084'}