Using Patsy in your library¶
Our goal is to make Patsy the de facto standard for describing models in Python, regardless of the underlying package in use – just as formulas are the standard interface to all R packages. Therefore we’ve tried to make it as easy as possible for you to build Patsy support into your libraries.
Patsy is a good houseguest:
Pure Python, no compilation necessary.
Exhaustive tests (>98% statement coverage at time of writing) and documentation (you’re looking at it).
No dependencies besides numpy.
Tested and supported on every version of Python since 2.5. (And 2.4 probably still works too if you really want it, it’s just become too hard to keep a working 2.4 environment on the test server.)
So you can be pretty confident that adding a dependency on Patsy won’t create much hassle for your users.
And, of course, the fundamental design is very conservative – the formula mini-language in S was first described in Chambers and Hastie (1992), more than two decades ago. It’s still in heavy use today in R, which is one of the most popular environments for statistical programming. Many of your users may already be familiar with it. So we can be pretty certain that it will hold up to real-world usage.
Using the high-level interface¶
If you have a function whose signature currently looks like this:
def mymodel2(X, y, ...):
...
or this:
def mymodel1(X, ...):
...
then adding Patsy support is extremely easy (though of course like any other API change, you may have to deprecate the old interface, or provide two interfaces in parallel, depending on your situation). Just write something like:
def mymodel2_patsy(formula_like, data={}, ...):
y, X = patsy.dmatrices(formula_like, data, 1)
...
or:
def mymodel1_patsy(formula_like, data={}, ...):
X = patsy.dmatrix(formula_like, data, 1)
...
(See dmatrices() and dmatrix() for details.) This won’t
force your users to switch to formulas immediately; they can replace
code that looks like this:
X, y = build_matrices_laboriously()
result = mymodel2(X, y, ...)
other_result = mymodel1(X, ...)
with code like this:
X, y = build_matrices_laboriously()
result = mymodel2((y, X), data=None, ...)
other_result = mymodel1(X, data=None, ...)
Of course in the long run they might want to throw away that
build_matrices_laboriously() function and start using formulas,
but they aren’t forced to just to start using your new interface.
Working with metadata¶
Once you’ve started using Patsy to handle formulas, you’ll probably
want to take advantage of the metadata that Patsy provides, so that
you can display regression coefficients by name and so forth. Design
matrices processed by Patsy always have a .design_info
attribute which contains lots of information about the design: see
DesignInfo for details.
Predictions¶
Another nice feature is making predictions on new data. But this requires that we can take in new data, and transform it to create a new X matrix. Or if we want to compute the likelihood of our model on new data, we need both new X and y matrices.
This is also easily done with Patsy – first fetch the relevant
DesignInfo objects by doing input_data.design_info, and
then pass them to build_design_matrices() along with the new
data.
Example¶
Here’s a simplified class for doing ordinary least-squares regression, demonstrating the above techniques:
Warning
This code has not been validated for numerical correctness.
import numpy as np
from patsy import dmatrices, build_design_matrices
class LM(object):
"""An example ordinary least squares linear model class, analogous to R's
lm() function. Don't use this in real life, it isn't properly tested."""
def __init__(self, formula_like, data={}):
y, x = dmatrices(formula_like, data, 1)
self.nobs = x.shape[0]
self.betas, self.rss, _, _ = np.linalg.lstsq(x, y)
self._y_design_info = y.design_info
self._x_design_info = x.design_info
def __repr__(self):
summary = ("Ordinary least-squares regression\n"
" Model: %s ~ %s\n"
" Regression (beta) coefficients:\n"
% (self._y_design_info.describe(),
self._x_design_info.describe()))
for name, value in zip(self._x_design_info.column_names, self.betas):
summary += " %s: %0.3g\n" % (name, value[0])
return summary
def predict(self, new_data):
(new_x,) = build_design_matrices([self._x_design_info],
new_data)
return np.dot(new_x, self.betas)
def loglik(self, new_data):
(new_y, new_x) = build_design_matrices([self._y_design_info,
self._x_design_info],
new_data)
new_pred = np.dot(new_x, self.betas)
sigma2 = self.rss / self.nobs
# It'd be more elegant to use scipy.stats.norm.logpdf here, but adding
# a dependency on scipy makes the docs build more complicated:
Z = -0.5 * np.log(2 * np.pi * sigma2)
return Z + -0.5 * (new_y - new_x) ** 2/sigma2
And here’s how it can be used:
In [1]: from patsy import demo_data
In [2]: data = demo_data("x", "y", "a")
# Old and boring approach (but it still works):
In [3]: X = np.column_stack(([1] * len(data["y"]), data["x"]))
In [4]: LM((data["y"], X))
Out[4]:
Ordinary least-squares regression
Model: y0 ~ x0 + x1
Regression (beta) coefficients:
x0: 0.677
x1: -0.217
# Fancy new way:
In [5]: m = LM("y ~ x", data)
In [6]: m
Out[6]:
Ordinary least-squares regression
Model: y ~ 1 + x
Regression (beta) coefficients:
Intercept: 0.677
x: -0.217
In [7]: m.predict({"x": [10, 20, 30]})
Out[7]:
array([[-1.48944498],
[-3.65620297],
[-5.82296096]])
In [8]: m.loglik(data)
Out[8]:
array([[ -0.28884193, -1.46289596],
[ -2.64235743, -0.8254485 ],
[ -2.44930737, -2.36666465],
[ -0.90233651, -6.24317017],
[ -1.58762894, -5.56817766],
[ -0.65148056, -10.80114045]])
In [9]: m.loglik({"x": [10, 20, 30], "y": [-1, -2, -3]})
Out[9]:
array([[ -7.39939265, -215.51261221],
[ -16.29311998, -861.19721649],
[ -28.74433824, -1937.33822362]])
# Your users get support for categorical predictors for free:
In [10]: LM("y ~ a", data)
Out[10]:
Ordinary least-squares regression
Model: y ~ 1 + a
Regression (beta) coefficients:
Intercept: 0.33
a[T.a2]: 0.241
# And variable transformations too:
In [11]: LM("y ~ np.log(x ** 2)", data)
Out[11]:
Ordinary least-squares regression
Model: y ~ 1 + np.log(x ** 2)
Regression (beta) coefficients:
Intercept: 0.399
np.log(x ** 2): 0.148
Other cool tricks¶
If you want to compute ANOVAs, then check out
DesignInfo.term_name_slices, DesignInfo.slice().
If you support linear hypothesis tests or otherwise allow your users
to specify linear constraints on model parameters, consider taking
advantage of DesignInfo.linear_constraint().
Extending the formula syntax¶
The above documentation assumes that you have a relatively simple model that can be described by one or two matrices (plus whatever other arguments you take). This covers many of the most popular models, but it’s definitely not sufficient for every model out there.
Internally, Patsy is designed to be very flexible – for example,
it’s quite straightforward to add custom operators to the formula
parser, or otherwise extend the formula evaluation machinery. (Heck,
it only took an hour or two to repurpose it for a totally different
purpose, parsing linear constraints.) But extending Patsy in a
more fundamental way then this will require just a wee bit more complicated
API than just calling dmatrices(), and for this initial release,
we’ve been busy enough getting the basics working that we haven’t yet
taken the time to pin down a public extension API we can support.
So, if you want something fancier – please give us a nudge, it’s entirely likely we can work something out.