<jupyter_start><jupyter_text>Chapter 3 Applied Labs Q14 (dealing with collinearity)<jupyter_code>import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.stats.anova
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np

sns.set(style="whitegrid")

boston = sm.datasets.get_rdataset("Boston", "MASS").data
boston = sm.add_constant(boston)
rs = []
for predictor in boston.columns.drop(["crim", "const"]):
    print(predictor)
    model = smf.ols(f"crim ~ {predictor}", data=boston)
    results = model.fit()
    rs.append([predictor, results.rsquared])
    print(results.summary())
    
print(rs)<jupyter_output>[['zn', 0.04018790803211081], ['indus', 0.16531007043075163], ['chas', 0.0031238689633055206], ['nox', 0.17721718179269363], ['rm', 0.048069116716083604], ['age', 0.12442145175894637], ['dis', 0.14414937492539848], ['rad', 0.39125668674998926], ['tax', 0.3396142433788121], ['ptratio', 0.08406843894373672], ['black', 0.1482742394241312], ['lstat', 0.20759093253433558], ['medv', 0.15078046904975717]]<jupyter_text>All bar `chas` seem to be significant. The highest R-values come from `rad`, `tax`, and `lstat`<jupyter_code>sns.scatterplot("rad", "crim", data=boston) 
 sns.scatterplot("tax", "crim", data=boston) 
 sns.scatterplot("lstat", "crim", data=boston)<jupyter_output><empty_output><jupyter_text>(b) Fit a multiple regression model using all the predictors. For which predictors can you reject the hypothesis $\beta_j=0$ ?<jupyter_code>predictors = boston.columns.drop(["crim", "const"])
model = smf.ols(f"crim ~ {'+'.join(predictors)}" , data=boston)
results = model.fit()
print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                   crim   R-squared:                       0.454
Model:                            OLS   Adj. R-squared:                  0.440
Method:                 Least Squares   F-statistic:                     31.47
Date:                Fri, 17 May 2019   Prob (F-statistic):           1.57e-56
Time:                        17:29:45   Log-Likelihood:                -1653.3
No. Observations:                 506   AIC:                             3335.
Df Residuals:                     492   BIC:                             3394.
Df Model:                          13                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Only for `dis`, `rad`, and `medv` can we reject the null hypothesis at a p=0.01 level  (c) plot the multi-regression coefficients against the single-regression coefficients to show how they compare.<jupyter_code>single_betas = [] 
for predictor in boston.columns.drop(["crim", "const"]):
    model = smf.ols(f"crim ~ {predictor}", data=boston)
    results = model.fit()
    single_betas.append(results.params[1])
    
predictors = boston.columns.drop(["crim", "const"])
model = smf.ols(f"crim ~ {'+'.join(predictors)}" , data=boston)
results = model.fit()
multi_betas = results.params[1:]
fig, ax = plt.subplots(figsize=(15,10))
plt.plot(single_betas, multi_betas, 'o')
plt.xlabel("beta for single-factor regression")
plt.ylabel("beta for multi-factor regression")
plt.show()<jupyter_output><empty_output><jupyter_text>(d) Is there evidence of non-linear association between any of the predictors and the response? To answer this question, for each predictor X, fit a model of the form $Y = \beta_0 + \beta_1x + \beta_2 x^2 + \beta_3 x^3 + \epsilon $ for each predictor<jupyter_code>single_betas = [] 
for predictor in boston.columns.drop(["crim", "const"]):
    print(predictor)
    simple_model = smf.ols(f"crim ~ {predictor}", data=boston)
    simple_results = simple_model.fit()
    poly_model = smf.ols(f"crim ~ {predictor} + np.square({predictor}) + np.power({predictor}, 3)", data=boston)
    poly_results = poly_model.fit()
    # Get the ANOVA
    print(statsmodels.stats.anova.anova_lm(simple_results, poly_results))
    print("-"*80)<jupyter_output>zn
   df_resid           ssr  df_diff     ss_diff         F    Pr(>F)
0     504.0  35861.672261      0.0         NaN       NaN       NaN
1     502.0  35187.114178      2.0  674.558083  4.811821  0.008512
--------------------------------------------------------------------------------
indus
   df_resid           ssr  df_diff      ss_diff         F        Pr(>F)
0     504.0  31186.705132      0.0          NaN       NaN           NaN
1     502.0  27661.567803      2.0  3525.137328  31.98696  8.408754e-14
--------------------------------------------------------------------------------
chas
   df_resid           ssr  df_diff  ss_diff   F  Pr(>F)
0     504.0  37246.504181      0.0      NaN NaN     NaN
1     504.0  37246.504181     -0.0     -0.0 NaN     NaN
--------------------------------------------------------------------------------
nox
   df_resid           ssr  df_diff      ss_diff          F        Pr(>F)
0     504.0  30741.817087      0.0          NaN        NaN           NaN
1     50[...]
{'fraction_non_alphanumeric': 0.18999427152950163, 'fraction_numerical': 0.09299217109031889, 'mean_word_length': 2.176470588235294, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '280044', 'license': 'MIT', 'n_tokens_mistral': 1983, 'n_tokens_neox': 1659, 'n_words': 389, 'path': 'Chapter3/Labs/question15.ipynb', 'repo_name': 'wllgrnt/islr-python-examples', 'size': '5237'}