The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Getting Started with bartisan

Noah Greifer

2026-10-18

Introduction

bartisan fits Bayesian additive regression trees (BART) using a similar interface to that of glm(). One supplies a formula, a data frame, and a response family, and the model estimates the relationship between the predictors and the outcome without requiring manual specification of nonlinearities and interactions. Unlike most other machine learning models, BART yields a posterior that can be used for valid Bayesian inference of any prediction or function thereof from the model.

In this guide, we will work through a complete analysis at surface level. Other vignettes provide more detail on each step. First, we’ll fit a model and check both that the sampler converged and that the fit describes the data. Next, we’ll see which predictors the forest uses, read the effects off it, and predict for new observations. Finally, we’ll compare two models and note some limitations of the model. The guide assumes a reader comfortable with regression but assumes nothing about machine learning or Bayesian methods.

The main thing to take from this guide is that the default settings are meant to be used (in general). The priors, the number of trees, and the sampler settings are chosen to work across a wide range of problems, and tuning them is rarely where the gains are. Almost everything below is a single function call with no arguments beyond the formula and the data.

library(bartisan)

Parallelization and Progress Bars

Before we fit any models, it’s useful to know that using parallelization can speed up fitting bartisan() models with multiple chains (i.e., distinct runs of the same model). Parallelization is controlled through the future package. To request multi-session evaluation, one can simply call the following:

future::plan("multisession")

The future.apply package must also be installed (installing it also installs future).

Some operations in bartisan can be take some time, and one can request a progress bar using the progressr package. One can request progress bars globally for all functions that produce them using the following:

progressr::handlers(global = TRUE)

To request a progress bar for a single call, wrap the function evaluation in progressr::with_progress(), e.g.,

progressr::with_progress(
  bartisan(y ~ ., data = d)
)

The Data

We’ll use the rhc dataset that comes with bartisan to demonstrate its functionality. rhc includes a random sample of 1500 of the 5735 patients in the SUPPORT study and whether each received right heart catheterization, a monitoring procedure, within a day of arriving in intensive care (Connors et al. 1996). The question the study asked is whether the procedure helps or harms.

data(rhc)

str(rhc)
#> 'data.frame':    1500 obs. of  16 variables:
#>  $ rhc   : int  0 0 0 0 0 1 0 0 0 1 ...
#>  $ death : int  1 0 0 1 0 0 1 1 0 1 ...
#>  $ days  : int  37 235 189 13 202 203 663 15 239 22 ...
#>  $ age   : num  75.3 55 34.4 42.2 41.4 ...
#>  $ sex   : Factor w/ 2 levels "female","male": 1 2 2 1 2 2 2 2 2 2 ...
#>  $ race  : Factor w/ 3 levels "white","black",..: 1 1 1 1 2 1 1 1 1 1 ...
#>  $ edu   : num  9 14 15 16 11 ...
#>  $ aps   : int  48 29 21 55 60 68 26 89 59 105 ...
#>  $ meanbp: num  55 67 66 77 53 47 63 44 50 33 ...
#>  $ resp  : num  26 10 30 40 12 40 22 0 33 44 ...
#>  $ hema  : num  26.3 29 23.8 53 31 ...
#>  $ pafi  : num  157 149 202 171 390 ...
#>  $ paco2 : num  30 45 37 25 31 28 40 39 36 31 ...
#>  $ crea  : num  1.7 1 0.5 1.9 15 ...
#>  $ surv2m: num  0.441 0.339 0.846 0.672 0.777 ...
#>  $ card  : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 2 2 1 1 ...

The outcome comes in two forms. death is whether the patient died during follow-up, and days is how long that took. This vignette uses the binary form; vignette("survival") uses the other.

The remaining variables are the patient’s age, sex, race and years of education, seven physiological measurements taken on the first day, whether cardiovascular disease was among the diagnoses, and surv2m, the study’s own estimate of the patient’s chance of surviving two months. All were recorded before catheterization.

Fitting the Model

To fit the model, we simply call bartisan::bartisan() with the same model formula one would use for glm() or any other standard model-fitting function, additionally supplying the dataset to data and the response family to family. Because BART involves random processes, we must also set a seed using set.seed() to ensure reproducibility (though note the kind does not matter, even with parallelization).

set.seed(2026)

fit <- bartisan(
  death ~ rhc + age + sex + race + edu + aps + meanbp + resp + hema +
    pafi + paco2 + crea + surv2m + card,
  data = rhc, family = binomial()
)

