Graph-based multivariate conditional autoregressive models

ABSTRACT

The conditional autoregressive model is a routinely used statistical model for areal data that arise from, for instances, epidemiological, socio-economic or ecological studies. Various multivariate conditional autoregressive models have also been extensively studied in the literature and it has been shown that extending from the univariate case to the multivariate case is not trivial. The difficulties lie in many aspects, including validity, interpretability, flexibility and computational feasibility of the model. In this paper, we approach the multivariate modelling from an element-based perspective instead of the traditional vector-based perspective. We focus on the joint adjacency structure of elements and discuss graphical structures for both the spatial and non-spatial domains. We assume that the graph for the spatial domain is generally known and fixed while the graph for the non-spatial domain can be unknown and random. We propose a very general specification for the multivariate conditional modelling and then focus on three special cases, which are linked to well-known models in the literature. Bayesian inference for parameter learning and graph learning is provided for the focused cases, and finally, an example with public health data is illustrated.

KEYWORDS:

• Areal data
• disease mapping
• graphical model
• G-wishart distribution
• Markov random field
• reversible jump

1. Introduction

Areal data, sometimes called lattice data, are usually represented by an undirected graph where each vertex represents an areal unit and each edge represents a neighbouring relationship. A finite set of random variables on an undirected graph, where each vertex is a random variable, is called a Markov random field if it has the Markov property. Hence, the Markov random field models are often used for the areal data. The univariate conditional autoregressive (CAR) model, originated from Besag (1974), is a Gaussian Markov random field model, for which the joint distribution is multivariate Gaussian. Let $\mathit{u}=(u_1,\dots ,u_I{)}^T$ be a vector of random variables on I areal units (i.e., I vertices). The zero-centred conditional autoregressive model specifies full conditional Gaussian distributions $u_i\mid u_{-i}\sim \mathrm{N}\left(\sum _{i^{\prime }\ne i}{b}_{i{i}^{\prime }}{u}_{i^{\prime }},{\tau }_i^2\right),i=1,\dots ,I,$ where $u_{-i}$ is the collection of $u_{i^{\prime }}$ for $i^{\prime }\ne i$. The resulting joint distribution, derived using Brook's lemma, has a density function as follows, $f(\mathit{u}\mid {\mathit{T}}_{\mathrm{C}\mathrm{A}\mathrm{R}},{\mathit{B}}_{\mathrm{C}\mathrm{A}\mathrm{R}})\propto \exp \left\{-\frac{1}{2}{\mathit{u}}^T{\mathit{T}}_{\mathrm{C}\mathrm{A}\mathrm{R}}^{-1}(\mathit{I}-{\mathit{B}}_{\mathrm{C}\mathrm{A}\mathrm{R}})\mathit{u}\right\},$ where $\mathit{I}$ is an identity matrix; ${\mathit{B}}_{\mathrm{C}\mathrm{A}\mathrm{R}}$ is an $I\times I$ matrix whose off-diagonal entries are $b_{i{i}^{\prime }}$ and diagonal entries are zeros, and ${\mathit{T}}_{\mathrm{C}\mathrm{A}\mathrm{R}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{{\tau }_1^2,\dots ,{\tau }_I^2\}$. The joint distribution is multivariate Gaussian if and only if ${\mathit{T}}_{\mathrm{C}\mathrm{A}\mathrm{R}}^{-1}(\mathit{I}-{\mathit{B}}_{\mathrm{C}\mathrm{A}\mathrm{R}})$ is symmetric and positive definite. A further parameterization on ${\mathit{B}}_{\mathrm{C}\mathrm{A}\mathrm{R}}$ and ${\mathit{T}}_{\mathrm{C}\mathrm{A}\mathrm{R}}$ is needed to reduce the number of parameters in the model. Consider a so-called adjacency matrix ${\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}$ for the undirected graph, where the $i{i}^{\prime }$th entry $C_{i{i}^{\prime }}=1$ if unit i and unit $i^{\prime }$ are neighbours (denoted as $i\sim i^{\prime }$) and $C_{i{i}^{\prime }}=0$ otherwise. One popular parameterization is to let $b_{i{i}^{\prime }}=\rho C_{i{i}^{\prime }}/C_{i+}$ and ${\tau }_i^2={\sigma }^2/C_{i+}$, where $C_{i+}$ is the ith row sum of ${\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}$, representing the total number of neighbours of unit i. Let ${\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{C_{1+},\dots ,C_{I+}\}$. When ρ is strictly between the smallest and largest eigenvalues of ${\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}^{-1/2}{\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}{\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}^{-1/2}$, or sufficiently, when $|\rho |<1$, and ${\sigma }^2>0$, the joint distribution of $\mathit{u}$ is a zero-mean multivariate Gaussian distribution: $\mathit{u}\sim N\{\mathbf{0},{\sigma }^2({\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}-\rho {\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}{)}^{-1}\}$. This is called the proper conditional autoregressive model in the literature. When $\rho =1$, it is called the intrinsic conditional autoregressive model which is an improper distribution due to the singular covariance matrix.

Turning to the multivariate case, consider J responses (e.g., multiple diseases) on I areal units. Let $\mathit{U}$ be an $I\times J$ matrix-variate, where the ijth entry $u_{ij}$ is a random variable for the ith areal unit and jth response. Each column of $\mathit{U}$ is an areal vector for a single response and hence can be modelled by the univariate conditional autoregressive model. However, a multivariate model is desired for the matrix-variate $\mathit{U}$ in order to simultaneously model the dependence across responses. Initially proposed by Mardia (1988), the multivariate conditional autoregressive model specifies full conditional distributions on row vectors of $\mathit{U}$. Let ${\mathit{u}}_i$ be the ith row vector of $\mathit{U}$. Following Besag (1974), specify (1) ${\mathit{u}}_i\mid {\mathit{u}}_{-i}\sim \mathrm{N}\left(\sum _{i^{\prime }\ne i}{\mathit{B}}_{i{i}^{\prime }}{\mathit{u}}_{i^{\prime }},{\mathbf{\Sigma }}_i\right),i=1,\dots ,I,$(1) where ${\mathit{B}}_{i{i}^{\prime }}$ and ${\mathbf{\Sigma }}_i$ are $J\times J$ matrices needing a further parameterization. To make the joint distribution for $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{U}}^T\right)$ a multivariate Gaussian, ${\mathit{B}}_{i{i}^{\prime }}$ and ${\mathbf{\Sigma }}_i$ must satisfy certain conditions (Mardia, 1988). Gelfand and Vounatsou (2003) showed a convenient parameterization, ${\mathit{B}}_{i{i}^{\prime }}=\left(\rho C_{i{i}^{\prime }}/C_{i+}\right){\mathit{I}}_J$ and ${\mathbf{\Sigma }}_i=\mathbf{\Sigma }/C_{i+}$. When $|\rho |<1$ and $\mathbf{\Sigma }$ is positive definite, $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{U}}^T\right)$ has a zero-mean multivariate Gaussian distribution: $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{U}}^T\right)\sim N\{\mathbf{0},({\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}-\rho {\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}{)}^{-1}\otimes \mathbf{\Sigma }\}$. It is clear that this multivariate specification is a Kronecker-product formula where $({\mathit{D}}_{\mathrm{C}\mathrm{A}\mathrm{R}}-\rho {\mathit{C}}_{\mathrm{C}\mathrm{A}\mathrm{R}}{)}^{-1}$ models the covariance structure across rows of $\mathit{U}$ (spatial domain) and $\mathbf{\Sigma }$ models the covariance structure across columns of $\mathit{U}$ (response domain). From the modelling perspective, Mardia's specification has a difficulty with parameterization. It is usually difficult to have a meaningful parameterization for ${\mathit{B}}_{i{i}^{\prime }}$ and ${\mathbf{\Sigma }}_i$ unless one pursues a simple formulation. It is arguable that the Mardia's specification presents a conflict, where the between vector variation is specified through an inverse covariance matrix, but the within vector variation is specified through a covariance matrix. It seems more intuitive to either work with the joint covariance or the joint inverse covariance directly. Notice that most multivariate spatial models for point reference data focus on the joint covariance structure. In this paper, we focus on the joint inverse covariance structure of elements in the multivariate areal data. In particular, we consider the joint adjacency structure of the lattice based on graphical structures of both the spatial domain and the response domain. We build a framework for graph-based multivariate conditional autoregressive models and discuss parameterizations under this framework. The advantage is that this framework is very general and we demonstrate it through multiple case examples. Furthermore, we allow graph learning for multiple responses in such models, which is potentially useful for many modern applications.

We shall point out other recent work on multivariate conditional autoregressive models. Kim, Sun, and Tsutakawa (2001) and Jin, Carlin, and Banerjee (2005) proposed conditional autoregressive models for bivariate areal data. Multivariate models were considered by Gelfand and Vounatsou (2003), Jin, Banerjee, and Carlin (2007), MacNab (2011, 2016), Martinez-Beneito (2013), Martinez-Beneito, Botella-Rocamora, and Banerjee (2017) among many others. MacNab (2018) reviewed some recent developments on multivariate Gaussian Markov random field models. We will show that some of the earlier work can be reconstructed in our proposed framework and some can be extended to graphical models. The paper is organised as follows. Section 2 presents the general framework and three special parameterizations. Section 3 presents a real data example using the proposed models. Section 5 contains further discussions and remarks. Technical details are given in the Appendices.

