Priors for Random-Effects Covariance Matrices in Joint Models

Why Parameterize the Covariance Matrix?

In a joint model, subject-specific random effects describe variation in latent longitudinal trajectories across subjects. Their covariance matrix determines both the amount of variation in each random effect and the dependence among random effects. Those features can affect the trajectories used in the survival component, so a prior on the covariance matrix is a substantive modeling choice, not just a technical requirement.

Let $b_i$ be the $q$-dimensional vector of random effects for subject $i$. In the standard Gaussian random-effects specification, $b_i$ has mean zero and covariance matrix $D$. The matrix $D$ must be symmetric positive definite. Directly assigning a prior to all its entries can make it difficult to ensure that condition and to understand the prior implications. A useful alternative separates random-effect standard deviations from their correlations.

Definition: Write the random-effects covariance matrix as $$ D = S R S,\qquad S=\operatorname{diag}(\tau_1,\ldots,\tau_q), $$ where each $\tau_j>0$ is the standard deviation of random effect $j$, and $R$ is a $q\times q$ correlation matrix. Thus $R_{jj}=1$, its off-diagonal entries are correlations, and $R$ must be positive definite. This decomposition gives $D_{jk}=\tau_j\tau_k R_{jk}$.

A common prior specification assigns proper priors to the positive scales $\tau_j$, a prior to $R$, and makes the scales independent of $R$ a priori. This modular construction lets the modeler express scale knowledge in the units of each random effect and correlation knowledge separately. It does not imply that scales and correlations will be independent in the posterior.

Scale Priors and Units

Half-normal and half-$t$ distributions are standard choices for random-effect standard deviations. If $Z\sim N(0,A^2)$, then $|Z|$ has a half-normal distribution with scale $A$. A half-$t$ distribution is the distribution of the absolute value of a centered $t$ random variable with a specified degrees of freedom and scale. Both have support only on positive values. The half-$t$ has heavier tails, which can allow larger scales than a half-normal with a comparable central spread.

The scale parameter must be chosen in context. If a random intercept is measured in biomarker units, its prior scale is in biomarker units; if a random slope is measured in biomarker units per year, its prior scale is in those units per year. Changing the time unit changes the numerical scale of a slope and its covariance with an intercept. Therefore, a prior that seems weak on one measurement scale may be informative on another.

Worked Example: Interpreting a Half-Normal Scale Prior

Suppose a random intercept is measured in biomarker units and receives a half-normal prior with scale 5. The prior probability that its standard deviation $\tau_1$ exceeds 10 is $$ P(\tau_1>10)=P(|Z|>10),\qquad Z\sim N(0,5^2). $$ Standardizing gives $10/5=2$, so $$ P(\tau_1>10)=2\{1-\Phi(2)\}\approx 2(0.02275)=0.0455. $$

Thus this prior assigns about 4.55% probability to a random-intercept standard deviation above 10 biomarker units. That statement describes the prior, not the probability that the standard deviation exceeds 10 after observing data. Whether 5 is a sensible scale depends on plausible between-subject variation in this biomarker and on the outcome scale used in the longitudinal model.

Correlation Priors and Positive Definiteness

A correlation matrix is constrained: every pairwise correlation must lie between $-1$ and $1$, but in dimensions above two those pairwise bounds alone do not ensure that the whole matrix is positive definite. Priors on $R$ must respect the joint constraint. The LKJ family is a standard choice, with density proportional to a power of the determinant.

Definition: An LKJ prior with parameter $\eta>0$ on a $q\times q$ correlation matrix has density $$ p(R)\propto |R|^{\eta-1} $$ over positive-definite correlation matrices. Larger values of $\eta$ place more prior weight on matrices with larger determinant, generally favoring correlations closer to zero. When $q=2$, the single correlation $\rho$ has density proportional to $(1-\rho^2)^{\eta-1}$ for $-1<\rho<1$.

For $q=2$, $\eta=1$ makes $\rho$ uniform on $(-1,1)$. For larger $q$, $\eta=1$ is uniform over the space of valid correlation matrices; it does not mean that all pairwise correlations are mutually independent and uniform. For $\eta>1$, the prior increasingly favors correlations nearer zero. The prior remains capable of assigning weight to strong correlations, but the amount depends on $\eta$ and the dimension.

Worked Example: A Two-Effect LKJ Prior

