A Spatial Bayesian Semiparametric Mixture Model for Positive Definite Matrices with Applications to Diffusion Tensor Imaging

TL;DR

This paper introduces a Bayesian semiparametric mixture model for analyzing spatially dependent positive definite matrices in diffusion tensor imaging, improving inference by capturing spatial correlations and utilizing matrix-variate data.

Contribution

It extends spatial statistical methods to matrix-variate data, providing a novel Bayesian approach for DTI analysis that accounts for spatial dependence and preserves matrix information.

Findings

The proposed model outperforms non-spatial methods in simulations.

Application to brain data reveals insights into cocaine's effects on brain structure.

The method offers a fast and elegant Bayesian computational framework.

Abstract

Diffusion tensor imaging (DTI) is a popular magnetic resonance imaging technique used to characterize microstructural changes in the brain. DTI studies quantify the diffusion of water molecules in a voxel using an estimated 3x3 symmetric positive definite diffusion tensor matrix. Statistical analysis of DTI data is challenging because the data are positive definite matrices. Matrix-variate information is often summarized by a univariate quantity, such as the fractional anisotropy (FA), leading to a loss of information. Furthermore, DTI analyses often ignore the spatial association of neighboring voxels, which can lead to imprecise estimates. Although the spatial modeling literature is abundant, modeling spatially dependent positive definite matrices is challenging. To mitigate these issues, we propose a matrix-variate Bayesian semiparametric mixture model, where the positive definite matrices are distributed as a mixture of inverse Wishart distributions with the spatial dependence captured by a Markov model for the mixture component labels. Conjugacy and the double Metropolis-Hastings algorithm result in fast and elegant Bayesian computing. Our simulation study shows that the proposed method is more powerful than non-spatial methods. We also apply the proposed method to investigate the effect of cocaine use on brain structure. The contribution of our work is to provide a novel statistical inference tool for DTI analysis by extending spatial statistics to matrix-variate data.

Tables

Table 1. Table 1 : The simulation results. The true positive rate, false positive rate, false discovery rate, and typical computation time of the Potts Model and Random Ellipsoid Model are summarized.

Method	Potts			Random Ellipsoid
	K=10	K=50	K=100
TPR	0.99	0.97	0.98	0.29
FPR	0.013	0.010	0.009	0.00
FDR	0.025	0.025	0.024	0.002
Time (hours)	0.5	0.8	1.0	0.01

Table 2. (a) Data generated from mixture models

Method	Potts			Random Ellipsoid
	K=10	K=50	K=100
TPR	0.99	0.97	0.98	0.29
FPR	0.013	0.010	0.009	0.00
FDR	0.025	0.025	0.024	0.002
Time (hours)	0.5	0.8	1.0	0.01

Table 3. (b) Data generated from the spatial Cholesky process

Method	Potts			Random Ellipsoid
	K=10	K=50	K=100
TPR	0.79	0.80	0.77	0.50
FPR	0.01	0.01	0.01	0.00
FDR	0.03	0.03	0.02	0.03
Time (hours)	1.0	1.5	1.8	0.01

Table 4. Table 2 : The Rand index for measuring clustering similarities. The off-diagonals of the table are the Rand indices for any two K 𝐾 K .

K	100	200	300	400	500
100	.	0.92	0.91	0.92	0.92
200	.	.	0.94	0.95	0.94
300	.	.	.	0.96	0.95
400	.	.	.	.	0.96
500	.	.	.	.	.

Equations

$${\mathbf{A}}_v|{\mathbf{M}}_v,m\stackrel{indep.}{\sim }{\mathcal{I}\mathcal{W}}_p({\mathbf{M}}_v,m),$$

$${\mathcal{P}}_k=P(g_v=k|\beta ,{\eta }_k,{\mathbf{g}}_{-v})\propto \exp \left[{\eta }_k+\beta \sum _{u\in N_v}\mathcal{I}(g_u=k)\right],$$

$$\sum _{k=1}^K{\mathcal{P}}_k{\mathcal{I}\mathcal{W}}_p({\mathbf{A}}_v|{\mathbf{V}}_k,m),$$

$$\mathcal{V}(u,v)=\mathbb{E}||{\mathbf{A}}_u-{\mathbf{A}}_v|{|}_F^2=\gamma (m,\nu ,\mathbf{\Sigma })P(u,v|\beta ,\xi ),$$

$${\mathbf{A}}_{iv}|{\mathbf{M}}_{iv},m\stackrel{indep.}{\sim }{\mathcal{I}\mathcal{W}}_p({\mathbf{M}}_{iv},m),{\mathbf{M}}_{iv}:={\mathbf{V}}_{g_{iv}},{\mathbf{V}}_k\stackrel{i.i.d}{\sim }{\mathcal{W}}_p(\mathbf{\Sigma },\nu ).$$

$$P(g_{iv}=k|\alpha ,\beta ,\xi ,{\mathbf{g}}_{-iv},\mathbf{h})$$

$$P(h_{xv}=k|\alpha ,\beta ,{\mathbf{h}}_{-xv},\mathbf{g})$$

$$\sum _{k=1}^K\exp \left[-k^{\xi }+\alpha \mathcal{I}(h_{xv}=k)\right]{\mathcal{I}\mathcal{W}}_p({\mathbf{A}}_{iv}|{\mathbf{V}}_k,m),$$

$${\mathcal{H}}_{ov}:h_{0v}=h_{1v}$$

$${\mathcal{H}}_{av}:h_{0v}\ne h_{1v}.$$

$${\mathcal{V}}_{ij}(u,v)=\mathbb{E}||{\mathbf{A}}_{iu}-{\mathbf{A}}_{jv}|{|}_F^2=\gamma (m,\nu ,\mathbf{\Sigma })P_{ij}(u,v|\alpha ,\beta ,\xi )$$

$${\mathbf{A}}_{iv}|{\mathbf{M}}_{iv},m\stackrel{indep.}{\sim }{\mathcal{I}\mathcal{W}}_p({\mathbf{M}}_{iv},m),{\mathbf{M}}_{iv}:={\mathbf{V}}_{g_{iv}},{\mathbf{V}}_k\stackrel{i.i.d}{\sim }{\mathcal{W}}_p(\mathbf{\Sigma },\nu )$$

$$P(g_{iv}=k|\alpha ,\beta ,\xi ,{\mathbf{g}}_{-(iv)},\mathbf{h})\propto \exp \left[-k^{\xi }+\beta \sum _{u\in N_v}\mathcal{I}(g_{iu}=k)+\alpha \mathcal{I}(h_{x_{i}v}=k)\right]$$

