Multivariate model estimation with brms (2023)

Introduction

In this vignette, we want to discuss how to specify multilevel multivariate models usingbrms. we call it a modelmultivariateif it contains multiple response variables, each predicted by its own set of predictors. Consider an example from biology. Hadfield, Nutall, Osorio, and Owens (2007) analyzed data from the Eurasian blue tit (https://en.wikipedia.org/wiki/Eurasian_blue_tit). They predicted ithocklength, as well as theto returncolor of chicks. Half of the litter was placed in anothernest of adoption, while the other half stayed in the nest of their own childrenbarrier. This makes it possible to separate genetic from environmental factors. In addition, we have information aboutDate of birthmsexof offspring (the latter are known in 94% of animals).

data("BT Data",package = "MCMCglmm")head(BT data)

hind tarsus animal mother foster child hatch sex1 -1.89229718 1.1464212 R187142 R187557 F2102 -0.6874021 female2 1.13610981 -0.7596521 R1871504 R187154 0 .9 846894 6 0.1449373 R187341 R187568 A602 -0.4279814 Male4 0.37900806 0.2555847 R046169 R187518 A11465 -1 07525299 -0.3006992 R04 6161 R187528 A2602 -1.4656641 Female6 -1.13519543 1.5577219 R187409 R187945 C2302805 Female 0.3

Basic multivariate models

We start with a relatively simple multivariate normal model.

setting 1<- brm( mvbind(tarsus, Kostas)~sex+Date of birth+(1|Pi|adoption nest)+(1|q|barrier), data =BT data,chains = 2,cores = 2)

As can be seen from the template code, we usemvbindnotation for measurementbrmsthat bothhockmto returnare separate response variables. the term(1|p|first child)indicates an interrupt variablenest of adoption. Written|p|in the middle, we indicate that all variable effects ofnest of adoptionthey must be modeled as correlated. This makes sense, as we actually have two parts of the model, one forhockand one forto return. the indexPiis arbitrary and can be replaced by other symbols that come to mind (for details on the multi-level syntax ofbrms, edhelp ("brmsformula")mβινιέτα ("brms_multililevel")). Likewise, the term(1|q|dam)indicates correlated maternal genetic variable effects of neonates. Alternatively, we could also have modeled genetic similarities through pedigrees and matching paternity matrices, but this is not the focus of this vignette (seeβινιέτα ("brms_phylogenetics")). The results of the model are easily summarized

setting 1<- add_criterion(adjustment 1,"toilet")summary(adaptation 1)

Family: MV (Gaussian, Gaussian) Links: mu = identity; sigma = identity mu = identity; sigma = identity Type: tarsus ~ sex + date of birth + (1 | p | foster nest) + (1 | q | mother) back ~ sex + date of birth + (1 | p | foster nest) + (1 | q | mother ) Data: BTdata (Number of observations: 828) Plots: 2 strings, each with iter = 2000; heating = 1000; skinny = 1; Total Bonds After Heat = 2000 Group Level Effects: ~barrier (Number of Levels: 106) Rating Rating. Error l-95% CI u-95% CI Rhat Bulk_ESSsd(tarsus_Intercept) 0.48 0.05 0.40 0.59 1.2000 back . 0.10 0.39 1.00 364cor(tarsus_Intercept,back_Intercept) -0.52 0.22 -0.93 -0.08 1.00 425 Tail_ESSsd(tarsus_Intercept) 1675sd(back_Intercept,back_Intercept) (back_Intercept) of levels: 104 ) EstimateErrorEst l-95% CI u- 95% CI Rhat Bulk_ESSsd(tarsus_Intercept) 0.27 0.05 0 .17 0.38 1.00 885sd (back_Intercept) 0.35 0.06 0.24 0.47 1 ... CI u-95% CI Rhat Bulk_ESS Tail_ESStarsus_Intercept -0.41 0.07 -0.54 -0.28 1.00 1692 1592back_Intercept .0201 .701 . 10 13 85tarsus_sexMale 0.77 0.05 0.66 0.87 1.00 4820 1704tarsus_sexUNK 0.23 0.13 -0 .02 0.49 1.00 3290 1493tarsus_hatchdate -0.04 0.06 -0.15 0 .07 1.00 1684 1564back_sexMale 0.01 0.07 -0.131306013060. .1 4 0.15 -0.16 0.43 1.00 3837 1586back_hatchdate -0.09 0.05 -0.19 0.01 1.00 2244 1725Family Parameters: Estimate Est. Error l- 95% CI u- 95% CI Rhat Bulk_ESS Tail_ESSsigma_tarsus 0.76 0.02 0.72 0.80 1.00 2985 1455sigma_back 0.90 0.02 0.86 0.015031 .Error l-95%CI u-95%CI Rhat Bulk_ESS Tail_ESSrescor(tarsal ,back) - 0.05 0.04 -0.13 0.02 1.00 3757 1184 Specimens were sampled (NUTS). For each parameter, Bulk_ESS and Tail_ESS are effective measures of sample size, and Rhat is the potential downscaling factor in segregation chains (at convergence, Rhat = 1).

