Why doScikit Learn's regression coefficients only match formal equation with standardized data?

Question

I was playing around with LinearRegression in Scikit Learn and I found a peculiarity that I'm trying to make sense out of.

If you compare values from the coef_ attribute, they only match the values given from the equation:

if you standardize the X feature matrix.

As an example, let's look at the results from the boston housing data set:

import pandas as pd
import numpy as np
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression

boston        = load_boston()
X             = pd.DataFrame(boston.data)  
X.columns     = boston.feature_names
y             = boston.target

With non-standardized data the results between the two methods differ:

# Comparison of Coefficients Without Standardized Data
np_coeffs = (np.linalg.inv(X.T.dot(X))).dot(X.T.dot(y))

lreg      = LinearRegression()
sk_coeffs = lreg.fit(X, y).coef_

# Put results into a dataframe for easy comparison
coeff_comparison1 = pd.DataFrame({
  'Variable' : X.columns,
  'np_coeffs': np_coeffs,
  'sk_coeffs': sk_coeffs
})

Which yields the following results:

If we standardize the data though, the results match exactly:

# comparison for coefficients with standardized data
X_std = (X - X.mean()) / X.std()

np_std_coeffs = (np.linalg.inv(X_std.T.dot(X_std))).dot(X_std.T.dot(y))
sk_std_coeffs = lreg.fit(X_std, y).coef_

coeff_comparison2 = pd.DataFrame({
  'Variable' : X.columns,
  'np_std_coeffs': np_std_coeffs,
  'sk_std_coeffs': sk_std_coeffs
})

Which gives the following dataframe:

And as you can see the results match exactly.

I'm curious what Scikit Learn does that causes one set of results to match those of the matrix equation, and the others to diverge.

Thank you.

Answer 1

When you wrote

np_coeffs = (np.linalg.inv(X.T.dot(X))).dot(X.T.dot(y))

you were solving for the coefficients for a linear model with no intercept term since you did not prepend a column of 1's to the feature matrix. When you standardize your data beforehand, the resulting linear model does not have an intercept term, so your results agree in that case.

Try instead

X1 = np.c_[np.ones(len(X)), X]  # Append a column for the intercept term
np_coeffs = (np.linalg.inv(X1.T.dot(X1))).dot(X1.T.dot(y))

and then compare:

lreg = LinearRegression().fit(X, y)  # Fit with original X matrix
sk_coeffs = [lreg.intercept_] + list(lreg.coef_)

coeff_comparison1 = pd.DataFrame({
    'np_coeffs': np_coeffs,
    'sk_coeffs': sk_coeffs},
    index=['Intercept'] + list(X.columns))

print(coeff_comparison1)

The printout is

            np_coeffs  sk_coeffs
Intercept  36.459488  36.459488
CRIM       -0.108011  -0.108011
ZN          0.046420   0.046420
INDUS       0.020559   0.020559
CHAS        2.686734   2.686734
NOX       -17.766611 -17.766611
RM          3.809865   3.809865
AGE         0.000692   0.000692
DIS        -1.475567  -1.475567
RAD         0.306049   0.306049
TAX        -0.012335  -0.012335
PTRATIO    -0.952747  -0.952747
B           0.009312   0.009312
LSTAT      -0.524758  -0.524758

as you would expect.