---
title: "Getting Started with bartisan"
author: "Noah Greifer"
date: "`r Sys.Date()`"
output:
  html_vignette:
    toc: yes
vignette: >
  %\VignetteIndexEntry{Getting Started with bartisan}
  %\VignetteEngine{knitr::rmarkdown_notangle}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
link-citations: true
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5,
                      fig.height = 3.2, dpi = 110, fig.align = "center",
                      message = FALSE, warning = FALSE)
options(digits = 4, marginaleffects_safe = FALSE)

# marginaleffects summarizes posterior draws with collapse, so both are needed.
has_me <- rlang::is_installed(c("marginaleffects", "collapse"))
has_loo <- rlang::is_installed("loo")
has_bp <- rlang::is_installed("bayesplot")

# The fits below are too slow to run on CRAN's machines, so their output is
# saved ahead of time and replayed where the code is not run, as on CRAN; see
# saved-results.R.
source("saved-results.R", local = knitr::knit_global())
knitr::opts_chunk$set(saved = TRUE)
```

## 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.

```{r libs}
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*](https://CRAN.R-project.org/package=future) package. To request multi-session evaluation, one can simply call the following:

```{r, eval = FALSE}
future::plan("multisession")
```

The [*future.apply*](https://CRAN.R-project.org/package=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*](https://CRAN.R-project.org/package=progressr) package. One can request progress bars globally for all functions that produce them using the following:

```{r, eval = FALSE}
progressr::handlers(global = TRUE)
```

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

```{r, eval = FALSE}
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 [@connors1996]. The question the study asked is
whether the procedure helps or harms.

```{r data}
data(rhc)

str(rhc)
```

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).

```{r fit}
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
```

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:

```{r}
summary(fit)
```

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.

```{r diagnose}
diagnose(fit)
```

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.

```{r diagnose_effect}
diagnose(estimate_effect(fit, treat = "rhc"))
```

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.

```{r ppcheck, eval = has_bp && has_loo, fig.height = 3.0}
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.

```{r varimp}
variable_importance(fit)
```

`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.

```{r comparisons, eval = has_me}
marginaleffects::avg_comparisons(fit)
```

`Contrast` says what change was made.^[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")`.]
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:

```{r comparisons_iqr, eval = has_me}
marginaleffects::avg_comparisons(fit, variables = list(surv2m = "iqr"))
```

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:

```{r pdp, fig.height = 3.0}

pd <- partial_dependence(fit, ~ surv2m)

pd

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.

```{r effects}
eff <- estimate_effect(fit, treat = "rhc")

eff
```

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.

```{r pred}
new_patient <- rhc[1, ]
new_patient$rhc <- 1

predict(fit, newdata = new_patient)
```

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

```{r predint, eval = has_me}
marginaleffects::predictions(fit, newdata = new_patient)
```

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.

```{r loo, eval = has_loo}
library(loo)

loo(fit)
```

`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:

```{r loocompare, eval = has_loo}
set.seed(2026)
demographics <- bartisan(death ~ rhc + age + sex + race + edu, data = rhc,
                         family = binomial())

loo_compare(list(full = loo(fit),
                 demographics = loo(demographics)))
```

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