Recently I've been going back to the core of multiple linear regression through theory books, and rather than use standard packages, for example python package statsmodels, to calculate summary statistics, I've been trying to re-create those results using the underlying theory. The only problem is that while my results are similar, there is one still very different.
So, first, the following result is from using the statsmodels package:
My currrent method for calculating the coefficients is:
$$\hat{\beta}=(\tilde{X}^{T}\tilde{X})^{-1}\tilde{X}^{T}\tilde{Y}$$
With $X$ and $Y$ both centered and scaled, using: $$\tilde{Y} = \frac{Y-\bar{y}}{s_{y}}$$ $$\tilde{X}_{j} = \frac{X_{j}-\bar{x}_{j}}{s_{j}}$$ $$s_{j} = \sqrt{\frac{\sum_{i=1}^{n}(x_{ij} - \bar{x}_{j})^2}{n-1}}$$ $$s_{y} = \sqrt{\frac{\sum_{i=1}^{n}(y_{i} - \bar{y}_{j})^2}{n-1}}$$
My variance estimate being: $$\hat{\sigma}^{2} = \frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{n-p-1}$$ The standard error of $\beta_{j}$ being: $$C = (\tilde{X}^{T}\tilde{X})^{-1}$$ $$s.e.(\beta_{j}) = \hat{\sigma}\sqrt{c_{jj}}$$ The above $c_{jj}$ is the $j^{th}$ diagonal element of $C$.
And then the t-statistic calculated with $\beta_{j}^{0}=0$: $$t_{j}=\frac{\hat{\beta}_{j} - \beta_{j}^{0}}{s.e.(\hat{\beta}_{j})}$$
With the above the following table comparing statsmodels as sm and my using the theory above, the row with the * highlights the one with the massive difference.
| sm.beta | my.beta | sm.t | my.t
--------------------------------------
| 0.6234 | 0.6295 | 3.892 | 3.816
| -0.0585 | -0.0536 | -0.435 | -0.406
| 0.3450 | 0.2980 | 2.079 | 1.921
| 0.0991 | 0.0923 | 0.450 | 0.522
*| 0.1293 | 0.0271 | 1.180 | 0.232
| -0.2197 | -0.1808 | -1.236 | -1.262
I just can't understand why all the other results are so close to the standard package. The small differences are most likely due to a much more efficient linear algebra solution with less rounding errors in their source code statsmodels OLS source code.
The dataset is shown below:
columns = [
"Overall rating of job being done by supervisor",
"Handles employee complaints",
"Does not allow special privileges",
"Opportunity to learn new things",
"Raises based on performance",
"Too critical of poor performance",
"Rate of advancing to better jobs"
]
rows = [
[43, 51, 30, 39, 61, 92, 45],
[63, 64, 51, 54, 63, 73, 47],
[71, 70, 68, 69, 76, 86, 48],
[61, 63, 45, 47, 54, 84, 35],
[81, 78, 56, 66, 71, 83, 47],
[43, 55, 49, 44, 54, 49, 34],
[58, 67, 42, 56, 66, 68, 35],
[71, 75, 50, 55, 70, 66, 41],
[72, 82, 72, 67, 71, 83, 31],
[67, 61, 45, 47, 62, 80, 41],
[64, 53, 53, 58, 58, 67, 34],
[67, 60, 47, 39, 59, 74, 41],
[69, 62, 57, 42, 55, 63, 25],
[68, 83, 83, 45, 59, 77, 35],
[77, 77, 54, 72, 79, 77, 46],
[81, 90, 50, 72, 60, 54, 36],
[74, 85, 64, 69, 79, 79, 63],
[65, 60, 65, 75, 55, 80, 60],
[65, 70, 46, 57, 75, 85, 46],
[50, 58, 68, 54, 64, 78, 52],
[50, 40, 33, 34, 43, 64, 33],
[64, 61, 52, 62, 66, 80, 41],
[53, 66, 52, 50, 63, 80, 37],
[40, 37, 42, 58, 50, 57, 49],
[63, 54, 42, 48, 66, 75, 33],
[66, 77, 66, 63, 88, 76, 72],
[78, 75, 58, 74, 80, 78, 49],
[48, 57, 44, 45, 51, 83, 38],
[85, 85, 71, 71, 77, 74, 55],
[82, 82, 39, 59, 64, 78, 39]
]
And the calculation code is shown below, a lot of it can be better written, I just quickly wrote it to make sure that the calculations are coming out ok.
import typing as ty
import numpy as np
class MultipleLinearRegression():
def __init__(self, X, y):
self.X = X
self.y = y
self.identity_size = X.shape[1]
self.identity_matrix = np.zeros((self.identity_size, self.identity_size))
np.fill_diagonal(self.identity_matrix, 1)
self.number_of_observations, self.number_of_predictors = self.__get_X_meta_info()
self.n_m_p_m_1 = self.number_of_observations - self.number_of_predictors - 1
def summary(self) -> None:
""" Print summary to screen """
# Center and scale the variables
self.X = self.__center_and_scale_x()
self.y = self.__center_and_scale_y()
beta_coefficients = self.__get_beta_coefficients()
variance = self.__get_variance(beta_coefficients)
big_c = self.__get_big_c()
print("idx \t beta \t t")
for idx, beta_coefficient in enumerate(beta_coefficients):
se = self.__get_standard_error(variance, big_c, idx)
res = (beta_coefficient / se)
print(f"{idx}\t{beta_coefficient}\t{res}")
return None
def __center_and_scale_x(self) -> ty.List:
""" Centering and scaling should be done """
standardized_x = []
for idx in range(self.number_of_predictors):
values = self.X[:, idx]
avg = np.average(values)
top_sum = sum([(i-avg)**2 for i in values])
sd = np.sqrt(top_sum / (len(values)-1))
col = []
for value in values:
col.append((value-avg) / sd)
standardized_x.append(col)
standardized_x = np.array(standardized_x)
transpose = standardized_x.T
return transpose
def __center_and_scale_y(self) -> ty.List:
""" Centering and scaling should be done """
values = self.y
avg = np.average(values)
top_sum = sum([(i-avg)**2 for i in values])
sd = np.sqrt(top_sum / (len(values)-1))
col = []
for value in values:
col.append((value-avg) / sd)
standardized_col = np.array(col)
return standardized_col
def __get_big_c(self):
""" Calculate $(X^{T}X)^{-1}$ """
xTx = self.X.T.dot(self.X) + 1 * self.identity_matrix
XtX = np.linalg.inv(xTx)
return XtX
def __get_variance(self, beta_coefficients: ty.List) -> float:
""" Get variance for y estimation """
y_list = []
for idx, row in enumerate(self.X):
y_hat = beta_coefficients[0]*row[0] + \
beta_coefficients[1]*row[1] + \
beta_coefficients[2]*row[2] + \
beta_coefficients[3]*row[3] + \
beta_coefficients[4]*row[4] + \
beta_coefficients[5]*row[5]
y = self.y[idx]
s1 = (y - y_hat)**2
y_list.append(s1)
return sum(y_list) / self.n_m_p_m_1
def __get_X_meta_info(self) -> tuple:
""" Get number of observations and predictors """
number_of_observations = X.shape[0]
number_of_predictors = X.shape[1]
return (number_of_observations, number_of_predictors)
def __get_standard_error(self, var, big_c, j):
""" Get standard error for a coefficient """
return np.sqrt(var) * np.sqrt(big_c[j][j])
def __get_standard_deviation(self, column_idx: int) -> float:
""" Get standard deviation for a coefficient """
values = self.X[:, column_idx]
avg = np.average(values)
top_sum = sum([(i-avg)**2 for i in values])
sd = np.sqrt(top_sum / (len(values)-1))
nsd = np.std(values)
# It should be noted, that numpy standard deviation DOES NOT use n-1
# when dividing top sum.
return sd
def __get_beta_coefficients(self) -> ty.List:
""" Get the beta coefficients for the model """
xTx = self.X.T.dot(self.X) + 1 * self.identity_matrix
XtX = np.linalg.inv(xTx)
XtX_xT = XtX.dot(self.X.T)
theta = XtX_xT.dot(self.y)
return theta
x = []
y = []
for row in rows:
y.append(row[0])
x.append(row[1:])
X = np.array(x, np.int32)
y = np.array(y, np.int32)
mlr = MultipleLinearRegression(X, y)
mlr.summary()
"""
idx beta t
0 0.6295040961600624 3.8157747856899435
1 -0.05355985635512084 -0.40554514013311027
2 0.29803825790238525 1.9212304315894024
3 0.09230975543564345 0.5219481794579525
4 0.027097331841556226 0.23236578727430382
5 -0.18079045869224042 -1.2620479201518775
"""
XtX_xT = XtX.dot(self.X.T). This should usenp.linalg.solveinstead of explicitly inverting the matrix. $\endgroup$+ 1 * self.identity_matrixlooks like Ridge regression to me and not OLS. $\endgroup$