Spatial von-Mises Fisher Regression for Directional Data

A random $p$-dimensional directional data point, $\bm{E}_{iv}$, lies on $\mathcal{S}^{p-1}$, as described in Section 1. The von Mises-Fisher (vMF) distribution (Mardia and Jupp, 2009) is a commonly used probability distribution for handling directional data, and spherical regression models are built based on the vMF distribution. However, building spherical regression models is not straightforward since the inference might depend on an arbitrary coordinate system in which the responses are defined. Following the approaches in Scealy and Wood (2019); Paine et al. (2020), We assume that $\bm{Q}^{T}\bm{E}_{iv}$ follows a von Mises-Fisher (vMF) distribution (Mardia and Jupp, 2009), denoted as $\bm{Q}^{T}\bm{E}_{iv}\sim\text{vMF}\big{(}\bm{\mu}_{iv},\kappa\big{)}$ with the probability density function expressed as

$\displaystyle f\Big{[}\bm{Q}^{T}\bm{E}_{iv}|\bm{\mu}_{iv},\kappa\Big{]}=C_{p}(\kappa)\exp\Big{(}\kappa\bm{\mu}_{iv}^{T}\bm{Q}^{T}\bm{E}_{iv}\Big{)}=C_{p}(\kappa)\exp\Big{(}\kappa\mathcal{M}_{iv}^{T}\bm{E}_{iv}\Big{)},$

(1)

In the expression, $C_{p}(\kappa)$ is the normalizing constant such as $C_{p}(\kappa)={\frac{\kappa^{p/2-1}}{(2\pi)^{p/2}I_{p/2-1}(\kappa)}}$, where $I_{a}$ stands for the modified Bessel function of the first kind at order $a$. In this paper, we primarily focus on the case with $p=3$, due to the nature of the DTI data, and in that case, $C_{3}(\kappa)=\frac{\kappa}{2\pi(e^{\kappa}-e^{-\kappa})}$ and the density simplifies to $\frac{\kappa\exp(\kappa(\mu^{T}\mathbf{x}-1))}{2\pi(1-\exp(-2\kappa))}$ for a point $\mathbf{x}\in\mathcal{S}^{2}$. There are some pictorial illustrations available in Straub (2017). The term $\bm{\mu}_{iv}$ is a directional vector, i.e., $\|\bm{\mu}_{iv}\|=1$, where $\|\cdot\|$ represents the $\ell_{2}$ norm. $\kappa\geq 0$ is a concentration parameter. The term $\bm{Q}$ is set as an unknown parameter in our model and modeled as a special orthogonal matrix.

The density is maximized at $\bm{E}_{iv}=\bm{Q}\bm{\mu}_{iv}$, the mode direction parameter. Thus, we further introduce the notation

$$\mathcal{M}_{iv}=\bm{Q}\bm{\mu}_{iv},$$

(2)

and it is the mode direction for the original response $\bm{E}_{iv}$, and is named as direct mode direction to differentiate it from $\bm{\mu}_{iv}$ which is called rotated mode direction. In the next subsection, we propose a link function to relate predictors with $\bm{\mu}_{iv}$. Since the link function is known and pre-specified, the orthogonal matrix $\bm{Q}$ here provides additional flexibility by identifying an optimal coordinate system with respect to the link.

In accordance with the generalized linear regression (GLM) framework, we propose a relationship between a linear combination of the predictors and the response variables through a link function. It is also reasonable to allow the linear combination of the predictors to be unbounded, mirroring other GLM frameworks. Therefore, our proposed link function is a composition of two functions, each handling a different level of transformations. The following section provides detailed insight into the construction of the link function.

