Response families for generalized BART

bartisan() accepts the stats::family objects used by stats::glm(), so gaussian(), binomial("probit"), poisson() and stats::Gamma() all work unchanged. The functions documented here supply the additional families that have no glm() counterpart, in the same style, so that they can be passed to the family argument the same way.

One caveat carries over from base R: stats::Gamma() defaults to link = "inverse", which is the worst link for this sampler, so write Gamma("log"), or name the family as the string family = "Gamma", which is this package's own spelling and resolves to the log link; see Details.

Usage

negbin(link = "log", theta = NULL)

ordinal(link = "logit")

multinomial(link = "logit", reference = NULL, replicates = 200L)

dpm_aft(
  nu = 10,
  q = 0.95,
  k_s = 10,
  alpha = NULL,
  max_clusters = NULL,
  psi = 0.5
)

dpm(nu = 10, q = 0.95, k_s = 10, alpha = NULL, max_clusters = NULL, psi = 0.5)

weibull_aft()

loglogistic_aft()

lognormal_aft()

ph(num_bins = NULL, lambda_shape = 1, update_lambda = TRUE)

gaussian_ls(link = "identity")

Gamma_ls(link = "log")

zi_poisson(link = "log")

zi_negbin(link = "log", theta = NULL)

Beta(link = "logit", phi = NULL)

ordbeta(link = "logit", phi = NULL)

tweedie(link = "log", power = 1.5, phi = NULL)

custom_family(
  logdens,
  num_predictors = 1L,
  start = 0,
  derivatives = NULL,
  aux_names = NULL,
  aux_start = NULL,
  name = "custom"
)

Arguments

