<jupyter_start><jupyter_text>Regression DiagnosticsThe purpose of this assigment is introduce you to a new library for linear regression called statmodels which is much better suited for inferential modeling than sklearn. This assignment is also to familiarize yourself with some of most important procedures for improving the interpretability of regression coefficients. You will also perform important statistical tests that will help establish that whether or not important assumptions that safeguard the interpretability of OLS coefficients have been met. We will continue to use the Ames Housing Dataset so that you can focus on the techniques and not on cleaning/getting associated with a brand new dataset.<jupyter_code>import seaborn as sns
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import stats<jupyter_output><empty_output><jupyter_text>1.1 Choose an X and Y variable from your dataset and use them to create a Seaborn Regplot<jupyter_code># Read in dataset 
df = pd.read_csv('https://raw.githubusercontent.com/ryanleeallred/datasets/master/Ames%20Housing%20Data/train.csv')

y = df['SalePrice']
X = df['YearBuilt']
sns.regplot(X, y)<jupyter_output><empty_output><jupyter_text>1.2 Now using the X variables that you feel like will be the best predictors of y use statsmodel to run the multiple regression between these variables and Y. You don't need to use every X variable in your dataset, in fact it's probably better if you don't. Just pick ones that you have already cleaned that seem the most relevant to house prices.<jupyter_code>X = df[['YearBuilt','OverallQual','Fireplaces','GarageArea','LotArea','GrLivArea','PoolArea','FullBath','HalfBath','YearRemodAdd']]
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())
df.info()<jupyter_output><empty_output><jupyter_text>1.3 Identify the standard errors and P-Values of these coefficients in the output table. What is the interpretation of the P-values here?<jupyter_code># All but 1 P-values indicate that we reject the null hypothesis that the coef is 0
# Pool Area in this model has a P-value of 0.368 which means that this feature does not help in accurately infering the Sale Price of a home<jupyter_output><empty_output><jupyter_text>1.4 Remove outliers from your dataset and run the regression again. Do you see a change in some coefficients? Which seem to move the most?<jupyter_code>X['Price'] = y
print(X.shape)
X = X[(np.abs(stats.zscore(X)) < 3).all(axis=1)]
print(X.shape)
# This got rid if all the PoolArea entries that weren't 0 since most of the houses don't have pools
y = X['Price']
X = X[['YearBuilt','OverallQual','Fireplaces','GarageArea','LotArea','GrLivArea','PoolArea','FullBath','HalfBath','YearRemodAdd']]
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())
# Fireplaces coef increaed by about 1000
# All of the coef moved<jupyter_output><empty_output><jupyter_text>1.5 Create a new log(y) variable and use it to run a log-linear regression of your variables using statmodels<jupyter_code>y = np.log(y)
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                  Price   R-squared:                       0.842
Model:                            OLS   Adj. R-squared:                  0.841
Method:                 Least Squares   F-statistic:                     815.0
Date:                Wed, 01 May 2019   Prob (F-statistic):               0.00
Time:                        21:52:55   Log-Likelihood:                 688.51
No. Observations:                1387   AIC:                            -1357.
Df Residuals:                    1377   BIC:                            -1305.
Df Model:                           9                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t     [...]<jupyter_text>2.1 Run a test for high levels of collinearity in your dataset. Calculate the Variance Inflation Factor for each X variable. Do you see VIF values greater than ten? If so try omitting those X variables and run your regression again. Do the standard errors change? Do the coefficients change? Do the coefficients seem to have an interpretation that matches your intuition?<jupyter_code>from statsmodels.stats.outliers_influence import variance_inflation_factor