2. Graph-based multivariate conditional autoregressive models

2.1. General framework

Instead of specifying full conditional distributions on vectors like (1(1) ${\mathit{u}}_i\mid {\mathit{u}}_{-i}\sim \mathrm{N}\left(\sum _{i^{\prime }\ne i}{\mathit{B}}_{i{i}^{\prime }}{\mathit{u}}_{i^{\prime }},{\mathbf{\Sigma }}_i\right),i=1,\dots ,I,$(1) ), we approach this problem from an element-based perspective. Following Besag (1974), specify full conditional distributions for each element $u_{ij}$ in the matrix-variate $\mathit{U}$ as follows, $\begin{array}{rl}& u_{ij}\mid u_{-\left\{ij\right\}}\sim \mathrm{N}\left(\sum _{\left\{i^{\prime }{j}^{\prime }\right\}\ne \left\{ij\right\}}{b}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}{u}_{i^{\prime }{j}^{\prime }},{\tau }_{ij}^2\right),\\ & i=1,\dots ,I\mathrm{a}\mathrm{n}\mathrm{d}j=1,\dots ,J,\end{array}$ where $\left\{i^{\prime }{j}^{\prime }\right\}\ne \left\{ij\right\}$ means either $i^{\prime }\ne i$ or $j^{\prime }\ne j$. In fact, here we consider a lattice consisting of all elements in $\mathit{U}$. Using Brook's lemma, the resulting joint distribution for $\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}\right)$ is $\begin{array}{rl}& f\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(\mathit{U})\mid \mathit{B},\mathit{T})\propto \exp \\ & \times \left\{-\frac{1}{2}\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}{)}^T{\mathit{T}}^{-1}\right(\mathit{I}-\mathit{B}\left)\mathrm{v}\mathrm{e}\mathrm{c}\right(\mathit{U})\right\},\end{array}$ where $\mathit{I}$ is an $IJ\times IJ$ identity matrix, $\mathit{T}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{{\tau }_{11}^2,\dots ,{\tau }_{I1}^2,\dots ,{\tau }_{1J}^2,\dots ,{\tau }_{IJ}^2\}$ and $\mathit{B}$ can be expressed block-wisely, $\begin{array}{rl}\mathit{B}& =\left(\begin{array}{ccc}{\mathit{B}}_{11}& \cdots & {\mathit{B}}_{1J}\\ ⋮& \ddots & ⋮\\ {\mathit{B}}_{J1}& \cdots & {\mathit{B}}_{JJ}\end{array}\right)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ {\mathit{B}}_{j{j}^{\prime }}& =\left(\begin{array}{ccc}{b}_{\left\{1j\right\},\left\{1j^{\prime }\right\}}& \cdots & b_{\left\{1j\right\},\left\{I{j}^{\prime }\right\}}\\ ⋮& \ddots & ⋮\\ b_{\left\{Ij\right\},\left\{1j^{\prime }\right\}}& \cdots & b_{\left\{Ij\right\},\left\{I{j}^{\prime }\right\}}\end{array}\right),\end{array}$ and the diagonal elements $b_{\left\{ij\right\},\left\{ij\right\}}$ are zeros. The joint distribution for $\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}\right)$ is multivariate Gaussian if and only if ${\mathit{T}}^{-1}(\mathit{I}-\mathit{B})$ is symmetric and positive definite. It is desired that $\mathit{B}$ and $\mathit{T}$ are further parameterised to reduce the number of parameters in the model. We denote this general model $\mathrm{M}\mathrm{C}\mathrm{A}\mathrm{R}(\mathit{B},\mathit{T})$ for later use.

Consider the adjacency structure of the undirected graph for the lattice of $\mathit{U}$. In the univariate situation, the adjacency structure is determined by those geographical locations. Two areal units are connected by an edge if they are neighbours geographically. However, it is not obvious which elements should be neighbours in $\mathit{U}$. Consider that the J responses can be connected through an undirected graph. Let ${\mathit{C}}^{\left(s\right)}$ be the adjacency matrix for all I areal units and ${\mathit{C}}^{\left(r\right)}$ be the adjacency matrix for all J responses. Both the spatial graph and the response graph are then uniquely determined by ${\mathit{C}}^{\left(s\right)}$ and ${\mathit{C}}^{\left(r\right)}$, respectively. Let $\mathit{C}$ be the joint adjacency matrix for the lattice of $\mathit{U}$. A general construction of $\mathit{C}$ can be made through ${\mathit{C}}^{\left(s\right)}$ and ${\mathit{C}}^{\left(r\right)}$, (2) $\mathit{C}={\mathit{C}}^{\left(r\right)}\otimes {\mathit{C}}^{\left(s\right)}+{\mathit{C}}^{\left(r\right)}\otimes {\mathit{I}}_I+{\mathit{I}}_J\otimes {\mathit{C}}^{\left(s\right)}.$(2) This construction connects $u_{ij}$ with $u_{i^{\prime }\sim i,j}$, $u_{i,j^{\prime }\sim j}$ and $u_{i^{\prime }\sim i,j^{\prime }\sim j}$, meaning its spatial neighbour, response neighbour and interaction neighbour, respectively. One may add edges for secondary neighbours or drop edges in a specific modelling. For example, some reduced constructions would be (i) $\mathit{C}={\mathit{I}}_J\otimes {\mathit{C}}^{\left(s\right)}$ (independent conditional autoregressive models, no dependence between responses); (ii) $\mathit{C}={\mathit{C}}^{\left(r\right)}\otimes {\mathit{I}}_I$ (independent multivariate variables, no spatial dependence); and (iii) $\mathit{C}={\mathit{C}}^{\left(r\right)}\otimes {\mathit{I}}_I+{\mathit{I}}_J\otimes {\mathit{C}}^{\left(s\right)}$ (drop edges for interaction neighbours $u_{i^{\prime }\sim i,j^{\prime }\sim j}$).

Let $C_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}$ denote entries in $\mathit{C}$, analogous to the block-wise notation $b_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}$ for $\mathit{B}$. Let $d_j^{\left(r\right)}$ be the jth row sum in ${\mathit{C}}^{\left(r\right)}$ and $d_i^{\left(s\right)}$ be the ith row sum in ${\mathit{C}}^{\left(s\right)}$. Then the ijth row sum in $\mathit{C}$ is $d_{ij}=d_j^{\left(r\right)}{d}_i^{\left(s\right)}+d_j^{\left(r\right)}+d_i^{\left(s\right)}$. Let ${\mathit{D}}^{\left(r\right)}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{d_1^{\left(r\right)},\dots ,d_J^{\left(r\right)}\}$, ${\mathit{D}}^{\left(s\right)}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{d_1^{\left(s\right)},\dots ,d_I^{\left(s\right)}\}$ and $\mathit{D}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{d_{11},\dots ,d_{I1},\dots ,d_{1J},\dots ,d_{IJ}\}$. With the adjacency constructions and notations, we then explore further parameterisation on $\mathit{B}$ and $\mathit{T}$ in the following subsections, and specifically, we discuss three specifications made from this general framework, all of which are linked to well-known models in the literature.

2.2. Model 1: nonseparable multifold specification

Kim et al. (2001) developed a twofold conditional autoregressive model for bivariate areal data (J = 2), using different linkage parameters for different types of neighbours. Those linkage parameters, in their work, are called smoothing and bridging parameters, representing the strength of information sharing. If we extend their specification to an arbitrary J, we can parameterise $\mathit{B}$ and $\mathit{T}$ in the following way (assuming $i\ne i^{\prime }$ and $j\ne j^{\prime }$): $\begin{array}{rl}& b_{\left\{ij\right\},\left\{i^{\prime }j\right\}}=\frac{{\lambda }_j}{d_{ij}}{C}_{\left\{ij\right\},\left\{i^{\prime }j\right\}},b_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}}=\frac{{\psi }_{j{j}^{\prime }}}{d_{ij}}\sqrt{\frac{{\delta }_j}{{\delta }_{j^{\prime }}}}{C}_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}},\\ & b_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}=\frac{{\phi }_{j{j}^{\prime }}}{d_{ij}}\sqrt{\frac{{\delta }_j}{{\delta }_{j^{\prime }}}}{C}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}},{\tau }_{ij}^2=\frac{{\delta }_j}{d_{ij}},\end{array}$ where ${\lambda }_j$, ${\psi }_{j{j}^{\prime }}$ and ${\phi }_{j{j}^{\prime }}$ are linkage parameters and ${\delta }_j$ are variance components. Linkage parameters are for three types of neighbour: $u_{i^{\prime }\sim i,j}$, $u_{i,j^{\prime }\sim j}$ and $u_{i^{\prime }\sim i,j^{\prime }\sim j}$. Having this specification, the conditional mean of $u_{ij}$ essentially is $\begin{array}{rl}\mathrm{E}(u_{ij}\mid u_{-\left\{ij\right\}})& =\frac{1}{d_{ij}}\left({\lambda }_j\sum _{i^{\prime }\sim i}{u}_{i^{\prime }j}+{\psi }_{j{j}^{\prime }}\sum _{j^{\prime }\sim j}\sqrt{\frac{{\delta }_j}{{\delta }_{j^{\prime }}}}{u}_{i{j}^{\prime }}\right.\\ & +\left.{\phi }_{j{j}^{\prime }}\sum _{i^{\prime }\sim i}\sum _{j^{\prime }\sim j}\sqrt{\frac{{\delta }_j}{{\delta }_{j^{\prime }}}}{u}_{i^{\prime }{j}^{\prime }}\right),\end{array}$ which is a weighted average of all its neighbours in $\mathit{C}$. This specification generalises (Kim et al., 2001)' twofold model and hence could be called a multifold specification. It can be shown that, for this parameterisation, the joint precision matrix is (3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) where $\mathbf{\Delta }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{{\delta }_1,\dots ,{\delta }_J\}$, $\mathbf{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{{\lambda }_1,\dots ,{\lambda }_J\}$, $\mathbf{\Psi }$ and $\mathbf{\Phi }$ are $J\times J$ symmetric matrices with entries ${\psi }_{j{j}^{\prime }}$ and ${\phi }_{j{j}^{\prime }}$, respectively, and the operator ° means an element-wise product. A derivation of (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ) is given in Appendix 1. Note that only nonzero entries of $\mathbf{\Psi }$ and $\mathbf{\Phi }$ are parameters in the model, the number of which depends on ${\mathit{C}}^{\left(r\right)}$.

