A Graphical Framework for Interpretable Correlation Matrix Models

1]Anna Freni-Sterrantinomycorrespondingauthor
2]Denis Rustand
2]Janet van Niekerk
2]Elias Teixeira Krainski
2]Håvard Rue
[1]The Alan Turing Institute, 2QR, John Dodson House, 96 Euston Rd, NW1 2DB London, United Kingdom
[2]Statistics Program, CEMSE Division, King Abdullah University of Science and Technology Thuwal, 23955, Makkah, Saudi Arabia

[mycorrespondingauthor]Corresponding author. Email address<EMAIL_ADDRESS>
In this work, we present a new approach for constructing models for correlation matrices with a user-defined graphical structure. The graphical structure makes correlation matrices interpretable and avoids the quadratic increase of parameters as a function of the dimension. We suggest an automatic approach to define a prior using a natural sequence of simpler models within the Penalized Complexity framework for the unknown parameters in these models.

We illustrate this approach with three applications: a multivariate linear regression of four biomarkers, a multivariate disease mapping, and a multivariate longitudinal joint modelling. Each application underscores our method's intuitive appeal, signifying a substantial advancement toward a more cohesive and enlightening model that facilitates a meaningful interpretation of correlation matrices.

Complexity penalized priorsCorrelation Matrix Modeling Graphical StructureMultivariate Joint Modeling

§ INTRODUCTION

Estimating the covariance matrix and specifying a covariance prior pose challenges from a Bayesian perspective.
Traditionally, various natural prior conjugate distributions have been used, such as the Wishart <cit.>, inverse-Wishart <cit.> and scaled inverse-Wishart <cit.>.  However, the inverse-Wishart prior and its generalization are associated with certain issues. For instance, the uncertainty for all variances is controlled by a single degree of freedom parameter and the marginal distribution for the variances has a low density in a region near zero. There is also a priori dependence between correlations and variances, where larger variances are associated with correlations near $\pm$1, whereas smaller variances are associated with correlations near zero <cit.>.
In the era of computational frameworks, the attractiveness of preserving prior conjugacy has diminished, as Bayesian inference is mainly conducted using either sampling-based methods like Markov Chain Monte Carlo (MCMC), approximate methods such as integrated nested Laplace approximation (INLA) or variational Bayes.

The covariance matrix contains information on two distinct fronts: dependencies and marginal variances. The correlation matrix is calculated based on the covariance matrix from which the nature and strength of pairwise linear dependence can be extracted.
An alternative view is to consider two priors: one for the marginal variances and one for the correlation.
For the latter, two major challenges emerge (i) the unit diagonal and positive definite constraints; (ii)  the number of parameters grows quadratically as the
number of variables increases.  The typical approach assumes correlation among variables (or random effects), but these are not explicitly elicited.  Hence, how variables (or random effects) are correlated is unknown and unaccounted for in the prior choice.

Proposals for correlation priors include priors based on the spectral decomposition of a symmetric matrix <cit.> and priors on the singular value decomposition of an asymmetric matrix <cit.>. Both are characterized by excessive informative constraints on the correlation matrix elements, leading to shrinkage estimation towards the matrix diagonal.
Another proposal is the  LKJ  prior <cit.>, which is a weakly informative prior defined on a simplex. By tuning the strength of just one parameter allows to control how closely the sampled matrices resemble the identity matrices.

Nonetheless, these proposals' priors tend to shrink the estimation towards the identity matrix. This feature can be problematic as it excludes models that have varying degrees of correlation and are not limited to just being fully correlated or fully uncorrelated.
Therefore, there is a need for more interpretable models that can account for correlation matrices with different structures beyond the conventional fully correlated or uncorrelated matrices.
Additionally, high-dimensional models pose a challenge, thus requiring dimension-reduction techniques for estimable parameters in correlation matrices.

We develop a framework based on a user-defined graphical construction addressing two main challenges. Firstly, we reduce the dimension of the estimable parameters of the correlation matrix by formulating a latent model. Secondly, we use this graphical construction to elicit the structure of the correlation matrix, and subsequently, we develop an intuitive prior that shrinks to a sequence of simpler models.

We formulate a latent model to control the quadratic growth of parameters for high-dimensional models, as the covariance matrix $\mm{\Sigma}$ tends to become large and can involve many correlations and marginal variances to be estimated. Thus,  we can explain some of the correlations between response variables through their relationship to latent factors.  One class of latent models is the multivariate Gaussian Markov random field (MGMRF) model.
These models find their most popular application in spatial statistics, i.e. in multivariate disease mapping, using multivariate conditional autoregressive model <cit.>; in
longitudinal and survival analysis <cit.>, spatio-temporal datasets <cit.>, and describing nature's and systems' multivariate spatial and temporal dynamics <cit.>.

Based on the elicited correlation structure from the graphical construction, we adopt the penalizing complexity (PC) prior <cit.> framework to develop a prior based on a sequence of simpler correlation matrices models. This prior contracts toward a simpler model (base model), that is not necessarily the identity matrix regulated by a user-specified rate parameter.

The remainder of this paper is the following: in Section <ref>, we show how to derive correlation matrices from a graphical construction.  Then, in Section  <ref>, we describe how we derive a sequence of simpler models for correlation matrices using the graphical framework.  We present the penalized complexity priors principles in Section <ref>, and we define the prior correlation matrix. In Section <ref>, using the proposed modelling approach, we illustrate three data applications of multivariate analysis,  and concluding remarks are given in Section <ref>.

§ GRAPH-DERIVED CORRELATION MATRICES

§.§ Correlation and graphs

Correlation is interpreted as the association between two variables by quantification of the degree of joint variability. Often two variables are correlated because they share a common factor.
This perspective suggests that the observed correlation between two variables may arise from both variables being influenced by latent or unobservable factors. Within this theoretical framework, the observed correlation does not necessarily signify a direct causal relationship between the variables; rather, it underscores their collective responsiveness to shared influences.

To summarize the correlation between variables, the graphical Gaussian models (GGM) <cit.> offer a visual and computational framework to organize the conditional independence in sparse matrices.  GGMs focus on modelling the relationships between observed variables directly. Each observed variable corresponds to a node in the graphical representation, and edges in the graph represent conditional dependencies between variables. On the other hand, the absence of an edge between them implies that they are conditionally independent, given all other variables in the model, corresponding to a zero entry in the precision matrix

However, if the correlation is constructed as two variables sharing a common effect, GGMs are extended to include latent (unobserved) variables. These latent variables capture shared information among the observed variables, allowing the model to represent more complex relationships. Latent GGMs provide greater flexibility in capturing dependencies,
allowing for a more comprehensive representation of the underlying structure in the data. By capturing heterogeneity in the data by allowing for the presence of unobserved factors to contribute to variability. This is particularly useful in scenarios where subgroups or hidden patterns influence the relationships among variables.

Based on the Latent GGM, we define a tree graph where the correlations between two nodes are given by sharing a latent ancestor node, which will induce the children's correlation. Edges connect only a latent node to children or to other latent nodes, with no edges between children and only one edge to the latent parent,  with a rooted node as a latent factor.
From the graphical representation, we factorize the joint in terms of local conditional probabilities. This provides a sparse matrix for the precision matrix, with a dense block for the observed (children) nodes.
While sparse matrices come with their computational benefits,  the dense matrix obtained by this graphical structure leverages the advantages of latent GGMs and offers an interpretable framework for understanding the intricate relationships within the data. The inclusion of the parent latent nodes avoids specifying directly the correlation among children and guarantees a valid correlation matrix.

For interpretability and illustration purposes, we show how we can derive the correlation matrices from graphs with up to three latent factors (see Figures <ref>, <ref> and <ref>), but the procedure is trivially extended to a higher number of variables and latent factors, as shown in data applications.

§.§ A correlation matrix with one parameter

[shorten >=1pt,node distance=2cm,on grid,auto]
state=[shape=circle,thick,draw,minimum size=0.6cm]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
[statesq, fill=rgb:black,1;white,2] (b) $p_1$;
[state,below left of=b] (Z1) $c_1$;
[state,below of=b](Z2)  $c_2$;
[state,below right of=b] (Z3) $c_3$;
(b) edge node  (Z1)
(b) edge node  (Z2)
(b) edge node  (Z3);

One parameter graph with two variables.  Grey and white nodes represent latent factors and observed variables, respectively.