$$P(h_{xv}=k|\alpha ,\beta ,{\mathbf{h}}_{-(xv)},\mathbf{g})\propto \exp \left[\beta \sum _{u\in N_v}\mathcal{I}(h_{xu}=k)+\sum _{j:x_j=x}\alpha \mathcal{I}(g_{jv}=k)\right]$$

$${\mathbf{\Sigma }}_k\sim {\mathcal{W}}_3((k+1){\mathbf{I}}_3,30)\text{ for k=1,2,3,4};{\mathbf{\Sigma }}_5\sim {\mathcal{W}}_3(1.5{\mathbf{I}}_3,30).$$

$$f(\mathbf{X}|\mathbf{V},n)=\frac{1}{2^{np/2}{\left|\mathbf{V}/n\right|}^{n/2}{\Gamma }_p\left(\frac{n}{2}\right)}{\left|\mathbf{X}\right|}^{(n-p-1)/2}{e}^{-(1/2)\mathrm{tr}({[\mathbf{V}/n]}^{-1}\mathbf{X})}$$

$$f(\mathbf{X}|\mathbf{\Psi },\nu )=\frac{{\left|(\nu -p-1)\mathbf{\Psi }\right|}^{\nu /2}}{2^{\nu p/2}{\Gamma }_p(\frac{\nu }{2})}{\left|\mathbf{X}\right|}^{-(\nu +p+1)/2}{e}^{-\frac{1}{2}\mathrm{tr}((\nu -p-1)\mathbf{\Psi }{\mathbf{X}}^{-1})}$$

$$\mathbb{E}||{\mathbf{A}}_u-{\mathbf{A}}_v|{|}_F^2$$

$$=Tr(\mathbb{E}{\mathbf{A}}_u{\mathbf{A}}_u)+Tr(\mathbb{E}{\mathbf{A}}_v{\mathbf{A}}_v)-2Tr(\mathbb{E}{\mathbf{A}}_u{\mathbf{A}}_v)$$

$$=[Tr(\mathbb{E}{\mathbf{A}}_u{\mathbf{A}}_u)+Tr(\mathbb{E}{\mathbf{A}}_v{\mathbf{A}}_v)-2Tr(\mathbb{E}[{\mathbf{A}}_u{\mathbf{A}}_v|g_u\ne g_v])]\times P(g_u\ne g_v|\beta )$$

$${\mathbf{A}}_u{\mathbf{A}}_u:=[{\mathbf{A}}_u{\mathbf{A}}_v|g_u=g_v]$$

$$=\gamma (m,\nu ,\mathbf{\Sigma })P(g_u\ne g_v|\beta )$$

$$Tr(\mathbb{E}[{\mathbf{A}}_u{\mathbf{A}}_v|g_u\ne g_v])$$

$$=\mathbb{E}[Tr(\mathbb{E}[{\mathbf{A}}_v|{\mathbf{M}}_{g_v}]\mathbb{E}[{\mathbf{A}}_u|{\mathbf{M}}_{g_u}])]$$

$$=\sum _{i=1}^p{\sigma }_i(\mathbf{\Sigma }{)}^2=\lambda$$

$$\mathbb{E}[{\mathbf{A}}_u{\mathbf{A}}_u|{\mathbf{M}}_{g_u}={\mathbf{M}}_{g_v}=\mathbf{M}]$$

$$Tr(\mathbb{E}{\mathbf{A}}_u{\mathbf{A}}_u)=Tr(\mathbb{E}{\lambda }^{\ast })$$

$$+c_2(2{\nu }^{-1}\sum _{i=1}^p{\sigma }_i(\mathbf{\Sigma }{)}^2+\sum _{i=1}^p{\sigma }_i(\mathbf{\Sigma })\sum _{i=1}^p{\sigma }_i(\mathbf{\Sigma })))\}(m-p-1{)}^2$$

$$=\sigma$$

$$P(\mathbf{h},\mathbf{g})\propto \exp \left[\sum _{i=1}^N\sum _{v=1}^n\alpha \mathcal{I}(g_{iv}=h_{x_{i}v})+\sum _{x=0}^1\sum _{u\sim v}\beta \mathcal{I}(h_{xu}=h_{xv})+\sum _{i=1}^N\sum _{u\sim v}\left(\beta \mathcal{I}(g_{iu}=g_{iv})-g_{iu}^{\xi }\right)\right]$$

$$=\exp (U(\mathbf{g},\mathbf{h},\mathbf{\theta }))$$

$$\left[{\mathbf{A}}_{iv}|\{g_u:u\in N_v\},h_{x_{i}v},\{{\mathbf{V}}_k:k\},m\right]=\sum _{k=1}^{K}P(g_{iv}=k|.){\mathcal{I}\mathcal{W}}_p({\mathbf{V}}_k,m)$$

$$=\sum _{k=1}^K\underset{\begin{array}{c}Normalizing\\ Constant\end{array}}{\underbrace{C_1}}\exp \left[-k^{\xi }+\beta \sum _{u\in N_v}\mathcal{I}(g_{iu}=k)+\alpha \mathcal{I}(h_{x_{i}v}=k)\right]{\mathcal{I}\mathcal{W}}_p({\mathbf{V}}_k,m)$$

$$Q=\sum _{\begin{array}{c}{g}_{iu}=1:K,\\ u\in N_v\end{array}}\exp \left[\beta \sum _{u\in N_v}\mathcal{I}(g_{iu}=k)\right]\underset{\begin{array}{c}jointp.m.f\\ of{g}_{iu},u\in N_v\end{array}}{\underbrace{{\mathcal{P}}_{g_{iu},u\in N_v}}}$$

$$\left[{\mathbf{A}}_{iv}|h_{x_{i}v},\{{\mathbf{V}}_k:k\},m\right]=\sum _{k=1}^K\underset{{\Phi }_k(h_{x_{i}v})}{\underbrace{C_{1}Q\exp \left[-k^{\xi }+\alpha \mathcal{I}(h_{x_{i}v}=k)\right]}}{\mathcal{I}\mathcal{W}}_p({\mathbf{V}}_k,m)$$

Full text

Keywords: Bayesian semiparametric, Diffusion tensor imaging, Inverse Wishart distribution, Matrix-variate, Positive definite matrix, Spatial statistics

