NAMpy: Generalized Additive Models and Interpretable Deep Learning in Python

I’ve always found generalized additive models one of the most elegant ideas in statistical modeling. They sit in an interesting place between linear regression and the more flexible machine-learning models. You can learn fairly complicated nonlinear relationships from data, and you can still look at the individual functions and see what the model is doing.

At the same time, I’ve often felt that the ecosystem around these models is fragmented.

If you work with statistical GAMs, you probably know R’s mgcv. It is an incredibly powerful piece of software, with years of work behind its spline construction, smoothness selection, inference, and diagnostics. If you work primarily in Python, particularly with deep learning, there isn’t quite the same unified experience.

On the other side there has been a lot of interesting work on neural additive models. Instead of using splines to represent feature effects, these models learn the functions with neural networks. Papers on NAMs, Neural Basis Models, NodeGAM, and several other architectures have explored different ways of keeping a flexible model interpretable.

What struck me was that these approaches are closely related mathematically, while their implementations tend to live in separate places. Each comes with its own way of preparing data, training, predicting, and drawing the fitted functions.

That was one of the main reasons I started NAMpy.

NAMpy is an open-source Python library that puts classical statistical GAMs and neural additive models in one framework. It includes a Python-native implementation of much of the statistical machinery associated with mgcv, together with a collection of neural architectures, and methods for distributional regression, quantiles, shape constraints, and more.

The part I find most interesting is not really the number of methods. It is that the methods can share a common architecture, which makes it much easier to try a new idea without rewriting the surrounding machinery.

Why additive models?

The basic idea is surprisingly simple.

Suppose we are trying to predict a response \(Y\) from features \(X_1,\ldots,X_p\). In ordinary linear regression we write

Each feature contributes linearly. That is easy to understand, and often too rigid. The effect of temperature on electricity use might be U-shaped. An income effect might flatten. Age might rise, level off, and eventually fall.

A generalized additive model keeps the additive structure and replaces the linear terms with functions:

Here \(g\) is a link function, and each \(f_j\) is learned from the data.

It is a small change in notation, and a very useful one. We stop having to write down the exact shape of a feature effect, and we do not give up the ability to look at the pieces separately.

There are real limitations. An additive model with only main effects cannot represent an arbitrary interaction. A fitted feature effect is not automatically a causal effect. As a starting point for interpretable nonlinear modeling, though, I still think this is a remarkably good framework.

Why bring mgcv to Python?

One thing I wanted from NAMpy was to make serious statistical GAM modeling available directly in Python.

There are already Python libraries that fit splines and generalized additive models. mgcv goes considerably deeper than fitting a smooth curve. A large part of the work is deciding how smooth a function should be, and then being honest about how precisely we know that function.

Typically we represent a function through basis functions,

and estimate the coefficients by optimizing a penalized likelihood,

The penalty stops the fitted function from becoming unnecessarily wiggly. The smoothing parameter \(\lambda\) decides how strongly we penalize that complexity. Too little regularization and the curve follows noise. Too much and it misses structure that is actually there.

Which basis we use matters less often than people expect, once the penalty is doing its job. A cubic regression spline is what I reach for first on a single covariate. The basis is cheap, the penalty is the integrated squared second derivative, and \(k\) is the dimension I am willing to allow before the penalty takes over. A thin-plate spline is what I use when I do not want to place knots by hand. It is built from the covariates themselves, and the same construction extends to several covariates, with a truncated eigenbasis so that the thing remains computable. A P-spline feels different in a way I like. The basis is local, usually B-splines, and the penalty is a discrete difference penalty on neighbouring coefficients. That makes it easy to change what “smooth” means, and it is the route into the shape-constrained bases later.

On the left of the figure below I fit those three bases, with the same budget of basis functions, to one noisy curve. The estimates sit almost on top of each other. The interesting differences show up when the covariates stop being a single column of numbers. A tensor product, te, multiplies one-dimensional marginal bases and gives each margin its own smoothing parameter, so a time effect and a temperature effect can be smooth to different degrees. That is a different object from a thin-plate surface, which uses one notion of roughness for the whole plane. If I already have main-effect smooths and I only want the interaction left over, ti is the version that tries not to re-fit those main effects inside the surface. A random effect, re, is the other end of the same machinery: a penalty that shrinks a whole group of coefficients toward a common value. I spent a long time thinking of these as separate model classes. In the penalty formulation they are neighbours.

The right-hand panel is the case where the basis is no longer a minor choice. The true mean rises overall, but it does turn down for a stretch in the middle. An ordinary P-spline follows that dip. A monotone increasing P-spline, mpi, refuses it and waits until the data force the curve up again. That is what I want when the science says the effect cannot decrease. It is the wrong constraint when the dip is real. I like seeing both fits, because the constraint is an assumption I am making, not a property the data have offered me.