The summary output of multivariate models is very similar to univariate models, except that the parameters are now prefixed with the corresponding response variable. In mothers, tarsus length and back color appear to be negatively correlated, while in foster nestlings the opposite occurs. This indicates differential effects of genetic and environmental factors on these two traits. Moreover, the small residual correlationrescor (tarsus, back)at the bottom of the output suggests that there is little unmodeled dependence between tarsus length and back color. Although not necessary at this point, we have already calculated and stored its LOO information criterionsetting 1, which we will use for model comparisons. Next, let's take a look at some further predictive checks, which give us a first impression of the model's application.

pp_check(adjustment 1,resp = "hock")

pp_check(adjustment 1,resp = "to return")

This appears fairly stable, but we observe a slight unmodeled left skew in its distributionhock. We will come back to this later. Next, we want to investigate how much variation in the response variables can be explained by our model, and we use a Bayesian generalization\(R^2\)coefficient.

bayes_R2(adaptation 1)

Estimativa Est.Error Q2.5 Q97.5R2tarsus 0.4337432 0.02339576 0.3842017 0.4767824R2back 0.2002619 0.02780397 0.14042534

Clearly, there is a great deal of variation in the characteristics of both animals that we cannot explain, but apparently we can explain variation in tarsus length more than back color.

More complex multivariate models

Now suppose we just want to check forsexemhockbut not insideto returnand vice versa forDate of birth. Not that this makes much sense for the present example, but it allows us to show how to specify different types for different response variables. we can no longer usemvbindsyntax and therefore we need to use a more inclusive approach:

bf_tarsus<- friend(hock~sex+(1|Pi|adoption nest)+(1|q|barrier))bf_back<- friend(to return~Date of birth+(1|Pi|adoption nest)+(1|q|barrier))adapt 2<- brm(bf_tarsus+bf_back,data =BT data,chains = 2,cores = 2)

Note that we have literalsaddedthe two parts of the model through it+operator, which in this case is equivalent to writingmvbf(bf_tarsus, bf_back). Verhelp ("brmsformula")mhelp ("mvbrmsformula")for more details on this pension. Again, we summarize the model first.

adapt 2<- add_criterion(adjustment2,"toilet")summary(adaptation 2)

