I could not find a way to stepwise regression in scikit learn. I have checked all other posts on Stack Exchange on this topic. Answers to all of them suggests using f_regression.

But f_regression does not do stepwise regression but only give F-score and pvalues corresponding to each of the regressors, which is only the first step in stepwise regression.

What to do after 1st regressors with the best f-score is chosen?


Scikit-learn indeed does not support stepwise regression. That's because what is commonly known as 'stepwise regression' is an algorithm based on p-values of coefficients of linear regression, and scikit-learn deliberately avoids inferential approach to model learning (significance testing etc). Moreover, pure OLS is only one of numerous regression algorithms, and from the scikit-learn point of view it is neither very important, nor one of the best.

There are, however, some pieces of advice for those who still need a good way for feature selection with linear models:

  1. Use inherently sparse models like ElasticNet or Lasso.
  2. Normalize your features with StandardScaler, and then order your features just by model.coef_. For perfectly independent covariates it is equivalent to sorting by p-values. The class sklearn.feature_selection.RFE will do it for you, and RFECV will even evaluate the optimal number of features.
  3. Use an implementation of forward selection by adjusted $R^2$ that works with statsmodels.
  4. Do brute-force forward or backward selection to maximize your favorite metric on cross-validation (it could take approximately quadratic time in number of covariates). A scikit-learn compatible mlxtend package supports this approach for any estimator and any metric.
  5. If you still want vanilla stepwise regression, it is easier to base it on statsmodels, since this package calculates p-values for you. A basic forward-backward selection could look like this:


from sklearn.datasets import load_boston
import pandas as pd
import numpy as np
import statsmodels.api as sm

data = load_boston()
X = pd.DataFrame(data.data, columns=data.feature_names)
y = data.target

def stepwise_selection(X, y, 
                       threshold_out = 0.05, 
    """ Perform a forward-backward feature selection 
    based on p-value from statsmodels.api.OLS
        X - pandas.DataFrame with candidate features
        y - list-like with the target
        initial_list - list of features to start with (column names of X)
        threshold_in - include a feature if its p-value < threshold_in
        threshold_out - exclude a feature if its p-value > threshold_out
        verbose - whether to print the sequence of inclusions and exclusions
    Returns: list of selected features 
    Always set threshold_in < threshold_out to avoid infinite looping.
    See https://en.wikipedia.org/wiki/Stepwise_regression for the details
    included = list(initial_list)
    while True:
        # forward step
        excluded = list(set(X.columns)-set(included))
        new_pval = pd.Series(index=excluded)
        for new_column in excluded:
            model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included+[new_column]]))).fit()
            new_pval[new_column] = model.pvalues[new_column]
        best_pval = new_pval.min()
        if best_pval < threshold_in:
            best_feature = new_pval.argmin()
            if verbose:
                print('Add  {:30} with p-value {:.6}'.format(best_feature, best_pval))

        # backward step
        model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
        # use all coefs except intercept
        pvalues = model.pvalues.iloc[1:]
        worst_pval = pvalues.max() # null if pvalues is empty
        if worst_pval > threshold_out:
            worst_feature = pvalues.argmax()
            if verbose:
                print('Drop {:30} with p-value {:.6}'.format(worst_feature, worst_pval))
        if not changed:
    return included

result = stepwise_selection(X, y)

print('resulting features:')

This example would print the following output:

Add  LSTAT                          with p-value 5.0811e-88
Add  RM                             with p-value 3.47226e-27
Add  PTRATIO                        with p-value 1.64466e-14
Add  DIS                            with p-value 1.66847e-05
Add  NOX                            with p-value 5.48815e-08
Add  CHAS                           with p-value 0.000265473
Add  B                              with p-value 0.000771946
Add  ZN                             with p-value 0.00465162
resulting features:
['LSTAT', 'RM', 'PTRATIO', 'DIS', 'NOX', 'CHAS', 'B', 'ZN']
  • $\begingroup$ Does this mean that the scikit-learn point imply that p-values are useless? $\endgroup$ – Digio Jul 3 '18 at 13:36
  • $\begingroup$ "pure OLS is only one of numerous regression algorithms, and from the scikit-learn point of view it is neither very important, nor one of the best." I don't understand. OLS solves finds a closed-form unique solution to a convex problem. How can any other algorithm perform better than one which is already at the global optimum? $\endgroup$ – Digio Jul 4 '18 at 8:12
  • $\begingroup$ OLS minimizes MSE of a linear model on the train set. Other algorithms may: (1) use various regularizations, which increase MSE on training data, but hope to improve generalizing ability - such as Lasso, Ridge, or bayessian linear regression; (2) minimize other losses instead of MSE - e.g. MAE or Huber loss; (3) use a non-linear model, e.g. ensemble of decision trees, or a neural network. $\endgroup$ – David Dale Jul 4 '18 at 10:52
  • $\begingroup$ As far as I understand, p-values (1) are a very specific interpretation of a single OLS algorithm, and (2) are useful for inference (to decide whether a single predictor matters), but not so useful for prediction (model with lots of bad p-values may have good predictive power, and vice versa) $\endgroup$ – David Dale Jul 4 '18 at 10:56
  • $\begingroup$ Stating that OLS is just not good enough compared to other methods is misleading. For a linearly separable dataset where the Gauss-Markov assumptions are satisfied, OLS will be more efficient than any other linear or nonlinear method. It's more of a question of data and model structure than anything else. $\endgroup$ – Digio Jul 5 '18 at 10:56

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