fit
#> Generalized BART
#> 
#> Call:
#> bartisan(formula = death ~ rhc + age + sex + race + edu + aps + 
#>     meanbp + resp + hema + pafi + paco2 + crea + surv2m + card, 
#>     data = rhc, family = binomial())
#> 
#> Family: "binomial" with the "logit" link
#> Observations: 1500
#> Structure: 1 forest of 50 trees, soft decision rules
#> Draws: 800 kept after 200 warmup

Setting family = binomial() says the outcome is binary, exactly as in glm(). The family may be omitted, in which case it is read off the outcome and reported; naming it is clearer and silences the message.

The fit uses the default single chain, which is enough for the estimates and keeps this guide brief. In a real analysis, we would run several chains (e.g., by setting chains = 4) and possibly change the number of burn-in and retained draws before reading any results; vignette("diagnostics") explains how to set these fit controls and how to check that they were enough.

Choosing a family is the one modeling decision that usually matters more than any sampler setting. vignette("families") covers the choice, and vignette("survival") covers censored outcomes such as the days variable in rhc.

The other decision worth knowing about is the shape of the decision rules. By default, a rule is soft, so an observation near a split contributes to both sides of it and the fitted function comes out smooth rather than piecewise constant; the gate argument of bartisan_control() can be used to switch to the hard rules of standard BART, which fit faster but less accurately.

Calling summary() on the output gives a summary of the convergence and variable importance:

summary(fit)
#> Convergence and mixing
#> 
#> Log likelihood: R-hat 1.039, bulk ESS 95, tail ESS 187, over 1 chain
#> 
#> ℹ Use diagnose() (`?bartisan::diagnose`) to examine convergence and mixing
#>   diagnostics.
#> 
#> Variable importance
#> 
#> Predictors ranked by use; 5 of 14 shown.
#>  variable prop_used prop_splits splits
#>    surv2m         1       0.116    8.5
#>       age         1       0.089    6.5
#>     paco2         1       0.073    5.4
#>       rhc         1       0.070    5.1
#>      pafi         1       0.069    5.1
#> 
#> ℹ Use variable_importance() (`?bartisan::variable_importance`) to examine
#>   variable importance.
#> 
#> Further tools
#> 
#> ℹ Use loo() (`?bartisan::loo.bartisan_fit`) to compare this fit with others, or
#>   kfold() (`?bartisan::kfold.bartisan_fit`) if `loo()` reports many Pareto k
#>   values above 0.7.
#> ℹ Use partial_dependence() (`?bartisan::partial_dependence`) and plot()
#>   (`?bartisan::plot.bartisan_fit`) to view the partial dependence of the
#>   predictions on a predictor.

This output is not as meaningful as the full checks, which we describe below.

Checking the Model

In a Bayesian analysis, one must first determine whether the sampler converged before moving forward with interpreting a model’s results. Separately, one should assess whether the model is a good fit to the data. These steps go together because a poorly fitting model will often take longer to converge.

Convergence

diagnose() computes the convergence and mixing statistics, applies the conventional thresholds, and says what to change about whichever of them fall short.