link: string; the link function. Allowable options include "logit" (the default), "probit", and "cloglog" for ordinal() and Beta(), and "logit" (the default) and "probit" for multinomial(). The remaining families take one link each, which is therefore the default: "log" for negbin(), zi_poisson(), zi_negbin(), and Gamma_ls(), "logit" for ordbeta(), and "identity" for gaussian_ls(). Each family compiles the links for which the additive predictor is the natural unconstrained scale; any other link is applied from R, for the families where that is well defined. See Details.
theta: numeric; for negbin() and zi_negbin(), a fixed value for the dispersion parameter, which must be positive. Default is NULL to draw it along with everything else.
reference: for multinomial(), the response category to hold as the reference, given as a single value naming one of the response's levels. Default is NULL, which with the logit link fits one forest per category instead and leaves the model unidentified, which is what makes the prior symmetric in the categories; see Details. The probit link is always written as contrasts against a reference, so there the default is the first level.
replicates: numeric; for multinomial("probit"), how many simulation draws to use for the category probabilities, which have no closed form. Default is 200. More is always better but resulting calculations will take longer.
nu, q: numeric; for dpm() and dpm_aft(), the degrees of freedom of the baseline's inverse-chi-square prior on a component's variance and the quantile of that prior placed at a rough estimate of the residual standard deviation. Defaults are 10 and .95, following George et al. (2019), and are tighter than BART's own 3 and .90 because the mixture covers small errors with extra components rather than with one component's left tail.
k_s: numeric; for dpm() and dpm_aft(), how many units of the baseline's own scale the component means are allowed to reach out to. Must be positive. Default is 10, which places the marginal of a component mean so that it reaches the largest residual of a linear fit.
alpha: numeric; for dpm() and dpm_aft(), a fixed Dirichlet process concentration, which must be positive. Default is NULL to draw it.
max_clusters, psi: numeric; for dpm() and dpm_aft(), the largest number of mixture components thought plausible (2 or greater) and the shape of the taper toward it, which together set the prior on alpha. Defaults are NULL to use a tenth of the sample size, and .5.
num_bins: Advanced. numeric; for ph(), how many pieces the baseline hazard has, with the edges at evenly spaced quantiles of the observed times. Must be 2 or greater. Default is NULL to use about the cube root of the sample size, which is the order the Freedman-Diaconis rule gives for a histogram. This should not need to be set: the estimates are flat in it over a sixty-fold range, and it is here for checking that rather than for tuning. See Details.
lambda_shape: numeric; for ph(), the shape of the gamma prior on each bin's baseline hazard, which must be positive. Its rate is drawn. Default is 1.
update_lambda: logical; for ph(), whether to draw the baseline hazards. Default is TRUE; FALSE holds them at their prior mean, which is for diagnosis rather than analysis.
phi: numeric; for Beta() and ordbeta(), a fixed value for the beta precision, and for tweedie() a fixed value for the dispersion. Must be positive. Default is NULL to draw it.
power: numeric; for tweedie(), the variance power $p$ in $\mathrm{Var}(y) = \phi\mu^p$, strictly between 1 and 2. Default is 1.5, and the value is held fixed there rather than drawn, because the power is weakly identified from data of the sizes this package is used on and a badly determined power drags the dispersion around with it. Pass NULL to draw it, which is worth doing only with a large sample and a real interest in the shape rather than the mean.
logdens: for custom_family(), the log density, given as a function of the response and the additive predictors, function(y, eta), where y is a numeric vector of length n and eta an n by num_predictors matrix, returning a numeric vector of length n. With nuisance parameters it takes a third argument, function(y, eta, aux), where aux is a numeric vector of their current values. It is the log density of one unit of prior weight, so that weights behave as they do elsewhere, and terms free of eta may be dropped.
num_predictors: numeric; for custom_family(), how many additive predictors the density has (i.e., how many forests to fit). Default is 1.
start: numeric; for custom_family(), the value each additive predictor starts at, in place of the intercept-only fit the compiled families use. One value or one per predictor. Default is 0.
derivatives: for custom_family(), optional; a function(y, eta, h) returning a list with elements score and info, the first derivative of logdens with respect to the hth predictor and minus its second derivative, each a vector of length n. Default is NULL to take central differences of logdens. It covers the additive predictors only: a nuisance parameter is always differenced, which costs three calls per sweep rather than three per leaf.
aux_names: optional character; for custom_family(), the names of the nuisance parameters to draw, if any. Naming them is what declares them, because the names label the columns of fit$aux and are what summary() and diagnose() report them under. They must be distinct and non-empty. Default is NULL for none, unless aux_start is given, in which case the parameters are named by names(aux_start) when it carries names and positionally ("aux1", "aux2", and so on) when it does not.
aux_start: optional numeric; for custom_family(), the value each nuisance parameter starts at, given as one value or one per parameter. Default is NULL, which is 0 for each. The sampler will walk to the posterior from a poor start, so this need only be the right order of magnitude. Supplying it is a second way to declare the parameters, so aux_start = c(shape = 1) both names one and starts it at 1.
name: string; for custom_family(), a label used when printing the fit. Default is "custom".

Value

A <bartisan_family> object, which is a list containing at least the elements family and link and which inherits from family, so that bartisan() recognizes it wherever it recognizes an ordinary stats::family object.

Details

Every family reduces to a scalar additive predictor, or to several of them, together with the first two derivatives of the log density with respect to each. That is the whole interface the sampler needs, which is why the set of available families is not restricted to the conditionally conjugate ones.

vignette("families", package = "bartisan") covers all of this at length: how to choose a family, how to choose among the links a family offers, and what each family is and is not for. What follows is the short version, and the points that could lead to output being misread.

The supported families and links are:

Family	Links	Additive predictors	Drawn nuisance parameters
`gaussian()`	`identity`	1	residual standard deviation
`binomial()`	`logit`, `probit`, `cloglog`	1	none
`poisson()`	`log`	1	none
`negbin()`	`log`	1	dispersion
`Gamma("log")`	`log`; any other link is ignored	1	shape
`ordinal()`	`logit`, `probit`, `cloglog`	1	cutpoints
`multinomial()`	`logit`, `probit`	one per category, or per non-reference level	latent covariance, for the probit link
`weibull_aft()`, `loglogistic_aft()`, `lognormal_aft()`	none	1	scale
`ph()`	none	1	baseline hazard per bin
`dpm_aft()`	none	1	error mixture, concentration
`gaussian_ls()`	`identity`	2	none
`Gamma_ls()`	`log`	2	none
`zi_poisson()`	`log`	2	none
`zi_negbin()`	`log`	2	dispersion
`Beta()`	`logit`, `probit`, `cloglog`	1	precision
`ordbeta()`	`logit`	1	2 cutpoints, precision
`tweedie()`	`log`	1	dispersion, and the power if it is drawn
`dpm()`	`identity`	1	error mixture, concentration

A family with more than one additive predictor fits one forest per predictor. Nuisance parameters are drawn alongside the trees and reported in fit$aux.

Links the Engine Does Not Compile

The links listed above are the ones the sampler evaluates in compiled code. Any other link is accepted for gaussian(), binomial(), poisson() and Beta(), and applied from R by composing the caller's inverse link with the family's own, with the chain rule carrying the derivatives back. So binomial("cauchit") works, as does any link object of the kind stats::make.link() returns. It costs a call into R for every leaf the sampler visits, and the leaf prior scale is calibrated for the compiled link.

Two single-predictor families are exceptions. negbin() takes "log" alone, so a link given to it is an error rather than a composition. stats::Gamma() accepts any link and fits none of them but "log": every other link is dropped with a message, because the ones base R offers have inverses that go non-positive over part of the line and the additive predictor is unconstrained. custom_family() is the route to a gamma response on another link.

A composed link whose inverse has a restricted range (poisson("identity"), poisson("sqrt")) will give non-finite densities for some predictors. Those proposals are rejected rather than breaking the chain, but they are wasted work and the fit is worse for it, so bartisan() says so when it starts. Prefer links whose inverse is defined on the whole line.

The families with more than one additive predictor, or whose link enters somewhere other than a single mean (ordinal(), multinomial(), the accelerated failure time families, gaussian_ls(), the zero-inflated families and ordbeta()), take only their listed links. custom_family() is the way to reach anything else.

If No Family Is Given

family may be omitted, in which case it is read off the response's type and the choice reported with a message; bartisan() tabulates the lookup, since family is its argument. Two things about it are worth knowing here, because both are deliberate rather than oversights: a count is not given poisson() and a numeric response taking two values other than 0 and 1 is not given binomial(), either being a modeling decision rather than a reading of the response's type. dpm() cannot take prior weights, so a weighted fit with no family named is an error rather than a silent substitution.

What to Know Before Reading the Output

ordinal() accepts a numeric response as well as an ordered factor, taking its sorted unique values as the categories, and that is a method rather than a fallback: the cutpoints absorb the marginal distribution of the response and the forest explains only the ordering, so nothing is assumed about the error distribution and the model for $P(Y \le y \mid x)$ is invariant to any monotone transformation of the response. Predict with type = "mean". Bin the response onto twenty-odd quantiles first: one cutpoint per distinct value costs 73 seconds against 2.7 for gaussian() at 1000 observations, and twenty-five bins was both sixteen times faster and slightly more accurate. See the vignette.

ordinal() uses the cumulative-link parameterization of MASS::polr() , in which $P(Y \le k) = F(c_k - \eta)$, so larger values of the additive predictor shift mass toward higher categories. Only the differences $c_k - \eta_i$ are identified, so one location has to be pinned: with three or more categories the draws are reported in the chart where the additive predictor has mean zero over the fitted sample and every cutpoint is free, which is the chart polr() reports in when its predictors are centered. With exactly two categories the single boundary is folded into the intercept instead, so a two-category response is exactly binary regression with the matching link and on the same scale. cut1 is therefore a free parameter rather than a constant zero, which is a change from earlier versions.

