The conditional autoregressive model is one of the routines for applications with univariate areal data. However, many authors have shown that it is not trivial to extend the univariate specification to multivariate situations. The difficulties lie in many aspects, including validity, interpretability, flexibility and computational feasibility. We approach this problem from an element-based perspective which utilizes given graphical information and builds joint adjacency structures. Our framework is general but three special specifications are discussed in details with Bayesian computations, all of which are related to existing models in literature. We provide an example with public health data, not only to illustrate the implementation, but also to compare these specifications. \\ Keywords: \normalfont Bayesian analysis; Disease mapping; Graphical model; G-Wishart distribution; Markov random field; Multivariate areal data.
Macros.tex
% FONT DEFINITIONS
\def\reduce#1{{\small  #1}}
\def\c#1{\ensuremath{\mathcal{#1}}}
\def\foot#1{\footnote{#1}}

% Bold face difinitions 
%\newcommand{\nmathbf}{}
\newcommand{\nmathbf}{\bm}

\def\bfA{\nmathbf A}
\def\bfB{\nmathbf B}
\def\bfC{\nmathbf C}
\def\bfD{\nmathbf D}
\def\bfE{\nmathbf E}
\def\bfF{\nmathbf F}
\def\bfG{\nmathbf G}
\def\bfH{\nmathbf H}
\def\bfI{\nmathbf I}
\def\bfJ{\nmathbf J}
\def\bfK{\nmathbf K}
\def\bfL{\nmathbf L}
\def\bfM{\nmathbf M}
\def\bfN{\nmathbf N}
\def\bfO{\nmathbf O}
\def\bfP{\nmathbf P}
\def\bfQ{\nmathbf Q}
\def\bfR{\nmathbf R}
\def\bfS{\nmathbf S}
\def\bfT{\nmathbf T}
\def\bfU{\nmathbf U}
\def\bfV{\nmathbf V}
\def\bfW{\nmathbf W}
\def\bfX{\nmathbf X}
\def\bfY{\nmathbf Y}
\def\bfZ{\nmathbf Z}

\def\bfa{\nmathbf a}
\def\bfb{\nmathbf b}
\def\bfc{\nmathbf c}
\def\bfd{\nmathbf d}
\def\bfe{\nmathbf e}
\def\bff{\nmathbf f}
\def\bfg{\nmathbf g}
\def\bfh{\nmathbf h}
\def\bfi{\nmathbf i}
\def\bfj{\nmathbf j}
\def\bfk{\nmathbf k}
\def\bfl{\nmathbf l}
\def\bfm{\nmathbf m}
\def\bfn{\nmathbf n}
\def\bfo{\nmathbf o}
\def\bfp{\nmathbf p}
\def\bfq{\nmathbf q}
\def\bfr{\nmathbf r}
\def\bfs{\nmathbf s}
\def\bft{\nmathbf t}
\def\bfu{\nmathbf u}
\def\bfv{\nmathbf v}
\def\bfw{\nmathbf w}
\def\bfx{\nmathbf x}
\def\bfy{\nmathbf y}
\def\bfz{\nmathbf z}

\def\bfalpha  {\nmathbf \alpha}
\def\bfbeta   {\nmathbf \beta}
\def\bfgamma  {\nmathbf \gamma}
\def\bfdelta  {\nmathbf \delta}
\def\bfepsilon{\nmathbf \epsilon}
\def\bfzeta   {\nmathbf \zeta}
\def\bfeta    {\nmathbf \eta}
\def\bftheta  {\nmathbf \theta}
\def\bfiota   {\nmathbf \iota}
\def\bfkappa  {\nmathbf \kappa}
\def\bflambda {\nmathbf \lambda}
\def\bfmu     {\nmathbf \mu}
\def\bfnu     {\nmathbf \nu}
\def\bfxi     {\nmathbf \xi}
\def\bfomicron{\nmathbf \omicron}
\def\bfpi     {\nmathbf \pi}
\def\bfrho    {\nmathbf \rho}
\def\bfsigma  {\nmathbf \sigma}
\def\bftau    {\nmathbf \tau}
\def\bfupsilon{\nmathbf \upsilon}
\def\bfphi    {\nmathbf \phi}
\def\bfpsi    {\nmathbf \psi}
\def\bfchi    {\nmathbf \chi}
\def\bfomega  {\nmathbf \omega}

\def\bfAlpha  {\nmathbf \Alpha}
\def\bfBeta   {\nmathbf \Beta}
\def\bfGamma  {\nmathbf \Gamma}
\def\bfDelta  {\nmathbf \Delta}
\def\bfEpsilon{\nmathbf \Epsilon}
\def\bfZeta   {\nmathbf \Zeta}
\def\bfEta    {\nmathbf \Eta}
\def\bfTheta  {\nmathbf \Theta}
\def\bfIota   {\nmathbf \Iota}
\def\bfKappa  {\nmathbf \Kappa}
\def\bfLambda {\nmathbf \Lambda}
\def\bfMu     {\nmathbf \Mu}
\def\bfNu     {\nmathbf \Nu}
\def\bfXi     {\nmathbf \Xi}
\def\bfOmicron{\nmathbf \Omicron}
\def\bfPi     {\nmathbf \Pi}
\def\bfRho    {\nmathbf \Rho}
\def\bfSigma  {\nmathbf \Sigma}
\def\bfTau    {\nmathbf \Tau}
\def\bfUpsilon{\nmathbf \Upsilon}
\def\bfPhi    {\nmathbf \Phi}
\def\bfPsi    {\nmathbf \Psi}
\def\bfChi    {\nmathbf \Chi}
\def\bfOmega  {\nmathbf \Omega}

\newcommand{\ttheta}{\tilde{\theta}}
\newcommand{\bfzero}{{\nmathbf 0}}
\newcommand{\bfone}{{\nmathbf 1}}
\newcommand{\vareps}{\varepsilon}
\def\bfvareps{\nmathbf \varepsilon}
\newcommand{\tgamma}{\tilde\gamma}

%% Calligraphic 
\newcommand{\cfA}{\mbox{\c{A}}}
\newcommand{\cfB}{\mbox{\c{B}}}
\newcommand{\cfC}{\mbox{\c{C}}}
\newcommand{\cfD}{\mbox{\c{D}}}
\newcommand{\cfE}{\mbox{\c{E}}}
\newcommand{\cfF}{\mbox{\c{F}}}
\newcommand{\cfG}{\mbox{\c{G}}}
\newcommand{\cfH}{\mbox{\c{H}}}
\newcommand{\cfI}{\mbox{\c{I}}}
\newcommand{\cfJ}{\mbox{\c{J}}}
\newcommand{\cfK}{\mbox{\c{K}}}
\newcommand{\cfL}{\mbox{\c{L}}}
\newcommand{\cfM}{\mbox{\c{M}}}
\newcommand{\cfN}{\mbox{\c{N}}}
\newcommand{\cfO}{\mbox{\c{O}}}
\newcommand{\cfP}{\mbox{\c{P}}}
\newcommand{\cfQ}{\mbox{\c{Q}}}
\newcommand{\cfR}{\mbox{\c{R}}}
\newcommand{\cfS}{\mbox{\c{S}}}
\newcommand{\cfT}{\mbox{\c{T}}}
\newcommand{\cfU}{\mbox{\c{U}}}
\newcommand{\cfV}{\mbox{\c{V}}}
\newcommand{\cfX}{\mbox{\c{X}}}
\newcommand{\cfY}{\mbox{\c{Y}}}
\newcommand{\cfZ}{\mbox{\c{Z}}}