Family: MV (Gaussian, Gaussian) Links: mu = identity; sigma = identity mu = identity; sigma = identity Type: tarsus ~ sex + (1 | p | foster nest) + (1 | q | mother) back ~ date of birth + (1 | p | foster nest) + (1 | q | mother) Data: BTdata ( Number of observations: 828) Draws: 2 chains, each with iter = 2000; heating = 1000; skinny = 1; Total Bonds After Heat = 2000 Group Level Effects: ~barrier (number of levels: 106) Estimate Estimate. Error l-95% CI u-95% CI Rhat Bulk_ESSsd(tarsus_Intercept) 0.48 0.05 0.39 0.58 1.009 1.009 1.009 0.07 0.39 1.00 271cor(tarsus_Intercept,back_Intercept) -0.51 0.23 -0.93 -0.05 1.00 447 Tail_ESSsd(tarsus_Intercept) 1538sd(back_Intercept,back_Intercept,back_Intercepttar) of Levels: 104 ) ErrorEstimateEst l-95% CI u -95% CI Rhat Bulk_ESSsd(tarsus_Intercept) 0 .27 0.05 0 .16 0.37 1.00 833sd(back_Intercept) 0.35 0.06 0.23 0.46 1.00 563corp, 0.01 .01 0.01 .9 8 1.02 256 Tail_ESSsd( tarsus_Intercept) 1109sd(back_Intercept) 973cor(tarsus_Intercept, back_Intercept ) 605Population Level Effects: Estimate Estimate Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESStarsus_Intercept -0.41 0.07 -0.54 -0.28 1 .00 1866_01.01.000 back. 1.00 2281 15 00tarsus_sexMale 0.77 0.06 0.66 0.88 1.00 3657 1507tarsus_sexUNK 0.23 0.12 -0.01 0.47 1.00 3656 1588back_hatchdate -0.09 0.05 -0. ,19 0.02 1.00 1942 1747 urrors%lTimly 1747Es. CI Rhat Bulk_ESS Tail_ESSsigma_tarsus 0. 76 0.02 0. 72 0.79 1.00 3449 1647sigma_back 0.90 0.02 0.86 0.95 1.00 3053 1477 Residual Correlations : Estimate Error l-95% CI u -95% CI Rhat Bulk_ESS Tail_ESSrescor(tarsus,back) -0.05 0.03011-0.05 0.05 0.301 -00. sampling using sampling (NUTS). For each parameter, Bulk_ESS and Tail_ESS are effective measures of sample size, and Rhat is the potential downscaling factor in segregation chains (at convergence, Rhat = 1).

Let's find out how the model fit changed due to the exclusion of some effects from the original model:

toilet(slot1, slot2)

Output of model 'fit1': Calculated from 2000 with 828 log-likelihood matrix Estimate SEelpd_loo -2126.0 33.6p_loo 175.8 7.2looic 4252.1 67.1-------Montep Carlo diagnostareo SEelpd_loo . min. n_eff(-Inf, 0.5] (good) 812 98.1% 321 (0.5, 0.7] (ok) 13 1.6% 70 (0.7, 1] (bad) 3 0.4 % 138 (1, Inf) (very bad) 0 0.0%See help ('pareto-k-diagnostic') for details. Output of model 'fit2': Calculated from 2000 with 828 log probability matrix. counting. Min. n_eff(-Inf, 0.5] (good) 805 97.2% 330 (0.5, 0.7] (ok) 22 2.7% 90 (0.7, 1] (bad) 1 0.1 % 74 (1, Inf) (very bad) 0 0.0%See help ('pareto-k-diagnostic') for details. Model comparisons: elpd_diff se_difffit2 0.0 0.0 fit1 -0.2 1.2

Obviously, there is no noticeable difference in the application of the model. So we don't really need to modelsexmDate of birthfor both response variables, but it doesn't hurt to include them (so I probably would).

To give you a taste of its potentialbrmsmultivariate syntax, we change our model in many directions at once. Note its slight tilt to the lefthock, which we are now going to model using theskew_normalfamily instead ofgaussianfamily. As we no longer have a normal multivariate (or t-student) model, estimation of residual correlations is no longer possible. We make this explicit by usingset_rescormode. In addition, we investigated whether theto returnmDate of birthis indeed linear as we assumed earlier by fitting a non-linear spline to itDate of birth. In addition, we model separate residual variances ofhockfor male and female chicks.