Given the (directed) graph in Figure <ref>,  we assume that the children nodes $c_i, i=1,2,3$ representing the observed variables and one parent latent factor $p_1$,  are distributed as:

\begin{equation}
\label{eq:eqdis1}
\begin{split}
c_i|p_1 \stackrel{\text{iid}}{\sim} \mathcal{N}(p_1,1)\\
p_1 \sim \mathcal{N}(0,q_1^2).
\end{split}
\end{equation}

Using the law of total expectation, variance and covariance, we can integrate out the latent factor. We get that $\E(c_i) = \E(\E(c_i|p_1))=0$,
$\Var(c_i) = \Var(\E(c_i|p_1)) + \E(\Var(c_i|p_1)) = 1 + q_1^2$, and for $i\not=j$, $\Cov(c_i, c_j) = \Cov(\E(c_i|p_1), \E(c_j|p_1)) + \E(\Cov(c_i, c_j | p_1)) = q_1^2$. The correlation between $c_i$ and $c_j$, is then

\[
\Cor(c_i, c_j) = \frac{q_1^2} {1 + q_1^2}, \quad i\not=j.
\]

However, the idea is to use this graphical approach where correlations are introduced by sharing (random) latent factors in the graph.

It turns out it is easier to construct the correlation matrix from the joint (Gaussian) distribution. As the $c_i$ are independent given $p_1$, the  joint distribution for this Gaussian Markov Random Field (GMRF)<cit.> for a fixed $q_1^2$ is:

\[
\pi(c_1c_2,c_3,p_1)\propto \pi(c_1|p_1)\pi(c_2|p_1)\pi(c_3|p_1)\pi(p_1)\propto
\exp\left(-\frac{1}{2}\left(\frac{p_1^2}{q_1^2}+\sum_{i=1}^{2}(c_i-p_1)^2\right)\right).
\]

To obtain the precision matrix for this GMRF, we need to compute the Hessian of the negative logarithm of this distribution.  The Hessian matrix carries the conditional information of the random vector  $\mm{\theta}=(c_1,c_2,c_3,p_1)$, when all the other parameters are fixed.  This also follows by the definition of a GMRF (Theorem 2.2 in <cit.>) where the diagonal elements of the precision matrix express the conditional precision given the other variables; the off-diagonal elements correspond to pairwise correlation between the nodes of a graph (see Figure <ref>), conditional on the other nodes.
The precision matrix is then:

\begin{equation*}\label{eq:eqn2}%
\mm{\Sigma}^{-1}=\begin{bmatrix}{}
1 & . & . & -1 \\
. & 1 & . & -1 \\
. & . & 1 & -1 \\
-1 & -1 & -1  & 3+\frac{1}{q_1^2}
\end{bmatrix}.
\end{equation*}

Zero elements are indicated as '$.$' and $(3+1/q_1^2)$ represents conditional precision between the latent factor and $c_i$'s.
The variance-covariance matrix is:

\begin{equation*}\label{eq:eqn3}%
\mm{\Sigma}=\begin{bmatrix}{}
1+q_1^2 &q_1^2 &q_1^2  & q_1^2 \\
q_1^2 & 1+q_1^2 &  q_1^2& q_1^2\\
q_1^2 &  q_1^2& 1 +q_1^2&  q_1^2\\
q_1^2&  q_1^2& q_1^2 &  q_1^2
\end{bmatrix}
\end{equation*}

whose diagonal elements are marginal node variances, the correlation matrix is:

\begin{equation}
\label{eq:norm}
\mm{C}=\Diag(\mm{\Sigma})^{-1/2} \mm{\Sigma}~\Diag(\mm{\Sigma})^{-1/2}.
\end{equation}

Finally,  we extract  the left upper block of dimension $ 3 \times 3 $ of $C$ (i.e., the number of child nodes):

\begin{equation*}
\label{eq:C1}%
\mm{C}[1:3,1:3]=\begin{bmatrix}{}
1    &  \frac{q_1^2}{1+q_1^2}& \frac{q_1^2}{1+q_1^2} \\
\frac{q_1^2}{1+q_1^2}  & 1  & \frac{q_1^2}{1+q_1^2 } \\
\frac{q_1^2}{1+q_1^2}   &   \frac{q_1^2}{1+q_1^2 } & 1 \\
\end{bmatrix}
\end{equation*}

where $\rho$ the correlation parameter $\rho=\frac{q_1^2}{1+q_1^2}$ is a  function only of $q_1$. It is immediately obvious that is an exchangeable correlation matrix is generated by a graph with one parent and multiple children.

§.§ A correlation matrix with two parameters

[shorten >=1pt,node distance=2.5cm,on grid,auto]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
state=[shape=circle,thick,draw,minimum size=0.6cm]
[statesq, fill=rgb:black,1;white,2] (b1) $p_2$;
[statesq, above right of=b1, yshift=10mm ,fill=rgb:black,1;white,2] (b2) $p_1$;
[state, below right of=b1] (Z2) $c_2$;
[state,below left of=b1] (Z1) $c_1$;
[state, below right of=b2] (Z3) $c_3$;
(b1) edge node  (Z1)
(b1) edge node  (Z2)
(b2) edge node  (Z3)
(b2) edge node   (b1);

Two parameters graph with three variables. Grey and white nodes represent latent factors and observed variables.

With two latent factors, as depicted in Figure <ref>, the relationship is between variables $c_1,c_2,c_3$ and $p_1,p_2$. We specify the conditional distribution between latent factors and observed,  as follows:

\[
\begin{matrix}
c_i|p_2 \stackrel{\text{iid}}{\sim} \mathcal{N}(p_2,1), i=1,2\\
c_3|p_1 \sim \mathcal{N}(p_1,1)\\
p_1 \sim \mathcal{N}(0, q_1^2)\\
p_2|p_1\sim \mathcal{N}(p_1,q_2^2)\\
\end{matrix}
\]

while the marginal and conditional distributions are the same as stated in equation <ref>.
The Latent GGM has a joint distribution (for fixed $q_1^2,q_2^2$):

\[
\pi(c_1,c_2,c_3,p_2,p_1)\propto \exp\left(-\frac{1}{2}\left(\frac{p_1^2}{q_1^2}+\sum_{i=1}^{p=2}(c_i-p_1)^2+(c_3-p_2)^2+\frac{(p_2-p_1)^2}{q_2^2}
\right)\right)
\]

with precision matrix:

\begin{equation*}\label{eq:equ2}%
\mm{\Sigma}^{-1}=\begin{bmatrix}{}
\begin{array}{ccccc}
1 & 0 & 0 & -1 & 0
\\
0 & 1 & 0 & -1 & 0
\\
0 & 0 & 1 & 0 & -1
\\
-1 & -1 & 0 & 2+\frac{1}{q_{2}^{2}} & -\frac{1}{q_{2}^{2}}
\\
0 & 0 & -1 & -\frac{1}{q_{2}^{2}} & \frac{1}{q_{1}^{2}}+1+\frac{1}{q_{2}^{2}}
\end{array}
\end{bmatrix}
\end{equation*}

and, covariance matrix:

\begin{equation*}%
\mm{\Sigma}=\begin{bmatrix}{}
\begin{array}{ccccc}
q_{1}^{2}+q_{2}^{2}+1 & q_{1}^{2}+q_{2}^{2} & q_{1}^{2} & q_{1}^{2}+q_{2}^{2} & q_{1}^{2}
\\
q_{1}^{2}+q_{2}^{2} & q_{1}^{2}+q_{2}^{2}+1 & q_{1}^{2} & q_{1}^{2}+q_{2}^{2} & q_{1}^{2}
\\
q_{1}^{2} & q_{1}^{2} & q_{1}^{2}+1 & q_{1}^{2} & q_{1}^{2}
\\
q_{1}^{2}+q_{2}^{2} & q_{1}^{2}+q_{2}^{2} & q_{1}^{2} & q_{1}^{2}+q_{2}^{2} & q_{1}^{2}
\\
q_{1}^{2} & q_{1}^{2} & q_{1}^{2} & q_{1}^{2} & q_{1}^{2}
\end{array}
\end{bmatrix}.
\end{equation*}

Taking the upper left $3 \times 3$ matrix and the correlation matrix is:

