Executive Summary¶
This analysis seeks three objectives:¶
- (1) describe, visualize, and model a complex dataset using supervised machine learning techniques.
- (2) describe, visualize, and segment a complex dataset using unsupervised machine learning techniques.
- (3) extract business insights from applied SQL analysis.
First, this analysis surveyed the distributions with the following selected findings:¶
- x1 (independent variable):
- uniform distribution
space
- uniform distribution
- x2 (independent variable):
- normal distribution
space
- normal distribution
- x3 (independent variable):
- t-distribution
space
- t-distribution
- x5 (independent variable):
- exponential distribution
space
- exponential distribution
- x6 (independent variable):
- binomial distribution
space
- binomial distribution
- ya (dependent variable):
- chi2 distribution
space
- chi2 distribution
- yb (dependent variable):
- reverse log normal distribution
space
- reverse log normal distribution
- yc (dependent variable):
- logistic distribution
Regarding predicting the dependent variables, this analysis found the following:¶
- For ya:
- For yb:
- For yc:
- No model appropriate because yc does not appear to be related to the independent variables.
Second, this analysis found through unsupervised machine learning:¶
- that the unsupervised learning data set could be reduced from 10 dimensions to just 3 (depicted as spheroids below).
- via a Principal Component Analysis (PCA):
- that the observations could be grouped into 4 cohorts (depicted in hues below).
- via a K-Means Clustering Analysis:
Third, this analysis derived business insights through SQL:¶
- such as the following query finding the number of posts in 2017:
- such as the following query finding the distribution of views for users:
- such as the following query finding the top 10 tags associated with posts created in 2017:
- such as the following query finding the average response time in seconds from the time a 2017 question first post to when the official answer is accepted:
Credit to God, my Mother, family and friends.
All errors are my own.
Best,
George John Jordan Thomas Aquinas Hayward, Optimist
Data Scientist
September 30, 2019
Table of Contents¶
- Executive Summary
- Selected Data Visualizations
- Part I. - Supervised Machine Learning
- Part II. - Unsupervised Machine Learning
- Part III. - SQL Coding
- Thank You!
- #TheLordAlwaysDelivers
Selected Data Visualizations¶
Surveying the DS1 Distributions ¶
DS1 Pair Plot ¶
DS1 Correlation Matrix ¶
Using XGBoost for 'ya' Feature Selection ¶
Using XGBoost for 'yb' Feature Selection ¶
Multiple Linear Regression for 'ya' Prediction ¶
Simple Linear Regression for 'yb' Prediction ¶
Selecting Number of Components for DS2 PCA ¶
3D-Plotting DS2's K-Means Cohort Results After PCA¶
SQL Love-Language Queries¶
And it begins...¶
In [1]:
import pandas as pd
import numpy as np
import missingno as msno
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import warnings
warnings.filterwarnings("ignore")
from fitter import Fitter
from scipy import stats
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn import linear_model
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, r2_score, median_absolute_error, \
explained_variance_score, confusion_matrix, accuracy_score, precision_score, recall_score
import xgboost as xgb
import statsmodels.formula.api as sm
from statsmodels.api import add_constant
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import FunctionTransformer
import matplotlib as mpl
from mpl_toolkits import mplot3d
from sklearn.cluster import KMeans
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
import hypertools as hyp
Part I. - Supervised Machine Learning ¶
- Consider data set 1 (ds1.csv). The data set comprises features (the Five xs) along with three sequences that may or may not be generated from the features (3 ys).
Question 1 - Part I.A ¶
- Describe the data set in a few sentences. E.g. What are the distributions of each feature? Summary statistics?
👇🏾 We begin by reading in the data. 👇🏾¶
- We will check for NULLS.
- For some reason there is an index column in the raw data, but a Pandas dataframe comes with an index column by default so we will drop this 'extra' index column in the raw data.
- We will summarize the data.
- We will check the distributions.
In [2]:
#read in data
ds1 = pd.read_csv("ds1.csv")
In [3]:
#peek at at data
ds1.head(3)
Out[3]:
In [4]:
#drop "Unamed:0" column and peek at the data again:
ds1 = ds1.drop(ds1.columns[0], axis=1)
ds1.head(3)
Out[4]:
In [5]:
#check for NULLS
msno.matrix(ds1, color = (.0,.0,.2))
Out[5]:
💭👆🏾This indicates that there are no NULLS.👆🏾💭¶
👇🏾 We now describe the data. 👇🏾¶¶
- When we talk about describing the data, we look at numerical summaries. These include:
- the standard deviation to get a sense of the spread of the data
- the mean and median
- the min and max
- the null count
- Pandas has a nice 'describe' method, which we'll use here.
- Next we'll get to the distributions.
In [6]:
ds1.describe()
Out[6]:
👇🏾 We now will now inspect and fit the distributions. 👇🏾¶¶
- We will plot histograms to visualize the distributions.
- We will use the Python third party package "fitter" to check these distributions against some of the most common and useful statistical distributions.
- Fitter was developed by Thomas Cokelaer
🔬 Function 1: Visualizing & Fitting the Distribution 🔬
- I am about to write a function which will:
- first, plot a histogram for our five independent and three dependent variables
- second, use fitter to fit a distribution to our histogram
- third, just for good measure, add a probability plot
- I will loop it through each column in the data frame to do all this automatically.
- See you on the other side!
In [7]:
def distribution_visualizer_and_fitter(df, col_str):
#this part makes the histogram:
fig, axs = plt.subplots(figsize=(12,4))
plt.hist(df[col_str], bins=100)
plt.title('Histogram of {}'.format(col_str), fontweight = 'bold')
plt.ylabel('Frequency')
plt.xlabel('Values')
plt.show()
#this part fits the distribution to the histogram:
f = Fitter(df[col_str], distributions=['gamma', 'beta', 'rayleigh', 'norm', 'pareto', 'uniform', \
'logistic', 'expon', 'chi2', 'pearson3'], verbose = False, timeout = 10)
f.fit()
f.summary()
#this part makes the probability plot
fig, axs = plt.subplots(ncols=2, figsize=(12,4))
axs[1].set_axis_off()
stats.probplot(df[col_str], plot=axs[0])
plt.show()
In [8]:
for i in ds1.columns:
distribution_visualizer_and_fitter(ds1, i)
💭👆🏾Conclusion: What do now know about the distributions?👆🏾💭¶
- x1 (independent variable):
- uniform distribution
space
- uniform distribution
- x2 (independent variable):
- normal distribution
space
- normal distribution
- x3 (independent variable):
- t-distribution
- some will also call this a 'nearly normal' distribution
- can be modeled with beta distribution
space
- t-distribution
- x5 (independent variable):
- exponential distribution
- a kind of gamma distribution
space
- a kind of gamma distribution
- exponential distribution
- x6 (independent variable):
- binomial distribution
space
- binomial distribution
- ya (dependent variable):
- chi2 distribution
- can also say it's a positively skewed normal distribution, BUT:
- to me, the key is the long tail on the right, looks a lot like chi2 distribution
space
- chi2 distribution
- yb (dependent variable):
- reverse log normal distribution
- can be approximated with beta distribution
space
- yc (dependent variable):
- logistic distribution
- more unclear
- in my subsequent analysis, I'll discuss how I don't think this 'dependent' feature is actually based on the other independent feature
- however, from my fitting algorithm, it appears a logistic distribution is the closest
- logistic distribution
Question 2 - Part I.B ¶
- Try to come up with a predictive model, e.g. y = f(x_1 , … , x_n) for each y sequence. Describe your model and how you came up with them. What (if any) are the predictive variables? How good would you say each of your models is?
👇🏾 Modeling Gameplan 👇🏾¶¶
- First a few assumptions:
- The problem seems to ask for an equation as the description of the model
- I assume this means that the client needs a linear regression output
- Many clients like this because, even if there are more accurate models, a linear regression is highly interpretable
- Nevertheless, we will still want to begin with a decision tree model because we can use it to figure out which features are most important
- I will use the decision tree models as a kind of feature selection algorithm
- I will proceed in two steps:
- 1. Exploratory Data Analysis
- Check correlations
- Check pair plots
- Check if polynomial transformations or interactions are helpful here
- I am NOT going to do log transformations because there are too many negative numbers
- You can add a constant to correct this issue, but you need to really know the underlying data well to do so; and this is something regarding which the Data Scientist should discuss with the team.
- Why? Because there are some data sets that even though the shift by adding a constant will not alter the variance, it could alter the interpretability or the meaning of the numbers
- Since we don't know too much about these underlying numbers (as in what they actually represent in real life), I make the decision, for here, not to apply a log transformation
- I am NOT going to do log transformations because there are too many negative numbers
- To select features I will use:
- A XGBoost feature importance
- Backwards elimination regression, if needed
- Please note that since this is just EDA, there will be no train/test split...I'm just trying to get a sense of what variables 'matter most' across the dataset
- I will use XGBoost first.
- This is often my preferred model.
- However, I'm assuming that the client actually needs an interpretable linear-regression-esque equation
- Thus we cannot end with XGBoost
- The nature of the tree splits is to be less interpretable
- We know that something matters, but we wont' have a coefficient or a direction (positive/negative)
- 2. Modeling
- So after deploying XGBoost for selection, I will run Linear Regressions next for the official modeling
- This will help us with our ability to understand our features.
- This regression will NOT be scaled.
- Why?
- In a word, interpretability
- Once I scale the data, especially if a normalize it, then I lose the units
- To be fair, this data does not seem to have units, but I want to make this as real as possible
- Before we lose the units through scaling and then apply the regularization via Ridge or Lasso, we want to be sure we care more about prediction than inference
- Given how this problem is stated, I'm assuming it's an odd combination of 'interpretability prediction', which for here would have to be an interpretable linear regression equation
- 1. Exploratory Data Analysis
👇🏾 Question 2 - Part I.B - 1. Exploratory Data Analysis 👇🏾¶¶
In [9]:
#check correlations
#all credit due to Pedro Marcelino; https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python
corrmat = ds1.corr()
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True)
plt.savefig('ds1_correlation_matrix.png',dpi=300, bbox_inches='tight')
💭👆🏾Just a few thoughts...👆🏾💭¶
- We can already see that 'dependent' variable 'yc' doesn't seem to be correlated with any of the independent features.
- We also see that the the 'dependent' variables seem to correlate to one another. This may be be because they are somehow generated from the same independent variables.
- This is also going to come up later on, but we can see that x1 and x3 are correlated.
- We can begin to say:
In [10]:
#check the pair plot
#all credit due to Pedro Marcelino; https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python
sns.pairplot(ds1, height = 2.5)
plt.savefig('ds1_pair_plot_lower_res.png',dpi=150, bbox_inches='tight')
plt.show()
💭👆🏾Just a few thoughts...👆🏾💭¶
- As we scan the chart, we can see a few linear-ish or nearly linear relationships.
- These charts don't include interactions.
- We see some blobs, and some bizarre relationships.
- We will eliminate features in the feature selection process.
👇🏾 Using XGBoost & Regressions To See Which Features Seem Important Before Modeling👇🏾¶¶
- Now we shall use XGBoost to get an idea for the features that are most important.
- Before checking the polynomial relationships, we just want to get an idea of what's going on.
- After each XGBoost output we'll also put all the features into a regression against the given dependent variable and see how closely the regression follows the XGBoost insights.
- First, we will start by defining our dependent and independent dataframes.
In [11]:
#define the dependent and independent relationships
features = ds1[['x1','x2','x3','x5','x6']]
ya = ds1['ya']
yb = ds1['yb']
yc = ds1['yc']
🔬 Function 2: XGBoost Modeler 🔬¶
- This function will check all features against the dependent variables
- It will then plot the predicteds vs. actuals as well as the residuals.
- It will provide pertinent model metrics.
- It will then make a horizontal bar chart for feature importance.
In [12]:
def xgboost_modeler(X, y, cv_folds_int, x_left, x_right, y_bottom, y_top):
X_train_xgb, X_test, y_train_xgb, y_test = train_test_split\
(X, y, test_size = 0.2, random_state = 1)
xgb_model = xgb.XGBRegressor(objective="reg:squarederror", random_state=42)
xgb_model.fit(X_train_xgb,y_train_xgb)
y_pred_xgb = xgb_model.predict(X_test)
#linear regression standard
#print r^2 and plot linear regression
print("R^2:"+str(round(r2_score(y_test, y_pred_xgb),4)))
print("Mean Absolute Error (in original units): "+str(round(mean_absolute_error((y_test),\
(y_pred_xgb)),2)))
print("Median Absolute Error (in original units): "+str(round(median_absolute_error((y_test),\
(y_pred_xgb)),2)))
print("STD of Training Data (in original units): "+str(round(np.std((y_test)),2)))
print("Mean of Training Data (in original units): "+str(round(np.mean((y_test)),2)))
print("Residual Skew: "+str(round(stats.skew((y_pred_xgb)-(y_test)),2)))
print("Residual Kurtosis: "+str(round(stats.kurtosis((y_pred_xgb)-(y_test)),2)))
fig, axs = plt.subplots(ncols=(3), nrows=(1), figsize=(15,4))
plt.subplots_adjust(top = 0.95, bottom=0.01, hspace=0.4, wspace=0.16)
axs[0].scatter((y_test), (y_pred_xgb))
axs[0].set_title('Actuals (x) vs Predicteds (y)',fontweight = 'bold')
axs[0].set_xlim([x_left, x_right])
axs[0].set_ylim([y_bottom, y_top])
axs[1].scatter((y_test), ((y_pred_xgb)-(y_test)))
axs[1].set_title('Actuals (x) vs Residuals [pred. - actual] (y)',fontweight = 'bold')
axs[1].set_xlim([x_left, x_right])
axs[1].set_ylim([y_bottom, y_top])
axs[2].hist(((y_pred_xgb)-(y_test)), bins=100)
axs[2].set_title('Histogram of Residuals [predicted - actual]',fontweight = 'bold')
plt.savefig('xgboost_modeling_',dpi=300, bbox_inches='tight')
plt.show()
print("Cross-Validation Scoring for XGBoost")
print('Mean Mean Absolute Error: {}'.format(-1*round(cross_val_score(xgb_model, X_train_xgb, y_train_xgb, \
cv=cv_folds_int, scoring='neg_mean_absolute_error').mean(),2)))
print('Median Mean Absolute Error: {}'.format(-1*round(np.median(cross_val_score(xgb_model, X_train_xgb,\
y_train_xgb, \
cv=cv_folds_int, scoring='neg_mean_absolute_error')),2)))
print('Mean R^2: {}'.format(round(cross_val_score(xgb_model, X_train_xgb, y_train_xgb, cv=cv_folds_int,\
scoring='r2').mean(),3)))
print('Median R^2: {}'.format(round(np.median(cross_val_score(xgb_model,\
X_train_xgb, y_train_xgb, cv=cv_folds_int, scoring='r2')),3)))
print('Max R^2: {}'.format(round(np.max(cross_val_score(xgb_model,\
X_train_xgb, y_train_xgb, cv=cv_folds_int, scoring='r2')),3)))
print('Min R^2: {}'.format(round(np.min(cross_val_score(xgb_model,\
X_train_xgb, y_train_xgb, cv=cv_folds_int, scoring='r2')),3)))
top_xgb_features = pd.DataFrame(sorted(list(zip(X,xgb_model.feature_importances_))\
,key = lambda x: abs(x[1]),reverse=True)[:10], columns=['Feature', 'XGBoost Importance'])
top_xgb_features
bar_count = range(len(top_xgb_features.Feature))
fig, axs = plt.subplots(ncols=2, figsize=(14,4))
#using a subplot method coupled with an inline parameter to have high resolution
#note: "[::-1]" reverses the column in a pandas dataframe
axs[1].set_axis_off()
axs[0].barh(bar_count, top_xgb_features['XGBoost Importance'][::-1],\
align='center', alpha=1)
axs[0].set_xlabel('Values')
axs[0].set_yticks(bar_count)
axs[0].set_yticklabels(top_xgb_features.Feature[::-1], fontsize=10)
axs[0].set_xlabel('XGBoost Importance')
axs[0].set_title("XGBoost's Feature Importances",fontweight = 'bold')
extent = axs[0].get_window_extent().transformed(fig.dpi_scale_trans.inverted())
fig.savefig('xgb_features_importance_dependent_var_',dpi=300, bbox_inches=extent.expanded(1.5, 1.5))
plt.show()
👇🏾 First, we check dependent variable 'ya'👇🏾¶¶
In [13]:
xgboost_modeler(X=features, y=ya, cv_folds_int=4, x_left=-100, x_right=100, y_bottom=-100, y_top=100)
💭👆🏾Just a few thoughts...👆🏾💭¶
- This lets us know that for dependent variable ya it looks like:
- x5 and x6 are **not** important
- x2 is most important, followed by x1 and x3
- We'll also peek at a standard regression...
In [14]:
features_regression = add_constant(features)
regressor_OLS = sm.OLS(endog = ya, exog = features_regression).fit()
regressor_OLS.summary()
Out[14]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- This regression confirms our XGBoost insights regarding dependent variable ya.
👇🏾 Second, we check dependent variable 'yb'👇🏾¶¶
In [15]:
xgboost_modeler(X=features, y=yb, cv_folds_int=4, x_left=-5, x_right=5, y_bottom=-5, y_top=5)
💭👆🏾Just a few thoughts...👆🏾💭¶
- This is amazing!!!
- It looks like just one variable, **x1** decides what dependent variable 'yb' is going to be.
- We'll check this with a standard regression...
In [16]:
features_regression = add_constant(features)
regressor_OLS = sm.OLS(endog = yb, exog = features_regression).fit()
regressor_OLS.summary()
Out[16]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- Even though you could argue that x3 is also statistically significant, you can see that x1's t-score of over 500 shows that it is likely doing most of the work here, just like XGBoost said.
- Let's actually run this regression again...but only for x1 and see if we keep the R^2 at about the same level (meaning that the variance can be accounted for by x1).
💭👇🏾Let's also check that one feature, x1, which seems to be doing the brunt of the work here.👇🏾💭¶
In [17]:
features_regression_x1 = add_constant(ds1.x1)
regressor_OLS = sm.OLS(endog = yb, exog = features_regression_x1).fit()
regressor_OLS.summary()
Out[17]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- Yup, the R^2 didn't budge. It stayed right there at 0.895.
- Looks like we have the x1 independent variable predicting the yb dependent variable.
👇🏾 Third, we check dependent variable 'yc'👇🏾¶¶
In [18]:
xgboost_modeler(X=features, y=yc, cv_folds_int=4, x_left=-0.1, x_right=0.1, y_bottom=-0.1, y_top=0.1)
💭👆🏾Just a few thoughts...👆🏾💭¶
- We know from the correlation chart above that yc was going to be trouble because it did not show any correlation to any of the x's.
- Here we see more of the same. The R^2 in the test case was literally negative. Meaning the model is actually worse than guessing.
- This is a tell-tale sign that yc has nothing to do with the x's.
- We'll double-check this with a regression, and we'll give polynomial feature engineering a shot later on.
In [19]:
features_regression = add_constant(features)
regressor_OLS = sm.OLS(endog = yc, exog = features_regression).fit()
regressor_OLS.summary()
Out[19]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- Again, we see an R^2 of 0.
- And p-values don't seem too impressive either. (There are two around the 5% p-value level, but given everything else we've seen, I don't think yc relates to the given independent features)
- ...and if this were Hayward Data Science (my consulting business), I would not model this because we might get an even worse outcome then not knowing -- thinking we know when we are actually wrong.
💭👇🏾Brief 'Side Path' - Checking for Polynomial Features👇🏾💭¶
- We've already discussed that we are not going to perform log transforms due to the negative numbers and the lack of units (such that we don't know if adding a constant will hurt our interpretability even though it wouldn't affect variance.
- What we will want to do is train our XGBoost models on the polynomial features to see if there is marked improvement.
- Just for demonstration purposes I will also run a Backwards Elimination Regression Algorithm as well.
In [20]:
def backwardElimination(X, y, significance_level_float):
numCols = len(X.columns)
for i in range(numCols):
regressor_OLS = sm.OLS(y, X).fit()
maxVar = max(regressor_OLS.pvalues)
if maxVar > significance_level_float:
for j in range(len(regressor_OLS.pvalues)):
if (regressor_OLS.pvalues[j] == maxVar):
X = X.drop(X.columns[j], axis = 1)
regressor_OLS.summary()
return X
In [21]:
#make the polynomial features
poly = PolynomialFeatures(degree=3, include_bias=False)
features_w_poly = poly.fit_transform(features)
features_w_poly = pd.DataFrame(features_w_poly, columns=poly.get_feature_names(features.columns))
In [22]:
#adding in a constant to prepare for regression (we had bias as False above so we need to add the constant now)
features_poly_regression = add_constant(features_w_poly)
In [23]:
#XGboost polynomial test for ya
xgboost_modeler(X=features_w_poly, y=ya, cv_folds_int=4, x_left=-100, x_right=100, y_bottom=-100, y_top=100)
💭👆🏾Just a few thoughts...👆🏾💭¶
- Even with the new polynomial features, we still are focused on x1, and x2.
- Plus the R^2 actually dropped from 0.704 to 0703.
- We've added a ton of features, but we don't see a marked improvement.
- All things being equal, we want a simpler model, so I'll be sticking to the features we had when we move to the modeling stage for ya.
In [24]:
#backwards elimination regression algo for ya
#note this cell runs the Backwards Elimination, and this cell will takw a while to execute
Features_Poly_Pruned_ya = backwardElimination(features_poly_regression, ya, significance_level_float = 0.00001)
Features_Poly_Pruned_ya = add_constant(Features_Poly_Pruned_ya)
pruned_regressor_OLS = sm.OLS(endog = ya, exog = Features_Poly_Pruned_ya).fit()
pruned_regressor_OLS.summary()
Out[24]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- The regression here seems to generally map on to what I just said about the XGBoost run above.
- We added a lot of features, but we don't see much uptick in explanation of variance.
- I'm inclined to stick with my simpler features.
In [25]:
#XGboost polynomial test for yb
xgboost_modeler(X=features_w_poly, y=yb, cv_folds_int=4, x_left=-6, x_right=6, y_bottom=-6, y_top=6)
💭👆🏾Just a few thoughts...👆🏾💭¶
- This is truly amazing!!!!!!!!
- No matter what we throw at yb, it just wants to stick with x1.
- It's yb is the faithful, monogamous type 😊haahahah!
In [26]:
#backwards elimination regression algo for yb
#note this cell runs the Backwards Elimination, and this cell will takw a while to execute
Features_Poly_Pruned_yb = backwardElimination(features_poly_regression, yb, significance_level_float = 0.00001)
Features_Poly_Pruned_yb = add_constant(Features_Poly_Pruned_yb)
pruned_regressor_OLS = sm.OLS(endog = yb, exog = Features_Poly_Pruned_yb).fit()
pruned_regressor_OLS.summary()
Out[26]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- Again, it's all x1 and variations thereof.
- We're defintiey just going to stick to x1 for the modeling.
- ❤️ x1 and yb are a true data love story. ❤️
In [27]:
#XGboost polynomial test for yc
xgboost_modeler(X=features_w_poly, y=yc, cv_folds_int=4, x_left=-0.1, x_right=0.1, y_bottom=-0.1, y_top=0.1)
💭👆🏾Just a few thoughts...👆🏾💭¶
- Again, it's not looking good for modeling yc.
- At first blust it looks like XGBoost may have finally found something, but:
- We see that we are essentially now interacting x1, x5 and x6 together in one mega column, and we still are only at about 0.10 gain.
- Since we've seen such limited correlation in our EDA phase, I'm inclined to be quite skeptical here and say that if you throw enough things against the wall, then something is bound to stick.
- But we are trying to follow the data, not jam the data into our goals.
- I would be inclined to not model yc at all here.
In [28]:
#backwards elimination regression algo for yc
#note this cell runs the Backwards Elimination, and this cell will takw a while to execute
Features_Poly_Pruned_yc = backwardElimination(features_poly_regression, yc, significance_level_float = 0.00001)
Features_Poly_Pruned_yc = add_constant(Features_Poly_Pruned_yc)
pruned_regressor_OLS = sm.OLS(endog = yc, exog = Features_Poly_Pruned_yc).fit()
pruned_regressor_OLS.summary()
Out[28]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- And here we have, as my late Dad would say, 'more trouble in River City.'
- Namely, in the Backwards Elimination Regression Round none of the independent variables were able to meet the threshold level needed to say in the regression.
- The data is screaming at us not to model yc.
👇🏾 Question 2 - Part I.B - 2. Modeling 👇🏾¶¶
🔬 Function 3: Multiple Linear Regression Modeler 🔬¶
In [29]:
def multiple_linear_regression_modeler(X, y, cv_folds_int, x_left, x_right, y_bottom, y_top):
X_train1, X_test, y_train1, y_test = train_test_split\
(X, y, test_size = 0.2, random_state = 1)
model = LinearRegression()
model.fit(X_train1,y_train1)
y_slr_predict = model.predict(X_test)
#linear regression standard
#print r^2 and plot linear regression
print("R^2:"+str(round(r2_score(y_test, y_slr_predict),4)))
print("Mean Absolute Error (in original units): "+str(round(mean_absolute_error((y_test),\
(y_slr_predict)),2)))
print("STD of Training Data (in original units): "+str(round(np.std((y_test)),2)))
print("Mean of Training Data (in original units): "+str(round(np.mean((y_test)),2)))
print("Residual Skew: "+str(round(stats.skew((y_slr_predict)-(y_test)),2)))
print("Residual Kurtosis: "+str(round(stats.kurtosis((y_slr_predict)-(y_test)),2)))
fig, axs = plt.subplots(ncols=(3), nrows=(1), figsize=(15,4))
plt.subplots_adjust(top = 0.95, bottom=0.01, hspace=0.25, wspace=0.16)
axs[0].scatter((y_test), (y_slr_predict))
axs[0].set_title('Actuals (x) vs Predicteds (y)',fontweight = 'bold')
axs[0].set_xlim([x_left, x_right])
axs[0].set_ylim([y_bottom, y_top])
axs[1].scatter((y_test), ((y_slr_predict)-(y_test)))
axs[1].set_title('Actuals (x) vs Residuals [pred. - actual] (y)',fontweight = 'bold')
axs[1].set_xlim([x_left, x_right])
axs[1].set_ylim([y_bottom, y_top])
axs[2].hist(((y_slr_predict)-(y_test)), bins=100)
axs[2].set_title('Histogram of Residuals [predicted - actual]',fontweight = 'bold')
plt.savefig('linear_regression_modeling_',dpi=300, bbox_inches='tight')
plt.show()
print("Cross-Validation Scoring for Standard Linear Regression")
print('Mean Absolute Error: {}'.format(-1*round(cross_val_score(model, X_train1, y_train1, \
cv=cv_folds_int, scoring='neg_mean_absolute_error').mean(),2)))
print('Median Mean Absolute Error: {}'.format(-1*round(np.median(cross_val_score(model, X_train1, y_train1, \
cv=cv_folds_int, scoring='neg_mean_absolute_error')),2)))
print('Mean R^2: {}'.format(round(cross_val_score(model, X_train1, y_train1, cv=cv_folds_int, \
scoring='r2').mean(),2)))
print('Median R^2: {}'.format(round(np.median(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
print('Max R^2: {}'.format(round(np.max(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
print('Min R^2: {}'.format(round(np.min(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
linear_regression_feature_list = []
linear_regression_coef_list = []
for i in X.columns:
linear_regression_feature_list.append(i)
for i in range(len(model.coef_)):
linear_regression_coef_list.append(model.coef_[i])
top_10_linear_regression_features = pd.DataFrame(sorted(list(zip(linear_regression_feature_list,\
linear_regression_coef_list))\
,key = lambda x: abs(x[1]),reverse=True)[:10], columns=['Feature', 'Linear Regression Coefficient'])
top_10_linear_regression_features
print(top_10_linear_regression_features)
print('Intercept: {}'.format(model.intercept_))
🔬 Function 4: Single Linear Regression Modeler 🔬¶
In [30]:
def single_linear_regression_modeler(X, single_feature_str, y, cv_folds_int, x_left, x_right, y_bottom, y_top):
X_train1, X_test, y_train1, y_test = train_test_split(X, y, test_size = 0.2, random_state = 1)
model = LinearRegression()
model.fit(X_train1,y_train1)
y_slr_predict = model.predict(X_test)
print("R^2:"+str(round(r2_score(y_test, y_slr_predict),4)))
print("Mean Absolute Error (in original units): "+str(round(mean_absolute_error((y_test),\
(y_slr_predict)),2)))
print("STD of Training Data (in original units): "+str(round(np.std((y_test)),2)))
print("Mean of Training Data (in original units): "+str(round(np.mean((y_test)),2)))
print("Residual Skew: "+str(round(stats.skew((y_slr_predict)-(y_test)),2)))
print("Residual Kurtosis: "+str(round(stats.kurtosis((y_slr_predict)-(y_test)),2)))
fig, axs = plt.subplots(ncols=(3), nrows=(1), figsize=(15,4))
plt.subplots_adjust(top = 0.95, bottom=0.01, hspace=0.25, wspace=0.16)
axs[0].scatter((y_test), (y_slr_predict))
axs[0].set_title('Actuals (x) vs Predicteds (y)',fontweight = 'bold')
axs[0].set_xlim([x_left, x_right])
axs[0].set_ylim([y_bottom, y_top])
axs[1].scatter((y_test), ((y_slr_predict)-(y_test)))
axs[1].set_title('Actuals (x) vs Residuals [pred. - actual] (y)',fontweight = 'bold')
axs[1].set_xlim([x_left, x_right])
axs[1].set_ylim([y_bottom, y_top])
axs[2].hist(((y_slr_predict)-(y_test)), bins=100)
axs[2].set_title('Histogram of Residuals [predicted - actual]',fontweight = 'bold')
plt.savefig('linear_regression_modeling_',dpi=300, bbox_inches='tight')
plt.show()
print("Cross-Validation Scoring for Standard Linear Regression")
print('Mean Absolute Error: {}'.format(-1*round(cross_val_score(model, X_train1, y_train1, \
cv=cv_folds_int, scoring='neg_mean_absolute_error').mean(),2)))
print('Median Mean Absolute Error: {}'.format(-1*round(np.median(cross_val_score(model, X_train1, y_train1, \
cv=cv_folds_int, scoring='neg_mean_absolute_error')),2)))
print('Mean R^2: {}'.format(round(cross_val_score(model, X_train1, y_train1, cv=cv_folds_int, \
scoring='r2').mean(),2)))
print('Median R^2: {}'.format(round(np.median(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
print('Max R^2: {}'.format(round(np.max(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
print('Min R^2: {}'.format(round(np.min(cross_val_score(model,\
X_train1, y_train1, cv=cv_folds_int, scoring='r2')),3)))
linear_regression_feature_list = []
linear_regression_coef_list = []
for i in list([single_feature_str]):
linear_regression_feature_list.append(i)
for i in range(len(model.coef_)):
linear_regression_coef_list.append(model.coef_[i])
top_10_linear_regression_features = pd.DataFrame(sorted(list(zip(linear_regression_feature_list,\
linear_regression_coef_list))\
,key = lambda x: abs(x[1]),reverse=True)[:10], columns=['Feature', 'Linear Regression Coefficient'])
top_10_linear_regression_features
print(top_10_linear_regression_features)
print('Intercept: {}'.format(model.intercept_))
👇🏾 First, let's build a multiple linear regression model for dependent variable 'ya.'👇🏾¶
In [31]:
multiple_linear_regression_modeler(X=ds1[['x1','x2','x3']],
y=ya, cv_folds_int=4, x_left=-100, x_right=100, y_bottom=-100,y_top=100)
👇🏾 Conclusion: Equation for Predicting 'ya.'👇🏾¶
$ \begin{align} ya = 0.251(x1) + 6.236(x2) + 0.956(x3) + 7.58 \end{align} $
👇🏾 Second, let's build a single linear regression model for dependent variable 'yb.'👇🏾¶¶
In [32]:
single_linear_regression_modeler(X=ds1.x1.values.reshape(-1, 1), single_feature_str='x1',\
y=yb, cv_folds_int=4, x_left=-6, x_right=6, y_bottom=-6, y_top=6)
👇🏾 Conclusion: Equation for Predicting 'yb.' 👇🏾¶
$ \begin{align} yb = 0.253(x1) + 0.8422 \end{align} $
👇🏾 Finally, let's build a multiple linear regression model for dependent variable 'yc.' We don't expect it to work 👇🏾¶
In [33]:
multiple_linear_regression_modeler(X=ds1,y=yc, cv_folds_int=4, x_left=-0.1, x_right=0.1, y_bottom=-0.1, y_top=0.1)
👇🏾 Conclusion: No Equation for Predicting 'yc.' 👇🏾¶
Part II. - Unsupervised Machine Learning ¶
- Consider data set 2 (ds2.csv). The data asset comprises a set of observations.
👇🏾 Just as we did in Part 1, Question 1, we begin by reading in the data. 👇🏾¶
- We will check for NULLS.
- For some reason there is an index column in the raw data, but a Pandas dataframe comes with an index column by default so we will drop this 'extra' index column in the raw data.
- We will summarize the data.
- We will also check the distributions.
In [34]:
ds2 = pd.read_csv('ds2.csv')
In [35]:
#drop the extra index column since a Pandas dataframe has one built in
ds2 = ds2.drop(ds2.columns[0], axis=1)
In [36]:
#peek at the data
ds2.head()
Out[36]:
In [37]:
#check for NULLS
msno.matrix(ds2, color = (.0,.0,.2))
Out[37]:
👇🏾 Just as we did in Part 1, Question 1, we begin by reading in the data. 👇🏾¶
In [38]:
#let's describe the data
ds2.describe()
Out[38]:
In [39]:
#all credit due to: Pedro Marcelino; https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python
#correlation matrix
corrmat = ds2.corr()
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True)
plt.savefig('ds1_correlation_matrix.png',dpi=300, bbox_inches='tight')
💭👆🏾Just a few thoughts...👆🏾💭¶
In [40]:
#all credit due to: Pedro Marcelino, https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python
sns.pairplot(ds2, size = 2.5)
plt.show()
In [41]:
#lets also check the distributions again
for i in ds2.columns:
distribution_visualizer_and_fitter(ds2, i)
💭👆🏾Just a few thoughts...👆🏾💭¶
- One thing that immediately jumps out at you is how bi-modal the data is.
Question 4 - Part II.B ¶
- How would you visualize this data set? Can you make an interesting visualization?
👇🏾 Gameplan for Principal Component Analysis (PCA)👇🏾¶
- The difficult with visualizing data like this, or doing any type of serious unsupervised learning is the high dimensionality.
- As I read somewhere on StackOverflow once, 'do you know anyone who can see in 10D?" 😊
- So what we need to is figure out how to reduce the dimensionality.
- We want to see how few components we need to approximate the cumulative variance.
- Once we reduce the components a bit, we'll be able to tease out insight.
- Once we tease out the insight, we'll be able to visualize it.
- And, not to give away the next section, but then we can begin to cluster it with K-Means.
👇🏾 First thing's first: how many PCA components do we actually need?👇🏾¶
- The first thing we do is figure out how many components we need to explain most of the variance in this high dimensional dataset.
- I'll show you two methods.
- One through a standard PCA method. Special thanks to Bartosz Mikulski for blogging about it.
- One through a special Python third-party package. Special thanks to Andrew Heusser for writing it.
In [42]:
#it's important to scale the data so you can compare apples-to-apples how much each component contributes to variance
scaler = MinMaxScaler()
data_rescaled = scaler.fit_transform(ds2)
In [43]:
#we'll need to set a variance threshold. I choose 95%; some others my choose 99%; it's up to you.
pca = PCA(n_components = 0.95)
pca.fit(data_rescaled)
reduced = pca.transform(data_rescaled)
In [44]:
#all credit here to Bartosz Mikulski, https://www.mikulskibartosz.name/pca-how-to-choose-the-number-of-components/
pca = PCA().fit(data_rescaled)
plt.rcParams["figure.figsize"] = (12,6)
fig, ax = plt.subplots()
xi = np.arange(1, 11, step=1)
y = np.cumsum(pca.explained_variance_ratio_)
plt.ylim(0.0,1.1)
plt.plot(xi, y, marker='o', linestyle='--', color='b')
plt.xlabel('Number of Components')
plt.xticks(np.arange(0, 11, step=1)) #change from 0-based array index to 1-based human-readable label
plt.ylabel('Cumulative variance (%)')
plt.title('The number of components needed to explain variance')
plt.axhline(y=0.95, color='r', linestyle='-')
plt.text(0.5, 0.85, '95% cut-off threshold', color = 'red', fontsize=16)
ax.grid(axis='x')
plt.show()
💭👆🏾Just a few thoughts...👆🏾💭¶
- Reading a chart like this, I'd say is both an art and a science.
- You are sort of looking for the 'elbow' point.
- As we look above we can see that after 3 components we've already gotten to above 80% of the cumulative variance.
- That's a major win for us because it reduces our data set to just 3 components.
- To get to the 95% cumulative variance level we'd need 7 components.
- But that doesn't really help us so much because 7 components is still a lot to interpret. If we, say, reduced a 3000 component set to 7 then that would be saying something. But if we take 10 components and make it 7, in my opinion we haven't really helped ourselves.
- Thus, I will go with 3 components as the number for our reduction.
In [45]:
#there is also a third-party Python plug-in called 'Hypertools' which has a similiar procedure
#it looks to correlation
#the package was written by Andrew Heusser,
#https://hypertools.readthedocs.io/en/latest/auto_examples/plot_describe.html#sphx-glr-auto-examples-plot-describe-py
hyp.describe(ds2)
Out[45]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- So, we see that 3 components seems to be the magic elbow point.
- Let's see what happens when we graph these 3 components....
👇🏾 Second thing's second: let's look at our data in 3D, as opposed to 10D.👇🏾¶
- Now we will take a look at our first visualization.
- We will set our PCA algorithm to have n_components of 3 after all the work we just did to arrive at that number.
- It appears that 3 really is the magic number.
In [46]:
pca_execution = PCA(n_components=3)
pca_doing_its_work = pca_execution.fit_transform(data_rescaled)
print("original shape: ", data_rescaled.shape)
print("transformed shape:", pca_doing_its_work.shape)
In [47]:
pca_conjecture = pd.DataFrame(pca_doing_its_work, columns=['component_1','component_2', 'component_3'])
In [48]:
pca_conjecture.head()
Out[48]:
In [49]:
#all credit here to Jake Vanderplas,
#https://jakevdp.github.io/PythonDataScienceHandbook/04.12-three-dimensional-plotting.html
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter3D(pca_conjecture.component_1, pca_conjecture.component_2, pca_conjecture.component_3)
plt.savefig('pca_3D_graph.png',dpi=300, bbox_inches='tight')
💭👆🏾Just a few thoughts...👆🏾💭¶
- This is a really cool 3D plot of what how we can group the data.
- Although there are 10(!) dimensions, we can reduce those to basically 3 dimensions.
- Very, very conveniently, the human mind can picture 3 dimensions on a 3D plot.
- It's very possible to have reduced to components which still would not be easily interpreted, so I think we really lucked out here! That's why I say 'The Lord Always Delivers.' 😊🙏
- Key point: It's tempting to think that this means that all the observations (i.e., the rows) can be split into 3 fields. But we're not there just yet.
- We have basically reduced the columns in this gigantic data set from 10 to 3, but we now will want to figure out how to apply K-Means clustering to find out more.
- Note: K-Means begins to break down at high dimensions because of what's called 'the curve of dimensionality,'
- Basically, anything that uses Euclidean distance is going to get a lot less meaningful because dimensions will begin to make everything potentially close to everything else when you combine it all to one vector. Therefore, it helps to reduce the components first.
- Also, please note: There are other clustering methods that can perform better with high dimensionality, but that's outside the scope of this notebook. If you are curious, more info here.
- Basically, anything that uses Euclidean distance is going to get a lot less meaningful because dimensions will begin to make everything potentially close to everything else when you combine it all to one vector. Therefore, it helps to reduce the components first.
Question 5 - Part II.C ¶
- Someone suggests that the observations are really from multiple different files and were accidentally joined into one larger data set. Does anything about the data set suggest this? If so, how many different sources/file do you think there are?
👇🏾 We begin our K-Means Clustering analysis now we've completed the PCA round.👇🏾¶
- Let the K-Means games begin!
- The first rule of thumb is to use an algorithm to choose you n clusters for K means.
- It's tempting to just jump right in and use the number 3 from before, but we will want to evaluate an **'elbow'** graph again.
- By the way, in case you are curious, I actually tested K-means on the PCA-generated dataframe framework offline, and 3 was **not** the optimum number of K clusters.
In [50]:
#all credit to Dmitriy Kavyazin, http://bit.ly/2mYDWv7
ks = range(1, 10)
inertias = []
for k in ks:
# Create a KMeans instance with k clusters: model
model = KMeans(n_clusters=k)
# Fit model to samples
model.fit(pca_conjecture) #.iloc[:,:3])
# Append the inertia to the list of inertias
inertias.append(model.inertia_)
plt.plot(ks, inertias, '-o', color='black')
plt.xlabel('number of clusters, k')
plt.ylabel('inertia')
plt.xticks(ks)
plt.show()
💭👆🏾Just a few thoughts...👆🏾💭¶
- So we can actually see the elbow point comes at **4** clusters!!! ✅🙏
- Let's go ahead and plot that next!!
In [51]:
#here comes K-Means at long last!
kmeans = KMeans(n_clusters=4)
kmeans.fit(pca_conjecture)
choice_kmeans = kmeans.predict(pca_conjecture)
In [52]:
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter3D(pca_conjecture.component_1, pca_conjecture.component_2, pca_conjecture.component_3,\
c=choice_kmeans, cmap='Dark2')
plt.savefig('kmeans_on_pca.png',dpi=300, bbox_inches='tight')
💭👆🏾Just a few thoughts...👆🏾💭¶
- And voilà: now we have a great visualization of the 10D data.
- We can see that there are basically **3 principal components (the blobs)** and **4 cohorts (the colors)** the actual observations can be sorted into.
- In future analysis, we'll want to study why 2 cohorts seem to get jammed into 1 component, yet they are not themselves considered the same cohort.
- We'll want to better understand what factors K-Means() is using to differentiate them.
- In future analysis, we'll want to study why 2 cohorts seem to get jammed into 1 component, yet they are not themselves considered the same cohort.
Question 6 - Part II.D - Bonus ¶
- Bonus points: If you think there are more than one source in ds2, can you assign each observation to the right source (based on the number of sources you identified in 2c)?
👇🏾 Applying our K-means work from the last question...👇🏾¶
- We can simply take the same K-Means classifications we used to make the hues in the 3D plot above, and apply them to the ds2 dataframe.
In [53]:
#just to be safe, we can check the shape so we know that they'll match
print(ds2.shape)
print(choice_kmeans.shape)
In [54]:
ds2['cohort_guess'] = choice_kmeans
In [55]:
ds2.sample(10)
Out[55]:
💭👆🏾Just a few thoughts...👆🏾💭¶
- And there we have it!
- I took a sample of the dataframe instead of the head() method, so you can get an idea of the different 'cohort guesses'.
- You can see that there are 4 different K-Means cohorts now.
- 0, 1, 2, and 3
- I would argue that these are where the disparate files came from before they were joined in ds2.
Part III. - SQL Coding ¶
- Stack Overflow provides a tool at https://data.stackexchange.com/stackoverflow/query/new that allows SQL queries to be run against their data. After reviewing the database schema provided on their site, please answer the questions below by providing both your answer and the query used to derive it.
👇🏾 Full answers and clickable/runable results can be found on the Stack Overflows query permalink.👇🏾¶¶
http://bit.ly/2mJ5nct¶
In [56]:
query_3_a_results = pd.read_csv('stackoverflow_query_3-a_results.csv')
query_3_a_results
Out[56]:
Question 8 - Part III.B ¶
- What is the distribution of views for users?
- Stack Overflow Query Permalink:
In [57]:
query_3_b_results = pd.read_csv('stackoverflow_query_3-b_results.csv')
query_3_b_results.head(10)
Out[57]:
In [58]:
ax = sns.regplot(x="View_Level", y="Frequency", data=query_3_b_results, fit_reg=False)
ax.set_ylabel('Frequency', fontweight = 'bold')
ax.set_xlim(0,40)
ax.set_xlabel("View Level", fontweight = 'bold')
ax.yaxis.set_major_formatter(mpl.ticker.StrMethodFormatter('{x:,.0f}'))
ax.set_title("Histogram of Views Generated by Above SQL Query", fontweight = 'bold')
plt.show()
In [59]:
plt.figure(figsize=(16,5))
plt.bar(query_3_b_results.View_Level, query_3_b_results.Frequency)
plt.ylim(0,300000)
plt.xlim(0,80)
plt.title("Zoomed-In Snapshot of Stack Overflow Views Distribution", fontweight='bold')
plt.ylabel("Frequency", fontweight = 'bold')
plt.xlabel("View Level", fontweight = 'bold')
plt.show()
Question 9 - Part III.C ¶
- What are the top 10 tags associated with posts created in 2017?
- Stack Overflow Query Permalink:
In [60]:
query_3_c_results = pd.read_csv('stackoverflow_query_3-c_results.csv')
query_3_c_results
Out[60]:
Question 10 - Part III.D - Bonus ¶
- Bonus points: For the questions created in 2017, what was the average time (in seconds) between when the question was created and when the accepted answer was provided?
- Stack Overflow Query Permalink:
In [61]:
query_3_d_bonus_results = pd.read_csv('stackoverflow_query_3-d_results.csv')
query_3_d_bonus_results
Out[61]:
Thank you so much for your time and consideration!¶
Best,¶
George John Jordan Thomas Aquinas Hayward, Optimist
