Hierarchical Bayesian modeling of random and residual variance–covariance matrices in bivariate mixed effects models

Bivariate mixed effects models are often used to jointly infer upon covariance matrices for both random effects ( u ) and residuals ( e ) between two different phenotypes in order to investigate the architecture of their relationship. However, these (co)variances themselves may additionally depend upon covariates as well as additional sets of exchangeable random effects that facilitate borrowing of strength across a large number of clusters. We propose a hierarchical Bayesian extension of the classical bivariate mixed effects model by embedding additional levels of mixed effects modeling of reparameterizations of u‐ level and e ‐level (co)variances between two traits. These parameters are based upon a recently popularized square‐root‐free Cholesky decomposition and are readily interpretable, each conveniently facilitating a generalized linear model characterization. Using Markov Chain Monte Carlo methods, we validate our model based on a simulation study and apply it to a joint analysis of milk yield and calving interval phenotypes in Michigan dairy cows. This analysis indicates that the e ‐level relationship between the two traits is highly heterogeneous across herds and depends upon systematic herd management factors.