1 Introduction

Measurement of signal attenuation from water diffusion, often considered one of the most important magnetic resonance contrast mechanisms (Alexander et al., 2007), is usually achieved via diffusion tensor imaging (DTI) that maps and characterizes the 3-D diffusion of water molecules as a function of the spatial location (Basser et al., 1994). The diffusion process in the brain reflects interactions with many obstacles, such as fibers, thereby revealing microscopic details about the underlying tissue architecture. Unlike ordinary images where scalars are summarized for each voxel, a distinguishing feature of DTI is each voxel is associated with a $3\times 3$ symmetric positive definite matrices which can be interpreted as the covariance matrix of a 3-D Gaussian distribution modeling the local Brownian motion of the water molecules (Schwartzman et al., 2008). These positive definite matrices are also called the diffusion tensors (DTs). One important clinical application of the DTI is to detect regions of local differences in the brain between two groups (i.e., normal versus disease), revealing anatomical structural differences (Lo et al., 2010). For example, the motivating data for this paper comes from a clinical DTI study (Ma et al., 2017), where the scientific objective is to detect regions of differences between cocaine users and non-cocaine users.

Statistical analysis of DTI data is challenging due to the difficulty of modeling matrix-variate responses. One option is to project the DTs into fractional anisotropy (FA), a scalar describing the degrees of anisotropy of a diffusion process. However, some information is lost because different positive definite matrices may produce the same FA (Ennis and Kindlmann, 2006). Matrix-variate methods potentially avoid information loss. There are relatively few matrix-variate methods available to analyze DTI data, and they can be broadly classified into the (inverse) Wishart matrix methods (Dryden et al., 2009) and the random ellipsoid models (Schwartzman et al., 2008; Lee and Schwartzman, 2017). However, these voxel-level models ignore information from neighboring voxels that may have similar neuronal activity (see Figure 1), despite recommendations of incorporating this non-negligible spatial association to achieve efficient and valid inference (Spence et al., 2007) as well as studies revealing that the disease status at proximally-located/neighboring voxels can be similar (see Wu et al., 2013; Xue et al., 2018). This motivates us to develop an improved spatial statistical model which (a) utilizes full matrix information, (b) captures spatial dependence, and (c) can be implemented via fast and elegant computing.

The spatial neuroimaging toolbox for univariate responses is considerably rich: Woolrich et al. (2004b) proposed a fully Bayesian model for spatiotemporal imaging data; Kang et al. (2011) implemented spatial point processes for meta-analysis of imaging data; To select essential biological features, Musgrove et al. (2016) introduced spatial Bayesian variable selection for neuroimaging data; Recently, Reich et al. (2018) proposed spectral methods for ADNI data to provide computational benefits. All of these studies demonstrated an improvement in the precision of estimates by properly accounting for spatial dependence.

In this vein, the Potts model, a generalization of the Ising model in statistical mechanics, has also been successfully applied to imaging (Johnson et al., 2013; Li et al., 2018). A desirable property of the Potts model is that it avoids smoothing over abrupt changes in the image intensity (Johnson et al., 2013), and this makes it more attractive than available Gaussian kernel methods. To this end, we assume the positive definite DTs follow a mixture of inverse Wishart distributions, with the mixture component labels modeled via a (spatial) Potts model, representing a discrete Markov random field. This semiparametric mixture specification refers to a class of flexible mixture distributions with a finite number of components (Lindsay and Lesperance, 1995).

Besides spatial modeling, another important topic in neuroimaging is detecting regions of differences between two groups. Previous attempts at identifying regions of differences between two groups were formulated through voxel-wise hypothesis testing (Schwartzman et al., 2008; Lee and Schwartzman, 2017). An alternative option is to construct multilevel hierarchical modeling accounting both subjective-level and group-level variation and use the group-level parameters for voxel-wise hypothesis testing (Woolrich et al., 2004a; Liu et al., 2014). In this paper, we use the latter approach via extending the latent classic Potts model into a hierarchical two-way framework, allowing hypothesis testing via group-level parameters and inter-subject variability simultaneously.

Our proposal is implemented using the Bayesian approach, accounting for the uncertainty of model parameters in all levels of the hierarchy. However, the Bayesian approach is often problematic in neuroimaging because of its heavy computational burden (Cohen et al., 2017). Although the associated Markov chain Monte Carlo (MCMC) algorithm is mostly composed of computationally tractable Gibbs steps that can be paralleled, a major drawback of the Potts model is the intractable normalizing constant, creating a bottleneck for hyperparameter updates. In this paper, it is resolved via the double Metropolis-Hastings algorithm (Liang, 2010) recommended by Park and Haran (2018).

To the best of our knowledge, this is the first work on exploring spatial associations in modeling positive definite matrix-variate data under a Bayesian semiparametric framework, with applications to DTI. In the rest of the paper, we first introduce the single-subject and multi-subject model, and the group hypothesis testing framework in Section 2. Relevant MCMC computational details appear in Section 3. To demonstrate the improvement in performance compared to plausible alternatives, we perform simulation studies in Section 4. In Section 5, we present the application to the motivating cocaine data set. Finally, Section 6 concludes with a discussion.

2 Model

In this section, we introduce the spatial Bayesian semiparametric mixture model for positive definite matrices. We introduce the single-subject model first in Section 2.1 and then extend to the multi-subject model in Section 2.2.

2.1 Single-Subject Model

Let $\bm{A}_{v}$ be the $p\times p$ DT at voxel $v\in\{1,2,..,n\}$. To ensure $\bm{A}_{v}$ is symmetric and positive definite, it is usually parameterized as a (inverse) Wishart matrix (Dryden et al., 2009) or a Gaussian symmetric matrix-variate distribution (Schwartzman et al., 2008). In this paper, we assume that $\bm{A}_{v}$ follows an inverse Wishart distribution as

[TABLE]

where $\mathcal{IW}_{p}(\bm{M}_{v},m)$ is the inverse Wishart distribution parameterized (Appendix A) to have mean $\bm{M}_{v}$ and degrees of freedom $m>p+1$, and the DTs are independently distributed across $v$ given the mean matrices $\bm{M}_{v}$ and the degrees of freedom $m$. The mean matrices are modeled as a finite mixture of Wishart distributions, denoted as $[\bm{M}_{v}|g_{v}=k]:=\bm{V}_{k}$ where $g_{v}\in\{1,2,...,K\}$ is the latent cluster label. The prior of $\bm{V}_{k}$ is $\bm{V}_{k}\mathbin{\overset{i.i.d}{\kern 0.0pt\sim}}\mathcal{W}_{p}(\bm{\Sigma},\nu)$ where $\mathcal{W}_{p}(\bm{\Sigma},\nu)$ is the Wishart distribution parameterized (Appendix A) to have mean $\bm{\Sigma}$ and degrees of freedom $\nu>p$.

