Statistical Inference of Auto-correlated Eigenvalues with Applications to Diffusion Tensor Imaging

Diffusion tensor imaging (DTI) is a prevalent neuroimaging tool in analyzing the anatomical structure of a tissue. The rich use of DTI in brain science triggers many meaningful interdisciplinary studies. Due to the nature of interdisciplinarity, providing a tutorial of DTI, which is intuitive to both statisticians and clinicians, is key in interdisciplinary researches of DTI. Based on Figure 1 which is modified from figures in (Zhou, 2010, Figure 1.6), we give a short tutorial of DTI as follows. DTI is a neuroimage measuring the white matter tracts within a tissue (e.g., human brain) voxel-by-voxel. Let the structure at the upper panel of Figure 1 be the anatomical structure measured at a voxel. The movement of water molecules in different environments is measured by giving their respective diffusion ellipsoids at the lower panel of Figure 1. The eigenvectors are the diffusion directions $[\bm{e}_{1}\ \bm{e}_{2}\ \bm{e}_{3}]$ and the squared roots of the eigenvalues $[l_{1}\ l_{2}\ l_{3}]$ are the corresponding semi-axis lengths (Zhou, 2010, Section 1.2.3). The corresponding diffusion tensor $D=[\bm{e}_{1}\ \bm{e}_{2}\ \bm{e}_{3}]\text{diag}(l_{1},l_{2},l_{3})[\bm{e}_{1}\ \bm{e}_{2}\ \bm{e}_{3}]^{T}$ describes the anatomical structure at the voxel. Combining the diffusion tensors, which are positive definite over all the voxels, we have an image whose voxel-wise variables are positive definite diffusion tensors, describing the anatomical structure of a tissue.

The use of DTI is extensive (e.g., fiber tracking, clinical studies, etc.). Our paper focuses on understanding how the subject’s covariates (e.g., age, gender, disease status) drive the brain’s anatomical structures in terms of DTI. Many brain scientists raised this scientific objective. The relevant researches are now divided into two paths. The statisticians are interested in developing methodologies that handle the statistical randomness of the whole diffusion tensor. There are many representative works (e.g., Schwartzman et al., 2008; Zhu et al., 2009; Yuan et al., 2012; Lee and Schwartzman, 2017; Lan et al., 2021). The clinicians are more interested in using convenient summary scalar to tell the changes of diffusion tensor quickly and easily. For example, the fractional anisotropy, denoted as ${\displaystyle{\text{FA}}={\sqrt{\frac{1}{2}}}{\frac{\sqrt{(l_{1}-l_{2})^{2}+(l_{2}-l_{3})^{2}+(l_{3}-l_{1})^{2}}}{\sqrt{l_{1}^{2}+l_{2}^{2}+l_{3}^{2}}}},}$ is a frequently used summary scalar. Via implementing basic statistical tools, i.e., t-test, simple linear regression with these quantities as data, the researchers can easily tell the changes of the tissue statistically (e.g., Ma et al., 2009; Lane et al., 2010).

In light of Figure 1, we propose that analyzing the DTI’s eigenvalues and eigenvectors separately can be another promising research direction. For statisticians, this path does not preclude any valuable available data information. For clinicians, the scientific explanations in terms of the statistical uncertainties of eigenvalues and/or eigenvectors are more transparent than the other ways, e.g., coefficients based on Cholesky decomposition (Zhu et al., 2009; Lan et al., 2021) or local polynomial regression (Yuan et al., 2012). In our paper, we focus on the statistical analysis of eigenvalues because the variations of eigenvalues are primarily interested by the clinicians given the current clinical reports (e.g., Ma et al., 2009; Lane et al., 2010). In this paper, the analysis of eigenvectors is left as future works.

Compared to other statistical models regarding the eigenvalues of DTI (e.g., Zhu et al., 2010, 2011; Jin et al., 2017), our paper may be a very first paper discussing the statistical inference of random eigenvalues. It is an important fact that the eigenvalues are stochastically ordered, and thus the probability density function of the eigenvalues must ensure the inequality of the eigenvalues stochastically. This inequality makes many classic models such as multivariate Gaussian processes not be applied to this problem straightforwardly. In this paper, we use the density function of a Wishart distributed matrix’s eigenvalues as the basic probability distribution, referred to as distribution of Wishart’s eigenvalues. The distribution of Wishart’s eigenvalues ensures the inequality of the eigenvalues stochastically.