\begin{equation*}\label{eq:C2}
\left[\begin{array}{ccc}
1 & \frac{q_{1}^{2}+q_{2}^{2}}{q_{1}^{2}+q_{2}^{2}+1} & \frac{q_{1}^{2}}{\sqrt{q_{1}^{2}+q_{2}^{2}+1}\, \sqrt{q_{1}^{2}+1}}
\\
\frac{q_{1}^{2}+q_{2}^{2}}{q_{1}^{2}+q_{2}^{2}+1} & 1 & \frac{q_{1}^{2}}{\sqrt{q_{1}^{2}+q_{2}^{2}+1}\, \sqrt{q_{1}^{2}+1}}
\\
\frac{q_{1}^{2}}{\sqrt{q_{1}^{2}+q_{2}^{2}+1}\, \sqrt{q_{1}^{2}+1}} & \frac{q_{1}^{2}}{\sqrt{q_{1}^{2}+q_{2}^{2}+1}\, \sqrt{q_{1}^{2}+1}} & 1
\end{array}\right]
\end{equation*}

\[
\rho_1=\frac{q_1^2+q_2^2}{q_1^2+q_2^2+1}
\]

\[
\rho_2= \frac{q_{1}^{2}}{\sqrt{q_{1}^{2}+q_{2}^{2}+1}\, \sqrt{q_{1}^{2}+1}}.
\]

It follows that the correlation matrix elements are functions of $q_1$ and $q_2$. With $\rho_1$ indicating the correlation between $c_1$ and $c_2$, and $\rho_2$ indicates correlation of $c_3$ with the children for $p_2$.

§.§ A correlation matrix  with three parameters

[shorten >=1pt,node distance=2.5cm,on grid,auto]
state=[shape=circle,thick,draw,minimum size=0.6cm]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
[statesq,fill=rgb:black,1;white,2] (c1) $p_1$;
[statesq,fill=rgb:black,1;white,2, below left of=c1, yshift=5mm] (b1) $p_2$;
[statesq,fill=rgb:black,1;white,2, below right of=c1, yshift=5mm] (b2) $p_3$;
[state, below of=b2] (Z4) $c_4$;
[state,below left of=b1] (Z1) $c_1$;
[state, below of=b1] (Z2) $c_2$;
[state, below right of=b1] (Z3) $c_3$;

(c1) edge node  (b1)
(c1) edge node  (b2)
(b2) edge node  (Z4)
(b1) edge node  (Z1)
(b1) edge node  (Z2)
(b1) edge node  (Z3) ;

Three parameters graph with three variables. Grey and white nodes represent latent factors and observed variables, respectively.

Finally, the third case is represented by a graph in Figure <ref>. The latent factor $p_1$ induces correlation between the four variables $c_i$ with $i=1,..,4$  and on the two latent factors $p_2$ and $p_3$. We consider the conditional distributions:

\[
\begin{matrix}
c_i|p_2 \stackrel{\text{iid}}{\sim} \mathcal{N}(p_2,1), i=1,2,3\\
c_4|p_3 \sim \mathcal{N}(p_3,1)\\
p_1 \sim \mathcal{N}(0, q_1^2)\\
p_2|p_1\sim \mathcal{N}(p_1,q_2^2)\\
p_3|p_1 \sim \mathcal{N}(p_1,q_3^2)\\
\end{matrix}
\]

therefore, the joint (for fixed $q_1^2,q_2^2,q_3^2$) is:

\[
\pi(c_1,c_2,c_3,c_4,p_2,p_3,p_1)\propto\\ \exp\left (-\frac{1}{2} \left (\frac{p_1^2}{q_1^2} +\sum_{i=1}^{p=3}(c_i-p_2)^2+(c_4-p_3)^2+
\frac{(p_2-p_1)^2}{q_2^2} + \frac{(p_3-p_1)^2}{q_3^2}\right) \right)
\]

it follows that the precision matrix is:

\begin{equation*}\label{equsub2}%
\mm{\Sigma}^{-1}=\begin{bmatrix}{}
1 & . & . & . & -1 &.      & .  \\
. & 1 & . & . & -1 & .      & .  \\
. & . & 1 & . & -1 &       &.  \\
. & . & . & 1 & .   & -1  & . \\
-1 & -1 &-1 & . & 3+1/q_2^2 & .&-1/q_2^2 \\
.  & . &. &-1 &. & 1+1/q_3^2 &-1/q_3^2 \\
.  & . & . & .  & -1/q_2^2& -1/q_3^2 & 1/q_1^2+1/q_2^2+1/q_3^2\\
\end{bmatrix}
\end{equation*}

giving  rise to the  correlation matrix

\begin{equation*}\label{eqc2}%
1        &  \rho_{1}  & \rho_{1}    &  \rho_{2}   \\
\rho_{1}  & 1 &  \rho_{1}     &  \rho_{2}    \\
\rho_{1}  & \rho_{1}    & 1  &  \rho_{2}  \\
\rho_{2} &  \rho_{2}  &   \rho_{2}  &  1 \\
\end{bmatrix}
\end{equation*}

\[
\rho_1=\frac{q_1^2+q_2^2}{q_1^2+q_2^2+1}
\]

\begin{equation}
\label{eqrho2}
\rho_2=\frac{q_1^2}{\sqrt{q_1^2+q_2^2+1}\sqrt{q_1^2+q_3^2+1}}
\end{equation}

$\rho_1$ indicate the correlation between $c_1$ and $(c_2,c_3)$ and $\rho_2$ indicate  correlation between $c_4$ with $p_2$'s children.
In the case of three latent nodes, we observed that the correlation among $(c_1,c_2,c_3)$  is given by $\rho_1=f(q_1^2, q_2^2)$, while $\rho_2=f(q_1^2,q_2^2,q_3^2)$. The latter correlation also accounts for some variability that is in $p_2$.

Extending these models to more variables is straightforward,
It is easy to see that it is sufficient to express the joint distribution of the graph models and that the derived correlation matrix is a function of the latent nodes' variances.

§.§ An explicit result

Using the law of total expectation, variance and covariance, we can directly compute $\Cor(c_i, c_j)$, $i\not=j$ from the graph. This is a simple consequence of the calculations at the beginning of Section <ref>. As an example let us compute $\rho_2 = \Cor(c_1, c_4)$ in Figure <ref>. We get that $\E(c_i)=0$,  $\Var(c_1) = 1 + q_2^2 + q_1^2$, $\Var(c_4) = 1 + q_3^2 + q_1^2$ simply by following the reverse paths in the graph from $c_1$ and $c_4$ to $p_1$. For the covariance, then $\Cov(c_1, c_4) = q_1^2$, as $p_1$ is the first parent node where the reverse paths from $c_1$ and $c_4$ to $p_1$ meet, and $\Var(p_1) = q_1^2$. These results gives $\rho_2$ in equation <ref>.

§ CONTRACTIONS FROM USER-DEFINED CORRELATION MATRICES TO IDENTITY MATRIX

The correlation coefficients are functions of the parent's variances. By further exploiting the graph, we can automatically generate a sequence of simpler correlation matrices by removing one parent at a time.  Based on this approach,  there are as many simpler models as there are parent nodes that will represent parsimonious models.
It is clear why the identity matrix could be a reasonable choice only when the correlation matrix is an exchangeable matrix.  When the correlation matrix is unstructured, choosing an uncorrelated as a simpler model might cause too much shrinkage. Hence, a solution is to have the user-defined graph to guide the definition of a sequence of nested models that represents a simpler version of the initial model, where every simpler model has one parameter less than the previous model.

The process automatically produces a sequence of correlation matrices by nullifying one parent's variance at a time.
In Figure <ref>$a$  (shown later), we show an example of how, from the graph with eight children and seven parents,  by removing, in order, the parent node, we generate a sequence of seven graphs,  with the last two devising an exchangeable and an uncorrelated matrix.
We illustrate this in the following example.

§.§ Contacting correlation matrices from three parameters to one parameter

The graph construction defines the complex model correlation matrix, and we can use it to define a sequence of simpler models. A sequence of nested graphs is automatically generated by removing a parent at each step and setting the parent node's variance equal to zero. Children whose parents are removed, become children of the antecedent parents.
The parents' order governs the reverse removal ordering, with parent $p_1$ being the last to be removed. Hence, this ordering impacts the magnitude of the correlation among children.
Here, we compare the correlation between three sets of children between graphs $e$, the complex model, and $f,g$, the simpler models, from Figure <ref>. The correlation equations are the following:

