I am using the answer given to this SO question to try to recover Standard Errors for a large linear regression from the LinearRegression() method in the python package sklearn.

I believe I am doing everything correctly, except I am using sparse matrices as I have a large predictor matrix with many zeros, for many fixed effects coded as dummy variables.

Using the sparse predictor matrix dataX (that already includes a column of 1's for the intercept) and the outcome variable y, I run:

from scipy import sparse
from sklearn.linear_model import LinearRegression
import pandas as pd
import numpy as np

mod = LinearRegression(fit_intercept=False).fit(dataX, y)
y_hat = mod.predict(dataX)
residuals = y - y_hat
N, p = dataX.shape
residual_sum_of_squares = residuals.T @ residuals
sigma_squared_hat = residual_sum_of_squares.iloc[0][0] / (N - p)
var_beta_hat = sparse.linalg.inv(sparse.csc_matrix(dataX.T @ dataX)) * sigma_squared_hat
diag = var_beta_hat.diagonal()
se = diag ** 0.5

coefs = pd.DataFrame({'names': labs, 'coefs': mod.coef_.tolist()[0], 'se': se})

But this yields many nan values. Upon inspecting the diag, which contains the variances, I am seeing negative values.

I have read that this may be due to serial correlation in my fixed effects. Does this mean that there is something wrong with my model? Can I still make valid inferences based on the coefficients, and the non-negative variances?


Your Answer

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

Browse other questions tagged or ask your own question.