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.
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:
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:
To request a progress bar for a single call, wrap the function
evaluation in progressr::with_progress(), e.g.,
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.
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 warmupSetting 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.
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.
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.
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.
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.
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.
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.
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: responseContrast 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 - Q1Moving 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.
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.
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.
predict() returns the posterior mean prediction.
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: responseThis 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.
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.00The 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.
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.
| 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.
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").↩︎