3
$\begingroup$

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:

result from 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
"""
$\endgroup$
  • $\begingroup$ XtX_xT = XtX.dot(self.X.T). This should use np.linalg.solve instead of explicitly inverting the matrix. $\endgroup$ – Matthew Drury Apr 7 at 17:02
  • $\begingroup$ Are you using standardized data also with statsmodels. It does not add a constant automatically. Your + 1 * self.identity_matrix looks like Ridge regression to me and not OLS. $\endgroup$ – Josef Apr 9 at 0:34
3
$\begingroup$

Your problem is the algorithm by which you're calculating the coefficients and the variance-covariance matrix.

While it's (usually) algebraically fine, it's not numerically stable.

You will need to use a reasonably stable algorithm to do the calculations -- the way you're doing it really can be very inaccurate; it's quite possible for some values to be close, and for others to be nowhere near the correct least squares solution.

A common choice in a lot of regression routines is via QR decomposition of the X matrix; other choices make more sense in particular circumstances.

$\endgroup$
  • $\begingroup$ I have looked into using QR-decomposition for the calculation of co-efficients and agree that it seems your statement about my theory method being numerically unstable is true. Could you please point me in the right direction for a better way to calculate the variance-covariance matrix? $\endgroup$ – jupiar Apr 9 at 1:10
  • $\begingroup$ stackoverflow.com/questions/39568978/… $\endgroup$ – Glen_b Apr 9 at 1:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.