Effects, Curves, and Interactions

Introduction

A forest has no coefficients: there is no table of slopes to read, and no standard error to put beside one. This is the part of the workflow that changes most when moving from parametric models to BART. At the same time, this is the very reason to switch to BART: the absence of coefficients accompanies to the extreme flexibility BART has in modeling relationships without the analyst having to specify their form.

Though this may seem like a loss, modern approaches to characterizing variable relationships can support a model-agnostic workflow (Rohrer and Arel-Bundock 2026). The replacement for coefficients is to ask the fitted model questions about predictions. What does it predict for these people? What would it predict if this variable were different? How much does the answer differ between groups? The marginaleffects package (Arel-Bundock et al. 2024) answers all of them, and returns posterior intervals with those answers.

Two of those questions bartisan answers itself: estimate_effect() considers one binary or factor treatment and reports the contrast as the average treatment effect or as one effect per unit, optionally within subgroups; partial_dependence() averages the fitted surface over the sample at each value of one or two predictors. vignette("causal") contains the worked version of the first. Everything else in this vignette uses marginaleffects: slopes, contrasts between covariate values, a continuous predictor, hypotheses comparing one estimate to another, and any grid that fixes the other covariates at chosen values rather than averaging over them.

This vignette covers the questions worth asking and how to phrase them. vignette("bartisan") is the shorter tour, and this expands its section on interpreting the fit.

In this guide, we will start from the three kinds of question the package answers and then work through them in turn. First, we’ll consider average effects of categorical and numeric predictors, along with the prior setting (sparsity) that can quietly attenuate a contrast when it is turned on. Next we’ll split those effects by subgroup and test a moderation with the difference of the two, then plot the shape of a fitted relationship and choose the scale on which an effect is reported, and close with what these estimates can and cannot be taken to mean.

Below, we load in the rhc dataset (see vignette("bartisan") or help("rhc", package = "bartisan") for details) and fit a BART logistic regression model predicting death from the available predictors. (Normally we would use multiple chains, perform diagnostics on mixing, and possibly modify the fitting parameters; we skip those here for brevity.)

library(bartisan)
library(marginaleffects)

data("rhc")

set.seed(2026)

# Fit a BART logistic regression model
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

Predictions, Comparisons, and Slopes

Everything below is one of the following:

A prediction is the outcome the model expects for a set of covariate values. predictions() gives one per unit or covariate profile, avg_predictions() averages them.

A comparison is the difference between two predictions that differ in one variable. This is the closest thing to a regression coefficient, and we often want to report it.

A slope is the derivative of the prediction with respect to a numeric variable. These are often also used to report but come with some limitations due to some peculiarities of using a tree-based model.

Average Effects (avg_comparisons())

avg_comparisons() with no variables argument gives every predictor at once, which for this model is a long table. A few at a time is easier to read, and often only a few effects are of primary interest anyway (Keele et al. 2020):

comp <- avg_comparisons(fit, variables = c("rhc", "age", "card"))

comp
#> 
#>  Term Contrast Estimate    2.5 %  97.5 %
#>  age  +1        0.00309  0.00153 0.00467
#>  card yes - no  0.04164 -0.00868 0.08973
#>  rhc  1 - 0     0.05968  0.01182 0.11118
#> 
#> Type: response

This reads as a coefficient table. Each estimate is an average difference in predicted probability, holding everything else at each patient’s own values: for a factor, between the levels named in the Contrast column, and for a numeric predictor, for an increase of one unit. rhc is coded 0 and 1, so its one-unit contrast is the treatment effect.

The comparison for age, for example, means that increasing all units’ age by 1, keeping all other predictors at their observed values, would yield an increase in the overall probability of death of 0.309 percentage points (a very small effect on this scale), with a 95% credible interval excluding 0. This isn’t a causal estimate; it’s just what the model predicts would happen.

avg_comparisons(fit, variables = list(aps = 10))
#> 
#>  Estimate   2.5 % 97.5 %
#>    0.0157 0.00171 0.0309
#> 
#> Term: aps
#> Type: response
#> Comparison: +10

Effects for Subgroups (by)

Supplying by to avg_comparisons() and other functions in marginaleffects splits the average by a grouping variable.

avg_comparisons(fit, variables = "rhc", by = "card")
#> 
#>  card Estimate  2.5 % 97.5 %
#>   no    0.0616 0.0115  0.116
#>   yes   0.0558 0.0037  0.110
#> 
#> Term: rhc
#> Type: response
#> Comparison: 1 - 0

The two subgroup estimates are close, and both intervals exclude zero by a small margin. A common mistake is to stop here and conclude that the effect differs between groups; that comparison is not a test. The question is whether the two effects differ from each other, which needs the difference of the two (i.e., a difference of differences) with an interval of its own. This can be requested by specifying hypothesis = ~pairwise.