Spatial dependence of the DTs is achieved through the dependence of the mean matrices $\bm{M}_{v}$. We induce spatial dependence via the latent cluster labels that follow a weighted Potts model, specified via the full conditional distributions:

[TABLE]

where $\bm{g}_{-v}$ is the full set $\bm{g}=\{g_{1},g_{2},...,g_{n}\}$ excluding $g_{v}$, $N_{v}$ is a set of indices of the neighboring voxels of $v$, and $\mathcal{I}[\mathcal{A}]=1$ if event $\mathcal{A}$ is true and $\mathcal{I}[\mathcal{A}]=0$ otherwise. Given $\bm{g}_{-v}$ but marginal over $g_{v}$, the distribution of $\bm{A}_{v}$ is the mixture of $K$ inverse Wishart distributions:

[TABLE]

where $\mathcal{IW}_{p}(\bm{A}|\bm{V},m)$ is the inverse Wishart density function of $\bm{A}$ with the mean matrix $\bm{V}$ and the degrees of freedom $m$. Therefore, this semiparametric mixture model spans a rich class of density functions.

Via the Potts model, an image can be considered as a network whose nodes are the voxels. In this network, every voxel is connected to its neighboring voxels. The full conditional distribution of $g_{v}$ depends only on the voxels in the neighboring set $N_{v}$ and therefore the process is Markovian. Since the spatial parameter $\beta$ is the coefficient of the neighboring term $\sum_{u\in N_{v}}\mathcal{I}(g_{u}=k)$, the spatial parameter $\beta$ controls the dependence on the neighboring voxels.

Unlike the classic Potts model (Wu, 1982), the terms $\eta_{k}$ are added as offset terms controlling the overall mass on each cluster. We set $\eta_{k}=-k^{\xi}$ so the parameter $\xi>0$ is the concentration parameter controlling the homogeneity of the latent cluster labels. It is problematic to pre-specify the number of components $K$ in a mixture model (McCullagh et al., 2008) but the offset terms provide more weight on the key components and fewer weights on the trivial components. We fit the model by setting $K$ to be an upper bound on the number of active clusters and allow the data to determine the number of active clusters via estimation of $\xi$: if $\xi\rightarrow 0$, there are several active clusters; if $\xi$ is large, there are a few active clusters. As a result, the model is less sensitive to the number of components $K$ when the offset term $\eta_{k}$ is included. This claim is verified in the simulation studies (Section 4) and the real data application (Section 5), where we find similar results for different $K$.

Quantifying spatial dependence is a vital issue in spatial statistics and neuroimaging. Since this model is for matrix-variate data, we use the expected squared Frobenius norm to measure dependence. The dependence between matrices $\bm{A}$ and $\bm{B}$ can be summarized as $\mathbb{E}||\bm{A}-\bm{B}||_{F}^{2}=\mathbb{E}Tr[(\bm{A}-\bm{B})^{T}(\bm{A}-\bm{B})]$. The norm increases as dependence decreases. If $\bm{A}$ and $\bm{B}$ are $1\times 1$, the expected squared Frobenius norm is the classic variogram (Cressie, 1992) of spatial statistics. In this regard, the expected squared Frobenius norm can be treated as the variogram for matrix-variate data and useful in measuring spatial dependence. In the rest of the paper, we simply call the expected squared Frobenius norm the variogram.

For the Potts model described above, the variogram is

[TABLE]

where $P(u,v|\beta,\xi)$ is the marginal (over all other cluster labels $\bm{g}$) probability of $g_{u}\not=g_{v}$, and $\gamma(m,\nu,\bm{\Sigma})$ is a measure of the variability in $\bm{A}_{v}|\bm{M}_{v}$ and variability of $\bm{V}_{k}$ across $K$. Therefore, the multivariate spatial dependence structure is separable (Cressie, 1992) in that the dependence is the product of a non-spatial term $\gamma(m,\nu,\bm{\Sigma})$ that controls cross dependence and a spatial term $P(u,v|\beta,\xi)$ that controls spatial dependence.

We give the expression of $\gamma(m,\nu,\bm{\Sigma})$ in Appendix B. When $p=3$ and $\bm{\Sigma}=\bm{I}$, the non-spatial term has an expression that is $\frac{12(m+\nu-4)(2m-7)}{\nu(m-3)(m-6)}$ where $m>6$ and $\nu>3$. Therefore, in this special case, the cross dependence decreases if $m$ or $\nu$ is larger (See Figure 2).

The spatial term $P(u,v|\beta,\xi)$ is intractable and so we use a Monte Carlo approximation to study the function. In Figure 3(a), the function is computed under the scenario that the image is a 1-D grid with $K=100$ and $\xi=0$. The spatial term $P(u,v|\beta,\xi)$ increases with distance and larger spatial parameter $\beta$ leads to the stronger spatial dependence. We also compute the function value with different $K$ in Figure 3(b). We have $\lim_{|u-v|\rightarrow\infty}P(u,v|\beta,\xi)=1-\frac{1}{K}$, where relevant result can be found in studies of extreme value analysis (see Reich and Shaby, 2018). Increasing $K$ leads to smaller spatial dependence. Hence, we fix $K$ to be large to eliminate long-range dependence (i.e., $P(u,v|\beta,\xi)<1$ for large $|u-v|$) and estimate $\beta$ to capture local dependence.

For a more intuitive understanding of this model, we also simulate the DTs and visualize the DTs as ellipsoids in a $40\times 40$ grid. In these simulations, we use $\xi=0$, $m=4$, and $\nu=30$. In Figure 4(a), the DTs within the same latent cluster label are similar to each other, indicating that spatial dependence of the DTs can be achieved by the latent cluster labels following the Potts model. In Figure 4(b), larger spatial parameter $\beta$ leads to a realization with more dependence on their neighbors. Figure 4 also illustrates that the Potts model allows sharp breaks, which is desirable if neighboring voxels are in different tracts.

2.2 Muti-Subject Model