multinomial() by default fits one forest per category and leaves the model unidentified, since adding any function of the predictors to every category's forest leaves the probabilities alone. This is the parameterization of Murray (2021), whose point is that the prior is then symmetric in the categories; every identified quantity is recovered from the draws. Passing reference instead pins that category at zero and fits one fewer forest, giving log odds against it.

multinomial("probit") lets the latent utilities correlate, which a multinomial logit cannot express at all. $\Sigma$ is normalized by the trace constraint $\mathrm{tr}(\Sigma) = C$ (Burgette and Nordheim 2012), and its lower triangle appears in fit$aux as sigma11, sigma21 and so on. Those correlations are weakly identified: read the fitted probabilities, not the covariance. They enter the likelihood only through orthant probabilities of a distribution whose location is a sum of trees, so a flexible mean absorbs much of the dependence they are meant to measure. At 900 observations a true correlation of zero came back as -0.57, and posterior intervals ran up to 1.07 wide on a parameter confined to $(-1, 1)$; at 3000 observations the posterior tracks the truth to within about 0.2. Two further consequences: the likelihood has no closed form, so it and every category probability are simulated with replicates draws, and augment does not apply, because the latent variables are the model rather than a rewriting of it. The sampler is Algorithm P2 of Xu et al. (2025).

dpm() is not a distribution but a way of not choosing one. It is DPMBART (George et al. 2019): a numeric response with the sum of trees for its mean, as gaussian() has, and a Dirichlet process mixture of normals for its errors instead of a single normal, so the error distribution comes out as whatever mixture the data ask for. error_density() gives that density, which is the object the method exists to produce.

It is the family to reach for by default on a numeric response, because it does not pay for its flexibility: on normal errors, where gaussian() is exactly right, it came out slightly ahead on both held-out error and log score at the same time to one decimal place, and on heavy-tailed, skewed and bimodal errors it was ahead by a great deal (on bimodal errors at a thousand observations, 0.050 against 0.154 in held-out RMSE, a factor of three, at the same time to a tenth of a second). So it is the family a numeric response gets when none is named. The reasons to prefer gaussian() are not statistical: it takes prior weights, which dpm() refuses, and it reports one interpretable sigma where dpm() has a mixture. It is also faster, by 1.4 times at a thousand observations. The vignette has the comparison.

Two things to know about dpm() itself. It does not buy heteroskedasticity: the error distribution is flexible but it is the same distribution at every $x$, and gaussian_ls() is the family for a spread that depends on the predictors. And the additive predictor is the conditional mean, as it is for gaussian(): nothing in the model forces the mixture to be centered, so the sampler works in a chart where only the sum of the predictor and the error mean is identified, but reporting is done in the chart where the mixture has mean zero and the whole conditional mean sits on the predictor. type = "link" and type = "response" therefore agree exactly, and fit$aux reports the shift that was taken out as center rather than an error mean, which is zero by construction. Prior weights are refused, since a weight would have to be a multiplicity in the Dirichlet process.

The gamma family puts the forest on the log mean and draws the shape, which acts as the inverse dispersion; it does not regress the shape on the predictors. negbin() and ordbeta() take theta and phi to fix their equivalents, and the gamma shape has no such argument because a caller who knows it is rare. The link is where the care is needed. Only log is compiled, and the base R default of inverse is the worst case for this sampler: its inverse maps a negative predictor to a negative mean, whose log is not a number, so the proposal is rejected, dozens of times per fit. Measured on 600 observations and 50 trees, stats::Gamma() took 7.2 seconds against 3.8 for Gamma("log") and fitted the mean slightly worse. So write the link, or name the family as the string "Gamma", which resolves to the log link; base R's own function is left as base R defines it, so that attaching bartisan cannot change what glm() does. Any composed link whose inverse has a restricted range is reported when the fit starts.

The accelerated failure time families expect a right-censored response, supplied either as a survival::Surv() object or as a two-column matrix of times and event indicators. They model $\log T = \eta + \sigma\epsilon$ with $\epsilon$ a standard smallest extreme value, logistic or normal variate respectively, giving Weibull, log-logistic and log-normal survival times.