In order to make (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ) positive definite, constraints on ${\lambda }_j$, ${\psi }_{j{j}^{\prime }}$ and ${\phi }_{j{j}^{\prime }}$ are needed, assuming that ${\delta }_j>0$. In general, it is difficult to find a sufficient and necessary condition for the positive definiteness of (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ). Kim et al. (2001)'s solution to this problem was a sufficient condition: $\text{max}\{|{\lambda }_j|,|{\psi }_{j{j}^{\prime }}|,|{\phi }_{j{j}^{\prime }}|;\mathrm{\forall }j,j^{\prime }\}<1$, under which the matrix (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ) is diagonally dominant and hence is positive definite. Though their proof was under J = 2, it is true for any J by the same arguments. The advantage of this condition is that it is simple and implementable. However, this is not a necessary condition meaning that it is impossible to reach all possible positive definite structures for the model under such a condition. In a Bayesian model, priors on parameters ${\lambda }_j$, ${\psi }_{j{j}^{\prime }}$ and ${\phi }_{j{j}^{\prime }}$ can be chosen based on their actual constraints. In our case, a uniform prior $\text{Unif}(-1,1)$ is adequate for these linkage parameters. Priors on the variance components ${\delta }_j$ can be weakly-informative inverse gamma priors $\text{IG}(a_j,b_j)$. Inference and computation under this model are given in Appendix 2.

2.3. Model 2: separable specification with homogeneous spatial smoothing

Gelfand and Vounatsou (2003)'s Kronecker-product model is a convenient parameterisation of Mardia (1988)'s model. In our framework, this specification can be obtained and extended by having the following parameterisation for $\mathit{B}$ and $\mathit{T}$ (assuming $i\ne i^{\prime }$ and $j\ne j^{\prime }$): $\begin{array}{rl}& b_{\left\{ij\right\},\left\{i^{\prime }j\right\}}=\frac{\rho }{d_i^{\left(s\right)}}{C}_{\left\{ij\right\},\left\{i^{\prime }j\right\}},b_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}}=-\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}{C}_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}},\\ & b_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}=\frac{\rho {\omega }_{j{j}^{\prime }}}{d_i^{\left(s\right)}{\omega }_{jj}}{C}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}},{\tau }_{ij}^2=\frac{1}{d_i^{\left(s\right)}{\omega }_{jj}},\end{array}$ where ρ and ${\omega }_{j{j}^{\prime }}$ are linkage parameters, and $1/{\omega }_{jj}$ are variance components. This parameterisation does not seem straightforward, but is much clearer in the form of conditional mean: (4) $\begin{array}{rl}& \mathrm{E}\left(u_{ij}-\frac{\rho }{d_i^{\left(s\right)}}\sum _{i^{\prime }\sim i}{u}_{i^{\prime }j}|u_{-\left\{ij\right\}}\right)\\ & =-\sum _{j^{\prime }\sim j}\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}\left(u_{i{j}^{\prime }}-\frac{\rho }{d_i^{\left(s\right)}}\sum _{i^{\prime }\sim i}{u}_{i^{\prime }{j}^{\prime }}\right).\end{array}$(4) Note that for a single response, the univariate conditional autoregressive model specifies $\text{E}\left(u_i|u_{-i}\right)=\rho \sum _{i^{\prime }\sim i}{u}_{i^{\prime }}/d_i^{\left(s\right)}$. In the multivariate setting, $\rho \sum _{i^{\prime }\sim i}{u}_{i^{\prime }j}/d_i^{\left(s\right)}$ is no longer the conditional mean for $u_{ij}\mid u_{-\left\{ij\right\}}$ and their conditional difference is regressed on other differences through ${\omega }_{j{j}^{\prime }}$. This parameterisation yields the joint precision matrix (5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) where $\mathbf{\Omega }$ is a symmetric $J\times J$ matrix with entries ${\omega }_{j{j}^{\prime }}$. A derivation of (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ) is given in Appendix 1. The linkage parameter ρ is interpreted as a spatial smoothing parameter and $\mathbf{\Omega }$ controls the dependence across J responses. It is noteworthy that only nonzero entries in $\mathbf{\Omega }$ are parameters in the model and we denote ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}=\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})$ for simplicity. The notation ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$, commonly used in graphical models, means the precision matrix restricted by graph ${\mathit{C}}^{\left(r\right)}$. The model (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ) is a natural extension of Gelfand and Vounatsou (2003)'s model. When ${\mathit{C}}^{\left(r\right)}$ is the complete graph (any two vertices are connected), ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ is free of zero entries. Then let $\mathbf{\Sigma }={\mathbf{\Omega }}^{-1}$ and (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ) is equivalent to Gelfand and Vounatsou (2003)'s specification. We call this a completely separable specification because the Kronecker product completely separates the spatial domain and the response domain. This complete separation is often not desirable because it makes the spatial smoothing common for all j. We call this homogeneous spatial smoothing because the linkage ρ is the same for any i and $i^{\prime }$ which distinguishes Model 2 from Model 3 in the next section.