avg_comparisons(fit, variables = "rhc", by = "card",
                hypothesis = ~pairwise)
#> 
#>    Hypothesis Estimate   2.5 % 97.5 %
#>  (yes) - (no) -0.00259 -0.0611 0.0248
#> 
#> Type: response

The difference is small with an interval covering zero, so there is no evidence here that the effect of catheterization varies by cardiovascular disease status. This is how to test an interaction in a model that never had an interaction term to test, and the interval on that difference is the only thing that separates a real interaction from two subgroup estimates that merely look different.

avg_predictions() does the same thing for predictions rather than differences, which is useful for describing groups:

avg_predictions(fit, by = "card")
#> 
#>  card Estimate 2.5 % 97.5 %
#>   no     0.630 0.604  0.661
#>   yes    0.699 0.660  0.733
#> 
#> Type: response

These correspond to the fitted predictions among those in each cardiovascular disease status group. That answers a different question from whether the model predicts that changing units’ disease status changes the probability of death, which corresponds to running the following:

avg_predictions(fit, variables = "card")
#> 
#>  card Estimate 2.5 % 97.5 %
#>   no     0.639 0.613  0.670
#>   yes    0.682 0.639  0.719
#> 
#> Type: response

In the latter, avg_predictions() computes predictions for all units, first setting card to "no" and then to "yes".

The Shape of a Relationship (partial_dependence())

Averages can hide the shape of the relationship between a predictor and the outcome. To see the fitted function we plot predictions against one predictor while the others are averaged over. partial_dependence() does this: it sets the predictor to each value of a grid in turn, predicts for every unit, and averages within each posterior draw, so the interval that comes back is on the average prediction rather than on any one patient’s. partial_dependence(fit, ~ x) builds the grid, and plot() draws it.

library(ggplot2)

pd <- partial_dependence(fit, ~aps)

pd
#> Partial dependence
#> 
#> Predictor: "aps"
#> Averaged over 1500 units, on the "response" scale
#> 
#>     aps estimate lower upper
#>    4.00    0.592 0.491 0.657
#>    9.72    0.593 0.495 0.657
#>   15.44    0.593 0.499 0.657
#>   21.16    0.597 0.513 0.656
#>   26.88    0.606 0.543 0.658
#>   --- 16 rows omitted ---
#>  124.12    0.724 0.651 0.811
#>  129.84    0.724 0.651 0.811
#>  135.56    0.724 0.651 0.811
#>  141.28    0.724 0.651 0.812
#>  147.00    0.724 0.651 0.812
#> 
#> ℹ 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) +
  labs(x = "APACHE III score on day 1", y = "Fitted probability of death")

# Same thing:
## plot(fit, ~aps)

The probability of death rises with the illness score, and the rise is not a straight line on any scale the model was told about; nothing was specified to find the shape. plot() returns a ggplot2 object, so the labels above are added the way they would be to any other plot.

The band is a credible interval and is wide at the top, where few patients were that sick, so the ends of a curve deserve more caution than the middle: with little data out there the forest shrinks its predictions toward the overall mean, which flattens the curve at both edges of the predictor’s range.

Naming a second predictor shows how the shape differs across groups.

pd2 <- partial_dependence(fit, ~ aps + rhc)

pd2
#> Partial dependence
#> 
#> Predictors: "aps" and "rhc"
#> Averaged over 1500 units, on the "response" scale
#> 
#>     aps rhc estimate lower upper
#>    4.00   0    0.568 0.460 0.638
#>    9.72   0    0.568 0.462 0.637
#>   15.44   0    0.569 0.472 0.637
#>   21.16   0    0.573 0.493 0.637
#>   26.88   0    0.581 0.518 0.640
#>     --- 42 rows omitted ---
#>  124.12   1    0.758 0.680 0.840
#>  129.84   1    0.758 0.680 0.840
#>  135.56   1    0.759 0.680 0.841
#>  141.28   1    0.759 0.680 0.841
#>  147.00   1    0.759 0.680 0.841
#> 
#> ℹ 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(pd2) +
  labs(x = "APACHE III score on day 1", y = "Fitted probability of death",
       color = "Catheterized", fill = "Catheterized")

The two curves run close together; if they diverged, that would be a moderation worth reporting, and the difference of differences above is how to put a number on it.

plot(fit, ~ meanbp + aps) +
  labs(x = "Mean arterial blood pressure", y = "Fitted probability of death",
       color = "APACHE III", fill = "APACHE III")
