Logistic Regression Scikit-learn vs Statsmodels

What’s the difference between Statsmodels and Scikit-learn? Both have ordinary least squares and logistic regression, so it seems like Python is giving us two ways to do the same thing. Statsmodels offers modeling from the perspective of statistics. Scikit-learn offers some of the same models from the perspective of machine learning.

So we need to understand the difference between statistics and machine learning! Statistics makes mathematically valid inferences about a population based on sample data. Statistics answers the question, “What is the evidence that X is related to Y?” Machine learning has the goal of optimizing predictive accuracy rather than inference. Machine learning answers the question, “Given X, what prediction should we make for Y?”

In the example below, we’ll create a fake dataset with predictor variables and a binary Y variable. Then we’ll perform logistic regression with scikit-learn and statsmodels. We’ll see that scikit-learn allows us to easily tune the model to optimize predictive power. Statsmodels will provide a summary of statistical measures which will be very familiar to those who’ve used SAS or R.

If you need an intro to Logistic Regression, see this Finxter post.

Create Fake Data for the Logistic Regression Model

I tried using some publicly available data for this exercise but didn’t find one with the characteristics I wanted. So I decided to create some fake data by using NumPy! There’s a post here that explains the math and how to do this in R.

import numpy as np
import pandas as pd

#The next line is setting the seed for the random number generator so that we get consistent results
rg = np.random.default_rng(seed=0)
#Create an array with 500 rows and 3 columns
X_for_creating_probabilities = rg.normal(size=(500,3))

Create an array with the first column removed. The deleted column can be thought of as random noise, or as a variable that we don’t have access to when creating the model.

X1 = np.delete(X_for_creating_probabilities,0,axis=1)
X1[:5]
"""
array([[-0.13210486,  0.64042265],
       [-0.53566937,  0.36159505],
       [ 0.94708096, -0.70373524],
       [-0.62327446,  0.04132598],
       [-0.21879166, -1.24591095]])
"""

Now we’ll create two more columns correlated with X1. Datasets often have highly correlated variables. Correlation increases the likelihood of overfitting. Concatenate to get a single array.

X2 = X1 + .1 * np.random.normal(size=(500,2))
X_predictors = np.concatenate((X1,X2),axis=1)

We want to create our outcome variable and have it be related to X_predictors. To do that, we use our data as inputs to the logistic regression model to get probabilities. Then we set the outcome variable, Y, to True when the probability is above .5.

P = 1 / (1 + np.e**(-np.matmul(X_for_creating_probabilities,[1,1,1])))
Y = P > .5
#About half of cases are True
np.mean(Y)
#0.498


Now divide the data into training and test data. We’ll run a logistic regression on the training data, then see how well the model performs on the training data.

#Set the first 50 rows to train the model
X_train = X_predictors[:50]
Y_train = Y[:50]

#Set the remaining rows to test the model
X_test = X_predictors[50:]
Y_test = Y[50:]

print(f"X_train: {len(X_train)} X_test: {len(X_test)}")
#X_train: 50 X_test: 450

Logistic regression with Scikit-learn

We’re ready to train and test models.

As we train the models, we need to take steps to avoid overfitting. A machine learning model may have very accurate results with the data used to train the model. But this does not mean it will be equally accurate when making predictions with data it hasn’t seen before. When the model fails to generalize to new data, we say it has “overfit” the training data. Overfitting is more likely when there are few observations to train on, and when the model uses many correlated predictors.

How to avoid overfitting? By default, scikit-learn‘s logistic regression applies regularization. Regularization balances the need for predictive accuracy on the training data with a penalty on the magnitude of the model coefficients. Increasing the penalty reduces the coefficients and hence reduces the likelihood of overfitting. If the penalty is too large, though, it will reduce predictive power on both the training and test data.

from sklearn.linear_model import LogisticRegression
scikit_default = LogisticRegression(random_state=0).fit(X_train, Y_train)
print(f"intecept: {scikit_default.intercept_} coeficients: {scikit_default.coef_}")
print(f"train accuracy: {scikit_default.score(X_train, Y_train)}")
print(f"test accuracy: {scikit_default.score(X_test, Y_test)}")
"""
Results will vary slightly, even when you set random_state.
intecept: [-0.44526823] coeficients: [[0.50031563 0.79636504 0.82047214 0.83635656]]
train accuracy: 0.8
test accuracy: 0.8088888888888889
"""

