Healthy Algorithms

I haven’t seen a question like this before today. I wonder what the answers have been like.

does this blog still work?

I’ve been appreciating the Python debugger lately, and I want everyone who does data science work in Python to have a chance to appreciate it, too.

When I went looking for a good place to refer colleagues who want to learn this tool, I decided this extensive tutorial might be best: https://realpython.com/python-debugging-pdb/ But if you have something else you think they could use, please let me know.

I thought that there was a Software Carpentry lesson on this as well, but I’m not sure where the definitive source is (perhaps this is close?), and I think Wes McKinney’s book has a very practical section on it, but I can’t confirm that without a trip to my office (it looks like his new book will have it!).

Also I use pdb++ as much as I can now, and this inspired me to read the docs about some cool things it does that I haven’t used yet so maybe I’m going to appreciate debugging in Python even more soon.

A paper I helped with is now in print, Comfort et al, Association Between Subnational Vaccine Coverage, Migration, and Incident Cases of Measles, Mumps, and Rubella in Iraq, 2001–2016.

It is a good chance to test out a new python package for regression modeling that I have been excited about, the BAyesian Model-Building Interface (Bambi).

In the past, it has sometimes been too much work to include random effects in a regression model in Python. The heart of the methods section in this paper, for example, is this “In our linear mixed effects regression model, we set vaccine coverage as the independent variable and disease incidence per 100,000 as the dependent variable. We also included governorate as a random effect term to control for any correlation in incidence within governorates.”

This sort of model is standard enough that there are multiple R and Stata methods to apply it in one line once you have your data all prepped and loaded. But before Bambi, I didn’t know an easy way to do it in Python. I could code up a custom PyMC model, but that seems like more work than it should be.

So I was pleasantly surprised when I went to replicate some of Haley’s findings that the Bambi code to do this is super-simple:

model = bmb.Model("incidence ~ lagged_coverage + (1|gov)", df, dropna=True)
results = model.fit()

If you want to run it yourself, here is an IPython Notebook that includes the data wrangling and shows that code in context: https://gist.github.com/aflaxman/3fad36937c2d082abf99314061e16db1

I have had my eye on a python package called Bambi for a while now, because I often need a regression model that is a little more complicated that sklearn.linear_model.LinearRegression but not complicated enough to make a whole new PyMC model.

Here is a minimal example (adapted from Awesome Open Source):

import bambi as bmb
import numpy as np, pandas as pd, arviz as az

data = pd.DataFrame({
    "y": np.random.normal(size=50),
    "x1": np.random.normal(size=50),
    "x2": np.random.normal(size=50)
})

model = bmb.Model("y ~ x1 + x2", data)
results = model.fit()
az.summary(results)

One cool thing about Bambi is that while it is simpler that writing a whole new PyMC model, it is a lot like writing a PyMC model. For example, if I need to add an informative prior, that is pretty easy:

priors = {'x1': bmb.Prior("Uniform", lower=0, upper=.05)}
model = bmb.Model("y ~ x1 + x2", data,
                 priors=priors)
results = model.fit()
az.summary(results)

And if I need a more complex distribution on that prior, Bambi exposes a “potential” parameter that puts additional terms in the posterior distribution, just like PyMC:

potentials = [
 (('x1', 'x2'),
 lambda x1, x2: bmb.math.switch((x1+x2 <= 0.0), 0, -999)),
]

model = bmb.Model("y ~ x1 + x2", data,
                 potentials=potentials)
results = model.fit()
az.summary(results)

I’m guessing that the syntax will continue to evolve, which is just one more reason Bambi is a python package that I am going to continue to watch.