X = sm.add_constant(X)
vif = [variance_inflation_factor(X.values, i) for i in range(len(X.columns))]
pd.Series(vif, X.columns)
X = X.drop('PoolArea', axis=1)
# PoolArea 
y_log = np.log(y)
model = sm.OLS(y_log, sm.add_constant(X))
results = model.fit()
print(results.summary())
# After Getting rid of PoolArea a few of the coef changed<jupyter_output><empty_output><jupyter_text>2.2 Variables that have high levels of multicollinearity should also be highly correlated with each other. Calculate your X matrix's correlation matrix to check if the variables highlighted by the VIF test truly are highly correlated.<jupyter_code>X.corr()
# over all qual had the highest VIF and it also has the most correlation with the other features 
# YearBuilt and overall qual make sense that they're are correlated because newer houses or most likey of higher quality 
# After adding more features to the model GrLivArea generated the largest VIF<jupyter_output><empty_output><jupyter_text>2.3 If you have variables with high Variance Inflation Factors, try excluding them from your regression. Do your standard errors improve? (get smaller). If high levels of multicollinearity are removed, the precision of the dataset should increase.<jupyter_code># After getting rid of PoolArea some of the atd error values went down to almost 0 but others shot up dramaticly! I guess they're compensating for the super low std errors?
X = X.drop('GrLivArea', axis=1)
y_log = np.log(y)
model = sm.OLS(y_log, sm.add_constant(X))
results = model.fit()
print(results.summary())
# Jb has decreased after omitting GrLivArea<jupyter_output><empty_output><jupyter_text>2.4 Recalculate your regression using Robust Standard Errors? What happens to your standard errors?<jupyter_code>model = sm.OLS(y_log, X)
results = model.fit(cov_type='HC3')
print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                  Price   R-squared:                       0.806
Model:                            OLS   Adj. R-squared:                  0.805
Method:                 Least Squares   F-statistic:                     511.7
Date:                Wed, 01 May 2019   Prob (F-statistic):               0.00
Time:                        21:58:57   Log-Likelihood:                 3989.2
No. Observations:                1387   AIC:                            -7960.
Df Residuals:                    1378   BIC:                            -7913.
Df Model:                           8                                         
Covariance Type:                  HC3                                         
================================================================================
                   coef    std err          z     [...]<jupyter_text>2.5 Use scatterplots or Seaborn's pairplot functionality to perform an eyeball test for potential variables that would be candidates for generating polynomial regressors.<jupyter_code>X.columns
features = X.columns
target = y_log
for feature in features:
  sns.scatterplot(X[feature], target)
  plt.show()<jupyter_output><empty_output><jupyter_text>2.6 Use seaborn's residplot to plot the distribution of each x variable's residuals. Does these plots indicate any other features that would be potential candidates for polynomial features.<jupyter_code>for feature in features:
    sns.residplot(X[feature], y, lowess=True, line_kws=dict(color='r'))
    plt.show()<jupyter_output>/usr/local/lib/python3.6/dist-packages/statsmodels/nonparametric/smoothers_lowess.py:165: RuntimeWarning: invalid value encountered in true_divide
  res = _lowess(y, x, frac=frac, it=it, delta=delta)
/usr/local/lib/python3.6/dist-packages/numpy/lib/function_base.py:3405: RuntimeWarning: Invalid value encountered in median
  r = func(a, **kwargs)
/usr/local/lib/python3.6/dist-packages/statsmodels/nonparametric/smoothers_lowess.py:165: RuntimeWarning: invalid value encountered in greater_equal
  res = _lowess(y, x, frac=frac, it=it, delta=delta)<jupyter_text>2.6 Feature Engineer the appropriate polynomial features from your analysis above and include them in one final log-polynomial, robust standard error, regression. Do the coefficients of this most advanced regression match your intuition better than the coefficients of the very first regression that we ran with the Statmodels library?<jupyter_code>X['YearBuilt**2'] = X['YearBuilt']**2
model = sm.OLS(y_log, X)
results = model.fit(cov_type='HC3')
print(results.summary())
X['LotArea**2'] = X['LotArea']**2
X['LotArea**3'] = X['LotArea']**3
model = sm.OLS(y_log, X)
results = model.fit(cov_type='HC3')
print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                  Price   R-squared:                       0.814
Model:                            OLS   Adj. R-squared:                  0.813
Method:                 Least Squares   F-statistic:                 1.259e+06
Date:                Wed, 01 May 2019   Prob (F-statistic):               0.00
Time:                        22:34:11   Log-Likelihood:                 4017.3
No. Observations:                1387   AIC:                            -8013.
Df Residuals:                    1376   BIC:                            -7955.
Df Model:                          10                                         
Covariance Type:                  HC3                                         
================================================================================
                   coef    std err          z     [...]
{'fraction_non_alphanumeric': 0.12107493739163938, 'fraction_numerical': 0.021286842612213447, 'mean_word_length': 2.8641607740975066, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '105502', 'license': 'MIT', 'n_tokens_mistral': 2982, 'n_tokens_neox': 2739, 'n_words': 984, 'path': 'module3-regression-diagnostics/regression-diagnostics-assignment.ipynb', 'repo_name': 'PWalis/DS-Unit-2-Sprint-2-Regression', 'size': '10382'}