In this section, we focus on our link function for 3-dimensional spherical coordinates, i.e., when $\bm{E}_{iv}$ lies on $\mathcal{S}^{p-1}$ and $p=3$. Let $(x,y,z)$ represent the Cartesian coordinates of a point on the sphere such that $x^{2}+y^{2}+z^{2}=1$. This point can be represented in spherical coordinates using two parameters, as depicted in Figure 3(a). The azimuth angle $\theta\in(-\pi,\pi]$ represents the counterclockwise angle in the x-y plane from the positive x-axis, and the elevation angle $\phi\in(-\frac{\pi}{2},\frac{\pi}{2}]$ is the angle from the x-y plane. The one-to-one transformation from Cartesian to spherical coordinates is as follows:

$$u:(x,y,z)\rightarrow\Big{(}\text{atan2}(y,x),\text{atan2}(z,\sqrt{1-z^{2}})\Big{)}=(\theta,\phi),$$

where $\text{atan2}(y,x)$ is defined by the limit expression $\lim_{z\rightarrow x^{+}}(\frac{y}{z})+\frac{\pi}{2}\text{sgn}(y)\text{sgn}(x)(\text{sgn}(x)-1)$. While the spherical coordinates $(\theta,\phi)$ are located in Euclidean space, they are bounded. Hence, we first scale the spherical coordinates $(\theta,\phi)$ to fall within the interval $(0,1]$, then apply the logit function. The transformation is:

$$g:(\theta,\phi)\rightarrow\Big{(}\text{logit}\big{(}\frac{\theta+\pi}{2\pi}\big{)},\text{logit}\big{(}\frac{\phi+\pi/2}{\pi}\big{)}\Big{)}=(\widetilde{\theta},\widetilde{\phi})\in\mathbb{R}\times\mathbb{R}.$$

Here, $\text{logit}(x)=\log\left(\frac{x}{1-x}\right)$. Although the logit transformation ill-behaves on the boundary (i.e., at $x=1$), the regression model built on this transformation is still adequately flexible, as the set of such degenerate cases has measure zero. This transformation serves to map the bounded spherical coordinates to a new set of parameters $(\widetilde{\theta},\widetilde{\phi})$, which are unbounded. For clarity, the overall mapping from $(x,y,z)$ to $(\widetilde{\theta},\widetilde{\phi})$ is denoted by $\ell(\cdot)$, which is a composite of $g(\cdot)$ and $u(\cdot)$, i.e., $\ell(\cdot):=u(\cdot)\circ g(\cdot)$. Thus, $\ell(x,y,z)=[\widetilde{\theta},\widetilde{\phi}]$ represents this innovative function that projects the directional data from $\mathcal{S}^{2}$ to $\mathbb{R}\times\mathbb{R}$.

A point on a $p$ dimensional hypersphere, $\mathcal{S}^{p-1}$, can be expressed using $p-1$ angles. The first $p-2$ angles (denoted as ${\bm{\phi}}^{(p-2)}=[{{\phi}}_{1}^{(p-2)},\ldots,{{\phi}}_{p-2}^{(p-2)}]$) are within the interval $(-\pi/2,\pi/2]$, while the last angle ${\theta}$ is within $(-\pi,\pi]$. The mapping function, when extended to the general dimension, is formulated as:

		$\displaystyle g:(\theta,\bm{\phi}^{(p-2)})\rightarrow\Big{(}\text{logit}\big{(}\frac{\theta+\pi}{2\pi}\big{)},\big{[}\text{logit}\big{(}\frac{{\phi}_{1}^{(p-2)}+\pi/2}{\pi}\big{)},\ldots,\text{logit}\big{(}\frac{{\phi}_{p-1}^{(p-2)}+\pi/2}{\pi}\big{)}\big{]}\Big{)}$
		$\displaystyle=(\widetilde{\theta},\big{[}\widetilde{{\phi}}_{1}^{(p-2)},\ldots,\widetilde{{\phi}}_{p-2}^{(p-2)}\big{]})\in\mathbb{R}\times\mathbb{R}^{p-2}.$

In our proposed generalized regression model, the link function $\ell(\cdot)$ plays a key role. The term $\ell(\bm{\mu}_{iv})=(\widetilde{\theta},\big{[}\widetilde{{\phi}}_{1}^{(p-2)},\ldots,\widetilde{{\phi}}_{p-2}^{(p-2)}\big{]})$ is considered as the prediction terms and is formulated as a function of the covariate vector $\bm{X}_{i}$. Specifically, we have:

$$\mathbb{E}(\widetilde{\theta},\big{[}\widetilde{{\phi}}_{1}^{(p-2)},\ldots,\widetilde{{\phi}}_{p-2}^{(p-2)}\big{]})=(\bm{X}_{i}\bm{\alpha}_{v},[\bm{X}_{i}\bm{\beta}^{(p-2)}_{v,1},\ldots,\bm{X}_{i}\bm{\beta}^{(p-2)}_{v,p-2}]).$$

Lemma 1 illustrates that the link function $\ell(\cdot)$ is a bijective function. This provides us with a direct mapping between the coordinates $(\widetilde{\theta},\big{[}\widetilde{{\phi}}_{1}^{(p-2)},\ldots,\widetilde{{\phi}}_{p-2}^{(p-2)}\big{]})$ and the mode direction $\bm{\mu}_{iv}$. This mapping is crucial in facilitating the interpretation of the regression coefficients.

For the remainder of this article, in view of our DTI application and to simplify the exposition, we set $p=3$, i.e.,

$$\mathbb{E}[\widetilde{\theta}_{iv},\widetilde{\phi}_{iv}]^{T}=[\bm{X}_{i}\bm{\alpha}_{v},\bm{X}_{i}\bm{\beta}_{v}]^{T},$$

and drop the superscript $(p-2)$. Despite this simplification, it should be noted that it is straightforward to generalize all the steps we take here to accommodate any general $p$.

Spatial dependence is an essential component in many applications. In this section, we further illustrate the model specification which accounts for spatial dependence. By preserving $\mathbb{E}[\widetilde{\theta}_{iv},\widetilde{\phi}_{iv}]^{T}=[\bm{X}_{i}\bm{\alpha}_{v},\bm{X}_{i}\bm{\beta}_{v}]^{T}$, we let $[\widetilde{\theta}_{iv},\widetilde{\phi}_{iv}]^{T}=[\bm{X}_{i}\bm{\alpha}_{v}+\epsilon_{iv},\bm{X}_{i}\bm{\beta}_{v}+\xi_{iv}]^{T}$, where $\bm{\epsilon}_{i}=(\epsilon_{1i},\ldots,\epsilon_{iv},\ldots,\epsilon_{iV})$ and $\bm{\xi}_{i}=(\xi_{1i},\ldots,\xi_{iv},\ldots,\xi_{iV})$ are distributed as the mean-zero Gaussian distribution with covariance matrices $\bm{\Sigma}^{(\epsilon)}$ and $\bm{\Sigma}^{(\xi)}$, respectively.

The focus of our applications is principal directional data where the spatial dependence is along a given fiber streamline. Suppose we have $K$ fiber streamlines, $\mathbb{Z}_{k}$ is a set of consecutive indices within the $k$-th fiber streamline. We let the directional statistics are auto-correlated along a given streamline. For completeness, we first state the definition and notations related to the Gaussian autoregressive modeling in Definition 1.

We thus have $\bm{\epsilon}_{i}=(\epsilon_{1i},\ldots,\epsilon_{iv},\ldots,\epsilon_{iV})$ and $\bm{\xi}_{i}=(\xi_{1i},\ldots,\xi_{iv},\ldots,\xi_{iV})$ follow a multivariate normal distribution with an autoregressive covariance matrix, denoted as

	$\displaystyle\bm{\epsilon}_{ik}=\{\epsilon_{iv}:v\in\mathbb{Z}_{k}\}\sim\text{AR}_{P}\Big{(}[\Phi_{\epsilon 1},\ldots,\Phi_{\epsilon P}],\tau_{\epsilon}^{2}\Big{)}$		(3)
	$\displaystyle\bm{\xi}_{ik}=\{\xi_{iv}:v\in\mathbb{Z}_{k}\}\sim\text{AR}_{P}\Big{(}[\Phi_{\xi 1},\ldots,\Phi_{\xi P}],\tau_{\xi}^{2}\Big{)},$