A contrast in the predictor is a log time ratio in all of them, and in dpm_aft() too, because with $\epsilon$ independent of $x$ every quantile and both means of $T$ scale by $e^{\Delta\eta}$ whatever shape the error has. What differs is what $e^{\eta}$ is on its own, each family pinning its error's location differently: the median of $T$ for loglogistic_aft() and lognormal_aft(), the geometric mean for dpm_aft(), and the Weibull scale for weibull_aft(). Contrasts are unaffected by any of that; vignette("survival") tabulates the levels and measures the difference between them.

weibull_aft() is also the one family whose predictor carries a log hazard ratio, of $-\Delta\eta/\sigma$, alongside its log time ratio. That is a property of the smallest extreme value error rather than of the structural part it shares with the other three: it is the only error making an accelerated failure time model a proportional hazards model as well.

ph() is the proportional hazards alternative, with a piecewise-constant baseline: $\lambda(t \mid x) = \lambda_0(t)\exp(r(x))$, so its predictor is a log hazard ratio and the baseline is free to take any shape rather than the monotone one a Weibull imposes. num_bins sets how many pieces, with the edges at evenly spaced quantiles of the observed times; the default is about the cube root of the sample size. The bin hazards are drawn from their exact gamma conditionals and reported as lambda1, lambda2, ... in fit$aux, together with the rate of their own prior. The predictor and the baseline are identified only jointly, so the baseline carries the level and the predictor is reported centered on it.

Cox's partial likelihood is what cannot be used here: it couples observations through risk sets and so does not decompose into a sum over the observations reaching a leaf. The full likelihood of the piecewise-exponential model does decompose, and it approaches the partial likelihood as the bins shrink, which is how ph() reaches proportional hazards without it.

num_bins is not a modeling decision, and its default should be left alone. It is exposed for checking that rather than for tuning: swept over a sixty-fold range the estimates move by less than the spread between replicates, with no trend. What the bin count does change is the effective number of parameters, which grows with it, and that is what makes the default matter for loo() and waic() rather than for the estimates: one parameter per event time would leave each observation's density inflated by a parameter only it informs, and leave-one-out unable to do its job. vignette("survival") has the sweep.

The three differ in cost, though not enough to decide a model on. lognormal_aft() and loglogistic_aft() impute each censored failure time above its censoring time, which makes their targets quadratic and is worth a large multiple of the speed; weibull_aft() needs no imputation because its likelihood already has a form the sampler can collapse to a single pass, but only under hard rules, which makes it the slowest of the three at the default gate.

dpm_aft() is the accelerated failure time model with the error distribution estimated rather than assumed: $\log T = m(x) + W$ with $W$ a Dirichlet process mixture of normals constrained to mean zero, and censored log-times imputed. It is dpm()'s error model with censoring, so its predictor is the conditional mean of $\log T$ and a contrast in it is a log time ratio exactly as for the three parametric families above, error_density() reports the fitted error density, and prior weights are refused for the same reason dpm() refuses them. Following Henderson, Louis, Rosner and Varadhan (2020).

Reach for it when the shape of the error is in doubt and there is no reason to assert one. Measured against a two-component error it was worth 210 held-out log points and a third of the error in $S(t \mid x)$ over the best fixed-error family; against a log-normal error, where lognormal_aft() is correctly specified, the two were within 0.1 log points of each other. So it gains where the assumption would have been wrong and costs nothing where it would have been right, which is the property dpm() has against gaussian().

gaussian_ls() regresses the mean and the log standard deviation of a normal response on separate forests, so the variance is an unrestricted function of the predictors.

Gamma_ls() does the same for a gamma response: the first forest is the log mean, exactly as Gamma("log")'s is, and the second is the log dispersion, so the shape is exp(-log_dispersion) at each observation rather than one value drawn for the whole sample. What it relaxes is the assumption that the coefficient of variation is constant, which is what a gamma with a single shape asserts, and giving its second forest an intercept-only formula puts that assumption back; see "Several additive predictors" above.

