Suppose you have data on $y_{it}$: outcome $y$ of unit $i$ at time $t$. You fit a varying intercept (hierarchical) model: $$y_{it} \sim Normal(\mu_{i}, \sigma_y)\\ \mu_i \sim Normal(0, \psi_{\mu})\\ \sigma_y \sim Cauchy(0,1)\\ \psi_{\mu} \sim Cauchy(0,1)\\ $$

I often see specifications where the random effects $\mu_i$ follow some multivariate distribution. The stated rationale for this is typically that the random effects are correlated in some way. For example, https://cran.r-project.org/web/packages/brms/vignettes/brms_multivariate.html talks about the need to correlate a bunch of random intercepts and proposes a multivariate distribution for this.

This puzzles me since typically the justification for the varying intercepts is that the outcomes $y_{it}$ might be correlated between time periods $t$ but are uncorrelated within units $i$. Accordingly, including the random effects in the model now means that $y$ are independent of each other. If we thought the $u_i$ were correlated with each other, then we could add another set of intercepts to model those higher level correlations. However, it's typically not the case!

So, why don't people just use univariate normal likelihoods/priors instead of multivariate?

  • $\begingroup$ Analogous to multilevel/mixed-effects models in a frequentist framework, random effects on the same level of measurement are often correlated, e.g., intercepts and time-slopes (co)varying across subjects. The assumption of uncorrelated errors would be violated if the random-effects model components are not correctly specified. $\endgroup$
    – Terrence
    Commented Oct 3, 2022 at 9:44
  • $\begingroup$ @Terrence 1. I'm not sure I understand -in the above model, the unit-specific intercepts $\mu_i$ are drawn from a distribution with a fixed mean (0) and common variance parameter. In this case, the correlation in the errors for each $it$ (unit time/level 1) is modeled by including the intercepts, so those errors are now independent. At the second level (the unit level $i$) we might have correlations -- units might not be independent. However, when we e.g. estimate a model with OLS we usually assume independent errors at the "top level." Wouldn't we do the same here? $\endgroup$ Commented Oct 3, 2022 at 19:56
  • $\begingroup$ In other words, isn't the marginal dependence broken by simply including the random intercept so that now we have conditional independence? Accordingly, since an identity matrix for the covariance of the random effects is the same as independent univariate normals, no issue...? $\endgroup$ Commented Oct 3, 2022 at 20:08
  • $\begingroup$ Conditioning on random intercepts will make errors independent if that is the only thing that varies across subjects in the population. But I think this actually does not answer your question. The link you provided is about multivariate models, not merely repeated measures of a single variable. I'll post an answer now that I see what the misunderstanding is. $\endgroup$
    – Terrence
    Commented Oct 4, 2022 at 11:48

1 Answer 1


Your intuition works for univariate multilevel/hierarchical models, but the link you provide describes multivariate multilevel models. So rather than

data on $𝑦_{𝑖𝑡}$: outcome $𝑦$ of unit $𝑖$ at time $𝑡$

the outcomes are a vector, say $[y_{it}, z_{it}]^{'}$ (variables tarsus and back in the posted example). Thus, there are 2 random intercepts (for $y$ and $z$) rather than 1 (for $y$ only), say:

  • $u_i$: the degree to which subject $i$'s average score on $y$ is higher/lower than the grand mean $\bar{y}$
  • $v_i$: the degree to which subject $i$'s average score on $z$ is higher/lower than the grand mean $\bar{z}$

Conditioning on these random intercepts can indeed make the errors $[e^{(y)}_{it}, e^{(z)}_{it}]^{'}$ conditionally independent across occasions (e.g., $t$ with $t-1$) within each subject $i$. However, that has nothing to do with correlations across variables at the within- or between-subject levels.

  • Correlation of $u_i$ with $v_i$ indicates whether subjects above average on $y$ tend to be above/below average on $z$
  • Correlation of $e^{(y)}_{it}$ with $e^{(z)}_{it}$ indicates whether subject $i$ being higher than usual on $y$ at Time $t$ also tends to be higher/lower than usual on $z$ at Time $t$

So, whereas a univariate MLM decomposes $y_{it}$ into subject-level ($\sigma^2_i$) and occasion-level ($\sigma^2_t$) variance components, a multivariate MLM decomposes the entire covariance matrix $\Sigma$ among multiple outcomes ($y$ and $z$) into a subject-level $\Sigma_i$ and an occasion-level $\Sigma_t$. This is the basis for multilevel SEM, if you are familiar with such models, e.g.:


So to clarify: across-variable correlations of random effects (at Level 2) or errors (at Level 1) are not a violation of conditional independence. Conditional on the random effects, the error-vectors for rows of data are still (conditionally) independent, although there can still be correlation across columns of data (i.e., correlation of residuals of variable $y$ with residuals of variable $z$).

  • $\begingroup$ Great, thanks. I'm just left somewhat confused why then here they discuss multivariate priors for hierarchical models mc-stan.org/docs/stan-users-guide/… but then the multilevel 2pl model shown here doesn't mc-stan.org/docs/stan-users-guide/item-response-models.html $\endgroup$ Commented Oct 4, 2022 at 14:17
  • $\begingroup$ Similarly, the way you describe it I'd expect for e.g. 3 outcomes there to be a 3x3 covariance matrix for the item intercepts, but I don't see this. Similarly, with 1 item intercept and 1 factor loading per outcome, should I expect 1 bivariate normal per item (since there are 2 parameter per item) or is is the dimensionality larger? , I'd expect the item-level discrimination and difficulty to be MVN but am not sure about the dimensionality $\endgroup$ Commented Oct 4, 2022 at 14:23
  • 1
    $\begingroup$ The IRT example you posted is several indicators of a single person-level trait, analogous to a single random intercept in a MLM. Actually, the observations are cross-classified under persons and items, but the item factor is fixed (directly estimating item-specific discriminations/loadings) rather than random. Skrondal & Rabe-Hasketh's book "Generalized Latent Variable Modeling" sheds a lot of light on a the relation between random-effect and latent-trait / common-factor models. Anyway, IRT with multiple latent traits would require multinormal priors for person parameters. $\endgroup$
    – Terrence
    Commented Oct 5, 2022 at 17:10
  • $\begingroup$ That book is terrific but unfortunately the notation at times is a little confusing since it's often not clear what is a matrix or vector. I also found De Boeck's Explanatory Item Response Model's book helpful but it doesn't have good 2-PL/CFA examples to extend to multivariate longitudinal models with only 2 time periods (what I'm working with). I think your post here and revisiting the "adventures in covariance" section of Statistical Rethinking has cleared this up a lot. Many thanks! $\endgroup$ Commented Oct 5, 2022 at 22:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.