$\overbrace{\rule{4.5cm}{0pt}}^{\text{Graph } e} $      $\overbrace{\rule{4.5cm}{0pt}}^{\text{Graph } f~(q_3^2=0)}$       $\overbrace{\rule{4.5cm}{0pt}}^{\text{Graph } g~(q_2^2=0)} $

$\rho_{children 12}= \frac{q_1^2+q_2^2}{1+q_1^2+q_2^2} $     $\rho_{children 12}= \frac{q_1^2+q_2^2}{1+q_1^2+q_2^2} $

$\rho_{children 15} = \frac{q_1^2}{\sqrt{1+q_1^2+q_2^2}\sqrt{1+q_1^2+q_3^2}} $     $\rho_{children 15} = \frac{q_1^2}{\sqrt{1+q_1^2+q_2^2}\sqrt{1+q_1^2}}$           $\rho_{Exchangeable}=\frac{q_1^2}{q_1^2+1}$

$\rho_{children 56} = \frac{q_1^2+q_3^2}{1+q_1^2+q_3^2}$     $\rho_{children 56} = \frac{q_1^2}{1+q_1^2}$

where $\rho_{children 56}$, form graph $f$ equals $\rho_{Exchangeable}$.
In Figure <ref>, we  plotted the correlations   values for $q_1^2=c(0.5,1,2,5,10)$, $0\leq q_2^2\leq 10$ and $q_3^2=2$.
In the complex model panels (graph $e$),  correlation $\rho_{children 12}$ and $\rho_{children 56}$  dominate over $\rho_{children 15}$.

In the simpler model panels (graph $f$), when $q_3^2=0$,  $\rho_{children 12}$  continues to dominate over $\rho_{children 15}$, while the correlation  $\rho_{children 56}$ depends only by $q_1^2$ , and increases to match the correlation with $\rho_{children 12}$ for higher values of $q_1^2$. In the graph denoted as $g$, correlations exhibit uniformity and increase with higher values of $q_1^2$. The last simpler model is always the identity matrix for $q_1^2=0$.
This result provides a way to assess how the variables (children nodes) can be grouped and what correlation magnitudes should be expected by the graph construction; hence, it can be reorganised to account for that. In a high-dimensional setting, when the graphs describe the correlation between random slope and intercepts, it helps to grasp the correlation structure.
In summary, the parent removal process gradually contracts toward the uncorrelated matrix.  This insight serves as a foundation for constructing a prior for these model components within the framework of a penalized complexity prior, which are going to be introduced next.

Panel shows the correlation values for 3 sets of children between graph $e$ representing the complex model and graph $f$  representing the simpler model, from Figure <ref>). In both models, the correlation between children of the same parent (solid) $\rho_[children 12]$ dominates the correlation between the sets of children from different parents $\rho_[children 15]$ (dotted), for $q_1^2=0.5, 1,2,5,10$ and  $0\leq q_2^2, \leq 10$.  The correlation of $\rho_[children 56]$, in the complex model is similar to  $\rho_[children 12]$ as $q_1^2$ increases and $q_3^2=2$ . In the simpler model $q_3^2=0$, hence the correlation for $\rho_[children 56]$ is a function only of $q_1^2$.

§ PENALIZED COMPLEXITY PRIOR

The core of the framework relies on specifying a base model (or simple),  that represents the simple version and the complex or flexible model that can be viewed as an extension of the parsimonious base model.
For example, the skew-normal is an extension of the normal model, a random effects model is an extension of a fixed effects model, a Weibull model is an extension of the exponential model and so on. Penalized complexity priors (PC-priors)  are a class of weakly informative priors that penalize departures from the complex model to the base model, as introduced by Simpson et al. (<cit.>).
Briefly, we present the four principles of a PC prior:

* The Occam's razor parsimony. Unless there is evidence to move toward a complex model, the prior should penalise departure from the base model.

* Measure of complexity. To compute the measure of complexity between the base and complex model, the Kullback-Leibler Divergence (KLD)

\[
KLD(\pi(\mm{x}|\xi)|| \pi(\mm{x}|\xi=0))=\int \pi(\mm{x}|\xi)\log\left(\frac{\pi(\mm{x}|\xi)}{ \pi(\mm{x}|\xi=0)}\right)\emph{d}\mm{x},
\]

is used. The KLD measures the information lost when the complex model is approximated with the base model. For two normal distribution $f_0$ and $f_1$  where
$f_i \sim \mathcal{N}(0,\mm{\Sigma}_i)$ $i = 0,1$, while $n$ is the covariance matrix dimension , the KLD is:

\begin{equation*}
\label{eq:kld}
KLD(f_0 || f_1)=\frac{1}{2}\left(tr(\mm{\Sigma}_0^{-1} \mm{\Sigma}_1)-n-\log \left( \frac{|\mm{\Sigma}_1|}{|\mm{\Sigma}_0|}\right)\right).
\end{equation*}

The KLD is transformed to a unidirectional distance measure,

\[
d(f_0 || f_1)=\sqrt{2KLD(f_0 || f_1)}.
\]

By  Pinsker's inequality, this is also an upper bound of the total variation distance between $f_0$ and $f_1$.  The prior is then defined as a function of the distance $\pi(d)$, with a high mass on areas where replacing the flexible model with the base model will not compromise too much loss of information, i.e. where $d=0$.

* Constant rate of penalisation. The prior distance should be chosen in a way that the mode should be located at the base model, while the density decays as the distance from the base model
increases. We chose a constant rate of penalization that satisfies the following:

\[
\frac{\pi_d(d+\delta)}{\pi_d(d)}=r^{\delta},
\]

with $\delta, d \geq 0$ and $0<r<1$. This means that the prior is independent of the actual distance. The resulting is an exponential prior is formulated for the distance

\[
\pi(d)=\lambda \exp(-\lambda d),
\]

where $r=\exp(-\lambda)$, corresponding to a change of variable in the prior for $\xi$ as:

\begin{equation*}
\label{eq:prior}
\pi(\xi)=\lambda \exp(-\lambda d(\xi)) \left|\frac{\partial d(\xi)}{\partial \xi}\right|.
\end{equation*}

* User-defined scaling. The PC-prior has one user-specified rate parameter:$\lambda$. This parameter controls the probability mass at the tail, and the user can thus intuitively specify this value using the following formulation:

\[
Prob(\mathcal{Q}(\xi) > U)=\alpha,
\]

Where $\mathcal{Q}(\xi)$ is a transformation of the flexibility parameter into an interpretable quantity,  $U$ is user-defined upper bound for the tail-event, and $\alpha$ is the weight assigned to the event.

§.§ Penalized complexity prior for the correlation matrix

Despite a correlation matrix being characterized by multiple parameters,  we can formulate a multivariate prior by combining individual univariate conditional priors corresponding to each simpler model (see <cit.> for an example with autoregressive time series).  To illustrate the derivation of the multivariate prior, that is based on the relationships between a  sequence of flexible-base model pairs, we examine the graphical structure outlined in Figure <ref>, as detailed in Section <ref>.

Given the parameters ($q_1^2, q_2^2, q_3^2$), we define a multivariate prior as a product of univariate PC priors based on the distance between models after removing one parent at each step.
The PC prior for the correlation matrix becomes:

$$\pi(q_1^2, q_2^2, q_3^2) = \pi(q_3^2| q_2^2, q_1^2) \cdot \pi(q_2^2| q_1^2) \cdot \pi(q_1^2).$$

This equates to sequentially setting the parents' variance equal to zero. The first step evaluates the distance between the flexible (3 parents) and base model (2 parents), conditional on the first 2 parents. Hence, it evaluates the prior for the parameter $q_3^2$ conditional on $q_1^2$ and $q_2^2$. In the next step, the flexible model is defined by the 2 remaining parents (parent 1 with variance $q_1^2$ and parent 2 with variance $q_2^2$), and the base model only includes parent 1, thus evaluating the prior for $q_2^2$ conditional on $q_1^2$. The last step is based on the distance between the model with the last parent remaining and the independence between the children.

To define the penalized complexity prior for the correlation matrix,
we rely on the relationship between the  Kullback-Leibler divergence and Fisher information computed at the base model for exponential families <cit.>. For each model couple in the sequence, we approximate the KLD as:

