2
$\begingroup$

I am learning about both the statsmodel library and sklearn. I am trying to construct a logistic model for both libraries trained on the same dataset.

In sklearn, the following works:

# import the data
from sklearn.datasets import load_breast_cancer

data = load_breast_cancer()

X_df = pd.DataFrame(data.data, columns=data.feature_names)
y_df = pd.DataFrame(data.target, columns=['target'])

# split into train and test datasets
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X_df, y_df, test_size=0.2, random_state=42)

print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)

# fit the model and make a prediction
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
Xscaled  = scaler.fit_transform(X_train)

clf1     = LogisticRegression()
clf1.fit(Xscaled, y_train.values.ravel())

y_pred = clf1.predict(scaler.fit_transform(X_test))

accuracy_score(y_test.values.ravel(), y_pred)*100)

This works and gives a result of

98.24561403508771

Now I want to do something similar in the statsmodel library

I do the following (continuing in the same notebook):

import statsmodels.api as sm

Xs = sm.add_constant(Xscaled)
res = sm.Logit(y_train, Xs).fit()

But this gives an error:

LinAlgError: Singular matrix

What is causing the discrpancy between sklearn and statsmodel?

$\endgroup$

3 Answers 3

8
$\begingroup$

I suspect the reason is that in scikit-learn the default logistic regression is not exactly logistic regression, but rather a penalized logistic regression (by default ridge-regresion i.e. with a L2-penalty). This has the result that it can provide estimates etc. even in case of perfect separation (e.g. some predictors have all 1 or all 0) or situations where some combination of predictors results in "perfect" prediction, while standard non-penalized logistic regression runs into problems (either you think this is a legitimate infite estimate, e.g. 100% of stones thrown into the air fall to the ground, or think of this as a problem e.g. 100% of 10 people that fell out of plane died, but occasionally people will survive).

$\endgroup$
2
$\begingroup$

In statsmodels, GLM may be more well developed than Logit. If you fit the model as below with GLM, it fails with a perfect separation error, which is exactly as it should. It is also possible to use fit_regularized to do L1 and/or L2 penalization to get parameter estimates in spite of the perfect separation.

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm
import pandas as pd

data = load_breast_cancer()

X_df = pd.DataFrame(data.data, columns=data.feature_names)
y_df = pd.DataFrame(data.target, columns=['target'])

X_train, X_test, y_train, y_test = train_test_split(X_df, y_df, test_size=0.2, random_state=42)

scaler = StandardScaler()
Xscaled = scaler.fit_transform(X_train)

Xs = sm.add_constant(Xscaled)

# Below will fail with a perfect separation error
#res = sm.GLM(y_train, Xs, family=sm.families.Binomial()).fit()

res2 = sm.GLM(y_train, Xs, family=sm.families.Binomial()).fit_regularized(L1_wt=0.0, alpha=0.1)
params = res2.params
$\endgroup$
0
$\begingroup$

In your data, the reason statsmodel crashes is due to highly correlated variables. If you look at your predictors:

import seaborn as sns
sns.clustermap(X_train.corr())

enter image description here

You can see this bunch of variables on diagonal left, perimeter,mean radius etc that is highly correlated. So if we include all of them into the fit, it is going to be hard to estimate these coefficients. We can remove these first, and fit:

import statsmodels.api as sm
excl = ['mean radius','mean perimeter','mean area',
        'worst area','worst perimeter','worst radius',
        'radius error','perimeter error','area error']
Xs = scaler.fit_transform(X_train.drop(excl,axis=1))
res = sm.Logit(y_train, Xs).fit()

Optimization terminated successfully.
         Current function value: 0.073780
         Iterations 12

This fits ok, we can check the accuracy and it's lower than what you have with scikit-learn, and I tried adding back the area, perimeter variables using PCA but it doesn't get much better:

y_scores = res.predict(scaler.fit_transform(X_test.drop(excl,axis=1)))
y_pred = (y_scores >0.5).astype(int)
accuracy_score(y_test.values.ravel(), y_pred)

0.9210526315789473

So fitting a model via maximum likelihood like in statsmodel is going to be highly unstable with high degree of collinearity. As @Björn pointed out, scikit-learn's logistic regression works because it by default uses a L2 penalty.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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