<jupyter_start><jupyter_text>Final Project Submission: Regression Models - KC Housing DataPlease fill out:* Student name: Princess '<NAME>* Student pace: Self Paced* Scheduled project review date/time: March 09, 2021 / 5:00PM* Instructor name: <NAME>* Blog post URL: [Visualizations for Data Science](https://seunotusanya.medium.com/visualization-for-data-science-ae6cf75c9514)  ObjectiveIn this Jupyter notebook we will finally be working on Linear Regression and finding the best fit model for our business case. The order of business is as follows: - Linear Regression- Feature Selection- Train-Test Split- Cross Validation Import neccesary Libraries and Load Data<jupyter_code># Let's start by importing nessesary libraries
import pandas as pd
import numpy as np

# For visualizations
import seaborn as sns  
import matplotlib.pyplot as plt
%matplotlib inline

# For regression
import statsmodels.api as sm
import scipy.stats as stats
import statsmodels.formula.api as smf
from sklearn.feature_selection import RFE
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler
# Import new and improved data and drop extra default index
kc_reg = pd.read_csv('Data/kc_reg.csv')
kc_reg.drop(['Unnamed: 0'], axis=1, inplace=True)
# Print top 5 rows to ensure data loads properly
kc_reg.head()<jupyter_output><empty_output><jupyter_text>RecapA couple of things to remember from the previous notebook before getting started on Regression:* No variable is 'perfectly' normal but 'sqft_living' comes pretty close.* There are 2 predictors that are highly Correlated: 'sqft_living' and 'sqft_above'.* Dummy variables have been created for _Norminal_ categorical variables.  Linear RegressionLet's perform a series of simple regressions between the target variable and some predictors. Step 1: Check DistributionWe already checked for this in the EDA notebook and if we recall the results were:- No variable is "perfectly" normal, but some do tend to follow an overall normal pattern, i.e, 'price', 'sqft_living', 'sqft_living15'.- The target variable **'Price'** is moderately skewed to the right.- **'Sqft_living' and 'Sqft_living15** are more normally distributed than the other features, this can be due to a high amount of variance in other features.- The rest of the features are worse of in terms of normal distribution.  Step 2: Test for Linearity AssumptionsThere are 3 assumptions for Linear Regression: _Linearity_, _Normality_, and _Homoscedasticity_. We can check for them by performing the following steps:- Run a simple regression in Statsmodel wtih the dependent and independent variables- Plot a prediction line on a scatter plot with the dependent and independent variables- Visualize the error term for variance and heteroscedasticity- Check the normality assumptions by creating a QQ-plot> Let's check to see if we have any practical Linear Regressions.<jupyter_code># Build a formula. 'sqft_living' is the most normally distributed variable so we'll start with that
f = 'price~sqft_living'
# Create a fitted model
model = smf.ols(formula=f, data=kc_reg).fit()
# Get regression model summary
model.summary()
# Create a new Dataframe with the minimum and maximum value of 'sqft_living'
X_new = pd.DataFrame({'sqft_living': [kc_reg.sqft_living.min(), kc_reg.sqft_living.max()]});
# Make predictions on the x values and store them 
preds = model.predict(X_new)
# Plot the observed data and the least squares line
kc_reg.plot(kind='scatter', x='sqft_living', y='price');
plt.plot(X_new, preds, c='red', linewidth=2);
plt.show()
# Visualize error terms for variance and heteroscedasticity
fig = plt.figure(figsize=(15,8))
fig = sm.graphics.plot_regress_exog(model, 'sqft_living', fig=fig)
plt.show()
# Show distribution of residuals
sns.distplot(model.resid)
# Create qq plot to check normality assuption
residuals = model.resid
fig = sm.graphics.qqplot(residuals, dist=stats.norm, line='45', fit=True)<jupyter_output><empty_output><jupyter_text>Insights:- The _R-squared_ value is _0.305_ i.e, _%_ of the variance in the target variable can be explained by 'sqft_living'. - R-squared(Coefficient of determination) also tells us that as a predictor, 'sqft_living' is not a very good fit as only _30%_ of the data points fall within the regression line.- If we _increase_ 'sqft_living' by 1 unit we can expect the price of a house to _increase_ by _158_.- The _P-value_ tells us that we can reject the null hypothesis that the _Coefficient = 0_. In other words a changes in 'sqft_living' are related to changes in 'price'.- A normally distributed data has a Skewness between -0.5 and 0.5 and a Kurtosis of 3. The model tells us that we have a Skewness of 0.61 and a Kurtosis of 3.3 this tell us that the data still has quite a few outliers but is pretty normal overall.<jupyter_code># Build a formula. 
f = 'price~distance_from_bellevue'
# Create a fitted model
model = smf.ols(formula=f, data=kc_reg).fit()
# Get regression model summary
model.summary()
# Create a new Dataframe with the minimum and maximum value of 'distance_from_bellevue'
X_new = pd.DataFrame({'distance_from_bellevue': [kc_reg.distance_from_bellevue.min(),
                                                 kc_reg.distance_from_bellevue.max()]});