Motivated by the cocaine users data set (Ma et al., 2017) that includes 11 cocaine users and 11 non-cocaine users, we extend the single-subject model to the multi-subject setting. The clinical objective is to analyze the brain’s physical structure for differences between the two groups. The objective can be statistically formulated as finding regions in the brain where the distribution of the DTs across the subjects is different between cocaine users and non-cocaine users.

Let $\bm{A}_{iv}$ be the DT and $g_{iv}$ be the cluster label for voxel $v\in\{1,2,...,n\}$ and subject $i\in\{1,2,...,N\}$. By extending $g_{v}$ to $g_{iv}$, the subject-level cluster labels not only model intra-subject spatial dependence but also allow inter-subject variability. As in the single-subject model, the DTs are conditionally independent given the random matrices $\bm{M}_{iv}$ following a finite mixture model:

[TABLE]

However, the latent Potts model is generalized to account for multiple subjects. We define $x_{i}$ as the binary group indicator of subject $i$. In the motivating data, cocaine users have $x_{i}=1$ and non-cocaine users have $x_{i}=0$. To model intra-subject spatial dependence within a group, we extend the latent cluster process by introducing the group-level cluster labels $h_{xv}$ for group $x\in\{0,1\}$ and voxel $v\in\{1,2,...,n\}$. Both $h_{xv}$ and $g_{iv}$ are also spatially dependent with full conditional distributions

[TABLE]

where $\bm{g}_{-(iv)}$ is the set on $\bm{g}_{i}=\{g_{i1},...,g_{in}\}$ excluding $g_{iv}$, $\bm{h}_{-(xv)}$ is the set $\bm{h}_{x}=\{h_{x1},...,h_{xn}\}$ excluding ${h}_{xv}$, $\bm{g}$ is the set on $\{\bm{g}_{1},...,\bm{g}_{N}\}$, and $\bm{h}$ is the set on $\{\bm{h}_{0},\bm{h}_{1}\}$. The joint probability mass function (PMF) of $\{\bm{g}_{1},\bm{g}_{2},...,\bm{g}_{N}\}\cup\{\bm{h}_{0},\bm{h}_{1}\}$ is given in Appendix C. Since the conditional densities in (4) satisfy the conditions of the Hammersley-Clifford Theorem (Clifford, 1990), the existence of joint distribution of $\{\bm{g}_{1},\bm{g}_{2},...,\bm{g}_{N}\}\cup\{\bm{h}_{0},\bm{h}_{1}\}$ is guaranteed (Appendix C).

A graphical representation of this latent Potts model is provided in Figure 5. Cluster labels $\bm{g}_{i}$ and $\bm{h}_{x}$ can be understood as the spatial pattern of subject $i$ and the general spatial pattern of subjects in group $x$, respectively. In comparison to the single-subject model, the group-clustering parameter $\alpha$ is introduced for modeling multiple subjects. If $\alpha=0$, $\bm{A}_{iv}$ is independently distributed over subjects; otherwise, the subject-level cluster label $g_{iv}$ depends on the group-level cluster label $h_{x_{i}v}$, leading to the smaller inter-subject variability of spatial dependence pattern within one group.

To further understand the role of $\alpha$ and $h_{xv}$, we inspect the density of $\bm{A}_{iv}$ conditioned on $h_{xv}$ and marginal over all other labels (Appendix D). The conditional density of $\bm{A}_{iv}$ given $x_{i}=x$ and $h_{xv}$ is the mixture of inverse Wishart distributions proportional to

[TABLE]

where the term $\exp\left[-k^{\xi}+\alpha\mathcal{I}(h_{xv}=k)\right]$ is the proportional weight of cluster $k$. Therefore, $h_{xv}=k$ elevates the mass on mixture component $k$ at voxel $v$ for all subjects with $x_{i}=x$. Assuming $\alpha>0$, the conditional density (5) depends on $x_{i}$ if and only if $h_{0v}\neq h_{1v}$.

The clinical objective is to find regions of differences between two groups, which can be formulated into finding voxels for which the distribution of $\bm{A}_{iv}$ is different for $x_{i}=0$ or $x_{i}=1$. As shown in the conditional density (5), the test can be simplified to the test that

[TABLE]

Bayesian inference provides estimates of the posterior probabilities of the hypotheses, which is further discussed in Section 3.

We investigate spatial dependence as in the single-subject model. To measure spatial dependence within and across subjects, we propose the variogram

[TABLE]

with $i=j$ for individual variogram and $i\not=j$ for inter-subject variogram, respectively. For inter-subject variogram, we also compare the within-group variogram for subjects with $x_{i}=x_{j}$ and between-group variogram for subjects with $x_{i}\neq x_{j}$. Both individual variogram and inter-subject variogram are also separable (Cressie, 1992). $\gamma(m,\nu,\bm{\Sigma})$ is the non-spatial term which has been discussed in Section 2.1. The spatial term $P_{ij}(u,v|\alpha,\beta,\xi)$ is the marginal probability of $g_{iu}\not=g_{jv}$.

In Figure 6(a), the function $P_{ii}(u,v|\alpha,\beta,\xi)$ is computed using a Monte Carlo approximation under the scenario that the image is a 1-D grid with $N=5$ and $K=100$. The spatial parameter $\beta$ largely controls the within-subject dependence. In Figure 6(b) plotting $P_{ij}(u,v|\alpha,\beta,\xi)$, larger $\alpha$ leads to more dependence in the within-group variogram. Therefore, $\alpha$ controls the dependence of subjects within one group. Since $g_{iu}$ and $g_{ju}$ are assumed to be independent, the spatial term $P_{ij}(u,v|\alpha,\beta,\xi)$ is a constant $1-\frac{1}{K}$ in between-group variogram (Figure 6(c)).

3 Computation

We use MCMC to fit the model described in Section 2 (Appendix E). The codes are written in hybrid R and C++ codes. The final model is

[TABLE]

which is referred to as the Potts Model in the rest of the paper. Using the moment method (Robert, 2007)[Section 3.2.4], we set $\bm{\Sigma}$ as the sample mean of all observed DTs. The priori information brought by $\bm{\Sigma}$ has a little impact if the number of observations is large. We put uniform prior for the degrees of freedom $m$ and $\nu$ on $[5,50]\times[4,50]$. Following Liang (2010), we put a uniform prior for $\bm{\theta}=\{\alpha,\beta,\xi\}$ on $[0,20]\times[0,20]\times[0,1]$, denoted as $\pi(\bm{\theta})$. Below we describe the updating rule for each parameter.