diagnose(fit)
#> Convergence and mixing
#> 
#>                             quantity  rhat rhat_late ess_bulk ess_tail
#>                               loglik 1.039     1.049       95      187
#>                           splits.eta 1.073     1.152       23      134
#>  eta.eta (average over observations) 1.000     1.003      435      723
#>   eta.eta (worst 5% of observations) 1.038     1.069       77      197
#>       bandwidth (average over trees) 1.022     1.037      172      299
#>        bandwidth (worst 5% of trees) 1.064     1.111       50       59
#> 
#> ✖ Only one chain, so R-hat can only compare it with itself; set `chains = 4`
#> ✖ R-hat is above 1.01 for loglik
#> ✖ That R-hat rests on only 95 effective draws, where a single chain averages
#>   1.011 even when its two halves agree
#> ℹ A longer warmup is not the fix: R-hat stays high on the second half of the
#>   draws alone as well
#> ✖ The chains disagree about how many splitting rules the forest has (R-hat
#>   1.07)
#> ✖ The chains disagree about how wide the decision rules are (R-hat 1.06)
#> ✖ Bulk ESS is 77 for eta.eta (worst 5% of observations), below 400
#> ✖ Tail ESS is 187 for loglik, below 400
#> 
#> What to do
#> 
#> • Refit with `chains = 4`. R-hat compares chains against each other, and one
#>   chain can only be compared with itself, so nothing below is reliable until
#>   there are several. With future installed the chains run in parallel.
#> • Raise `num_draws`, which was 800. R-hat is above the threshold for a quantity
#>   that carries too few effective draws for the threshold to mean anything: with
#>   this many chains it would sit about where it does even if the chains agreed
#>   exactly, as the check above reports. A larger effective sample size makes it
#>   readable, and that grows with the total number of draws; using fewer chains
#>   lowers the bar as well, since R-hat's null rises with the number of chains
#>   being compared.
#> • If that does not settle it, reduce `num_trees`, which was 50. A smaller
#>   forest has fewer ways to represent the same fit, so the sampler has less room
#>   to move between them.
#> • Then check the family. A likelihood that fits the data badly can give a
#>   posterior with no single place to be; `bayesplot::pp_check()` is the
#>   diagnostic.
#> • The forest's own size is a different kind of failure from the others and does
#>   not take the same advice. A sum of trees represents one function through many
#>   different partitions, so two chains can agree about every fitted value while
#>   disagreeing about how many rules they used to get there, and the quantities a
#>   fit reports are integrals over that structure. Raising `num_draws` moves this
#>   row slowly and may not clear the threshold at any affordable length. Act on
#>   it when split counts are themselves what gets reported --
#>   `variable_importance()` and `vignette("importance")` -- and not when fitted
#>   values, predictions or effects are.

Because the fit has one chain, the first line of the report says to refit with several, and in a real analysis we would do that before reading the rest. rhat compares variation between chains to variation within them; with one chain it compares the chain’s two halves, which detects drift but cannot detect chains that have settled in different places. rhat_late is the same statistic on the second half of the draws alone, which distinguishes a warmup that ended too early from chains that have each settled somewhere different. ess_bulk and ess_tail are effective sample sizes, and count how many independent draws the correlated ones are worth, in the middle of the distribution and in the tails.

One should read the rows that correspond to the quantities being reported. eta.eta is the fitted function and appears twice, once averaged over the observations and once over the worst 5% of them; the second describes a prediction for one observation. Forests mix slowly on their fitted values, so a figure above 1.01 on the worst-5% row is ordinary rather than alarming; vignette("diagnostics") explains what to do about it and when to worry.

No row here describes an effect, though, and the average row in particular should not be read as if it did. An effect is a contrast, and a contrast can mix badly where the function it is a contrast of mixes well. estimate_effect() produces estimates of the effect of a categorical variable (described more below), and diagnose() can be applied to its output if the effect is of interest.

diagnose(estimate_effect(fit, treat = "rhc"))
#> Convergence and mixing
#> 
#>     quantity  rhat rhat_late ess_bulk ess_tail
#>  Y[1] - Y[0] 1.008     0.998      492      393
#>         Y[0] 1.002     1.002      496      736
#>         Y[1] 1.007     1.007      564      615
#> 
#> ✖ Only one chain, so R-hat can only compare it with itself; set `chains = 4`
#> ✔ R-hat is below 1.01 for every reported quantity
#> ✔ Warmup was long enough, since R-hat is already fine
#> ✔ Bulk ESS is at least 492 for every reported quantity, above 400
#> ✖ Tail ESS is 393 for Y[1] - Y[0], below 400
#> 
#> What to do
#> 
#> • Refit with `chains = 4`. R-hat compares chains against each other, and one
#>   chain can only be compared with itself, so nothing below is reliable until
#>   there are several. With future installed the chains run in parallel.
#> • Raise `num_draws`, which was 800. The chains agree and are stationary, so
#>   they simply have not run long enough. Do not reach for `num_thin`: thinning
#>   discards draws already paid for and lowers the effective sample size per unit
#>   of time.
#> • The tail is the binding constraint, so a posterior mean is already fine and
#>   an interval endpoint is not. Raise `num_draws`, which was 800 if intervals
#>   are what gets reported.

The contrast’s tail effective sample size falls short where the average fitted value’s in the table above did not, which the fit’s own table could not have shown, since the other predictors keep the fitted function moving even where the contrast mixes slowly. The table is in diagnose(.)$table when the values are wanted rather than the report.

Fit

