stopifnot(require(knitr)) options(width = 90) opts_chunk$set( comment = NA, message = FALSE, warning = FALSE, eval = if (isTRUE(exists("params"))) params$EVAL else FALSE, dev = "png", dpi = 150, fig.asp = 0.8, fig.width = 5, out.width = "60%", fig.align = "center" ) library(brms) ggplot2::theme_set(theme_default())
In the present vignette, we want to discuss how to specify multivariate multilevel models using brms. We call a model multivariate if it contains multiple response variables, each being predicted by its own set of predictors. Consider an example from biology. Hadfield, Nutall, Osorio, and Owens (2007) analyzed data of the Eurasian blue tit (https://en.wikipedia.org/wiki/Eurasian_blue_tit). They predicted the
tarsus length as well as the
back color of chicks. Half of the brood were put into another
fosternest, while the other half stayed in the fosternest of their own
dam. This allows to separate genetic from environmental factors. Additionally, we have information about the
sex of the chicks (the latter being known for 94\% of the animals).
data("BTdata", package = "MCMCglmm") head(BTdata)
We begin with a relatively simple multivariate normal model.
fit1 <- brm( mvbind(tarsus, back) ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam), data = BTdata, chains = 2, cores = 2 )
As can be seen in the model code, we have used
mvbind notation to tell
brms that both
back are separate response variables. The term
(1|p|fosternest) indicates a varying intercept over
fosternest. By writing
|p| in between we indicate that all varying effects of
fosternest should be
modeled as correlated. This makes sense since we actually have two model parts,
tarsus and one for
back. The indicator
p is arbitrary and can be
replaced by other symbols that comes into your mind (for details about the
multilevel syntax of brms, see
vignette("brms_multilevel")). Similarly, the term
correlated varying effects of the genetic mother of the chicks. Alternatively,
we could have also modeled the genetic similarities through pedigrees and
corresponding relatedness matrices, but this is not the focus of this vignette
vignette("brms_phylogenetics")). The model results are readily
fit1 <- add_criterion(fit1, "loo") summary(fit1)
The summary output of multivariate models closely resembles those of univariate
models, except that the parameters now have the corresponding response variable
as prefix. Within dams, tarsus length and back color seem to be negatively
correlated, while within fosternests the opposite is true. This indicates
differential effects of genetic and environmental factors on these two
characteristics. Further, the small residual correlation
on the bottom of the output indicates that there is little unmodeled dependency
between tarsus length and back color. Although not necessary at this point, we
have already computed and stored the LOO information criterion of
we will use for model comparisons. Next, let's take a look at some
posterior-predictive checks, which give us a first impression of the model fit.
pp_check(fit1, resp = "tarsus") pp_check(fit1, resp = "back")
This looks pretty solid, but we notice a slight unmodeled left skewness in the
tarsus. We will come back to this later on. Next, we want to
investigate how much variation in the response variables can be explained by our
model and we use a Bayesian generalization of the $R^2$ coefficient.
Clearly, there is much variation in both animal characteristics that we can not explain, but apparently we can explain more of the variation in tarsus length than in back color.
Now, suppose we only want to control for
tarsus but not in
vice versa for
hatchdate. Not that this is particular reasonable for the
present example, but it allows us to illustrate how to specify different
formulas for different response variables. We can no longer use
and so we have to use a more verbose approach:
bf_tarsus <- bf(tarsus ~ sex + (1|p|fosternest) + (1|q|dam)) bf_back <- bf(back ~ hatchdate + (1|p|fosternest) + (1|q|dam)) fit2 <- brm(bf_tarsus + bf_back, data = BTdata, chains = 2, cores = 2)
Note that we have literally added the two model parts via the
which is in this case equivalent to writing
mvbf(bf_tarsus, bf_back). See
help("mvbrmsformula") for more details about this
syntax. Again, we summarize the model first.
fit2 <- add_criterion(fit2, "loo") summary(fit2)
Let's find out, how model fit changed due to excluding certain effects from the initial model:
Apparently, there is no noteworthy difference in the model fit. Accordingly, we
do not really need to model
hatchdate for both response variables,
but there is also no harm in including them (so I would probably just include
To give you a glimpse of the capabilities of brms' multivariate syntax, we
change our model in various directions at the same time. Remember the slight
left skewness of
tarsus, which we will now model by using the
family instead of the
gaussian family. Since we do not have a multivariate
normal (or student-t) model, anymore, estimating residual correlations is no
longer possible. We make this explicit using the
set_rescor function. Further,
we investigate if the relationship of
hatchdate is really linear as
previously assumed by fitting a non-linear spline of
hatchdate. On top of it,
we model separate residual variances of
tarsus for male and female chicks.
bf_tarsus <- bf(tarsus ~ sex + (1|p|fosternest) + (1|q|dam)) + lf(sigma ~ 0 + sex) + skew_normal() bf_back <- bf(back ~ s(hatchdate) + (1|p|fosternest) + (1|q|dam)) + gaussian() fit3 <- brm( bf_tarsus + bf_back + set_rescor(FALSE), data = BTdata, chains = 2, cores = 2, control = list(adapt_delta = 0.95) )
Again, we summarize the model and look at some posterior-predictive checks.
fit3 <- add_criterion(fit3, "loo") summary(fit3)
We see that the (log) residual standard deviation of
tarsus is somewhat larger
for chicks whose sex could not be identified as compared to male or female
chicks. Further, we see from the negative
alpha (skewness) parameter of
tarsus that the residuals are indeed slightly left-skewed. Lastly, running
conditional_effects(fit3, "hatchdate", resp = "back")
reveals a non-linear relationship of
hatchdate on the
back color, which
seems to change in waves over the course of the hatch dates.
There are many more modeling options for multivariate models, which are not
discussed in this vignette. Examples include autocorrelation structures,
Gaussian processes, or explicit non-linear predictors (e.g., see
vignette("brms_multilevel")). In fact, nearly all the
flexibility of univariate models is retained in multivariate models.
Hadfield JD, Nutall A, Osorio D, Owens IPF (2007). Testing the phenotypic gambit: phenotypic, genetic and environmental correlations of colour. Journal of Evolutionary Biology, 20(2), 549-557.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.