The joint precision matrix (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ) is positive definite if $|\rho |<1$ and ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ is positive definite. Let $M^{+}\left({\mathit{C}}^{\left(r\right)}\right)$ be the cone of symmetric positive definite matrices restricted by ${\mathit{C}}^{\left(r\right)}$ and then ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\in M^{+}\left({\mathit{C}}^{\left(r\right)}\right)$. In a Bayesian model, a widely used prior on ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ is the G-Wishart distribution (Atay-Kayis & Massam, 2005; Letac & Massam, 2007). The G-Wishart distribution is a conjugate family for the precision matrix of a Gaussian graphical model, whose density function is given by $\begin{array}{rl}& p({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid b,\mathit{V})=I_{{\mathit{C}}^{\left(r\right)}}(b,\mathit{V}{)}^{-1}{\left|{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\right|}^{\left(b-2\right)/2}\exp \\ & \times \left\{-\frac{1}{2}\mathrm{t}\mathrm{r}\left(\mathit{V}{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\right)\right\}{1}_{{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\in M^{+}\left({\mathit{C}}^{\left(r\right)}\right)},\end{array}$ where b>2 is the number of degrees of freedom; $\mathit{V}$ is the scale matrix and $I_{{\mathit{C}}^{\left(r\right)}}(\cdot )$ is the normalising constant. It is practically attractive because of its conjugacy. That said, for a prior distribution $\text{GWis}(b,\mathit{V})$ and a given sample covariance matrix $\mathit{S}$ of sample size n, the posterior distribution of ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ is $\text{GWis}(b+n,\mathit{V}+\mathit{S})$. Inference and computation under this model are given in Appendix 2.

2.4. Model 3: separable specification with heterogenous spatial smoothing

Dobra, Lenkoski, and Rodriguez (2011) introduced a multivariate lattice model by giving Kronecker-product G-Wishart priors to the matrix-variate $\mathit{U}$. In our framework, $\mathit{B}$ and $\mathit{T}$ can be parameterised in the following way (assuming $i\ne i^{\prime }$ and $j\ne j^{\prime }$): $\begin{array}{rl}& b_{\left\{ij\right\},\left\{i^{\prime }j\right\}}=-\frac{{\omega }_{i{i}^{\prime }}^{\left(s\right)}}{{\omega }_{ii}^{\left(s\right)}}{C}_{\left\{ij\right\},\left\{i^{\prime }j\right\}},b_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}}=-\frac{{\omega }_{j{j}^{\prime }}^{\left(r\right)}}{{\omega }_{jj}^{\left(r\right)}}{C}_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}},\\ & b_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}=-\frac{{\omega }_{i{i}^{\prime }}^{\left(s\right)}{\omega }_{j{j}^{\prime }}^{\left(r\right)}}{{\omega }_{ii}^{\left(s\right)}{\omega }_{jj}^{\left(r\right)}}{C}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}},{\tau }_{ij}^2=\frac{1}{{\omega }_{ii}^{\left(s\right)}{\omega }_{jj}^{\left(r\right)}},\end{array}$ which is equivalent to the version of conditional mean (6) $\begin{array}{rl}& \mathrm{E}\left(u_{ij}-\frac{1}{{\omega }_{ii}^{\left(s\right)}}\sum _{i^{\prime }\sim i}{\omega }_{i{i}^{\prime }}^{\left(s\right)}{u}_{i^{\prime }j}|u_{-\left\{ij\right\}}\right)\\ & =-\sum _{j^{\prime }\sim j}\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}\left(u_{i{j}^{\prime }}-\frac{1}{{\omega }_{ii}^{\left(s\right)}}\sum _{i^{\prime }\sim i}{\omega }_{i{i}^{\prime }}^{\left(s\right)}{u}_{i^{\prime }{j}^{\prime }}\right).\end{array}$(6) Comparing (6(6) $\begin{array}{rl}& \mathrm{E}\left(u_{ij}-\frac{1}{{\omega }_{ii}^{\left(s\right)}}\sum _{i^{\prime }\sim i}{\omega }_{i{i}^{\prime }}^{\left(s\right)}{u}_{i^{\prime }j}|u_{-\left\{ij\right\}}\right)\\ & =-\sum _{j^{\prime }\sim j}\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}\left(u_{i{j}^{\prime }}-\frac{1}{{\omega }_{ii}^{\left(s\right)}}\sum _{i^{\prime }\sim i}{\omega }_{i{i}^{\prime }}^{\left(s\right)}{u}_{i^{\prime }{j}^{\prime }}\right).\end{array}$(6) ) with (4(4) $\begin{array}{rl}& \mathrm{E}\left(u_{ij}-\frac{\rho }{d_i^{\left(s\right)}}\sum _{i^{\prime }\sim i}{u}_{i^{\prime }j}|u_{-\left\{ij\right\}}\right)\\ & =-\sum _{j^{\prime }\sim j}\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}\left(u_{i{j}^{\prime }}-\frac{\rho }{d_i^{\left(s\right)}}\sum _{i^{\prime }\sim i}{u}_{i^{\prime }{j}^{\prime }}\right).\end{array}$(4) ), instead of a homogeneous spatial smoothing with ρ, it has a heterogeneous specification with ${\omega }_{i{i}^{\prime }}^{\left(s\right)}$. This is hence more flexible in the spatial domain. The resulting joint precision matrix is (7) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\\ & \otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}},\end{array}$(7) where ${\mathbf{\Omega }}^{\left(r\right)}$ is a symmetric $J\times J$ matrix with entries ${\omega }_{j{j}^{\prime }}^{\left(r\right)}$ and ${\mathbf{\Omega }}^{\left(s\right)}$ is a symmetric $I\times I$ matrix with entries ${\omega }_{i{i}^{\prime }}^{\left(s\right)}$. A derivation of (7(7) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\\ & \otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}},\end{array}$(7) ) is given in Appendix 1. We again use ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ and ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$ for simplicity. In model (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ), the spatial part is the conventional conditional autoregressive model while in model (7(7) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\\ & \otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}},\end{array}$(7) ), it is modelled by a more flexible one ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$.

The precision matrix (7(7) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\\ & \otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}},\end{array}$(7) ) is positive definite if both ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ and ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$ are positive definite. In a Bayesian model, both can have G-Wishart priors. The specification has an obvious problem of identification: ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}=z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes \left(1/z\right){\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$, where z is an arbitrary constant scalar. Following Wang and West (2009), one can impose a constraint ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)},11}=1$ and add an auxiliary variable z. Then specify a joint prior on $(z,z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}})$: (8) $p(z,z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid b^{\left(r\right)},{\mathit{V}}^{\left(r\right)})\propto p_{GWis}(z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid b^{\left(r\right)},{\mathit{V}}^{\left(r\right)})\cdot 1,$(8) where $p_{GWis}(\cdot )$ is the density of G-Wishart distribution. Transform this joint density to $(z,{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}})$ and we obtain the desired joint prior. There is no additional constraint imposed on ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$ and let ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}\sim \text{GWis}(b^{\left(s\right)},{\mathit{V}}^{\left(s\right)})$. Inference and computation under this model are given in Appendix 2.

2.5. Priors for the graph

The two types of graphs used in this modelling framework should be treated differently. On one hand, the spatial graph should be treated known and fixed because the geographical locations and their neighbouring structure is fixed in most scenarios. On the other hand, the response graph should be treated unknown because we often know little about the relationship between multiple responses. In the literature of Gaussian graphical model determination, usually the unknown graph is assumed random and a prior on the graph is assigned. The Markov chain Monte Carlo (MCMC) sampling scheme, such as the reversible jump MCMC (Green, 1995), is often used to sample graphs from the posterior distribution. In this paper, we adopt and slightly modify existing MCMC algorithms for the graph determination (Dobra et al., 2011; Wang & Li, 2012), with computational details given in Appendix 2, for each aforementioned model. For the prior choice of ${\mathit{C}}^{\left(r\right)}$, consider (9) $\begin{array}{rl}P\left({\mathit{C}}^{\left(r\right)}\right)& \propto B(a+\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}({\mathit{C}}^{\left(r\right)}),b+m\\ & -\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)\left)/B\right(a,b),\end{array}$(9) where $B(\cdot ,\cdot )$ is the beta function, m is the total number of possible edges $\left(\genfrac{}{}{0ex}{}{J}{2}\right)$, $\text{size}\left({\mathit{C}}^{\left(r\right)}\right)\in \{0,1,\dots ,m\}$, and a and b are given hyperparameters. More details about this prior can be found in Scott and Berger (2006) and Scott and Carvalho (2009). The following prior is often used as well (Dobra et al., 2011): (10) $P\left({\mathit{C}}^{\left(r\right)}\right)\propto {\pi }^{\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)}(1-\pi {)}^{m-\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)},$(10) where $\pi \in (0,1)$ is a given hyperparameter. Sparser graphs can be favoured by choosing a small value for π. The prior (9(9) $\begin{array}{rl}P\left({\mathit{C}}^{\left(r\right)}\right)& \propto B(a+\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}({\mathit{C}}^{\left(r\right)}),b+m\\ & -\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)\left)/B\right(a,b),\end{array}$(9) ) can be obtained by integrating π out with a hyperprior $\text{Beta}(a,b)$ on π.

3. An application

We illustrate the proposed models with a real example of disease mapping. It is known that smoking is linked with multiple diseases in the population, of which leading diseases include lung diseases and heart diseases. The dataset under consideration here includes six variables, among which four variables are related to the smoke exposure and the other two are diseases. Obtained from the 2011 Missouri County Level Survey, the four smoke exposure variables are: Current Cigarette Smoking, Current Smokeless Tobacco Use, Current Other Tobacco Use, Exposure to Secondhand Smoke. Data are binary responses to the survey questionnaires (Yes or No), aggregated to each county level. The other two variables, obtained from the Surveillance, Epidemiology and End Results (SEER) program, are the Lung Cancer Mortality and the Heart Diseases Mortality, both of which are counts for each county within a specified time period. To summarise, we have I = 115 counties and J = 6 response variables. Let $n_{i1},\dots ,n_{i4}$ be the numbers of respondents in the survey and let $E_{i5}$ and $E_{i6}$ be the age-adjusted expected mortality for the two diseases. Then, the proportions $y_{ij}/n_{ij},j=1,\dots ,4$ are empirical estimates of the prevalences of the survey variables, and the proportions $y_{ij}/E_{ij},j=5,6$ are standardised mortality ratios of the diseases.