Note that the mean forest is quadratic in neither family, but it takes the same cheap exponential form Gamma("log")'s does, while the dispersion forest takes the general path, so the cost sits in the second forest and it is worth giving that one fewer trees.

zi_poisson() and zi_negbin() are zero-inflated counts. Both parts get their own forest: the first predictor is the log mean of the count component and the second the log odds that an observation is a structural zero, so the excess-zero mechanism is free to depend on the predictors. The two are reported as the count and zero predictors.

Beta() is beta regression for a response strictly inside the unit interval: a forest on the link of the mean, and a precision drawn alongside it. A response at either endpoint has no beta density, so it is an error rather than something to nudge inward.

ordbeta() is the ordered beta regression of Kubinec (2023), for a response on the closed unit interval with point masses at zero and one. One predictor drives both the probability of landing on an endpoint, through a pair of cutpoints as in an ordinal model, and the mean of the beta density in between. Because the predictor also enters the beta mean it is identified, so unlike ordinal() both cutpoints are drawn.

Choose between them on whether the response can reach a boundary, not on whether it happens to in the sample: the two ask different questions, and ordbeta() fitted to a response with no boundary observations leaves its cutpoints with nothing to identify them.

tweedie() is the compound Poisson-gamma, for a non-negative response with a point mass at zero and a continuous positive part, which is the shape of spending, rainfall, insurance claims and earnings. It is the analogue of ordbeta() at the other end: one predictor again drives both parts, but through the mean rather than through a cutpoint, since $\mu = \exp(\eta)$ and $\mathrm{Var}(y) = \phi\mu^p$ together fix the probability of a zero at $\exp(-\mu^{2-p}/(\phi(2-p)))$. That is what makes it a single process and is also its restriction: the share of zeros has no level of its own, so a response whose zeros are more or less common than its mean implies wants a two-part model instead, which zi_poisson() and zi_negbin() are for counts and which custom_family() can supply for anything else.

Two things follow from the mean being $\exp(\eta)$ exactly. A counterfactual mean through marginaleffects needs nothing beyond the forest, unlike a two-part model where it has to be recombined across predictors; and the fit is comparable with a poisson() or Gamma("log") fit of the same response, since all three put the same quantity on the same scale.

Several Additive Predictors

Most families model one parameter with one forest. Some model several, and then every argument that could mean something different for each of them may be given per forest, keyed by the names below or positionally; bartisan_control() states the recycling rule and lists which arguments it covers, and formula is among them, so a forest can have predictors of its own.

The first forest is always the main parameter, the one a single-forest family would have on its own. This table is the canonical list of the names, which vignette("families") reproduces:

Family	Forests, in order
`gaussian_ls()`	`mean`, `log_sd`
`Gamma_ls()`	`mean`, `log_dispersion`
`zi_poisson()`, `zi_negbin()`	`count`, `zero`
`custom_family(num_predictors = k)`	`eta1` ... `etak`
`multinomial("logit")`	one per level, or per non-reference level
`multinomial("probit")`	one per non-reference level, named for its contrast
everything else	`eta`

mean is the mean and log_sd is the logarithm of the standard deviation, which is the scale the forest works on. Gamma_ls()'s mean is the log mean, as Gamma("log")'s single forest is, and its log_dispersion is the logarithm of the dispersion, so that in both location-scale families a larger second predictor means more spread. count is the linear predictor of the count component and zero that of the inflation component. A custom family's nuisance parameters are not on this list: they are carried as trailing forests pinned to a single leaf, and nothing about them is set per forest.

