Should I use statsmodels or scikit-learn for regression?

The short version

They serve different goals. statsmodels is built for inference: it reports coefficients with standard errors, p-values, confidence intervals and diagnostic tests, so you can explain a relationship. scikit-learn is built for prediction: it gives pipelines, cross-validation, hyperparameter search and a consistent interface across many models. Choose by the question you are answering.

Statsmodels for explaining

import statsmodels.formula.api as smf

model = smf.ols("sales ~ temp", data=df).fit()
print(model.summary())

The summary gives you everything you need to talk about the relationship: the coefficient, whether it is distinguishable from zero, the plausible range, and how well the model fits. If your output is a sentence like "each extra degree is associated with 1.3 more sales, plus or minus 0.2", this is the right tool.

Scikit-learn for predicting

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

model = LinearRegression()
scores = cross_val_score(model, X, y, cv=5)

Scikit-learn deliberately does not give p-values, because its purpose is to make good predictions on unseen data, measured by held-out performance. If your output is a number predicted for a new row, and you care about how accurate it is, this is the right tool.

The tell-tale difference

Scikit-learn's LinearRegression has no summary() and no p-values, which surprises people coming from R or Stata. That is not an omission; it reflects a different philosophy. Statistical inference requires assumptions about the data-generating process, and scikit-learn's model-agnostic, cross-validated approach deliberately avoids them.

The intercept trap in statsmodels

The formula interface adds an intercept automatically, but the array interface does not.

import statsmodels.api as sm
X = sm.add_constant(X)      # required, or you fit a line through the origin
sm.OLS(y, X).fit()

Forgetting add_constant silently fits a model with no intercept, which usually gives badly wrong coefficients. The formula interface (smf.ols) avoids this, which is one reason to prefer it.

Using both

There is no conflict in using both on the same problem. Fit with statsmodels to understand which predictors matter and how, then build a scikit-learn pipeline to make and validate predictions. They answer complementary questions, and a good analysis often needs both.

A caution about p-values in prediction models

If you select variables by their p-values and then report performance on the same data, you have overfitted and your p-values are no longer meaningful. Keep the two activities separate: use inference to understand, use held-out data to measure prediction, and do not let one contaminate the other.