Consider a Bayesian hierarchical model for $y_{ij}$. We use the binomial-logit model and the Poisson-lognormal model (Banerjee, Gelfand, & Carlin, 2004) for $y_{i,1-4}$ and $y_{i,5-6}$, respectively, i.e., $\begin{array}{rl}& y_{ij}\sim \mathrm{B}\mathrm{i}\mathrm{n}(n_{ij},p_{ij}),\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left(p_{ij}\right)={\text{β}}_j+u_{ij},\\ i& =1,\dots ,115\mathrm{a}\mathrm{n}\mathrm{d}j=1,\dots ,4;\\ & y_{ij}\sim \text{Poi}\left(E_{ij}{\eta }_{ij}\right),\mathrm{l}\mathrm{o}\mathrm{g}\left({\eta }_{ij}\right)={\text{β}}_j+u_{ij},\\ i& =1,\dots ,115\mathrm{a}\mathrm{n}\mathrm{d}j=5,6.\end{array}$ For simplicity, we do not consider other covariates in this example. The primary interest here is to model the random effects $u_{ij}$, which are expected to be correlated in both the spatial domain and the response domain. To complete the model specification, specify a weakly-informative normal prior for the intercepts ${\text{β}}_j$ and a multivariate conditional autoregressive model $\mathrm{M}\mathrm{C}\mathrm{A}\mathrm{R}(\mathit{B},\mathit{T})$ for the random effects $\mathit{U}=\left\{u_{ij}\right\}$. We apply the three proposed versions of $\mathrm{M}\mathrm{C}\mathrm{A}\mathrm{R}(\mathit{B},\mathit{T})$ here. Hyperparameters for prior distributions are specified as follows. For the graph, noticing that the choice of π in (10(10) $P\left({\mathit{C}}^{\left(r\right)}\right)\propto {\pi }^{\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)}(1-\pi {)}^{m-\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)},$(10) ) can influence the posterior inference, we consider the prior (10(10) $P\left({\mathit{C}}^{\left(r\right)}\right)\propto {\pi }^{\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)}(1-\pi {)}^{m-\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\left({\mathit{C}}^{\left(r\right)}\right)},$(10) ) with both $\pi =0.2$ in favour of a sparse graph and $\pi =0.5$ as no preference. All other priors are chosen to be only weakly-informative and have little impact on the posterior inference. In Model 1, we specify hyperparameters in the inverse gamma prior as $a_j=b_j=0.5$. In Model 2, we specify hyperparameters in the G-Wishart prior as b = 3 and $\mathit{V}=\mathit{I}$. In Model 3, we specify hyperparameters in the two G-Wishart priors as $b=3,\mathit{V}=\mathit{I},b^{\left(s\right)}=24$ and ${\mathit{V}}^{\left(s\right)}=(b^{\left(s\right)}-2)(\mathit{D}-0.95\mathit{C}{)}^{-1}$, which implies a prior mode at a proper conditional autoregressive model. For each model, we perform the Markov chain Monte Carlo for 150,000 iterations with a burn-in size of 50,000. Posterior results are based on the remaining samples. Figure 1 shows the convergence of the log-joint-posterior and notice that they all converge quickly.

Figure 1. Data analysis: convergence of the log-joint posterior $\log p(\mathit{\theta }\mid \text{data})$ under the three models (first 50,000 iterations), where $\mathit{\theta }$ represents the collection of all parameters.

Table 1 shows the posterior edge inclusion probabilities for the response graph ${\mathit{C}}^{\left(r\right)}$. First, all three models seem to agree on the link between Cigarette Smoking and Secondhand Smoke Exposure, as well as the link between Lung Diseases Mortality and Heart Diseases Mortality. There is a moderate agreement on the links between Secondhand Smoke Exposure and Lung Diseases Mortality, and between Cigarette Smoking and Lung Diseases Mortality. In general, Model 1 tends to be a sparser graph, which is possibly due to the diagonal dominance condition. Model 2 is the simplest model as reflected by its pD, the effective number of parameters but has the largest DIC. Model 3 is the most flexible model among the three, and the inferred graph tends to be denser than the other two. It is as expected that its pD is larger but the overall criterion DIC is much smaller than the other two. Second, the edge inclusion probabilities are in general higher when $\pi =0.5$, as expected, but it has little material impact on the final inferred graph. The DIC has little change with different π values. Lastly, Figures 2–4  show the maps of spatial random effects $u_{ij}$ for the three models, respectively, and for a problem of disease mapping, this is often the eventual output for practitioners.

Figure 2. Data analysis: posterior means of spatial random effects $u_{ij}$ under Model 1.

Figure 3. Data analysis: posterior means of spatial random effects $u_{ij}$ under Model 2.

Figure 4. Data analysis: posterior means of spatial random effects $u_{ij}$ under Model 3.

Table 1. Data analysis: Posterior edge inclusion probabilities for the response graph and deviance information criterion for Models 1–3, respectively.

	Cigarette	Smokeless	Other	Secondhand	Lung	Heart	pD	DIC
	Model 1
Cigarette		0	0.165	0.963	0.276	0	460.8	5766.5
Smokeless	0		0.181	0	0	0.002
Other	0.258	0		0.032	0.014	0.005
Secondhand	0.983	0	0.062		0.453	0
Lung	0.301	0	0.040	0.544		0.533
Heart	0.005	0	0.039	0.007	0.772		460.2	5765.0
	Model 2
Cigarette		0.213	0.919	1	0.756	0.373	445.1	5813.8
Smokeless	0.340		0.207	0.281	0.249	0.706
Other	0.872	0.316		0.390	0.346	0.898
Secondhand	1	0.527	0.475		0.732	0.343
Lung	0.584	0.407	0.412	0.871		1
Heart	0.352	0.942	0.838	0.400	1		447.2	5812.6
	Model 3
Cigarette		0.268	0.328	0.827	0.541	0.436	543.3	5630.4
Smokeless	0.249		0.417	0.692	0.515	0.726
Other	0.317	0.403		0.534	0.453	0.510
Secondhand	0.819	0.694	0.531		0.936	0.801
Lung	0.521	0.500	0.435	0.941		0.997
Heart	0.411	0.732	0.508	0.805	0.998		543.1	5629.2

Notes: Inclusion probabilities higher than 0.5 are in bold. The upper-right panel is for $\pi =0.2$ and the lower-left panel is for $\pi =0.5$.

4. Simulation

To validate the proposed algorithms, we perform a simulation study on a $7\times 7$ regular grid (I = 49 area units) with J = 4 response variables. Consider the true response graph with two edges $C_{13}^{\left(r\right)}$ and $C_{24}^{\left(r\right)}$. In this simulation study, we do not consider the scenario with misspecified models, and therefore, data are generated under each of the three models and the correct model is then used for inference. The parameter settings are given as follows. For Model 1, ${\lambda }_j=0.95$, ${\phi }_{j{j}^{\prime }}={\psi }_{j{j}^{\prime }}=0.9$, ${\delta }_j=1$ and ${\text{β}}_j=1$. For Model 2, $\rho =0.9$, ${\omega }_{jj}=4$, ${\omega }_{13}={\omega }_{24}=-3.2$ and ${\text{β}}_j=1$. For Model 3, parameters are the same as Model 2 but ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$ is generated from $\text{GWis}(10,8(\mathit{D}-0.9\mathit{C}{)}^{-1})$. We repeat the simulation and inference process for L = 50 times and for each time, the MCMC iteration number is 5000. We consider three measures for validating and comparing the three algorithms. The first measure is the mean inclusion probability matrix with standard deviations. We call the second measure the error rate of mis-identified edges. If we use 0.5 as the threshold for identifying an edge in the graph, for each replication, we obtain an inferred graph and then compare with the true graph to record a proportion of wrong edges/non-edges. The error rate is the average proportion of L replications. The third measure is the mean absolute error (MAE) of random effects in the model, $\mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{L}\frac{1}{J}\frac{1}{I}\sum _l\sum _j\sum _i\left|\frac{{\hat{u}}_{ijl}-u_{ijl}}{u_{ijl}}\right|$ where $u_{ijl}$ is the true value and ${\hat{u}}_{ijl}$ is the posterior mean.

Simulation results are given in Table 2. For all three models, the algorithms can correctly identify the true edges. The algorithm for Model 1 appears to be unstable as the standard deviation is large and tends to underestimate inclusion probabilities, while the algorithm for Model 3 tends to overestimate inclusion probabilities for non-edges. The algorithm for Model 2 presents the smallest error rate and MAE. Note that this simulation study validates the proposed algorithms under correct model specifications and hence the result cannot imply that Model 2 is the best model for a real dataset. In fact, as shown in the data analysis, Model 2 is the simplest specification and is the least preferred model in that case according to DIC.

Table 2. Simulation: Mean posterior edge inclusion probabilities (standard deviations in parentheses), error rates of mis-identified edges, and the mean absolute errors of random effects.

	Var 2	Var 3	Var 4	Error Rate	MAE
	Model 1
Var 1	0.128 (0.214)	0.664 (0.358)	0.088 (0.107)	0.15	3.167
Var 2		0.100 (0.148)	0.615 (0.352)
Var 3			0.087 (0.155)
	Model 2
Var 1	0.211 (0.108)	1 (0)	0.200 (0.086)	0.02	2.153
Var 2		0.213 (0.134)	1 (0)
Var 3			0.199 (0.107)
	Model 3
Var 1	0.246 (0.037)	0.961 (0.051)	0.292 (0.054)	0.07	2.218
Var 2		0.438 (0.092)	0.996 (0.008)
Var 3			0.426 (0.083)

Note: True edges are in bold.

5. Further discussion