where $(\Phi_{\epsilon,1},\ldots,\Phi_{\epsilon,P},\tau_{\epsilon}^{2})$ and $(\Phi_{\xi,1},\ldots,\Phi_{\xi,P},\tau_{\xi}^{2})$ are the corresponding parameters. Note that our characterization of the AR$(P)$ process is based on the partial correlations but not the autoregressive coefficients directly.

We now put priors on the model parameters to proceed with our Bayesian inference. The matrix $\bm{Q}$ is first reparametrized using the Cayley transformation as $\bm{Q}=(\bm{I}-\bm{A})(\bm{I}+\bm{A})^{-1}$, where $\bm{A}$ is a skew-symmetric matrix. This parametrization is useful for efficient posterior computation. Setting $\bm{A}=\begin{bmatrix}0&a_{1}&a_{2}\\ -a_{1}&0&a_{3}\\ -a_{2}&-a_{3}&0\end{bmatrix}$, the priors for the three parameters $a_{1},a_{3},a_{3}$ are set as mean-zero normal with variance $100$. For the regression coefficients, we also assign AR-$P$ process priors, denoted as

	$\displaystyle\bm{\alpha}_{k}(c)=\{\alpha_{v}(c):v\in\mathbb{Z}_{k}\}\sim\text{AR}_{P}\Big{(}[\Phi_{\alpha,1},\ldots,\Phi_{\alpha,P}],\sigma_{\alpha}^{2}\Big{)}$		(4)
	$\displaystyle\bm{\beta}_{k}(c)=\{\beta_{v}(c):v\in\mathbb{Z}_{k}\}\sim\text{AR}_{P}\Big{(}[\Phi_{\beta,1},\ldots,\Phi_{\beta,P}],\sigma_{\beta}^{2}\Big{)},$

where $\bm{\alpha}_{v}=\Big{(}{\alpha}_{v}(0),\ldots,{\alpha}_{v}(c),\ldots,{\alpha}_{v}(C)\Big{)}$ and $\bm{\beta}_{v}=\Big{(}{\beta}_{v}(0),\ldots,{\beta}_{v}(c),\ldots,{\beta}_{v}(C)\Big{)}$. The terms ${\alpha}_{v}(0)$ and ${\beta}_{v}(0)$ are the coefficients for the intercepts. The terms ${\alpha}_{v}(c)$ and ${\beta}_{v}(c)$ are the coefficients for the covariate $c$. The variances are assigned with weakly informative inverse-gamma priors with shape and rate parameters set to 0.1, denoted as $\tau_{\epsilon}^{-2},\tau_{\xi}^{-2},\sigma_{\alpha}^{-2},\sigma_{\beta}^{-2}\sim\mathcal{GA}(0.1,0.1)$. We put a diffuse prior for the concentration parameter, $\kappa$, denoted as $\kappa\propto 1$. The concentration parameter $\kappa$ can be viewed as a nugget effect. An autoregressive process with the partial autocorrelation parameters restricted to $(-1,1)$ can ensure stationarity (Barnett et al., 1996). Therefore, we assign uniform $(-1,1)$ prior on the partial autocorrelation parameters. Specifically, $\rho_{\alpha,p},\rho_{\beta,p},\rho_{\epsilon,p},\rho_{\xi,p}\sim\mathcal{U}(-1,1)$ for $p\in\{1,2,\ldots,P\}$, where $\rho_{\alpha,p},\rho_{\beta,p},\rho_{\epsilon,p},\rho_{\xi,p}$ are the corresponding partial autocorrelation parameters associated with the autoregressive process priors for $\bm{\alpha}_{v}(c)$’s, $\bm{\beta}_{v}(c)$’s, and the autoregressive random effects $\bm{\epsilon}_{ik}$’s, and $\bm{\xi}_{ik}$’s, respectively. Applying the Durbin-Levinson recursion on the partial autocorrelations ($\rho_{\alpha,p},\rho_{\beta,p},\rho_{\epsilon,p},\rho_{\xi,p}$), we get back the autoregressive coefficients ($\Phi_{\alpha,p},\Phi_{\beta,p},\Phi_{\epsilon,p},\Phi_{\xi,p}$) using the R function inla.ar.pacf2phi() from the R package inla (Lindgren and Rue, 2015). The reverse transformation is applied using inla.ar.phi2pacf() from the same package. We utilize a Markov chain Monte Carlo (MCMC) algorithm, which includes Metropolis-Hastings-within-Gibbs steps, to generate the posterior samples needed for our Bayesian inference. The codes and usages of our codes are stated in Section B.1. of the supplementary materials.