#> ℹ Grouping by `aps` at 41, 54, and 68, three of its values near its quartiles.
#> ℹ Set `values` to choose them yourself.

plot(fit, ~ meanbp + aps,
     values = list(aps = function(x) quantile(x, c(.05, .5, .95)))) +
  labs(x = "Mean arterial blood pressure", y = "Fitted probability of death",
       color = "APACHE III", fill = "APACHE III")

at <- quantile(rhc$aps, c(.1, .5, .9))

# Uses posterior mean
partial_dependence(fit, ~ aps, values = list(aps = at))
#> Partial dependence
#> 
#> Predictor: "aps"
#> Averaged over 1500 units, on the "response" scale
#> 
#>   aps estimate lower upper
#>  29.9    0.610 0.552 0.657
#>  54.0    0.656 0.623 0.695
#>  83.0    0.703 0.649 0.765
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.

# Uses posterior median by default
avg_predictions(fit, variables = list(aps = at))
#> 
#>   aps Estimate 2.5 % 97.5 %
#>  29.9    0.611 0.552  0.657
#>  54.0    0.656 0.623  0.695
#>  83.0    0.701 0.649  0.765
#> 
#> Type: response

plot_predictions(fit, condition = list(aps = at), draw = FALSE)
#>   rowid estimate conf.low conf.high  df   age card  crea   edu  hema meanbp
#> 1     1   0.6299   0.4969    0.7594 Inf 61.42   no 2.132 11.64 31.68     78
#> 2     2   0.6817   0.5604    0.7867 Inf 61.42   no 2.132 11.64 31.68     78
#> 3     3   0.7386   0.6018    0.8416 Inf 61.42   no 2.132 11.64 31.68     78
#>   paco2  pafi  race resp rhc  sex surv2m  aps
#> 1 38.85 217.4 white   28   0 male  0.587 29.9
#> 2 38.85 217.4 white   28   0 male  0.587 54.0
#> 3 38.85 217.4 white   28   0 male  0.587 83.0

Slopes and the Predictor Transform (avg_slopes())

marginaleffects::avg_slopes() reports the derivative of the regression function with respect to a predictor; this is closest quantity to a regression slope. However, whether that derivative means anything depends on gate and x_transform in bartisan_control(). Under a soft gate, the fit is a smooth function of the (transformed) predictor; with gate = "hard", the fit is piecewise constant under any transform and has no derivative worth taking.

Because a numeric predictor is transformed to the interval \([0, 1]\) according to x_transform, how its effect translates into a slope depends on the value supplied. Writing the fit as \(f(T(x))\), a slope is \(f'(T(x))\,T'(x)\), where \(T(\cdot)\) is the transform function specified in x_transform: an affine \(T(\cdot)\) has a known constant derivative, so only \(f'\) is estimated.

The default, "smoothcdf", is differentiable, so avg_slopes() works, but any error in estimating the density of \(X\) yields additional error in the estimation of the slope. Under "range", which is an affine map from the covariate’s observed range to \([0, 1]\), slopes are more accurate. "quantile" maps each predictor through its empirical distribution function, which is a step function, so the fitted function is a step function of the original predictor and the difference quotient grows without bound as the step shrinks. There is no derivative there to estimate, and the number avg_slopes() returns is a property of the step size rather than of the fit.

So, when a slope is the quantity being reported, it is best to fit the model with x_transform = "range", though the default can be okay as well, especially with a large sample (which yields less error in the estimation of the predictor distribution). For everything else, it makes most sense to use avg_comparisons() with a step we can interpret, as with aps = 10 above, which evaluates the fit at two points a substantive distance apart rather than dividing by a vanishing one. This is documented at ?bartisan-marginaleffects.

Descriptions of the Fitted Model

Everything here is a description of the fitted model. avg_comparisons() reports what the model predicts would differ between two versions of the data, which is a causal quantity only if the model contains enough covariates to account for confounding. Patients were not randomized to catheterization, so for this fit that is a strong assumption. vignette("causal") covers what is needed.

The intervals are posterior credible intervals under the model: they cover the uncertainty in the fitted function, and they do not cover the possibility that the model is missing a confounder, that the outcome is measured with bias, or that the sample is not the population of interest.

Further Reading

vignette("varying") covers vc(), which writes an effect into the model as a parameter with a prior of its own rather than reading it off as a contrast, and is where to go when the effect is the point of the model. vignette("importance") covers which predictors the forest uses, which is a different question from how much they move the outcome. vignette("diagnostics") covers whether the fit can be trusted before any of this is read. ?bartisan-marginaleffects documents which marginaleffects functions are supported and the arguments that are specific to this package, and ?estimate_effect and ?partial_dependence document the two questions answered without it.