\begin{equation}
\label{eq:eqkld}
KLD(\pi(\mm{x}|\xi)|| \pi(\mm{x}|\xi=0))=\frac{1}{2}I(0)\xi^2+\text{ higher order terms} ,
\end{equation}

where $I(0)$ is the Fisher information computed at the base model where $\xi=0$, represents the variance associated  with the parent removed. The prior is then:

\[
\pi(\xi)=\lambda \sqrt{I(0)} \exp(-\lambda \sqrt{I(0)}\xi).
\]

Hence, the PC-prior is a  tilted  Jeffrey's prior for $\xi$. Thus, the  PC prior for the correlation matrix follows from a numerical approximation to the Hessian matrix of the KLD and
equation <ref>, in one dimension.
We assume a common rate parameter, $\lambda$, for each model pair in the sequence. Different parameters could be considered, although the interpretation and specification by the user would be cumbersome. We advocate using one rate parameter that is common for each step.
The effect of certain choices of $\lambda$ is investigated in the next section.

§.§ Choosing the value for $\lambda$

The prior depends on the  $\lambda$ parameter, which corresponds to the rate of the exponential prior of the pairwise distance between the base model and its corresponding complex model. We proceed numerically to define a reasonable value for $\lambda$ by sampling distances and mapping them to correlations. The most appropriate value for $\lambda$ is the value that contracts correlation towards 0. In Figure  <ref>, we illustrate the range of correlations obtained for different $\lambda$ values and correlation, generated with only one parent and a simpler model made of single children, i.e. independent and no parents.
The histograms show that we reach a reasonable contraction around $\lambda=5$, as too small $\lambda$ values do not contract sufficiently towards zero.

In Figure <ref>, we illustrate the correlations obtained from sampled distances for various $\lambda$ values for a  base model that has 1 parent (as opposed to the base model from Figure <ref>, which had 0 parents), hence a complex model with 2 parents.  In this case, we need to fix the value of the standard deviation of this parent to compute the correlations. We chose three different values for standard deviation $(0.1, 1, 2)$  to illustrate the impact of the choice of $\lambda$.  Still, $\lambda=5$ is a reasonable choice for the prior as it contracts towards 0, without being too restrictive.   The chosen $\lambda $  behaves in a similar way for a correlation derived by graphs with 2 parents and the two simpler models with exchangeable and identity correlation matrices, respectively.

Correlations obtained when sampling distances between the flexible model (1 parent) and the base model (0 parents) for a given lambda.

Correlations obtained when sampling distances between the flexible model (2 parents) and the base model (1 parent) for a given lambda and a given standard deviation of the base model.

Graph sequence that generates correlation matrices, from the complex $a$ to the simplest $h$ (identity correlation matrix), for 8 observed variables or random effects (circle nodes) and 7 latent nodes (square grey nodes).

§ APPLICATIONS

This section provides three applications for the proposed approach for correlation matrices derived from graphs. We show how the methods can be used to identify a potential correlation structure, how it extends to a multivariate conditional autoregressive model for disease mapping and how it applies to longitudinal simulated datasets with multiple random effects.
We define the model, as follows.  For response variables $\pmb y_{ij}, i = 1,2,\dots,n \quad \text{and} \quad j = 1,2,\dots,p$ define the joint distribution through the univariate conditionals as

\[
\pi(\pmb y|\pmb \eta, \mm{\theta}) \sim \prod_{i,j} p(y_{ij}|\eta_{ij}, \mm{\theta})
\]

where the linear predictors $\pmb\eta$ connects the responses $\pmb y$ to the covariates $\pmb Z$ and $\pmb u$ through fixed effects $\pmb\beta$ and random effects $\pmb A$, respectively, as follows:

\begin{equation*}
\pmb \eta = \pmb \beta\pmb Z + \pmb A\pmb u,
\end{equation*}

where we assume Gaussian priors for $\pmb\beta$ and $\pmb u$ with sparse precision matrices. Consider the stacked vector $\mm{x} = (\pmb\beta, \pmb u)$ as the random latent field, then $\pi(\mm{x})$ is a multivariate Gaussian density with a precision matrix composed by the structure of the various effects.  All the models were estimated with R-INLA<cit.>.

§.§ Simple multivariate model

We illustrate the simple case of multivariate penalized prior by analysing a subset of data from the Irish Longitudinal Study on Ageing (TILDA) cohort <cit.>.  The TILDA cohort investigates the physical, mental health, and cognitive measures of Irish residents over 50 years old <cit.>.
The multivariate outcome is given by four epigenetic biomarkers:  Horvath, Hannum, PhenoAge and GrimAge.  The epigenetic clock measures an individual's chronological and biological age. Individuals with good lifestyles (no smoking, limited drinking, active lifestyle), are usually genetically younger than their chronological age, and vice versa. Horvath and Hannum are described as the “first generation clock” as both have been trained on blood samples and aimed at clocks to predict human chronological age. In contrast,  PhenoAge and GrimAge are known as the “second generation clock”, as both were trained using other biomarkers and mortality data to reflect ageing-related physiological conditions.

[shorten >=1pt,node distance=2.5cm,on grid,auto]
state=[shape=circle,thick,draw,minimum size=0.6cm]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
[statesq,fill=rgb:black,1;white,2] (p1) $p_1$;
[statesq,fill=rgb:black,1;white,2, below left of=p1, yshift=5mm] (p2) $p_2$;
[statesq,fill=rgb:black,1;white,2, below right of=p1, yshift=5mm] (p3) $p_3$;
[state, below left of=p2, yshift=5mm] (c1) $c_1$;
[state, below of=p2, yshift=12mm] (c2) $c_2$;
[state,below of=p3, yshift=12mm] (c3) $c_3$;
[state, below right of=p3, yshift=5mm] (c4) $c_4$;
(p1) edge node  (p2)
(p1) edge node  (p3)
(p3) edge node  (c3)
(p3) edge node  (c4)
(p2) edge node  (c1)
(p2) edge node  (c2);

Graph for the model fitted on epigenetic data in Section <ref> with 3 parents and 4 children.

Based on this prior knowledge, we can design the correlation structure of these 4 outcomes using a graph as described in Figure <ref>. The first-generation clock Horvath (children 1) and Hannum (children 2) share some variability captured by their common ancestors' parent 2 ($p_2$) and parent 1 ($p_1$), while the second-generation clock shares some variability captured by their common ancestors' parent 3 ($p_3$) and parent 1. The correlation between first and second-generation clocks is constrained to be lower because they only share parent 1 as a common ancestor. This structure reduces the number of parameters to be estimated from 6 correlations to only 3 compared to an unspecified covariance structure, as there is only one correlation parameter between the two groups of children. As described in Section <ref>, the prior is based on the sequential distance from the complex model described in Figure <ref> until reaching the model, assuming all children are independent. The observed biomarkers show a correlation of $0.85$ and $0.49$ for first and second generation, respectively (see Table <ref>).
We  fit a simple model with no covariates to compare with the correlations provided by our model:

\[
y_{ij}=\alpha+x_{ij}, \quad i=1,...,215, j=1,...,4
\]

Empirical epigenetic biomarkers correlation.

Horvath     Hannum     PhenoAge     GrimAge

Horvath     1

Hannum     0.85     1

PhenoAge     0.08     0.12     1

GrimAge     0.05     0.14     0.49     1

The estimated correlations are reported in Table <ref>, and they are close to the empirical correlations computed from the raw data.

Posterior mean and 95% credible interval for the correlations between the epigenetic clocks.

0.848  [0.807 - 0.882]     Hannum

2c:2*$0.170\: [0.092 - 0.285]$     PhenoAge

0.446 [0.371 - 0.572]     GrimAge

§.§ Multivariate disease mapping

Our second application is the multivariate disease mapping of male lung, oral, oesophagus and larynx cancer mortality data in 544 districts of Germany from 1986 to 1990 <cit.>.
We define a correlation structure based on graphs to model the dependence across diseases in a Multivariate Conditional autoregressive model (MCAR).

The diseases $y_{ij}$  where $j=1,\dots,j$ are distributed as counts over the area  $i=1,\dots,n$. We use a  Poisson  distribution with mean $\theta_{ij}= E_{ij} r_{ij}$, where $E_{ij}$ is the expected number of cases (computed based on the demographic characteristics of a reference population) and $r_{ij}$ is the relative risk:

\[
y_{ij} |E_{ij} r_{ij} \sim \textrm{Pois}(E_{ij} r_{ij}).
\]