For a random unit vector $\mathcal{M}$, the angular expectation is defined as

$$\mathbb{A}\mathcal{M}=\arg\min_{\bm{m}\in\mathcal{S}^{p-1}}\int\delta(\bm{m},\mathcal{M})\pi(\mathcal{M})d\mathcal{M},$$

where $\mathbb{A}$ is introduced as an operator returning the angular expectation, $p$ is the dimension of $\mathcal{M}$, and $\pi(\mathcal{M})$ stands for the distribution of $\mathcal{M}$. The angular expectation $\mathbb{A}$ is the main tool for our proposed Bayesian angular inference scheme. For a new subject with the covariate vector $\bm{X}_{\text{new}}$, the posterior predictive distribution of the direct mode direction $\bm{\mu}_{\text{new},v}$ for the $v$-th index (Equation 2), is $\Big{[}\mathcal{M}_{\text{new},v}|\bm{X}_{\text{new}};\text{data}\Big{]}$. For clarity, we use subscript “new” to replace the original subscript “$i$” in the model for referring to a new subject. For summarizing above posterior of $\mathcal{M}_{\text{new},v}$, we compute the posterior angular expectation using the MCMC samples $\{\mathcal{M}_{\text{new},v}^{(t)}:t=1:T\}$ as

$\displaystyle\mathbb{A}\Big{[}\mathcal{M}_{\text{new},v}|\bm{X}_{\text{new}};\text{data}\Big{]}\approx\arg\min_{\bm{m}\in\mathcal{S}^{p-1}}\frac{1}{T}\sum_{t=1}^{T}\delta(\bm{m},\mathcal{M}_{\text{new},v}^{(t)})=\frac{\sum_{t=1}^{T}\mathcal{M}_{\text{new},v}^{(t)}/T}{\|\sum_{t=1}^{T}\mathcal{M}_{\text{new},v}^{(t)}/T\|}.$

Relating to the real applications, the above term profiles the expected fiber direction after the model-based adjustment for the covariate vector $\bm{X}_{\text{new}}$.

A transparent inference framework illustrating the covariate effects is essential for any scientific study. In our proposed model, the coefficients $\bm{\alpha}_{v}$ and $\bm{\beta}_{v}$ do not offer a straightforward illustration of the covariate effects on the mode directions. This is also a limitation of many regression models for complicated responses. To overcome this, we incorporate the tangent-normal decomposition (Mardia and Jupp, 2009, Page 169) to summarize the covariate effects on directional data as follows. The tangent-normal decomposition is an integrated approach, but can be decomposed into the following three steps for a concise illustration.

Our proposed procedure starts by defining the following two covariate vectors: a) $\bm{X}_{\text{T}}=[1,{X}_{\text{T}1},\ldots,{X}_{\text{T}c},\ldots,{X}_{\text{T}C}]$: Covariate vector for the typical condition; b) ${\bm{X}}_{\text{S}}=[1,{X}_{\text{S}1},\ldots,{X}_{\text{S}c},\ldots,{X}_{\text{S}C}]$: Covariate vector for the specific condition. Thus, we replace the original subscript “$i$” with “T” and “S” in our notation, to represent a subject in the typical and specific condition, respectively.