The MCMC algorithm is a combination of Gibbs and Metropolis-Hastings steps. The latent mean matrices and cluster labels are updated via Gibbs steps. Their full conditional distributions are

• •

$\bm{V}_{k}|.\sim\mathcal{W}_{p}((\bm{\Sigma}^{-1}\nu+(m-p-1)\sum_{i}\sum_{v:g_{iv}=k}\bm{A}_{iv}^{-1})^{-1}(Nn_{k}m+\nu),Nn_{k}m+\nu)$

• •

$P(g_{iv}=k|.)\propto\mathcal{IW}_{p}(\bm{A}_{iv}|\bm{V}_{g_{iv}},m)\exp\left[-k^{\xi}+\beta\sum_{u\in N_{v}}\mathcal{I}(g_{iu}=k)+\alpha\mathcal{I}(h_{x_{i},v}=k)\right]$

• •

$P(h_{xv}=k|.)\propto\exp\left[\beta\sum_{u\in N_{v}}\mathcal{I}(h_{xu}=k)+\alpha\sum_{j:x_{j}=x}\mathcal{I}(g_{jv}=k)\right]$

where $n_{k}=\sum_{i,v}\mathcal{I}(g_{iv}=k)$. In addition, $P(g_{iv}=k|.)$ and $P(h_{xv}=k|.)$ can be updated in parallel over $i$ and $x$, respectively. Since the uniform prior is not conjugate, we have to sample $[m|.]$ and $[\nu|.]$ via Metropolis-Hastings sampling with log-normal random walk as proposal distribution.

To select regions of differences via Bayesian hypothesis testing, we reject the null hypothesis in (6) if $P(h_{0v}\not=h_{1v}|.)<P(h_{0v}=h_{1v}|.)$. The posterior probabilities can be estimated through MCMC samples that $h_{0v}=h_{1v}$ or $h_{0v}\not=h_{1v}$.

Updating the Potts hyperparameters $\alpha$, $\beta$, and $\xi$ is problematic because the normalizing constant in the joint distribution function of the cluster labels is intractable (see the joint PMF in Appendix C). A simple approach is to estimate the parameters outside of MCMC. The plug-in values can be obtained from cross-validation (i.e., Goldsmith et al., 2014), pseudo-likelihood comparison (i.e., Zhao et al., 2014; Lan et al., 2016), or by comparing empirical and model-based variograms (i.e., Figure 7 and 9). However, these methods fail to account for uncertainty about these imputed parameters and so we update them using the double Metropolis-Hastings algorithm (Liang, 2010).

Park and Haran (2018) review several Monte Carlo methods for models with intractable normalizing constants and recommend the double Metropolis-Hastings algorithm proposed by Liang (2010) because of its ease of implementation and computational efficiency. Li et al. (2018) combine the double Metropolis–Hastings algorithm with usual Bayesian tools for implementing the Potts model. The double Metropolis-Hastings update for $\bm{\theta}$ begins with a candidate $\bm{\theta}^{\prime}$ drawn from $q(\bm{\theta}|\widetilde{\bm{\theta}})$ where $\widetilde{\bm{\theta}}$ is the current value and $q(\bm{\theta}|\widetilde{\bm{\theta}})$ is a log-normal random walk transitional probability centered at $\widetilde{\bm{\theta}}$. Given the candidate $\bm{\theta}^{\prime}$, we draw labels $\bm{g}_{i}^{\prime}=\{g_{iv}^{\prime},...,g_{in}^{\prime}\}$ and $\bm{h}_{x}^{\prime}=\{h_{xv}^{\prime},...,h_{xv}^{\prime}\}$ using Gibbs sampling for each $i$ and $x$, respectively. The candidate $\bm{\theta}^{\prime}$ is accepted with the probability $\min(1,r)$ where $r=\frac{\pi(\bm{\theta}^{\prime})\mathcal{P}(\bm{g}^{\prime},\bm{h}^{\prime}|\widetilde{\bm{\theta}})\mathcal{P}(\widetilde{\bm{g}},\widetilde{\bm{h}}|{\bm{\theta}^{\prime}})}{\pi(\widetilde{\bm{\theta})}\mathcal{P}(\widetilde{\bm{g}},\widetilde{\bm{h}}|\widetilde{\bm{\theta}})\mathcal{P}(\bm{g}^{\prime},\bm{h}^{\prime}|{\bm{\theta}^{\prime}})}$ where $\mathcal{P}(\bm{g},\bm{h}|{\bm{\theta}})$ is the likelihood of $\{\bm{g}_{1},\bm{g}_{2},...,\bm{g}_{N}\}\cup\{\bm{h}_{0},\bm{h}_{1}\}$ conditioned on $\bm{\theta}$. $\widetilde{\bm{g}}_{i}$ and $\widetilde{\bm{h}}_{x}$ are current values. Due to the concern that the double Metropolis-Hastings algorithm is not an exact sampling (Liang, 2010; Park and Haran, 2018), the MCMC convergence of $\bm{\theta}$ in the simulation studies (Section 4) and the real data analysis (Section 5) are monitored by Heidelberger and Welch’s convergence diagnostic (Heidelberger and Welch, 1981).

4 Simulation

In this section, we illustrate the performance of our method using two simulation studies under different scenarios for synthetic data. We compare our method to the non-spatial DTI inference method the Gaussian symmetric matrix model (Schwartzman et al., 2008) referred to as the Random Ellipsoid Model. The Random Ellipsoid Model is a non-spatial matrix-variate method and assumes that $\bm{A}_{iv}$ follows a Gaussian symmetric random matrix distribution with the probability density function (PDF) as $f(\bm{A}_{iv}|x_{i}=x,\bm{\Sigma}_{v},\sigma^{2})=H(\bm{A}_{iv})\exp\left[\frac{1}{2\sigma^{2}}Tr(2\bm{\Sigma}_{xv}\bm{A}_{iv}-\bm{\Sigma}_{xv}\bm{\Sigma}_{xv})\right]$, where $\bm{\Sigma}_{xv}$ is the DT’s population mean at voxel $v$ of group $x$ and $\sigma^{2}$ is the nuisance parameter. The Random Ellipsoid Model selects regions of differences via testing $\bm{\Sigma}_{0v}=\bm{\Sigma}_{1v}$, where the test statistics are constructed by maximum likelihood estimations. the Potts Model has $8,000$ MCMC samples with $3,000$ discarded as burn-in. Methods are evaluated in terms of true positive rate (TPR), false positive rate (FPR), false discovery rate (FDR), and typical computation time.