The multivariate disease model that smooths the risks is  modelled as follows:

\begin{equation*}
% \label{eqn:eqn1}
\log(r_{ij})=\alpha_j + x_{ij}
\end{equation*}

where $\alpha_j$ is the log baseline risk for the $j$-th disease,
and the  term $x_{ij}$ is a random effect to capture
the variability in $i$ region for $j$-th disease. Thus, let   $\mm X$ be the  random effect matrix, with $X_{i\mathord{\cdot}}, (i=1\dots,n)$ indicating the  $i-th$ row, and $X_{\mathord{\cdot}j}, (j=1,\dots,J)$ the $j-th$ column, the $\text{Vec}({\mm X})=(X_{\mathord{\cdot}1}^T, \cdots, X_{\mathord{\cdot}j}^T)^T$. The $\text{vec(\mm)}$ is defined using a Normal
distribution with zero mean and a highly structured precision matrix $\mm Q$.  This representation is  a multivariate Gaussian vector, with a sparse precision matrix

\begin{equation}
\label{eq:eqmcar1}
\text{Vec}(\mm X) \sim \mathcal{N}(0, \mm Q)
% \pi(\mm{x}|\mm  ) \propto (\vert \mm Q \vert^*)^{1/2}\exp \left(-\frac{1}{2}\tilde{\mm{x}}^T(\mm{Q}_{dis} \otimes \mm{R}_{sp})\tilde{\mm{x}}\right)
\end{equation}

where $\mm Q =\mm{R}_{dis} \otimes \mm{R}_{sp}$, are the
$j \times j$  disease inverse correlation and the $n \times n$ spatial correlation matrix respectively, and  $\otimes$ indicating the Kronecker product.
For spatial component represented by $\mm R_{sp}$ in equation  <ref>,
we considered a precision matrix structure from an
intrinsic CAR model<cit.>.
The precision structure  is usually defined as