Next, we obtain the corresponding predictive posterior distributions of the mode directions with the given typical or specific covariate vector. They are referred as to the typical and specific mode direction, respectively, expressed as follows: a) Typical Mode Direction: $\overline{\mathcal{M}}_{v}=\Big{[}\mathcal{M}_{\text{T}v}|{\bm{X}}_{\text{T}};\text{data}\Big{]}$; b) Specific Mode Direction: $\widetilde{\mathcal{M}}_{v}=\Big{[}\mathcal{M}_{\text{S}v}|{\bm{X}}_{\text{H}};\text{data}\Big{]}$. They can also be approximately learned from the MCMC output through the $\mathbb{A}$ operator defined above.

To quantify the covariate effect, we employ the tangent-normal decomposition on the specific mode direction $\widetilde{\mathcal{M}}_{v}$ as in Figure 3(b) where the typical mode direction $\overline{\mathcal{M}}_{v}$ is the tangent direction and the term $\bm{R}_{v}$ is the normal direction. Applying this to the principal diffusion directions, the term $\bm{R}_{v}$ can be interpreted as the deviation direction that is caused by the difference from the specific condition to the typical condition. The scalar $m_{v}\in(0,1)$, which controls the magnitude of this effect, can be used to quantify the importance of this covariate: the larger value of $m_{v}$ indicates the stronger deviation between the specific condition to the typical condition. Computationally, the tangent-normal decomposition can be applied to each MCMC sample $t$. This provides us with the posterior samples of $\bm{R}_{v}$ and $m_{v}$, denoted as $\{\bm{R}_{v}^{(t)}:t=1:T\}$ and $\{m_{v}^{(t)}:t=1:T\}$. We use those samples to obtain posterior estimates of $\mathbb{A}[\bm{R}_{v}|\text{data}]$ to spatially profile the deviation directions caused by the difference between typical condition and specific condition. Similarly, the posterior mean estimate of $m_{v}$ can also be computed and used to quantify the magnitude. We use the following notation to express the tangent-normal decomposition:

$\displaystyle\mathcal{D}_{TN}(\overline{\mathcal{M}}_{v}||\widetilde{\mathcal{M}}_{v})=\{\bm{R}_{v},m_{v}\}$

(5)

Our above illustration provides the general description of tangent-normal decomposition. As long as it is feasible to specify the typical and specific conditions, the may use $\{\bm{R}_{v},m_{v}\}$ to infer the covariate effects caused by the difference between the typical and specific conditions. The below are two important remarks which are important for the general use of our approach.

While we illustrate the use of our approach as a generic use, interpreting one covariate effect at a time is usually a common manner. Therefore, we further provide two examples below to illustrate how to interpret one covariate effect at a time. In summary, the key strategy is to have typical and specific conditions differ only in the covariate $c^{\prime}$ of our primary interest.

