brmsformula
Set up a model formula for use in the brms package
Set up a model formula for use in the brms package allowing to define (potentially nonlinear) additive multilevel models for all parameters of the assumed response distribution.
Usage
brmsformula(formula, ..., nonlinear = NULL)
Arguments
 formula
 An object of class
formula
(or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under 'Details'.  ...
 Additional
formula
objects to specify predictors of special model parts and auxiliary parameters. Formulas can either be named directly or contain names on their lefthand side. Currently, the following names are accepted:sigma
(residual standard deviation of thegaussian
andstudent
families);shape
(shape parameter of theGamma
,weibull
,negbinomial
and related zeroinflated / hurdle families);nu
(degrees of freedom parameter of thestudent
family);phi
(precision parameter of thebeta
andzero_inflated_beta
families);zi
(zeroinflation probability);hu
(hurdle probability). All auxiliary parameters are modeled on the log or logit scale to ensure correct definition intervals after transformation.  nonlinear
 An optional list of formuluas, specifying
linear models for nonlinear parameters. If
NULL
(the default)formula
is treated as an ordinary formula. If notNULL
,formula
is treated as a nonlinear model andnonlinear
should contain a formula for each nonlinear parameter, which has the parameter on the left hand side and its linear predictor on the right hand side. Alternatively, it can be a single formula with all nonlinear parameters on the left hand side (separated by a+
) and a common linear predictor on the right hand side. More information is given under 'Details'.
Details
The formula
argument accepts formulae of the following syntax:
response  addition ~ Pterms + (Gterms  group)
The Pterms
part contains effects that are assumed to be the
same across obervations. We call them 'populationlevel' effects
or (adopting frequentist vocabulary) 'fixed' effects. The optional
Gterms
part may contain effects that are assumed to vary
accross grouping variables specified in group
. We
call them 'grouplevel' effects or (adopting frequentist
vocabulary) 'random' effects, although the latter name is misleading
in a Bayesian context (for more details type vignette("brms")
).
Multiple grouping factors each with multiple grouplevel effects
are possible. Instead of 
you may use 
in grouping terms
to prevent correlations from being modeled.
Alternatively, it is possible to model different grouplevel terms of
the same grouping factor as correlated (even across different formulae,
e.g. in nonlinear models) by using 
instead of 
.
All grouplevel terms sharing the same ID will be modeled as correlated.
If, for instance, one specifies the terms (1+x2g)
and
(1+z2g)
somewhere in the formulae passed to brmsformula
,
correlations between the corresponding grouplevel effects
will be estimated.
Smoothing terms can modeled using the s
and t2
functions of the mgcv package
in the Pterms
part of the model formula.
This allows to fit generalized additive mixed models (GAMMs) with brms.
The implementation is similar to that used in the gamm4 package.
For more details on this model class see gam
and gamm
.
The Pterms
part may contain two nonstandard types
of populationlevel effects namely monotonic and category specific effects,
which can be specified using terms of the form monotonic(
and cse(
respectively. The latter can only be applied in
ordinal models and is explained in more detail in the package's vignette
(type vignette("brms")
). The former effect type is explained here.
A monotonic predictor must either be integer valued or an ordered factor,
which is the first difference to an ordinary continuous predictor.
More importantly, predictor categories (or integers) are not assumend to be
equidistant with respect to their effect on the response variable.
Instead, the distance between adjacent predictor categories (or integers)
is estimated from the data and may vary across categories.
This is realized by parameterizing as follows:
One parameter takes care of the direction and size of the effect similar
to an ordinary regression parameter, while an additional parameter vector
estimates the normalized distances between consecutive predictor categories.
A main application of monotonic effects are ordinal predictors that
can this way be modeled without (falsely) treating them as continuous
or as unordered categorical predictors.
The third exception is the optional addition
term, which may contain
multiple terms of the form fun(variable)
seperated by +
each
providing special information on the response variable. fun
can be
replaced with either se
, weights
, disp
, trials
,
cat
, cens
, or trunc
. Their meanings are explained below.
For families gaussian
and student
, it is
possible to specify standard errors of the observation, thus allowing
to perform metaanalysis. Suppose that the variable yi
contains
the effect sizes from the studies and sei
the corresponding
standard errors. Then, fixed and random effects metaanalyses can
be conducted using the formulae yi  se(sei) ~ 1
and
yi  se(sei) ~ 1 + (1study)
, respectively, where
study
is a variable uniquely identifying every study.
If desired, metaregression can be performed via
yi  se(sei) ~ 1 + mod1 + mod2 + (1study)
or yi  se(sei) ~ 1 + mod1 + mod2 + (1 + mod1 + mod2study)
, where
mod1
and mod2
represent moderator variables.
For all families, weighted regression may be performed using
weights
in the addition part. Internally, this is
implemented by multiplying the logposterior values of each
observation by their corresponding weights.
Suppose that variable wei
contains the weights
and that yi
is the response variable.
Then, formula yi  weights(wei) ~ predictors
implements a weighted regression.
The addition argument disp
(short for dispersion) serves a
similar purpose than weight
. However, it has a different
implementation and is less general as it is only usable for the
families gaussian
, student
,
lognormal
, Gamma
, weibull
, and negbinomial
.
For the former four families, the residual standard deviation
sigma
is multiplied by the values given in
disp
, so that higher values lead to lower weights.
Contrariwise, for the latter three families, the parameter shape
is multiplied by the values given in disp
. As shape
can be understood as a precision parameter (inverse of the variance),
higher values will lead to higher weights in this case.
For families binomial
and zero_inflated_binomial
,
addition should contain a variable indicating the number of trials
underlying each observation. In lme4
syntax, we may write for instance
cbind(success, n  success)
, which is equivalent
to success  trials(n)
in brms
syntax. If the number of trials
is constant across all observation (say 10
),
we may also write success  trials(10)
.
For all ordinal families, addition
may contain a term
cat(number)
to specify the number categories (e.g, cat(7)
).
If not given, the number of categories is calculated from the data.
With the expection of categorical
and ordinal families,
left and right censoring can be modeled through
yi  cens(censored) ~ predictors
.
The censoring variable (named censored
in this example) should
contain the values 'left'
, 'none'
, and 'right'
(or equivalenty 1, 0, and 1) to indicate that the corresponding observation is
left censored, not censored, or right censored.
With the expection of categorical
and ordinal families, the response
distribution can be truncated using the trunc
function in the addition part.
If the response variable is truncated between, say, 0 and 100, we can specify this via
yi  trunc(lb = 0, ub = 100) ~ predictors
.
Instead of numbers, variables in the data set can also be passed allowing
for varying truncation points across observations.
Defining only one of the two arguments in trunc
leads to onesided truncation.
Mutiple addition
terms may be specified at the same time using
the +
operator, for instance formula = yi  se(sei) + cens(censored) ~ 1
for a censored metaanalytic model.
For families gaussian
and student
,
multivariate models may be specified using cbind
notation.
In brms 1.0.0, the multvariate 'trait' syntax was removed
from the package as it repeatedly confused users, required much
special case coding, and was hard to maintain. Below the new
syntax is described.
Suppose that y1
and y2
are response variables
and x
is a predictor.
Then cbind(y1,y2) ~ x
specifies a multivariate model,
The effects of all terms specified at the RHS of the formula
are assumed to vary across response variables (this was not the
case by default in brms < 1.0.0). For instance, two parameters will
be estimated for x
, one for the effect
on y1
and another for the effect on y2
.
This is also true for grouplevel effects. When writing, for instance,
cbind(y1,y2) ~ x + (1+xg)
, grouplevel effects will be
estimated separately for each response. To model these effects
as correlated across responses, use the ID syntax (see above).
For the present example, this would look as follows:
cbind(y1,y2) ~ x + (1+x2g)
. Of course, you could also use
any value other than 2
as ID. It is not yet possible
to model terms as only affecting certain responses (and not others),
but this will comebe implemented in the future.
Categorical models use the same syntax as multivariate
models. As in most other implementations of categorical models,
values of one category (the first in brms) are fixed
to identify the model. Thus, all terms on the RHS of
the formula correspond to K  1
effects
(K
= number of categories), one for each nonfixed category.
Grouplevel effects may be specified as correlated across
categories using the ID syntax.
As of brms 1.0.0, zeroinflated and hurdle models are specfied
in the same way as as their noninflated counterparts.
However, they have additional auxiliary parameters
(named zi
and hu
respectively)
modeling the zeroinflation / hurdle probability depending on which
model you choose. These parameters can also be affected by predictors
in the same way the response variable itself. See the end of the
Details section for information on how to accomplish that.
Parameterization of the populationlevel intercept
The populationlevel intercept (if incorporated) is estimated separately
and not as part of populationlevel parameter vector b
.
also have to be specified separately
(see set_prior
for more details).
Furthermore, to increase sampling efficiency, the fixed effects
design matrix X
is centered around its column means
X_means
if the intercept is incorporated.
This leads to a temporary bias in the intercept equal to
, where <,>
is the scalar product.
The bias is corrected after fitting the model, but be aware
that you are effectively defining a prior on the temporary
intercept of the centered design matrix not on the real intercept.
This behavior can be avoided by using the reserved
(and internally generated) variable intercept
.
Instead of y ~ x
, you may write
y ~ 0 + intercept + x
. This way, priors can be
defined on the real intercept, directly. In addition,
the intercept is just treated as an ordinary fixed effect
and thus priors defined on b
will also apply to it.
Note that this parameterization may be a bit less efficient
than the default parameterization discussed above.
Formula syntax for nonlinear models
Using the nonlinear
argument, it is possible to specify
nonlinear models in brms. Contrary to what the name might suggest,
nonlinear
should not contain the nonlinear model itself
but rather information on the nonlinear parameters.
The nonlinear model will just be specified within the formula
argument. Suppose, that we want to predict the response y
through the predictor x
, where x
is linked to y
through y = alpha  beta * lambda^x
, with parameters
alpha
, beta
, and lambda
. This is certainly a
nonlinear model being defined via
formula = y ~ alpha  beta * lambda^x
(addition arguments
can be added in the same way as for ordinary formulas).
Now we have to tell brms the names of the nonlinear parameters
and specfiy a (linear mixed) model for each of them using the nonlinear
argument. Let's say we just want to estimate those three parameters
with no further covariates or random effects. Then we can write
nonlinear = alpha + beta + lambda ~ 1
or equivalently
(and more flexible) nonlinear = list(alpha ~ 1, beta ~ 1, lambda ~ 1)
.
This can, of course, be extended. If we have another predictor z
and
observations nested within the grouping factor g
, we may write for
instance nonlinear = list(alpha ~ 1, beta ~ 1 + z + (1g), lambda ~ 1)
.
The formula syntax described above applies here as well.
In this example, we are using z
and g
only for the
prediction of beta
, but we might also use them for the other
nonlinear parameters (provided that the resulting model is still
scientifically reasonable).
Nonlinear models may not be uniquely identified and / or show bad convergence.
For this reason it is mandatory to specify priors on the nonlinear parameters.
For instructions on how to do that, see set_prior
.
Formula syntax for predicting auxiliary parameters
It is also possible to predict auxiliary parameters of the response
distribution such as the residual standard deviation sigma
in gaussian models or the hurdle probability hu
in hurdle models.
The syntax closely resembles that of a nonlinear
parameter, for instance sigma ~ x + s(z) + (1+xg)
.
All auxiliary parameters currently supported by brmsformula
have to positive (a negative standard deviation or precision parameter
doesn't make any sense) or are bounded between 0 and 1 (for zeroinflated /
hurdle proabilities).
However, linear predictors can be positive or negative, and thus
the log link (for positive parameters) or logit link (for probability parameters)
are used to ensure that auxiliary parameters are within their valid intervals.
This implies that effects for auxiliary parameters are estimated on the
log / logit scale and one has to apply the inverse link function to get
to the effects on the original scale.
Value

An object of class
brmsformula
, which inherits
from class formula
but contains additional attributes.
Examples
# multilevel model with smoothing terms
brmsformula(y ~ x1*x2 + s(z) + (1+x11) + (1g2))
# additionally predict 'sigma'
brmsformula(y ~ x1*x2 + s(z) + (1+x11) + (1g2),
sigma ~ x1 + (1g2))
# use the shorter alias 'bf'
(formula1 < brmsformula(y ~ x + (xg)))
(formula2 < bf(y ~ x + (xg)))
# will be TRUE
identical(formula1, formula2)
# incorporate censoring
bf(y  cens(censor_variable) ~ predictors)
# define a nonlinear model
bf(y ~ a1  a2^x, nonlinear = list(a1 ~ 1, a2 ~ x + (xg)))
# correlated grouplevel effects across parameters
bf(y ~ a1  a2^x, nonlinear = list(a1 ~ 1 + (12g), a2 ~ x + (x2g)))
# define a multivariate model
bf(cbind(y1, y2) ~ x * z + (1g))
# define a zeroinflated model
# also predicting the zeroinflation part
bf(y ~ x * z + (1+xID1g), zi ~ x + (1ID1g))
# specify a predictor as monotonic
bf(y ~ mono(x) + more_predictors)
# specify a predictor as category specific
# for ordinal models only
bf(y ~ cse(x) + more_predictors)