Linear Regression with statsmodels

What you’ll learn

• What statsmodels is for, and how it differs from a prediction-only library like scikit-learn
• Fitting OLS with the array API: sm.add_constant + sm.OLS(y, X).fit()
• Reading .summary() — coefficients, std err, t, p-values, and R-squared
• The formula API smf.ols("y ~ x1 + x2", data=df).fit() and basic Patsy formula syntax
• Predicting on new data with .predict()

What statsmodels is for

statsmodels is a Python library for fitting statistical models, running statistical tests, and exploring data. Where a library like scikit-learn is built to answer “what will this predict?”, statsmodels is built to answer “why, and how confident are we?” — it reports standard errors, t-statistics, p-values, and confidence intervals alongside every coefficient, not just the fitted numbers.

It covers a lot of ground: ordinary least squares and other linear models, ANOVA, time series models, and more. This page focuses on the workhorse of the bunch — OLS (ordinary least squares) linear regression — using both ways statsmodels lets you describe a model: plain arrays, and R-style text formulas.

statsmodels imports

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

statsmodels.api (sm) is the array-based interface. statsmodels.formula.api (smf) is the formula-based interface, powered by Patsy under the hood.

The array API: sm.OLS

A linear model almost always needs an intercept term — a column of 1s so the fitted line doesn’t have to pass through the origin. sm.add_constant adds it for you:

ols_array_api.py

import numpy as np
import statsmodels.api as sm
 
rng = np.random.default_rng(seed=12345)
 
def dnorm(mean, variance, size):
    return mean + np.sqrt(variance) * rng.standard_normal(size)
 
N = 100
X = np.column_stack([
    dnorm(0, 0.4, N),
    dnorm(0, 0.6, N),
    dnorm(0, 0.2, N),
])
eps = dnorm(0, 0.1, N)
beta = [0.1, 0.3, 0.5]   # the "true" coefficients we're trying to recover
y = X @ beta + eps
 
X_model = sm.add_constant(X)     # prepend a column of 1s for the intercept
model = sm.OLS(y, X_model)
results = model.fit()
 
print(results.params)

Output

[-0.0208  0.0658  0.269   0.4494]

model.fit() returns a results object, not just a plain array of numbers. That object is where statsmodels’ inference focus shows up — it’s carrying around standard errors, residuals, and everything .summary() needs, in addition to results.params.

Reading .summary()

ols_summary.py

print(results.summary())

Output (trimmed)

                     coef    std err       t      P>|t|      [0.025   0.975]
------------------------------------------------------------------------------
const             -0.0208      0.032    -0.653     0.516     -0.084    0.042
x1                 0.0658      0.054     1.220     0.226     -0.041    0.173
x2                 0.2690      0.043     6.312     0.000      0.184    0.354
x3                 0.4494      0.068     6.567     0.000      0.314    0.585
==============================================================================
R-squared:                       0.470

A few things worth reading every time:

• coef — the estimated slope for that variable, holding the others fixed.
• std err — how much that coefficient would jitter if you re-sampled the data.
• t and P>|t| — the t-statistic and its p-value, testing “is this coefficient really different from 0?” Small p-values (x2, x3 above) mean strong evidence the variable matters; a large one (const, x1) means we can’t rule out that its true effect is 0.
• R-squared — the fraction of variance in y explained by the model, from 0 (no fit) to 1 (perfect fit).

The formula API: smf.ols and Patsy

Patsy formulas describe a model as a short string instead of hand-built arrays. The basic shape is "y ~ x1 + x2" — read the ~ as “explained by” and + as “and also include this term” (not arithmetic addition!). When you use the formula API with a DataFrame, statsmodels adds the intercept and attaches real column names automatically — no add_constant needed:

ols_formula_api.py

import pandas as pd
import statsmodels.formula.api as smf
 
data = pd.DataFrame(X, columns=["col0", "col1", "col2"])
data["y"] = y
 
results = smf.ols("y ~ col0 + col1 + col2", data=data).fit()
print(results.params)

Output

Intercept    -0.020799
col0          0.065813
col1          0.268970
col2          0.449419
dtype: float64

Common Patsy building blocks:

Syntax	Meaning
y ~ x1 + x2	fit x1 and x2 as separate terms, intercept included by default
y ~ x1 + x2 + 0	drop the intercept
y ~ x1 * x2	shorthand for x1 + x2 + x1:x2 (both main effects plus their interaction)
y ~ x1 : x2	just the interaction term between x1 and x2
y ~ C(category_col)	treat a numeric column as categorical, expanding it into dummy variables

Non-numeric columns are turned into dummy variables automatically, and C(...) forces the same treatment for a column that’s stored as numbers but is really a category (like a store ID or a survey rating).

Predicting with .predict()

Once a model is fit, .predict() scores new rows using the fitted coefficients:

ols_predict.py

new_rows = data[["col0", "col1", "col2"]].head(3)
print(results.predict(new_rows))

Output

0   -0.592959
1   -0.531160
2    0.058636
dtype: float64

statsmodels vs. scikit-learn

Both fit models to data, but they’re built for different jobs:

• statsmodels — you care about why: coefficients, p-values, confidence intervals, formal statistical tests. Slower to fit, richer to interpret.
• scikit-learn — you care about predicting well: a consistent .fit() / .predict() API across dozens of algorithms (including ones with no p-values at all, like random forests), built-in cross-validation, and pipelines.

A common workflow: explore and validate with statsmodels, then move to scikit-learn once you’ve settled on which features matter and just need fast, repeatable predictions.

From data to inference

diagram OLS fitting pipeline mermaid

  flowchart LR
  A["Raw data
(DataFrame or arrays)"] --> B["Design matrix
(add_constant, or a Patsy formula)"]
  B --> C["sm.OLS(y, X).fit()
or smf.ols(formula, data).fit()"]
  C --> D["results.summary()
(coef, std err, t, p-values, R-squared)"]
  C --> E["results.predict(new_data)"]

How a DataFrame becomes a design matrix, a fitted model, and finally a readable summary or a prediction.

Visualize it

An OLS fit is the line that makes the residuals — the vertical gaps between each point and the line — as small as possible overall. Watch how the fitted line sits right through the middle of the cloud, and how each residual is just the leftover distance the line didn’t explain:

sketch OLS fit and residuals p5.js Pause Restart

The fitted line minimizes the sum of squared vertical distances (residuals) from every point.

From the book

McKinney is careful to point out that a + b inside a Patsy formula does not mean addition — it means “these are two separate terms in the design matrix.” If you actually want to add two columns together as one term, wrap them in I(...), like y ~ I(x0 + x1). He also notes the formula API’s biggest convenience: when you fit against a pandas DataFrame with smf.ols, statsmodels adds the intercept and attaches the real column names to results.params for you — no manual add_constant or bookkeeping required.

🧪 Try It Yourself

The exercises below fit OLS by hand with numpy.linalg.lstsq — the same least-squares math statsmodels runs under the hood — so you can see exactly where those coefficients come from.

Exercise 1 – Build a design matrix and fit OLS

pch.exerciseTag pch.exerciseTitle

Exercise 2 – Predict a new observation

pch.exerciseTag pch.exerciseTitle

Exercise 3 – R-squared: how much variance is explained

pch.exerciseTag pch.exerciseTitle

Next

That wraps up the statistics module. Continue to Introduction to SQL for Data Analytics (Phase 9) to start pulling and aggregating the data you’ll analyze.