Following many neuroimaging studies, the spatial dependence of the eigenvalues among voxels is encouraged to be incorporated. In the previous analysis (e.g., Zhu et al., 2010, 2011; Jin et al., 2017), the spatial dependence can be easily achieved via the spatial Gaussian process model. Many fiber tracking reports (e.g., Wong et al., 2016) assume that the diffusion tensors are smooth along with a fiber tract (see Figure 2) whose diffusion tensors have similar directions. In light of this assumption, we assume that the eigenvalues are auto-correlated. The model based on auto-correlated eigenvalues captures the spatial correlation and produces an analytic expression of the density function of the spatially correlated eigenvalues (James et al., 1964).

Combining all the above, we construct a hierarchical spatial model for the random auto-correlated eigenvalues. This model is both statistically and scientifically proper regarding the applications to diffusion tensor imaging. Given the hierarchical structure, we implement the Monte-Carlo Expectation-Maximization (EM) algorithm (Levine and Casella, 2001, Section 2.1) for the parameter estimation and implement parametric bootstrapping for uncertainties quantification.

Both simulation studies and real data applications are used to demonstrate our proposed method. We compare our proposal to a non-linear regression model treating the data as normally distributed. The simulation studies indicate that our model gives a more accurate parameter point and confidence interval estimation. The real data analysis based on the IXI data-set (https://brain-development.org/ixi-dataset/) provided by Center for the Developing Brain, Imperial College London, an open data-set made available under the Creative Commons CC BY-SA 3.0 license, shows more clinically meaningful results after treating the data as random eigenvalues.

In the rest of the paper, we first give the hierarchical spatial model of random auto-correlated eigenvalues in Section 2. The estimation of the proposed model is provided in Section 3. We give simulation studies and real data analysis in Section 4 and Section 5, respectively. Finally, we give a conclusion in Section 6. The theoretical details in our paper are in Appendix A. The codes for implementing our method are summarized in Appendix B.

In this section, we give the model construction step-by-step. First, we provide the basic model setup in Section 2.1, where the probability distribution of the eigenvalues and the covariate effect specification are introduced. Under the basic model, we further give a latent layer of eigenvalues that induce spatial correlation and provide the details in Section 2.2. To wrap up, we provide a model summary and introduce its application to the statistical inference of eigenvalues of DTI in Section 2.3.

Our model is based on the the parameterized Wishart distribution, and the parameterized Wishart distribution is slightly different from the Wishart distribution introduced in the classic textbooks. We give the construction of the parameterized Wishart distribution as follows. $\bm{Z}_{j}$ independently follows a mean-zero normal distribution with the covariance matrix $\bm{\Sigma}$ over $j\in\{1,...,U\}$ , denoted as $\bm{Z}_{j}\stackrel{{\scriptstyle ind}}{{\sim}}\mathcal{N}(\bm{0},\bm{\Sigma})$. $\bm{W}=\sum_{j=1}^{U}\bm{Z}_{j}\bm{Z}_{j}^{T}/U$ follows a parameterized Wishart distribution with the mean matrix $\bm{\Sigma}$ and the degrees of freedom $U$, denoted as $\bm{W}\sim\mathcal{W}_{p}(\bm{\Sigma},U)$. The density function of $\bm{W}$ is given in Appendix A.1.

In this section, we introduce the statistical estimation of Model 6. We use $\bm{l}_{i,g_{t}}=[l_{i,g_{t}}^{(1)},l_{i,g_{t}}^{(2)},l_{i,g_{t}}^{(3)}]$ to denote the observed eigenvalues of the random eigenvalues $\bm{L}_{i,g_{t}}=[L_{i,g_{t}}^{(1)},L_{i,g_{t}}^{(2)},L_{i,g_{t}}^{(3)}]$. We also give $\bm{\xi}_{i,g_{t}}=[\xi_{i,g_{t}}^{(1)},\xi_{i,g_{t}}^{(2)},\xi_{i,g_{t}}^{(3)}]$, $\bm{L}=\{\bm{L}_{i,g_{t}}:i=1,2,...,N;g=1,2,...,G;t=1,2,...,T_{g}\}$, and $\bm{\xi}=\{\bm{\xi}_{i,g_{t}}:i=1,2,...,N;g=1,2,...,G;t=1,2,...,T_{g}\}$. The parameters to be estimated are $\bm{\theta}=\{M,m,p,\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\beta}^{(3)}\}$ where $\bm{\beta}^{(k)}=\{\bm{B}_{g_{t}}^{(k)}=[{\beta}_{0,g_{t}}^{(k)},{\bm{\beta}_{g_{t}}^{(k)}}^{T}]^{T}:g=1,2,...,G;t=1,2,...,T_{g}\}$.

