statsmodels
statsmodels is what you reach for when you need to explain a relationship rather than predict a label. It fits classical statistical models (OLS and GLS regression, logit and probit, Poisson and negative binomial GLMs, ARIMA and full state space time series, survival models, mixed effects) and then prints the things scikit-learn deliberately withholds: coefficients with standard errors, t-statistics, p-values, confidence intervals, R-squared, AIC, and a long list of specification and diagnostic tests. There are two front doors, an array-based API in the spirit of Stata and a formula API using R's y ~ x1 + x2 syntax, and the same engine powers seaborn's regression plots and most econometrics coursework written in Python.
The only serious home for classical statistics and econometrics in Python, and worth the install the moment anyone needs a standard error rather than a prediction. Budget a day for the API's quirks: no automatic intercept, silent NaN propagation, and a 0.x version number that has never matched how load-bearing the library actually is.
Use it if
- Someone will ask you whether a coefficient is significant, how wide the confidence interval is, or whether the residuals violate an assumption; the summary() table answers all three and scikit-learn answers none of them
- You are porting analysis from R or Stata and want the formula syntax, the same default contrasts for categorical variables, and comparable output tables
- You need classical time series done properly: ARIMA and SARIMAX, VAR and VECM, state space models with Kalman filtering, Markov switching, and STL or seasonal decomposition
- You need the statistics grab bag around the models: multiple-testing corrections, variance inflation factors, Breusch-Pagan and Durbin-Watson, ANOVA tables, power analysis, and cluster or heteroskedasticity-consistent covariance estimators
- Your goal is prediction accuracy on a large dataset: statsmodels fits in memory, has no pipeline or cross-validation machinery, and no regularization path search. scikit-learn does that job and does it faster
- You only need a handful of hypothesis tests: scipy.stats gives you t-tests, chi-square, Mann-Whitney and correlation without adding a compiled 30MB dependency and its numpy, scipy, pandas and patsy chain
- You want panel data or instrumental variables as first-class citizens: statsmodels can be pushed there, but linearmodels was written for fixed effects, random effects, IV/2SLS and system estimation and gives cleaner output
- Deep Bayesian work: outside the Bayesian mixed GLM there is no MCMC story here, and PyMC or a Stan interface will be less painful than bending frequentist models into a posterior
- You dislike surprises from a 0.x project: the package has been on 0.x for well over a decade, still ships the Development Status :: 4 - Beta classifier, and does remove things across minors (the old statsmodels.tsa.arima_model.ARIMA was dropped in favour of statsmodels.tsa.arima.model.ARIMA)
- You need fast issue turnaround: 2673 open issues (2881 counting PRs) against a small core team, with releases arriving roughly twice a year
Setup reality
pip install statsmodels pulls numpy, scipy, pandas, patsy and packaging, so a slim container gains a few hundred megabytes. Wheels cover the common platforms, but the package contains compiled Cython extensions, so the first weeks after a new CPython release usually mean a source build with Cython and a C compiler. Two behaviours catch nearly everyone on day one. First, sm.OLS does not add an intercept: fit a model without sm.add_constant and you silently get a regression through the origin with a flattering R-squared. The formula API adds one for you, which makes the two APIs disagree on the same data. Second, the default missing='none' means NaNs are not checked, so a single missing value turns the entire results table into NaN with no error; pass missing='drop' or clean the frame first. Finally, patsy handles formula parsing and is itself in low-activity maintenance, so exotic formula transforms are unlikely to gain features.
Patterns
Fit OLS with the array APIfit-ols-arrays
import numpy as np
import statsmodels.api as sm
X = sm.add_constant(df[["income", "age"]]) # REQUIRED for an intercept
y = df["spend"]
model = sm.OLS(y, X, missing="drop")
results = model.fit()
print(results.summary())
print(results.params["income"], results.pvalues["income"])Without sm.add_constant you fit a line through the origin and the reported R-squared is computed differently, which usually makes a bad model look excellent. missing='drop' matters too: the default does not check for NaN and the whole table comes back NaN.
Fit the same model with R-style formulasfit-ols-formula
import statsmodels.formula.api as smf
results = smf.ols(
"spend ~ income + age + C(region) + income:age",
data=df,
).fit()
print(results.summary())The formula API adds the intercept automatically, so the two APIs give different results on identical data unless you are careful. C(region) creates treatment-coded dummies dropping the first level; use C(region, Treatment(reference='north')) to choose the baseline.
Logistic regression with odds ratios and marginal effectslogistic-regression
import numpy as np
import statsmodels.formula.api as smf
fit = smf.logit("converted ~ visits + C(plan)", data=df).fit()
print(fit.summary())
odds_ratios = np.exp(fit.params)
ci = np.exp(fit.conf_int())
print(fit.get_margeff(at="mean").summary())fit.predict() returns probabilities, not 0/1 labels; threshold them yourself. Perfect separation raises a PerfectSeparationError or produces enormous coefficients with enormous standard errors, which is the model telling you a predictor encodes the outcome.
Heteroskedasticity-consistent and clustered standard errorsadjust-standard-errors
hc = smf.ols("spend ~ income", data=df).fit(cov_type="HC3")
clustered = smf.ols("spend ~ income", data=df).fit(
cov_type="cluster",
cov_kwds={"groups": df["store_id"]},
)
print(clustered.bse)Only the standard errors, t-statistics and p-values change; the coefficients are identical. Clustered errors need a decent number of clusters (rules of thumb start around 30 to 50) or they are biased downward and overstate significance.
Predict new data with confidence and prediction intervalspredict-with-intervals
new = pd.DataFrame({"income": [50_000, 90_000], "age": [30, 45]})
pred = results.get_prediction(new)
frame = pred.summary_frame(alpha=0.05)
print(frame[["mean", "mean_ci_lower", "mean_ci_upper",
"obs_ci_lower", "obs_ci_upper"]])mean_ci is the interval for the average response and obs_ci is the interval for a single new observation; the second is always wider and is the one people actually want when forecasting. With the array API you must add the constant column to new data too.
Count data with a Poisson GLM and an offsetpoisson-glm
import statsmodels.api as sm
import numpy as np
X = sm.add_constant(df[["promo", "weekend"]])
fit = sm.GLM(
df["orders"],
X,
family=sm.families.Poisson(),
offset=np.log(df["visits"]),
).fit()
print(fit.summary())
print("dispersion:", fit.pearson_chi2 / fit.df_resid)The offset must be on the log scale because that is the link function. If the dispersion ratio is far above 1 the Poisson standard errors are too small; switch to sm.families.NegativeBinomial or refit with cov_type='HC0'.
Fit ARIMA or SARIMAX and forecastarima-forecast
from statsmodels.tsa.arima.model import ARIMA
series = df.set_index("date")["sales"].asfreq("D")
fit = ARIMA(series, order=(2, 1, 2), seasonal_order=(1, 0, 1, 7)).fit()
print(fit.summary())
fc = fit.get_forecast(steps=14)
print(fc.predicted_mean)
print(fc.conf_int())Import from statsmodels.tsa.arima.model; the old statsmodels.tsa.arima_model path was removed and most tutorials online still use it. The index needs a frequency (asfreq or a proper DatetimeIndex) or forecasts come back with integer positions instead of dates.
Split a series into trend, season and remainderdecompose-seasonality
from statsmodels.tsa.seasonal import STL, seasonal_decompose
stl = STL(series, period=7, seasonal=13).fit()
trend, season, resid = stl.trend, stl.seasonal, stl.resid
# quick classical version
classic = seasonal_decompose(series, model="additive", period=7)STL handles a changing seasonal shape and outliers better than seasonal_decompose, which just uses moving averages. The seasonal argument must be an odd number of at least 7, and period is in observations, not days.
Check regression assumptionsdiagnostic-tests
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.stattools import durbin_watson
lm, lm_p, f, f_p = het_breuschpagan(results.resid, results.model.exog)
print("Breusch-Pagan p:", lm_p)
print("Durbin-Watson:", durbin_watson(results.resid))
vifs = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]variance_inflation_factor expects the design matrix including the constant column, and the VIF for that constant is meaningless, so ignore index 0. A VIF above roughly 5 to 10 is the usual flag for collinearity.
Test a series for stationarityunit-root-test
from statsmodels.tsa.stattools import adfuller, kpss
stat, pvalue, lags, nobs, crit, icbest = adfuller(series.dropna())
print(pvalue, crit["5%"]) # p < 0.05 rejects a unit root
kpss_stat, kpss_p, kpss_lags, kpss_crit = kpss(series.dropna(), regression="c")The two tests have opposite null hypotheses: ADF's null is non-stationary, KPSS's null is stationary. Run both, and when they disagree the series usually needs differencing or detrending rather than a verdict.
Fit a random-intercept mixed modelmixed-effects-model
import statsmodels.formula.api as smf
fit = smf.mixedlm(
"score ~ hours + C(method)",
data=df,
groups=df["school_id"],
re_formula="~hours", # random slope for hours too
).fit()
print(fit.summary())
print(fit.random_effects["school_3"])Convergence warnings are common with random slopes; drop re_formula back to a random intercept before believing anything. The fit uses REML by default, so likelihood ratio tests comparing different fixed-effect structures need reml=False.
Get results out as a DataFrameresults-to-dataframe
coefs = results.summary2().tables[1] # DataFrame of estimates
print(coefs.columns.tolist())
tidy = pd.DataFrame({
"coef": results.params,
"se": results.bse,
"p": results.pvalues,
}).join(results.conf_int().rename(columns={0: "lo", 1: "hi"}))summary() returns a printable object, not data; parsing its text is a trap people fall into. summary2().tables[1] is a real DataFrame, and building one from .params, .bse and .pvalues is more stable across versions.
Alternatives
| Package | Registry | Pick it when |
|---|---|---|
| scikit-learn | PyPI | You care about out-of-sample prediction and want pipelines, cross-validation and regularization rather than p-values. |
| scipy | PyPI | You need a few standalone hypothesis tests or distributions and no model fitting at all. |
| linearmodels | PyPI | Your data is panel or you need instrumental variables, 2SLS, GMM or system estimation as the primary use case. |
| pymc | PyPI | You want posterior distributions and priors instead of point estimates and p-values. |