Solved – Wildly different $R^2$ between linear regression in statsmodels and sklearn

pythonregressionscikit learnstatsmodels

My question is related to:

Difference between statsmodel OLS and scikit linear regression

I essentially have the same problem, except my results are even more substantially different. Performing the following simple linear regressions, I get almost completely opposite results for the coefficient of determination:

import statsmodels.api as sm
from sklearn import linear_model

    x1 = [26.0, 31.0, 47.0, 51.0, 50.0, 49.0, 37.0, 33.0, 49.0, 54.0, 31.0, 49.0, 48.0, 49.0, 49.0, 47.0, 44.0, 48.0, 35.0, 43.0]
    y1 = [116.0, 94.0, 100.0, 102.0, 116.0, 116.0, 68.0, 118.0, 91.0, 104.0, 78.0, 116.0, 90.0, 109.0, 116.0, 118.0, 108.0, 119.0, 110.0, 102.0]

# Fit and summarize statsmodel OLS model
model_sm = sm.OLS(x1, y1)
result_sm = model_sm.fit()
print(result_sm.summary())


# Create sklearn linear regression object
ols_sk = linear_model.LinearRegression(fit_intercept=True)

# fit model
model_sk = ols_sk.fit(pd.DataFrame(x1), pd.DataFrame(y1))

# sklearn coefficient of determination
coefofdet = model_sk.score(pd.DataFrame(x1), pd.DataFrame(y1))

print('sklearn R^2: ' + str(coefofdet))

Statsmodels give me an $R^2$ of 0.962, while sklearn gives me an $R^2$ of 0.0584069073664.

What is causing such a drastic difference?

Best Answer

In your scikit-learn model, you included an intercept using the fit_intercept=True method. This fit both your intercept and the slope.

In statsmodels, if you want to include an intercept, you need to run the command x1 = stat.add_constant(x1) in order to create a column of constants. Then running the sm.OLS() command would yield an R-squared value of around 0.056.

It's also important to note that when constructing a model in statsmodels, you want to put your y1 first and x1 second rather than x1, then y1. The arguments are reversed from statsmodels to scikit-learn.

Related Solutions

Regression – Difference Between Statsmodel OLS and Scikit-Learn Linear Regression

First in terms of usage. You can get the prediction in statsmodels in a very similar way as in scikit-learn, except that we use the results instance returned by fit

predictions = results.predict(X_test)

Given the predictions, we can calculate statistics that are based on the prediction error

prediction_error = y_test - predictions

There is a separate list of functions to calculate goodness of prediction statistics with it, but it's not integrated into the models, nor does it include R squared. (I've never heard of R squared used for out of sample data.) Calculating those requires a bit more work by the user and statsmodels does not have the same set of statistics, especially not for classification or models with a binary response variable.

To your other two points:

Linear regression is in its basic form the same in statsmodels and in scikit-learn. However, the implementation differs which might produce different results in edge cases, and scikit learn has in general more support for larger models. For example, statsmodels currently uses sparse matrices in very few parts.

The most important difference is in the surrounding infrastructure and the use cases that are directly supported.

Statsmodels follows largely the traditional model where we want to know how well a given model fits the data, and what variables "explain" or affect the outcome, or what the size of the effect is. Scikit-learn follows the machine learning tradition where the main supported task is chosing the "best" model for prediction.

As a consequence, the emphasis in the supporting features of statsmodels is in analysing the training data which includes hypothesis tests and goodness-of-fit measures, while the emphasis in the supporting infrastructure in scikit-learn is on model selection for out-of-sample prediction and therefore cross-validation on "test data".

This points out the distinction, there is still quite a lot of overlap also in the usage. statsmodels also does prediction, and additionally forecasting in a time series context. But, when we want to do cross-validation for prediction in statsmodels it is currently still often easier to reuse the cross-validation setup of scikit-learn together with the estimation models of statsmodels.

Logistic Regression – Comparing Scikit Learn and Statsmodels

Your clue to figuring this out should be that the parameter estimates from the scikit-learn estimation are uniformly smaller in magnitude than the statsmodels counterpart. This might lead you to believe that scikit-learn applies some kind of parameter regularization. You can confirm this by reading the scikit-learn documentation.

There is no way to switch off regularization in scikit-learn, but you can make it ineffective by setting the tuning parameter C to a large number. Here is how that works in your case:

# module imports
from patsy import dmatrices
import pandas as pd
from sklearn.linear_model import LogisticRegression
import statsmodels.discrete.discrete_model as sm

# read in the data & create matrices
df = pd.read_csv("http://www.ats.ucla.edu/stat/data/binary.csv")
y, X = dmatrices('admit ~ gre + gpa + C(rank)', df, return_type = 'dataframe')

# sklearn output
model = LogisticRegression(fit_intercept = False, C = 1e9)
mdl = model.fit(X, y)
model.coef_

# sm
logit = sm.Logit(y, X)
logit.fit().params

UPDATE: As correctly pointed out in the comments below, now you can switch off the relularization in scikit-learn by setting penalty='none' (see the docs).

Best Answer

Related Solutions

Regression – Difference Between Statsmodel OLS and Scikit-Learn Linear Regression

Logistic Regression – Comparing Scikit Learn and Statsmodels

Related Question