In this paper, we proposed a modelling framework for multivariate areal data from a graphical model perspective. We rebuilt three well-known models in our framework and developed Bayesian inference tools for the proposed models. It is our perspective that this framework is very general and can contain other models that are beyond the cases discussed in the paper. For example, Jin et al. (2007) specified a co-regionalised areal data model, in which their Case 3 is a very general specification. We show that this specification can be reproduced and extended in our framework. Consider the Cholesky decomposition $\mathbf{\Sigma }=\mathit{A}{\mathit{A}}^T$. Jin et al. (2007)'s Case 3 specification of the joint covariance matrix is $(\mathit{A}\otimes {\mathit{I}}_I)({\mathit{I}}_J\otimes {\mathit{D}}^{\left(s\right)}-\mathbf{\Phi }\otimes {\mathit{C}}^{\left(s\right)}{)}^{-1}(\mathit{A}\otimes {\mathit{I}}_I{)}^T$ whose inverse is then (11) $(\mathit{A}{\mathit{A}}^T{)}^{-1}\otimes {\mathit{D}}^{\left(s\right)}-({\mathit{A}}^{-1}{)}^T\mathbf{\Phi }{\mathit{A}}^{-1}\otimes {\mathit{C}}^{\left(s\right)}$(11) where $\mathbf{\Phi }$ is a symmetric $J\times J$ matrix. Let $\mathbf{\Omega }={\mathbf{\Sigma }}^{-1}=(\mathit{A}{\mathit{A}}^T{)}^{-1}$ and $\mathit{Q}=({\mathit{A}}^{-1}{)}^T\mathbf{\Phi }{\mathit{A}}^{-1}$. Obviously it is one-to-one from $(\mathit{A},\mathbf{\Phi })$ to $(\mathbf{\Omega },\mathit{Q})$. Specification (11(11) $(\mathit{A}{\mathit{A}}^T{)}^{-1}\otimes {\mathit{D}}^{\left(s\right)}-({\mathit{A}}^{-1}{)}^T\mathbf{\Phi }{\mathit{A}}^{-1}\otimes {\mathit{C}}^{\left(s\right)}$(11) ) is hence equivalent to (12) $\mathbf{\Omega }\otimes {\mathit{D}}^{\left(s\right)}-\mathit{Q}\otimes {\mathit{C}}^{\left(s\right)},$(12) where $\mathit{Q}$ is a symmetric $J\times J$ matrix with entries $q_{j{j}^{\prime }}$. To reproduce this specification in our framework, parameterise $\mathit{B}$ and $\mathit{T}$ as follows (assuming $i\ne i^{\prime }$ and $j\ne j^{\prime }$): $\begin{array}{rl}& b_{\left\{ij\right\},\left\{i^{\prime }j\right\}}=\frac{q_{jj}}{d_i^{\left(s\right)}{\omega }_{jj}}{C}_{\left\{ij\right\},\left\{i^{\prime }j\right\}},b_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}}=-\frac{{\omega }_{j{j}^{\prime }}}{{\omega }_{jj}}{C}_{\left\{ij\right\},\left\{i{j}^{\prime }\right\}},\\ & b_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}=\frac{q_{j{j}^{\prime }}}{d_i^{\left(s\right)}{\omega }_{jj}}{C}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}},{\tau }_{ij}^2=\frac{1}{d_i^{\left(s\right)}{\omega }_{jj}}.\end{array}$ This parameterisation leads to the joint precision matrix (13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) The expression (13(13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) ) reduces to (12(12) $\mathbf{\Omega }\otimes {\mathit{D}}^{\left(s\right)}-\mathit{Q}\otimes {\mathit{C}}^{\left(s\right)},$(12) ) which is equivalent to Jin et al. (2007)'s (11(11) $(\mathit{A}{\mathit{A}}^T{)}^{-1}\otimes {\mathit{D}}^{\left(s\right)}-({\mathit{A}}^{-1}{)}^T\mathbf{\Phi }{\mathit{A}}^{-1}\otimes {\mathit{C}}^{\left(s\right)}$(11) ) when ${\mathit{C}}^{\left(r\right)}$ is a complete graph. A derivation of (13(13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) ) is given in Appendix 1. The validity of this model relies on the positive definiteness of (13(13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) ). Jin et al. (2007) showed that it is positive definite if $\mathbf{\Omega }$ is positive definite and eigenvalues of $\mathbf{\Phi }={\mathit{A}}^T\mathit{Q}\mathit{A}$ are between $1/{\xi }_{\text{min}}$ and $1/{\xi }_{\text{max}}$, reciprocals of the smallest and largest eigenvalues of ${\mathit{D}}^{\left(s\right)-1/2}{\mathit{C}}^{\left(s\right)}{\mathit{D}}^{\left(s\right)-1/2}$, which are known constants. The graphical version (13(13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) ) must also satisfy this condition, that is, $1/{\xi }_{\text{min}}\le \lambda \left({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}^{-1}{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\right)\le 1/{\xi }_{\text{max}}$, where $\lambda \left(\mathit{M}\right)$ is any eigenvalue of $\mathit{M}$. Considering that both ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ and ${\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}$ are restricted by the underlying graph, the eigenvalue condition is not easy to implement in computations. This matter is worth investigating in the future.

In general, flexible models are desired for modelling multivariate areal data because overly simplistic models may misspecify the true underlying covariance structure. However, there is almost always a trade-off between the simplicity and the flexibility of a model. It is probably reasonable to allow certain flexibilities for specific purposes, such as in this paper, for learning a graphical relationship between multiple responses. It is usually the practitioner's choice whether a more flexible but complicated model is needed for the problem at hand, especially when the performance improvement is negligible.

Acknowledgments

The author would like to thank two editors and two anonymous referees for their helpful comments which lead to a much improved paper.

Disclosure statement

No potential conflict of interest was reported by the author.

ORCID

Ye Liang  http://orcid.org/0000-0002-6513-8962

Additional information

Notes on contributors

Ye Liang

Ye Liang is an associate professor in the Department of Statistics at Oklahoma State University. He received his Ph.D. in statistics from the University of Missouri.

Appendices

Appendix 1. Derivations

A.1. Derivation of Equation (3)

With the parameterisation in Model 1, we have $\mathit{T}={\mathit{D}}^{-1}(\mathbf{\Delta }\otimes {\mathit{I}}_I)$ and $\begin{array}{rl}\mathit{B}& =({\mathbf{\Delta }}^{1/2}\otimes {\mathit{I}}_I){\mathit{D}}^{-1}\left\{\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}+(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I\right.\\ & +\left.(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I).\end{array}$ Then immediately we have expression (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ) for ${\mathit{T}}^{-1}(\mathit{I}-\mathit{B})$.

A.2. Derivation of Equation (5)

With the parameterisation in Model 2, we have $\mathit{T}=\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\mathbf{\Omega })\otimes {\mathit{D}}^{\left(s\right)}{)}^{-1}$ and $\begin{array}{rl}\mathit{B}& =\rho {\mathit{I}}_J\otimes {\mathit{D}}^{\left(s\right)-1}{\mathit{C}}^{\left(s\right)}-\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{\Omega }{)}^{-1}\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)}\right)\otimes {\mathit{I}}_I\\ & +\rho \left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{\Omega }{)}^{-1}\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)}\right)\otimes {\mathit{D}}^{\left(s\right)-1}{\mathit{C}}^{\left(s\right)}.\end{array}$ Then $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathbf{\Omega }\right)\otimes {\mathit{D}}^{\left(s\right)}-\rho \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathbf{\Omega }\right)\otimes {\mathit{C}}^{\left(s\right)}\\ & +(\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{D}}^{\left(s\right)}-\rho (\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\\ & =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes \rho {\mathit{C}}^{\left(s\right)}\\ & =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$ which is expression (5(5) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)})\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}),\end{array}$(5) ).

A.3. Derivation of Equation (7)

With the parameterisation in Model 3, we have $\mathit{T}=\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right({\mathbf{\Omega }}^{\left(r\right)})\otimes \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\Omega }}^{\left(s\right)}){)}^{-1}$ and $\begin{array}{rl}\mathit{B}& =-{\mathit{I}}_J\otimes \left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\Omega }}^{\left(s\right)}{)}^{-1}{\mathbf{\Omega }}^{\left(s\right)}\circ {\mathit{C}}^{\left(s\right)}\right)\\ & -\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\Omega }}^{\left(r\right)}{)}^{-1}{\mathbf{\Omega }}^{\left(r\right)}\circ {\mathit{C}}^{\left(r\right)}\right)\otimes {\mathit{I}}_I\\ & -\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\Omega }}^{\left(r\right)}{)}^{-1}{\mathbf{\Omega }}^{\left(r\right)}\circ {\mathit{C}}^{\left(r\right)}\right)\\ & \otimes \left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\Omega }}^{\left(s\right)}{)}^{-1}{\mathbf{\Omega }}^{\left(s\right)}\circ {\mathit{C}}^{\left(s\right)}\right).\end{array}$ Then $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathbf{\Omega }}^{\left(r\right)}\right)\otimes \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathbf{\Omega }}^{\left(s\right)}\right)+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathbf{\Omega }}^{\left(r\right)}\right)\\ & \otimes ({\mathbf{\Omega }}^{\left(s\right)}\circ {\mathit{C}}^{\left(s\right)})\\ & +({\mathbf{\Omega }}^{\left(r\right)}\circ {\mathit{C}}^{\left(r\right)})\otimes \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathbf{\Omega }}^{\left(s\right)}\right)+({\mathbf{\Omega }}^{\left(r\right)}\circ {\mathit{C}}^{\left(r\right)})\\ & \otimes ({\mathbf{\Omega }}^{\left(s\right)}\circ {\mathit{C}}^{\left(s\right)})\\ & =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\},\end{array}$ which is expression (7(7) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{{\mathbf{\Omega }}^{\left(r\right)}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\\ & \otimes \left\{{\mathbf{\Omega }}^{\left(s\right)}\circ ({\mathit{I}}_I+{\mathit{C}}^{\left(s\right)})\right\}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}},\end{array}$(7) ).