We first investigate the performance of our method when the data are generated from a mixture model. We use a $40\times 40$ grid with spacing 1 between adjacent grid points as an image. Each simulated data set consists of 5 subjects in the control group ($x_{i}=0$) and 5 subjects in the treatment group ($x_{i}=1$). For the control group ($x_{i}=0$), we equally partition the graph into 4 parts by rectangular regions so that $g_{iv}\in\{1,2,3,4\}$, ordered by right-to-left. Thus each region is a $40\times 10$ region. The treatment group has the same partition as the control group, except a $10\times 10$ region at the middle of the second region where $g_{iv}=5$. This simulates the brain with a small region of difference between the two groups. For each simulation, $\bm{\Sigma}_{k}$ is generated based on the model

[TABLE]

Given the simulated $\bm{\Sigma}_{k}$, the data are generated based on the model $\bm{A}_{iv}|\bm{\Sigma}_{g_{iv}}\sim\mathcal{IW}_{3}(\bm{\Sigma}_{g_{iv}},5)$.

Our model with $K=10,50,100$ is compared to the Random Ellipsoid Model. The results averaged over 50 data sets are summarized in Table 1(a). Our model has significantly improved performance in terms of the TPR, FPR, FDR in comparison to the Random Ellipsoid Model. The small number of subjects might be one of the causes that the alternative produces low TPR. Schwartzman et al. (2008) discusses that the accuracy of maximum likelihood estimates of the Random Ellipsoid Model is dependent on the number of subjects. Since the choice of $K$ does not affect selection accuracy, the simulation results also support the claim that the model can be less sensitive to the number of clusters if $K$ is larger than the true clusters.

To determine robustness to model misspecification, we also simulate data from the spatial Cholesky process described as follows: The DT matrix for subject $i$ at voxel $v$ is determined by six independent spatial Gaussian processes $U_{ivk}$ ($k\in\{1,2,...,6\}$). These spatial Gaussian processes are arranged in the lower triangular matrix $\bm{L}_{iv}$ with $\bm{L}_{iv}=\begin{bmatrix}e^{U_{iv1}}&0&0\\ U_{iv4}&e^{U_{iv2}}&0\\ U_{iv5}&U_{iv6}&e^{U_{iv3}}\end{bmatrix}$. The responses $\bm{A}_{iv}$ are then constructed as $\bm{A}_{iv}=\bm{L}_{iv}\bm{L}_{iv}^{T}$, thereby introducing spatial dependence and guaranteeing positive definite of responses. We again use a $40\times 40$ grid with spacing 1 between adjacent grid points as an image. There are 10 subjects in the control group and 10 subjects in the treatment group. The six spatial Gaussian processes are simulated with variance $\tau^{2}=0.1$ and exponential correlation function with range parameter $\rho=2$. The mean of the six Gaussian processes are all [math] except for treatment subjects’ $10\times 10$ region in the center of the image where $U_{ivk}$ has mean $0.5$ for $k\leq 3$ and $0.25$ for $k>3$. This simulates the brain with a small region of difference between the two groups. We compare our model with $K=10,50,100$ to the Random Ellipsoid Model. The results averaged over 50 simulations are summarized in Table 1(b). The results demonstrate that our model maintains good performance, indicating the Potts Model is robust to this form of misspecification. In addition, under the spatial dependence assumption, the spatial models produce an overall better performance than the non-spatial model.

A problematic issue to the use of Bayesian methods in neuroimaging data is their heavy computational burden. In both simulations, the Potts Model has a computational speed within a few hours. The Random Ellipsoid Model avoids the expensive MCMC, but since the performance of the Random Ellipsoid Model is too conservative, the Potts Model is a reasonable trade-off.

5 Real Data Application

We apply this model to the data set of cocaine users (Ma et al., 2017) described in Section 2. The data are provided by the Institute for Drug and Alcohol Studies of Virginia Commonwealth University (VCU). The study recruited 11 cocaine users and 11 controls to test for microstructural changes of the brain through DTI. In the data analysis, we focus on the corpus callosum containing 15,273 voxels because this region plays important roles such as transferring motor, sensory, and cognitive information between the brain hemispheres (Ma et al., 2009). Conventionally, studies on cocaine use focus on this region (i.e., Ma et al., 2009, 2017; Lane et al., 2010).

Before model fitting, we examine the fit of the proposed model to the cocaine users data via empirical estimates of variograms. We denote $\hat{\mathcal{V}}_{ij}(d)=\frac{1}{N_{d}}\sum_{|u-v|=d}||\bm{A}_{iu}-\bm{A}_{jv}||_{F}^{2}$ as the empirical variogram value of subjects $i$ and $j$ at distance $d$, where $N_{d}$ is the number of pairs with $|u-v|=d$. We plot these empirical variograms of the motivating data in Figure 7. The DTs have a strong within-subject spatial dependence (Figure 9(a)). The empirical within-group variogram (Figure 9(b)) also increases with distance, indicating inter-subject dependence within a group, however, the empirical variogram is almost flat in the between-group variogram (Figure 9(c)), which suggests that the subjects are independent if they are in different groups. Since these empirical variograms perfectly match the theoretical variograms in Figure 6, the the hierarchical Potts model assumptions about spatial dependence are reasonable.

We sample $11,000$ MCMC samples with $3,000$ discarded as burn-in and it typically takes $5$ hours using a CPU with 3.4 GHz Intel Core i5. To study the sensitivity to $K$, we fit the model with $K$ as $100$, $200$, $300$, $400$ and $500$ and use the Rand index (Rand, 1971) for measuring the similarities of regions of differences detection with different $K$. The Rand index measures clustering similarity: if the two clusterings are almost identical, the index is close 1; otherwise, the index is close 0. In Table 2, the Rand indices of any two $K$ are close to 1, hence the selection is not sensitive to $K$. For a concise illustration, we use the result of $K=100$ in the rest of this section.