For two random effects, let $\eta=2$. The density of their correlation is proportional to $1-\rho^2$. Its normalizing constant is obtained from $$ \int_{-1}^{1}(1-\rho^2)\,d\rho =2\left(1-\frac{1}{3}\right)=\frac{4}{3}, $$ so $p(\rho)=\frac{3}{4}(1-\rho^2)$ for $-1<\rho<1$.

The prior probability of a correlation between $-0.5$ and $0.5$ is therefore $$ P(-0.5<\rho<0.5)=\frac{3}{4}\int_{-0.5}^{0.5}(1-\rho^2)\,d\rho =\frac{3}{4}\left(1-\frac{1}{12}\right)=0.6875. $$ Thus this prior assigns 68.75% probability to that interval. It is more concentrated near zero than the uniform prior, which assigns 50% probability to the same interval. This calculation describes the two-effect case; it should not be transferred unchanged to a higher-dimensional LKJ prior.

Building and Checking a Covariance Matrix

The decomposition makes it straightforward to construct a covariance matrix once the standard deviations and correlations are specified. For any valid positive-definite $R$ and positive scales, $D=SRS$ is positive definite: for any nonzero vector $x$, $x^\mathsf{T}Dx=(Sx)^\mathsf{T}R(Sx)>0$, because $S$ is invertible and $R$ is positive definite. This is one reason to parameterize a prior through scales and a valid correlation matrix rather than independent entries of $D$.

Worked Example: From Scales and Correlation to Covariance

Suppose a hypothetical model has a random intercept measured in biomarker units and a random slope measured in biomarker units per year. Let their standard deviations be $\tau_1=10$ biomarker units and $\tau_2=2$ biomarker units per year, with correlation $\rho=-0.4$. Then $$ S=\begin{pmatrix}10&0\\0&2\end{pmatrix}, \qquad R=\begin{pmatrix}1&-0.4\\-0.4&1\end{pmatrix}. $$ The covariance matrix is $$ D=SRS= \begin{pmatrix}100&-8\\-8&4\end{pmatrix}. $$

The off-diagonal entry is $(-0.4)(10)(2)=-8$, in units of biomarker squared per year. The diagonal entries are $10^2=100$ biomarker-units squared and $2^2=4$ biomarker-units squared per year squared. The determinant is $100(4)-(-8)^2=336>0$, and the first diagonal entry is positive, so this two-dimensional symmetric matrix is positive definite. The negative correlation represents a prior or model feature in which subjects with higher random intercepts tend to have lower random slopes; it does not by itself establish a causal relationship.

Why an Inverse-Wishart Prior Can Behave Differently

A traditional alternative is to assign an inverse-Wishart prior directly to $D$. Under one common parameterization, $D\sim\operatorname{IW}_q(\Psi,\nu)$ has density proportional to $$ |D|^{-(\nu+q+1)/2}\exp\left\{-\frac{1}{2}\operatorname{tr}(\Psi D^{-1})\right\}, $$ where $\Psi$ is positive definite and $\nu>q-1$. Under this parameterization, the mean exists for $\nu>q+1$ and is $E(D)=\Psi/(\nu-q-1)$.

The inverse-Wishart prior is compact to specify, but it does not independently set priors for standard deviations and correlations. Its parameters jointly shape these quantities. In particular, changing $\nu$ to alter the concentration of the prior also changes the induced marginal distributions of the variances and correlations. A diagonal $\Psi$ does not make the standard deviations and correlations independent a priori.

Worked Example: Reading an Inverse-Wishart Prior

Consider a two-dimensional random-effects covariance matrix with $D\sim\operatorname{IW}_2(\Psi,5)$ and $$ \Psi=\begin{pmatrix}2&0\\0&0.5\end{pmatrix}. $$ Since $5>2+1$, the prior mean exists and is $$ E(D)=\frac{\Psi}{5-2-1} =\frac{1}{2}\begin{pmatrix}2&0\\0&0.5\end{pmatrix} =\begin{pmatrix}1&0\\0&0.25\end{pmatrix}. $$

The mean prior variances are therefore 1 and 0.25 in the corresponding squared units. They are not prior means of the standard deviations: in general $E(\sqrt{D_{jj}})\ne\sqrt{E(D_{jj})}$. Moreover, the diagonal $\Psi$ does not assert that the random-effect correlation is fixed at zero. The prior remains a distribution over covariance matrices, including matrices with nonzero off-diagonal entries.