A.4. Derivation of Equation (13)

With the parameterisation in Model 4, we have $\mathit{T}=\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\mathbf{\Omega })\otimes {\mathit{D}}^{\left(s\right)}{)}^{-1}$ and $\begin{array}{rl}\mathit{B}& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathbf{\Omega }{)}^{-1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\mathit{Q})\otimes {\mathit{D}}^{\left(s\right)-1}{\mathit{C}}^{\left(s\right)}\\ & -\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{\Omega }{)}^{-1}\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)}\otimes {\mathit{I}}_I\\ & +\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{\Omega }{)}^{-1}\mathit{Q}\circ {\mathit{C}}^{\left(r\right)}\otimes {\mathit{D}}^{\left(s\right)-1}{\mathit{C}}^{\left(s\right)}.\end{array}$ Then $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathbf{\Omega }\right)\otimes {\mathit{D}}^{\left(s\right)}-\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{Q}\right)\otimes {\mathit{C}}^{\left(s\right)}\\ & +\mathbf{\Omega }\circ {\mathit{C}}^{\left(r\right)}\otimes {\mathit{D}}^{\left(s\right)}-\mathit{Q}\circ {\mathit{C}}^{\left(r\right)}\otimes {\mathit{C}}^{\left(s\right)}\\ & =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)},\end{array}$ which is expression (13(13) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =\left\{\mathbf{\Omega }\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{D}}^{\left(s\right)}\\ & -\left\{\mathit{Q}\circ ({\mathit{I}}_J+{\mathit{C}}^{\left(r\right)})\right\}\otimes {\mathit{C}}^{\left(s\right)}\\ & ={\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{D}}^{\left(s\right)}-{\mathit{Q}}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}.\end{array}$(13) ).

Appendix 2: Bayesian computations

A.5. A hierarchical generalised linear model

For illustration, we now assume a full Bayesian hierarchical model and give computational details for this model. Assume binomial counts $y_{ij}/n_{ij}$ for J responses and I areal units. Specify a Bayesian model as follows, for $i=1,\dots ,I$ and $j=1,\dots ,J$, $\begin{array}{rl}& y_{ij}\mid p_{ij}\sim \mathrm{B}\mathrm{i}\mathrm{n}(n_{ij},p_{ij}),\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left(p_{ij}\right)={\text{β}}_j+u_{ij},\\ & {\text{β}}_j\sim \mathrm{N}(0,{\tau }_0^2),\mathit{U}\mid \mathit{B},\mathit{T}\sim \mathrm{M}\mathrm{C}\mathrm{A}\mathrm{R}(\mathit{B},\mathit{T}),\end{array}$ where ${\tau }_0^2$ is a given constant, $\mathit{U}$ is the matrix-variate of $u_{ij}$, and priors for $\mathit{B}$ and $\mathit{T}$ depend on the specific parameterisation. This section is organised as follows: we first give details of updating effects parameters ${\text{β}}_j$ and $u_{ij}$, and then, separately for each model, details of updating parameters of $\mathrm{M}\mathrm{C}\mathrm{A}\mathrm{R}$ and updating the random response graph ${\mathit{C}}^{\left(r\right)}$.

A.6. Updating effects parameters

Our experience has shown that the convergence is poor if we directly update ${\text{β}}_j$ and $u_{ij}$. We apply the hierarchical centring technique (Gelfand, Sahu, & Carlin, 1995) and block sampling. Let ${\gamma }_{ij}={\text{β}}_j+u_{ij}$ and $\mathit{\gamma }=\mathrm{v}\mathrm{e}\mathrm{c}\left[\right({\gamma }_{ij}{)}_{I\times J}]$ has a non-centred MCAR prior. We update $({\gamma }_{ij},{\text{β}}_j)$ instead of $(u_{ij},{\text{β}}_j)$. The full conditional distribution of ${\gamma }_{ij}$ is $\begin{array}{rl}& p({\gamma }_{ij}\mid \cdot )\propto \frac{e^{y_{ij}{\gamma }_{ij}}}{(1+e^{{\gamma }_{ij}}{)}^{n_{ij}}}\exp \left\{-\frac{1}{2{\tau }_{ij}^2}\left({\gamma }_{ij}-{\text{β}}_j\right.\right.\\ & -\left.\left.\sum _{\left\{i^{\prime }{j}^{\prime }\right\}\ne \left\{ij\right\}}{b}_{\left\{ij\right\},\left\{i^{\prime }{j}^{\prime }\right\}}({\gamma }_{i^{\prime }{j}^{\prime }}-{\text{β}}_{j^{\prime }})\right)\right\}.\end{array}$ We use Metropolis-Hastings algorithm to sample ${\gamma }_{ij}$ from this conditional density. We block sample $\text{β}$ in the following way. For now denote $\mathit{M}={\mathit{T}}^{-1}(\mathit{I}-\mathit{B})$, the joint precision matrix. Let ${\mathit{\gamma }}^{\ast }=\mathit{M}\mathit{\gamma }$ and ${\mathit{\gamma }}^{\ast \ast }$ be a $J\times 1$ vector such that ${\gamma }_1^{\ast \ast }$ is the sum of the first I elements in ${\mathit{\gamma }}^{\ast }$, ${\gamma }_2^{\ast \ast }$ is the sum of the second I elements in ${\mathit{\gamma }}^{\ast }$ and so on. Partition $\mathit{M}$ into $J\times J$ blocks and define $\mathit{H}=\left(\begin{array}{ccc}{\mathbf{1}}^T{\mathit{M}}_{11}\mathbf{1}& \cdots & {\mathbf{1}}^T{\mathit{M}}_{1J}\mathbf{1}\\ ⋮& \ddots & ⋮\\ {\mathbf{1}}^T{\mathit{M}}_{J1}\mathbf{1}& \cdots & {\mathbf{1}}^T{\mathit{M}}_{JJ}\mathbf{1},\end{array}\right)$ where $\mathbf{1}$ is the all-one vector. Then the full conditional distribution for the vector $\text{β}$ is $(\text{β}\mid \cdot )\sim \mathrm{N}\left[(\mathit{H}+1/{\tau }_0^2\mathit{I}{)}^{-1}{\mathit{\gamma }}^{\ast \ast },(\mathit{H}+1/{\tau }_0^2\mathit{I}{)}^{-1}\right].$

A.7. Model 1: updating ${\delta }_j$, ${\lambda }_j$, ${\psi }_{j{j}^{\prime }}$, ${\phi }_{j{j}^{\prime }}$ and ${\mathit{C}}^{\left(r\right)}$

Given the current graph ${\mathit{C}}^{\left(r\right)}$, parameters are updated through Gibbs sampling. Recall priors on these parameters: ${\delta }_j\sim \text{IG}(a_j,b_j)$ and ${\lambda }_j,{\psi }_{j{j}^{\prime }},{\phi }_{j{j}^{\prime }}\sim \text{Unif}(-1,1)$. Let ${\mathit{u}}_j$ be the jth column vector of $\mathit{U}$, $j=1,\dots ,J$ and ${\mathit{D}}_j$ be the jth diagonal block of $\mathit{D}$. The full conditional distribution of ${\delta }_j$ is given by $\begin{array}{rl}& p({\delta }_j\mid \cdot )\propto {\delta }_j^{-I/2-a_j-1}\exp \left\{-\frac{1}{2{\delta }_j}{\mathit{u}}_j^T({\mathit{D}}_j-{\lambda }_j{\mathit{C}}^{\left(s\right)}){\mathit{u}}_j\right.\\ & +\left.\sum _{j^{\prime }\sim j}\frac{1}{\sqrt{{\delta }_j{\delta }_{j^{\prime }}}}{\mathit{u}}_j^T({\psi }_{j{j}^{\prime }}{\mathit{I}}_I+{\phi }_{j{j}^{\prime }}{\mathit{C}}^{\left(s\right)}){\mathit{u}}_{j^{\prime }}-\frac{b_j}{{\delta }_j}\right\}.\end{array}$ It can be shown that the transformed one $(\sqrt{1/{\delta }_j}\mid \cdot )$ is log-concave when $I+2a_j-1>0$. Thus, we use the adaptive rejection sampling to update ${\delta }_j$.