We first use the R package brainR (Muschelli et al., 2014) for 3-D visualizing the regions of differences. To investigate if the performance is improved by introducing spatial dependence, we also compare it to the Random Ellipsoid Model. We give the confidence level $0.9$ for the Random Ellipsoid Model. The selected regions of differences are displayed in Figure 8. The Rand index of the two clusterings is $0.86$ and so the results of the two analyses are similar but the Potts Model finds more spatial contiguous regions. Our study shows that a region of difference is detected in the splenium, which is consistent with previous clinical studies on cocaine use (see Lane et al., 2010). The splenium is a component located on the posterior end of the corpus callosum with an essential role on cognition. Since many studies revealed that the disease status at proximally-located/neighboring voxels can be similar (see Wu et al., 2013; Xue et al., 2018), our method may be more clinically meaningful than other studies because the detected regions are spatially contiguous. Furthermore, our test is able to find more regions of differences compared to alternative clinical studies on investigating the effect of cocaine use. For example, Ma et al. (2017) did not find such regions of differences by using the same data set.

The MCMC trace plots and posterior densities of the Potts spatial dependence parameters $\bm{\theta}$ are in Figure 9. The concentration parameter $\xi$ has $95\%$ credible region $[0.884,0.888]$, indicating there are a few active clusters. The group-clustering parameter $\alpha$ and spatial parameter $\beta$ control the within and inter subject spatial dependence and have $95\%$ credible regions $[0.323,0.327]$ and $[18.698,18.703]$, respectively. The dependence information revealed by the two credible regions is identical to the information obtained from the empirical variograms (Figure 7). Thus similar to the usage of the classic variogram, the generalized empirical variograms may also be a tool for obtaining the plug-in values of hyperparameters as an alternative (see Reich and Shaby, 2018).

6 Discussion

Although the spatial statistics literature on models and tools for matrix-variate data is sparse, the usage of positive definite matrix-variate is broad, which includes multiple-input and multiple-output (MIMO) systems (Smith and Garth, 2007) and computer vision (Cherian et al., 2016). Our major contribution is to present a spatial statistics formulation to model matrix-variate data via a Bayesian semiparametric mixture model. Our formulation retains the original data structure and accounts for spatial dependence by a computationally elegant model. In simulation studies, our model produces significantly improved performance compared to the non-spatial alternatives. The application to the DTI data set of cocaine users demonstrates the novelty of this model for detecting clinically meaningful regions of differences.

The current work primarily focuses on finding between-region differences in the brain at a single time-point (baseline). Temporally dependent matrix-variate data are also studied (Smith and Garth, 2007), and corresponding spatiotemporal extensions of our model are possible, though non-trivial. Extensions to incorporate covariates (i.e., socio-demographics, such as age, gender, etc.) may be possible via including a regression term in the full conditional distribution of the cluster labels.

Acknowledgements

We thank Institute for Drug and Alcohol Studies of Virginia Commonwealth University (VCU) for providing the cocaine users data set (Ma et al., 2017) and data manipulation instructions.

SUPPLEMENTAL MATERIALS

Appendix A: Density Functions

The PDFs of Parameterized Wishart and Inverse Wishart Distribution are given as

The PDF of Wishart Distribution $\bm{X}\sim\mathcal{W}_{p}(\bm{V},n)$:

[TABLE]

The PDF of Inverse Wishart Distribution $\bm{X}\sim\mathcal{IW}_{p}({\mathbf{\Psi}},\nu)$:

[TABLE]

Appendix B: Variograms

In this section, we give the details of derivations of variogram. We first give the variogram in the Single-Subject Model below:

[TABLE]

Obviously, $\mathbb{E}||\bm{A}_{iu}-\bm{A}_{iv}||_{F}^{2}$ can be derived in the same way. Next, we give the the explicit expression of $\gamma(m,\nu,\bm{\Sigma})$. We first have

[TABLE]

where $\sigma_{i}(.)$ returns the $i$-th eigenvalue of the input function.

The term $Tr(\mathbb{E}\bm{A}_{u}\bm{A}_{u})$ is complex. Gupta and Nagar (1999)[Section 3.3.6] provides trace moments of Wishart and inverse Wishart distribution.

Theorem 1**.**

Let $\bm{S}\sim\mathcal{W}_{p}(\bm{\Sigma},m)$ or $\bm{W}\sim\mathcal{IW}_{p}(\bm{M},m)$. Also, $\bm{S}^{-1}=\bm{W}$ and $\bm{M}=\bm{\Sigma}^{-1}$. $p$ is the matrix dimension. Then we have

1) $\mathbb{E}\bm{W}\bm{W}=(c_{1}+c_{2})\bm{M}\bm{M}(m-p-1)^{2}+c_{2}Tr(\bm{M})\bm{M}(m-p-1)^{2}$

2) $\mathbb{E}\bm{S}\bm{S}=\frac{m+1}{m}\bm{\Sigma}\bm{\Sigma}+\frac{1}{m}Tr(\bm{\Sigma})\bm{\Sigma}$,

3) $\mathbb{E}Tr(\bm{S})\bm{S}=2\frac{1}{m}\bm{\Sigma}\bm{\Sigma}+Tr(\bm{\Sigma})\bm{\Sigma}$

where $c_{1}=(m-p-2)c_{2}$ and $c_{2}=\frac{1}{(m-p)(m-p-1)(m-p-3)}$

Proof.

See Gupta and Nagar (1999)[Section 3.3.6] ∎

Then we first can obtain

[TABLE]

Next, we have Then we can obtain

[TABLE]

In summary, $\gamma(m,\nu,\bm{\Sigma})=2(\sigma-\lambda)$.

Appendix C: The Joint Probability Mass Density (PMF) of $\{\bm{g}_{1},\bm{g}_{2},...,\bm{g}_{N}\}\cup\{\bm{h}_{0},\bm{h}_{1}\}$

In this section, we show the expression of the PMF and validate that the PMF is valid:

[TABLE]

where $u\sim v$ means $u$ and $v$ are connected. It is obvious that the probability is positive and satisfies the pairwise Markov property stated in the Hammersley and Clifford Theorem. The normalizing constant $Z(\alpha,\beta,\xi)=\sum_{\bm{g},\bm{h}}\exp(U(\bm{g},\bm{h},\bm{\theta}))$ is intractable. Since the summation is over finite and discrete indices, we have that $0<Z(\alpha,\beta,\xi)<\infty$, revealing that $P(\bm{h},\bm{g})$ is proper.

Appendix D: The Statistical Role of $h_{xv}$

Step 1: Marginalizing $g_{iv}$

[TABLE]

Step 2: Marginalizing $\{g_{u}:u\in N_{v}\}$

[TABLE]

Appendix E: Codes

The codes and example scripts are available at https://github.com/ZhouLanNCSU/Potts_DTI.