\[
\mm R_{sp_{ig}} = \left\{ \begin{array}{rl}
n_i & \textrm{if } i = g \\
-1 & \textrm{if } i\sim g \textrm{ (area $i$ is neighbor to area $g$)} \\
0 & \textrm{otherwise.}
\end{array}
\right.
\]

Since this gives a marginal variance that depends on the neighbourhood
structure, we scaled this matrix, as described by <cit.>.

In contrast to the usual approach that leaves the disease correlation to be revealed only after the model fit, in our
approach, $\mm{R}_{dis}$ is defined by a graph for the four types of cancer based on previous work in the field <cit.>. In Figure <ref>, we depict the expected relationships, knowing that Oesophagus (Osph), Larynx(Lary), and Oral cancer exhibit higher correlations, indicative of shared risk factors or underlying biological mechanisms. Lung cancer, however, has been observed to have a lower correlation with the other three cancer types.  We assume parent 1 ($p_1$) as the common ancestor to all children (i.e., cancer types), including Lung ($c_4$) as its direct children. Larynx ($c_3$), Oesophagus ($c_2$) and Oral ($c_1$) share an additional common ancestor ($p_2$), making the correlation between them higher compared to their correlation with Lung cancer. Finally, Oral cancer has an additional latent factor ($p_3$), allowing for a lower correlation with Lung cancer compared to the other two cancer types as observed in <cit.>. We also admit a simplified graph that allows Oral cancer ($c_1$) to have the same correlation to Lung cancer as Larynx ($c_3$) and Oesophagus ($c_2$), by removing parent 3.

[shorten >=1pt,node distance=2.5cm,on grid,auto]
state=[shape=circle,thick,draw,minimum size=0.6cm]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
[statesq,fill=rgb:black,1;white,2] (p1) $p_1$;
[statesq,fill=rgb:black,1;white,2, below left of=p1, yshift=10mm] (p2) $p_2$;
[statesq,fill=rgb:black,1;white,2, below left of=p2, yshift=10mm] (p3) $p_3$;
[state, below left of=p3, yshift=5mm] (c1) $c_1$;
[state, below of=p2, yshift=5mm] (c2) $c_2$;
[state, below right of=p2, yshift=5mm] (c3) $c_3$;
[state, below right of=p1, yshift=5mm] (c4) $c_4$;
(p3) edge node  (c1)
(p1) edge node  (p2)
(p1) edge node  (c4)
(p2) edge node  (p3)
(p2) edge node  (c2)
(p2) edge node  (c3);

Graph to model the correlation between the four disease risks for Model D.

We fitted five models with various levels of flexibility to illustrate how our new approach offers a trade-off between the model assuming independence and the model with full correlations, both in terms of flexibility and computational burden, see Table <ref>.
First, we fitted the Besag-York-Mollié model <cit.>
for each vector in $\mm x$ corresponding to each one of the diseases,  independently, as model $M_A$.
In the second model $M_B$, the disease part is modelled with an unstructured correlation matrix, for which we used a PC prior for
the precision parameters and LKJ prior <cit.> for the correlation parameters with $\eta=1$, and an independent identically distributed random effect for the spatial areas. While model $M_C$ is the same, except the spatial part is an MCAR. The last two models, $M_D, M_E$  present a correlation derived from the graph for the disease and an MCAR for the spatial part, respectively.
These  models assume  the same spatial smoothing
for each disease and that the correlation among diseases
is the same in the different parts of the study area.
Finally,  to  evaluate and compare the five fitted models, we  reported for each the  Deviance Information Criteria (DIC), the Widely applicable Bayesian information criterion (WAIC), and the Conditional predictive ordinate using either leave-one-out (CPO) or leave-group-out (GCPO) cross-validation,
with level set $m=5$, which means that at least $5$ data points  are removed to cross-validate each observation,
see <cit.> for details.

Description of the fitted models.

Model     Disease     Spatial     Parameters

$M_A$     Independent     Correlated (BYM)     8

$M_B$    Unstructured correlation     Independent     10

$M_C$      Unstructured correlation     Correlated (Besag)     10

$M_D$     Graph model as in Figure <ref>     Correlated (Besag)     7

$M_E$     Simplified graph model (2 parents)     Correlated (Besag)     6

Overall, the models considering spatial and between-diseases correlation ($M_C$, $M_D$, $M_E$) performed better
compared to models assuming independence between diseases ($M_A$)
or spatial independence ($M_B$).  Furthermore, models based on our approach ($M_D$ and $M_E$) performed well compared to the model with the highest flexibility ($M_C$).
On the performance indexes, there are minor differences in the DIC, WAIC, CPO and GCPO
between $M_D$, $M_E$ and $M_C$ with no clear model outperforming the others. Our approach, hence, manages to provide a fit similar to the most complex model at a lower cost. While the difference in the number of parameters is limited, this difference increases a lot when the number of diseases to model increases, making our approach scalable where the flexible model with unstructured correlation quickly reaches limitations in terms of the computational burden.

The DIC, WAIC, CPO and GCPO with the five models fitted.

Model      DIC     WAIC     CPO     GCPO

$M_A$    13811.73     13736.50     7057.66     7033.07

$M_B$    14042.71     13927.67     7461.09     7178.57

$M_C$     13687.23     13643.21     6946.31     6868.24

$M_D$     13688.10     13643.03     6944.93     6888.27

$M_E$    13689.65     13638.85     6945.60     6885.80

Figure <ref> displays the four maps of the raw observed and
estimated SMR based on $M_D$.
As expected, the estimated SMR maps are a smoothed version of the observed ones
highlighting the spatial pattern, which is similar among Oral, Oesophagus and Larynx
cancer risk and different for Lung cancer.

Observed SMR (top maps) and estimated SMR from $M_D$ (bottom maps).

The estimated correlation between the log risk of each disease considering $M_D$ is shown in Table <ref>. The highest correlation is between the Oesophagus and the Larynx, similar to Oral- Oesophagus and Oral-Larynx, with the lowest observed for  Oral-Lung, in line with the graph in Figure <ref>.

Posterior mean and 95% credible interval for the correlations between the logarithm of the relative risk.

2*0.73 (0.59 - 0.82)    Oesophagus

0.74 (0.62 - 0.83)    Larynx

0.49 (0.37 - 0.58)     2c0.49 (0.38 - 0.59)     Lung

The simplified model $M_E$ assumes only two distinct correlations:
(Oral, Oesophagus), (Oral, Larynx) and (Oesophagus, Larynx)
share the same correlation of 0.65 (0.51 - 0.76),
and the correlation between Lung and the others 0.50 (0.40, 0.60). From these values and the goodness-of-fit metrics, one can see that the simplified model returns a similar fit and similar correlations at a lower complexity.

§.§ Multivariate longitudinal analysis

During clinical trials, a group of patients undergoes regular visits, during which data such as blood-measured markers, various tests, and patient-reported outcomes are systematically collected. Our focus lies in modelling several longitudinal markers of interest. Each individual exhibits a distinct deviation from the mean distribution of these markers (random effects), and we posit that there may be potential correlations among the longitudinal markers. Typically, a recognized hierarchical structure exists for the correlation between measurements, such as repeated measurements within an individual or a subset of individuals. Consequently, describing the correlation structure of random effects becomes straightforward through a graphical approach.
Let $y_{ijk}$ denote the value of longitudinal marker $k^{th} (k=1,2)$ for individual $i^{th} (i=1, ..., N_k)$ measured at time points $t_{ijk}$ with $j^{th}$ occasion $(j=1, ..., n_{ik})$. We can describe the distribution of this marker with a mixed effects model:

\[
\textrm{E}(y_{ijk}|\ldots) = \eta_{ik}(t_{ijk}) =\boldsymbol{X}_{ik}(t_{ijk})^\top \boldsymbol{\beta}_k + \boldsymbol{Z}_{ik}(t_{ijk})^\top \boldsymbol{b}_{ik},
\]

where $\eta_{ik}(t_{ijk})$ is the linear predictor defined by fixed effects $\boldsymbol{\beta}_k$ of covariates $\boldsymbol{X}_{ik}(t_{ijk})$ and random effects $\boldsymbol{b}_{ik}$ of covariates $\boldsymbol{Z}_{ik}(t_{ijk})$. The residual error is assumed to follow a Gaussian distribution.

For the sake of simplicity, we simulate two datasets assuming  Gaussian data.  Extension to generalized linear modelling follows trivially.  The first data set includes two longitudinal outcomes modelled with fixed and random intercepts and linear slopes. The second data set shows how this approach easily accommodates more complex hierarchical models with fixed and random effects for the intercepts and linear, quadratic, and cubic slopes.

§.§.§ Longitudinal joint  model with linear effect

The model is defined  for Gaussian data and a linear temporal effect as:

\begin{equation*}
\begin{split}
\textrm{E}(y_{ijk} | \ldots) &= \beta_{0k} + b_{0ik} + (\beta_{1k} + b_{1ik})t_{ijk}.
\end{split}
\end{equation*}

Here $\beta_{0k}$ and $\beta_{1k}$ are the fixed intercept and slope for marker $k$, respectively, while $b_{0ik}$ and $b_{1ik}$ are the corresponding random intercept and slope. For the sake of simplicity, the residual error is assumed negligible. The model includes 4 random effects (random intercept and slope for each marker). We assume the structure described in Figure <ref>.

[shorten >=1pt,node distance=2.5cm,on grid,auto]
state=[shape=circle,thick,draw,minimum size=0.6cm]
statesq=[shape=rectangle,thick,draw,minimum size=0.6cm]
[statesq,fill=rgb:black,1;white,2] (p1) $p_1$;
[statesq,fill=rgb:black,1;white,2, below left of=p1, yshift=5mm] (p2) $p_2$;
[statesq,fill=rgb:black,1;white,2, below right of=p1, yshift=5mm] (p3) $p_3$;
[state, below left of=p2, xshift=5mm, yshift=5mm] (c1) $c_1$;
[state, below right of=p2, xshift=-5mm, yshift=5mm] (c2) $c_2$;
[state,below left of=p3, xshift=5mm, yshift=5mm] (c3) $c_3$;
[state, below right of=p3, xshift=-5mm, yshift=5mm] (c4) $c_4$;
(p1) edge node  (p2)
(p1) edge node  (p3)
(p2) edge node  (c1)
(p2) edge node  (c2)
(p3) edge node  (c3)
(p3) edge node  (c4);

Graph for the model fitted on longitudinal data in Section <ref> with 3 parents and 4 children. Children 1 and 2 correspond to random intercept and slope for the first marker, while children 3 and 4 correspond to random intercept and slope for the second marker.

In this graph, $c_1$ and $c_2$ represent the random intercept and slope of the first longitudinal marker, and $c_3$ and $c_4$ represent the random intercept and slope of the second marker. One can deduce a hierarchy in the correlations between these random effects as $c_1$ and $c_2$ share the common ancestors $p_2$ and $p_1$, and similarly, the other effect shares $p_1$ and $p_3$. The correlation between the two subgroups of children is quantified by their common ancestor $p_1$ and is lower than the correlations between random intercepts and slopes. The correlation structure is then defined as follows:

\begin{equation*}
\left.
\begin{NiceArray}{ll:ll}
\sigma_{c1} & \\
\rho_1 & \sigma_{c2} \\
\hdashline
\multicolumn{2}{c}{\multirow{2}{*}{$\rho_3$}} & \sigma_{c3}\\
& & \rho_2 & \sigma_{c4}
\end{NiceArray} \right.
\end{equation*}

We simulated 200 individuals with 11 repeated measurements of the 2 longitudinal markers for each individual. The value of the fixed intercept and slope for each longitudinal marker is defined as $\beta_{0k}=0.2$ and $\beta_{1k}=-0.1$, respectively, while we fixed the correlation between the subgroups of children to 0.9 and the correlation between the two longitudinal markers to 0.8.  The results of the fitted model given this graph structure are presented in Table <ref>.

Simulations results with 3 parents and 4 children.

Parameter:      True value     Mean Estimate     Credible Interval

$\sigma_{c1}$     1     1.009     [0.916 ; 1.106]

$\sigma_{c2}$     0.2     0.201     [0.183 ; 0.220]

$\sigma_{c3}$     0.1     0.101     [0.092 ; 0.112]

$\sigma_{c4}$     0.5     0.495     [0.451 ; 0.544]

$\rho_1$     0.9     0.912     [0.892 ; 0.932]

$\rho_2$     0.9     0.894     [0.869 ; 0.918]

$\rho_3$     0.8     0.798     [0.737 ; 0.831]

$\beta_{01}$     0.2     0.120     [-0.020 ; 0.259]

$\beta_{11}$     -0.1     -0.113     [-0.141 ; -0.086]

$\beta_{02}$     0.2     0.193     [0.179 ; 0.207]

$\beta_{12}$     -0.1     -0.152     [-0.22 ; -0.083]

The model structure has two advantages. Firstly, it provides a balance between an unspecified covariance structure and a fully independent structure by enabling the definition of single correlation parameters between groups of random effects. In the given example, the correlation parameter $\rho_3$ measures the correlation between two sets of correlated random effects, corresponding to the random intercept and slope for each longitudinal marker. An unspecified correlation would necessitate four correlation parameters instead of one, accounting for the relationship between the two longitudinal markers.

Secondly, this structure facilitates a more intuitive interpretation of the correlation between markers, as it is quantified by a single parameter. In contrast, an unspecified structure involves four terms that could have opposite signs, making the interpretation of the correlation between markers challenging.

§.§.§ Longitudinal joint  model with non-linear effect

It is common in medical research to have flexible trajectories to properly capture the complex nature of biological markers dynamics (e.g., including quadratic or cubic slopes in addition to the linear slope, with associated individual deviations captured by random effects), therefore leading to complex models with many random effects when studying the joint distribution of multiple markers.
Suppose we observe two longitudinal markers ($k=2$); we want to define a model for each marker that includes fixed effects for a given binary ($X_{bin} \sim Binomial(p=0.5)$) and a continuous covariate ($X_{con}$ are independent Gaussians with unit mean and variance 0.5). Additionally, each model includes fixed and random effects for the intercept and the linear, quadratic and cubic slopes. The model is defined as:

\begin{equation*}
\begin{split}
y_{ijk} &= \beta_{0k} + b_{0ik} + (\beta_{1k} + b_{1ik})t_{ijk} + (\beta_{2k} + b_{2ik})t_{ijk}^2   \\
&+(\beta_{3k} + b_{3ik})t_{ijk}^3  +
\beta_{4k} X_{bin} + \beta_{5k} X_{con} + \sigma_{\varepsilon k}, \ k=1,2.
\end{split}
\end{equation*}

This model includes eight random effects; given contextual prior knowledge, we can represent the relationship among those random effects with a graphical structure as displayed as graph $a$ in Figure <ref>.
Following the derivation described in Section <ref>, the parents define the correlation between each random effect, giving the following correlation structure:

\begin{equation*}
\left.
\begin{NiceArray}{cc:cc:cc:cc}
\sigma_{c1} & \\
\rho_1 & \sigma_{c2} \\
\hdashline
\multicolumn{2}{c}{\multirow{2}{*}{$\rho_5$}} & \sigma_{c3}\\
& & \rho_2 & \sigma_{c4}\\
\hdashline
\multicolumn{2}{c}{\multirow{2}{*}{$\rho_7$}} & \multicolumn{2}{c}{\multirow{2}{*}{$\rho_9$}} & \sigma_{c5}\\
& & & & \rho_3 & \sigma_{c6}\\
\hdashline
\multicolumn{2}{c}{\multirow{2}{*}{$\rho_8$}} & \multicolumn{2}{c}{\multirow{2}{*}{$\rho_{10}$}} & \multicolumn{2}{c}{\multirow{2}{*}{$\rho_6$}} & \sigma_{c7}\\
& & & & & & \rho_4 & \sigma_{c8}\\
\end{NiceArray} \right.
\end{equation*}

We simulated observations of the two longitudinal biomarkers for 500 individuals, each with 11 repeated measurements; the fitted values are presented in Table <ref>. In this example, only 15 parameters were estimated, compared with the 36 parameters (28 correlation terms and 8 random-effects specific standard deviation terms) given by an unstructured correlation.

Simulations results with 7 parents and 8 children.

Parameter:      True value     Mean Estimate     Credible Interval

$\sigma_{c1}$     1     0.976     [0.921 ; 1.040]

$\sigma_{c2}$     0.2     0.195     [0.184 ; 0.207]

$\sigma_{c3}$     0.1     0.1     [0.094 ; 0.106]

$\sigma_{c4}$     0.5     0.503     [0.476 ; 0.535]

$\sigma_{c5}$     3     3.049     [2.873 ; 3.234]

$\sigma_{c6}$     4     3.932     [3.683 ; 4.178]

$\sigma_{c7}$     0.5     0.526     [0.495 ; 0.560]

$\sigma_{c8}$     0.2     0.208     [0.195 ; 0.222]

$\rho_1$     0.9     0.904     [0.889 ; 0.917]

$\rho_2$     0.9     0.899     [0.883 ; 0.913]

$\rho_3$     0.9     0.888     [0.872 ; 0.904]

$\rho_4$     0.9     0.913     [0.900 ; 0.925]

$\rho_5$     0.8     0.814     [0.785 ; 0.838]

$\rho_6$     0.7     0.716     [0.674 ; 0.751]

$\rho_7$     0.6     0.578     [0.518 ; 0.624]

$\rho_8$     0.6     0.592     [0.533 ; 0.640]

$\rho_9$     0.5     0.510     [0.452 ; 0.560]

$\rho_{10}$     0.5     0.523     [0.468 ; 0.571]

$\beta_{01}$     0.2     0.165     [0.049 ; 0.280]

$\beta_{11}$     -0.1     -0.108     [-0.125 ; -0.091]

$\beta_{21}$     -0.1     -0.100     [-0.109 ; -0.092]

$\beta_{31}$     0.1     0.098     [0.054 ; 0.142]

$\beta_{41}$     -0.2     -0.169     [-0.238 ; -0.099]

$\beta_{51}$     0.1     0.073     [0.001 ; 0.145]

$\beta_{02}$     0.2     -0.028     [-0.401 ; 0.346]

$\beta_{12}$     -0.1     -0.075     [-0.421 ; 0.272]

$\beta_{22}$     -0.1     -0.089     [-0.136 ; -0.042]

$\beta_{32}$     0.1     0.105     [0.087 ; 0.123]

$\beta_{42}$     -0.2     0.084     [-0.152 ; 0.319]

$\beta_{52}$     0.1     0.240     [-0.002 ; 0.483]

§ CONCLUSIONS

Using an intuitive graphical construction,  this paper introduced an approach to defining correlation matrices for observed variables or random effects. The essence of our approach is rooted in the interpretation of the model structure, unlike the conventional methods, which offer an unstructured and often ambiguous interpretation of correlations.  Each element of the correlation matrix is no longer an isolated parameter but a piece of a coherent, inspired narrative that describes the intricate relationships among variables.

This new approach offers several notable advantages. Firstly, it aligns the statistical modelling process more closely with prior knowledge and empirical observations than covariance models. Researchers and practitioners can map their theoretical constructs and hypotheses directly onto the statistical model, facilitating the union of theory and practical analysis. Hence, this approach enhances the interpretability and transparency of the correlation structure.
Every estimated correlation is embedded within a framework that makes the results more accessible and interpretable for a broader audience, including those without advanced statistical training.

Secondly,  our approach improves interpretability and eases computational burden by streamlining the models and decreasing the necessary number of parameters. This introduces a balance between the excessively simplistic independence assumption and the computationally demanding unstructured correlation structures commonly found in conventional practices.

Lastly, the multivariate penalized priors derived are proper by construction, and posterior proprieties are guaranteed.  It also provides a stable choice of prior if no other expert choice is available since the prior ensures contraction to a less complex model if the data lacks support for the complex model.

We have restricted our interpretation to a subclass of directed acyclic graphs to avoid non-suitable correlation matrices. Our starting point was to have a drawing (a graph) based on how researchers hypothesize the multiple variables are related.
We acknowledge that this is a convenient approach and privileges simple graphs, being aware that similar correlation matrices could have been generated by other types of graphs, but for which we could not have a better understanding.

In conclusion, the graphical framework provides interpretable correlation matrices, and the automatic prior construction uses a sequence of simpler models embedded in the penalized complexity prior framework.

§ TECHNICAL DETAILS

§.§ Computation of the PC prior for a complex exchangeable matrix

We derive first the exact PC prior for an exchangeable correlation matrix (one latent factor and 2-children)correlation matrix is $\mm C_1$ from equation <ref>
and the simpler model with an identity matrix $\mm C_0= \mm I$.
The KLD is:

\[
\text{KLD}= \frac{1}{2} \left\{\text{tr}(\mm C_0^{-1}\mm C_1) - p - \log\left(\frac{|\mm C_1|}{|\mm C_0|}\right)\right\}
\]

where $p=2$ (dimension of $\mm C_1$), $\text{tr}(\mm C_0^{-1}\mm C_1)=2$, $|\mm C_0|=1$ and $|\mm C_1|=\frac{2q_1^2+1}{(1+q_1^2)^2}$,

\[
\text{KLD} =\frac{1}{2}\left(2 - 2 - \log\left(\frac{2q_1^2+1}{(1+q_1^2)^2}\right)\right).
\]

Then, the distance is reduced to:

\[
\text{Dist} = \sqrt{2 \cdot \text{KLD}} =\sqrt{-\log\left(\frac{2q_1^2+1}{(1+q_1^2)^2}\right)}
\]

Hence, the prior for $\theta=\log(q_1^2)$ is:

\[
\pi(\theta) = \lambda \exp\left(-\lambda \cdot \sqrt{-\log\left(\frac{2q_1^2+1}{(1+q_1^2)^2}\right)}\right) \left|\frac{\partial \text{Dist}}{\partial \theta}\right|
\]

Using the second order expansion instead, we get $\text{Dist}= q_1^2$, so that $\pi(\theta) = \lambda \exp(-\lambda \exp(\theta)) \exp(\theta)$.

§ ACKNOWLEDGMENTS

AFS is grateful to Prof. Cathal McCrory for TILDA  epigenetic biomarkers data.

§ CONFLICT OF INTEREST

§ FUNDING

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

§ CREDIT AUTHORSHIP CONTRIBUTION STATEMENT

Anna Freni Sterrantino:  Conceptualization,  Methodology, Writing – original draft, Writing – review & editing, Visualization.

Denis Rustand: Software, Methodology, Writing – original draft, Writing – review & editing, Visualization.

Janet van Niekerk: Methodology, Writing – original draft,  Writing – review & editing.

Elias Teixeira Krainski:Software, Methodology, Writing – original draft,Visualization.

Håvard Rue: Conceptualization,  Methodology,  Supervision, Writing – review & editing.