Varying Coefficients and Random Intercepts
Noah Greifer
2026-10-18
Source:vignettes/varying.Rmd
varying.RmdIntroduction
A BART model is a sum of trees, and a sum of trees has no
coefficients. It predicts, and an effect is read off it as a contrast
between two predictions, as estimate_effect() and the
marginaleffects functions in vignette("effects")
do. That is the right reading for most predictors. For some, though, we
want the effect to be a parameter: something with a prior of
its own, a posterior of its own, and a value we can report as a
coefficient rather than reconstruct as a difference. The treatment in a
causal analysis is the usual case, and so is any predictor whose effect
is the whole point of the model.
vc() is how a predictor is given a coefficient in
bartisan, and the coefficient it gets is a function rather than
a number: a forest of its own that says how the effect moves with the
other predictors. This is the varying-coefficient model of Deshpande et al. (2026), of which the Bayesian causal
forest of Hahn et al. (2020) is the case of one binary
covariate and Woody et al. (2020) the case of one continuous one.
bcf() fits that special case with the settings the causal
problem wants, and this guide is where the general model is
described.
In this guide we will write the model down and say what each of its
parts is, then declare a varying coefficient with vc() and
work through its three arguments: which predictors the coefficient
varies with, where the control function is read, and what happens when
the coefficient is not allowed to vary at all, which is how a linear
term is written into an otherwise nonparametric model. Next we’ll read
the fit, both as coefficients and as effects, and see how the settings
of a fit are keyed by forest once there is more than one. Then we’ll
show how bcf() writes its model in these terms. Finally,
because they are declared in the same formula and explained nowhere
else, we’ll cover random intercepts.
Throughout we use the rhc dataset (see
vignette("bartisan") or
help("rhc", package = "bartisan")), predicting
death from the covariates and from rhc,
whether the patient received right heart catheterization. The fits use
the default single chain; vignette("diagnostics") covers
checking convergence, which we skip here for brevity.
The Model
A model with varying coefficients is
\[g(\mu_i) = f_0(Z_i) + \sum_j (X_{ij} - c_j)\, f_j(Z_i)\]
where \(g\) is the link, \(Z_i\) is the vector of predictors, and each \(X_{ij}\) is a predictor that has been given a coefficient. \(f_0\) is a forest we will call the control function, following Hahn et al. (2020): the surface for the outcome when every \(X_j\) is at its centering value \(c_j\). Each \(f_j\) is a forest of its own, and it is the coefficient on \(X_j\): how much the linear predictor moves per unit of \(X_j\), as a function of whatever \(f_j\) is allowed to split on. With a single binary \(X\), \(f_1\) is the conditional treatment effect on the link scale, one value per covariate profile.
All of the forests are fitted at once. That gives the coefficient a prior of its own. Were \(X_j\) simply one predictor among many in a single forest, its effect would be whatever difference that forest happened to produce, regularized by the same prior that regularizes everything else, and Hahn et al. (2020) show that shrinking the prognostic surface then shrinks the effect along with it. With \(f_j\) separate, it can have a prior of its own, and in particular a tighter one, since the ways an effect varies are usually simpler than the surface it sits on.
Declaring a Varying Coefficient (vc())
vc() is a marker written inside the formula. It is never
evaluated as a function (calling it directly is an error) but is read
out of the formula by bartisan(), which builds the forests
it asks for.
set.seed(2026)
fit_vc <- bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card +
vc(rhc),
data = rhc, family = binomial())
fit_vc
#> Generalized BART
#>
#> Call:
#> bartisan(formula = death ~ age + sex + race + edu + aps + meanbp +
#> resp + hema + pafi + paco2 + crea + surv2m + card + vc(rhc),
#> data = rhc, family = binomial())
#>
#> Family: "binomial" with the "logit" link
#> Observations: 1500
#> Structure: 2 forests of 50 trees, soft decision rules
#> Draws: 800 kept after 200 warmupThe print shows two forests where a fit without vc()
would show one. These forests are the control function, named
(Intercept), and the coefficient, named for its covariate.
Both are fitted with the same number of trees here, which is the default
and is taken up below.
A predictor that has been given a coefficient should be kept out of
the control function. Were rhc among \(f_0\)’s predictors as well as multiplying
\(f_1\), the two would not be
separately identified, since any function of rhc can move
from one to the other without changing a single prediction.
We can examine the structure of the forest for the rhc
coefficient using variable_importance(), which gives a
table for each forest in the model:
variable_importance(fit_vc)
#> Variable importance
#>
#> predictor variable prop_used prop_splits splits
#> (Intercept) surv2m 1.000 0.117 8.9
#> (Intercept) age 1.000 0.097 7.4
#> (Intercept) paco2 1.000 0.084 6.4
#> rhc aps 1.000 0.082 6.2
#> (Intercept) pafi 1.000 0.081 6.1
#> rhc race 1.000 0.079 6.0
#> rhc meanbp 1.000 0.079 6.0
#> rhc resp 1.000 0.078 5.9
#> rhc card 1.000 0.077 5.8
#> (Intercept) race 1.000 0.067 5.1
#> rhc hema 0.999 0.083 6.3
#> (Intercept) aps 0.999 0.082 6.3
#> rhc surv2m 0.999 0.075 5.7
#> (Intercept) edu 0.999 0.073 5.6
#> (Intercept) card 0.999 0.069 5.2
#> rhc age 0.998 0.077 5.9
#> (Intercept) hema 0.998 0.071 5.4
#> (Intercept) meanbp 0.998 0.071 5.4
#> (Intercept) sex 0.998 0.064 4.9
#> rhc pafi 0.996 0.076 5.7
#> rhc sex 0.996 0.074 5.6
#> rhc paco2 0.996 0.073 5.5
#> rhc edu 0.995 0.075 5.7
#> rhc crea 0.995 0.073 5.5
#> (Intercept) crea 0.989 0.065 4.9
#> (Intercept) resp 0.984 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.Under the default splitting prior every predictor is used in nearly
every draw of both forests, so the prop_used column says
little here and the splits column carries the ranking:
surv2m and age take the most rules in the
control function ((Intercept)) forest, and
hema in the varying coefficient (rhc)
forest.
Effect Modifiers (modifiers)
By default, a coefficient is allowed to vary with every predictor in
the model except the covariate itself. Supplying an argument to
modifiers, the second argument to vc() taking
a one-sided formula, allows one to control which predictors enter the
coefficient forest.
set.seed(2026)
fit_aps <- bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card +
vc(rhc, ~ aps + age),
data = rhc, family = binomial())Now the effect of catheterization may differ by severity of illness
(aps), age, and nothing else, which is a
statement about the model and not a restriction the data can undo.
Restricting the modifiers is worth doing when there is reason to think
the effect depends only on a few specific predictors, since a forest
asked to search fourteen predictors for heterogeneity that lives in two
will spend some of its prior on the other twelve.
Calling variable_importance() on fit_aps,
we can see that the only variables used in the varying coefficient
forest are indeed aps and age, as
requested:
variable_importance(fit_aps)
#> Variable importance
#>
#> predictor variable prop_used prop_splits splits
#> rhc aps 1.000 0.500 37.7
#> rhc age 1.000 0.500 37.6
#> (Intercept) surv2m 1.000 0.122 9.2
#> (Intercept) age 1.000 0.096 7.2
#> (Intercept) paco2 1.000 0.081 6.1
#> (Intercept) aps 0.999 0.082 6.2
#> (Intercept) hema 0.999 0.073 5.5
#> (Intercept) edu 0.999 0.069 5.2
#> (Intercept) pafi 0.998 0.074 5.6
#> (Intercept) meanbp 0.996 0.073 5.5
#> (Intercept) crea 0.996 0.070 5.3
#> (Intercept) sex 0.996 0.065 4.9
#> (Intercept) card 0.994 0.066 5.0
#> (Intercept) resp 0.993 0.064 4.8
#> (Intercept) race 0.990 0.064 4.8
#> rhc sex 0.000 0.000 0.0
#> rhc race 0.000 0.000 0.0
#> rhc edu 0.000 0.000 0.0
#> rhc meanbp 0.000 0.000 0.0
#> rhc resp 0.000 0.000 0.0
#> rhc hema 0.000 0.000 0.0
#> rhc pafi 0.000 0.000 0.0
#> rhc paco2 0.000 0.000 0.0
#> rhc crea 0.000 0.000 0.0
#> rhc surv2m 0.000 0.000 0.0
#> rhc card 0.000 0.000 0.0
#>
#> ℹ 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.When the modified variable is continuous, modifiers
decides whether the effect is linear in the covariate. As written, the
model is \(z\, f_1(Z)\), which for a
binary \(z\) is no assumption at all,
there being only two values, and for a continuous one says the effect is
conditionally linear. That is a real restriction, and a numeric
covariate is allowed to modify its own coefficient to remove it:
# The effect of `aps` may itself change across `aps`, so the dose response is
# a curve rather than a line through the origin
bartisan(death ~ age + sex + vc(aps, ~ . + aps),
data = rhc, family = binomial())This fits \(\text{aps} \cdot
f_1(\text{aps},\, \text{age},\, \text{sex})\), so the slope moves
across the range of aps, and it is worth reaching for
whenever the effect of a continuous predictor might not be proportional
to it. A categorical covariate is removed from its own forests instead:
a level’s indicator is nonzero only on the rows where that level holds,
and the variable is constant on exactly those rows, so a split on it
would separate rows that contribute from rows that contribute
nothing.
Centering the Covariate (center)
The \(c_j\) in the model are the
values at which the control function is read, and the
center argument chooses them. This is a reparameterization
of \(f_0\) alone: every coefficient and
every estimand is identical under any choice, and the choice changes
only what the control function means and how well the two forests mix.
The default, "auto", picks by the covariate. A
0/1 covariate is left at zero, so \(f_0\) is the surface among the binary
predictor’s reference level, which is a quantity with a meaning of its
own. Any other numeric covariate is centered at its mean, because zero
may be nowhere near the data and a control function read there would be
an extrapolation. "zero", "mean",
"mid" (the midpoint of the range), or a specific number
override that.
A factor is handled differently. It is always fitted mean-centered
and gets one forest per level, coded symmetrically rather than as
contrasts against whichever level happened to sort first, and
coef() recenters the coefficients to sum to zero across the
levels. That coding leaves the reference level a reporting choice: the
center argument names the level to report against, and no
refit is needed to change it.
Setting center = "estimate" is different in kind. Rather
than subtract a number from the covariate, it gives each of the
covariate’s values a coefficient of its own and draws it, so
that every contrast carries the same prior whatever the number of values
and no value is a reference. At two values it restricts nothing, and
bcf() uses it for a binary treatment: with it, the answer
stops depending on which arm was written as 1. Above two values every
contrast becomes one shared shape times a scalar, where the symmetric
coding gives each level its own, so it is the parsimonious model against
a general one, and the one to reach for when the levels plausibly differ
in degree rather than in kind. It needs a covariate with between two and
twenty distinct values, and a family whose leaf target is quadratic in
the predictor it feeds, which is why bcf() falls back to
the fixed coding for a family that is not.
A Constant Coefficient (vc(z, ~ 1))
Naming no modifier at all leaves the coefficient’s forest nothing to split on. Every tree in it is then a stump, the coefficient is a single drawn number, and the covariate enters the model as a linear term while everything else stays nonparametric. That is the semiparametric case of the General BART model of Tan and Roy (2019), a forest and a parametric part taking disjoint sets of predictors, and it is the shape stan4bart and flexBART fit1.
Comparing the constant coefficient against the varying one is a
question neither variable importance nor a dropped predictor answers:
not whether rhc matters, but whether what it does depends
on anything else. loo_compare() from loo answers
it; vignette("comparison") covers what it is doing.
set.seed(2026)
constant <- bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card +
vc(rhc, ~ 1),
data = rhc, family = binomial())
constant
#> Generalized BART
#>
#> Call:
#> bartisan(formula = death ~ age + sex + race + edu + aps + meanbp +
#> resp + hema + pafi + paco2 + crea + surv2m + card + vc(rhc,
#> ~1), data = rhc, family = binomial())
#>
#> Family: "binomial" with the "logit" link
#> Observations: 1500
#> Structure: 2 forests of 50 trees, soft decision rules
#> Draws: 800 kept after 200 warmup
library(loo)
loo_compare(list(constant = loo(constant), varying = loo(fit_vc)))
#> model elpd_diff se_diff p_worse diag_diff diag_elpd
#> constant 0.0 0.0 NA
#> varying -0.5 1.0 0.68 |elpd_diff| < 4The constant coefficient is not beaten, so there is no evidence here that the effect of catheterization on the log odds of death varies with the other covariates. That buys an interpretation the varying fit cannot offer, because one number describes the whole sample:
b <- coef(constant, draws = TRUE)[["rhc"]][, 1]
c(log_OR = mean(b),
OR = exp(mean(b)),
lower = quantile(exp(b), .025, names = FALSE),
upper = quantile(exp(b), .975, names = FALSE)) |>
round(3)
#> log_OR OR lower upper
#> 0.338 1.402 1.059 1.799An odds ratio of about 1.4, conditional on the predictors in the forest. A constant conditional odds ratio is a much more clinically meaningful effect summary than a marginal effect measure because it describes the estimated effect for each patient profile as determined by the included covariates, as opposed to an average effect across a whole population. This is the same interpretation a coefficient in a logistic regression would have, but fit with BART, we can be reasonably confident that any nonlinearities and interactions among covariates are adjusted for, whereas a logistic regression must impose the model structure a researcher believes to be there.
It is important to remember that the coefficient is drawn under the leaf prior of its own forest and not under a prior written for a regression coefficient, so it is shrunk toward zero and is a regularized log odds ratio rather than a maximum likelihood one; because the leaf scale is divided by \(\sqrt{m}\) for a forest of \(m\) trees, that prior does not change with the number of trees. Also, the comparison between a fixed and a varying coefficient is far better powered on some outcomes than others. On a Gaussian outcome with a coefficient truly ranging from 0.3 to 2.5, the varying coefficient was easily preferred at \(n = 1000\); on a binary outcome with the same coefficients, both models were equally preferred. A null result on a binary outcome at a few hundred observations says little, and should not be read as evidence that an effect is constant.
Reading a Fit
The Coefficients (coef())
coef() returns the varying coefficients and nothing
else. The control function is a prediction rather than a coefficient,
and predict() reports it.
head(coef(fit_vc))
#> rhc
#> [1,] 0.3886
#> [2,] 0.1962
#> [3,] 0.2041
#> [4,] 0.2841
#> [5,] 0.2481
#> [6,] 0.2463coef() returns one row per observation and one column
per coefficient, since a coefficient allowed to vary has a value for
each observation: for the first patient, or more precisely for a patient
with the first patient’s covariate profile, the log odds of death move
by 0.389 under catheterization. The values are on the link scale, so for
a binary outcome they are differences in log odds rather than in
probability. Supplying newdata evaluates the coefficients
at other covariate profiles, and setting draws = TRUE
returns every posterior draw rather than the mean, as a list with one
draws-by-observations matrix per coefficient, which an interval on any
one patient’s coefficient needs.
draws <- coef(fit_vc, draws = TRUE)[["rhc"]]
# A 95% interval on the coefficient for the first three patients
t(apply(draws[, 1:3], 2L, quantile, c(.025, .975)))
#> 2.5% 97.5%
#> [1,] -0.002961 0.9387
#> [2,] -0.240273 0.6338
#> [3,] -0.214174 0.6488summary() on the fit reports the forests as it does for
any fit, and prior_summary() writes out the prior each one
was drawn under, including the leaf scale of the coefficient forest,
which is where the regularization that separates the effect from the
surface lives.
Effects From the Coefficients (estimate_effect())
A coefficient on the link scale is one step from an effect on the
response scale, and estimate_effect() makes it easier to
report meaningful contrasts. With treat naming the
covariate, it predicts outcomes under each of the predictor’s values and
contrasts them. Setting estimand = "CATE" returns one
effect per unit, here on the risk difference scale, and
plot() on that draws them ordered by their estimate with
the marginal effect as a separate interval past the right edge.
cate <- estimate_effect(fit_vc, treat = "rhc", estimand = "CATE")
plot(cate)
The marginal effect, effects by subgroup, or the effect as a ratio
are all the same call with different arguments. For example, setting
comparison = "or" displays the effects on the conditional
odds ratio scale. ?estimate_effect documents these
arguments, and vignette("causal") is where the assumptions
live that turn any of them into a causal effect.
Settings That Vary by Forest
Every forest has its own prior, and every setting in
bartisan_control() that could mean something different for
one forest than for another may be given once, to apply to all, or once
per forest. Given per forest, the values are either positional in the
order the forests are listed in the fit’s print, or keyed by the forest
names. Both of these ask for fifty trees in the control function and
twenty-five in the coefficient:
bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card +
vc(rhc),
data = rhc, family = binomial(), num_trees = c(50L, 25L))
bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card +
vc(rhc),
data = rhc, family = binomial(),
num_trees = c("(Intercept)" = 50L, rhc = 25L))A forest a named argument does not mention keeps that argument’s own
default rather than borrowing another forest’s value. The settings this
covers are num_trees, k,
sigma_mu, sparsity, split_prior,
bandwidth, gamma, beta, the four
alpha arguments, and the three update_ flags;
?bartisan_control lists them.
Fewer trees for a coefficient than for the control function is the
setting worth knowing about, and it is the one bcf() makes
by default. A coefficient forest describes how an effect varies, and
that is usually a simpler surface than the outcome’s, so it needs less
capacity; and with less capacity it competes less with the control
function for the same variation, which keeps the two well separated in
the posterior. When a fit with a varying coefficient mixes badly, a
coefficient forest that is too large is the first thing to try.
The sparsity argument is the other setting that reads
differently on a coefficient forest. On the control function it selects
among the predictors the way it does in any BART fit. On a coefficient
forest what it selects among is the modifiers, and dropping all
of them leaves an effect that does not vary rather than an effect that
is zero, since nothing can drop the covariate the forest multiplies.
That is the shrinkage a heterogeneity model wants, and it makes the
sparsity prior safe to turn on for a coefficient forest where it is not
for a single forest with the treatment among its predictors.
bcf() keeps the default of off for both forests, and
setting sparsity = c(FALSE, TRUE) turns it on for the
coefficient forest alone; vignette("effects") covers the
single-forest case.
Families With Several Additive Predictors
gaussian_ls(), zi_poisson() and the other
families with more than one additive predictor fit one control function
per parameter, and each parameter’s formula carries its own
vc() terms. The forests are then two-dimensional, a control
function and its coefficients for each parameter, and are named
accordingly:
# Forests: mean, mean:z, log_sd
bartisan(list(mean = y ~ x1 + x2 + vc(z),
log_sd = ~ x1 + x2),
data = d,
family = gaussian_ls())
# Forests: mean, mean:z, log_sd, log_sd:z. One formula reaches every parameter,
# which is the rule every per-forest argument follows too
bartisan(y ~ x1 + x2 + vc(z), data = d, family = gaussian_ls())
# Each coefficient with modifiers of its own
bartisan(list(mean = y ~ x1 + x2 + vc(z, ~ x2),
log_sd = ~ x1 + x2 + vc(z, ~ x1)),
data = d, family = gaussian_ls())So the same covariate may have a coefficient on more than one
parameter. z shifting the mean and z widening
the spread are different questions, and both are answered at once, and
coef() returns one column per coefficient, named
mean:z and log_sd:z for the forests they come
from.
Those names are also how the per-forest settings are keyed, and every
forest can be set at once. With four forests, a control function and a
coefficient for each parameter, num_trees takes four
values, either positionally in the order the forests are listed in the
fit’s print or by name:
# The mean's control function and coefficient, then the log standard
# deviation's, in the order the print lists them
bartisan(y ~ x1 + x2 + vc(z), data = d, family = gaussian_ls(),
num_trees = c(50L, 20L, 15L, 10L))
# The same, by name; the order then does not matter
bartisan(y ~ x1 + x2 + vc(z), data = d, family = gaussian_ls(),
num_trees = c("mean" = 50L, "mean:z" = 20L,
"log_sd" = 15L, "log_sd:z" = 10L))A forest a named vector leaves out keeps the argument’s own default
rather than taking its parameter’s value: setting
num_trees = c(mean = 50L, log_sd = 15L) gives both
coefficient forests the default of 50, not 50 and 15. Name the
coefficient forests too when they are meant to differ.
Bayesian Causal Forests (bcf())
The Bayesian causal forest of Hahn et al. (2020) is a varying coefficient model
with one binary covariate, the treatment, and a handful of choices made
for the causal problem. bcf() implements this model with
these modifications to support treatment effect estimation. Given a
formula for the covariates and the treatment named in
treat, it writes a bartisan() call that could
have been written directly:
set.seed(2026)
fit_bcf <- bcf(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card,
treat = ~ rhc, data = rhc, family = binomial())
fit_bcf
#> Generalized BART
#>
#> Call:
#> bcf(formula = death ~ age + sex + race + edu + aps + meanbp +
#> resp + hema + pafi + paco2 + crea + surv2m + card, treat = ~rhc,
#> data = rhc, family = binomial())
#>
#> Family: "binomial" with the "logit" link
#> Observations: 1500
#> Structure: 2 forests of 50 and 25 trees, soft decision rules
#> Draws: 800 kept after 200 warmup
#>
#> Posterior means: b.rhc.0 = -0.0312, b.rhc.1 = 0.0635
#>
#> Treatment: "rhc"
#> Effect moderators: "age", "sex", "race", "edu", "aps", "meanbp", "resp", "hema", "pafi", "paco2", "crea", "surv2m", and "card"
#> ℹ `estimate_effect()` reports the treatment effect, with the average potential
#> outcomes beside it; `plot()` draws the conditional ones.bcf() makes the following modifications to a standard
bartisan() call:
- The treatment gets a
vc()term, so its effect is a forest with a prior of its own, which is everything above. - The estimated propensity score, fitted by a logistic BART model on
the same covariates, goes into the control function as a predictor named
.propensityand not into the effect forest, so that the control function absorbs the selection into treatment while the effect stays free of it; that is the remedy Hahn et al. (2020) propose for what they call regularization-induced confounding. - The effect forest gets fewer trees than the control function, 25 against 50, for the reason given in the section above.
- A binary treatment’s coding is drawn rather than fixed, which is
center = "estimate", so the answer is the same whichever arm was written as 1.
The sparsity argument is left at its default of off for
both forests. Because the treatment is the coefficient rather than a
predictor a splitting proportion could drop, the prior can be turned on
for the coefficient forest alone with
sparsity = c(FALSE, TRUE), as the section above
describes.
The moderators argument is vc()’s
modifiers under the name the causal literature uses, and
the propensity model is bartisan() with
propensity_args for its settings.
Written out by hand, and leaving aside that the propensity score would have to be fitted and added to the data first, the model is
bartisan(death ~ age + sex + race + edu + aps +
meanbp + resp + hema + pafi +
paco2 + crea + surv2m + card + .propensity +
vc(rhc, center = "estimate"),
data = rhc, family = binomial(), num_trees = c(50L, 25L))and every method that works on a bartisan() fit works on
a bcf() fit unchanged, since the second is the first with
the treatment recorded. What the record adds is that
estimate_effect() needs no treat,
print() names the treatment and its moderators, and
plot() on the fit is
plot(estimate_effect(fit, estimand = "CATE")):
plot(fit_bcf)
The reasons to write the model by hand rather than through
bcf() are the ones this guide is about: a coefficient on
something other than a treatment, more than one of them, a continuous
covariate whose effect should not be linear in it, or a family with
several additive predictors. vignette("causal") covers what
has to be true of the design for any of these effects to be read
causally, which no model supplies.
Random Intercepts ((1 | group))
A (1 | group) term in the formula adds an intercept per
level of group, drawn from a common mean-zero normal whose
standard deviation is itself drawn, under the same half-Cauchy prior the
leaf scale uses (Gelman
2006; Polson and Scott 2012).
Several grouping factors are allowed, and (1 | a/b) expands
to nesting as it does in lme4.
The rhc data has no grouping structure, so we simulate
one: forty groups of about fifteen observations each, a nonlinear effect
of one covariate, and a group effect with a standard deviation of one
half.
set.seed(2026)
n_groups <- 40L
d <- data.frame(g = factor(sample(n_groups, 600L, replace = TRUE)),
x = runif(600L))
u <- rnorm(n_groups, 0, 0.5)
d$y <- sin(2 * pi * d$x) + u[d$g] + rnorm(600L, 0, 0.5)
fit_re <- bartisan(y ~ x + (1 | g), data = d)
fit_re
#> Generalized BART
#>
#> Call:
#> bartisan(formula = y ~ x + (1 | g), data = d)
#>
#> Family: "dpm" with the "identity" link
#> Observations: 600
#> Structure: 1 forest of 50 trees, soft decision rules
#> Draws: 800 kept after 200 warmup
#> Random intercepts: g (40 levels)
#>
#> Posterior means: alpha = 3.37, clusters = 15.7, center = -0.0182, error_sd = 0.521The intercepts are in fit_re$ranef and the standard
deviation they were drawn under in fit_re$tau, one matrix
of draws per additive predictor with a column per grouping factor, and
summary() reports that scale beside the leaf scale.
ranef() is the generic from nlme, which
lme4 re-exports, so either qualification reaches the method;
setting draws = TRUE gives the draws, which an interval on
a group needs.
# How much of the variation is between groups: the posterior of the group scale
quantile(fit_re$tau[[1L]][, "g"], c(.025, .5, .975))
#> 2.5% 50% 97.5%
#> 0.4220 0.5221 0.6656
# The intercepts, as posterior means, against the values they were drawn from
re <- nlme::ranef(fit_re)$g
cor(re[, "(Intercept)"], u)
#> [1] 0.9613
# With intervals, for the first four groups
apply(nlme::ranef(fit_re, draws = TRUE)$g[["(Intercept)"]][, 1:4], 2L,
quantile, c(.025, .975))
#> 1 2 3 4
#> 2.5% 0.3092 0.09548 -0.5970 -0.7757
#> 97.5% 0.8688 0.72409 0.0372 -0.2062A posterior mean is the wrong summary for a level with few observations, which is the case a group intercept exists for, and a level whose interval covers zero is one the data had little to say about. The intercepts are shrunk toward zero by their prior and are deviations from the additive predictor, so they come out approximately centered without being constrained to sum to zero exactly; nothing is lost by that, since the level of the fitted function is the additive predictor’s.
There are a few details about this term worth knowing. Only intercepts are supported, and a random slope is refused rather than ignored. A random intercept is a scalar entering the predictor with weight one for the observations in its level, exactly as a leaf is once its gate is removed, so the sampler’s leaf machinery handles it, closed forms included; a slope is a scalar multiplying a covariate, which is a different shape of parameter. A variable whose effect varies by group belongs in the fixed part of the formula, where a tree can split on the group and on the variable together and get an interaction of any shape. A grouping factor can also go in the fixed part, where a tree splits on it like anything else, and with few large groups that is the better choice: the group means are well determined without pooling, and a split can interact the group with the covariates. The random intercept wins where there are many small groups, which is where partial pooling earns its keep. A level not seen at fitting time is given the prior mean of zero when predicting, with a warning.
A family with several additive predictors gets a separate set of
intercepts for each, independent of one another, so a zero-inflated
count model has a group effect on the count part and another on the
inflation part, and the columns of ranef() are named for
the predictors rather than (Intercept). A group intercept
reaches every control function and no coefficient, since a group-varying
coefficient would be a random slope. prior_summary()
reports the prior the scale was drawn under.
Further Reading
?vc is the reference for the arguments used here, and
?bcf for what it writes into a bartisan() call
and how it handles a categorical or continuous treatment.
vignette("causal") covers the assumptions under which a
coefficient on a treatment is a causal effect, and
vignette("effects") the contrasts and curves that read an
effect off a fit without a coefficient.
vignette("comparison") covers loo_compare(),
on which the constant-against-varying comparison above rests.
?bartisan documents the (1 | group) syntax,
and ?ranef.bartisan_fit what comes back.