This is Bob.
Are neural nets the future of regression?
Andrew was visiting Flatiron last Friday (really last Friday, not six months ago), and I was asking the question that’s been on my mind lately: will neural networks put regression modelers out of work?
Andrew hired me and Matt Hoffman in 2010 to work out how to specify and fit hierarchical regression models with interactions. He wanted to create a system that automatically added interactions, non-linearities, etc., guided by a vaguely conceived “topology of models.” This is a combinatorial nightmare, even with a handful of covariates and non-linearities, and not even considering continuous variation in things like priors.
Black box non-linear function approximations
Fast forward 15 years and regressions from neural networks are ubiquitous. Rather than specifying interactions, non-linearities, etc., we just let a highly overparameterized deep neural net sort it out. This idea of black-box, non-linear function approximation is not new. I first saw it with random forests (the Bayesian analogue of which is Bayesian additive regression trees) and more recently, gradient-boosted decision trees (the go-to method in Kaggle competitions).
Do we have enough data?
The only thing holding us back from using neural networks everywhere is limited data. It’s clear as our data sets get bigger that neural network regression works very well (see, e.g., LLMs, image recognition, and image generation systems, all of which fit largely black-box deep neural network models).
Uncertainty quantification
As we were talking about this, Andrew kept returning to uncertainty quantification. I somehow couldn’t convince him that we can do exactly the same thing as we are currently doing. There’s no fundamental difference between using a neural network and using bespoke hand-tooled covariate combinations—just different functions mapping the covariates to expected values.
Here’s a document explaining the connection
I didn’t have time to explain this to Andrew at the board, so I wrote it up as a document. This goes over how you can take a GLM and swap out the linear component for a neural network and then proceed as usual. It contains an example of a two-hidden layer perceptron model coded in Stan.
- Bob Carpenter. 2025. It’s all neural nets to me. Technical note.
I’m always happy to get feedback if people have comments or suggestions. Keep in mind the purpose here is not a publication, but just explaining how a neural network can be swapped in for the linear function in a generalized linear model.
Interesting stuff, Bob, thanks for posting! I think for me this is maybe more about interpretability than uncertainty quantification. Your Stan code is very helpful for this BTW. So I see how we are getting the same uncertainty quantification via MCMC on your alphas and betas for the 3-layer net. I’m just unsure what to make of those terms! Part of this is I haven’t really worked with these models, and so its not clear to me how the conformal arrangement of output from the ‘relu’ function is really doing. In my mind, this is sort of like fitting a big flexible GP or Fourier Series or something, basically a “black box”. By contrast, the linear model is very easy to interpret (well…comparatively easy depending on interactions and whatnot :)). So the uncertainty quantification is easier to map to something I can appreciate there, or with a bespoke non-linear model that has some scientific or mathematical derivation that makes to me in connection with the analysis. OTOH, we can also consider linear models to simply be extremely stripped down Taylor Series’ approximations to some ‘true’ underlying model. From that perspective, swapping out the linear model for one of these neural nets really is just about familiarity rather than any kind of conceptual difference in the framework. Thanks for sharing your perspective!
Interesting thoughts.
Today I was messing around with an interaction model fit using brms, and to really understand the result I had to plot it using the conditional_effects() function. It wasn’t easy for me to parse the printed coefficient weights or mentally derive their implications–it was much easier for me to just check the marginal model predictions.
When I’ve fit neural networks I’ve found myself doing a similar thing. For example, I’ve checked what probabilities my neural net spits out for two people who vary on a couple theoretically important parameters (e.g. what does my model predict for a black male who’s 24, versus a white man who’s 24, versus another white man who’s 58).
I’m not sure I approach neural networks and linear models too differently when the linear models are complex enough (which for me is not very complex: two-way interactions are usually enough for me abandon individual parameter estimates and start plotting model predictions instead).
Agreed Peter, even with simple linear models it is usually the case we need to work with MCMC draws of various combinations of coefficients, rather than the coefficients themselves, to get at the questions of interest.
The parameters don’t have natural interpretations. In a GLM, a coefficient beta[n] can be interpreted as the partial derivative of the linear predictor with respect to the covariate x[n] evaluated at x[n], with all other coefficients and covariates held fixed. It gets a little trickier with non-linearities, but it’s the same effect on the linear predictor.
Nevertheless, in a linear model and a deep net, I think it’s easiest to understand the uncertainty in terms of its effect on posterior predictive inference. Suppose I observe data x and y, then want to sample the outcome for new covariates x_new, i.e., we want to draw from the posterior predictive distribution y’ ~ p(y’ | x, y, x’). We can sample y’ in two states. First, we sample parameters from the posterior, alpha, beta ~ p(alpha, beta | x, y), then we sample the outcome from the sampling distribution, y'[n] ~ p(y'[n] | alpha(m), beta(m), x'[n]). Uncertainty in alpha, beta goes down with more data x, y, but the uncertainty from the sampling distribution cannot be eliminated. With draws of y’, you can characterize predictive uncertainty. This lets you do posterior predictive checks, cross-validation, etc.
Yes, this was exactly the point I was trying to make. In practice, you’d probably want something on the order of 10 layers deep and at least 20 neurons wide, which leads to overparameterization—I’m going to go answer the question about that in a reply to another comment.
> Uncertainty in alpha, beta goes down with more data x, y, but the uncertainty from the sampling distribution cannot be eliminated.
To ask for some clarification: are you making the point here that with sufficient data, uncertainty quantification is only really concerned about sampling uncertainty rather than parameter uncertainty? So despite the non-identifiability of an overparameterized NN, we can more or less content ourselves that the model captures the mean behavior arbitrarily well and all talk of uncertainty quantification really comes down to effectively modeling the sampling noise?
Thanks Bob, the timing is perfect because I’ve just started diving into this stuff and can grok it way better with some Stan code. Great write up!
I feel like this glosses over some of the differences introduced by the severe overparameterization. Imo both multimodality and extremely wide isosurfaces in the likelihood represent a theoretical challenge to bayesian sampling, not just computational (which to me is just the log-likelihood taking a long time to evaluate). In addition, mini-batch learning is important not just for computation but plays a central algorithmic role in regularization (and usually outperforms explicit L^p regularization in practice). Even with infinite VRAM and GPU cores, full-batch fitting does not give desirable results. Given the mini-batch algorithms, I’m also not sure that that the UQ we get with bayesian neural networks have a particularly well defined relationship to the same network fit with “normal” SGD based methods.
I’ve recently gotten into fitting embedding vectors to model data clusters where I might otherwise have used a bayesian hierarchical model. I have to say that the behavior can be kind of futzy with low data in the group. And I haven’t had a lot of luck with the sort of spatiotemporal tasks where gaussian processes tend to do well. You can make them work there, but there’s kind of a stone soup quality where all the actual magic sauce is in the manual work you do massaging data representations and the neural network provides a tiny bit of smoothing.
You’re right—I’m glossing over a world of differences to try to highlight the underlying similarity in structure.
Overparameterization: these models are usually not identified due to overparameterization. Ages ago I reported on this with Stan for evaluations I was doing on logistic regression speed. This relates to what the ML folks call “double descent,” but you can see it in regression as the Wikipedia says (I stumbled on this doing a grid evaluation of Stan vs. R’s glm function about 10 yeas ago).
Minibatches: This ties into the optimization that is being used, stochastic gradient. You’re right that nobody’s really figured out how to do this robustly with other optimizers. They seem to get trapped in local optima in a way that SGD with momentum like Adam does not. If you just try to use Stan’s built-in L-BFGS to fit, you can get a local optimum, but it’s not going to be a good fit. I don’t think it’s well understood why these optimizers work for neural networks. We have a postdoc, Jeremy Cohen, who specializes in just this problem, see, for example, Gradient Descent on Neural Networks Typically Occurs at the Edge of Stability, which has some really nice plots.
Stability: You can use HMC to fit—that’s what Radford Neal did in his thesis (supervised by Geoff Hinton!). But you’re right that it won’t be easily relateable to the fits you get with SGD. But then I don’t know that two fits we get with SGD have a well-defined relationship, either. As Matt Hoffman joked (I don’t think he made it up), the most important parameter in an ML model is the random seed.
Data representation: That was one of the things Andrew was concerned about for hierarchical models. I think that no matter what one does in data science, one is dogged by the need to deal with data representations.
Network architecture: The other thing I swept under the rug is that there are a bajillion potential neural network architectures, activation functions, tricks to training, etc. It’s all very fussy, despite the papers usually describing everything as a simpel and general black box.
> It’s all very fussy
This has been an issue for me recently. Even simple models with lots of dropout — fit it a few times and get wildly different results.
I expect that if I spent more time working on the models I’d end up with something more stable, but it was at least as frustrating as sampler convergence problems.
So the overparametrized DL models with a suitable optimizer are great for prediction. I’m actually trying them for that, though I try to tip my toes there gradually (with JAX that is kind of doable).
But given that you then don’t have a proper posterior, neither a model but a “system” that blends a model, an optimiser and a bit of the skills of the developer, no interpretability except maybe afterwards with some extra techniques, and uncertainty with calibration only, this seems like a totally separate landscape from traditional bayesian work.
Or is there a middle ground?
I got good results with an MLP, strong L2 and the Lion optimiser after some heavy experimentation. True, L-BFGS didn’t make it. “Full batch”, dropout didn’t work. The likelihood is a distribution, not sure about calibration yet. But quite different from doing the same in Stan with a more parametric approach, and even further away from the initial approach of having an interpretable model that could even be sampled. (I’m after predictions in my project, so no problem with that as long as the modelling process is manageable and somewhat generalizable to new data.)
Typos: equation above section 3.1 should be exp(xn…; first sentence in 4 “liner”–>”linear”.
I think I’m buying what you’re selling, but I have one immediate question (and forgive my nontechnical language…):
Classical uncertainty statements of statistical models are of the form “assuming the relationship between the predictors and the predictand has this particular model structure, the uncertainty in value y for a given set of predictor values xn is such-and-such”. And you can bootstrap to get the uncertainties in the predictions, which are a function of the uncertainties in the model parameters, given the model.
In a simple linear model, outside the range of the predictor data, the uncertainties grow linearly. But if we expand the model from a + bnxn to a + bnxn + cnxn^2 + dnxn^3… znxn^m, the uncertainties get huge very quickly outside the range of the predictors, and bootstrapping can reveal large uncertainties with the range of the predictors if the dimensionality of the model starts to rival the dimensionality of the predictors.
How does the uncertainty of a deep neural net behave outside (or even within) the range of predictor data?
I’m not trying to sell anything. I’ve never published anyhting on neural networks or done anything interesting with them. I do see lots of talks about them and have lots of colleagues working on them.
I wish I knew the answer to your question about extrapolation uncertainty. Someone must have thought about this more deeply than I have, which is not at all. I hadn’t even thought much about the linear case closely enough in terms of extrapolation. Thinking through it, the linear fits wobble around the data, so it makes sense the uncertianty grows outside the data range. Won’t it also grow out from the “center” even within the range of the observed covariates? I guess the intercept and slope can vary, so it’s not all pivoted around a single point.
I think the point is that in high dimensions, everything behaves like extrapolation.
As the saying goes, In high-dimensional space, no one can hear you scream. That’s why we need models, which are sometimes explicit and sometimes implicit in the fitting algorithm.
I believe it depends strongly on your activation function, and with a typical feedforward and the fashionable “ELU” family of activations today, the result is usually close to multilinear in each input continuing along the trajectory from the last seen datapoint.
Getting it to do what you want is usually a matter of encoding your prior knowledge in the data representations; periodic transformations, linearizing transformations, etc
Bob:
I agree with your general point that all sorts of different parameterizations can work for regression models. Neural nets as implemented with regularization seem to have the desired “unfolding flower” property that the effective size of the model increases as more data are included in the fit.
One important application area of linear regression is multilevel regression and poststratification (MRP). I’ve argued that we should more generally be talking about regularized regression and poststratification (RRP), because what’s important is (a) the poststratification, which puts the onus on the model to (b) adjust for many variables and their interactions, which requires (c) regularization to fit the model to data. Multilevel regression is just one particular form of regularized regression.
The other thing I wanted to add, and this is also completely consistent with your above post, is Hal Stern’s principle that the most important aspect of a statistical method is not what it does with the data, but what data it uses. I’ve been frustrated over the years with social scientists who seem to work hard to avoid using all the data: they’ll obsess over distributional issues or concerns about how to model ordered categorical data, but then they leave tons of information on the table because they’re too focused on whatever particular model they’re fitting. Or they’ll design a randomized experiment and not collect pre-treatment variables because they think they don’t need to. Or they’ll performa a regression discontinuity analysis and not adjust for obviously important pre-treatment variables that are sitting right in front of them. Or they’ll get stuck aggregating data from different places or different time periods because they don’t know about hierarchical modeling.
A big advantage of multilevel modeling (for some problems) and deep nets (for lots of problems) is that they automatically include regularization, so they are well suited for including lots of information. As you say, deep nets have the advantage of working really well on big problems. Multilevel modeling has some advantages of its own, for example allowing modeling of measurement error and correlation structure, but as you say it typically requires more user modeling effort, so it will often make sense to use more fully-automatic methods to the extent that they are available.
I should have mentioned that you can do what you’re calling regularized regression and poststratification.
Rather than engineering your own hierarchical model, you’d just code up all the inputs and let the neural network figure out the structure. But then I don’t see how that’d work to let you build something like our SARS-CoV-2 PCR testing model where we’d want to get predictions for new groups, because we’d presumably use a one-hot encoding of group membership.
I’ve also been surprised in working with domain scientists that their first inclination is to clean the data rather than try to model how it got that way.
Bob:
Yes, I consider the “throw the data at a neural network so that you can get predictions of y given X for any X” to be an example of regularized regression. Then poststratification is immediate, as long as your prediction algorithm spits out a matrix of N x S predictions, where N indexes the cases being predicted in the poststratification and S is some number of simulation draws. Indeed, such examples exist in the literature, for example there is a paper on BARP, which is MRP using Bayesian additive regression trees rather than multivariate linear or logistic regression for the “regularized regression” step.
Regarding predictions for new groups: Yeah, that’s why we’re always talking about developing a method that combines machine learning with multilevel modeling. It’s also possible to use bootstrapping to hack a prediction for new groups. For example, if you have data from 50 counties and you want to make a prediction for county 51, you’d include a bunch of county-level predictors (in machine-learning vocab, they would be “features”) and then 50 county indicators, and then when you want to make the prediction for county 51, you’d use its features and then, for the county indicators, pick one of the 50 at random.
It’s always possible to hack something, but at some point it might be easier to set the epicycles aside and just fit a model, for example when we’re trying to figure out dosing in pharmacology or when we need to account for nonrandom, nonrepresentative sampling.
Nice writeup. My only comment is that K is depth in Section 2, but I believe K becomes something else in Section 3.2.
I am a health economist, not a statistician, so a lot of this goes over my head. I mostly work with medical data, clinical prediction models, and the clinical-statistical intersection where models inform decision-making. However, I have been reluctant to develop models with a high risk of overfitting to sort things out for me because a) of the capability for overfitting, b) similar to point a, that I’m not just noise mining by throwing in another 100 irrelevant “predictors”, and c) because they’re opaque and totally lack any way to check face validity.
Re: point a) in particular, I think it’s worth recognising that in many cases, there are idiosyncrasies in the data that seem to get isolated in highly flexible models. Having more data helps, but you will always get weird edge cases with humans. I guess I’d be more confident when I could see the parameters and understand how things might be misspecified. Maybe this is Luddism?
Ack, that should read “develop models with a high degree of flexibility”
Interesting stuff Bob.
You say that it’s not clear that it Stan can handle it. Did you try?
Thanks putting this together Bob. I’ve been advocating the GLM-ification of our transformer-based neural-networks back at work, so this writeup reassures me a little.
However, while this all works very nicely in theory, I do think you’re glancing over nonnegligible technical hurdles that have actually been hindering adoption in practice. Especially the discussion around uncertainty quantification is overly simplified.
For one, bootstrapping is infeasible for neural networks. Large models that take many many hours (if not days/weeks) to fit cannot be repeatedly refitted 100s/1000s of times on bootstrap samples.
And as for Bayesian deep learning, we’ve already been there. Not just with Radford Neal / David MacKay in the 90s, but the 2015-2019 renaissance of Bayesian deep learning (coupled with the rise of deep learning programming frameworks and probabilistic languages. I mainly remember lots of good work coming out of David Blei’s lab).
but MCMC didn’t seem to scale to such models and we were left with only ELBO-based variational inference. And even those cases, the aggregation of so many parameters (i.e., DNN weights) have often caused poorly-calibrated uncertainty estimates (of the posterior predictive), pulling the rug under the entire reason for Bayesian deep learning.
Eventually, the bandwagon moved forward (maybe the same researchers joined the post-2016 attention/transformer-based architectures for LLMs), and it feels like there’s less research in that area today. Furthermore, it doesn’t seem any of that late 2010s research sticked an got actually implemented and used in practice.
My hunch is that a solution lies in combining deterministic DNNs with probabilistic Bayesian “heads”. Something like a regular neural network up until the last layer, after reducing and compressing the original data into a denser embedding, and then applying Bayesian technics on that last layer for prediction. The problem with this approach, as far as I understand it, will be scaling/parallelizing this approach.
HI Bob:
Nice article as always. Sorry, but I was too shy to introduce myself at ISBA 2024 ;o( Anyways, the big difference between DNN and other machine learning algorithms like ensembles, is that DNNs do not work well with tabular data like age, sex, race, etc. Unfortunately, that is a very big issue.
Also for consideration, at the intersection of “Bayesian” and “transformer-based deep learning neural networks”, are Prior-data Fitted Networks (PFNs) from Prof. Frank Hutter’s group at the University of Freiburg:
* YouTube, AutoML Freiburg Education, “11.05 Prior Fitted Networks (PFNs)”, Jan. 7, 2025: https://www.youtube.com/watch?v=0Pi9ARZjIGg
* Nature, paper on TabPFN, “Accurate predictions on small data with a tabular foundation model”, Jan. 8, 2025: https://www.nature.com/articles/s41586-024-08328-6
* Perplexity page summarizing the group: https://www.perplexity.ai/page/hutter-group-bayesian-transfor-VX7TP7XiQxqjdkkp1SX5pQ
I would agree that Generalized Linear Models (GLMs) are overused for a variety of reasons and that we should be thinking about nonlinear relationships between predictors and outcomes beyond that induced by a nonlinear inverse link function. But my problem with neural networks and similar schemes that allow for almost any continuous (or perhaps differentiable) relationship is that they put positive probability on stuff that I would want to put negligible probability on (and GLMs rule out). It is the analogue of trying to use the weakest possible prior distribution on a parameter, instead of the prior distribution that favors what you actually believe and tolerates a reasonable amount of departure from that.
For example, I cannot think of anything off the top of my head in the social sciences where the outcome depends on cos(x), but a neural network will happily entertain functional forms that are even crazier than cosine. If you draw from the prior predictive distribution of a neural network, you get a massively overdispersed distribution, and it is not easy to choose prior distributions on the weights and biases to weed out the contributions that are inconsistent with any causal model I could imagine.
Proponents of neural networks are usually thinking of using them in contexts where the number of data points is large and the signal-to-noise ratio is high. They usually do not bother to draw from a prior predictive distribution because they presume that the data will induce a reasonable posterior predictive distribution for any non-dogmatic prior distribution. And they usually do not even bother to draw from a posterior distribution because it is difficult to do and because they presume the predictions will not be that different regardless of whatever local optimum they land on (or optima if they do it multiple times or somewhere within shouting distance of some local optima if they employ early stopping). I think all that is questionable even in the intended applications, dubious for a lot of science, and just wrong for most of social science.
If, for example, I were predicting vote turnout with age (and other predictors), something like a Restricted Cubic Spline seems like a reasonable thing because it would allow one or two smooth turns in the function but not thousands.
The argument of this post seems to be: GLMs can be viewed as regularized function approximators; deep neural networks can also be viewed as regularized function approximators. Since deep networks perform much better, why not replace the former with the latter?
But this reasoning overlooks a crucial point (as pointed out by other commenters): deep networks are difficult to interpret—at best. Yes, you can check predictions and swap covariates to observe changes in outputs, but that does not mean you understand the model. In science, interpretability is not just a “nice to have”—it is a fundamental requirement. A scientific model should have parameters that are directly linked to a well-defined theory, and these parameters should come with meaningful uncertainty estimates. Deep networks cannot provide this.
While deep learning is undeniably powerful for non-tabular data such as images, speech, and text, its value in scientific modeling remains limited. From a scientific perspective, deep networks are often uninformative, as they lack the structure and transparency needed to generate theoretical insights. Their strength lies in prediction, not in understanding.
Exactly! My take is that tabular data (unlike image etc) is designed by a hypothesis generating process (our brain’s superior neural processing or what we call a scientist’s brain). No AI model is yet to reach that level. Out of thousands of people only a few combine brain+curated knowledge (many many years of training and honed intuition) to ask the right hypotheses before collecting/using tabular data. Scientists discard (almost everything except the problem at hand) and then bootstrap on the tiny bit of relevant data, NNs average and love bloat.
“The only thing holding us back from using neural networks everywhere is limited data. It’s clear as our data sets get bigger that neural network regression works very well (see, e.g., LLMs, image recognition, and image generation systems, all of which fit largely black-box deep neural network models).”
Also the EU AI Act’s explainability requirement, under which black-boxing is banned in judicial contexts.