# Make predictions on the x values and store them 
preds = model.predict(X_new)
# Plot the observed data and the least squares line
kc_reg.plot(kind='scatter', x='distance_from_bellevue', y='price');
plt.plot(X_new, preds, c='red', linewidth=2);
plt.show()
# Visualize error terms for variance and heteroscedasticity
fig = plt.figure(figsize=(15,8))
fig = sm.graphics.plot_regress_exog(model, 'distance_from_bellevue', fig=fig)
plt.show()
# Show distribution of residuals
sns.distplot(model.resid)
# Create qq plot to check normality assuption
residuals = model.resid
fig = sm.graphics.qqplot(residuals, dist=stats.norm, line='45', fit=True)<jupyter_output><empty_output><jupyter_text>Insights:- The _R-squared_ value is _0.242_ i.e, _24%_ of the variance in the target variable can be explained by 'distance_from_bellevue'. - R-squared(Coefficient of determination) also tells us that as a predictor, 'distance_from_bellevue' is also not a very good fit. - If we _increase_ 'diatance_from_bellevue' by 1 unit we can expect the price of a house to _decrease_ by _18670_.- The _P-value_ tells us that we can reject the null hypothesis that the _Coefficient = 0_. In other words a changes in 'distance_from_bellevue' are related to changes in 'price'.- There is a Skewness of 0.65 which tells us that the data is moderately skewed to the right. The _QQ plot_ looks pretty good so the normality assumption seems fulfilled aswell.<jupyter_code># Build a formula.
f = 'price~grade'
# Create a fitted model
model = smf.ols(formula=f, data=kc_reg).fit()
# Get regression model summary
model.summary()
# Create a new Dataframe with the minimum and maximum value of 'grade'
X_new = pd.DataFrame({'grade': [kc_reg.grade.min(), kc_reg.grade.max()]});
# Make predictions on the x values and store them 
preds = model.predict(X_new)
# Plot the observed data and the least squares line
kc_reg.plot(kind='scatter', x='grade', y='price');
plt.plot(X_new, preds, c='red', linewidth=2);
plt.show()
# Visualize error terms for variance and heteroscedasticity
fig = plt.figure(figsize=(15,8))
fig = sm.graphics.plot_regress_exog(model, 'grade', fig=fig)
plt.show()
# Show distribution of residuals
sns.distplot(model.resid)
# Create qq plot to check normality assuption
residuals = model.resid
fig = sm.graphics.qqplot(residuals, dist=stats.norm, line='45', fit=True)<jupyter_output><empty_output><jupyter_text>Insights:- The _R-squared_ value is _0.307_ i.e, _31%_ of the variance in the target variable can be explained by 'grade'.- Just like 'sqft_living' and 'distance_from_bellevue', 'grade' is not a very good fit.- If we _increase_ 'grade' by 1 unit we can expect the price of a house to _increase_ by _127900_.- The _P-value_ tells us that we can reject the null hypothesis that the _Coefficient = 0_. In other words changes in 'grade' are related to changes in 'price'.- Based on the _QQ plot_ and residual distribution the normality assumption seems fufilled. > **We don't seem to have any reasonable single Linear Regressions, so we can move on to Multiple regression.**  Multiple RegressionIn order to construct a well functioning Multiple regression model we have to go through some prerequisites:- Deal with **Categorical data**. Which we've already did in the EDA notebook.- Handle **Multicollinearity of Features**. If we recall from the EDA notebook we found out the 'sqft_living' and 'sqft_above' are highly correlated. The best option is to get rid of one of the features but I've decided to keep them both as long as they are not used in the same model as this could affect our results.- **Feature Selection** is another step we can apply to our model when necessary. All features are important to a certain degree but not all of them are absolutely necessary so this step will help us weed out the features that may not be significant to our model.- **Feature Scaling and Normalization** is another important step in the process. We have a variety of Independent variables that happen to be on different scales, for example, we have 'bedrooms' ranging from 2 to 5 and on the other hand we have 'sqft_living' ranging roughly around 500 to 4000. 'sqft_living' is clearly on a larger scale and can be interpreted as more significant to the target variable than 'bedrooms'. Scaling and normalization can help avoid that.- **Train-Test Split** to make sure we are not _overfitting_ or _underfitting_ the model to the dataset.  Modeling for InferenceWhen modeling for _Inference_ we want to know how a change in _X(independent)_ affects and changes _Y(dependent)_. We are more concerned with the _Coefficients_ of the features than the predictive accuracy of the model.<jupyter_code># Let's first build a model with all the features.
# Define X and Y
X = kc_reg.drop(['price'], axis=1)
y = kc_reg['price']