First, we aim to give the explicit expression of the log-likelihood of the data $\bm{l}$ and the latent variables $\bm{\xi}$, in order to proceed with the next parameter estimation. The relevant proofs and derivations in this section are provided in Appendix A.2.

The log-likelihood of the data $\bm{l}$ and the latent variables $\bm{\xi}$ is expressed as

where the term $f_{\bm{L}_{i,g_{t}}|\bm{\xi}_{i,g_{t}}}(\bm{l}_{i,g_{t}}|\bm{\xi}_{i,g_{t}};m,\bm{\beta})$ is the conditional density of $\bm{L}_{i,g_{t}}$ given $\bm{\xi}_{i,g_{t}}$ and the term $f_{\bm{\xi}_{g}}(\bm{\xi}_{i,g_{1}},...,\bm{\xi}_{i,g_{t}},...,\bm{\xi}_{i,g_{T_{g}}};M,p)$ is the joint density function of $[\bm{\xi}_{i,g_{1}},...,\bm{\xi}_{i,g_{t}},...,\bm{\xi}_{i,g_{T_{g}}}]$.

Muirhead (2009, Theorem 3.2.18) gives the expression of $f_{\bm{L}_{i,g_{t}}|\bm{\xi}_{i,g_{t}}}(\bm{l}_{i,g_{t}}|\bm{\xi}_{i,g_{t}};m,\bm{\beta})$. The density function is difficult to compute due to the presence of an integral which is difficult to evaluate. Given the assumption that $m$ is assumed to be large to make $\mathbb{E}_{L_{it}^{(k)}|\xi_{it}^{(k)}}L_{it}^{(k)}$ dominated by $\lambda_{i,g_{t}}^{(k)}=\xi_{i,g_{t}}^{(k)}e^{\beta_{0,g_{t}}^{(k)}+\bm{X}_{i}\bm{\beta}_{g_{t}}^{(k)}}$, we give Lemma 3.1 borrowed from Srivastava (1979, Corollary 9.4.2).

Therefore, to simplify the computation, we approximate $f_{\bm{L}_{i,g_{t}}|\bm{\xi}_{i,g_{t}}}(\bm{l}_{i,g_{t}}|\bm{\xi}_{i,g_{t}};m,\bm{\beta})$ as a product of probability density functions of normal distributed random variables, expressed as

where $\lambda_{i,g_{t}}^{(k)}=\xi_{i,g_{t}}^{(k)}e^{\beta_{0,g_{t}}^{(k)}+\bm{X}_{i}\bm{\beta}_{g_{t}}^{(k)}}$ and $\phi(l;\mu,\sigma)$ is the probability density function of a normal distribution with mean $\mu$ and standard error $\sigma$. In practice, this simplification reduces the computational cost significantly but preserves the accurate parameter estimation.

Next we handle the probability density function $f_{\bm{\xi}_{ig}}(\bm{\xi}_{i,g_{1}},...,\bm{\xi}_{i,g_{t}},...,\bm{\xi}_{i,g_{T_{g}}};M,p)$. In the following paragraphs, we derive the expression of the probability density function. First, Lemma 3.2 below allows us to handle only two density functions: $f_{\bm{\xi}_{i,g_{1}}}(\bm{\xi}_{i,g_{1}};M)$ and $f_{\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}};M,p)$.

Because $\bm{\xi}_{i,g_{1}}\sim\mathcal{E}([1,1,1],M)$, Muirhead (2009, Corollary 3.2.19) directly gives the expression of the probability density function $f_{\bm{\xi}_{i,g_{1}}}(\bm{\xi}_{i,g_{1}};M)$, expressed

Given James et al. (1960, Equation 7) and James et al. (1964, Equation 68), we further derives the expression of $f_{\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}};M,p)$ (see Appendix A.2), expressed as

where $\ {}_{0}F_{1}(;)$ is a hypergeometric function of a matrix argument and its value can be numerically evaluated (Koev and Edelman, 2006). When $p=0$, $\bm{\xi}_{i,g_{t}}$ is independent of $\bm{\xi}_{i,g_{t-1}}$. Furthermore, either $f_{\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}};M,p)=f_{\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}}|\bm{\xi}_{i,g_{t-1}};M,p)f_{\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t-1}};M)$ or $f_{\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}};M,p)=f_{\bm{\xi}_{i,g_{t-1}}|\bm{\xi}_{i,g_{t}}}(\bm{\xi}_{i,g_{t-1}}|\bm{\xi}_{i,g_{t}};M,p)f_{\bm{\xi}_{i,g_{t}}}(\bm{\xi}_{i,g_{t}};M)$ produces the joint density function of $\bm{\xi}_{i,g_{t}}$ and $\bm{\xi}_{i,g_{t-1}}$, denoted as $f_{\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}}}(\bm{\xi}_{i,g_{t}},\bm{\xi}_{i,g_{t-1}};M,p)$.