The question of whether the model fits the data well can be investigated using a posterior predictive check, such as those implemented in bayesplot::pp_check(), which is compatible with the output of bartisan(). The default pp_check() compares the distribution of simulated outcomes with the observed distribution. This check only makes sense with continuous outcomes; with binary outcomes, the question worth asking is whether the predicted probabilities mean what they say: among the patients the model gave a 30% chance of dying, did about 30% die? A calibration plot, such as that produced by bayesplot::pp_check(., type = "loo_calibration"), answers it.

bayesplot::pp_check(fit, type = "loo_calibration")

The line should follow the diagonal, and here it does across the whole range. vignette("diagnostics") covers what a departure from the diagonal looks like and what else to check.

Variable Importance

A bartisan() fit produces no model coefficients; to see which covariates played the biggest role in the fit, we can use variable_importance() on its output. This produces statistics about the use of each predictor in the trees.

variable_importance(fit)
#> Variable importance
#> 
#>  variable prop_used prop_splits splits
#>    surv2m     1.000       0.116    8.5
#>       age     1.000       0.089    6.5
#>     paco2     1.000       0.073    5.4
#>       rhc     1.000       0.070    5.1
#>      pafi     1.000       0.069    5.1
#>       aps     0.999       0.082    6.1
#>       edu     0.998       0.062    4.6
#>      hema     0.996       0.065    4.8
#>      card     0.996       0.060    4.4
#>      resp     0.990       0.064    4.7
#>       sex     0.988       0.068    5.0
#>    meanbp     0.988       0.063    4.7
#>      race     0.986       0.059    4.4
#>      crea     0.980       0.060    4.5
#> 
#> ℹ Fitted with `sparsity = FALSE` (the default), so every predictor keeps a
#>   share of the rules and prop_used is near 1 throughout. Refit with `sparsity =
#>   TRUE` to read it as a selection rule.
#> ℹ splits_lower and splits_upper hold the 95% interval, not shown above.

splits is the average number of splitting rules the forest spends on each predictor per draw, prop_splits is its share of all the rules, and prop_used is the proportion of draws in which the predictor received any rule at all. Under the default splitting prior every predictor keeps a share of the rules, so prop_used sits near 1 for all of them, as the note below the table says, and the ranking is in the other two columns.

surv2m takes the most rules, which is unsurprising: it is a prognostic score built to predict survival. age and aps follow. The rest sit close together, at four to five rules per draw, and this table alone cannot say whether that spread means anything.

It’s important to remember that usage is not effect size: a predictor can be split on constantly and still move the prediction very little, and the comparisons in the next section are the better guide to that. Also, when predictors are correlated, the usage distributes among them more or less arbitrarily.

vignette("importance") covers variable importance and selection, including how to tell whether a difference in this table means anything.

Interpreting the Fit

Because a BART model has no table of coefficients to read, a fit is interpreted by asking questions of it: what happens to the prediction when a predictor is changed, and what shape does the prediction trace as that predictor varies. The functions below answer versions of that question, and they differ in what they average over rather than in what they are asking.

A Table of Average Comparisons

marginaleffects::avg_comparisons() is the closest thing to a table of regression slopes. For each predictor in turn, it changes that predictor, leaves the others as they are, predicts every patient twice, and averages the difference over the sample.

marginaleffects::avg_comparisons(fit)
#> 
#>    Term      Contrast  Estimate     2.5 %    97.5 %
#>  age    +1             0.003094  1.53e-03  0.004670
#>  aps    +1             0.001602  1.82e-04  0.003181
#>  card   yes - no       0.041645 -8.68e-03  0.089733
#>  crea   +1             0.005024 -1.80e-02  0.034841
#>  edu    +1            -0.006093 -1.84e-02  0.002854
#>  hema   +1            -0.000788 -3.99e-03  0.002115
#>  meanbp +1             0.000620 -4.11e-04  0.001858
#>  paco2  +1             0.003676  8.53e-04  0.006552
#>  pafi   +1             0.000256  5.17e-05  0.000528
#>  race   black - white  0.006047 -5.00e-02  0.059730
#>  race   other - white  0.024162 -5.21e-02  0.105626
#>  resp   +1             0.000480 -1.39e-03  0.002773
#>  rhc    1 - 0          0.059675  1.18e-02  0.111183
#>  sex    male - female  0.018854 -2.81e-02  0.059946
#>  surv2m +1            -0.177557 -2.60e-01 -0.114545
#> 
#> Type: response

