Frequently Asked Questions

What is BART?

BART stands for “Bayesian Additive Regression Trees”. It’s a Bayesian machine learning method used to model the relationship between predictors and outcomes. It’s similar to GBM/xgboost, extraTrees, and random forests in that the fit is built out of many small trees rather than one big one. BART sums their contributions, as boosting does, where random forests and extraTrees average theirs. This lets it approximate complex and nonlinear functions without you needing to specify their form.

How/When should I use BART?

Use BART exactly the same way you would use any regression method or any machine learning method.

How can a machine learning method be Bayesian?

Each parameter of the model, including which variables are used to split the trees, how deep each tree is, the predicted value in the leaf, the residual variance, etc., is assigned a prior, and through running the model on the data, each model output gains a posterior. In traditional machine learning models, these parameters are given an individual value, or the values are tuned across a grid. The posterior enables Bayesian inference on BART’s predictions.

How can a Bayesian method be machine learning?

Instead of specifying a specific functional form as you would in a Bayesian generalized linear model (GLM), you specify priors governing a tree structure (e.g., the depth of trees, the variables to split on), and that structure is used to model the regression function. It’s still Bayesian because all components of the model have a prior and posterior, but rather than putting priors on GLM coefficients, you put priors on aspects of the machine learning model.

How can BART be used for inference?

Because each prediction from the model has a posterior, any quantity derived from the predictions, like average marginal effects, has one too. These can be used as part of a modern regression workflow through tools like marginaleffects (Arel-Bundock et al. 2024), which produce interpretable estimates from regression models in model-agnostic ways (Rohrer and Arel-Bundock 2026). These tools make the black-box nature of machine learning models no longer problematic for arriving at interpretable model summaries, and the posteriors of these quantities allow for Bayesian inference on them.

How do I get estimates out of the model?

predict() gives you the fitted values or the full posterior of them. estimate_effect() turns those into a treatment effect: the ATE, the ATT, or one effect per unit, with a credible interval and on the response scale, which is the scale on which averaging unit-level contrasts gives the marginal effect. partial_dependence() shows how the prediction moves with one or two predictors, and plot() on a fit draws it. variable_importance() reports how often each predictor is split on, which is a rough guide to what the model is using but not a test of anything.

For the quantities those do not cover, marginaleffects::avg_comparisons() and its relatives work on a fit directly, including arbitrary contrasts between covariate values, hypotheses comparing one estimate to another, and slopes. Because every draw goes through the same machinery either way, the result carries a posterior rather than a point estimate with a delta-method standard error. See vignette("effects") for the worked versions and ?estimate_effect for the estimands.

What if I have missing data?

Missing values in the predictors are fine and nothing is dropped for them, though predict() accepts a missing value only in a column that had one at fitting time, since only those columns’ rules carry an answer for it. Each splitting rule carries its own answer for what to do with a missing value, drawn from the prior alongside the variable and the cutpoint, so the missingness is part of the model rather than something to impute first. This is missingness incorporated in attributes, and it means a covariate can be informative through whether it is observed as well as through its value. Missingness in the outcome is different: those rows are dropped, with a warning saying how many.

Can I supply my own likelihood?

Yes, through custom_family(), which is the feature the package is named for. You give it a function of the response and the additive predictors, function(y, eta), returning one log density per observation, and the sampler does the rest; you do not have to derive a conjugate update or write any C++. The trade-off is that a log density says nothing about how to draw from the distribution it describes, so a custom family has no posterior predictive draws, and everything built on them, including pp_check() and posterior_predict(), is unavailable. Fitting, prediction on the link scale, and loo() all work, though predict(type = "response") hands back the additive predictors rather than a fitted mean, since the package cannot know where the mean of a supplied density sits.

Why is my model slow, and what should I change first?

