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.