Contrast says what change was made.1 For a categorical predictor, it is a difference between two levels, and for a numeric one it is an increase of one unit, which is a default rather than anything the data suggested. The comparisons are on the probability scale, so rhc reads as an increase of about six percentage points in the probability of death and age as three tenths of a point per year of age. A logistic regression would report each of these as one slope on the log-odds scale; these are averages over the sample of a quantity the model allows to differ from patient to patient.

The one-unit default should be reconsidered whenever a predictor does not span a unit. surv2m is a probability, so an increase of one is wider than its whole observed range, and the model holds its prediction flat past the edge of that range rather than continuing any trend. Asking instead for a change the data contains gives a more useful answer:

marginaleffects::avg_comparisons(fit, variables = list(surv2m = "iqr"))
#> 
#>  Estimate  2.5 % 97.5 %
#>    -0.282 -0.355 -0.207
#> 
#> Term: surv2m
#> Type: response
#> Comparison: Q3 - Q1

Moving a patient from the first quartile of the prognostic score to the third lowers the predicted probability of death by about .28, against the .18 the one-unit contrast reported.

The Shape of a Relationship

One number for a predictor hides the shape of the relationship behind it, and the two surv2m contrasts show it happening when the shape matters: the answer depends on which change is asked about. Plotting the fit against one predictor shows the whole curve, averaging over the other predictors at each value. This can be done using partial_dependence() and plotting its output:


pd <- partial_dependence(fit, ~ surv2m)

pd
#> Partial dependence
#> 
#> Predictor: "surv2m"
#> Averaged over 1500 units, on the "response" scale
#> 
#>  surv2m estimate lower upper
#>  0.0000    0.859 0.806 0.907
#>  0.0376    0.858 0.805 0.905
#>  0.0752    0.858 0.805 0.904
#>  0.1128    0.857 0.806 0.903
#>  0.1504    0.856 0.807 0.903
#>   --- 16 rows omitted ---
#>  0.7896    0.495 0.434 0.554
#>  0.8272    0.483 0.419 0.548
#>  0.8648    0.477 0.404 0.547
#>  0.9024    0.475 0.395 0.546
#>  0.9400    0.474 0.393 0.546
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ `n_print` in `print()` (`?bartisan::print.bartisan_partial()`) sets how many
#>   rows are shown, half from each end; `print(., n_print = Inf)` shows all of
#>   them.

plot(pd) +
  ggplot2::labs(x = "Estimated probability of surviving two months",
                y = "Fitted probability of death")

The fitted probability of death falls from about .85 to just under .5 as the prognostic score rises, and the fall is not a straight line, which is why the two contrasts above disagree about its size. A logistic regression reports one slope on the log-odds scale for the whole range. Nothing had to be specified to find the shape.

The band is a credible interval on the average prediction at each value, not on any one patient’s, and it widens at the top where few patients were that healthy. partial_dependence() returns the same numbers without drawing them, and marginaleffects::plot_predictions() is the one to reach for when the grid or what is conditioned on needs more control.

vignette("effects") covers effects, curves, and interactions.

The Effect of a Treatment

rhc has a row in the table above, so in one sense its effect has already been reported. estimate_effect() asks the same question in the vocabulary of a causal analysis: it takes the treatment as an argument rather than as one predictor among many, it names the estimand being averaged, and it reports the two potential outcomes the contrast is a difference of.

eff <- estimate_effect(fit, treat = "rhc")

eff
#> Average treatment effect (difference)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> 
#>     contrast estimate  lower upper    n
#>  Y[1] - Y[0]   0.0608 0.0118 0.111 1500
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.631 0.602 0.661
#>      Y[1]    0.692 0.656 0.733
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".

Catheterization is associated with an increase of about six percentage points in the probability of death. The interval runs from about one point to roughly eleven, so the direction is reasonably clear and the size is not. The estimate differs a little from the rhc row of the comparisons table because it is the posterior mean of the same draws rather than their median.

Because the outcome is binary, this is a difference in probability, which is interpretable without reference to the model, and it is usually the number to report. The two probabilities it is a difference of are printed below it: about 63% of patients would be expected to die without catheterization and 69% with it, averaging over the covariates as they actually occur in this sample.

The treatment has to be categorical, since the contrast is taken between its levels. The effect of a continuous treatment is a slope or a dose-response curve instead, and estimate_effect() says so and names the functions that produce one.