A forest whose formula names no predictor is a constant. ~ 1 leaves that forest nothing to split on, so every tree in it is a stump and the parameter is one drawn scalar. Every family here that takes more than one formula accepts that, which is what makes the distinction between a nuisance parameter and an empty forest a thin one: gaussian_ls() with ~ 1 on its scale is gaussian(), Gamma_ls() with ~ 1 is Gamma("log"), and zi_poisson() with ~ 1 on its inflation part is the zero-inflated Poisson with a single structural-zero probability. Note that the scalar is drawn under the leaf prior rather than under the prior the corresponding built-in family puts on its nuisance parameter, so the two agree closely rather than exactly. The multinomial families are the exception, for the reason given below.

So, for a location-scale model with a smaller scale forest and a restricted set of predictors for it:

bartisan(list(y ~ x1 + x2 + x3, log_sd = ~ x1), data = d,
         family = gaussian_ls(), num_trees = c(mean = 50, log_sd = 10))

The multinomial families are the exception. Their forests are the levels of one categorical parameter and act together rather than describing separate components of the response distribution, so there is nothing a caller could mean by giving one level a different prior or a different set of predictors from another. Every argument applies to all of their forests at once, and more than one value is an error rather than a silent recycling.

Supplying a Likelihood

custom_family() takes the log density itself, as an R function, and fits the model that goes with it. Nothing else about the sampler changes: the leaf-level Laplace proposal needs the first two derivatives of the log density with respect to each additive predictor and nothing more, and central differences of the supplied function produce both.

# A Poisson model written out by hand. Terms free of eta may be dropped;
# they cancel from every acceptance ratio.
pois <- custom_family(function(y, eta) y * eta[, 1] - exp(eta[, 1]),
                      start = log(mean(d$y)))

# Two predictors: a mean and a log standard deviation.
ls <- custom_family(function(y, eta) dnorm(y, eta[, 1], exp(eta[, 2]),
                                           log = TRUE),
                    num_predictors = 2, start = c(0, 0))

The function is called once per leaf per Fisher-scoring step with the observations reaching that leaf, so it must be vectorized over y and the rows of eta; it must not be vectorized within an observation, and it must return exactly one value per row. Supplying derivatives cuts three calls to one and removes the differencing error.

Nuisance parameters are drawn alongside the trees when aux_names names them, and logdens then takes a third argument holding their current values:

# A Gaussian written out by hand, with its scale drawn rather than fixed.
by_hand <- custom_family(
  logdens = function(y, eta, aux) dnorm(y, eta[, 1], exp(aux[1]), log = TRUE),
  aux_names = "log_sigma", aux_start = 0)

They are reported in fit$aux under those names, and covered by summary() and diagnose() like any other family's. There is no prior argument and no bounds argument, because a nuisance parameter here is carried as an additive predictor whose forest is pinned at depth zero (one tree that can never split, so the forest is a single scalar), and it is drawn by the same Laplace-plus-Metropolis step as any leaf, under that step's Gaussian leaf prior. So a parameter with a restricted range is handled the way it would be for a real predictor, by writing the transform into logdens: the exp() above is what keeps the scale positive.

What custom_family() does not do: the response must be numeric, so a factor has to be coded first; and since the package cannot know what the mean of the density is, predict(type = "response") returns the additive predictors rather than a fitted mean.

References

Burgette, L. F., & Nordheim, E. V. (2012). The trace restriction: an alternative identification strategy for the Bayesian multinomial probit model. Journal of Business & Economic Statistics, 30(3), 404–410. doi:10.1080/07350015.2012.680416

George, E., Laud, P., Logan, B., McCulloch, R., & Sparapani, R. (2019). Fully nonparametric Bayesian additive regression trees. In Topics in Identification, Limited Dependent Variables, Partial Observability, Experimentation, and Flexible Modeling: Part B (Advances in Econometrics, vol. 40B, pp. 89–110). Emerald Publishing. doi:10.1108/S0731-90532019000040B006

Henderson, N. C., Louis, T. A., Rosner, G. L., & Varadhan, R. (2020). Individualized treatment effects with censored data via fully nonparametric Bayesian accelerated failure time models. Biostatistics, 21(1), 50–68. doi:10.1093/biostatistics/kxy028

