This appendix is a short tour of the formula interface to statsmodels, which lets you write a regression model as a string such as 'y ~ x1 + x2' instead of assembling a design matrix by hand.
In the regression lectures we built the design matrix \(X\) ourselves: pick the columns, add a column of ones for the intercept, one-hot encode anything categorical, apply any transformations, and only then call sm.OLS(y, X).
The statsmodels.formula.api sub-module (conventionally imported as smf) lets us do all of that by typing a formula instead. A formula can
smf is built on the patsy package, a Python implementation of the formula mini-language used in R and S. If you have seen lm(y ~ x, data) in R, you already know the syntax.
We start with the Longley macroeconomic data set used in the regularization lecture: 16 years of US data, predicting total employment (TOTEMP) from five economic indicators. With the formula interface the response and the predictors come from the same data frame, so we put TOTEMP in as a column.
from statsmodels.datasets.longley import load_pandas
longley = load_pandas().exog
longley['TOTEMP'] = load_pandas().endog
longley.head()| GNPDEFL | GNP | UNEMP | ARMED | POP | YEAR | TOTEMP | |
|---|---|---|---|---|---|---|---|
| 0 | 83.0 | 234289.0 | 2356.0 | 1590.0 | 107608.0 | 1947.0 | 60323.0 |
| 1 | 88.5 | 259426.0 | 2325.0 | 1456.0 | 108632.0 | 1948.0 | 61122.0 |
| 2 | 88.2 | 258054.0 | 3682.0 | 1616.0 | 109773.0 | 1949.0 | 60171.0 |
| 3 | 89.5 | 284599.0 | 3351.0 | 1650.0 | 110929.0 | 1950.0 | 61187.0 |
| 4 | 96.2 | 328975.0 | 2099.0 | 3099.0 | 112075.0 | 1951.0 | 63221.0 |
mod = smf.ols(formula='TOTEMP ~ GNPDEFL + GNP + UNEMP + ARMED + POP', data=longley)
res = mod.fit()
with warnings.catch_warnings():
warnings.simplefilter('ignore')
print(res.summary()) OLS Regression Results
==============================================================================
Dep. Variable: TOTEMP R-squared: 0.987
Model: OLS Adj. R-squared: 0.981
Method: Least Squares F-statistic: 156.4
Date: Thu, 15 Oct 2026 Prob (F-statistic): 3.70e-09
Time: 21:51:33 Log-Likelihood: -117.83
No. Observations: 16 AIC: 247.7
Df Residuals: 10 BIC: 252.3
Df Model: 5
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 9.246e+04 3.52e+04 2.629 0.025 1.41e+04 1.71e+05
GNPDEFL -48.4628 132.248 -0.366 0.722 -343.129 246.204
GNP 0.0720 0.032 2.269 0.047 0.001 0.143
UNEMP -0.4039 0.439 -0.921 0.379 -1.381 0.573
ARMED -0.5605 0.284 -1.975 0.077 -1.193 0.072
POP -0.4035 0.330 -1.222 0.250 -1.139 0.332
==============================================================================
Omnibus: 1.572 Durbin-Watson: 1.248
Prob(Omnibus): 0.456 Jarque-Bera (JB): 0.642
Skew: 0.489 Prob(JB): 0.725
Kurtosis: 3.079 Cond. No. 1.21e+08
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.21e+08. This might indicate that there are
strong multicollinearity or other numerical problems.The formula
'TOTEMP ~ GNPDEFL + GNP + UNEMP + ARMED + POP'
is an R-style formula string that specifies the model.
~ is the dependent variable, here TOTEMP.~ separates the dependent variable from the independent variables.GNPDEFL + GNP + UNEMP + ARMED + POP are the independent variables; the + means “include this column”, not arithmetic addition.Intercept row in the table. We never had to add a column of ones.The large condition number in the notes under the table is the multicollinearity warning discussed in the lecture.
To fit a smaller model, list fewer variables. To fit a model without an intercept, add - 1 to the formula: the - sign removes a term, and 1 stands for the intercept.
mod = smf.ols(formula='TOTEMP ~ GNPDEFL + GNP + UNEMP - 1', data=longley)
res = mod.fit()
with warnings.catch_warnings():
warnings.simplefilter('ignore')
print(res.summary()) OLS Regression Results
=======================================================================================
Dep. Variable: TOTEMP R-squared (uncentered): 1.000
Model: OLS Adj. R-squared (uncentered): 1.000
Method: Least Squares F-statistic: 1.127e+04
Date: Thu, 15 Oct 2026 Prob (F-statistic): 1.92e-22
Time: 21:51:33 Log-Likelihood: -137.20
No. Observations: 16 AIC: 280.4
Df Residuals: 13 BIC: 282.7
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
GNPDEFL 871.0961 25.984 33.525 0.000 814.961 927.231
GNP -0.0532 0.007 -8.139 0.000 -0.067 -0.039
UNEMP -0.8333 0.496 -1.679 0.117 -1.905 0.239
==============================================================================
Omnibus: 0.046 Durbin-Watson: 1.422
Prob(Omnibus): 0.977 Jarque-Bera (JB): 0.274
Skew: -0.010 Prob(JB): 0.872
Kurtosis: 2.359 Cond. No. 2.92e+04
==============================================================================
Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 2.92e+04. This might indicate that there are
strong multicollinearity or other numerical problems.The dependent variable is still TOTEMP; the independent variables are now only GNPDEFL, GNP and UNEMP; and there is no Intercept row. Forcing the fit through the origin is rarely what you want for economic data like this, but it is the standard way to get the design matrix \(X\) without a column of ones.
To look at categorical variables we switch to another famous data set.
The Guerry data set is a collection of historical data used in support of André-Michel Guerry’s 1833 Essay on the Moral Statistics of France.
André-Michel Guerry’s (1833) Essai sur la Statistique Morale de la France was one of the foundation studies of modern social science. Guerry assembled data on crimes, suicides, literacy and other “moral statistics,” and used tables and maps to analyze a variety of social issues in perhaps the first comprehensive study relating such variables.
Guerry’s results were startling for two reasons. First he showed that rates of crime and suicide remained remarkably stable over time, when broken down by age, sex, region of France and even season of the year; yet these numbers varied systematically across departements of France. This regularity of social numbers created the possibility to conceive, for the first time, that human actions in the social world were governed by social laws, just as inanimate objects were governed by laws of the physical world.
Source: “A.-M. Guerry’s Moral Statistics of France: Challenges for Multivariable Spatial Analysis”, Michael Friendly. Statistical Science 2007, Vol. 22, No. 3, 368–399.
We keep four columns: Lottery (per-capita wager on the Royal Lottery), Literacy, Wealth, and Region – a text label for one of five regions of France.
guerry = sm.datasets.get_rdataset("Guerry", "HistData").data
guerry = guerry[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()
guerry.head()| Lottery | Literacy | Wealth | Region | |
|---|---|---|---|---|
| 0 | 41 | 37 | 73 | E |
| 1 | 38 | 51 | 22 | N |
| 2 | 66 | 13 | 61 | C |
| 3 | 80 | 46 | 76 | E |
| 4 | 79 | 69 | 83 | E |
Here is a regression of Lottery on Literacy, Wealth, and Region.
mod = smf.ols(formula='Lottery ~ Literacy + Wealth + Region', data=guerry)
res = mod.fit()
print(res.summary()) OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.338
Model: OLS Adj. R-squared: 0.287
Method: Least Squares F-statistic: 6.636
Date: Thu, 15 Oct 2026 Prob (F-statistic): 1.07e-05
Time: 21:51:33 Log-Likelihood: -375.30
No. Observations: 85 AIC: 764.6
Df Residuals: 78 BIC: 781.7
Df Model: 6
Covariance Type: nonrobust
===============================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 38.6517 9.456 4.087 0.000 19.826 57.478
Region[T.E] -15.4278 9.727 -1.586 0.117 -34.793 3.938
Region[T.N] -10.0170 9.260 -1.082 0.283 -28.453 8.419
Region[T.S] -4.5483 7.279 -0.625 0.534 -19.039 9.943
Region[T.W] -10.0913 7.196 -1.402 0.165 -24.418 4.235
Literacy -0.1858 0.210 -0.886 0.378 -0.603 0.232
Wealth 0.4515 0.103 4.390 0.000 0.247 0.656
==============================================================================
Omnibus: 3.049 Durbin-Watson: 1.785
Prob(Omnibus): 0.218 Jarque-Bera (JB): 2.694
Skew: -0.340 Prob(JB): 0.260
Kurtosis: 2.454 Cond. No. 371.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.The Region column holds text strings, so patsy treated it as a categorical variable: it created one indicator (0/1) column per region, Region[T.E], Region[T.N], …, and dropped one level (C, the first alphabetically) as the reference category that the intercept represents. Each Region coefficient is the difference in mean Lottery between that region and region C, holding Literacy and Wealth fixed.
Wrapping a variable in C() forces it to be treated as categorical. This matters when a category is coded as a number – an institution rank of 1 to 4, a zip code, a month – and would otherwise be used as a continuous predictor.
Combining C() with - 1 removes the intercept, and then every region gets its own indicator column: the coefficients become the regions’ adjusted means instead of differences from a reference level.
res = smf.ols(formula='Lottery ~ Literacy + Wealth + C(Region) - 1', data=guerry).fit()
print(res.summary()) OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.338
Model: OLS Adj. R-squared: 0.287
Method: Least Squares F-statistic: 6.636
Date: Thu, 15 Oct 2026 Prob (F-statistic): 1.07e-05
Time: 21:51:33 Log-Likelihood: -375.30
No. Observations: 85 AIC: 764.6
Df Residuals: 78 BIC: 781.7
Df Model: 6
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
C(Region)[C] 38.6517 9.456 4.087 0.000 19.826 57.478
C(Region)[E] 23.2239 14.931 1.555 0.124 -6.501 52.949
C(Region)[N] 28.6347 13.127 2.181 0.032 2.501 54.769
C(Region)[S] 34.1034 10.370 3.289 0.002 13.459 54.748
C(Region)[W] 28.5604 10.018 2.851 0.006 8.616 48.505
Literacy -0.1858 0.210 -0.886 0.378 -0.603 0.232
Wealth 0.4515 0.103 4.390 0.000 0.247 0.656
==============================================================================
Omnibus: 3.049 Durbin-Watson: 1.785
Prob(Omnibus): 0.218 Jarque-Bera (JB): 2.694
Skew: -0.340 Prob(JB): 0.260
Kurtosis: 2.454 Cond. No. 653.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.Same fit, same predictions, different parameterization: compare the C(Region)[C] coefficient here with the Intercept in the previous table.
In the logistic regression lecture the ordinal rank (1–4) was used as a number, which assumes each step in rank shifts the log-odds by the same amount. The formula 'admit ~ gre + gpa + C(rank)' with smf.logit drops that assumption at the cost of two more coefficients – try it, and compare the three rank coefficients with the single slope the lecture reports.
We can also apply vectorized functions to the variables in our model, directly in the formula. Here we apply the natural log to the Literacy column and use that single transformed column to predict Lottery.
res = smf.ols(formula='Lottery ~ np.log(Literacy)', data=guerry).fit()
print(res.summary()) OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.161
Model: OLS Adj. R-squared: 0.151
Method: Least Squares F-statistic: 15.89
Date: Thu, 15 Oct 2026 Prob (F-statistic): 0.000144
Time: 21:51:33 Log-Likelihood: -385.38
No. Observations: 85 AIC: 774.8
Df Residuals: 83 BIC: 779.7
Df Model: 1
Covariance Type: nonrobust
====================================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------------
Intercept 115.6091 18.374 6.292 0.000 79.064 152.155
np.log(Literacy) -20.3940 5.116 -3.986 0.000 -30.570 -10.218
==============================================================================
Omnibus: 8.907 Durbin-Watson: 2.019
Prob(Omnibus): 0.012 Jarque-Bera (JB): 3.299
Skew: 0.108 Prob(JB): 0.192
Kurtosis: 2.059 Cond. No. 28.7
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.Any function visible in the calling namespace works (np.log, np.sqrt, your own), and I() protects ordinary arithmetic: I(Wealth**2) adds a squared term, whereas a bare Wealth**2 would be read as formula syntax.
| Syntax | Meaning |
|---|---|
y ~ a + b |
regress y on a and b, with an intercept |
y ~ a + b - 1 |
the same, without an intercept |
y ~ C(a) |
treat a as categorical (one indicator column per level, one level as reference) |
y ~ np.log(a) |
apply a function to a column |
y ~ a + I(a**2) |
add a transformed column that uses arithmetic operators |
The same formulas work with smf.logit for logistic regression and with the other models in statsmodels.formula.api. For the full language – interactions (a:b, a*b), nested terms, and custom coding – see the patsy documentation.