The settings below matter most, in the order worth trying them. Cost is close to linear in num_trees, and a second forest can usually be smaller than the first, so setting num_trees = c(mean = 50, log_sd = 10) on a location-scale model is much cheaper than 50 for both. Setting gate = "hard" is several times faster than the soft default and costs nothing if the relationship is well approximated by a step function. The bandwidth_every argument controls how often the soft-rule bandwidth is resampled, and raising it keeps soft rules while recovering some of that speed; it is not free, though, since the bandwidth update lets the rules sharpen toward a step, so it costs some mixing and some accuracy on a mean function with jumps. Setting update_bandwidth = FALSE stops resampling it altogether, which is faster again and can actually be more accurate on smooth functions but much worse on nonsmooth ones. Beyond those, run the chains in parallel with a future plan, which costs nothing in draws.

What if I’m a frequentist?

You don’t have to use the full posterior for inference; you can use the predictions from a BART model (e.g., the posterior mean for each prediction) in a frequentist estimator like double/debiased machine learning (DML) (Ahrens et al., n.d.; Chernozhukov et al. 2018) or targeted maximum likelihood estimation (TMLE) (Schuler and Rose 2017; Laan and Rubin 2006). These methods only require predictions from models and allow you to perform valid frequentist inference on specific estimands. That said, some also treat the resulting credible intervals as confidence intervals, but the frequentist operating characteristics of these intervals are not guaranteed.

Is BART good?

Yes! BART or a related method won the American Causal Inference Conference Data Competition in both of the years whose results have been written up (Dorie et al. 2019; Hahn et al. 2019). This is a competition to see which method can estimate treatment effects as accurately as possible (including with accurate inference) across a wide variety of data-generating processes. For general prediction, BART has been shown to do as well or better than popular methods like GBM and random forests (Chipman et al. 2010).

How can BART be used for causal inference?

The only step BART can be used for in causal inference is causal effect estimation, i.e., estimating causal effects given causal assumptions on the data. It works like any other regression method as part of a causal effect estimator. The award-winning use of BART often involves BART-based g-computation (Hill 2011; Snowden et al. 2011; Dorie et al. 2019; Carnegie 2019), which involves modeling the relationship between the outcome, treatment, and confounders, and using that model to predict the counterfactual outcomes under each treatment for each unit. See vignette("causal") for a walkthrough with bartisan. BART can also be used to estimate propensity scores (Hill et al. 2011), which, together with a BART outcome model, can be used in doubly robust estimators like DML and TMLE. BART can also be used in instrumental variables analysis (McCulloch et al., n.d.) and regression discontinuity (Alcantara et al. 2024).

What are Bayesian Causal Forests?

Bayesian Causal Forests (BCF) are a modification of BART used to estimate treatment effects (Hahn et al. 2020). In addition to flexibly modeling the relationship between the outcome and the covariates, it also flexibly models the relationship between the magnitude of the treatment effect and the covariates. In this way, BCF is a varying-coefficient model (Deshpande et al. 2026). Compared to traditional BART, BCF tends to have better calibrated intervals and regularization of heterogeneous treatment effects. BCF tends to outperform traditional BART in most causal inference contexts. It also tends to outperform generalized random forests (GRF), which serve a similar function. bcf() fits it; vignette("varying") covers the varying-coefficient model it is a case of and what bcf() sets on top of it.

How do I get a treatment effect?

estimate_effect() reports it, with the two average potential outcomes printed beneath. A fit from bcf() carries its own treatment, so nothing else has to be named:

fit <- bcf(y ~ x1 + x2, treat = ~ z, data = d)

estimate_effect(fit)

The estimand is an argument: setting estimand = "ATT" averages over the treated instead of over everyone, setting estimand = "CATE" returns one effect per unit, supplying by = ~ g gives subgroup effects, and setting comparison = "ratio" reports a ratio rather than a difference. On a fit from bartisan() rather than bcf() it works the same way once the treat argument names the column. plot() draws whatever was asked for, and diagnose() on the result says whether the estimand itself has mixed, which is not implied by the fit’s own diagnostics.

See vignette("causal") for a worked version, including the assumptions that turn the estimate into a causal effect, which no model supplies.

What papers should I read to better understand BART?

BART has a growing literature spread across multiple fields. We recommend reading the original BART paper by Chipman et al. (2010), the paper introducing BART for causal inference by Hill (2011), an accessible paper on BART for social scientists Green and Kern (2012), and the Annual Reviews paper on BART by Hill et al. (2020).