bf_tarsus<- friend(hock~sex+(1|Pi|adoption nest)+(1|q|barrier))+ with(Sigma~ 0 +sex)+ skew_normal()bf_back<- friend(to return~ small(hatch date)+(1|Pi|adoption nest)+(1|q|barrier))+ gaussian()ajuste3<- brm(bf_tarsus+bf_back+ set_rescor(FALSE), data =BT data,chains = 2,cores = 2, control = list(adapt_delta = 0,95))

Again, we summarize the model and consider some further predictive checks.

ajuste3<- add_criterion(just3,"toilet")summary(just3)

Family: MV(skew_normal, gaussian) Links: mu = identity; sigma = write; alpha = identity mu = identity; sigma = identity Type: tarsus ~ sex + (1 | p | foster nest) + (1 | q | mother) sigma ~ 0 + sex back ~ s(birthdate) + (1 | p | foster nest) + (1 | q|mother) Data: BTdata (Number of observations: 828) Plots: 2 strings, each with iter = 2000; heating = 1000; skinny = 1; total bonds after heat = 2000 Smooth terms: Estimate Est. Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSsds(back_shatchdate_1) 2.08 1.09 0.35 4.72 1.00 411 ~ 400 Level level Estimate Estimate Error l-95% CI u-95% CI Rhat Bulk_ESSsd(tarsus_Intercept) 0.47 0.05 0.38 0.57 1.00 552sd(back_Intercept) 0 .24 0.07 0.07 0.10_Intercept) 0.11003 -0.52 0.2 ​​2 -0.92 -0.08 1.00 363 Tail_ESSsd(tarsus_Intercept ) 839sd(back_Intercept) 552cor(tarsus_Intercept,back_Intercept) 515~fosternest (Number of levels: 104) Error estimate l-95% CI u_Intercept. 0.15 0.37 1.00 408d (back_Intercept) 0.31 0.06 0.19 0.43 1.00 345cor(tarsus_Intercept,back_Intercept) 0.65 0.22 0.17 0.98 1.02 205 Tail_ESSsd(tarsus_Intercept,back_6) Intercept) 312 Population-Level Effects: Estimation Est. Error l-95% CI u- 95% CI Rhat Bulk_ESS Tail_ESStarsus_Intercept -0.42 0.07 -0.54 -0.28 1.00 675 849back_Intercept 0.00 0.05 -0.10 0.10 1.00 1228 0.54 -0.28 . 9 1.00 2565 1723tarsus_sexUNK 0.21 0.12 -0.03 0.44 1.00 1828 1347sigma_tarsus_sexFemale -0.30 0.04 -0.38 -0.21 1.23100 . 6.57 1.00 793 1197 Family parameters: Estimate Estimate Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSsigma_back 0.90 0.02 0.85 0.95 1.00 1818 1527alpha_tarsus -1.25 0. 41 -1.90 -0.24 1.00 1275 568 Plans were sampled using NUTS sampling). For each parameter, Bulk_ESS and Tail_ESS are effective measures of sample size, and Rhat is the potential downscaling factor in segregation chains (at convergence, Rhat = 1).

We see that the residual (log) standard deviation ofhockit is slightly higher for chicks whose sex cannot be determined compared to male or female chicks. In addition, we see from the negativesalpha(deformation) parameter ofhockthat the residuals are actually slightly skewed to the left. Finally, running

conditional_effect(just3,"hatch date",resp = "to return")

reveals a non-linear relationship ofDate of birthnoto returncolor, which seems to change in waves during birth dates.

There are many more modeling options for multivariate models that are not discussed in this vignette. Examples include autocorrelation structures, Gaussian processes, or explicit non-linear predictions (for example, seehelp ("brmsformula")theβινιέτα ("brms_multililevel")). In fact, almost all the flexibility of univariate models is retained in multivariate models.

bibliographical references

Hadfield JD, Nutall A, Osório D, Owens IPF (2007). Testing the phenotypic gambit: phenotypic, genetic, and environmental correlates of color.Journal of Evolutionary Biology, 20(2), 549-557.