Kubinec, R. (2023). Ordered beta regression: a parsimonious, well-fitting model for continuous data with lower and upper bounds. Political Analysis, 31(4), 519–536. doi:10.1017/pan.2022.20

Murray, J. S. (2021). Log-linear Bayesian additive regression trees for multinomial logistic and count regression models. Journal of the American Statistical Association, 116(534), 756–769. doi:10.1080/01621459.2020.1813587

Xu, Y., Hogan, J., Daniels, M., Kantor, R., & Mwangi, A. (2025). Augmentation samplers for multinomial probit Bayesian additive regression trees. Journal of Computational and Graphical Statistics, 34(2), 498–508. doi:10.1080/10618600.2024.2388605

Examples

data("rhc")
set.seed(123)

# A right-censored response, given as the time and the event indicator,
# with the error distribution estimated rather than assumed
fit <- bartisan(cbind(days, death) ~ ., data = rhc, family = dpm_aft(),
                num_trees = 10, num_burn = 50, num_draws = 50)

# The shape the errors came out
head(error_density(fit))
#>          at         mean        lower        upper
#> 1 -7.951284 2.231688e-06 1.035634e-08 1.725141e-05
#> 2 -7.871771 2.740017e-06 1.400678e-08 2.073086e-05
#> 3 -7.792258 3.363112e-06 1.895640e-08 2.489060e-05
#> 4 -7.712745 4.126868e-06 2.629499e-08 2.985934e-05
#> 5 -7.633233 5.063076e-06 3.642054e-08 3.578917e-05
#> 6 -7.553720 6.210788e-06 4.584312e-08 4.285956e-05

# The same response under proportional hazards, whose predictor is a log
# hazard ratio and whose baseline is free to take any shape
bartisan(cbind(days, death) ~ ., data = rhc, family = ph(),
         num_trees = 10, num_burn = 50, num_draws = 50)
#> Generalized BART
#> 
#> Call:
#> bartisan(formula = cbind(days, death) ~ ., data = rhc, family = ph(), 
#>     num_trees = 10, num_burn = 50, num_draws = 50)
#> 
#> Family: "ph" with the "log" link
#> Observations: 1500
#> Structure: 1 forest of 10 trees, soft decision rules
#> Draws: 50 kept after 50 warmup
#> 
#> Posterior means: lambda1 = 0.0165, lambda2 = 0.0246, lambda3 = 0.0164, lambda4 = 0.0109, lambda5 = 0.00464, lambda6 = 0.00203, lambda7 = 0.00156, lambda8 = 0.00103, lambda9 = 0.00207, lambda10 = 0.00309, lambda11 = 0.00247, lambda12 = 0.00285, lambda_rate = 146

# An unordered response, with one forest per category and a prior that is
# symmetric in them
bartisan(race ~ . - days - death, data = rhc, family = multinomial(),
         num_trees = 10, num_burn = 50, num_draws = 50)
#> Generalized BART
#> 
#> Call:
#> bartisan(formula = race ~ . - days - death, data = rhc, family = multinomial(), 
#>     num_trees = 10, num_burn = 50, num_draws = 50)
#> 
#> Family: "multinomial" with the "logit" link
#> Observations: 1500
#> Structure: 3 forests of 10 trees, soft decision rules
#> Draws: 50 kept after 50 warmup

# A link the engine does not compile, applied from R
bartisan(death ~ . - days, data = rhc, family = binomial("cauchit"),
         num_trees = 10, num_burn = 50, num_draws = 50)
#> Generalized BART
#> 
#> Call:
#> bartisan(formula = death ~ . - days, data = rhc, family = binomial("cauchit"), 
#>     num_trees = 10, num_burn = 50, num_draws = 50)
#> 
#> Family: "binomial" with the "cauchit" link (supplied from R)
#> Observations: 1500
#> Structure: 1 forest of 10 trees, soft decision rules
#> Draws: 50 kept after 50 warmup