In this section, we carry out numerical analyses on both the ADNI data and synthetic principal diffusion direction data. The purpose of these studies is to demonstrate the superior performance of our proposed vMF regression relative to other methods. Here, we consider three competing methods based on multivariate Gaussian distribution: 1) Gaussian Regression 1: a multivariate regression model for the principal diffusion directions $\bm{E}_{iv}$, 2) Gaussian Regression 2: a multivariate regression model for transformed means $\ell\Big{(}\bm{E}_{iv}\Big{)}$, 3) Non-spatial vMF regression: $\bm{Q}^{T}\bm{E}_{iv}\sim\text{vMF}\big{(}\mathcal{M}_{iv},\kappa\big{)}$ and $\ell\Big{(}\bm{\mu}_{iv}\Big{)}=[\widetilde{\theta}_{iv},\widetilde{\phi}_{iv}]^{T}=[\bm{X}_{i}\bm{\alpha}_{gv},\bm{X}_{i}\bm{\beta}_{gv}]^{T}$. Specific details are described as follows. The Gaussian Regression 1 simply treats the principal diffusion directions $\bm{E}_{iv}$ as normally distributed random variables, such that $\bm{E}_{iv}\sim\mathcal{N}(\mathcal{M}_{iv},\bm{\Sigma}_{1})$, where $\mathcal{M}_{iv}^{(p)}=\bm{X}_{i}\bm{u}_{v}^{(p)}$, $\mathcal{M}_{iv}^{(p)}$ is the $p$-th element of $\mathcal{M}_{iv}$ for $p\in\{1,2,3\}$, i.e., $\mathcal{M}_{iv}=(\mathcal{M}_{iv}^{(1)},\mathcal{M}_{iv}^{(2)},\mathcal{M}_{iv}^{(3)})$. Alternatively, in case of Gaussian Regression 2:, the terms $\ell\Big{(}\bm{E}_{iv}\Big{)}=[A_{iv},B_{iv}]^{T}$ are assumed to follow a normal distribution, denoted as $\ell\Big{(}\bm{E}_{iv}\Big{)}=[A_{iv},B_{iv}]^{T}\sim\mathcal{N}([\tilde{\theta}_{iv},\tilde{\phi}_{iv}]^{T},\bm{\Sigma}_{2})$, where $\tilde{\theta}_{iv}=\bm{X}_{i}\bm{a}_{v}$ and $\tilde{\phi}_{iv}=\bm{X}_{i}\bm{b}_{v}$. In case of Non-spatial vMF Regression:, we simply remove the spatially dependent randomness from the proposed model.

We take the Bayesian route for inference for all these methods and put weakly informative priors on the unknown parameters. Specifically, we give normal priors with mean $\bm{0}$ and covariance matrix $1000\times\bm{I}$ on the coefficients, denoted as $\bm{u}_{v}^{(1)},\bm{u}_{v}^{(2)},\bm{u}_{v}^{(3)},\bm{a}_{v},\bm{b}_{v}\sim\mathcal{N}(\bm{0},1000\times\bm{I})$. We further put inverse Wishart prior with degrees of freedom $5$ and scale matrix $\bm{I}$ on the covariance matrices, denoted as $\bm{\Sigma}_{1}^{-1},\bm{\Sigma}_{2}^{-1}\sim\mathcal{W}(5,\bm{I})$. Our inference is based on the posterior samples, generated by applying the MCMC algorithm, where we collect $3,000$ post-burn samples after discarding the first $2,000$. The MCMC convergence is monitored by Heidelberger and Welch’s convergence diagnostic (Heidelberger and Welch, 1981).

We generate synthetic principal diffusion direction data using our proposed vMF regression model. We design a synthetic fiber tract that includes $K=12$ fibers and a total of $1,256$ voxels. We also have covariate vectors for 20 subjects for training and 20 subjects for validation. For each data replication, we fit our model to the training data to obtain the parameter estimates and these parameter estimates are used to predict the response in the validation data. The discrepancies between the predict the responses and true responses are used to measure model accuracy.

The tractography atlas, values of covariate vectors and voxel-wise coefficients utilized are stored in the R file synthetic.Rdata, as detailed in Section B.2 of the supplementary materials. For other parameters, we set $P=5$ and generate partial correlation coefficients from $\mathcal{U}(-1,1)$. The variances are specified as $\tau_{\epsilon}^{2}=\tau_{\xi}^{2}=\sigma_{\alpha}^{2}=\sigma_{\beta}^{2}=1$, and orthogonal matrix is specified as $\bm{Q}=\bm{I}$. We further give $\kappa\in\{20,30,40\}$, respectively. A total of 50 replicated data sets for each setting.