We can set turn off regularization by setting penalty as none. Applying regularization reduces the magnitude of the coefficients. Setting the penalty to none will increase the coefficients. Notice that the accuracy on the test data decreases. This indicates our model has overfit the training data.

from sklearn.linear_model import LogisticRegression
scikit_no_penalty = LogisticRegression(random_state=0,penalty='none').fit(X_train, Y_train)
print(f"intecept: {scikit_no_penalty.intercept_} coeficients: {scikit_no_penalty.coef_}")
print(f"train accuracy: {scikit_no_penalty.score(X_train, Y_train)}")
print(f"test accuracy: {scikit_no_penalty.score(X_test, Y_test)}")
"""
intecept: [-0.63388911] coeficients: [[-3.59878438  0.70813119  5.10660019  1.29684873]]
train accuracy: 0.82
test accuracy: 0.7888888888888889
"""


C is 1.0 by default. Smaller values of C increase the regularization, so if we set the value to .1 we reduce the magnitude of the coefficients.

from sklearn.linear_model import LogisticRegression
scikit_bigger_penalty = LogisticRegression(random_state=0,C=.1).fit(X_train, Y_train)
print(f"intecept: {scikit_bigger_penalty.intercept_} \
    coeficients: {scikit_bigger_penalty.coef_}")
print(f"train accuracy: {scikit_bigger_penalty.score(X_train, Y_train)}")
print(f"test accuracy: {scikit_bigger_penalty.score(X_test, Y_test)}")
"""
intecept: [-0.13102803]     coeficients: [[0.3021235  0.3919277  0.34359251 0.40332636]]
train accuracy: 0.8
test accuracy: 0.8066666666666666
"""


It’s nice to be able to adjust the smoothing coefficient, but how do we decide the optimal value? Scikit-learn’s GridSearchCV provides an effective but easy to use method for choosing an optimal value. The “Grid Search” in GridSearchCV means that we supply a dictionary with the parameter values we wish to test. The model is fit with all combinations of those values. If we have 4 possible values for C and 2 possible values for solver, we will search through all 4X2=8 combinations.

GridSearchCV Searches Through This Grid

The “CV” in GridSearchCV stands for cross-validation. Cross-validation is the method of segmenting the training data. The model is trained on all but one of the segments and the remaining segment validate the model.

GridSearch and cross-validation work in combination. GridsearchCV iterates through values of C and solver for different test and training segments. The algorithm selects the best estimator based performance on the validation segments.

Doing this allows us to determine which values of C and solver work best for our training data. This is how scikit-learn helps us to optimize predictive accuracy.

Let’s see it in action.

from sklearn.model_selection import GridSearchCV
parameters = {'C':[.01, .1, 1, 10],'solver':['newton-cg','lbfgs']}
Logistic = LogisticRegression(random_state=0)
scikit_GridSearchCV = GridSearchCV(Logistic, parameters)
scikit_GridSearchCV.fit(X_train, Y_train)
print(f"best estimator: {scikit_GridSearchCV.best_estimator_}")
#best estimator: LogisticRegression(C=0.1, random_state=0, solver='newton-cg')

Use the score method returns the mean accuracy on the given test data and labels. Accuracy is the percent of observations correctly predicted.

print(f"train accuracy: {scikit_GridSearchCV.score(X_train, Y_train)}")
print(f"test accuracy: {scikit_GridSearchCV.score(X_test, Y_test)}")
"""
train accuracy: 0.82
test accuracy: 0.8133333333333334
"""

Logistic regression with Statsmodels

Now let’s try the same, but with statsmodels. With scikit-learn, to turn off regularization we set penalty='none', but with statsmodels regularization is turned off by default. A quirk to watch out for is that Statsmodels does not include an intercept by default. To include an intercept, we use the sm.add_constant method.

import statsmodels.api as sm

#adding constant to X
X_train_with_constant = sm.add_constant(X_train)
X_test_with_constant = sm.add_constant(X_test)

# building the model and fitting the data
sm_model_all_predictors = sm.Logit(Y_train, X_train_with_constant).fit()

