In [24]:
import matplotlib.pyplot as plt
from pylab import mpl, plt
import numpy as np
from sklearn import datasets
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import datetime as dt
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy import stats
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'DejaVu Sans'
mpl.rcParams["savefig.dpi"] = 500
np.set_printoptions(precision=5, suppress=True, formatter={"float": lambda x: f"{x:6.3f}"})
%matplotlib inline
In [ ]:
In [6]:
# Load the diabetes dataset
diabetes_X, diabetes_y = datasets.load_diabetes(return_X_y=True)
diabetes_X
Out[6]:
array([[ 0.03807591, 0.05068012, 0.06169621, ..., -0.00259226,
0.01990842, -0.01764613],
[-0.00188202, -0.04464164, -0.05147406, ..., -0.03949338,
-0.06832974, -0.09220405],
[ 0.08529891, 0.05068012, 0.04445121, ..., -0.00259226,
0.00286377, -0.02593034],
...,
[ 0.04170844, 0.05068012, -0.01590626, ..., -0.01107952,
-0.04687948, 0.01549073],
[-0.04547248, -0.04464164, 0.03906215, ..., 0.02655962,
0.04452837, -0.02593034],
[-0.04547248, -0.04464164, -0.0730303 , ..., -0.03949338,
-0.00421986, 0.00306441]])
In [7]:
# Use only one feature
diabetes_X = diabetes_X[:, np.newaxis, 2]
# Split the data into training/testing sets
diabetes_X_train = diabetes_X[:-20]
diabetes_X_test = diabetes_X[-20:]
# Split the targets into training/testing sets
diabetes_y_train = diabetes_y[:-20]
diabetes_y_test = diabetes_y[-20:]
In [25]:
# Create linear regression object
regr = LinearRegression()
# Train the model using the training sets
regr.fit(diabetes_X_train, diabetes_y_train)
# Make predictions using the testing set
diabetes_y_pred = regr.predict(diabetes_X_test)
In [14]:
type(regr.coef_)
Out[14]:
numpy.ndarray
In [19]:
np.ndarray.round(regr.coef_,4)
Out[19]:
array([938.2379])
In [21]:
# The coefficients
print("Coefficients: \n", np.ndarray.round(regr.coef_,4))
# The mean squared error
print("Mean squared error: %.2f" % mean_squared_error(diabetes_y_test, diabetes_y_pred))
# The coefficient of determination: 1 is perfect prediction
print("Coefficient of determination: %.2f" % r2_score(diabetes_y_test, diabetes_y_pred))
Coefficients: [938.2379] Mean squared error: 2548.07 Coefficient of determination: 0.47
In [4]:
# Plot outputs
plt.scatter(diabetes_X_test, diabetes_y_test, color="black")
plt.plot(diabetes_X_test, diabetes_y_pred, color="blue", linewidth=3)
plt.xticks(())
plt.yticks(())
plt.show()
In [36]:
import statsmodels.api as sm
X_ols = sm.add_constant(diabetes_X_train)
model = sm.OLS(diabetes_y_train, X_ols).fit()
print(model.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.335
Model: OLS Adj. R-squared: 0.334
Method: Least Squares F-statistic: 211.7
Date: Sun, 29 May 2022 Prob (F-statistic): 3.98e-39
Time: 01:49:43 Log-Likelihood: -2346.5
No. Observations: 422 AIC: 4697.
Df Residuals: 420 BIC: 4705.
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 152.9189 3.069 49.830 0.000 146.887 158.951
x1 938.2379 64.484 14.550 0.000 811.487 1064.989
==============================================================================
Omnibus: 12.587 Durbin-Watson: 1.839
Prob(Omnibus): 0.002 Jarque-Bera (JB): 7.334
Skew: 0.142 Prob(JB): 0.0255
Kurtosis: 2.420 Cond. No. 21.0
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [38]:
X_ols
Out[38]:
array([[ 1.000, 0.062],
[ 1.000, -0.051],
[ 1.000, 0.044],
[ 1.000, -0.012],
[ 1.000, -0.036],
[ 1.000, -0.041],
[ 1.000, -0.047],
[ 1.000, -0.002],
[ 1.000, 0.062],
[ 1.000, 0.039],
[ 1.000, -0.084],
[ 1.000, 0.018],
[ 1.000, -0.029],
[ 1.000, -0.002],
[ 1.000, -0.026],
[ 1.000, -0.018],
[ 1.000, 0.042],
[ 1.000, 0.012],
[ 1.000, -0.011],
[ 1.000, -0.018],
[ 1.000, -0.057],
[ 1.000, -0.022],
[ 1.000, -0.004],
[ 1.000, 0.061],
[ 1.000, 0.036],
[ 1.000, -0.013],
[ 1.000, -0.077],
[ 1.000, 0.060],
[ 1.000, -0.021],
[ 1.000, -0.006],
[ 1.000, 0.044],
[ 1.000, -0.065],
[ 1.000, 0.125],
[ 1.000, -0.050],
[ 1.000, -0.063],
[ 1.000, -0.031],
[ 1.000, 0.023],
[ 1.000, 0.011],
[ 1.000, 0.071],
[ 1.000, 0.014],
[ 1.000, -0.008],
[ 1.000, -0.068],
[ 1.000, -0.011],
[ 1.000, -0.023],
[ 1.000, 0.068],
[ 1.000, -0.035],
[ 1.000, -0.012],
[ 1.000, -0.073],
[ 1.000, -0.042],
[ 1.000, 0.014],
[ 1.000, -0.007],
[ 1.000, 0.016],
[ 1.000, -0.009],
[ 1.000, -0.016],
[ 1.000, 0.025],
[ 1.000, -0.049],
[ 1.000, 0.041],
[ 1.000, -0.063],
[ 1.000, -0.064],
[ 1.000, -0.026],
[ 1.000, -0.004],
[ 1.000, 0.005],
[ 1.000, -0.007],
[ 1.000, -0.037],
[ 1.000, -0.026],
[ 1.000, -0.025],
[ 1.000, -0.018],
[ 1.000, -0.015],
[ 1.000, -0.030],
[ 1.000, -0.046],
[ 1.000, -0.070],
[ 1.000, 0.034],
[ 1.000, -0.004],
[ 1.000, -0.020],
[ 1.000, 0.002],
[ 1.000, -0.031],
[ 1.000, 0.028],
[ 1.000, -0.036],
[ 1.000, -0.058],
[ 1.000, -0.037],
[ 1.000, 0.012],
[ 1.000, -0.022],
[ 1.000, -0.035],
[ 1.000, 0.010],
[ 1.000, -0.040],
[ 1.000, 0.071],
[ 1.000, -0.075],
[ 1.000, -0.006],
[ 1.000, -0.041],
[ 1.000, -0.048],
[ 1.000, -0.026],
[ 1.000, 0.052],
[ 1.000, 0.005],
[ 1.000, -0.064],
[ 1.000, -0.017],
[ 1.000, -0.058],
[ 1.000, 0.010],
[ 1.000, 0.089],
[ 1.000, -0.005],
[ 1.000, -0.064],
[ 1.000, 0.018],
[ 1.000, -0.045],
[ 1.000, 0.028],
[ 1.000, 0.041],
[ 1.000, 0.065],
[ 1.000, -0.032],
[ 1.000, -0.076],
[ 1.000, 0.050],
[ 1.000, 0.046],
[ 1.000, -0.009],
[ 1.000, -0.032],
[ 1.000, 0.005],
[ 1.000, 0.021],
[ 1.000, 0.014],
[ 1.000, 0.110],
[ 1.000, 0.001],
[ 1.000, 0.058],
[ 1.000, -0.021],
[ 1.000, -0.011],
[ 1.000, -0.047],
[ 1.000, 0.005],
[ 1.000, 0.018],
[ 1.000, 0.081],
[ 1.000, 0.035],
[ 1.000, 0.024],
[ 1.000, -0.008],
[ 1.000, -0.061],
[ 1.000, -0.002],
[ 1.000, -0.062],
[ 1.000, 0.016],
[ 1.000, 0.096],
[ 1.000, -0.070],
[ 1.000, -0.021],
[ 1.000, -0.054],
[ 1.000, 0.043],
[ 1.000, 0.056],
[ 1.000, -0.082],
[ 1.000, 0.050],
[ 1.000, 0.111],
[ 1.000, 0.062],
[ 1.000, 0.014],
[ 1.000, 0.048],
[ 1.000, 0.012],
[ 1.000, 0.006],
[ 1.000, 0.047],
[ 1.000, 0.129],
[ 1.000, 0.060],
[ 1.000, 0.093],
[ 1.000, 0.015],
[ 1.000, -0.005],
[ 1.000, 0.070],
[ 1.000, -0.004],
[ 1.000, -0.001],
[ 1.000, -0.044],
[ 1.000, 0.021],
[ 1.000, 0.061],
[ 1.000, -0.011],
[ 1.000, -0.033],
[ 1.000, -0.065],
[ 1.000, 0.043],
[ 1.000, -0.062],
[ 1.000, 0.064],
[ 1.000, 0.030],
[ 1.000, 0.072],
[ 1.000, -0.019],
[ 1.000, -0.067],
[ 1.000, -0.060],
[ 1.000, 0.069],
[ 1.000, 0.060],
[ 1.000, -0.027],
[ 1.000, -0.020],
[ 1.000, -0.046],
[ 1.000, 0.071],
[ 1.000, -0.079],
[ 1.000, 0.010],
[ 1.000, -0.039],
[ 1.000, 0.020],
[ 1.000, 0.027],
[ 1.000, -0.008],
[ 1.000, -0.016],
[ 1.000, 0.005],
[ 1.000, -0.043],
[ 1.000, 0.006],
[ 1.000, -0.035],
[ 1.000, 0.024],
[ 1.000, -0.018],
[ 1.000, 0.042],
[ 1.000, -0.055],
[ 1.000, -0.003],
[ 1.000, -0.067],
[ 1.000, -0.013],
[ 1.000, -0.042],
[ 1.000, -0.031],
[ 1.000, -0.005],
[ 1.000, -0.059],
[ 1.000, 0.025],
[ 1.000, -0.046],
[ 1.000, 0.003],
[ 1.000, 0.054],
[ 1.000, -0.045],
[ 1.000, -0.058],
[ 1.000, -0.056],
[ 1.000, 0.001],
[ 1.000, 0.030],
[ 1.000, 0.007],
[ 1.000, 0.047],
[ 1.000, 0.026],
[ 1.000, 0.046],
[ 1.000, 0.040],
[ 1.000, -0.018],
[ 1.000, 0.014],
[ 1.000, 0.037],
[ 1.000, 0.003],
[ 1.000, -0.071],
[ 1.000, -0.033],
[ 1.000, 0.094],
[ 1.000, 0.036],
[ 1.000, 0.032],
[ 1.000, -0.065],
[ 1.000, -0.042],
[ 1.000, -0.040],
[ 1.000, -0.039],
[ 1.000, -0.026],
[ 1.000, -0.023],
[ 1.000, -0.067],
[ 1.000, 0.033],
[ 1.000, -0.046],
[ 1.000, -0.030],
[ 1.000, -0.013],
[ 1.000, -0.016],
[ 1.000, 0.071],
[ 1.000, -0.031],
[ 1.000, 0.000],
[ 1.000, 0.037],
[ 1.000, 0.039],
[ 1.000, -0.015],
[ 1.000, 0.007],
[ 1.000, -0.069],
[ 1.000, -0.009],
[ 1.000, 0.020],
[ 1.000, 0.075],
[ 1.000, -0.008],
[ 1.000, -0.023],
[ 1.000, -0.046],
[ 1.000, 0.054],
[ 1.000, -0.035],
[ 1.000, -0.032],
[ 1.000, -0.082],
[ 1.000, 0.048],
[ 1.000, 0.061],
[ 1.000, 0.056],
[ 1.000, 0.098],
[ 1.000, 0.060],
[ 1.000, 0.034],
[ 1.000, 0.056],
[ 1.000, -0.065],
[ 1.000, 0.161],
[ 1.000, -0.056],
[ 1.000, -0.025],
[ 1.000, -0.036],
[ 1.000, -0.008],
[ 1.000, -0.042],
[ 1.000, 0.127],
[ 1.000, -0.077],
[ 1.000, 0.028],
[ 1.000, -0.026],
[ 1.000, -0.062],
[ 1.000, -0.001],
[ 1.000, 0.089],
[ 1.000, -0.032],
[ 1.000, 0.030],
[ 1.000, 0.009],
[ 1.000, 0.007],
[ 1.000, -0.020],
[ 1.000, -0.025],
[ 1.000, -0.012],
[ 1.000, 0.026],
[ 1.000, -0.059],
[ 1.000, -0.036],
[ 1.000, -0.025],
[ 1.000, 0.019],
[ 1.000, -0.090],
[ 1.000, -0.005],
[ 1.000, -0.053],
[ 1.000, -0.022],
[ 1.000, -0.020],
[ 1.000, -0.055],
[ 1.000, -0.006],
[ 1.000, -0.017],
[ 1.000, 0.055],
[ 1.000, 0.077],
[ 1.000, 0.019],
[ 1.000, -0.022],
[ 1.000, 0.093],
[ 1.000, -0.031],
[ 1.000, 0.039],
[ 1.000, -0.061],
[ 1.000, -0.008],
[ 1.000, -0.037],
[ 1.000, -0.014],
[ 1.000, 0.074],
[ 1.000, -0.025],
[ 1.000, 0.034],
[ 1.000, 0.035],
[ 1.000, -0.039],
[ 1.000, -0.040],
[ 1.000, -0.002],
[ 1.000, -0.031],
[ 1.000, -0.046],
[ 1.000, 0.001],
[ 1.000, 0.065],
[ 1.000, 0.040],
[ 1.000, -0.023],
[ 1.000, 0.053],
[ 1.000, 0.040],
[ 1.000, -0.020],
[ 1.000, 0.014],
[ 1.000, -0.034],
[ 1.000, 0.007],
[ 1.000, 0.005],
[ 1.000, 0.030],
[ 1.000, 0.052],
[ 1.000, 0.062],
[ 1.000, -0.007],
[ 1.000, 0.006],
[ 1.000, 0.054],
[ 1.000, -0.008],
[ 1.000, 0.115],
[ 1.000, 0.067],
[ 1.000, -0.056],
[ 1.000, 0.030],
[ 1.000, -0.026],
[ 1.000, 0.105],
[ 1.000, -0.006],
[ 1.000, -0.047],
[ 1.000, -0.048],
[ 1.000, 0.085],
[ 1.000, -0.013],
[ 1.000, -0.033],
[ 1.000, -0.007],
[ 1.000, -0.014],
[ 1.000, 0.060],
[ 1.000, 0.022],
[ 1.000, 0.019],
[ 1.000, -0.012],
[ 1.000, -0.003],
[ 1.000, 0.018],
[ 1.000, -0.030],
[ 1.000, -0.020],
[ 1.000, -0.058],
[ 1.000, 0.061],
[ 1.000, -0.041],
[ 1.000, -0.072],
[ 1.000, -0.056],
[ 1.000, 0.046],
[ 1.000, -0.009],
[ 1.000, -0.033],
[ 1.000, 0.050],
[ 1.000, -0.085],
[ 1.000, 0.006],
[ 1.000, 0.021],
[ 1.000, -0.007],
[ 1.000, 0.105],
[ 1.000, -0.025],
[ 1.000, -0.006],
[ 1.000, -0.039],
[ 1.000, 0.137],
[ 1.000, 0.171],
[ 1.000, 0.002],
[ 1.000, 0.038],
[ 1.000, -0.058],
[ 1.000, -0.009],
[ 1.000, -0.023],
[ 1.000, -0.011],
[ 1.000, -0.034],
[ 1.000, -0.003],
[ 1.000, 0.068],
[ 1.000, 0.010],
[ 1.000, 0.002],
[ 1.000, -0.039],
[ 1.000, 0.026],
[ 1.000, -0.089],
[ 1.000, 0.061],
[ 1.000, -0.029],
[ 1.000, -0.030],
[ 1.000, -0.019],
[ 1.000, -0.041],
[ 1.000, 0.015],
[ 1.000, -0.025],
[ 1.000, 0.001],
[ 1.000, 0.069],
[ 1.000, -0.070],
[ 1.000, -0.030],
[ 1.000, -0.046],
[ 1.000, 0.019],
[ 1.000, 0.001],
[ 1.000, -0.031],
[ 1.000, -0.004],
[ 1.000, 0.015],
[ 1.000, 0.023],
[ 1.000, 0.046],
[ 1.000, -0.045],
[ 1.000, -0.033],
[ 1.000, 0.097],
[ 1.000, 0.054],
[ 1.000, 0.123],
[ 1.000, -0.081],
[ 1.000, 0.093],
[ 1.000, -0.050],
[ 1.000, -0.012],
[ 1.000, -0.028],
[ 1.000, 0.058],
[ 1.000, 0.085],
[ 1.000, -0.001],
[ 1.000, 0.007],
[ 1.000, 0.009],
[ 1.000, 0.080],
[ 1.000, 0.071],
[ 1.000, -0.025],
[ 1.000, -0.055],
[ 1.000, -0.036],
[ 1.000, 0.016]])
In [35]:
mod = sm.OLS(diabetes_y_train, sm.add_constant(diabetes_X_train))
res = mod.fit()
print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.335
Model: OLS Adj. R-squared: 0.334
Method: Least Squares F-statistic: 211.7
Date: Sun, 29 May 2022 Prob (F-statistic): 3.98e-39
Time: 01:49:37 Log-Likelihood: -2346.5
No. Observations: 422 AIC: 4697.
Df Residuals: 420 BIC: 4705.
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 152.9189 3.069 49.830 0.000 146.887 158.951
x1 938.2379 64.484 14.550 0.000 811.487 1064.989
==============================================================================
Omnibus: 12.587 Durbin-Watson: 1.839
Prob(Omnibus): 0.002 Jarque-Bera (JB): 7.334
Skew: 0.142 Prob(JB): 0.0255
Kurtosis: 2.420 Cond. No. 21.0
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [39]:
import pandas as pd
beta=np.linalg.inv(X_ols.T.dot(X_ols)).dot(X_ols.T.dot(diabetes_y_train))
pd.Series(beta)
Out[39]:
0 152.918862 1 938.237861 dtype: float64