Prediction errors are quantified in terms of a) the separation angle b) the root mean squared error. The separation angle is defined as $\sum_{i,v}\delta\Big{(}\widehat{\mathcal{M}}_{iv},\bm{E}_{iv}\Big{)}/N_{iv}$ for vMF regression and $\sum_{i,v}\delta\Big{(}\widehat{\mathcal{M}}_{iv},\bm{E}_{iv}\Big{)}/N_{iv}$ for Gaussian regression methods. The root mean squared error is defined as $\sum_{i,v}\|\widehat{\mathcal{M}}_{iv}-\hat{\bm{Q}}^{T}\bm{E}_{iv}\|/N_{iv}$ for vMF regression and $\sum_{i,v}\|\widehat{\mathcal{M}}_{iv}/\|\widehat{\mathcal{M}}_{iv}\|-\bm{E}_{iv}\|/N_{iv}$ for Gaussian case. We give the prediction errors spanning simulation replications in Figure 4. Our proposed method substantially outperforms Gaussian alternatives and the non-spatial vMF regression model based on these prediction errors. The spatial vMF regression with $P=5$ demonstrates comparatively better performance among $P\in\{2,\ldots,7\}$. This aligns with the fact that true data is generated by setting $P=5$. It suggests that the optimal lag can be chosen based on the predictive performances for our real data analysis.

In this section, we utilize our Bayesian angular inference method to analyze data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI). We illustrate the effects of covariates via our proposed tangent-normal decomposition. Our techniques shed fresh light on the neurological progression of Alzheimer’s disease. To maintain a specific emphasis in our study, we concentrate on examining the influence of the APOE-4 allele on brain anatomical structure. This analysis takes into account other covariates such as age, gender, and scores from the Mini-Mental State Examination (MMSE). We focus on fornix in this section.

To achieve this goal, the typical covariate vector is given as follows

$\displaystyle\bar{\bm{X}}_{T}^{(g)}=[1,\bar{X}_{age}^{(g)},\bar{X}_{T,gender}=0,\bar{X}_{MMSE}^{(g)},X_{T,APOE}=0]$

(8)

or

$\displaystyle\bar{\bm{X}}_{T}^{(g)}=[1,\bar{X}_{age}^{(g)},\bar{X}_{T,gender}=1,\bar{X}_{MMSE}^{(g)},X_{T,APOE}=0]$

(9)

where Equations (10) and (11) are the typical covariate vectors of male and female, respectively. We consider the wild type of APOE-4 (i.e., $X_{T,APOE}=0$) as typical. The terms $\bar{X}_{age}^{(g)}$ and $\bar{X}_{MMSE}^{(g)}$ are the averages of age and MMSE within group $g$. We decompose the gender specification into male and female, in order to provide the gender-specific inference. The corresponding specific covariate vector is given as

$\displaystyle\bar{\bm{X}}_{S}^{(g)}=[1,\bar{X}_{age}^{(g)},\bar{X}_{T,gender}=0,\bar{X}_{MMSE}^{(g)},X_{S,APOE}=1]$

(10)

or

$\displaystyle\bar{\bm{X}}_{S}^{(g)}=[1,\bar{X}_{age}^{(g)},\bar{X}_{s,gender}=1,\bar{X}_{MMSE}^{(g)},X_{S,APOE}=1],$

(11)

where the APOE-4 status is considered as mutated (i.e., $X_{S,APOE}=1$).

To incorporate group-level analysis into our study, we use the superscript $(g)$ to denote group-specific terms. In this section, we focused on the investigation of the effect of APOE-4. The posterior $[m_{v}^{(g)}|\text{data}]$ implies the strength of the covariate effect and we consider $m_{v}^{(g)}>0.65$ as large impact. Therefore, we calculate the posterior probabilities for $Pr(m_{v}^{(g)}-0.65|\text{data})$ As an exploratory analysis, we consider the voxels whose posterior z-value larger than the 9-th quantile of the posterior z-values over voxels and group as impacted voxels. The visualization of the ADNI group is in Figure 5. Compared to other clinical groups (see Figure A.2 in supplementary materials), the ADNI group’s fornix anatomical stricture is more impacted by the factor of APOE. Also, male has more laterality of the APOE effect, compared to women. Male’s left component of fornix is more impacted by the factor of APOE.