# printing the summary table
print(sm_model_all_predictors.params)
"""
Optimization terminated successfully.
         Current function value: 0.446973
         Iterations 7
[-0.57361523 -2.00207425  1.28872367  3.53734636  0.77494424]
"""

If you’re used to doing logistic regression in R or SAS, what comes next will be familiar. Once we have trained the logistic regression model with statsmodels, the summary method will easily produce a table with statistical measures including p-values and confidence intervals.

sm_model_all_predictors.summary()

There’s a lot here, but we’ll focus on the second table with the coefficients.

The first column shows the value for the coefficient. The fourth column, with the heading P>|z|, shows the p-values. A p-value is a probability measure, and p-values above .05 are frequently considered, “not statistically significant.” None of the predictors are considered statistically significant! This is because we have a relatively small number of observations in our training data and because the predictors are highly correlated. Some statistical packages like R and SAS have built-in methods to select the features to include in the model based on which predictors have low (significant) p-values, but unfortunately, this isn’t available in statsmodels.

If we try again with just x1 and x2, we’ll get a completely different result, with very low p-values for x1 and x2, meaning that the evidence for a relationship with the dependent variable is statistically significant. We’re cheating, though – because we created the data, we know that we only need x1 and x2.

sm_model_x1_x2 = sm.Logit(Y_train, X_train_with_constant[:,:3]).fit()
sm_model_x1_x2.summary()

Now we see x1 and x2 are both statistically significant.

Statsmodels doesn’t have the same accuracy method that we have in scikit-learn. We’ll use the predict method to predict the probabilities. Then we’ll use the decision rule that probabilities above .5 are true and all others are false. This is the same rule used when scikit-learn calculates accuracy.

all_predicted_train = sm_model_all_predictors.predict(X_train_with_constant)>.5
all_predicted_test = sm_model_all_predictors.predict(X_test_with_constant)>.5

x1_x2_predicted_train = sm_model_x1_x2.predict(X_train_with_constant[:,:3])>.5
x1_x2_predicted_test = sm_model_x1_x2.predict(X_test_with_constant[:,:3])>.5

#calculate the accuracy
print(f"train: {(Y_train==all_predicted_train).mean()} and test: {(Y_test==all_predicted_test).mean()}")
print(f"train: {(Y_train==x1_x2_predicted_train).mean()} and test: {(Y_test==x1_x2_predicted_test).mean()}")
"""
train: 0.8 and test: 0.8066666666666666
train: 0.8 and test: 0.8111111111111111
"""

Summarizing The Results

Let’s create a DataFrame with the results. The models have identical accuracy on the training data, but different results on the test data. The models with all the predictors and without smoothing have the worst test accuracy, suggesting that they have overfit on the training data and so do not generalize well to new data.

Even if we use the best methods in creating our model, there is still chance involved in how well it generalizes to the test data.

lst = [['scikit-learn','default', scikit_default.score(X_train, Y_train),scikit_default.score(X_test, Y_test)],
       ['scikit-learn','no penalty', scikit_no_penalty.score(X_train, Y_train),scikit_no_penalty.score(X_test, Y_test)],
       ['scikit-learn','bigger penalty', scikit_bigger_penalty.score(X_train, Y_train),scikit_bigger_penalty.score(X_test, Y_test)],
       ['scikit-learn','GridSearchCV', scikit_GridSearchCV.score(X_train, Y_train),scikit_GridSearchCV.score(X_test, Y_test)],
       ['statsmodels','include intercept and all predictors', (Y_train==all_predicted_train).mean(),(Y_test==all_predicted_test).mean()],
       ['statsmodels','include intercept and x1 and x2', (Y_train==x1_x2_predicted_train).mean(),(Y_test==x1_x2_predicted_test).mean()]
      ]
df = pd.DataFrame(lst, columns =['package', 'setting','train accuracy','test accuracy'])
df

Scikit-learn vs Statsmodels

Upshot is that you should use Scikit-learn for logistic regression unless you need the statistics results provided by StatsModels.

Here’s a table of the most relevant similarities and differences:

That’s it for now! Please check out my other work at learningtableau.com and my new site datasciencedrills.com.