Packages like mgcv choose \(\lambda\) by restricted maximum likelihood, along with a substantial amount of numerical linear algebra. REML treats the penalty as a Gaussian prior on the coefficients and integrates those coefficients out. What is left is a marginal likelihood for the smoothing parameters, and an outer Newton method maximises it. The practical consequence is on the right of the next figure. With almost no penalty the cubic regression spline picks up a little extra wiggle that is not in the true mean. With a very heavy penalty the same basis is pulled toward a flatter curve and misses the top and the bottom of the wave. REML, on the left, lands close to the function I actually simulated.

Once \(\lambda\) is chosen, the coefficients still have uncertainty. The object I care about is the penalized information for \(\beta\): the expected curvature of the log-likelihood, plus the penalty,

In the usual penalized least-squares calculation this is the weighted Gram matrix of the basis plus \(S_\lambda\). Its inverse is the Bayesian covariance of the coefficients, conditional on the smoothing parameter that we just selected. Take the row of the basis that evaluates the smooth at a new \(x\), and the quadratic form of that row with this covariance is a pointwise variance for the fitted function. The shaded band on the left is that calculation, drawn as an approximate 95% interval. It gets wider toward the edges, where the data say less. It is a statement about the estimated mean, not about a new observation, and it is pointwise. It is not a band that covers the whole curve at once.

There is a further, smaller uncertainty that this band leaves out. Because \(\lambda\) was estimated, a band that conditions on the selected value is a little too narrow. Propagating that extra uncertainty uses the curvature of the REML criterion in the log smoothing parameters, which is the Fisher information for those parameters, and folds it back into the coefficient covariance. NAMpy can carry that correction when the family supports it. I try not to draw one shaded region and call it the uncertainty. The usual picture conditions on the smoothness we chose.

I wanted these pieces in Python without an R process behind the scenes. Fitting a GAM with automatic smoothness selection looks like this:

from nampy.models import GAMRegressor

model = GAMRegressor(
    formula="y ~ s(x1, bs='cr', k=12) + s(x2, k=12)",
    smoothing_method="reml",
).fit(data)

print(model.summary())
predictions = model.predict(data)

The syntax should look familiar if you have used mgcv.

Something that has mattered a lot during development is numerical validation. Implementing a spline basis is one thing. Reproducing the statistical behavior of a mature package is another. For many of the supported methods, NAMpy tests fitted predictions, coefficients, smoothing parameters, effective degrees of freedom, and covariance estimates against results from R.

There are still differences and unsupported combinations, and I try to document them rather than paper over them. Shape-constrained smooths are one example: their smoothing parameters are selected by GCV, not by the REML outer iteration used above. The goal is not to claim that NAMpy is a drop-in replacement for every corner of mgcv. It is to bring as much of that statistical functionality as I can into a Python-native implementation, and to keep it connected to the rest of the Python machine-learning stack.

And then there are neural additive models

Additive models do not have to use splines.

In a classical GAM the functions \(f_j\) are spline bases with smoothness penalties. We can also give each feature its own neural network,

That is the idea behind neural additive models. Each feature has a small network, and the outputs are summed. The networks can learn nonlinear functions, and the model still has an additive decomposition.

I find the connection interesting because the two approaches are trying to learn the same kind of object with very different control of complexity. In one case I choose a basis and a penalty, and REML picks the smoothing parameter from the marginal likelihood. In the other I choose a network, a learning rate, and whatever regularization is in the optimizer. Both are function estimation. They fail in different ways, which is why I want them in the same place.

NAMpy includes the original NAM idea and several related architectures, including Neural Basis Models, NodeGAM, IGANN, GPNAM, and SIAN. Some use a network per feature. Others share a representation. Some are explicit about interactions. The interface is deliberately similar. A neural additive model can be fitted with:

from nampy.models import NAMRegressor

model = NAMRegressor().fit(X_train, y_train)
predictions = model.predict(X_test)
figures = model.plot_terms(X_test, center=True, reference_X=X_train)

What I like about this is that I can fit a statistical GAM and a neural additive model to the same table and compare the functions, not only a score. Two models with similar predictive error can still disagree about the shape of a relationship. On the simulated example below they mostly do not disagree. Both recover the sine in the first feature and the bowl in the second. The spline sits a little closer to the top of the sine. The neural fit is slightly flatter through the bottom of the quadratic. I would not turn one draw into a verdict about which class of model is better. I would want to see this comparison every time the sample size, the noise, or the correlation between predictors changes.

Interpretability gets more subtle here. A model that returns a separate contribution for each feature has not thereby shown that the contribution is unique or stable. Once an architecture has interactions, or a component that is shared across features, I have to be careful about what I call a feature effect. I have tried to keep those distinctions explicit in NAMpy, rather than treating every decomposition as the same kind of object.

Modeling more than the mean

Most regression models are built to estimate a conditional mean,

A mean is not the whole response. Two households can have the same expected electricity use and very different variability. In a loss distribution, the tail can matter more than the center.

Suppose

and both parameters have additive predictors,

Then a feature can move the level of the response and, separately, its spread. This is the idea behind GAMLSS, generalized additive models for location, scale, and shape.

NAMpy supports distributional regression on the statistical side and on the neural side. For a normal response the formulas can be written separately:

from nampy.models import GAMLSS

model = GAMLSS(
    family="normal",
    formula={
        "mu": "y ~ s(x1, k=12) + x2",
        "sigma": "~ s(x1, k=8)",
    },
).fit(data)

parameters = model.predict_parameters(data)

The sigma formula is on the log scale. predict_parameters returns the mean and the standard deviation themselves. In the figure I simulated a sine wave for the mean and a standard deviation that grows with \(x\). The fit recovers both shapes. The mean is a bit damped at the peaks, which is what a smoothness penalty does when it is not completely sure, and the scale comes out a little low while still rising. The same idea carries over to neural distributional models, where the additive functions are networks rather than splines.

There are other families in the library, and methods for quantiles, expectiles, censored outcomes, and survival. I think of these as the same modeling habit pointed at a different summary of the response. Once I can learn an interpretable function for a conditional mean, I do not see a good reason to stop there.

Constraints that come from the problem

Sometimes I already know something about the shape, even when I do not know the function. A dose effect might have to increase. I may not know whether it increases linearly or rises quickly and then flattens. I do know it should not decrease. That is a constraint on the derivative,

and there are analogous constraints for convexity, or for a function that has to stay positive.

NAMpy includes shape-constrained smooths in the tradition of scam. The constraint is part of the fit, not a cleanup applied to an unconstrained curve afterward. The right-hand panel of the first figure is the picture I keep in mind. Because these bases are parameterized through a transform of the coefficients, smoothness selection currently goes through GCV rather than REML:

from nampy.models import GAMRegressor

model = GAMRegressor(
    formula="y ~ s(x, bs='mpi', k=12)",
    smoothing_method="gcv",
    smoothing_optimizer="bfgs",
).fit(data)

mpi asks for a smooth that is monotone increasing. I like this style of modeling when the constraint really does come from the problem. The most flexible model is not always the one I want. Sometimes I want a function that is flexible enough, and that refuses shapes I am not willing to entertain.

The part I enjoy most

Adding statistical and neural methods has been interesting. The part of NAMpy I am most attached to is the internal split of responsibilities.

It is tempting to implement every paper as its own stack. A network, a training loop, a preprocessor, a prediction method, a plotting helper. After a few models, most of that code is solving the same problems again.

So I started separating them. On the neural side, the architecture that represents the functions is separate from the objective, the response distribution, and the training procedure. Data preparation, prediction, saving, and the interpretation tables are shared. On the statistical side, the formula parser, the smooth constructors, the penalty matrices, the model compiler, the optimizers, and the prediction code are separate components.

That split has been more useful than I expected. A new architecture for feature effects can reuse the training and the plotting. A quantile GAM or an expectile GAM can reuse the basis construction and change the criterion that is optimized. Not every combination is meaningful. Different objectives really do need different assumptions, and the inference that is valid for one of them can be nonsense for another. Having the pieces named makes those mismatches easier to refuse explicitly, instead of discovering them as a confusing error three layers down.

This is the direction I want for NAMpy as a research library. I would like experimenting with a new additive method to mean implementing the part that is actually new.

Where I want to take it

There is a lot left that I want to understand, not only to implement.

I want cleaner comparisons between penalized splines and neural networks as ways of learning an additive function. When does the spline recover the relationship more reliably? How much of that depends on the sample size, the noise, and the correlation among predictors? What changes when an interaction is really there? The small simulation above is a reminder that agreement on one example is cheap.

I am also interested in models that sit between the two, with a structured spline for the parts of the predictor we know how to penalize, and a neural piece for the rest. Distributional regression has the same appeal. A lot of recent work makes the predictive distribution more flexible. I think there is still a lot to be gained by keeping an additive structure while doing that, so that a change in spread can be read off a function rather than only off a score.

And I keep coming back to the uncertainty question from the Fisher information. Regularization, variable selection, and the stability of a learned effect are tied up with the fact that the smoothing parameter, or the optimization path in a neural model, was itself estimated. A band that pretends otherwise is tidier than the experiment we actually ran.

NAMpy already contains a fairly broad set of methods. I do not think of the project as finished. The point of the generic architecture was to make it easier to keep asking these questions, and to keep the numerical checks, the benchmarks, and the comparisons with existing implementations growing as the methods do.

Trying it

NAMpy is open source. You can install it with:

pip install nampy

The documentation goes from introductory GAMs and neural additive models through smoothness selection, distributional regression, shape constraints, interactions, and more specialized statistical methods.

If you know R’s mgcv, much of the statistical interface should feel recognizable. If you are coming from PyTorch or scikit-learn, the neural models should be a familiar kind of object.

The project is on GitHub, and the documentation is at nampy.readthedocs.io.

I would like to hear from people who work with GAMs, distributional regression, interpretable machine learning, or new additive architectures. Suggestions, issues, and contributions are welcome.

More than anything, I hope NAMpy makes it a little easier to move between statistical modeling and machine learning without dropping the useful ideas from either side. I think there is still a lot to learn from keeping the two closer together.