BART — Bayesian Additive Regression Trees

An Ensemble With a Posterior

Every other ensemble in this section returns a number. BART returns a distribution. It keeps the additive-tree structure of boosting — a sum of many shallow trees — but puts a prior on it: one that keeps each tree small, so no single tree can dominate, and a prior on the leaf values that shrinks them toward zero. Fitting is MCMC rather than greedy search, so the output is a posterior over the whole function, and uncertainty becomes a first-class result instead of an afterthought.

That is the payoff, and it is a real one — measurable rather than rhetorical. On the credit data the posterior pins low- and high-risk clients down tightly while the mid-range carries wide bands: interval width runs from 0.081 at a predicted risk of 0.06 to 0.344 at 0.61, and the correlation between width and closeness to the 0.5 decision boundary is +0.95. BART tells you which predictions to trust, which no point-prediction forest can. The probabilities themselves are close to honest — expected calibration error 0.025 across deciles of predicted risk. Out-of-sample AUC is 0.759 with a mean credible-interval width of 0.188; R's dbarts reaches 0.783 with a tighter 0.158.

Credible is not predictive

One distinction deserves emphasis, because getting it wrong halves your stated uncertainty. The posterior over the mean function is not a predictive interval — to cover an actual observation you must add the observation noise $\sigma$. Measured on the same fit in the R notebook: the credible interval on the mean covers just 50.2% of held-out values, while the predictive interval including $\sigma$ covers 91.7% against its nominal 90%. Both are correct answers — to different questions.

And why they over-cover

Why 96.2% and not 90%? An aggregate coverage number can hide a real failure — intervals too wide at one end of the range and too narrow at the other, the errors cancelling in the average. That is not what is happening here, and checking was worth the trouble. Split the test set into deciles of the prediction and coverage runs 0.884 to 1.000, at or above nominal in nine of ten deciles. The intervals are simply too wide almost everywhere, so the headline figure is an honest summary of them.

The mechanism shows up in the interval widths, which are nearly constant at 2.97 to 3.15 across the whole prediction range. They are dominated by a single estimated noise term, $\sigma = 0.690$, shared by every block group: that alone makes the band about 2.27 wide, against a response spanning roughly 0.2 to 5.0, and the posterior spread of the mean function — the genuinely Bayesian part — is a minor addition on top. So a homoskedastic likelihood lays a constant-width band over a target whose dispersion plainly varies with level, and over-covering is the arithmetic consequence. Coverage comes closest to nominal (0.880) exactly where the response is most spread out. The 96.2% is a modelling choice made visible, not a calibration failure hidden by averaging — and a heteroskedastic BART or a quantile variant is the fix.

And here the two engines part company in a way the headline numbers reverse. Run the same conditional check on dbarts and coverage falls almost monotonically with the predicted value — 0.993 in the cheapest deciles to 0.797 in the dearest, with only six of ten reaching nominal. Its intervals are too wide at the bottom of the range and too narrow at the top, and the two errors cancel to a healthy-looking 91.7%. That is the failure a single number cannot show, and it matters most where it appears: an interval covering 80% of outcomes in the expensive blocks is not delivering what it promises exactly where a valuation error costs the most.

So on the aggregate dbarts looks the better-calibrated implementation — 91.7% against PyMC-BART's 96.2%. Conditionally it is the one with the real problem. PyMC-BART over-covers, but it does so almost uniformly; dbarts lands closer on average because its narrower bands trade over-coverage at one end for under-coverage at the other. Neither notebook would have found that from the headline figure, and it is the same lesson as the RMSE gap below: a BART result is an implementation's result.

Two implementations, one disagreement

The two implementations do not agree, and the honest thing is to say so. (PyMC-BART's PGBART sampler is not seeded reproducibly, so every figure on this page drawn from it — RMSE, AUC, calibration error, interval widths — moves slightly between executions. The movement is in the third decimal and is far smaller than any comparison drawn here, but it is why these should be read as one run rather than as constants.) On the California regression, dbarts reaches RMSE 0.526 — right alongside the random forest's 0.523 — while PyMC-BART reaches 0.667 — 27% higher, the figure the R notebook now computes from the Python result directly rather than from a copied-in constant — with intervals that over-cover at 96.2% where dbarts sits at 91.7%. Three explanations were tested and ruled out: sampling (raising PyMC-BART to 4 chains × 1000 draws moved RMSE from 0.666 to 0.667), ensemble size (dbarts is flat between 50 and 200 trees, 0.532 against 0.531), and the prediction path (PyMC-BART's in-sample RMSE via its own posterior and via its prediction helper agree to a correlation of 0.9996). What remains is that PyMC-BART does not fit this data as closely in the first place — its in-sample RMSE is already ~0.64, worse than dbarts manages out of sample.

What to take from it

So BART is competitive with the best ensembles here and returns a posterior — but that is a statement about the method, not about any one library's number. A BART result is an implementation's result, and this arc only discovered that by running the same analysis twice in two languages. It is the same lesson the CART example reached from the other direction, where R and Python disagreed on whether a tree beats a logistic regression because they pruned differently.

The same method elsewhere

BART appears again on the joint scoreboard in Gaussian Processes & Splines — the Bayesian Kernel View, set beside Gaussian processes, additive splines and the frequentist ensembles on these same two datasets. The implementation gap diagnosed above is what decides its placing there: level with the random forest on credit, but 0.667 against 0.523 on California with pymc-bart, where dbarts reaches 0.526 and closes most of the gap.

Notebooks

Downloads

References

• Chipman, H. A., George, E. I. & McCulloch, R. E. (2010). BART: Bayesian additive regression trees. Annals of Applied Statistics 4(1), 266–298. — the model, its priors, and the backfitting MCMC
• Hill, J. L. (2011). Bayesian nonparametric modeling for causal inference. JCGS 20(1), 217–240. — why the posterior matters when the estimand is a treatment effect
• Dorie, V. (2023). dbarts: Discrete Bayesian Additive Regression Trees Sampler. — the R implementation used here
• Quiroga, M., Garay, P. G. et al. (2022). Bayesian additive regression trees for probabilistic programming. — PyMC-BART and its particle-Gibbs sampler
• Linero, A. R. (2018). Bayesian regression trees for high-dimensional prediction and variable selection. JASA 113(522), 626–636. — on BART's behaviour as dimension grows