# Show model summary table
X_int = sm.add_constant(X)
model = sm.OLS(y, X_int).fit()
model.summary()<jupyter_output>/usr/local/anaconda3/envs/learn-env/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)<jupyter_text>> The model has a $ R^2 $ of **0.821** or **82%**, which defines a pretty high _**Goodness of fit**_. The goal of this project is to find the features with the highest Coefficents. _Recursive Feature Elimination(RFE)_ lets us do just that.  Final ModelFor our final model we want the features with the highest Coefficents. There are currently 90 features in the dataset. Let's preselect some features to build our model on.<jupyter_code>predictors = kc_reg[['bedrooms', 'bathrooms', 'sqft_living', 'distance_from_bellevue', 'waterfront', 
                     'basement', 'renovated', 'yr_built', 'grade', 'condition', 'lat', 'long']]
predictors.head()
# Model with selected features
X_int = sm.add_constant(predictors)
model = sm.OLS(y, X_int).fit()
model.summary()<jupyter_output><empty_output><jupyter_text>> This model has a $ R^2 $ of **0.675** or **68%**, which defines a good enough _**Goodness of fit**_. With just 12 features we were able to get a $ R^2 $ just 12% short of the original model with all 89 features.  Feature Scaling - MinMaxScalerWhen working with data with different **datatypes** it is hard to guage the importance of each feature when they vary largely in magnitudes. For instance, take the 'sqft_living' and 'bedrooms' columns, the number of bedrooms ranges from about 1 to 4, sqft-living on the other hand ranges from about 1 to about 5000. It's easy to assume that sqft-living is more important as it is on a much grander scale, but this is not always the case. In other to get the most accurate coefficients for each feature we need to make sure that each feature is on the same grading scale and that can be done with the _**Sklearn's MinMaxScaler**_ function.<jupyter_code># Scale features with MinMaxScaler
scaler = MinMaxScaler()
scaler.fit(predictors)
scaled_data= scaler.transform(predictors)
predictors_scaled = pd.DataFrame(scaled_data)
# New dataframe with scaled predictors
predictors_scaled.columns = predictors.columns
predictors_scaled.head()<jupyter_output><empty_output><jupyter_text>_Recursive Feature Elimination(RFE)_Next we will be using RFE to determine which 6 of the 13 scaled features are more important.<jupyter_code># Use RFE to get features with the highest coefficients
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select=6)
selector = selector.fit(predictors_scaled, y)
# Boolean mask to pop out names of selected columns
selector.support_ * np.array(predictors_scaled.columns)
# Build a model with the selected features.
formula = 'price ~ sqft_living + distance_from_bellevue + waterfront + yr_built + grade + lat'
model = smf.ols(formula=formula, data=kc_reg).fit()
model.summary()
# Assign the final features to a new variable
final_feat = kc_reg[['sqft_living', 'distance_from_bellevue', 'waterfront', 'yr_built', 'grade', 'lat']]
# Print Coefficient of each feature in final model
linreg.fit(final_feat, y)
coefficients = dict(zip(final_feat,list(linreg.coef_)))
coefficients
# Check residual distribution
sns.distplot(model.resid)
plt.show()
# QQ Plot
fig = sm.graphics.qqplot(residuals, dist=stats.norm, line='45', fit=True)
# Visualize error terms for variance and heteroscedasticity
pred = model.predict(final_feat)
sns.scatterplot(pred, model.resid)
plt.axhline(0)
plt.show()<jupyter_output><empty_output><jupyter_text>> The final model has a $ R^2 $ of _**0.663**_, this tells us the the model explains _**66%**_ of the variability of the response data around it's mean. Looking at the displot graph ans QQ plot we can see that the data is moderately normal. The homoscedasticity assumption is fufilled.  Train-Test SplitNow that we've selected the features for our final model we will now test _overfitting_ and _underfitting_<jupyter_code># Split the data into Train-Test
X_train, X_test, y_train, y_test = train_test_split(final_feat, y, test_size=0.25, random_state=123)
# Check split sizes
print(len(X_train), len(X_test), len(y_train), len(y_test))
# Calculate and compare MSE of train and test data
model_test = linreg.fit(X_train, y_train)
y_hat_train = linreg.predict(X_train)
y_hat_test = linreg.predict(X_test)
# Calculate MSE with sklearn mean_squared_error() function
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
print('Train Mean Squarred Error:', train_mse)
print('Test Mean Squarred Error:', test_mse)

#Use .score to return R^2 of current test
R_squared = model_test.score(X_test, y_test)
print('This iteration has an R-squared value of:', R_squared)<jupyter_output>Train Mean Squarred Error: 12388530979.099497
Test Mean Squarred Error: 12371707262.100603
This iteration has an R-squared value of: 0.664077708172385
{'fraction_non_alphanumeric': 0.08095144227610224, 'fraction_numerical': 0.014061370655102176, 'mean_word_length': 5.336386344712739, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 42, 'https://': 1, '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': '340584', 'license': 'Apache-2.0', 'n_tokens_mistral': 4786, 'n_tokens_neox': 4452, 'n_words': 1901, 'path': 'Regression_Notebook.ipynb', 'repo_name': 'Princess-1/2_Module_Project', 'size': '15232'}