A single forest with the treatment among its predictors is not necessarily the best model for estimating a causal effect. Bayesian causal forests (BCF), which give the treatment effect a forest of its own, are described in vignette("causal"), which also covers what has to be true of the data before any of this can be read as an effect of the procedure.

Predicting New Observations

predict() returns the posterior mean prediction.

new_patient <- rhc[1, ]
new_patient$rhc <- 1

predict(fit, newdata = new_patient)
#> [1] 0.8366

For a prediction with an interval, use marginaleffects::predictions(), which reports the posterior median rather than the mean predict() reports:

marginaleffects::predictions(fit, newdata = new_patient)
#> 
#>  Estimate 2.5 % 97.5 %
#>     0.841 0.732  0.911
#> 
#> Type: response

This describes the probability that a patient with these characteristics dies. It is an interval for that probability, not a statement about which way any individual patient will go: the outcome itself is either 0 or 1, and a probability of .7 is entirely compatible with survival.

Comparing Models

Approximate leave-one-out cross-validation estimates how well a model predicts data it has not seen.

library(loo)

loo(fit)
#> 
#> Computed from 800 by 1500 log-likelihood matrix.
#> 
#>          Estimate   SE
#> elpd_loo   -847.6 17.3
#> p_loo        35.7  1.0
#> looic      1695.2 34.6
#> ------
#> MCSE of elpd_loo is 0.5.
#> MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.5]).
#> 
#> All Pareto k estimates are good (k < 0.66).
#> See help('pareto-k-diagnostic') for details.

elpd_loo is the estimated log predictive density on held-out data, where higher is better. p_loo is the effective number of parameters, which is a measure of how much of the data the forest is actually using. The Pareto k diagnostics are all good, meaning the approximation is trustworthy for this fit.

Two models can be compared directly. Here we ask whether everything recorded about how sick the patient was on arrival, the seven physiological measurements and the prognostic score among them, earns its keep over the demographics alone:

set.seed(2026)
demographics <- bartisan(death ~ rhc + age + sex + race + edu, data = rhc,
                         family = binomial())

loo_compare(list(full = loo(fit),
                 demographics = loo(demographics)))
#>         model elpd_diff se_diff p_worse diag_diff diag_elpd
#>          full       0.0     0.0      NA                    
#>  demographics     -67.9    11.2    1.00

The full model predicts better by around six times the standard error of the difference, as we would expect: how sick a patient is on arrival is the main thing that predicts whether they die.

vignette("comparison") covers model comparison, and the cases where leave-one-out fails.

Limitations

The model is flexible about the shape of the relationship and nothing else.

It does not make an association causal. Patients were not randomized to catheterization; sicker patients were more likely to receive it, which is exactly the kind of confounding that can produce an apparent harm. Whether the estimate above can be read as the effect of the procedure depends on whether the covariates account for that selection, which is a question about the study and not about the fit. vignette("causal") covers what is required, using this same data.

It does not extrapolate reliably. A predictor value outside the range of the training data is mapped to the edge of that range, so the prediction is held flat at whatever the fitted function is at the edge rather than continuing any trend.

It does not fix a badly chosen family. Getting the outcome distribution wrong matters more than any sampler setting.

Further Reading

Topic Vignette
How BART works, and the sampler vignette("implementation")
Choosing a family vignette("families")
Convergence and fit vignette("diagnostics")
Variable importance and selection vignette("importance")
Effects, Curves, and Interactions vignette("effects")
Varying coefficients and random intercepts vignette("varying")
Model comparison vignette("comparison")
Causal inference vignette("causal")
Censored and survival outcomes vignette("survival")
Frequently asked questions vignette("faq")

?bartisan-package has a longer map, organized by task and pointing at functions rather than vignettes.

References

Connors, Alfred F., Theodore Speroff, Neal V. Dawson, et al. 1996. “The Effectiveness of Right Heart Catheterization in the Initial Care of Critically Ill Patients.” JAMA 276 (11): 889–97. https://doi.org/10.1001/jama.1996.03540110043030.

  1. By default marginaleffects reports the posterior median as the point estimate, whereas predict() and estimate_effect() report the posterior mean, so the same quantity can differ slightly between them. Setting options(marginaleffects_posterior_center = mean) before the call switches it to the mean; we do this in vignette("causal").↩︎

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.