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.

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:

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)

[-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()

print(results.summary())

                     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:

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)

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

Common Patsy building blocks:

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:

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

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

  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)"]

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:

Check your understanding

🧪 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

Press Run to see what your code prints.

Exercise 2 – Predict a new observation

Press Run to see what your code prints.

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

Press Run to see what your code prints.

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.