\def\boldfacefake#1{\kern-4pt
   \hbox{ \mathsurround=0pt
   \hbox to 0.4pt{$#1$\hss}\hbox to 0.4pt{$#1$\hss}\hbox {$#1$}}}

\def\bfitI{\mbox{\boldfacefake{\it I}}}

\newcommand{\bcfA}{\boldsymbol{\mathcal{A}}}
\newcommand{\bcfB}{\boldsymbol{\mathcal{B}}}
\newcommand{\bcfC}{\boldsymbol{\mathcal{C}}}
\newcommand{\bcfD}{\boldsymbol{\mathcal{D}}}
\newcommand{\bcfE}{\boldsymbol{\mathcal{E}}}
\newcommand{\bcfF}{\boldsymbol{\mathcal{F}}}
\newcommand{\bcfG}{\boldsymbol{\mathcal{G}}}
\newcommand{\bcfH}{\boldsymbol{\mathcal{H}}}
\newcommand{\bcfI}{\boldsymbol{\mathcal{I}}}
\newcommand{\bcfJ}{\boldsymbol{\mathcal{J}}}
\newcommand{\bcfK}{\boldsymbol{\mathcal{K}}}
\newcommand{\bcfL}{\boldsymbol{\mathcal{L}}}
\newcommand{\bcfM}{\boldsymbol{\mathcal{M}}}
\newcommand{\bcfN}{\boldsymbol{\mathcal{N}}}
\newcommand{\bcfO}{\boldsymbol{\mathcal{O}}}
\newcommand{\bcfP}{\boldsymbol{\mathcal{P}}}
\newcommand{\bcfQ}{\boldsymbol{\mathcal{Q}}}
\newcommand{\bcfR}{\boldsymbol{\mathcal{R}}}
\newcommand{\bcfS}{\boldsymbol{\mathcal{S}}}
\newcommand{\bcfT}{\boldsymbol{\mathcal{T}}}
\newcommand{\bcfU}{\boldsymbol{\mathcal{U}}}
\newcommand{\bcfV}{\boldsymbol{\mathcal{V}}}
\newcommand{\bcfW}{\boldsymbol{\mathcal{W}}}
\newcommand{\bcfX}{\boldsymbol{\mathcal{X}}}
\newcommand{\bcfY}{\boldsymbol{\mathcal{Y}}}
\newcommand{\bcfZ}{\boldsymbol{\mathcal{Z}}}

% MATHEMATICAL NOTATION 

% Operators
\newcommand{\g}{\,\vert\,}
\newcommand{\p}{\mbox{P}}
\newcommand{\D}{\mbox{D}}
\newcommand{\E}{\mbox{E}}
\newcommand{\Mo}{\mbox{Mo}}
\newcommand{\Me}{\mbox{Me}}
\newcommand{\Cov}{\mbox{Cov}}
\newcommand{\Var}{\mbox{Var}}
\newcommand{\Corr}{\mbox{Corr}}
\newcommand{\Q}{\mbox{Q}}

% Distributions
\newcommand{\Bb}{\mbox{Bb}}
\newcommand{\Be}{\mbox{Be}}
\newcommand{\Bi}{\mbox{Bi}}
\newcommand{\Br}{\mbox{Br}}
\newcommand{\Ca}{\mbox{Ca}}
\newcommand{\Di}{\mbox{Di}}
\newcommand{\Ex}{\mbox{Ex}}
\newcommand{\Fs}{\mbox{Fs}}
\newcommand{\Ga}{\mbox{Ga}}
\newcommand{\Ge}{\mbox{Ge}}
\newcommand{\Gg}{\mbox{Gg}}
\newcommand{\Hy}{\mbox{Hy}}
\newcommand{\Ig}{\mbox{Ig}}
\newcommand{\Ip}{\mbox{Ip}}
\newcommand{\Lo}{\mbox{Lo}}
\newcommand{\Mu}{\mbox{Mu}}
\newcommand{\N}{\mbox{N}}
\newcommand{\Nb}{\mbox{Nb}}
\newcommand{\Ng}{\mbox{Ng}}
\newcommand{\Nw}{\mbox{Nw}}
\newcommand{\Pa}{\mbox{Pa}}
\newcommand{\Po}{\mbox{Po}}
\newcommand{\Pg}{\mbox{Pg}}
\newcommand{\Pn}{\mbox{Pn}}
\newcommand{\Ra}{\mbox{Ra}}
\newcommand{\St}{\mbox{St}}
\newcommand{\Un}{\mbox{Un}}
\newcommand{\Wi}{\mbox{Wi}}

% General Mathematics
\newcommand{\dd}[1]{\,d#1}
\newcommand{\barx}{\mbox{$\overline x$}}
\newcommand{\comb}[2]{{#1\choose#2}}
\newcommand{\ontop}[2]{{#1\atop#2}}
\newcommand{\h}{\hbox{$1\over2$}}
\newcommand{\ok}{\hfill\fbox{}}

% Abbreviations
\newcommand{\brow}[2]{\mbox{$\{{#1}_1,\ldots,{#1}_{#2}\}$}}
\newcommand{\prow}[2]{\mbox{$({#1}_1,\ldots,{#1}_{#2})$}}
\newcommand{\row}[2]{\mbox{${#1}_1,\ldots,{#1}_{#2}$}}
\newcommand{\data}{\row{x}{n}}
\newcommand{\bdata}{\prow{x}{n}}
\newcommand{\ie}{\emph{i.e.},\ }
\newcommand{\co}{\emph{cf.}\ }
\newcommand{\eg}{\emph{e.g.}, }
\newcommand{\etalc}{\emph{et al.},\ }
\newcommand{\etal}{\emph{et al.}\ }

% ARRAYS 
% Mathematical arrays
\newenvironment{mat}{\left(\begin{array}}{\end{array}\right)}

% Cases definition
\newcommand{\twomat}[5]{\begin{array}{ll}\displaystyle
              #1&\mbox{#2}\\[#5pt]#3&\mbox{#4}\end{array}}

\newcommand{\twocases}[6]
             {#1=\left\{\twomat{#2}{#3}{#4}{#5}{#6}\right.}

\newcommand{\mymatrix}[4]
           {\left(\begin{array}{ll}{#1} & {#2}\\
                                   {#3} & {#4}
                    \end{array} \right)}

% Tables
\newcommand{\btable}{\begin{table}[h]\centering}
\newcommand{\etable}{\end{table}}
\newcommand{\bt}{\begin{parag}\small \let\b=\nsb \let\sb=\nssb \begin{tabular}}
\newcommand{\et}{\end{tabular}\let\b=\nb \let\sb=\nsb\end{parag}}
\newcommand{\capt}[1]
         {\begin{quotation}\caption{{\small #1 }}\vs{-2}\end{quotation}}


% POSITIONING
\newenvironment{parag}{\par}{\par}
\newenvironment{dif}
    {\begin{parag}\small \let\b=\nsb \let\sb=\nssb \begin{parag}}
    {\let\b=\nb \let\sb=\nsb \end{parag}\end{parag}}

\newenvironment{proof}{\begin{dif} \noindent{\em Proof.~}}
            {\ok\vspace*{10pt}\end{dif}}
\newenvironment{exa}{\begin{list}{}
           {\setlength{\leftmargin}{10pt}
            \setlength{\rightmargin}{\leftmargin}}
           \item\begin{ex}\em}{\end{ex}\end{list}}

% Operative size command
\newcommand{\be}{\begin{eqnarray}}
\newcommand{\ee}{\end{eqnarray}}
\newcommand{\ba}{\begin{eqnarray*}}
\newcommand{\ea}{\end{eqnarray*}}

\newcommand{\go}{\rightarrow}
\newcommand{\goi}{\rightarrow \infty}
\newcommand{\ul}{\underline}
\newcommand{\ol}{\overline}
\newcommand{\fr}{\frac}
\newcommand{\pn}{\par\noindent}
\newcommand{\nc}{\nonumber\\}
\newcommand{\ssum}{\mbox{$\sum$}}
\newcommand{\hhline}{\hline\hline}

\newtheorem{theorem0}{Theorem}
\newtheorem{lemma0}{Lemma}
\newtheorem{remark0}{Remark}
\newtheorem{fact0}{Fact}
\newtheorem{example0}{Example}
\newtheorem{definition0}{Definition}
\newtheorem{corollary0}{Corollary}
\newtheorem{proposition0}{Proposition}
\newtheorem{algorithmY}{Algorithm}

\newenvironment{theorem}{\begin{theorem0} \mbox{} }{\end{theorem0}}
\newenvironment{lemma}{\begin{lemma0} \mbox{}}{\end{lemma0}}
\newenvironment{remark}{\begin{remark0} \mbox{}}{\end{remark0}}
\newenvironment{fact}{\begin{fact0} \mbox{}}{\end{fact0}}
\newenvironment{example}{\begin{example0} }{\end{example0}}
\newenvironment{definition}{\begin{definition0} \mbox{}}{\end{definition0}}
\newenvironment{corollary}{\begin{corollary0} \mbox{} }{\end{corollary0}}
\newenvironment{proposition}{\begin{proposition0}\mbox{} }{\end{proposition0}}
\newenvironment{algorithm1}{\begin{algorithmY}\mbox{} }{\end{algorithmY}}

\newcommand{\reals}{\mbox{\rm I\kern-.20em R}}
\newcommand{\sreals}{\mbox{\small \rm I\kern-.20em R}}
\newcommand{\mylinel}{\renewcommand{\baselinestretch}{1.8}\tiny\small}
\newcommand{\goto}{\rightarrow}
\newcommand{\expect}{\E}

\newcommand{\bdfn}{\begin{dfn}}
\newcommand{\edfn}{\end{dfn}}
\newcommand{\bteo}{\begin{teo}}
\newcommand{\eteo}{\end{teo}}
\newcommand{\bexa}{\begin{exa}}
\newcommand{\eexa}{\end{exa}}
\newcommand{\bdif}{\begin{dif}}
\newcommand{\edif}{\end{dif}}
\newcommand{\bpro}{\begin{proof}}
\newcommand{\epro}{\end{proof}}



main.tex
\documentclass[12pt]{article}
\usepackage{natbib,bm}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage[dvips]{epsfig}
\usepackage{lscape}

\voffset=-1.0in
\hoffset=-1.0in
\oddsidemargin=1.0truein
\evensidemargin=1.0truein
\topmargin=1.0truein

\textwidth=6.5truein
\textheight=8.0truein
\renewcommand{\baselinestretch}{1.1}
\tiny\normalsize
\input Macros
\setlength{\parskip}{1.0mm}

\makeatletter
\def\vec{\mathop{\operator@font vec}\nolimits}
\makeatother

\DeclareMathOperator{\diag}{diag}

\begin{document}
\input text 

\renewcommand{\baselinestretch}{1.1}

\input{main.bbl}

%\tiny\normalsize
%\bibliographystyle{agsm}
%\bibliography{ref}

%\newpage
%\input table
%\newpage
%\input graph

\end{document}

text.tex

\begin{center}
 {\large \bf A Graph-Based Multivariate Conditional Autoregressive Model}

\bigskip
 Ye Liang \\
 {\small \it Department of Statistics, Oklahoma State University, Stillwater, Oklahoma 74078, U.S.A.} \\
 {\small ye.liang@okstate.edu}

\end{center}
\medskip

\begin{abstract}
The conditional autoregressive model is one of the routines for applications with univariate
areal data. However, many authors have shown that it is not trivial to extend the univariate specification
to multivariate situations. The difficulties lie in many aspects, including validity, interpretability, flexibility 
and computational feasibility. We approach this problem from an element-based perspective
which utilizes given graphical information and builds joint adjacency structures. Our framework is general 
but three special specifications are discussed in details with Bayesian computations, all of which are 
related to existing models in literature. We provide an example with public health data, not only 
to illustrate the implementation, but also to compare these specifications. \\

\bf Keywords: \normalfont Bayesian analysis; Disease mapping; Graphical model; G-Wishart distribution; 
	Markov random field; Multivariate areal data.	  
\end{abstract}

%\newpage

\section{Introduction} \label{sec:intro}
%
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 neighboring 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 a Markov property. Hence, one often considers areal data as a Markov random field.
The univariate conditional autoregressive (\textsc{CAR}) model originated from \cite{Besag:1974} is 
a Gaussian Markov random field model, for which the joint distribution is multivariate Gaussian. Let $\bfu=(u_1,\ldots,u_I)^T$
be a vector of random variables on $I$ areal units (i.e. $I$ vertices). The zero-centered conditional autoregressive model
specifies full conditional Gaussian distributions
\ba
	u_i\mid u_{-i} \sim \mbox{N}\left( \sum_{i'\neq i} b_{ii'}u_{i'}, \tau_{i}^2 \right), ~~~ i=1,\ldots,I,
\ea
where $u_{-i}$ is the collection of $u_{i', i'\neq i}$. Its resulting joint distribution, by the Brook's lemma, 
has a density form
\ba
	f(\bfu\mid \bfT_{\textsc{CAR}},\bfB_{\textsc{CAR}})\propto \mbox{exp} \left\{ -\frac{1}{2}\bfu^T \bfT^{-1}_{\textsc{CAR}}
	(\bfI-\bfB_{\textsc{CAR}}) \bfu \right\},
\ea
where $\bfI$ is an identity matrix; $\bfB_{\textsc{CAR}}$ is the $I \times I$ matrix with off-diagonal entries $b_{ii'}$
and diagonal entries zero and $\bfT_{\textsc{CAR}}=\diag\{\tau_{1}^2,\ldots,\tau_{I}^2\}$.
The joint distribution is multivariate Gaussian if and only if $\bfT^{-1}_{\textsc{CAR}} (\bfI-\bfB_{\textsc{CAR}})$ is symmetric and positive
definite. A proper parameterization on $\bfB_{\textsc{CAR}}$ and $\bfT_{\textsc{CAR}}$ is needed. 
Consider a so-called adjacency matrix $\bfC_{\textsc{CAR}}$ for the undirected graph, where the $ii'$th entry $C_{ii'}=1$ 
if unit $i$ and unit $i'$ are neighbors (denoted as $i\sim i'$) and $C_{ii'}=0$ otherwise.
One popular parameterization is to let $b_{ii'}=\rho C_{ii'}/C_{i+}$ and $\tau_i^2=\sigma^2/C_{i+}$, where
$C_{i+}$ is the $i$th row sum of $\bfC_{\textsc{CAR}}$, representing the total number of neighbors of unit $i$. 
Let $\bfD_{\textsc{CAR}}=\diag\{C_{1+},\ldots,C_{I+}\}$. When $\rho$ is strictly between the smallest and largest eigenvalues 
of $\bfD_{\textsc{CAR}}^{-1/2}\bfC_{\textsc{CAR}}\bfD_{\textsc{CAR}}^{-1/2}$ (or sufficiently $|\rho|<1$) and $\sigma^2>0$, 
the joint distribution of $\bfu$ is a zero-mean multivariate Gaussian distribution with a covariance matrix 
$\sigma^2(\bfD_{\textsc{CAR}}-\rho\bfC_{\textsc{CAR}})^{-1}$. This is called the proper conditional autoregressive model.
When $\rho=1$, it is called the intrinsic conditional autoregressive model which is improper. 

Turning to the multivariate situation, consider $J$ responses (e.g. multiple diseases) on $I$ areal units. 
Let $\bfU$ be an $I \times J$ matrix-variate where the $ij$th entry $u_{ij}$ is a random variable for the  
$i$th areal unit and $j$th response. Each column of $\bfU$ is an areal vector for a single response and hence 
can be modeled by the univariate conditional autoregressive model. However, a multivariate model is desired for 
the matrix-variate $\bfU$ in order to simultaneously model the dependence across responses and areal units. 
Initially proposed by \cite{Mardia:1988}, the multivariate conditional autoregressive model specifies full 
conditional distributions on row vectors
of $\bfU$. Let $\bfu_i$ be the $i$th row vector of $\bfU$. In the spirit of \cite{Besag:1974}, specify
\be \label{mardia}
	\bfu_i\mid \bfu_{-i} \sim \mbox{N}\left(\sum_{i' \neq i} \bfB_{ii'}\bfu_{i'}, \bfSigma_i \right), ~~~ i=1,\ldots,I,
\ee
 where $\bfB_{ii'}$ and $\bfSigma_{i}$ are $J \times J$ matrices needing a further parameterization. 
 To make the joint distribution for $\vec(\bfU^T)$ a multivariate Gaussian, $\bfB_{ii'}$ and $\bfSigma_{i}$ must 
 satisfy certain conditions \citep{Mardia:1988}. \cite{Gelfand:2003} showed a convenient parameterization, 
 $\bfB_{ii'}=(\rho C_{ii'}/C_{i+})\bfI_J$ and $\bfSigma_i=\bfSigma/C_{i+}$. When $|\rho|<1$ and $\bfSigma$ is
 positive definite, $\vec(\bfU^T)$ has a zero-mean multivariate Gaussian distribution with a covariance matrix 
 $(\bfD_{\textsc{CAR}}-\rho\bfC_{\textsc{CAR}})^{-1}\otimes\bfSigma$. It is clear that this multivariate specification is a Kronecker product
 formula where $(\bfD_{\textsc{CAR}}-\rho\bfC_{\textsc{CAR}})^{-1}$ models the covariance structure across rows of $\bfU$
 (spatial domain) and $\bfSigma$ models the covariance structure across columns of $\bfU$ (response domain).
 
 From the modeling perspective, Mardia's specification has a difficulty with parameterization. 
 It is usually difficult to have a meaningful parameterization for $\bfB_{ii'}$ and $\bfSigma_i$ unless one 
 pursues a simple formulation. On the other hand, when we have graphical knowledge about the $J$ responses,
 it is also difficult to incorporate this information directly. It is arguable that Mardia's specification presents 
 a disconnect: the between vector variation is specified through a precision structure (the conditional autoregressive model) 
 but the within vector variation is specified through a conditional covariance structure ($\bfSigma_i$). To remedy this disconnect,
 one shall consider either the whole covariance structure or the whole precision structure. In this paper,
 we investigate the latter case as it is a natural extension of the univariate conditional autoregressive model and is particularly useful 
 when graphical information is available. 
 
 We shall point out other recent work on multivariate conditional autoregressive models. \cite{Kim:2001} developed a twofold model
 for the bivariate situation and \cite{Jin:2005} introduced a conditional specification also for the 
 bivariate situation. Multivariate situations were considered by \cite{Gelfand:2003}, \cite{Jin:2007}, 
 \cite{MacNab:2011}, \cite{Mart:2013} and many others. We will show the connection between our framework
 and some of these earlier work. The paper is organized as follows. Section \ref{sec:spec} is our general framework and three 
 special specifications. We present a real example with public health data in Section \ref{sec:app}. 
 Section \ref{sec:dis} is a closing section with discussions.
 
\section{A Graph-Based Model} \label{sec:spec}
%
\subsection{General framework}
Instead of specifying full conditional distributions on vectors like (\ref{mardia}), we approach this problem from 
an element-based perspective which seems more natural. This consideration also appeared in \citep{Sain:2011}. 
In the spirit of \cite{Besag:1974}, specify full conditional distributions for each element $u_{ij}$ in the matrix-variate $\bfU$,
\ba
	u_{ij}\mid u_{-\{ij\}} \sim \mbox{N}\left( \sum_{\{i'j'\}\neq \{ij\}} b_{\{ij\},\{i'j'\}}u_{i'j'}, \tau_{ij}^2 \right),
	~~~ i=1,\ldots,I ~\mbox{and}~j=1,\ldots,J,
\ea
where $\{i'j'\}\neq \{ij\}$ means either $i'\neq i$ or $j'\neq j$. In fact, here all elements in $\bfU$ form a lattice
system. By the Brook's lemma, the resulting joint distribution for $\vec(\bfU)$ has a density form
\ba
	f(\vec(\bfU)\mid \bfB,\bfT)\propto \exp\left\{ -\frac{1}{2}\vec(\bfU)^T \bfT^{-1}(\bfI-\bfB)\vec(\bfU) \right\},
\ea
where $\bfI$ is an $IJ \times IJ$ identity matrix, 
$\bfT=\diag\{\tau_{11}^2,\ldots,\tau_{I1}^2,\ldots,\tau_{1J}^2,\ldots,\tau_{IJ}^2 \}$ and $\bfB$
can be expressed block-wisely,
\ba
	\bfB=\begin{pmatrix}
		\bfB_{11} & \cdots & \bfB_{1J} \\
		\vdots & \ddots & \vdots \\
		\bfB_{J1} & \cdots & \bfB_{JJ} 
	\end{pmatrix}
	\mbox{~~~where~~~}
	\bfB_{jj'}=
	\begin{pmatrix}
		b_{\{1j\},\{1j'\}} & \cdots & b_{\{1j\},\{Ij'\}} \\
		\vdots & \ddots & \vdots \\
		b_{\{Ij\},\{1j'\}} & \cdots & b_{\{Ij\},\{Ij'\}} 
	\end{pmatrix},	
\ea
and the diagonal elements $b_{\{ij\},\{ij\}}$ are zero.
The joint distribution for $\vec(\bfU)$ is multivariate Gaussian if and only if $\bfT^{-1}(\bfI-\bfB)$ is symmetric
and positive definite. It is desired that $\bfB$ and $\bfT$ are parameterized. We denote this general model
$\textsc{MCAR}(\bfB,\bfT)$ for later use.

Consider the adjacency structure of the undirected graph for $\bfU$. In the univariate situation,
the adjacency structure is naturally determined by the spatial map. Two areal units 
are connected by an edge if they are neighbors geographically. However, it is not quite obvious that which elements 
should be neighbors in $\bfU$, which makes the topic somewhat subjective. Consider a situation in which
there is an undirected graph for the $J$ responses. Let $\bfC^{(s)}$ be the adjacency matrix for all areal units
and $\bfC^{(r)}$ be the adjacency matrix for all responses. Both spatial graph and response graph are then 
uniquely determined by $\bfC^{(s)}$ and $\bfC^{(r)}$, respectively. Let $\bfC$ be the joint adjacency matrix for
the lattice system $\bfU$. A general construction of $\bfC$ can be made through $\bfC^{(s)}$ and $\bfC^{(r)}$,
\be\label{adjacency}
	\bfC=\bfC^{(r)} \otimes \bfC^{(s)} + \bfC^{(r)} \otimes \bfI_I + \bfI_J \otimes \bfC^{(s)}.
\ee
This construction connects $u_{ij}$ with $u_{i'\sim i,j}$, $u_{i,j'\sim j}$ and $u_{i'\sim i, j'\sim j}$,
meaning its spatial neighbor, response neighbor and interaction neighbor, respectively. 
One may add edges for secondary neighbors or drop edges when needed. For example, some other 
special constructions would be: (i) $\bfC=\bfI_J \otimes \bfC^{(s)}$ (independent conditional autoregressive models, 
no dependence between responses); (ii) $\bfC=\bfC^{(r)} \otimes \bfI_I$ (independent multivariate variables, 
no spatial dependence); (iii) $\bfC=\bfC^{(r)} \otimes \bfI_I + \bfI_J \otimes \bfC^{(s)}$ (drop edges with the 
interaction neighbor $u_{i'\sim i, j'\sim j}$). The construction (\ref{adjacency}) here is analogous to the
concept of separability in spatio-temporal modeling, implying that the structure of the overall lattice system 
can be separated. We will focus on (\ref{adjacency}) throughout this paper though other constructions are possible. 

Let $C_{\{ij\},\{i'j'\}}$ denote entries of $\bfC$, analogous to the block-wise notation $b_{\{ij\},\{i'j'\}}$ for $\bfB$.
Let  $d_j^{(r)}$ be the $j$th row sum in $\bfC^{(r)}$ and $d_i^{(s)}$ be the $i$th row sum in $\bfC^{(s)}$.
Then the $ij$th row sum in $\bfC$ is $d_{ij}=d_j^{(r)}d_i^{(s)}+d_j^{(r)}+d_i^{(s)}$.
Let $\bfD^{(r)}=\diag\{d_1^{(r)},\ldots,d_J^{(r)}\}$, $\bfD^{(s)}=\diag\{d_1^{(s)},\ldots,d_I^{(s)}\}$ and 
$\bfD=\diag\{d_{11},\ldots,d_{I1},\ldots,d_{1J},\ldots,d_{IJ}\}$. Parameterization of $\bfB$ and $\bfT$ 
can be made on the basis of the joint adjacency structure. We mainly discuss three specifications
in following subsections, all of which are related to existing models in literature.  

\subsection{Model 1: multifold specification with different linkage parameters}
\cite{Kim:2001} developed a twofold conditional autoregressive model for bivariate areal data ($J=2$), 
using different linkage parameters for different types of neighbor. 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 parameterize $\bfB$ and $\bfT$ in the following way (assuming $i \neq i'$ and $j \neq j'$):
\ba
	&& b_{\{ij\},\{i'j\}} = \frac{\lambda_j}{d_{ij}} C_{\{ij\},\{i'j\}}, ~~~ 
		b_{\{ij\},\{ij'\}} = \frac{\psi_{jj'}}{d_{ij}} \sqrt{\frac{\delta_j}{\delta_{j'}}} C_{\{ij\},\{ij'\}},  \\
	&& b_{\{ij\},\{i'j'\}} = \frac{\phi_{jj'}}{d_{ij}} \sqrt{\frac{\delta_j}{\delta_{j'}}} C_{\{ij\},\{i'j'\}}, ~~~
		\tau_{ij}^2=\frac{\delta_j}{d_{ij}},
\ea
where $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$ are linkage parameters and $\delta_j$ are variance components.    
Linkage parameters are for three types of neighbor: $u_{i'\sim i,j}$, $u_{i,j'\sim j}$ and $u_{i'\sim i, j'\sim j}$.
Having this specification, the conditional mean of $u_{ij}$ essentially is
\ba
	\mbox{E}(u_{ij} \mid u_{-\{ij\}})=\frac{1}{d_{ij}}\left( \lambda_j\sum_{i' \sim i}u_{i'j}
	+\psi_{jj'}\sum_{j'\sim j}\sqrt{\frac{\delta_j}{\delta_{j'}}}u_{ij'}
	+\phi_{jj'}\sum_{i'\sim i}\sum_{j'\sim j}\sqrt{\frac{\delta_j}{\delta_{j'}}} u_{i'j'} \right),
\ea
which is a weighted average of all its neighbors in $\bfC$. This specification generalizes 
\cite{Kim:2001}' twofold model and hence could be called a multifold specification. 
It can be shown that, for this parameterization, the joint precision matrix is
\be \label{Model1}
	\bfT^{-1}(\bfI-\bfB) &=& (\bfDelta^{-\frac{1}{2}}\otimes\bfI_I ) \left\{\bfD-\bfLambda\otimes\bfC^{(s)} \right. \notag \\
	&& \left. -(\bfPsi\circ\bfC^{(r)})\otimes\bfI_I - (\bfPhi\circ\bfC^{(r)})\otimes\bfC^{(s)}\right\}
	(\bfDelta^{-\frac{1}{2}}\otimes\bfI_I ),
\ee
where $\bfDelta=\diag\{\delta_1,\ldots,\delta_J\}$, $\bfLambda=\diag\{\lambda_1,\ldots,\lambda_J\}$,
$\bfPsi$ and $\bfPhi$ are $J \times J$ symmetric matrices with entries $\psi_{jj'}$ and  $\phi_{jj'}$, respectively,
and the operator $\circ$ means an element-wise product. Derivation of (\ref{Model1}) is given in Appendix 1.
Note that only nonzero entries of $\bfPsi$ and $\bfPhi$ are parameters in the model, the number of which depends 
on $\bfC^{(r)}$. \cite{Jin:2007} mentioned a "dependent non-identical latent spatial process" which is also a special case 
of Model 1 when $\bfC=\bfC^{(r)} \otimes \bfC^{(s)} + \bfI_J \otimes \bfC^{(s)}$, implying that $\bfPsi$ is zero.

In order to make (\ref{Model1}) positive definite, constraints on $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$ are needed,
with $\delta_j>0$ assumed. In general, it is difficult to find a sufficient and necessary condition for the positive
definiteness of (\ref{Model1}). Even such a condition on $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$ can be found,
computational difficulties are expected in the inference procedure. \cite{Kim:2001}'s solution to this problem 
was a sufficient condition: $\mbox{max}\{|\lambda_j|, |\psi_{jj'}|, |\phi_{jj'}|; \forall j, j'\}<1$. 
The matrix (\ref{Model1}) under this condition is diagonally dominant and hence is positive definite. 
Though their proof was under $J=2$, it is generally true for any $J$ (by the same argument).
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 condition. We expect this to be a major deficiency of the multifold model 
and further evidences can be found in Section \ref{sec:app}. We adopt this solution for now due to lack
of better solutions.

In a Bayesian model, priors on parameters $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$ can be chosen based on
their actual constraints. In our case, a uniform prior $\mbox{Unif}(-1,1)$ is adequate for these linkage parameters. 
Priors on the variance components $\delta_j$ can be weakly-informative inverse-gamma priors $\mbox{IG}(a_j,b_j)$. 
Computational details are discussed in Appendix 2.

\subsection{Model 2: completely separable specification with homogeneous spatial smoothing}
\cite{Gelfand:2003}'s Kronecker-product model is a convenient parameterization of \cite{Mardia:1988}'s
model. In our framework, this specification can be achieved and extended by having the following 
parameterization for $\bfB$ and $\bfT$ (assuming $i \neq i'$ and $j \neq j'$):
\ba
	&& b_{\{ij\},\{i'j\}} = \frac{\rho}{d_{i}^{(s)}} C_{\{ij\},\{i'j\}}, ~~~ 
		b_{\{ij\},\{ij'\}} = -\frac{\omega_{jj'}}{\omega_{jj}} C_{\{ij\},\{ij'\}}, \\
	&& b_{\{ij\},\{i'j'\}} = \frac{\rho\omega_{jj'}}{d_{i}^{(s)}\omega_{jj}} C_{\{ij\},\{i'j'\}}, ~~~
		\tau_{ij}^2=\frac{1}{d_{i}^{(s)}\omega_{jj}},
\ea
where $\rho$ and $\omega_{jj'}$ are linkage parameters, and $1/\omega_{jj}$ are variance components.
This parameterization is not very straightforward at first glance, but is much clearer in the form
of conditional mean:
\be \label{cond_mean1}
	\mbox{E}\left( u_{ij}-\frac{\rho}{d_i^{(s)}}\sum_{i'\sim i}u_{i'j} ~\Bigg|~ u_{-\{ij\}} \right)
	=-\sum_{j'\sim j}\frac{\omega_{jj'}}{\omega_{jj}}\left( u_{ij'}-\frac{\rho}{d_i^{(s)}}\sum_{i'\sim i}u_{i'j'} \right).
\ee
Note that for a single response, the univariate conditional autoregressive model specifies 
$\mbox{E}(u_i | u_{-i})=\rho\sum_{i'\sim i} u_{i'}/d_i^{(s)}$. In the multivariate setting, 
$\rho\sum_{i'\sim i} u_{i'j}/d_i^{(s)}$ is no longer the conditional mean for $u_{ij}\mid u_{-\{ij\}}$
and their conditional difference is regressed on other differences through $\omega_{jj'}$. 
This parameterization yields the joint precision matrix
\be \label{Model2}
	\bfT^{-1}(\bfI-\bfB)= \left\{\bfOmega\circ(\bfI_J+\bfC^{(r)})\right\} \otimes (\bfD^{(s)}-\rho\bfC^{(s)}),
\ee 
where $\bfOmega$ is a symmetric $J\times J$ matrix with entries $\omega_{jj'}$. Derivation of (\ref{Model2})
is given in Appendix 1. The linkage parameter
$\rho$ is interpreted as a spatial smoothing parameter and $\bfOmega$ controls dependence across 
$J$ responses. It is noteworthy that only nonzero entries in $\bfOmega$ are parameters in the model and we denote
$\bfOmega_{\bfC^{(r)}}=\bfOmega\circ(\bfI_J+\bfC^{(r)})$ for simplicity. The notation $\bfOmega_{\bfC^{(r)}}$,
commonly used in graphical models, means the precision matrix for a graph $\bfC^{(r)}$. 
The model (\ref{Model2}) is a natural extension of \cite{Gelfand:2003}'s model.  When $\bfC^{(r)}$ is the complete graph 
(any two vertices are connected), $\bfOmega_{\bfC^{(r)}}$ is free of zero entries. Then let 
$\bfSigma=\bfOmega^{-1}$ and (\ref{Model2}) is equivalent to \cite{Gelfand: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 $\rho$ is the same for any $i$ and $i'$ which distinguishes Model 2 from Model 3.
The limitation of this model is apparently its inflexibility due to the common linkage $\rho$. 
 
The joint precision matrix (\ref{Model2}) is positive definite if $|\rho|<1$ and $\bfOmega_{\bfC^{(r)}}$ is positive
definite. In a Bayesian model, an appropriate prior on the graphical precision matrix $\bfOmega_{\bfC^{(r)}}$ is desired. A reasonable choice is the G-Wishart distribution \citep{Atay:2005,Letac:2007}. 
The G-Wishart distribution is a conjugate family for the precision matrix on a Gaussian graphical model, 
whose density function for $\bfOmega_{\bfC^{(r)}}$ is given by
\ba
	p(\bfOmega_{\bfC^{(r)}} \mid b,\bfV)=I_{\bfC^{(r)}}(b,\bfV)^{-1}
	\left| \bfOmega_{\bfC^{(r)}} \right|^{\frac{b-2}{2}}\mbox{exp}
	\left\{-\frac{1}{2}\mbox{tr}(\bfV\bfOmega_{\bfC^{(r)}})\right\}
	1_{\bfOmega_{\bfC^{(r)}}\in M^+({\bfC^{(r)}})},
\ea
where $b>2$ is the number of degrees of freedom; $\bfV$ is the scale matrix, analogous to the regular Wishart
distribution; $I_{\bfC^{(r)}}(\cdot)$ is the normalizing constant and $M^+({\bfC^{(r)}})$ is the cone
of possible symmetric positive definite precision matrices on the graph $\bfC^{(r)}$. The requirement of 
positive definiteness is naturally satisfied by using this G-Wishart distribution. It is also practically attractive
because of its conjugacy. That said, for a prior distribution $\mbox{GWis}(b,\bfV)$ and a given 
sample covariance matrix $\bfS$ of sample size $n$,
the posterior distribution of $\bfOmega_{\bfC^{(r)}}$ is $\mbox{GWis}(b+n,\bfV+\bfS)$.
Computational details are given in Appendix 2.

\subsection{Model 3: completely separable specification with heterogenous spatial smoothing}
\cite{Dobra:2011} introduced a multivariate lattice model by giving Kronecker product G-Wishart priors
to the matrix-variate $\bfU$. In our framework, $\bfB$ and $\bfT$ can be parameterized in the following way
(assuming $i \neq i'$ and $j \neq j'$):
\ba
	&& b_{\{ij\},\{i'j\}} = -\frac{\omega_{ii'}^{(s)}}{\omega_{ii}^{(s)}} C_{\{ij\},\{i'j\}}, ~~~
		b_{\{ij\},\{ij'\}} = -\frac{\omega_{jj'}^{(r)}}{\omega_{jj}^{(r)}} C_{\{ij\},\{ij'\}}, \\
	&& b_{\{ij\},\{i'j'\}} = -\frac{\omega_{ii'}^{(s)}\omega_{jj'}^{(r)}}{\omega_{ii}^{(s)}\omega_{jj}^{(r)}} C_{\{ij\},\{i'j'\}}, ~~~ 
		\tau_{ij}^2=\frac{1}{\omega_{ii}^{(s)}\omega_{jj}^{(r)}},
\ea
which is equivalent to the version of conditional mean
\be \label{cond_mean2}
	\mbox{E}\left( u_{ij}-\frac{1}{\omega_{ii}^{(s)}}\sum_{i'\sim i}\omega_{ii'}^{(s)}u_{i'j} ~\Bigg|~ u_{-\{ij\}} \right)
	=-\sum_{j'\sim j}\frac{\omega_{jj'}}{\omega_{jj}}\left( u_{ij'}-\frac{1}{\omega_{ii}^{(s)}}\sum_{i'\sim i}
	\omega_{ii'}^{(s)}u_{i'j'} \right).
\ee
Comparing (\ref{cond_mean2}) with (\ref{cond_mean1}), instead of a homogeneous spatial smoothing 
with $\rho$, we make it heterogeneous with $\omega_{ii'}^{(s)}$ for different $i$ and $i'$. 
This treatment makes the model more flexible in the spatial domain.     
The resulting joint precision matrix for this specification is
\be \label{Model3}
	\bfT^{-1}(\bfI-\bfB)= \left\{\bfOmega^{(r)}\circ(\bfI_J+\bfC^{(r)})\right\} \otimes 
		\left\{\bfOmega^{(s)}\circ(\bfI_I+\bfC^{(s)})\right\},
\ee
where $\bfOmega^{(r)}$ is a symmetric $J \times J$ matrix with entries $\omega^{(r)}_{jj'}$ and 
$\bfOmega^{(s)}$ is a symmetric $I \times I$ matrix with entries $\omega^{(s)}_{ii'}$. Derivation of (\ref{Model3}) 
is given in Appendix 1. We let $\bfOmega_{\bfC^{(r)}}=\bfOmega^{(r)}\circ(\bfI_J+\bfC^{(r)})$ and 
$\bfOmega_{\bfC^{(s)}}=\bfOmega^{(s)}\circ(\bfI_I+\bfC^{(s)})$ for simplicity. Again, we can see the difference 
between (\ref{Model3}) and (\ref{Model2}). In model (\ref{Model2}), the spatial part is the conventional 
conditional autoregressive model while in model (\ref{Model3}), it is modeled by $\bfOmega_{\bfC^{(s)}}$. 
\cite{Dobra:2011} showed that the G-Wishart prior for $\bfOmega_{\bfC^{(s)}}$ can contain the 
conventional version as a mode when $b$ and $\bfV$ are properly specified. 

The precision matrix (\ref{Model3}) is positive definite if both $\bfOmega_{\bfC^{(r)}}$ and $\bfOmega_{\bfC^{(s)}}$
are positive definite. In a Bayesian model, both can have G-Wishart priors. The specification has an obvious 
problem of identification: $\bfOmega_{\bfC^{(r)}} \otimes  \bfOmega_{\bfC^{(s)}}=z\bfOmega_{\bfC^{(r)}} \otimes  
(1/z)\bfOmega_{\bfC^{(s)}}$, where $z$ is an arbitrary constant scalar. Following \cite{Wang:2009}, we impose 
a constraint $\bfOmega_{\bfC^{(r)},11}=1$ and add an auxiliary variable $z$. Then specify a joint prior on
$(z,z\bfOmega_{\bfC^{(r)}})$:
\be \label{jointprior}
	p(z,z\bfOmega_{\bfC^{(r)}} \mid b^{(r)}, \bfV^{(r)})\propto p_{GWis}
	(z\bfOmega_{\bfC^{(r)}} \mid b^{(r)}, \bfV^{(r)})\cdot 1,
\ee
where $p_{GWis}(\cdot)$ is the density of G-Wishart distribution. Transform this joint density to $(z,\bfOmega_{\bfC^{(r)}})$
and we obtain the desired joint prior. Impose no constraint on $\bfOmega_{\bfC^{(s)}}$ and specify 
$\bfOmega_{\bfC^{(s)}}\sim \mbox{GWis}(b^{(s)},\bfV^{(s)})$. Computational details are given 
in Appendix 2.

\subsection{A brief summary of three models}
All three models are based on the joint graph (\ref{adjacency}) and our general framework. They are all 
related to existing multivariate models. We do not discuss other specifications in this paper,
but it is possible that there are better specifications under this framework. The joint precision matrix for 
Model 1, 2 and 3 are given by (\ref{Model1}), (\ref{Model2}) and (\ref{Model3}), respectively.
Constraints are imposed for all three models to ensure the positive definiteness of the joint precision matrix. 

We now briefly compare all three models. Let function $\nu(\cdot)$ be the total number of vertices and
edges in a graph. Model 1 specifies $\vec(\bfU)\mid \delta_j, \lambda_j,\psi_{jj'},\phi_{jj'},~j\neq j'$. 
The total number of parameters in the model is $2\nu(\bfC^{(r)})$, which depends on the response graph only.
The constraints on parameters are $\mbox{max}\{|\lambda_j|, |\psi_{jj'}|, |\phi_{jj'}|; \forall j, j'\}<1$ and $\delta_j>0$,
for all $j$.
Model 2 specifies $\vec(\bfU)\mid \rho,\bfOmega_{\bfC^{(r)}}$ whose total number of parameters is
$\nu(\bfC^{(r)})+1$, which depends on the response graph only. The constraints on parameters are 
$|\rho|<1$ and that $\bfOmega_{\bfC^{(r)}}$ is positive definite. Model 3 specifies 
$\vec(\bfU)\mid \bfOmega_{\bfC^{(r)}},\bfOmega_{\bfC^{(s)}}$ whose total number of parameters is 
$\nu(\bfC^{(r)})+\nu(\bfC^{(s)})$, which depends on both the response graph and spatial graph.
The constraints on parameters are that both $\bfOmega_{\bfC^{(r)}}$ and $\bfOmega_{\bfC^{(s)}}$
are positive definite. It is clear that Model 2 has the least number of parameters for any $J>1$
and $I>1$ (trivially satisfied). Model 3 often has more parameters than Model 1 because the size
of a spatial graph is often larger than a response graph. Not only can the complexity of the model affect
the computational feasibility, but also affect most model selection scores (penalized by a model selection criteria). 
However, on the other hand, a more complicated model usually has greater flexibility. We will not give an
ultimate recommendation on the three models because the question here is about a tradeoff. Model comparisons
with real data are discussed in Section \ref{sec:app}. 

\section{An Example in Multiple Disease Mapping} \label{sec:app}
%
We present an example with public health data in Missouri. In particular, we analyze a county level dataset with five response
variables, which implies that $I=115$ and $J=5$. Two of the response variables are survey counts in each county: $y_{i1}$ 
(smoking) and $y_{i2}$ (no health care), which are collected from the 2011 Missouri County Level Survey
Study. The other three response variables are mortality counts in each county: $y_{i3}$ (lung cancer), $y_{i4}$ (heart diseases) 
and $y_{i5}$ (COPD), which are collected from the Surveillance, Epidemiology and End Results (SEER) within a given time period. 
Let $n_{i1}$ and $n_{i2}$ be the numbers of respondents for the survey, and let $E_{i3}$, $E_{i4}$ and $E_{i5}$ be the age-adjusted 
expected numbers for the three diseases. Proportions $y_{ij}/n_{ij},~j=1,2$ give empirical estimates of the prevalence 
and proportions $y_{ij}/E_{ij},~j=3,4,5$ give standardized mortality ratios. 

We are now interested in jointly modeling all $y_{ij}$ in a Bayesian hierarchical framework. We use the binomial-logit model
and the Poisson-lognormal model \citep{Banerjee:2004} for $y_{i,1-2}$ and $y_{i,3-5}$, respectively, i.e.
\ba
	&& y_{ij} \sim \mbox{Bin}(n_{ij}, p_{ij}), ~~~ \mbox{logit}(p_{ij})=\beta_{j}+u_{ij}, ~~~ i=1,\ldots,I \mbox{ and } j=1,2; \\
	&& y_{ij} \sim \mbox{Poi}(E_{ij}\eta_{ij}), ~~~ \mbox{log}(\eta_{ij})=\beta_{j}+u_{ij}, ~~~ i=1,\ldots,I \mbox{ and } j=3,4,5.
\ea
Parameters $p_{ij}$ are underlying binomial prevalences and parameters $\eta_{ij}$ are underlying
relative risks. In our example, no covariates are included and hence $u_{ij}$ is interpreted as a spatial smoother. 
We specify a vague normal prior for $\beta_j$ and three different $\textsc{MCAR}(\bfB,\bfT$) for $\bfU=\{u_{ij}\}$. 
All computations for this model then follow Appendix 2.

Our discussion in this paper focuses on given response graphs though graph determination is possible with adjusted
computation steps \citep{Wang:2012, Dobra:2011}. To empirically determine the response graph for this example, 
we compute $r_{ij}=\mbox{logit}(y_{ij}/n_{ij}), ~j=1,2$ and $r_{ij}=\mbox{log}(y_{ij}/E_{ij}). ~j=3,4,5$. We treat $r_{ij}$
as 115 independent multivariate normal samples for the five response variables (ignoring the spatial dependence). 
For this small five-variable problem, we use algorithms in \cite{Wang:2012} to determine the graph and choose
the one with the highest posterior probability. This empirically assumed response graph is shown in Figure \ref{fig:graph},
along with the actual spatial graph and the corresponding joint precision structure. Note that, under adjacency construction
(\ref{adjacency}), this precision structure will be the same for all three models discussed before. We shall later
compare this assumed graph with the complete graph which is indeed the traditional setting without graphical considerations. 
We do not claim that the empirically assumed graph is the "best" one but it would be convenient, with a given graph, 
to compare our proposed models. 

\begin{figure}
\vspace*{-15mm}
\centering\includegraphics[scale=1,angle=270]{example.eps}
\caption{The assumed response graph ($J=5$), spatial graph ($I=115$) and structure of the corresponding 
	joint precision matrix ($575 \times 575$).}
\label{fig:graph}
\end{figure}

\begin{table}[t]
\caption{\textsc{DIC} of six models}
\label{tab:DIC}
\begin{center}
\begin{tabular}{c|ccc|ccc}
\hline
\hline
 & & Complete Graph $\bfC^{(r)}$ & & & Assumed Graph $\bfC^{(r)}$ & \\
\hline
 & \textsc{Dbar} & \textsc{pD} & \textsc{DIC} & \textsc{Dbar} & \textsc{pD} & \textsc{DIC} \\
\hline 
Model 1 & 4622.4 & 474.7 & 5097.1  & 4622.5 & 471.6 & 5094.1  \\
Model 2 & 4615.4 & 462.0 & 5077.4  & 4598.8 & 457.9 & 5056.7  \\
Model 3 & 4545.8 & 515.8 & 5061.6  & 4541.8 & 514.2 & 5055.9  \\
\hline
\end{tabular}
\end{center}
\end{table}

In total we fitted six models (Model 1-3 versus two graphs) for the data and performed Markov chain Monte Carlo for 100,000 cycles
with 50,000 burn-ins. Convergence and related issues have been carefully examined. We use the deviance information
criterion (\textsc{DIC}) to compare all six models. Note that $\textsc{DIC}=\textsc{Dbar}+\textsc{pD}$, giving the posterior mean deviance
penalized by the effective number of parameters. Since in this hierarchical model, usually the focus is underlying 
prevalences $p_{ij}$ and relative risks $\eta_{ij}$ which are mid-level parameters, \textsc{DIC} should be an appropriate measure 
for this type of comparisons \citep{Spieg:2002}. Comparison results are reported in Table \ref{tab:DIC}.  
When compare the complete graph with the assumed graph, the gain of considering a graphical structure
is substantial even though the assumed graph may not be the "best" one. We are also interested in comparing 
our proposed three models. Clearly, Model 1 is much worse than the other two and the reason, however, 
has been mentioned in Section \ref{sec:spec}: the unnecessary positivity condition seriously limits the precision structures
and the "true" structure may not be in this class. There is not much difference in \textsc{DIC} between Model 2 and 3,
indicating that both models are equally preferable for the data. However, it is interesting to notice the huge
difference in \textsc{Dbar} and \textsc{pD}. Model 3 is much more flexible than Model 2 in the spatial domain and hence
fits the data better (smaller \textsc{Dbar}), but suffers from a large number of parameters (greater \textsc{pD}).  

An intuitive way to see the dependence after modeling is to check posterior correlations of $u_{ij}$, as shown
in Figure 3. We also empirically calculated correlations using $r_{ij}$ for the five response variables 
(upper diagonal and by rows): $(0.07, 0.49, 0.31, 0.27)$, $(0.49, 0.35, 0.49)$, $(0.67, 0.74)$, $(0.63)$. 
The block patterns in Figure 3
generally reflect the empirical correlations. It can be seen that Model 2 and Model 3 are similar
while Model 1 is unsatisfactory. A map of posterior means of $u_{ij}$ under Model 2 is shown in Figure \ref{fig:map}
for illustration. A higher $u_{ij}$ implies a higher prevalence or relative risk. 
Notice that the three disease mortalities have similar spatial patterns, and "smoking" and "no health care" seem to be
positively correlated with those mortalities. 

%\begin{figure}
%\vspace*{-15mm}
%\centering\includegraphics[scale=.5]{corr_matrix.png}
%\caption{Posterior correlations of $u_{ij}$ ($575 \times 575$) under the assumed response graph.}
%\label{fig:corr}
%\end{figure}

\begin{figure}
\vspace*{-15mm}
\centering\includegraphics[scale=1]{map_wishart_gray.eps}
\caption{Posterior means of $u_{ij}$ under Model 2 with the assume response graph.}
\label{fig:map}
\end{figure}

\section{Discussion} \label{sec:dis}
%
The three models discussed in this paper mainly follow the adjacency structure (\ref{adjacency}) which is
constructed based on given response graph and spatial graph. We have already pointed out some other possible
choices that are all "separable". Indeed, the adjacency matrix describes the neighboring structure of the 
joint graph. It could be determined to serve a particular need. For instance, consider a spatio-temporal modeling
where the time domain determines the response graph. In this scenario, $\bfC^{(r)}=\bfC^{(t)}$ could be a directed
acyclic graph with $p$th order indicating a conditional dependence of $u_{it} \mid u_{i,t-1},\ldots,u_{i,t-p}$. 
Note that $\bfC$ will be asymmetric and $\bfB$ will be asymmetric as well. 
The construction (\ref{adjacency}) could also be used if the dynamic process is separable, otherwise more 
complicated constructions can be considered, for instance, time-varying spatial structures. 

Parameterization for $\bfB$ and $\bfT$ perhaps is the most challenging part for it is always associated with
the positivity condition. Under an arbitrary parameterization, it is difficult to seek a suitable condition for the positive
definiteness. As we have seen in Model 1, a sufficient condition seriously limits the richness of the precision structure,
but a sufficient and necessary condition is not obvious. On the other hand, perhaps the Kronecker product structure
is more straightforward and easier to deal with, i.e. a structure like $\bfOmega^{(r)} \otimes \bfOmega^{(s)}$,
because one only needs to ensure positive definite $\bfOmega^{(r)}$ and $\bfOmega^{(s)}$ separately.  
However, the Kronecker product structure is not always satisfactory due to its inflexibility. One option
is to decompose and modify both $\bfOmega^{(r)}$ and $\bfOmega^{(s)}$ and maintain their positive definiteness
\citep{Gelfand:2003, Mart:2013}, but a meaningful interpretation might be lost. With a somewhat arbitrary parameterization, 
such as Model 1, we are curious about its performance under a suitable condition. 


\section*{Appendix 1: Some Derivations} \label{app:a1}
%
\subsection*{Derivation of equation (\ref{Model1})} 
With the parameterization in Model 1, we have $\bfT=\bfD^{-1}(\bfDelta \otimes \bfI_I)$ and 
\ba
	\bfB=(\bfDelta^{\frac{1}{2}}\otimes\bfI_I)\bfD^{-1}\left\{\bfLambda\otimes\bfC^{(s)}+(\bfPsi\circ\bfC^{(r)})\otimes\bfI_I 
		+ (\bfPhi\circ\bfC^{(r)})\otimes\bfC^{(s)}\right\}(\bfDelta^{-\frac{1}{2}}\otimes\bfI_I).
\ea
Then immediately we have expression (\ref{Model1}) for $\bfT^{-1}(\bfI-\bfB)$. 
%It is positive definite if and only if
%$\bfD-\bfLambda\otimes\bfC^{(s)}-(\bfPsi\circ\bfC^{(r)})\otimes\bfI_I - (\bfPhi\circ\bfC^{(r)})\otimes\bfC^{(s)}$
%is positive definite, which has the following structure
%\ba
%	\begin{pmatrix}
%		\bfD_1-\lambda_1\bfC^{(s)} & -\psi_{12}\bfI_I-\phi_{12}\bfC^{(s)} & \cdots & -\psi_{1J}\bfI_I-\phi_{1J}\bfC^{(s)} \\
%		-\psi_{12}\bfI_I-\phi_{12}\bfC^{(s)} & \bfD_2-\lambda_2\bfC^{(s)} & \cdots & -\psi_{2J}\bfI_I-\phi_{2J}\bfC^{(s)} \\
%		\vdots & \vdots & \ddots & \vdots \\
%		-\psi_{1J}\bfI_I-\phi_{1J}\bfC^{(s)} & -\psi_{2J}\bfI_I-\phi_{2J}\bfC^{(s)} & \cdots & \bfD_J-\lambda_J\bfC^{(s)}
%	\end{pmatrix},
%\ea
%where some blocks are zeros if they are corresponding to missing edges in $\bfC^{(r)}$ and $\bfD=\mbox{Bdiag}(\bfD_1,\ldots,\bfD_J)$.
%For each row, the sum of absolute values of off-diagonal elements is less than the diagonal element if 
%$\mbox{max}\{|\lambda_j|, |\psi_{jj'}|, |\phi_{jj'}|; \forall j, j'\}<1$.

\subsection*{Derivation of equation (\ref{Model2})}
With the parameterization in Model 2, we have $\bfT=\left(\diag(\bfOmega) \otimes \bfD^{(s)}\right)^{-1}$ and
\ba
	\bfB=\rho\bfI_J\otimes\bfD^{(s)-1}\bfC^{(s)}-\left(\diag(\bfOmega)^{-1}\bfOmega\circ\bfC^{(r)}\right)\otimes\bfI_I
		+\rho\left(\diag(\bfOmega)^{-1}\bfOmega\circ\bfC^{(r)}\right)\otimes\bfD^{(s)-1}\bfC^{(s)}.
\ea
Then
\ba
	\bfT^{-1}(\bfI-\bfB) &=& \diag(\bfOmega)\otimes \bfD^{(s)}-\rho\diag(\bfOmega)\otimes\bfC^{(s)}
		+(\bfOmega\circ\bfC^{(r)})\otimes\bfD^{(s)}-\rho(\bfOmega\circ\bfC^{(r)})\otimes\bfC^{(s)} \\
		&=& \left\{\bfOmega\circ(\bfI_J+\bfC^{(r)})\right\}\otimes\bfD^{(s)}- \left\{\bfOmega\circ(\bfI_J+\bfC^{(r)})\right\}\otimes\rho\bfC^{(s)} \\
		&=&  \left\{\bfOmega\circ(\bfI_J+\bfC^{(r)})\right\}\otimes(\bfD^{(s)}-\rho\bfC^{(s)}),
\ea
which is expression (\ref{Model2}).

\subsection*{Derivation of equation (\ref{Model3})}
With the parameterization in Model 3, we have $\bfT=\left(\diag(\bfOmega^{(r)})\otimes\diag(\bfOmega^{(s)})\right)^{-1}$ and
\ba
	\bfB &=& -\bfI_J\otimes\left(\diag(\bfOmega^{(s)})^{-1}\bfOmega^{(s)}\circ\bfC^{(s)}\right)-
		\left(\diag(\bfOmega^{(r)})^{-1}\bfOmega^{(r)}\circ\bfC^{(r)}\right)\otimes\bfI_I \\
		&& -\left(\diag(\bfOmega^{(r)})^{-1}\bfOmega^{(r)}\circ\bfC^{(r)}\right)\otimes
		\left(\diag(\bfOmega^{(s)})^{-1}\bfOmega^{(s)}\circ\bfC^{(s)}\right).
\ea
Then
\ba
	\bfT^{-1}(\bfI-\bfB) &=& \diag(\bfOmega^{(r)})\otimes\diag(\bfOmega^{(s)})
		+ \diag(\bfOmega^{(r)})\otimes(\bfOmega^{(s)}\circ\bfC^{(s)}) \\
		&& + (\bfOmega^{(r)}\circ\bfC^{(r)})\otimes\diag(\bfOmega^{(s)})
			+ (\bfOmega^{(r)}\circ\bfC^{(r)})\otimes(\bfOmega^{(s)}\circ\bfC^{(s)}) \\
		&=& \left\{\bfOmega^{(r)}\circ(\bfI_J+\bfC^{(r)})\right\} \otimes 
			\left\{\bfOmega^{(s)}\circ(\bfI_I+\bfC^{(s)})\right\},
\ea
which is expression (\ref{Model3}).

\section*{Appendix 2: Bayesian Computations} \label{app:a2}
%
\subsection*{An example of hierarchical 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,\ldots,I$ and $j=1,\dots,J$,
\ba
	&& y_{ij} \mid p_{ij} \sim \mbox{Bin}(n_{ij}, p_{ij}), ~~~~~~~~ \mbox{logit}(p_{ij})=\beta_j+u_{ij}, \\
	&& \beta_j \sim \mbox{N}(0, \tau_0^2), ~~~~~~~~ \bfU \mid \bfB,\bfT \sim \textsc{MCAR}(\bfB,\bfT),
\ea
where $\tau_0^2$ is a given constant, $\bfU$ is the matrix-variate of $u_{ij}$,
and priors for $\bfB$ and $\bfT$ depend on the specific parameterization.
This section is organized as follows: we first discuss updating effects parameters $\beta_j$ and $u_{ij}$
and then discuss updating the parameterized $\bfB$ and $\bfT$ for Model 1, 2 and 3 separately. 

\subsection*{Updating effect parameters and hierarchical centering}
Our experience has shown that the convergence is poor if we directly update $\beta_j$ and $u_{ij}$. 
We apply the hierarchical centering technique \citep{Gelfand:1995} and block sampling, which work well
in our practice. Let $\gamma_{ij}=\beta_j+u_{ij}$ and $\bfgamma=\vec[(\gamma_{ij})_{I\times J}]$ has a non-centered 
\textsc{MCAR} prior. We update $(\gamma_{ij},\beta_j)$ instead of $(u_{ij},\beta_j)$. The full conditional distribution of $\gamma_{ij}$ is
\ba
	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}-\beta_j-\sum_{\{i'j'\}\neq \{ij\}} b_{\{ij\},\{i'j'\}}(\gamma_{i'j'}-\beta_{j'}) \right) \right\}.
\ea
We use Metropolis-Hastings algorithm to sample $\gamma_{ij}$ from this conditional density. 
We block sample $\bfbeta$ in the following way. For now denote $\bfA=\bfT^{-1}(\bfI-\bfB)$, the joint precision matrix.
Let $\bfgamma^*=\bfA\bfgamma$ and $\bfgamma^{**}$ be a $J \times 1$ "collapsed" vector of $\bfgamma^*$: $\gamma^{**}_1$
is the sum of the first $I$ elements in $\bfgamma^*$, $\gamma^{**}_2$ is the sum of the second $I$ elements 
in $\bfgamma^*$ and so on. Partition $\bfA$ into $J \times J$ blocks and define matrix $\bfH$ which is the "collapsed" $\bfA$:
\ba
	\bfH=\begin{pmatrix}
		{\bf1}^T\bfA_{11}{\bf1} & \cdots & {\bf1}^T\bfA_{1J}{\bf1} \\
		\vdots & \ddots & \vdots \\
		{\bf1}^T\bfA_{J1}{\bf1} & \cdots & {\bf1}^T\bfA_{JJ}{\bf1},
	\end{pmatrix}
\ea
where ${\bf1}$ is the all-one vector. Then the full conditional distribution for the vector $\bfbeta$ is
\ba
	(\bfbeta \mid \cdot) \sim \mbox{N}\left[ (\bfH+1/\tau_0^2 \bfI)^{-1}\bfgamma^{**}, (\bfH+1/\tau_0^2 \bfI)^{-1} \right].
\ea

\subsection*{Model 1: updating $\delta_j$, $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$}
Recall priors on these parameters: $\delta_j\sim \mbox{IG}(a_j,b_j)$ and $\lambda_j, \psi_{jj'}, \phi_{jj'}\sim \mbox{Unif}(-1,1)$.
Let $\bfu_j$ be the $j$th column vector of $\bfU$, $j=1,\ldots,J$ and $\bfD_j$ be the $j$th diagonal block of $\bfD$.
The full conditional distribution of $\delta_j$ is given by
\ba
	p(\delta_j\mid \cdot)\propto \delta_j^{-\frac{I}{2}-a_j-1} \exp\left\{ -\frac{1}{2\delta_j}\bfu_j^T
	(\bfD_j-\lambda_j\bfC^{(s)})\bfu_j + \sum_{j' \sim j} \frac{1}{\sqrt{\delta_j\delta_{j'}}}
	\bfu_j^T (\psi_{jj'}\bfI_I+\phi_{jj'}\bfC^{(s)}) \bfu_{j'} - \frac{b_j}{\delta_j} \right\}.
\ea
It can be shown that the transformed one $(\sqrt{1/\delta_j}\mid \cdot)$ is log-concave when $I+2a_j-1>0$ (trivially satisfied). 
Thus, we use the adaptive rejection sampling (ARS) to update $\delta_j$.

Let $\bfW$ be an $I \times J$ matrix, where $\vec(\bfW)=(\bfDelta^{-1/2}\otimes \bfI_I)\vec(\bfU)$ and 
$\bfw_j$ be the $j$th column vector of $\bfW$. Let $\bfM=\bfD-\bfLambda\otimes\bfC^{(s)}-(\bfPsi\circ\bfC^{(r)})\otimes\bfI_I 
- (\bfPhi\circ\bfC^{(r)})\otimes\bfC^{(s)}$ as in (\ref{Model1}). Then $\lambda_j$, $\psi_{jj'}$ and $\phi_{jj'}$ are 
sequentially updated through following full conditional distributions,
\ba
	p(\lambda_j \mid \cdot) &\propto& \left|\bfM(\lambda_j) \right|^{\frac{1}{2}}
		\exp\left( \frac{1}{2}\lambda_j\bfw_j^T \bfC^{(s)}\bfw_j \right), \\
	p(\psi_{jj'} \mid \cdot) &\propto& \left|\bfM(\psi_{jj'}) \right|^{\frac{1}{2}}
		\exp\left( \psi_{jj'}\bfw_j^T \bfw_{j'} \right), \\
	p(\phi_{jj'} \mid \cdot) &\propto& \left|\bfM(\phi_{jj'}) \right|^{\frac{1}{2}}
		\exp\left( \phi_{jj'}\bfw_j^T \bfC^{(s)} \bfw_{j'} \right).
\ea
We use Metropolis-Hastings algorithm to update these parameters. Note that evaluating $|\bfM|$ could be computationally
intensive because of its large dimension. However, $\bfM$ could often be sparse because it is a construction from $\bfC^{(r)}$
and $\bfC^{(s)}$. An efficient algorithm, usually based on the Cholesky decomposition, on sparse matrices is helpful in reducing the 
computation time. We recommend the R package "spam" \citep{Furrer:2010} for R users. 

\subsection*{Model 2: updating $\rho$ and $\bfOmega_{\bfC^{(r)}}$}
Recall priors on these parameters: $\rho\sim U(-1,1)$ and $\bfOmega_{\bfC^{(r)}}\sim \mbox{GWis}(b,\bfV)$.
Use Metropolis-Hastings algorithm to update $\rho$. It can be shown that the full conditional distribution for $\rho$ is
\ba
	p(\rho \mid \cdot) \propto \left| \bfD^{(s)}-\rho\bfC^{(s)} \right|^{\frac{J}{2}}
	\exp\left\{ \frac{\rho}{2}\vec(\bfU)^T (\bfOmega_{\bfC^{(r)}} \otimes \bfC^{(s)}) \vec(\bfU) \right\}.
\ea
Efficiency of evaluating $\left| \bfD^{(s)}-\rho\bfC^{(s)} \right|$ depends on the sparsity of the matrix.

Let $\bfW_1$ be an $I \times J$ matrix, where $\vec(\bfW_1)=[\bfI_J \otimes (\bfD^{(s)}-\rho\bfC^{(s)})^{1/2}]\vec(\bfU)$ and 
$\bfw_{1,j}$ be the $j$th column vector of $\bfW_1$. Let $\bfS$ be an $J\times J$ matrix with $s_{jj'}=\bfw_{1,j}^T \bfw_{1,j'}$. 
Then the full conditional distribution of $\bfOmega_{\bfC^{(r)}}$ is 
\ba
	p(\bfOmega_{\bfC^{(r)}} \mid \cdot)\propto \left| \bfOmega_{\bfC^{(r)}} \right|^{\frac{b+I-2}{2}} 
	\exp\left[ -\frac{1}{2}tr\left\{ \bfOmega_{\bfC^{(r)}}(\bfV+\bfS) \right\} \right]
	1_{\bfOmega_{\bfC^{(r)}}\in M^+({\bfC^{(r)}})},
\ea
which is GWis($b+I, \bfV+\bfS$). For sampling from the G-Wishart distribution, we use the block Gibbs sampler
(given the set of maximum cliques) by \cite{Wang:2012}.

\subsection*{Model 3: updating $\bfOmega_{\bfC^{(r)}}$ and $\bfOmega_{\bfC^{(s)}}$}
Recall that we impose a constraint and use a joint prior (\ref{jointprior}) on $(z,z\bfOmega_{\bfC^{(r)}})$ and have 
$\bfOmega_{\bfC^{(s)}}\sim \mbox{GWis}(b^{(s)},\bfV^{(s)})$. Let both $\bfW^{(r)}$ and $\bfW^{(s)}$ be $I \times J$ matrices,
where $\vec(\bfW^{(r)})=(\bfI_J \otimes \bfOmega_{\bfC^{(s)}}^{1/2})\vec(\bfU)$ and 
$\vec(\bfW^{(s)})=(\bfOmega_{\bfC^{(r)}}^{1/2} \otimes \bfI_I)\vec(\bfU)$. Let $\bfw_j^{(r)}$ be the $j$th column vector of $\bfW^{(r)}$
and $\bfw_i^{(s)}$ be the $i$th row vector of $\bfW^{(s)}$. Then let $\bfS^{(r)}$ be $J\times J$ with $s_{jj'}^{(r)}=\bfw_j^{(r)T}\bfw_j^{(r)}$
and $\bfS^{(s)}$ be $I\times I$ with $s_{ii'}^{(s)}=\bfw_i^{(s)T}\bfw_i^{(s)}$. With these notations, we have
\ba
	(z\mid \cdot) \sim Ga(a_z,b_z),
\ea 
where $a_z=J(b-2)/2+\nu(\bfC^{(r)})$ and $b_z=tr(\bfOmega_{\bfC^{(r)}}\bfV^{(r)})/2$;
\ba
	(\bfOmega_{\bfC^{(r)}}\mid \cdot) \sim \mbox{GWis}(b^{(r)}+I, z\bfV^{(r)}+\bfS^{(r)})
\ea
and 
\ba
	(\bfOmega_{\bfC^{(s)}}\mid \cdot) \sim \mbox{GWis}(b^{(s)}+J, \bfV^{(s)}+\bfS^{(s)}).
\ea

main.bbl
\begin{thebibliography}{xx}

\harvarditem{Atay-Kayis \harvardand\ Massam}{2005}{Atay:2005}
Atay-Kayis, A. \harvardand\ Massam, H.  \harvardyearleft 2005\harvardyearright
  , `A {M}onte {C}arlo method for computing the marginal likelihood in
  nondecomposable {G}aussian graphical models', {\em Biometrika} {\bf
  92},~317--335.

\harvarditem[Banerjee et~al.]{Banerjee, Gelfand \harvardand\
  Carlin}{2004}{Banerjee:2004}
Banerjee, S., Gelfand, A. \harvardand\ Carlin, B.  \harvardyearleft
  2004\harvardyearright , {\em Hierarchical modeling and analysis for spatial
  data}, CRC Press/Chapman \& Hall.

\harvarditem{Besag}{1974}{Besag:1974}
Besag, J.  \harvardyearleft 1974\harvardyearright , `Spatial interaction and
  the statistical analysis of lattice systems', {\em Journal of the Royal
  Statistical Society. Series B} {\bf 36},~192--236.

\harvarditem[Dobra et~al.]{Dobra, A. \harvardand\ Rodriguez}{2011}{Dobra:2011}
Dobra, A., A., L. \harvardand\ Rodriguez, A.  \harvardyearleft
  2011\harvardyearright , `Bayesian inference for general {G}aussian graphical
  models with applications to multivariate lattice data', {\em Journal of the
  American Statistical Association} {\bf 106},~1418--1433.

\harvarditem{Furrer \harvardand\ Sain}{2010}{Furrer:2010}
Furrer, R. \harvardand\ Sain, S.  \harvardyearleft 2010\harvardyearright ,
  `spam: A sparse matrix {R} package with emphasis on {MCMC} methods for
  {G}aussian {M}arkov random fields', {\em Journal of Statistical Software}
  {\bf 36},~1--25.

\harvarditem[Gelfand et~al.]{Gelfand, Sahu \harvardand\
  Carlin}{1995}{Gelfand:1995}
Gelfand, A., Sahu, S. \harvardand\ Carlin, B.  \harvardyearleft
  1995\harvardyearright , `Efficient parameterizations for normal linear mixed
  models', {\em Biometrika} {\bf 82},~479--488.

\harvarditem{Gelfand \harvardand\ Vounatsou}{2003}{Gelfand:2003}
Gelfand, A. \harvardand\ Vounatsou, P.  \harvardyearleft 2003\harvardyearright
  , `Proper multivariate conditional autoregressive models for spatial data
  analysis', {\em Biostatistics} {\bf 4},~11--25.

\harvarditem[Jin et~al.]{Jin, Banerjee \harvardand\ Carlin}{2007}{Jin:2007}
Jin, X., Banerjee, S. \harvardand\ Carlin, B.  \harvardyearleft
  2007\harvardyearright , `Order-free co-regionalized areal data models with
  application to multiple-disease mapping', {\em Journal of Royal Statistical
  Society. Series B} {\bf 69},~817--838.

\harvarditem[Jin et~al.]{Jin, Carlin \harvardand\ Banerjee}{2005}{Jin:2005}
Jin, X., Carlin, B. \harvardand\ Banerjee, S.  \harvardyearleft
  2005\harvardyearright , `Generalized hierarchical multivariate {CAR} models
  for areal data', {\em Biometrics} {\bf 61},~950--961.

\harvarditem[Kim et~al.]{Kim, Sun \harvardand\ Tsutakawa}{2001}{Kim:2001}
Kim, H., Sun, D. \harvardand\ Tsutakawa, R.  \harvardyearleft
  2001\harvardyearright , `A bivariate {B}ayes method for improving estimates
  of mortality rates with a twofold conditional autoregressive model', {\em
  Journal of the American Statistical Association} {\bf 96},~1506--1521.

\harvarditem{Letac \harvardand\ Massam}{2007}{Letac:2007}
Letac, G. \harvardand\ Massam, H.  \harvardyearleft 2007\harvardyearright ,
  `Wishart distributions for decomposable graphs', {\em Annals of Statistics}
  {\bf 35},~1278--1323.

\harvarditem{MacNab}{2011}{MacNab:2011}
MacNab, Y.  \harvardyearleft 2011\harvardyearright , `On {G}aussian {M}arkov
  random fields and {B}ayesian disease mapping', {\em Statistical Methods in
  Medical Research} {\bf 20},~49--68.

\harvarditem{Mardia}{1988}{Mardia:1988}
Mardia, K.  \harvardyearleft 1988\harvardyearright , `Multidimentional
  multivariate {G}aussian {M}arkov random fields with application to image
  processing', {\em Journal of Multivariate Analysis} {\bf 24},~265--284.

\harvarditem{Martinez-Beneito}{2013}{Mart:2013}
Martinez-Beneito, M.  \harvardyearleft 2013\harvardyearright , `A general
  modeling framework for multivariate disease mapping', {\em Biometrika} {\bf
  100},~539--553.

\harvarditem[Sain et~al.]{Sain, Furrer \harvardand\ Cressie}{2011}{Sain:2011}
Sain, S., Furrer, R. \harvardand\ Cressie, N.  \harvardyearleft
  2011\harvardyearright , `A spatial analysis of multivariate output from
  regional climate models', {\em The Annals of Applied Statistics} {\bf
  5},~150--175.

\harvarditem[Spiegelhalter et~al.]{Spiegelhalter, Best, Carlin \harvardand\
  van~der Linde}{2002}{Spieg:2002}
Spiegelhalter, D., Best, N., Carlin, B. \harvardand\ van~der Linde, A.
  \harvardyearleft 2002\harvardyearright , `Bayesian measures of model
  complexity and fit (with discussion)', {\em Journal of the Royal Statistical
  Society Series B} {\bf 64},~583--639.

\harvarditem{Wang \harvardand\ Li}{2012}{Wang:2012}
Wang, H. \harvardand\ Li, S.  \harvardyearleft 2012\harvardyearright ,
  `Efficient {G}aussian graphical model determination under {G}-{W}ishart prior
  distributions', {\em Electronic Journal of Statistics} {\bf 6},~168--198.

\harvarditem{Wang \harvardand\ West}{2009}{Wang:2009}
Wang, H. \harvardand\ West, M.  \harvardyearleft 2009\harvardyearright ,
  `Bayesian analysis of matrix normal graphical models', {\em Biometrika} {\bf
  96},~821--834.

\end{thebibliography}