Putting Equations 8-10 together, we have the explicit expression of the full log-likelihood $\ell_{full}(\bm{\theta}|\bm{l},\bm{\xi})$ (Equation 7). Standing on the above results, we next proceed the statistical inference of correlated eigenvalues. Our goal is to obtain the maximum likelihood estimators of $\bm{\theta}$. In particular, to answer the clinical question about covariate effect quantification, the maximum likelihood estimate of $\bm{\beta}_{g_{t}}^{(k)}$ quantifies the covariate effects, and its confidence interval describes the uncertainties of covariate effects. Therefore, we have to maximize the likelihood whose the latent variables $\bm{\xi}$ are integrated out. The log-likelihood $\ell(\bm{\theta}|\bm{l})$ whose latent variables are integrated out is expressed as

The hurdle caused by the unobserved latent variables $\bm{\xi}$ can be easily overcome through the EM algorithm. The iterative scheme of the algorithm is given as follows. We start with the current parameter $\bm{\theta}^{(r)}$, thus the full likelihood at the current stage is $\ell_{Full}(\bm{\theta}^{(r)}|\bm{l},\bm{\xi})$. We first compute the expected value of the likelihood function with respect to the current conditional distribution of $\bm{\xi}$ given the data $\bm{l}$ and the current parameter estimates $\bm{\theta}^{(r)}$, denoted as $Q(\bm{\theta}|\bm{\theta}^{(r)})=\mathbb{E}_{\bm{\xi}|\bm{L},\bm{\theta}^{(r)}}\ell_{Full}(\bm{\theta}|\bm{l},\bm{\xi})$. This step is called E-step since we take the expectation. Next we maximize the term $Q(\bm{\theta}|\bm{\theta}^{(r)})$ with respect to $\bm{\theta}$ to obtain the next estimates $\bm{\theta}^{(r+1)}$. This step is called M-step because we maximize the function. We iteratively repeat the two steps until convergence, in order to obtain the maximum likelihood estimator $\hat{\bm{\theta}}$ for the parameters $\bm{\theta}$.

In the M-step, the function $Q(\bm{\theta}|\bm{\theta}^{(r)})$ can be feasibly maximized using the Quasi-Newton method with the the constraints ${\beta}_{dt}^{(1)}>{\beta}_{dt}^{(2)}>{\beta}_{dt}^{(3)}$ for $d=0,1,2,...,D$. In the E-step, the term $\mathbb{E}_{\bm{\xi}|\bm{L},\bm{\theta}^{(r)}}\ell_{Full}(\bm{\theta}|\bm{l},\bm{\xi})$ has no analytic expression. Thus, the Monte-Carlo EM based on importance sampling (Levine and Casella, 2001, Section 2.1) is applied. We first generate $C$ Markov chain Monte-Carlo samples from the posterior distribution $\bm{\xi}|\bm{L},\bm{\theta}^{(0)}$, denoted as $\bm{\xi}^{(1)},...,\bm{\xi}^{(C)}$, where $\bm{\theta}^{(0)}$ is the initial values of the EM algorithm. For every step $r$, $Q(\bm{\theta}|\bm{\theta}^{(r)})$ is approximated as $Q(\bm{\theta}|\bm{\theta}^{(r)})=\mathbb{E}_{\bm{\xi}|\bm{L},\bm{\theta}^{(r)}}\ell_{Full}(\bm{\theta}|\bm{l},\bm{\xi})\approx\frac{1}{C}\sum_{c=1}^{C}\omega_{c}\ell_{Full}(\bm{\theta}|\bm{l},\bm{\xi}^{(c)})$. The weight $\omega_{c}$ is defined as $\omega_{c}=\frac{\omega_{c}^{\prime}}{\sum_{c=1}^{C}\omega_{c}^{\prime}}$ where $\omega_{c}^{\prime}=\exp\Big{[}\ell_{Full}(\bm{\theta}^{(r)}|\bm{l},\bm{\xi})-\ell_{Full}(\bm{\theta}^{(0)}|\bm{l},\bm{\xi})\Big{]}$. The correctness and relevant properties of the Monte-Carlo EM based on the importance sampling have already been provided in Levine and Casella (2001).