<jupyter_start><jupyter_text>Linear regression with statsmodels (OLS)<jupyter_code>import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tools.eval_measures as eval_measures
import seaborn as sns
import scipy.stats as stats
import pandas as pd
import numpy as np

from matplotlib import pyplot as plt
from sklearn import datasets, model_selection, metrics<jupyter_output><empty_output><jupyter_text>We load the boston house-prices dataset and `X` are our features and `y` is the target variable `medv` (Median value of owner-occupied homes in $1000s).<jupyter_code>boston = datasets.load_boston()
print(boston.DESCR)<jupyter_output>.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pu[...]<jupyter_text>Let's split the data in a test and training set.<jupyter_code>X = pd.DataFrame(boston.data, columns=boston.feature_names)
y = boston.target

X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, train_size=0.7)
df_train = pd.DataFrame(y_train, columns=['target'])
df_train['type'] = 'train'

df_test = pd.DataFrame(y_test, columns=['target'])
df_test['type'] = 'test'

df_set = df_train.append(df_test)

_ = sns.displot(df_set, x="target" ,hue="type", kind="kde", log_scale=False)<jupyter_output><empty_output><jupyter_text>Fitting models - the standard way Full model without an intercept<jupyter_code>model = sm.OLS(y_train, X_train)
result = model.fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")
result.pvalues < 0.05
predicted = result.predict(X_test)

print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))<jupyter_output>r2 score: 0.7423650859458641
mse: 22.426525918612626
rmse: 4.735665308973242
mae: 3.3603974168926563<jupyter_text>Full model with an intercept<jupyter_code>model = sm.OLS(y_train, sm.add_constant(X_train))
result = model.fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")

predicted = result.predict(sm.add_constant(X_test))

print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))<jupyter_output>r2 score: 0.7351303467796936
mse: 23.05629330095753
rmse: 4.801696918898311
mae: 3.5463377274856707<jupyter_text>Fitting models using R-style formulasWe can also fit a model with the R syntax `y ~ x_1 + x_2` and build some complexer models.<jupyter_code>dat = X_train.copy()
dat['MEDV'] = y_train
dat.head()<jupyter_output><empty_output><jupyter_text>Full model with an intercept<jupyter_code>result = smf.ols('MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B', data=dat).fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")
predicted = result.predict(X_test)

print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))<jupyter_output>r2 score: 0.7035468480293212
mse: 25.80556412834571
rmse: 5.079917728501684
mae: 3.7617323184841407<jupyter_text>Model with a polynomial and the target variable log transformed<jupyter_code>result = smf.ols('np.log(MEDV) ~ CRIM + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + B + pow(AGE, 2)', data=dat).fit()
print(result.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:           np.log(MEDV)   R-squared:                       0.678
Model:                            OLS   Adj. R-squared:                  0.669
Method:                 Least Squares   F-statistic:                     72.36
Date:                Wed, 10 Mar 2021   Prob (F-statistic):           2.60e-78
Time:                        21:16:10   Log-Likelihood:                 25.883
No. Observations:                 354   AIC:                            -29.77
Df Residuals:                     343   BIC:                             12.80
Df Model:                          10                                         
Covariance Type:            nonrobust                                         
===============================================================================
                  coef    std err          t      P[...]<jupyter_text>Let's plot the QQ-Plot for the residuals<jupyter_code>result.pvalues < 0.05
predicted = np.exp(result.predict(X_test))

print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))
_ = sm.qqplot(result.resid, fit=True, line="q")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1447686603256489, 'fraction_numerical': 0.04836369498629695, 'mean_word_length': 3.5955555555555554, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 32, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '211802', 'license': 'MIT', 'n_tokens_mistral': 2255, 'n_tokens_neox': 1989, 'n_words': 467, 'path': 'notebooks/supervised/regression/linear/multiple_linear_regression_statsmodels.ipynb', 'repo_name': 'rueedlinger/machine-learning-snippets-', 'size': '6203'}