Let $\mathit{W}$ be an $I\times J$ matrix, where $\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{W}\right)=({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}\right)$ and ${\mathit{w}}_j$ be the jth column vector of $\mathit{W}$. Let $\mathit{M}=\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}$ as in (3(3) $\begin{array}{rl}{\mathit{T}}^{-1}(\mathit{I}-\mathit{B})& =({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I)\left\{\mathit{D}-\mathbf{\Lambda }\otimes {\mathit{C}}^{\left(s\right)}\right.\\ & \left.-(\mathbf{\Psi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{I}}_I-(\mathbf{\Phi }\circ {\mathit{C}}^{\left(r\right)})\otimes {\mathit{C}}^{\left(s\right)}\right\}\\ & \times ({\mathbf{\Delta }}^{-1/2}\otimes {\mathit{I}}_I),\end{array}$(3) ). Then ${\lambda }_j$, ${\psi }_{j{j}^{\prime }}$ and ${\phi }_{j{j}^{\prime }}$ are sequentially updated through following full conditional distributions, $\begin{array}{rl}p({\lambda }_j\mid \cdot )& \propto {\left|\mathit{M}\left({\lambda }_j\right)\right|}^{1/2}\exp \left(\frac{1}{2}{\lambda }_j{\mathit{w}}_j^T{\mathit{C}}^{\left(s\right)}{\mathit{w}}_j\right),\\ p({\psi }_{j{j}^{\prime }}\mid \cdot )& \propto {\left|\mathit{M}\left({\psi }_{j{j}^{\prime }}\right)\right|}^{1/2}\exp \left({\psi }_{j{j}^{\prime }}{\mathit{w}}_j^T{\mathit{w}}_{j^{\prime }}\right),\\ p({\phi }_{j{j}^{\prime }}\mid \cdot )& \propto {\left|\mathit{M}\left({\phi }_{j{j}^{\prime }}\right)\right|}^{1/2}\exp \left({\phi }_{j{j}^{\prime }}{\mathit{w}}_j^T{\mathit{C}}^{\left(s\right)}{\mathit{w}}_{j^{\prime }}\right).\end{array}$ We use Metropolis-Hastings algorithm to update these parameters. Note that evaluating the sparse $|\mathit{M}|$ could be computationally intensive. An efficient algorithm, usually based on the Cholesky decomposition, on sparse matrices is helpful.

The graph ${\mathit{C}}^{\left(r\right)}$ is updated through a simple reversible jump MCMC algorithm. Propose a new graph ${\mathit{C}}^{\left(r\right)\ast }$ by only adding or deleting one edge from ${\mathit{C}}^{\left(r\right)}$. Without loss of generality, suppose that one edge $\{j0,k0\}$ is added to the new graph. Dimension has been changed by 2 from $({\mathit{C}}^{\left(r\right)},\mathbf{\Psi },\mathbf{\Phi })$ to $({\mathit{C}}^{\left(r\right)\ast },{\mathbf{\Psi }}^{\ast },{\mathbf{\Phi }}^{\ast })$. Propose $u_1\sim U(-1,1)$ and $u_2\sim U(-1,1)$, and let ${\psi }_{j0,k0}^{\ast }=u_1$ and ${\phi }_{j0,k0}^{\ast }=u_2$. The Jacobian from $(\mathbf{\Psi },\mathbf{\Phi },u_1,u_2)$ to $({\mathbf{\Psi }}^{\ast },{\mathbf{\Phi }}^{\ast })$ hence is 1. Choose a Bernoulli jump proposal with odds $q({\mathit{C}}^{\left(r\right)\ast },{\mathit{C}}^{\left(r\right)})/q({\mathit{C}}^{\left(r\right)},{\mathit{C}}^{\left(r\right)\ast })=p\left({\mathit{C}}^{\left(r\right)}\right)/p\left({\mathit{C}}^{\left(r\right)\ast }\right)$ and systematically scan through the graph for updating. Accept the move from ${\mathit{C}}^{\left(r\right)}$ to ${\mathit{C}}^{\left(r\right)\ast }$ with probability $min\{1,\alpha \}$ where $\alpha =\frac{|{\mathit{M}}^{\ast }{|}^{1/2}}{|\mathit{M}{|}^{1/2}}\exp \left\{-\frac{1}{2}\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{W}\right)({\mathit{M}}^{\ast }-\mathit{M})\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{W}\right)\right\}.$

A.8. Model 2: updating ρ, ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ and ${\mathit{C}}^{\left(r\right)}$

Given the current graph ${\mathit{C}}^{\left(r\right)}$, parameters are updated through Gibbs sampling. Recall priors on these parameters: $\rho \sim U(-1,1)$ and ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\sim \text{GWis}(b,\mathit{V})$. Use Metropolis-Hastings algorithm to update ρ. It can be shown that the full conditional distribution for ρ is $\begin{array}{rl}& p(\rho \mid \cdot )\propto {\left|{\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}\right|}^{J/2}\exp \\ & \times \left\{\frac{\rho }{2}\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}{)}^T\right({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\otimes {\mathit{C}}^{\left(s\right)}\left)\mathrm{v}\mathrm{e}\mathrm{c}\right(\mathit{U})\right\}.\end{array}$ Let ${\mathit{W}}_1$ be an $I\times J$ matrix, where $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{W}}_1\right)=[{\mathit{I}}_J\otimes ({\mathit{D}}^{\left(s\right)}-\rho {\mathit{C}}^{\left(s\right)}{)}^{1/2}\left]\mathrm{v}\mathrm{e}\mathrm{c}\right(\mathit{U})$ and ${\mathit{w}}_{1,j}$ be the jth column vector of ${\mathit{W}}_1$. Let $\mathit{S}$ be an $J\times J$ matrix with $s_{j{j}^{\prime }}={\mathit{w}}_{1,j}^T{\mathit{w}}_{1,j^{\prime }}$. Then the full conditional distribution of ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$ is $\begin{array}{rl}& p({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid \cdot )\propto {\left|{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\right|}^{\left(b+I-2\right)/2}\exp \\ & \times \left[-\frac{1}{2}\mathrm{t}\mathrm{r}\left\{{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}(\mathit{V}+\mathit{S})\right\}\right]1_{{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\in M^{+}\left({\mathit{C}}^{\left(r\right)}\right)},\end{array}$ which is GWis($b+I,\mathit{V}+\mathit{S}$). For sampling from the G-Wishart distribution, we use the block Gibbs sampler, given the set of maximum cliques, introduced by Wang and Li (2012).

The graph ${\mathit{C}}^{\left(r\right)}$ is updated using (Wang & Li, 2012)'s partial analytic structure algorithm (p. 188, Algorithm 2).

A.9. Model 3: updating ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}$, ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}$ and ${\mathit{C}}^{\left(r\right)}$

Given the current graph ${\mathit{C}}^{\left(r\right)}$, parameters are updated through Gibbs sampling. Recall that we impose a constraint and use a joint prior (8(8) $p(z,z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid b^{\left(r\right)},{\mathit{V}}^{\left(r\right)})\propto p_{GWis}(z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid b^{\left(r\right)},{\mathit{V}}^{\left(r\right)})\cdot 1,$(8) ) on $(z,z{\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}})$ and have ${\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}\sim \text{GWis}(b^{\left(s\right)},{\mathit{V}}^{\left(s\right)})$. Let both ${\mathit{W}}^{\left(r\right)}$ and ${\mathit{W}}^{\left(s\right)}$ be $I\times J$ matrices, where $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{W}}^{\left(r\right)}\right)=({\mathit{I}}_J\otimes {\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}^{1/2})\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}\right)$ and $\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{W}}^{\left(s\right)}\right)=({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}^{1/2}\otimes {\mathit{I}}_I)\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{U}\right)$. Let ${\mathit{w}}_j^{\left(r\right)}$ be the jth column vector of ${\mathit{W}}^{\left(r\right)}$ and ${\mathit{w}}_i^{\left(s\right)}$ be the ith row vector of ${\mathit{W}}^{\left(s\right)}$. Then let ${\mathit{S}}^{\left(r\right)}$ be $J\times J$ with $s_{j{j}^{\prime }}^{\left(r\right)}={\mathit{w}}_j^{\left(r\right)T}{\mathit{w}}_j^{\left(r\right)}$ and ${\mathit{S}}^{\left(s\right)}$ be $I\times I$ with $s_{i{i}^{\prime }}^{\left(s\right)}={\mathit{w}}_i^{\left(s\right)T}{\mathit{w}}_i^{\left(s\right)}$. With these notations, we have $(z\mid \cdot )\sim Ga(a_z,b_z),$ where $a_z=J(b-2)/2+\nu \left({\mathit{C}}^{\left(r\right)}\right)$ and $b_z=\mathrm{t}\mathrm{r}\left({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}{\mathit{V}}^{\left(r\right)}\right)/2$; $({\mathbf{\Omega }}_{{\mathit{C}}^{\left(r\right)}}\mid \cdot )\sim \mathrm{G}\mathrm{W}\mathrm{i}\mathrm{s}(b^{\left(r\right)}+I,z{\mathit{V}}^{\left(r\right)}+{\mathit{S}}^{\left(r\right)})$ and $({\mathbf{\Omega }}_{{\mathit{C}}^{\left(s\right)}}\mid \cdot )\sim \mathrm{G}\mathrm{W}\mathrm{i}\mathrm{s}(b^{\left(s\right)}+J,{\mathit{V}}^{\left(s\right)}+{\mathit{S}}^{\left(s\right)}).$ The graph ${\mathit{C}}^{\left(r\right)}$ is updated using (Wang & Li, 2012)'s partial analytic structure algorithm (p. 188, Algorithm 2).