What does the name “bartisan” mean?

bartisan is a portmanteau of BART and “artisan”. bartisan allows you to supply your own likelihood to fit the BART model of your choice without needing to derive or program the whole machinery to do so. In this way, you can build an artisanal model for your specific scenario. In addition, the package offers a high degree of flexibility and customization for its built-in models, including hard and smooth gates for the trees, random effects, varying coefficient models, many model families, and control of sparsity.

How does bartisan differ from other BART implementations in R?

There are many R packages that implement BART, each of which has its own strengths. bartisan aims to be highly general but with intelligent defaults that make it easy to simply plug in in place of a GLM or generalized additive model (GAM) as part of a data analysis workflow. bartisan offers many more model types and families than other packages. Often, R packages will be built to implement a few specific features of BART; for example, SoftBart implements soft decision trees, VCBART implements the varying coefficient model, dbarts implements Normal and probit BART regression with optional random intercepts. vignette("implementation") has a feature-by-feature comparison against the most popular BART packages, along with the relative cost of each family. bartisan was designed to incorporate as many recent advancements in BART theory and practice as possible, including models for different outcome types (e.g., multinomial, ordinal, \(\left[ 0, 1 \right]\)-bounded, zero-inflated), soft trees, sparsity priors, varying coefficients, random effects, and more, all in a single package¹. There are ways in which bartisan is inferior to these more specialized packages, but as a general-purpose tool, we hope you’ll find it effective.

What model family should I use?

This of course depends on the outcome and the assumptions you are willing to make about it. In general, the defaults are a good starting point, though for numeric outcomes, bounding (e.g., at 0 or 1, inclusive or exclusive, or integer or continuous) might suggest different families. Check out vignette("families") and ?bartisan-families for more details and help with this decision.

How should I select the values of the hyperparameters?

One of the great advantages of BART is that performance tends to be strong using the default hyperparameters, and tuning them with cross-validation tends to buy very little for the increase in computation. We recommend changing the default parameters only if you understand very well what each one does.

That said, there are a few hyperparameters that might be considered. The gate argument controls whether hard or soft decision trees are used; the default, soft decision trees, tend to perform better when modeling smooth relationships, but they cost substantially more computation, because every observation reaches every leaf with some weight rather than taking one side of each split. Setting gate = "hard" is the fast option when the relationship is expected to be a step function anyway. The sparsity argument controls whether a variable-selection prior is placed on the covariates chosen for splitting; it is off by default, so every covariate keeps a share of the splits. Setting sparsity = TRUE is worth considering when there are many predictors, when only a few of them are expected to matter, or when the goal is to identify the subset of predictors the model needs, but not when an effect is being estimated, since a predictor the prior drops has an effect of exactly zero in that draw. The split_prior argument can also be used instead of sparsity to manually decide which variables should be split on more often, though in general the Bayesian updates can discover this from the data anyway. See ?bartisan_control for a full list.

How should I choose how many chains and burn-in and posterior draws?

For all of these parameters, more is always better but requires more computation and more memory usage. It can be a good idea to build your model on synthetic data with the default number of burn-in and posterior draws, and then to use a much higher number for the final analysis on the real data. More chains come almost for free in wall-clock terms if you have parallel processing, as each chain can be run in parallel, whereas burn-in and posterior draws must be run sequentially. They are not free diagnostically, though, and it is worth knowing which way that cuts. Effective sample size depends on the total number of draws you keep and not on how they are split up, so four chains of 10,000 and eight of 5,000 give you about the same precision. But R-hat compares chains against each other, and the value it takes under a perfectly behaved sampler rises with the number of chains being compared: it sits near \(1 + m/S\) for \(m\) chains and effective sample size \(S\). Splitting a fixed budget into more chains therefore raises the bar R-hat has to clear without buying any more information, and a fit that passes at four chains can warn at twelve. Two to four chains is a good default, and when a quantity fails diagnose()’s R-hat threshold on too few effective draws for the comparison to mean anything, the report adds a line saying the number cannot be read yet rather than letting it stand on its own. See vignette("diagnostics") for more information on diagnosing and fixing poor mixing.