A further way to see why the decomposition matters is that a prior density depends on its parameterization. The transformation from $D$'s distinct entries to standard deviations and correlations has Jacobian $$ dD=2^q\left(\prod_{j=1}^{q}\tau_j^q\right) \left(\prod_{j=1}^{q}d\tau_j\right) \left(\prod_{j<k}dR_{jk}\right). $$ For $q=2$, this follows by differentiating $D_{11}=\tau_1^2$, $D_{22}=\tau_2^2$, and $D_{12}=\tau_1\tau_2\rho$; the absolute determinant of that derivative is $4\tau_1^2\tau_2^2$. Consequently, a density that is flat in one set of coordinates is not generally flat in the other. Priors should be specified on the quantities whose behavior is intended, rather than described vaguely as “noninformative.”

Practical Choices, Diagnostics, and Limitations

A useful workflow is to choose scales on interpretable random-effect units, select a correlation prior that allows scientifically plausible dependence, and then examine the implied covariance matrices. Prior simulation can reveal whether the chosen scales routinely imply implausibly large trajectory variation or whether the correlation prior concentrates more strongly near zero than intended. For each draw, the constructed $D$ should be symmetric positive definite. If a draw or parameterization produces an invalid matrix, the implementation is not respecting the intended support.

The main alternatives answer different modeling needs. A diagonal $D$ sets random-effect correlations to zero, reducing the number of covariance parameters but ruling out residual between-subject association among those effects. A full covariance matrix with separate scale priors and an LKJ prior allows more flexible dependence and makes the scale and correlation choices visible. An inverse-Wishart prior is a single conventional prior on the full matrix, but its induced marginal behavior can be less direct to calibrate. No option is universally best: the choice depends on the number of random effects, the available information, and meaningful scales for the application.

These priors do not remove the possibility that the data weakly identify some covariance components. With limited longitudinal information, especially for random slopes, posterior inferences may remain sensitive to scale and correlation choices. Strong posterior concentration at a boundary, markedly different results under reasonable scale priors, or simulated prior trajectories that are far outside plausible ranges are reasons to revisit the specification. Rescaling time or the biomarker may improve interpretability, but the prior scales must be adjusted consistently.

Mistake: assigning independent uniform priors to every off-diagonal correlation. In dimensions above two, arbitrary collections of pairwise correlations need not form a positive-definite matrix. Use a prior supported on valid correlation matrices, such as LKJ.
Mistake: treating a half-$t$ or half-normal scale as automatically weak. Its practical meaning depends on its scale relative to the units and plausible variation. Calibrate it by examining the implied random-effect variation.
Mistake: reading an inverse-Wishart scale matrix as a set of independent variance priors. The inverse-Wishart prior shapes the covariance components jointly. Examine the induced marginal behavior or use a standard-deviation/correlation decomposition when separate calibration is important.

Key takeaway: The representation $D=SRS$ separates the scales of random effects from their dependence while preserving a valid covariance matrix. Proper scale priors and a positive-definite correlation prior make the assumptions more interpretable; an inverse-Wishart prior is convenient but can impose less transparent joint behavior.

Check Your Understanding

Answer each question, then compare with the solutions below.

If a random intercept has standard deviation 3 and a random slope has standard deviation 0.5, with correlation 0.6, what is their covariance?
For two random effects under LKJ($\eta=1$), what is the prior probability that their correlation lies between $-0.25$ and $0.25$?
Why do independent priors on every off-diagonal entry fail to guarantee a valid correlation matrix when there are three or more random effects?
Under $\operatorname{IW}_2(\Psi,6)$, with $\Psi=\operatorname{diag}(3,1.5)$, calculate the prior mean of $D$.

Solutions. (1) The covariance is $0.6(3)(0.5)=0.9$, in the product units of the intercept and slope. (2) LKJ($\eta=1$) gives a uniform correlation on $(-1,1)$ when $q=2$, so the probability is the interval length divided by 2: $0.5/2=0.25$. (3) Pairwise values individually in $(-1,1)$ can still combine into a matrix with a nonpositive determinant, so they may violate positive definiteness. (4) Since $6>2+1$, $E(D)=\Psi/(6-2-1)=\operatorname{diag}(1,0.5)$, in the squared units of the two random effects.