[1]Rongke Lyu

1]Department of Statistics, Rice University, Houston, TX, United States

2]Department of Biostatistics, MD Anderson Cancer Center, Houston, TX, United
States

#  Bayesian Tensor Modeling for Image-based Classification of Alzheimer’s
Disease

<EMAIL_ADDRESS>Marina Vannucci  Suprateek Kundu  [  [

###### Abstract

Tensor-based representations are being increasingly used to represent complex
data types such as imaging data, due to their appealing properties such as
dimension reduction and the preservation of spatial information. Recently,
there is a growing literature on using Bayesian scalar-on-tensor regression
techniques that use tensor-based representations for high-dimensional and
spatially distributed covariates to predict continuous outcomes. However
surprisingly, there is limited development on corresponding Bayesian
classification methods relying on tensor-valued covariates. Standard
approaches that vectorize the image are not desirable due to the loss of
spatial structure, and alternate methods that use extracted features from the
image in the predictive model may suffer from information loss. We propose a
novel data augmentation-based Bayesian classification approach relying on
tensor-valued covariates, with a focus on imaging predictors. We propose two
data augmentation schemes, one resulting in a support vector machine (SVM)
type of classifier, and another yielding a logistic regression classifier.
While both types of classifiers have been proposed independently in
literature, our contribution is to extend such existing methodology to
accommodate high-dimensional tensor valued predictors that involve low rank
decompositions of the coefficient matrix while preserving the spatial
information in the image. An efficient Markov chain Monte Carlo (MCMC)
algorithm is developed for implementing these methods. Simulation studies show
significant improvements in classification accuracy and parameter estimation
compared to routinely used classification methods. We further illustrate our
method in a neuroimaging application using cortical thickness MRI data from
Alzheimer’s Disease Neuroimaging Initiative, with results displaying better
classification accuracy throughout several classification tasks.

###### keywords:

Alzheimer’s disease; Bayesian tensor modeling; logistic regression; support
vector machines; neuroimaging analysis.

##  1 Introduction

Neuroimaging studies stand as a cornerstone in contemporary neuroscience,
fundamentally transforming our comprehension of the intricate structure and
function of the brain. These non-invasive visualization techniques have not
only enriched our understanding of neurological disorders but have also
pioneered new frontiers in mental health research. Within the realm of risk
prediction, neuroimaging studies have emerged as an invaluable tool for
identifying individuals susceptible to neurological and psychiatric
conditions. By discerning subtle abnormalities in brain structure and
connectivity, researchers can now predict the risk of mental disorders (such
as Alzheimer’s disease or dementia) with greater precision. This early
identification facilitates timely intervention and treatment, potentially
offering an opportunity to provide improved health outcomes and quality of
life related to these disorders.

For example, neuroimaging studies have transformed the field of Alzheimer’s
disease (AD) research, providing invaluable insights into the pathological
mechanisms underlying this devastating neurodegenerative disorder (Chouliaras
and O’Brien, 2023). By visualizing the brain’s structural and functional
changes, neuroimaging techniques have enabled researchers to track the
progression of AD, identify early signs of the disorder, and differentiate it
from other causes of dementia. In particular, structural neuroimaging studies
involving magnetic resonance imaging (MRI), have revealed the intricate
patterns of atrophy in spatially distributed brain regions involved in memory
and cognition that is a hallmark of AD (Frenzel et al., 2020). Neuroimaging
features such as brain volume or cortical thickness that are derived from MRI
scans can be used to provide quantitative assessments of disease severity and
monitor disease progression over time. Importantly, such neuroimaging features
can be embedded in machine learning algorithms in order to perform risk
prediction or early detection in AD. However, several major challenges are
encountered when analyzing neuroimaging data. For example, the brain imaging
data is spatially dependent, high-dimensional and noisy, and it is often
unclear how to identify suitable neurobiological markers for the mental
disorder in the presence of heterogeneity.

In order to model such complex types of imaging data emerging at a rapid pace,
several statistical and machine learning approaches have been proposed. Among
them, classification models using neuroimaging features have seen a rapid
development (Rathore et al., 2017, Arbabshirani et al., 2017, Falahati et al.,
2014). These approaches typically either vectorize the image, or extract
informative summary features from the image, to be used as covariates. For
example, Plant et al. (2010) extracted the low-level-feature extraction
algorithm with feature selection criterion to select most discriminating
features, that are then coupled with a clustering algorithm to group spatially
coherent voxels to predict Alzheimer’s disease status. Ben Ahmed et al. (2015)
proposed a multi-feature fusion algorithm used both extracted visual features
from the hippocampal region of interest (ROI) and the quantity of
cerebrospinal fluid (CSF) in the hippocampal region, and then applied a late
fusion scheme to perform the binary classification of Alzheimer’s disease
subjects using the MRI images. Going beyond AD classification, Griffis et al.
(2016) implemented a voxel-based Gaussian Naive Bayes Classification of
ischemic stroke lesions in individual T1-weighted MRI scans, where the authors
separately created two feature maps as predictor variables for missing and
abnormal tissue to avoid including highly redundant information. Alternate
types of shape-based image analysis that go beyond voxel-level analysis have
also been proposed for prediction (Wu et al., 2022).

The above approaches, while useful, did not explicitly account for the spatial
configuration of imaging voxels. Some exceptions include Markov random field
(MRF) based methods that have been proposed in the prediction context (Smith
and Fahrmeir, 2007, Lee et al., 2014). However, given that these are not
equipped to perform dimension reduction, they may not be fully scalable to
high-dimensional images with tens of thousands of voxels, and their
performance in classification problems is unclear. In order to tackle the
spatial information in the image in the context of multi-class classification,
Pan et al. (2018) proposed a penalized linear discriminant analysis (LDA)
model using scalar and tensor covariates. Unfortunately, there is limited, if
any, literature on Bayesian classification approaches based on imaging
features that account for the spatial information in the image. This is
surprising, given the utility of Bayesian methods that can predict the
probability of an observation belonging to two or more classes, which can be
useful in the presence of measurement error or uncertainty regarding class
labels in medical imaging studies (Morales et al., 2013, Behler et al., 2022).
Existing Bayesian classification approaches that use vectorized features can
be not readily adapted to our problem of interest involving Bayesian image-
based classification, since it ignores the spatial structure of the image
resulting in information loss and potentially poor model performance.
Additionally, simply vectorizing the imaging features without an appropriate
lower dimensional representation also introduces the curse of dimensionality
since the number of voxels in the image is typically tens of thousands.
Alternate approaches that rely on first extracting lower dimensional features
from the image and subsequently using these features for classification, may
involve an additional layer of information loss resulting from the feature
extraction step, resulting in potential loss in accuracy.

Recently, there has been a growing literature on tensor analysis in
statistical modeling for imaging data that addresses some of the above
concerns. Guhaniyogi et al. (2017) proposed a Bayesian tensor regression with
a scalar response on scalar and tensor covariates. Other tensor models include
Bayesian response regression models that model the image outcome as a tensor
object. Guhaniyogi and Spencer (2021) implemented a Bayesian Tensor response
on scalar regression with an application of neuronal activation detection in
fMRI experiments using both tensor-valued brain images and scalar covariates.
Kundu et al. (2023) proposed a longitudinal Bayesian tensor response
regression model for mapping neuroplasticity across longitudinal visits.
Billio et al. (2023) proposed a novel linear autoregressive tensor process
model that introduces dynamics in linear tensor regression and allows for both
tensor-valued covariates and outcomes. Under the frequentist approach, Lock
(2018) proposed a penalized tensor on tensor regression using a (L2) Ridge
penalty. Zhou et al. (2013) proposed a tensor regression model by extending
generalized linear model to include tensor-structured covariates, where the
rank-R PARAFAC decomposition is assumed on the tensor parameters, with
adaptive lasso penalties applied on the tensor margins. With the exception of
the penalized GLM approach in Zhou et al. (2013), most existing tensor-based
approaches in literature have mainly focused on linear regression models that
can not be readily used for Bayesian classification.

In this article, we propose a data augmentation-based Bayesian classification
approach that models binary outcomes based on imaging covariates using a
tensor-based representation. We consider two different data augmentation
schemes resulting in two distinct Bayesian classifiers: a support vector
machine (SVM) and a logistic regression model. While these classifiers have
been extensively used in literature, the focus has been on using non-
structured covariates that ignore the spatial structure embedded in the image.
Our specific interest is in Bayesian classification based on imaging
predictors, where the images are registered across samples. Such a set-up is
routinely used in neuroimaging studies. Our main contribution is to develop a
Bayesian classification methodology based on high-dimensional tensor-valued
predictors via low-rank decompositions of the coefficient matrix and using
data augmentation. The low-rank PARAFAC decomposition assumed by the tensor
model is able to preserve the spatial configuration of imaging voxels, while
overcoming the challenges arising from the high dimensionality of the image
that can often contain tens of thousands of voxels. This results in
considerable improvements in classification accuracy, as illustrated via
rigorous numerical examples. In contrast to existing feature extraction
approaches that first use a tensor decomposition or alternate schemes to
obtain low level features to be subsequently used in modeling (Sen and Parhi,
2021), the proposed approach uses the full image as is in the classification
model, but employs a low rank PARAFAC decomposition to model the high-
dimensional tensor model coefficients. This ensures no information loss due to
feature extraction, while simultaneously allowing for dimension reduction. We
adopt the multiway shrinkage prior from Guhaniyogi et al. (2017) to model the
tensor margins of the assumed rank-R PARAFAC decomposition, which shrinks non-
significant parameters to near zero while inducing a minimal shrinkage effect
on the significant parameters. We develop efficient Markov chain Monte Carlo
(MCMC) algorithms for posterior inference that use data augmentation
techniques. Simulation studies show significant improvements in classification
accuracy, parameter estimation and feature selection compared to routinely
used classification methods that use vectorized images. We further illustrate
the advantages under our method via a detailed neuroimaging application using
voxel-wise cortical thickness features data from Alzheimer’s Disease
Neuroimaging Initiative (ADNI) study, with the proposed Bayesian classifiers
displaying better out-of-sample accuracy consistently throughout several
classification tasks. We leverage cortical thickness as our neuroimaging
feature of choice, since it is known to be a highly sensitive imaging
biomarker for modeling neurodegeneration in AD (Weston et al., 2016, Fjell et
al., 2015).

The rest of the paper is organized as follows: In section 2, we propose the
framework of the Bayesian tensor classification model with two types of data
augmentation, specify the prior and hyperparameter choices, and list the
posterior computation steps. In Section 3 we study the model performances
through comprehensive simulation studies. In Section 4 we provide the results
from the data analysis using the ADNI dataset. We conclude the manuscript with
a discussion.

##  2 Methods

###  2.1 Brief Introduction to Tensors

Tensor-based models have gained recognition as a promising way to model
neuroimaging data, due to their multifold advantages. Tensors naturally
inherit a multidimensional structure to represent complex data structures such
as spatial features of a brain region. Additionally, tensor-based techniques
achieve dimension reduction, which is particularly useful with neuroimaging
data to tackle the challenges of $p>>n$ in statistical modeling. A tensor is a
multi-dimensional array, with the order being the number of its dimensions.
For example, a one-way or first-order tensor is a vector, and a second-order
tensor is a matrix. A fiber resembles the idea of matrix rows and columns in
higher dimensionality, which is obtained by fixing every dimension of a tensor
except one. Similarly, a slice is defined by fixing every order of the tensor
except two. Tensor decomposition is a mathematical technique that expresses a
high dimensional tensor into the combination of lower dimensional factors. One
type of tensor decomposition is the Tucker decomposition (Kolda and Bader,
2009), which decomposes a tensor into a core tensor and a set of matrices, one
along each mode. It can be denoted as follows:

$\displaystyle\bm{B}$
$\displaystyle=\Lambda\times_{1}A\times_{2}B\times_{3}\cdots\times_{D}D$
$\displaystyle=\sum_{r_{1}=1}^{R_{1}}\cdots\sum_{r_{D}=1}^{R_{D}}\lambda_{r_{1},\ldots,r_{D}}\bm{a}^{\left(r_{1}\right)}\circ\cdots\circ\bm{d}^{\left(r_{D}\right)},$
(1)

where $\Lambda$ is the core tensor, and $A$, $B$, $D$ are the factor matrices.
The PARAFAC decomposition is a special case of the Tucker decomposition, where
the core tensor $\Lambda$ is restricted to be diagonal, and
$R_{1}=R_{2}=\ldots=R_{D}=R$. The rank-R PARAFAC model can then be expressed
as

$\displaystyle\bm{B}=\sum_{r=1}^{R}\bm{\beta}_{1}^{(r)}\circ\cdots\circ\bm{\beta}_{D}^{(r)},$
(2)

where $\bm{\beta}_{1},\ldots,\bm{\beta}_{D}$, known as tensor margins, are
vectors of length $p_{1},\ldots,p_{D}$, and where
$\bm{\beta}_{1}\circ\cdots\circ\bm{\beta}_{D}$ is a D-way outer product of
dimension $p_{1}\times p_{2}\times\ldots\times p_{D}$. It is essential to
recognize that tensor margins can only be uniquely identified up to a
permutation and a multiplicative constant unless we introduce additional
constraints. However, the lack of identifiability in tensor margins does not
create any complications for our scenario. This is because the tensor product
is fully identifiable, which suffices for our primary objective of estimating
coefficients. Consequently, we refrain from imposing extra identifiability
conditions on the tensor margins, aligning with the principles in the Bayesian
tensor modeling literature (Guhaniyogi, 2020). Furthermore, the PARAFAC
decomposition dramatically reduces the number of coefficients from
$p_{1}\times\ldots\times p_{D}$ to $R(p_{1}+\ldots+p_{D})$, which grows
linearly with the tensor rank $R$ and results in significant dimension
reduction. The appropriate tensor rank can vary depending on the specific
application context and can be chosen using a goodness-of-fit approach.

Prior to applying the tensor model, the image’s voxels are transformed onto a
regularly spaced grid, making them more suitable for a tensor-based approach.
This mapping conserves the spatial arrangements of the voxels, offering
notable advantages over a univariate voxel-wise analysis or a multivariable
analysis that vectorizes the voxels without considering their spatial
arrangements. While the grid mapping might not preserve exact spatial
distances between voxels, this has limited impact, as it can still capture
correlations between neighboring elements in the tensor margins. Moreover, the
tensor construction has the advantageous ability to estimate voxel-specific
coefficients by leveraging information from neighboring voxels through the
estimation of tensor margins with their inherent low-rank structure. This
feature results in brain maps which are more consistent and robust to missing
voxels and image noise. Additionally, it can be conveniently used for reliable
imputation of imaging features related to missing voxels. In contrast, voxel-
wise analysis lacks the capacity to share information among neighboring
voxels, treating voxel coefficients as independent entities regardless of
their spatial arrangement. For more details on the advantages and
characteristics of Bayesian tensor models, we refer readers to Kundu et al.
(2023).

###  2.2 Loss-function Based Classification

A loss function is a mathematical tool to quantify the difference between
predicted values under a model and the actual observed data values. Loss
functions are often used for finding optimal estimates of model parameters in
machine learning and statistical modeling tasks for regression,
classification, and more. In the Bayesian paradigm, loss functions translate
to different types of likelihood that are combined with additional priors on
the model parameters embedded in the loss function to obtain posterior
distributions that are subsequently used for estimation and uncertainty
quantification. Although we motivate our approach by drawing connections with
loss functions, we note that our work is distinct compared to decision
theoretic Bayesian approaches that use loss functions as a post-processing
step after Markov chain Monte Carlo (MCMC) to derive optimal estimates. See,
for example, Hahn and Carvalho (2015) and Kundu et al. (2019). In this article
we focus, in particular, on two types of commonly used loss functions: the
hinge loss represented as the support vector machine classifier, and the
logistic regression loss, both of which employ high-dimensional images as
covariates for classification.

Figure 1: From left to right: hinge loss and sigmoid loss function

Support vector machine (SVM): SVMs play a pivotal role in classification
tasks, and their significance stems from their ability to handle complex
decision boundaries with remarkable efficiency. SVMs work by finding an
optimal hyperplane that separates the data points into two classes. This
hyperplane is chosen to maximize the margin, which is the distance between the
hyperplane and the closest data points. This helps to ensure that the SVM
model is generalizable to new data, and protects against overfitting. SVMs are
well-suited for both linear and non-linear classification, via using suitable
kernel functions. This versatility makes SVMs applicable across various
domains, from image recognition and natural language processing to
bioinformatics - see Cervantes et al. (2020) for a review. Their robust
performance, ability to manage high-dimensional data, and capacity to handle
intricate relationships between features underscore their importance in
tackling diverse and challenging classification problems in machine learning.

Most SVM-based classifiers rely on point estimates with penalized approaches
to tackle high-dimensional covariates (Peng et al., 2016, Dedieu, 2019). The
SVM classifier uses the hinge loss function that takes the form

$\displaystyle\mathcal{L}(y\mid\beta)=\frac{1}{\sigma^{2}}max(1-yf({\bf
x};{\bm{\beta}}),0)+R,$ (3)

where $y\in\\{-1,1\\}$ is the binary outcome, $f(\cdot)$ is a linear or non-
linear function of (potentially high-dimensional) covariates ${\bf x}$, with
corresponding unknown parameters ${\bm{\beta}}$ that need to be estimated from
the data, and additional tuning parameters $\sigma^{2}$. See Figure 1 for a
visualization of the hinge loss function. Recently, Ma and Kundu (2022)
generalized the hinge loss to a smooth hinge loss and derived provably
flexible estimators in the presence of noisy high-dimensional covariates in
the SVM framework. Bayesian approaches for SVM based on a pseudo-likelihood
approach were proposed originally by Polson and Scott (2011) and subsequently
adopted for biomedical applications in Sun et al. (2018). In particular, the
pseudo-likelihood can be represented as a location-scale mixture of normals
with latent variable $\rho$ as

$\displaystyle L=\prod_{i=1}^{n}L_{i}(y_{i}\mid{\bf
x}_{i},{\bm{\beta}},{\sigma^{2}})=\prod_{i=1}^{n}\big{\\{}\frac{1}{\sigma^{2}}\exp\\{-\frac{2}{\sigma^{2}}\max(1-y_{i}f({\bf
x};{\bm{\beta}}),0)\\}\big{\\}}$
$\displaystyle=\int_{0}^{\infty}\prod_{i=1}^{n}\frac{1}{{\sigma}\sqrt{2\pi\rho_{i}}}\exp(-\frac{(1+\rho_{i}-y_{i}f({\bf
x};{\bm{\beta}}))^{2}}{2\rho_{i}{\sigma^{2}}})d\rho_{i},$ (4)

where $L_{i}$ represents the contribution corresponding to the $i$th sample.
The above representation essentially uses data augmentation techniques,
introducing a latent variable $\rho$ which, when marginalized over, gives back
the hinge loss function. Such a latent variable representation enables an
efficient Gibbs sampler for posterior inference, which will be described in
Section 2.4 below.

Logistic regression classifier: Logistic regression is a powerful and
versatile machine learning algorithm that is widely used in diverse
applications, including medical diagnosis, fraud detection, customer
segmentation, marketing campaigns, and recommendation systems. It is
characterized by simplicity in the implementation and interpretability of the
results, as the logistic regression coefficients can be understood in terms of
odd-ratios. In order to make this model scalable for high dimensional
biomedical applications, penalized versions of logistic regression models have
been proposed (Doerken et al., 2019, Devika et al., 2016).

The logistic loss function specifies a sigmoid type loss that takes the
analytical form

$\displaystyle\mathcal{L}(y=1\mid{\bm{\beta}})=\exp\\{f({\bf
x};{\bm{\beta}})\\}/\big{(}1+\exp\\{f({\bf x};{\bm{\beta}})\\}\big{)},$ (5)

where $f({\bf x};{\bm{\beta}})$ represents the contribution of the covariates
to the logistic loss that is quantified via the unknown parameters
${\bm{\beta}}$, to be estimated from the data. See Figure 1 for a
visualization of the sigmoid loss function. The binary outcome variables
$y\in\\{0,1\\}$ are assumed to follow a Bernoulli distribution, with the
corresponding probability function resembling the form
$\mathcal{L}(y=1\mid{\bm{\beta}})$. Typically $f(\cdot)$ represents linear
functions of covariates that facilitate straightforward interpretations for
the model parameters, although non-linear logistic regression models have also
been proposed (Tokdar and Ghosh, 2007). In the Bayesian paradigm, model
fitting and inference for logistic regression is often achieved using Polya-
Gamma latent variables (Nicholas G. Polson and Windle, 2013). A random
variable $X$ has a Polya–Gamma distribution with parameters $b>0$ and
$c\in\mathcal{R}$, denoted as $X\sim PG(b,c)$ if

$\displaystyle X\stackrel{{\scriptstyle
D}}{{=}}\frac{1}{2\pi^{2}}\sum_{k=1}^{\infty}\frac{g_{k}}{(k-1/2)^{2}+c^{2}/\left(4\pi^{2}\right)},$
(6)

where $g_{k}$ follows a Gamma distribution $Ga(b,1)$. The introduction of the
Polya-Gamma latent variable allows one to represent the binomial likelihood as
mixtures of Gaussians. It can be shown that the logistic loss function can be
recovered by marginalizing out the latent Polya-Gamma latent variable using
the following relationship

$\displaystyle\frac{\big{(}e^{f(\cdot)})\big{)}^{y}}{\big{(}1+e^{f(\cdot)}\big{)}^{b}}=2^{-b}e^{\kappa\psi}\int_{0}^{\infty}e^{-\omega\psi^{2}/2}p(\omega)d\omega,\mbox{
}b>0,\kappa=y-b/2,$ (7)

where $\omega\sim PG(b,0)$, and $f({\bf x};{\bm{\beta}})$ is a linear
predictor in most cases. The above identity facilitates conjugate updates
under a Gaussian prior, conditional upon the latent Polya-Gamma variable, as
detailed in Nicholas G. Polson and Windle (2013). The full data augmented
likelihood is given by the following expression:

$\displaystyle
L=\prod_{i=1}^{n}\frac{\left(e^{f_{i}}\right)^{y_{i}}}{\left(1+e^{f_{i}}\right)}=\prod_{i=1}^{n}2^{-1}e^{\kappa_{i}\psi_{i}}\int_{0}^{\infty}e^{-\omega_{i}f_{i}^{2}/2}p(\omega_{i})d\omega_{i},$
(8)

where $\kappa_{i}=y_{i}-1/2$, $b=1$, $\omega_{i}\sim PG(1,0)$.

###  2.3 Priors

For our implementation, we consider the widely used linear predictor using the
relationship

$\displaystyle
f_{i}=\langle\bm{X}_{i},\bm{B}\rangle+\bm{z}_{i}^{\prime}{\bm{\gamma}},$ (9)

where $\bm{X}_{i}$ and $\bm{z}_{i}$ denote the imaging predictors and the
supplemental (eg: demographic/clinical) features, respectively, for the $i$th
sample, $\langle\cdot,\cdot\rangle$ denotes the inner product operator,
$\bm{B}$ denotes the tensor-valued coefficient matrix that quantifies the
effect of the image on the classification model, and ${\bm{\gamma}}$ is a
vector of dimension $p_{z}+1$ capturing the effects of the supplemental
covariates. Furthermore, we assume
$\bm{B}\in\otimes^{D}_{j=1}\mathbb{R}^{p_{j}}$, that is modeled under a
PARAFAC decomposition as in (2).

For the prior choice on the tensor margins, we adopt the multiway Dirichlet
generalized double Pareto (M-DGDP) prior from Guhaniyogi et al. (2017), which
shrinks small coefficients towards zero while minimizing shrinkage of large
coefficients. The prior can be expressed in hierarchical form on the tensor
coefficient margins $\bm{\beta}_{j}^{(r)}$, $j=1,\ldots,D$ and $r=1,\ldots,R$
as:

$\displaystyle\bm{\beta}_{j}^{(r)}\sim\mathrm{N}\left(0,\left(\phi_{r}\tau\right)\bm{W}_{jr}\right),\mbox{
}w_{jr,k}\sim\operatorname{Exp}\left(\lambda_{jr}^{2}/2\right),$ (10)

where $\tau\sim Ga(a_{\tau},b_{\tau})$ is a global scale parameter,
$\Phi=\left(\phi_{1},\ldots,\phi_{R}\right)$ follows a Dirichlet distribution
that encourages shrinkage to lower rank in the assumed PARAFAC decomposition
with $\left(\phi_{1},\ldots,\phi_{R}\right)\sim
Dirichlet\left(\alpha_{1},\ldots,\alpha_{R}\right)$, and where
$\bm{W}_{jr}=\mbox{Diag}(w_{jr,1},\ldots,w_{jr,p_{j}})$ are scale parameters
that are margin-specific for each element and modeled under an Exponential
distribution as $w_{jr,k}\sim Exp(\lambda^{2}_{jr}/2)$ with $\lambda_{jr}$
unknown and modeled as
$\lambda_{jr}\sim\operatorname{Ga}\left(a_{\lambda},b_{\lambda}\right)$.
Additionally, one can obtain

$\beta_{j,k}^{(r)}\mid\lambda_{jr},\phi_{r},\tau\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{DE}\left(\lambda_{jr}/\sqrt{\phi_{r}\tau}\right),1\leq
k\leq p_{j},$ (11)

after marginalizing out the scale parameters $w_{jr,k}$, that is, prior (10)
induces a GDP prior on the individual margin coefficients which in turn has
the form of an adaptive Lasso penalty as in Armagan et al. (2013). Overall,
the flexibility in estimating $\bm{B}_{r}=\left\\{\bm{\beta}_{j}^{(r)};1\leq
j\leq D\right\\}$ is accommodated by component-specific scaling parameter
$w_{jr,k}$ and common rate parameter $\lambda_{jr}$, which shares information
between margin elements and encourages shrinkage at the local scale. We
complete the prior specification by assuming a $N(0,\Sigma_{0\gamma})$ as the
prior on $\mathbf{\gamma}$.

###  2.4 MCMC Algorithms for Posterior Inference

For posterior inference, we implemented efficient MCMC algorithms that take
advantage of data augmentation techniques. Algorithm 1 outlines the MCMC
updates for the proposed Bayesian Tensor SVM model (BT-SVM), and Algorithm 2
illustrates the updating scheme for the the proposed Bayesian Tensor logistic
regression model (BT-LR).

1:Update $\rho_{i}$ from inverse Gaussian: $\rho_{i}^{-1}\sim
IN(\mu_{i},\lambda_{i})$, where
$\mu_{i}=|1-y_{i}(<X_{i},B>+z^{\prime}\gamma)|^{-1}$ and
$\lambda_{i}=1/\sigma^{2}$

2:Update hyperparameters $[\alpha,\Phi,\tau\mid\bm{B},\bm{W}]$ compositionally
as $[\alpha\mid\bm{B},\bm{W}][\Phi,\tau\mid\alpha,\bm{B},\bm{W}]$, as
described in Guhaniyogi et al. (2017).

3:Sample $\\{\beta_{j}^{(r)},\omega_{jr},\lambda_{jr}\\}$ using a back-fitting
procedure to produce a sequence of draws from the margin-level conditional
distributions across components.

  * (a)

Draw
$\left[w_{jr},\lambda_{jr}\mid\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]=\left[w_{jr}\mid\lambda_{jr},\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]\left[\lambda_{jr}\mid\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]$.

    1. 1.

Draw
$\lambda_{jr}\sim\mathrm{Ga}\left(a_{\lambda}+p_{j},b_{\lambda}+\left\|\bm{\beta}_{j}^{(r)}\right\|_{1}/\sqrt{\phi_{r}\tau}\right)$;

    2. 2.

Draw
$w_{jr,k}\sim\operatorname{giG}\left(\frac{1}{2},\lambda_{jr}^{2},\beta_{j,k}^{2(r)}/\left(\phi_{r}\tau\right)\right)$
independently for $1\leq k\leq p_{j}$

  * (b)

Draw $\beta_{j}^{(r)}$ from multivariate normal distribution:
$\beta_{j}^{(r)}\sim N(\mu_{jr},\Sigma_{jr})$, where
$\mu_{jr}=\frac{\Sigma_{jr}H_{j}^{(r)}\tilde{\mathbf{y}}}{\sigma^{2}}$ and
$\Sigma_{jr}={(\frac{H_{j}^{(r)T}H_{j}^{(r)}}{\sigma^{2}}+\frac{W_{jr}^{-1}}{\phi_{r}\tau})}^{-1}$,
with

$h_{i,j,k}^{(r)}=\sum_{d_{1}=1,\ldots,d_{D}=1}^{p_{1},\ldots,p_{D}}I\left(d_{j}=k\right)x_{d_{1},\ldots,d_{D}}\left(\prod_{l\neq
j}\beta_{l,i_{l}}^{(r)}\right),$
$\bm{H}_{i,j}^{(r)}=\left(h_{i,j,1}^{(r)}/\sqrt{\rho_{i}},\ldots,h_{i,j,p_{j}}^{(r)}/\sqrt{\rho_{i}}\right)^{\prime},$
$\tilde{y}_{i}=\frac{y_{i}}{\sqrt{\rho_{i}}}(\rho_{i}+1-y_{i}(z_{i}^{\prime}\gamma+\sum_{l\neq
r}<X_{i},B_{l}>))$

4:Update $\gamma$ from a conjugate normal conditional distribution $\gamma\sim
N(\mu_{\gamma},\Sigma_{\gamma})$ where
$\mu_{\gamma}=\Sigma_{\gamma}Z^{T}(\tilde{y}y/\rho)$ and
$\Sigma_{\gamma}=(G^{T}G/\sigma^{2}+\Sigma_{0\gamma})^{-1}$, with
$\bm{G}_{i,p_{z}}=Z_{i,p_{z}}/\sqrt{\rho_{i}}$ and
$\tilde{y_{i}}=\rho_{i}+1-y_{i}<X_{i},B>$.

Algorithm 1 MCMC steps for BT-SVM

1:Update $\omega_{i}$ from Polya-gamma distribution:
$\omega_{i}\sim\mathrm{PG}(1,<X_{i},\bm{B}>+z_{i}^{T}\gamma)$

2:Update hyperparameters $[\alpha,\Phi,\tau\mid\bm{B},\bm{W}]$ compositionally
as $[\alpha\mid\bm{B},\bm{W}][\Phi,\tau\mid\alpha,\bm{B},\bm{W}]$, as
described in Guhaniyogi et al. (2017).

3:Sample $\\{\beta_{j}^{(r)},\omega_{jr},\lambda_{jr}\\}$ using a back-fitting
procedure to produce a sequence of draws from the margin-level conditional
distributions across components.

  * (a)

Draw
$\left[w_{jr},\lambda_{jr}\mid\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]=\left[w_{jr}\mid\lambda_{jr},\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]\left[\lambda_{jr}\mid\bm{\beta}_{j}^{(r)},\phi_{r},\tau\right]$.

    1. 1.

Draw
$\lambda_{jr}\sim\mathrm{Ga}\left(a_{\lambda}+p_{j},b_{\lambda}+\left\|\bm{\beta}_{j}^{(r)}\right\|_{1}/\sqrt{\phi_{r}\tau}\right)$;

    2. 2.

Draw
$w_{jr,k}\sim\operatorname{giG}\left(\frac{1}{2},\lambda_{jr}^{2},\beta_{j,k}^{2(r)}/\left(\phi_{r}\tau\right)\right)$
independently for $1\leq k\leq p_{j}$

  * (b)

Draw $\beta_{j}^{(r)}$ from multivariate normal distribution:
$\beta_{j}^{(r)}\sim N(\mu_{jr},\Sigma_{jr})$, with
$\mu_{jr}=\Sigma_{jr}(\Omega(H_{j}^{(r)})^{T}\tilde{y})$ and
$\Sigma_{jr}={((H_{j}^{(r)})^{T}\Omega
H_{j}^{(r)}+W_{jr}^{-1}/(\phi_{r}\tau))}^{-1},$ where
$\tilde{y}=\kappa/\omega$,
$\kappa=\left(y_{1}-N_{1}/2,\ldots,y_{n}-N_{n}/2\right)$,
$N_{1}=\ldots=N_{n}=1$, and $\Omega$ a diagonal matrix with diagonal elements
$\omega_{i}$’s.

4:Update $\gamma$ from a conjugate normal conditional distribution $\gamma\sim
N(\mu_{\gamma},\Sigma_{\gamma})$, where
$\mu_{\gamma}=\Sigma_{\gamma}Z^{T}(\tilde{y}\omega)$ and
$\Sigma_{\gamma}=(G^{T}G+\Sigma_{0\gamma})^{-1}$, with
${G}_{i,p_{z}}=Z_{i,p_{z}}*\sqrt{\omega_{i}}$ and
$\tilde{y_{i}}=\kappa_{i}/\omega_{i}-<X_{i},B>$.

Algorithm 2 MCMC steps for BT-LR

##  3 Simulation Study

###  3.1 Data Generation

We illustrate the performance of our methods using several simulation settings
and perform comparisons with competitive approaches, based on various types of
generated data sets involving several types of functional signals and with
data generated from SVM and logistic loss functions. We considered four
different types of signals for the tensor coefficient $\bm{B}$ to generate the
binary outcome, as defined below.

Scenario 1: In this setting, the tensor $\bm{B}$ is constructed from rank-R
PARAFAC decomposition with rank $R_{0}=3$ and dimension $p=c(48,48)$. Each
beta margin $\beta_{j}^{(r)}$ is generated from the independent binomial
distribution $Binomial(2,0.2)$. After the construction of tensor, we set the
maximal value of the tensor $\bm{B}$ cells to be 1.

Scenario 2: The tensor image is simulated by a rank-R PARAFAC decomposition
with rank $R_{0}=3$. Here, instead of generating the tensor margin from a
known distribution, we manually set up each value of $\beta_{j}^{(r)}$.

Scenario 3: Instead of generating the 2D tensor images from a PARAFAC
decomposition, the tensor coefficient $\bm{B}$ is set to be 1 for a
rectangular area and 0 otherwise. The non-zero elements cover approximately 30
percent of the area.

Scenario 4: The tensor coefficient $\bm{B}$ is set to be 1 for the circular
area and 0 otherwise. The non-zero elements cover approximately 10 percent of
the area.

The top panel in Figure 2 shows the true 2D tensor images, for the different
scenarios. For each setting, we generated the tensor covariates $\bm{X}$ from
standard normal distribution $N(0,1)$. For simplicity, we did not include non-
tensor covariates in our simulation settings, i.e. we assumed the true
$\gamma=(0,\ldots,0)^{\prime}$. Finally for each scenario, the binary outcome
$Y$ was generated according to both the SVM and logistic loss function as
follows. Denoting the linear prediction as $\psi=<\bm{X}_{i},\bm{B}>$, the
binary outcome is generated as $Y_{i}=1$ if $\psi>0$ and $Y_{i}=-1$ otherwise,
under the SVM loss, and from a Bernoulli distribution with probability
$p=1/(1+exp(-\psi)$ under a logistic loss. Therefore a total of 8 scenarios
are considered in our simulation set-up, and for each of these scenarios 10
replicates were generated.

###  3.2 Parameter Settings

We choose suitable values of hyperparameters in the prior distributions that
yield good overall performance. For example, we set the parameters of the
hyperprior on the global scale $\tau$ to $a_{\tau}=1$ and $b_{\tau}=\alpha
R^{(1/D)}$, where $R$ is the rank in the assumed PARAFAC decomposition, and
set $\alpha_{1}=\ldots=\alpha_{R}=1/R$. For the common rate parameter
$\lambda_{jr}$, we set $a_{\lambda}=3$ and
$b_{\lambda}=\sqrt[2D]{a_{\lambda}}$. Note that under the SVM loss, the
scaling parameter $\sigma^{2}$ is a fixed parameter that can be manually
adjusted for maximal model performance. Several values of the tuning parameter
$\sigma^{2}$ from 0.1 through 10 are tested and we choose $\sigma^{2}=6$. In
order to decide the rank of the fitted model, we fit the proposed model using
ranks 2-5 and choose the rank that minimizes the Deviance Information
Criterion (DIC) scores. DIC measures the goodness-of-fit of a set of Bayesian
hierarchical models adjusting for model complexity in a manner that penalizes
more complex models.

###  3.3 Performance Evaluation

We report estimation accuracy in terms of relative error(RE), Root Mean
Squared Error (RMSE) and correlation coefficient for point estimation of cell-
level tensor coefficients. We also illustrate the classification accuracy by
calculating misclassification error and the F1 score. Feature selection
performances are evaluated by Sensitivity, Specificity, F1 score, and Matthews
correlation coefficient (MCC). These metrics are defined as follows. Let
$\theta_{j},~{}j=1,\ldots,J$, be the vectorized tensor coefficients, with
$J=\prod_{k=1}^{D}p_{k}$ the total number of cells in the tensor coefficient
$\bm{B}$. Further, define the following terms related to classification
performance under the SVM classifier: (a) TP is the true positive, i.e. the
number of predictions where the classifier correctly predicts the positive
class as positive; (b) FP is the false positive, i.e. the number of
predictions where the classifier incorrectly predicts the negative class as
positive; (c) TN is the true negative, i.e. the number of predictions where
the classifier correctly predicts a negative class as negative; and (d) FN is
a false negative, i.e. the number of predictions where the classifier
incorrectly classifies a positive class as negative. For logistic
classification, the above definitions also hold, but after replacing the
negative class with zero class. The above definitions of TP/FP/TN/FN can also
be adopted for feature selection performance such that the positive class
corresponds to non-zero coefficients, and the negative/zero class refers to
absent or zero coefficients.

Metrics for evaluating coefficient estimation performance: These metrics
include: (i) Relative error of $\theta$, denoted as
$RE(\theta)=\frac{\sum_{i=1}^{J}|\hat{\theta_{i}}-\theta_{i}|}{\sum_{i=1}^{J}|\theta_{i}|}$,
that measures the scaled absolute deviations between the true and estimated
parameters, with $\hat{\theta_{i}}$ and $\theta_{i}$ the estimated and true
coefficients, respectively. We note that the values of relative error can be
larger than 1, since the discrepancy between estimated and true values of
tensor cells can be greater than zero when a large portion of the true values
is exactly 0. (ii) Root-mean-square error of $\theta$, denoted as
$RMSE(\theta)=\sqrt{\frac{\sum_{i=1}^{J}(\hat{\theta_{i}}-\theta_{i})^{2}}{J}}$,
that provides another measure of estimation accuracy; (iii) correlation
coefficient between the true and estimated coefficients.

Metrics for evaluating classification performance: These metrics include: (i)
Mis-classification rate, defined as $\frac{(FP+FN)}{TP+TN+FP+FN}$; and (ii)
F1-score, defined as the harmonic mean between precision (i.e. TP/(TP+FP)) and
recall or sensitivity (TP/(TP+FN)). The expression of F1-score is given as
$\frac{TP}{TP+(FP+FN)/2}$.

Metrics for evaluating feature selection performance: These metrics include:
(i) Sensitivity = $\frac{TP}{TP+FN}$; (ii) Specificity = $\frac{TN}{TN+FP}$;
and (iii) MCC = $\frac{TP\times TN-FP\times
FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$.

Results reported below were obtained by randomly splitting the data into
training and test sets in the ratio 70:30. The metrics for point estimation
and feature selection performances were calculated using the training set
data, while the metrics for classification performance were calculated based
on the test set. We report below averaged values of the selected metrics
across 10 replicates. Two state-of-the-art classification methods are used as
competing methods. The first is a penalized logistic regression model with
lasso penalty, which is available in the R package glmnet (Friedman et al.,
2010). The second competing method is the L1-norm SVM model from the R package
penalizedSVM (Becker et al., 2009, Bradley and Mangasarian, 1998). Both
methods use a vectorization approach, where the tensor covariates are
vectorized into scalar variables and then fit into the statistical models.
Therefore they do not respect the spatial information in the image.
Additionally, we use a grid search algorithm and cross-validation to select
the best tuning parameters prior to model fitting.

###  3.4 Results

We ran MCMC chains for 3,000 iterations, with 1,000 burn-in iterations. The
computation time varies depending on the rank and is expected to increase with
higher ranks. It took around 19 minutes to run a single MCMC chain with rank
2, and around 36 minutes with rank 4. A Geweke diagnostic (Geweke, 1991) was
applied to examine for signs of non-convergence parameters. We obtain the
z-score from Geweke for each element of the coefficient matrix $\bm{B}$. For
the proposed Bayesian tensor SVM model (BT-SVM), we observe the z-scores lie
in the range $(-1.96,1.96)$ for 91 percent of the coefficient matrix elements.
For the proposed Bayesian tensor logistic regression model (BT-LR), the
z-scores lie in the range $(-1.96,1.96)$ for 79 percent, indicating that most
chains have reached ergodicity.

Table 1: Point Estimation and Out-of-sample Classification results for the four 2D tensor images portrayed in Figure 2 top panel; Y generated from SVM loss Scenarios | Methods | RE | RMSE | Corr.Coef. | Mis. Class. | F1-score
---|---|---|---|---|---|---
Scenario 1 | LR w/ lasso | 1.000 | 0.558 | 0.071 | 0.54 | 0.147
| L1norm-SVM | 1.110 | 0.687 | 0.027 | 0.513 | 0.522
| BT-SVM | 0.986 | 0.489 | 0.442 | 0.260 | 0.764
| BT-LR | 1.055 | 0.504 | 0.383 | 0.342 | 0.689
Scenario 2 | LR w/ Lasso | 1.000 | 0.382 | 0.055 | 0.46 | 0
| L1norm-SVM | 1.226 | 0.533 | 0.092 | 0.48 | 0.502
| BT-SVM | 0.853 | 0.246 | 0.874 | 0.149 | 0.837
| BT-LR | 1.006 | 0.277 | 0.794 | 0.204 | 0.778
Scenario 3 | LR w/ Lasso | 1.001 | 0.541 | 0.074 | 0.507 | 0.191
| L1norm-SVM | 1.032 | 0.537 | 0.109 | 0.533 | 0.512
| BT-SVM | 0.836 | 0.412 | 0.810 | 0.197 | 0.823
| BT-LR | 0.909 | 0.439 | 0.733 | 0.238 | 0.781
Scenario 4 | LR w/ Lasso | 1.005 | 0.330 | 0.122 | 0.533 | 0.2
| L1norm-SVM | 1.335 | 0.535 | 0.053 | 0.44 | 0.565
| BT-SVM | 1.013 | 0.225 | 0.773 | 0.200 | 0.819
| BT-LR | 1.447 | 0.262 | 0.614 | 0.272 | 0.752
|  |  |  |  |  |
Table 2: Feature selection results for the four 2D tensor images portrayed in Figure 2 top panel; Y generated from SVM loss Scenarios | Methods | Sens. | Spec. | F1-score | MCC
---|---|---|---|---|---
Scenario 1 | LR w/ lasso | 0.031 | 0.984 | 0.057 | 0.046
| L1norm-SVM | 0.058 | 0.950 | 0.1 | 0.017
| BT-SVM | 0.070 | 0.999 | 0.130 | 0.213
| BT-LR | 0.214 | 0.938 | 0.313 | 0.225
Scenario 2 | LR w/ Lasso | 0.003 | 1 | 0.007 | 0.055
| L1norm-SVM | 0.021 | 0.985 | 0.037 | 0.016
| BT-SVM | 0.628 | 1 | 0.771 | 0.772
| BT-LR | 0.821 | 0.978 | 0.835 | 0.814
Scenario 3 | LR w/ Lasso | 0.025 | 0.981 | 0.047 | 0.019
| L1norm-SVM | 0.019 | 0.988 | 0.037 | 0.027
| BT-SVM | 0.192 | 1 | 0.32 | 0.377
| BT-LR | 0.517 | 0.987 | 0.668 | 0.626
Scenario 4 | LR w/ Lasso | 0.094 | 0.974 | 0.145 | 0.122
| L1norm-SVM | 0.020 | 0.984 | 0.034 | 0.009
| BT-SVM | 0.368 | 0.998 | 0.530 | 0.573
| BT-LR | 0.519 | 0.934 | 0.552 | 0.522
|  |  |  |  |
Table 3: Point Estimation and Out-of-sample Classification results for the four 2D tensor images portrayed in Figure 2 top panel; Y generated from logistic regression loss Scenarios | Methods | RE | RMSE | Corr.Coef. | Mis. Class. | F1-score
---|---|---|---|---|---|---
Scenario 1 | LR w/ lasso | 1 | 0.557 | 0.076 | 0.56 | 0.125
| L1norm-SVM | 1.017 | 0.562 | 0.054 | 0.48 | 0.55
| BT-SVM | 0.991 | 0.493 | 0.495 | 0.274 | 0.745
| BT-LR | 1.010 | 0.490 | 0.445 | 0.247 | 0.767
Scenario 2 | LR w/ Lasso | 1.000 | 0.382 | 0.085 | 0.46 | 0
| L1norm-SVM | 1.034 | 0.383 | 0.053 | 0.413 | 0.537
| BT-SVM | 0.875 | 0.245 | 0.859 | 0.179 | 0.804
| BT-LR | 1.196 | 0.260 | 0.732 | 0.252 | 0.734
Scenario 3 | LR w/ Lasso | 1.000 | 0.542 | 0.071 | 0.506 | 0.380
| L1norm-SVM | 1.068 | 0.587 | 0.041 | 0.433 | 0.586
| BT-SVM | 0.836 | 0.412 | 0.804 | 0.178 | 0.842
| BT-LR | 0.964 | 0.421 | 0.618 | 0.253 | 0.770
Scenario 4 | LR w/ Lasso | 1.002 | 0.331 | 0.131 | 0.467 | 0.557
| L1norm-SVM | 1.361 | 0.542 | 0.058 | 0.467 | 0.557
| BT-SVM | 1.017 | 0.226 | 0.766 | 0.203 | 0.813
| BT-LR | 1.202 | 0.255 | 0.658 | 0.227 | 0.788
|  |  |  |  |  |
Table 4: Feature selection results for the four 2D tensor images portrayed in Figure 2 top panel; Y generated from logistic regression loss Scenarios | Methods | Sens. | Spec. | F1-score | MCC
---|---|---|---|---|---
Scenario 1 | LR w/ lasso | 0.036 | 0.979 | 0.067 | 0.045
| L1norm-SVM | 0.023 | 0.982 | 0.044 | 0.020
| BT-SVM | 0.007 | 0.999 | 0.014 | 0.061
| BT-LR | 0.343 | 0.913 | 0.446 | 0.318
Scenario 2 | LR w/ Lasso | 0.017 | 0.999 | 0.034 | 0.098
| L1norm-SVM | 0.042 | 0.978 | 0.070 | 0.042
| BT-SVM | 0.590 | 0.999 | 0.736 | 0.742
| BT-LR | 0.801 | 0.948 | 0.764 | 0.736
Scenario 3 | LR w/ Lasso | 0.024 | 0.984 | 0.044 | 0.026
| L1norm-SVM | 0.010 | 0.981 | 0.020 | -0.031
| BT-SVM | 0.180 | 1 | 0.300 | 0.361
| BT-LR | 0.475 | 0.926 | 0.573 | 0.483
Scenario 4 | LR w/ Lasso | 0.086 | 0.972 | 0.133 | 0.102
| L1norm-SVM | 0.040 | 0.983 | 0.067 | 0.051
| BT-SVM | 0.409 | 0.999 | 0.577 | 0.612
| BT-LR | 0.609 | 0.965 | 0.645 | 0.610
|  |  |  |  |

We report results for estimation and classification accuracy, and feature
selection in Tables 1, 2, 3, 4. Specifically, Tables 1 and 2 reflect results
across all 4 scenarios when the binary outcome Y is generated from a SVM loss,
while Tables 3 and 4 reflect the logistic type of loss. These results
demonstrate that both the proposed methods (BT-SVM and BT-LR) consistently
outperform competing penalized methods in terms of coefficient estimation,
feature selection, and classification performance across all scenarios. When
the binary outcome data is generated from SVM loss, the BT-SVM approach has
superior coefficient estimation (as evident from lower RE/RMSE and higher
correlation coefficient in Table 1) and improved classification accuracy (as
evident from lower misclassification rate and higher F1-score in Table 1).
Even when the data is generated from a logistic loss, the same trends hold,
with the exception of Scenario 1 where the proposed BT-LR approach reports
improved classification accuracy and comparable coefficient estimation, as
evident from the results in Table 3.

Moreover, the proposed model with SVM loss (BT-SVM) generally performs worse
than the corresponding model with logistic loss (BT-LR) in terms of feature
selection, as evident from the results in Tables 2 and 4. In particular, the
BT-LR approach almost always has considerably higher sensitivity compared to
the BT-SVM model while having comparable or slightly lower specificity, even
when the outcome data is generated under the SVM loss. This results in the BT-
LR approach consistently having higher F1-score and higher or comparable MCC
values for feature selection, when compared with the BT-SVM approach, even
when the outcome data is generated under the SVM loss.

Figure 2: Row 1 from left: Simulated data with 48$\times$48 2D tensor images
from Scenario 1, Scenario 2, Scenario 3, and Scenario 4. Row 2: Recovered
images for the 48$\times$48 2D tensor images using BT-SVM corresponding to 4
scenarios in row 1. Row 3: Recovered images for the 48$\times$48 2D tensor
images using BT-LR corresponding to 4 scenarios in row 1.

Combining the above discussions, the BT-SVM model appears to generally show
improved coefficient estimation and classification performance over its
counterpart with logistic loss, regardless of which type of loss function was
used to generate the outcome data. On the contrary, the BT-LR model appears to
have improved sensitivity and comparable specificity compared to its
counterpart with SVM loss that translates to improved feature selection,
regardless of which loss function is used to generate the outcome data. Figure
2 presents the coefficient matrix cell estimation using the proposed methods
BT-SVM and BT-LR. From this Figure, it is evident that the proposed method is
able to broadly recover the shapes of the 2D tensor $\bm{B}$ regardless of
whether the underlying signal is generated using the PARAFAC decomposition or
not.

Figure 3: Row 1 from left: Simulated data with 48$\times$48 2D tensor images
from Scenario 1, Scenario 2, Scenario 3, and Scenario 4. Row 2: Recovered
images for the 48$\times$48 2D tensor images using competing method L1norm-SVM
for Scenario 1, Scenario 2, Scenario 3, and Scenario 4. Row 3: Recovered
images for the 48$\times$48 2D tensor images using logistic regression with
LASSO penalty.

In contrast, the penalized competing methods (logistic regression with LASSO
penalty and L1norm-SVM) illustrate poor estimation, feature selection, and
out-of-sample classification performance. This is evident from the results in
Tables 2, 4, where the competing approaches report notably low sensitivity,
almost approaching zero. This low sensitivity reflects their inability to
detect the true signals, resulting in poor coefficient estimation performance
as evident from Tables 1, 3, where the correlation coefficients between the
true and estimated coefficients are often close to zero. For a more detailed
perspective, Figure 3 visualizes the estimated coefficients under two
competing methods. The Figure reveals spatially disparate non-zero
coefficients, indicating that the absence of spatial smoothing hampers the
accurate estimation of true signals by these competing methods. Ultimately,
the substandard feature selection performance contributes to considerably
inferior classification results, as demonstrated in Tables 1, 3.

###  3.5 Sensitivity Analysis

Guhaniyogi et al. (2017) set up a series of default values for the prior
hyper-parameters for the Bayesian tensor regression model. Specifically,
hyperparameter of the Dirichlet component
$\left(\alpha_{1},\ldots,\alpha_{R}\right)=\alpha=(1/R,...,1/R)$,
$a_{\lambda}=3$, $b_{\lambda}=\sqrt[2D]{a_{\lambda}}$,
$a_{\tau}=\sum_{i=1}^{R}\alpha_{i}$, and $b_{\tau}=\alpha R^{(1/D)}$ are set
as defaulted values to control the cell-level variance on tensor $\bm{B}$. We
conducted a prior hyperparameter sensitivity analysis by tuning each
hyperparameter while fixing other hyperparameters. Table 5 displays the cell-
level RMSE values of the tensor coefficient matrix $\bm{B}$ of dimension $D=2$
and parafac rank-R, where for BT-SVM model, $\bm{B}$ is generated from
Scenario 4; and for BT-LR model, $\bm{B}$ is generated from Scenario 3. The
overall results reveal no strong sensitivity to the hyper-parameter choices
under the selected ranges that point to a robust performance.

Table 5: Prior hyperparameter sensitivity analysis | Hyperparameters | RMSE | Hyperparameters | RMSE | Hyperparameters | RMSE
---|---|---|---|---|---|---
BT-SVM |  $\alpha$ = 1/9 | 0.224 |  $a_{\lambda}$ = 3 | 0.225 |  $a_{\tau}$ = 1/3 | 0.225
$\alpha$= 1/6 | 0.228 |  $a_{\lambda}$ = 5 | 0.235 |  $a_{\tau}$ = 1/2 | 0.225
|  $\alpha$ = 1/3 | 0.226 |  $a_{\lambda}$ = 7 | 0.227 |  $a_{\tau}$ = 1 | 0.226
| $\alpha=3^{(-0.1)}$ | 0.227 |  $a_{\lambda}$ = 10 | 0.230 |  $a_{\tau}$ = 2 | 0.231
BT-LR |  $\alpha$ = 1/9 | 0.418 |  $a_{\lambda}$ = 3 | 0.413 |  $a_{\tau}$ = 1/3 | 0.428
$\alpha$= 1/6 | 0.422 |  $a_{\lambda}$ = 5 | 0.421 |  $a_{\tau}$ = 1/2 | 0.416
|  $\alpha$ = 1/3 | 0.419 |  $a_{\lambda}$ = 7 | 0.413 |  $a_{\tau}$ = 1 | 0.416
| $\alpha=3^{(-0.1)}$ | 0.421 |  $a_{\lambda}$ = 10 | 0.415 |  $a_{\tau}$ = 2 | 0.420
|  |  |  |  |  |

##  4 ADNI Data Analysis

###  4.1 Data Source and Pre-processing

This study utilized the data obtained from the Alzheimer’s Disease
Neuroimaging Initiative (ADNI), which is a longitudinal multicenter study
launched in 2004 for the early detection and tracking of Alzheimer’s disease
(AD). ADNI researchers collect multiple data types such as clinical,
behavioral and genetic data, along with neuroimaging measurements such as
magnetic resonance imaging (MRI), positron emission tomography (PET), and
biospecimens. We use data from the ADNI 1 study collected at baseline,
consisting of T1-weighted MRI scans, cognitive measurements in terms of Mini
Mental State examination (MMSE), and basic demographic data (age, gender, year
of education, APOE status) of 818 subjects. A more detailed description of the
ADNI data is provided in Table 6.

Table 6: Summary of demographic variables and cognitive measurements under study |  |  Overall (N = 818)
---|---|---
Age | Mean (SD) | 75.17 (6.83)
Median[Min, Max] | 75.45[54.40, 90.90]
Years of Education | Mean (SD) | 15.53 (3.05)
Median [Min, Max] | 16 [4, 20]
Gender | Female | 342 (41.8%)
Male | 476 (58.2%)
APOE4 | 0 | 418 (51.1%)
1 | 312 (31.2%)
2 | 88 (10.7%)
MMSE score | Mean(SD) | 26.74 (2.67)
Median[Min, Max] | 27 [18, 30]
log(ICV score) | Mean (SD) | 14.25 (0.11)
Median [Min, Max] | 24.25 [13.92, 14.56]
|  |

The T1-weighted MRI images underwent processing through the Advance
Normalization Tools (ANTs) registration pipeline (Tustison et al., 2014),
where all images were registered to a template image to ensure consistent
normalization of brain locations across participants. The population-based
template was constructed based on data from 52 normal control participants in
ADNI 1, originally from the ANTs group (Tustison et al., 2019). Notably, the
ANTs pipeline includes the N4 bias correction step, addressing intensity
discordance to inherently standardize intensity across samples (Tustison et
al., 2010). It also employs a symmetric diffeomorphic image registration
algorithm for spatial normalization, aligning each T1 image with a brain image
template to facilitate spatial comparability (Avants et al., 2008).
Subsequently, the pipeline utilized the processed brain images, estimated
brain masks, and template tissue labels for 6-tissue Atropos segmentation,
generating tissue masks for cerebrospinal fluid (CSF), gray matter (GM), white
matter (WM), deep gray matter (DGM), brain stem, and cerebellum. Finally,
cortical thickness measurements were derived using the DiReCT algorithm. The
3-D cortical thickness image was further downsampled to a dimension of
$48\times 48\times 48$ and divided into 48 2-D axial slices of dimension
$48\times 48$, and a subset of these 2D slices were used for our analysis. The
downsampling step reduces the dimension of the image, and also somewhat
alleviates the sparsity of the cortical thickness maps by consolidating
adjacent voxels and presenting the average cortical thickness. A reduction of
sparsity proves beneficial to fitting our Bayesian tensor model, and the same
should hold of other commonly used statistical models.

###  4.2 Analysis Outline

We apply the proposed approaches to perform various classification tasks using
data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) study.
Demographic information of age, gender, years of education, APOE4 allele
(0,1,2) and intracranial volume (ICV) are incorporated as scalar covariates,
while 2D cortical thickness slices (derived from the T1 weighted MRI scans)
are used as tensor covariates. Since the different 2D brain slices are
expected to contain varying amount of the brain cortical regions, it is
important to choose these slices carefully. In particular, we would like to
choose 2D slices such that they contain at least a certain portion of the
brain cortex, in order to contain enough information to perform
classification. Therefore, we present the analysis results from 7 different
axial slices (we will denote them slices 19 to 25), each of which has cortical
brain regions that cover at least 65% of the slice. We evaluate the
classification performance of the proposed approaches and compare with the
penalized logistic regression with lasso and L1-norm SVM as described in
Section 3.

In particular, we assess the ability of the proposed classifiers to
differentiate between different types of disease phenotypes using the 2D
imaging slices along with demographic information. In particular, we perform
the following types of classification tasks corresponding to disease
phenotypes: (i) normal control vs. AD patients; (ii) normal control vs. MCI
patients; (iii) MCI vs. AD patients. In addition to the above tasks, we also
perform (iv) gender classification (males vs females); and (v) intelligence
level classification based on high and low levels of MMSE scores. MMSE is
commonly used clinically for checking cognitive impairment with a low value
indicative of a cognitive decline. We stratify individuals into high vs low
MMSE categories, depending on whether their MMSE scores are above the 70th
percentile or below the 30th percentile of the MMSE distribution. The
corresponding sample sizes in the high and low intelligence categories were
280 and 249 respectively. For each classification task based on a given 2D
slice, the data is randomly split into training and test splits in the ratio
70:30, and 10 such splits are considered. The results are averaged over these
10 replicates and reported in Table 7.

###  4.3 Results

Table 7: Classification on Gender, Disease phenotype, and High/Low MMSE scores with Demographic covariates. | BT-SVM | BT-LR | LR-lasso | L1norm-SVM
---|---|---|---|---
| Mis. C.111Misclassification Rate. | F1-score | Mis. C. | F1-score | Mis. C. | F1-score | Mis. C. | F1-score
| Female vs. Male
slice 19 | 0.288 | 0.764 | 0.317 | 0.728 | 0.337 | 0.750 | 0.422 | 0.636
slice 20 | 0.280 | 0.762 | 0.292 | 0.748 | 0.361 | 0.729 | 0.398 | 0.654
slice 21 | 0.296 | 0.752 | 0.3 | 0.756 | 0.357 | 0.75 | 0.447 | 0.595
slice 22 | 0.284 | 0.779 | 0.333 | 0.723 | 0.361 | 0.729 | 0.504 | 0.550
slice 23 | 0.256 | 0.789 | 0.325 | 0.733 | 0.373 | 0.657 | 0.390 | 0.684
slice 24 | 0.280 | 0.771 | 0.325 | 0.720 | 0.398 | 0.611 | 0.447 | 0.618
slice 25 | 0.256 | 0.801 | 0.345 | 0.719 | 0.369 | 0.674 | 0.455 | 0.591
| NC. vs. AD
slice 19 | 0.327 | 0.669 | 0.277 | 0.715 | 0.388 | 0.595 | 0.452 | 0.564
slice 20 | 0.301 | 0.712 | 0.325 | 0.682 | 0.341 | 0.644 | 0.365 | 0.671
slice 21 | 0.341 | 0.656 | 0.309 | 0.697 | 0.341 | 0.638 | 0.404 | 0.564
slice 22 | 0.278 | 0.724 | 0.325 | 0.655 | 0.333 | 0.655 | 0.444 | 0.582
slice 23 | 0.277 | 0.720 | 0.309 | 0.698 | 0.293 | 0.654 | 0.436 | 0.444
slice 24 | 0.293 | 0.689 | 0.333 | 0.681 | 0.325 | 0.601 | 0.341 | 0.626
slice 25 | 0.269 | 0.746 | 0.309 | 0.677 | 0.301 | 0.660 | 0.396 | 0.510
| NC. vs MCI
slice 19 | 0.280 | 0.8 | 0.312 | 0.768 | 0.344 | 0.732 | 0.414 | 0.686
slice 20 | 0.291 | 0.791 | 0.380 | 0.707 | 0.302 | 0.751 | 0.418 | 0.691
slice 21 | 0.338 | 0.761 | 0.365 | 0.718 | 0.333 | 0.731 | 0.440 | 0.688
slice 22 | 0.317 | 0.771 | 0.349 | 0.736 | 0.333 | 0.722 | 0.458 | 0.664
slice 23 | 0.322 | 0.776 | 0.349 | 0.766 | 0.365 | 0.696 | 0.455 | 0.586
slice 24 | 0.328 | 0.760 | 0.349 | 0.750 | 0.360 | 0.704 | 0.433 | 0.620
slice 25 | 0.291 | 0.792 | 0.354 | 0.729 | 0.333 | 0.720 | 0.465 | 0.582
| AD vs. MCI
slice 19 | 0.296 | 0.805 | 0.355 | 0.746 | 0.334 | 0.763 | 0.395 | 0.690
slice 20 | 0.282 | 0.816 | 0.305 | 0.784 | 0.325 | 0.577 | 0.420 | 0.495
slice 21 | 0.287 | 0.815 | 0.338 | 0.781 | 0.350 | 0.759 | 0.412 | 0.691
slice 22 | 0.282 | 0.810 | 0.322 | 0.753 | 0.322 | 0.778 | 0.356 | 0.770
slice 23 | 0.271 | 0.825 | 0.344 | 0.751 | 0.384 | 0.723 | 0.401 | 0.702
slice 24 | 0.282 | 0.812 | 0.338 | 0.758 | 0.350 | 0.747 | 0.446 | 0.663
slice 25 | 0.291 | 0.808 | 0.389 | 0.721 | 0.339 | 0.765 | 0.463 | 0.613
| High vs Low MMSE scores
slice 19 | 0.283 | 0.723 | 0.327 | 0.675 | 0.333 | 0.693 | 0.496 | 0.606
slice 20 | 0.295 | 0.715 | 0.358 | 0.655 | 0.333 | 0.674 | 0.471 | 0.534
slice 21 | 0.289 | 0.722 | 0.352 | 0.654 | 0.345 | 0.667 | 0.484 | 0.549
slice 22 | 0.295 | 0.710 | 0.339 | 0.686 | 0.377 | 0.634 | 0.471 | 0.460
slice 23 | 0.295 | 0.712 | 0.345 | 0.645 | 0.327 | 0.653 | 0.421 | 0.645
slice 24 | 0.289 | 0.720 | 0.333 | 0.675 | 0.321 | 0.622 | 0.459 | 0.587
slice 25 | 0.301 | 0.707 | 0.327 | 0.662 | 0.327 | 0.653 | 0.471 | 0.590
|  |  |  |  |  |  |  |

Table 7 reports the misclassification rate and f1-score for slices 19 through
25. The proposed approach under both loss functions almost always results in
better classification accuracy consistently across slices and classification
tasks, compared to the penalized approaches. While both the two competing
approaches have inferior performance compared to their Bayesian counterparts,
the LR-Lasso generally performs slightly better than the L1norm-SVM. Among the
three disease phenotype classification tasks, the highest accuracy is achieved
for AD vs MCI classification under the BT-SVM approach, with f1-scores greater
than 0.8 across all 2D slices. The proposed approaches are also able to
perform well for the NC vs MCI classification task with F1-scores consistently
greater than 0.75, while the classification performance for NC vs AD is
slightly less impressive under all methods. We note that slice 23 provides the
highest accuracy for both AD vs MCI and NC vs MCI classification tasks. With
regards to the classification of non-imaging phenotypes, the gender
classification performance is generally improved compared to intelligence
classification, which is not surprising given the considerably lower training
sample size for classifying the high versus low intelligence cohorts.

Clearly, there exist fluctuations in classification performances across the
different 2-D slices. This is not surprising because cortical thickness
contained in each 2-D slice varies depending on the sections of the slice.
Further, different slices represent different brain regions that may have
differential effects on classifying phenotypic classes. Overall, these results
demonstrate that our proposed Bayesian tensor approach offers consistent
improvements for classifying Alzheimer’s disease and subjects’ demographic
information, by leveraging the spatial information in the images and
incorporating dimension reduction that potentially avoids overfitting.
Moreover, the unsatisfactory performance of the competing methods clearly
illustrates the perils of ignoring the spatial correlations in high-
dimensional images in classification tasks. This is not surprising, given that
Lasso-based approaches are known to be affected by multicollinearity. While it
is possible to explore alternate approaches such as principal component
regression (PCR), such approaches lead to loss of interpretability that is not
desirable in imaging studies due to various reasons including the inability to
perform feature selection.

##  5 Discussion

In this study, we proposed a Bayesian tensor classification approach based on
high-dimensional neuroimaging data and scalar predictors, via data
augmentation. The proposed approach essentially extends the literature on
Bayesian tensor regression to classification problems that has the ability to
perform inference and uncertainty quantification. The tensor structure is
particularly suitable for neuroimaging data since it respects the spatial
information of imaging voxels while keeping the number of parameters to be
estimated at a manageable level via a PARAFAC decomposition. With two data
augmentation schemes implemented, we demonstrated the superiority of the
proposed method via simulation and data application. In particular,
comprehensive simulation studies concretely illustrate the advantages under
the proposed approach, with each type of data augmentation having distinct
advantages in classification, coefficient estimation and/or feature selection.
We applied the proposed method to the ADNI dataset to classify disease
phenotype, as well as gender and intelligence levels. Both data augmentation
schemes showed consistently higher classification accuracy and improved
feature selection than other penalized methods that vectorized the image. This
showed the benefits of incorporating the spatial information in the image.

A potential limitation of the study is that we use 2D brain slices instead of
using the whole 3D brain image. While using the full image may potentially
result in improvements in accuracy, it can also cause computational
bottlenecks due to a massive number of model parameters. In future work, we
intend to propose more scalable versions of the proposed approach that can be
used to incorporate high-dimensional 3D images for classification. We also
propose to explore various prior choices for the tensor coefficients in order
to enable stronger shrinkage and thicker tails, with a view to improving the
feature selection performance. Another potential limitation is that we did not
consider potential heterogeneity across samples in our modeling scheme. In
future work, we plan to account for the heterogeneity across samples via a
mixture of tensors approach that is particularly relevant for AD studies.
While we plan to investigate the above issues in future research, we believe
that the proposed approach in this article provides a valuable Bayesian
classification tool based on imaging covariates that fills an important gap in
literature.

##  6 Declarations

Information Sharing Statement

We will make our code available on Github upon acceptance of the paper. The
sample data from ADNI can be accessed by accepting the Data Use Agreement and
submitting an online application form at https://adni.loni.usc.edu/data-
samples/access-data/.

Author Contributions S.K. designed the research and conceived the analysis
plan; R.L. performed model implementation and data analyses; S.K. and M.V.
provided overall supervision and direction of the work. R.L., S.K. and M.V.
wrote the paper.

Competing Interests The authors have no relevant financial or non-financial
interests to disclose. The authors have no competing interests to declare that
are relevant to the content of this article. All authors certify that they
have no affiliations with or involvement in any organization or entity with
any financial interest or non-financial interest in the subject matter or
materials discussed in this manuscript. The authors have no financial or
proprietary interests in any material discussed in this article.

## References

  * Arbabshirani et al. (2017) M. R. Arbabshirani, S. Plis, J. Sui, and V. D. Calhoun.  Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls.  _Neuroimage_ , 145:137–165, 2017.
  * Armagan et al. (2013) A. Armagan, D. B. Dunson, and J. Lee.  Generalized double Pareto shrinkage.  _Statistica Sinica_ , 23(1):119, 2013.
  * Avants et al. (2008) B. B. Avants, C. L. Epstein, M. Grossman, and J. C. Gee.  Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain.  _Medical image analysis_ , 12(1):26–41, 2008.
  * Becker et al. (2009) N. Becker, W. Werft, G. Toedt, P. Lichter, and A. Benner.  penalizedSVM: a R-package for feature selection SVM classification.  _Bioinformatics_ , 25(13):1711–1712, 2009.
  * Behler et al. (2022) A. Behler, H.-P. Müller, A. C. Ludolph, D. Lulé, and J. Kassubek.  A multivariate Bayesian classification algorithm for cerebral stage prediction by diffusion tensor imaging in amyotrophic lateral sclerosis.  _NeuroImage: Clinical_ , 35:103094, 2022.  10.1016/j.nicl.2022.103094.
  * Ben Ahmed et al. (2015) O. Ben Ahmed, J. Benois-Pineau, M. Allard, C. Ben Amar, G. Catheline, and A. D. N. Initiative.  Classification of Alzheimer’s disease subjects from MRI using hippocampal visual features.  _Multimedia Tools and Applications_ , 74:1249–1266, 2015.
  * Billio et al. (2023) M. Billio, R. Casarin, M. Iacopini, and S. Kaufmann.  Bayesian dynamic tensor regression.  _Journal of Business & Economic Statistics_, 41(2):429–439, 2023.
  * Bradley and Mangasarian (1998) P. S. Bradley and O. L. Mangasarian.  Feature selection via concave minimization and support vector machines.  In _ICML_ , volume 98, pages 82–90, 1998.
  * Cervantes et al. (2020) J. Cervantes, F. Garcia-Lamont, L. Rodríguez-Mazahua, and A. Lopez.  A comprehensive survey on support vector machine classification: Applications, challenges and trends.  _Neurocomputing_ , 408:189–215, 2020.
  * Chouliaras and O’Brien (2023) L. Chouliaras and J. T. O’Brien.  The use of neuroimaging techniques in the early and differential diagnosis of dementia.  _Molecular Psychiatry_ , pages 1–14, 2023.
  * Dedieu (2019) A. Dedieu.  Error bounds for sparse classifiers in high-dimensions.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 48–56. PMLR, 2019.
  * Devika et al. (2016) S. Devika, L. Jeyaseelan, and G. Sebastian.  Analysis of sparse data in logistic regression in medical research: A newer approach.  _Journal of postgraduate medicine_ , 62(1):26, 2016.  10.4103/0022-3859.173193.
  * Doerken et al. (2019) S. Doerken, M. Avalos, E. Lagarde, and M. Schumacher.  Penalized logistic regression with low prevalence exposures beyond high dimensional settings.  _PLoS One_ , 14(5):e0217057, 2019.  10.1371/journal.pone.0217057.
  * Falahati et al. (2014) F. Falahati, E. Westman, and A. Simmons.  Multivariate data analysis and machine learning in Alzheimer’s disease with a focus on structural magnetic resonance imaging.  _Journal of Alzheimer’s disease_ , 41(3):685–708, 2014.
  * Fjell et al. (2015) A. M. Fjell, H. Grydeland, S. K. Krogsrud, I. Amlien, D. A. Rohani, L. Ferschmann, A. B. Storsve, C. K. Tamnes, R. Sala-Llonch, P. Due-Tønnessen, et al.  Development and aging of cortical thickness correspond to genetic organization patterns.  _Proceedings of the National Academy of Sciences_ , 112(50):15462–15467, 2015.
  * Frenzel et al. (2020) S. Frenzel, K. Wittfeld, M. Habes, J. Klinger-Koenig, R. Buelow, H. Voelzke, and H. J. Grabe.  A biomarker for Alzheimer’s disease based on patterns of regional brain atrophy.  _Frontiers in psychiatry_ , 10:953, 2020.
  * Friedman et al. (2010) J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of statistical software_ , 33(1):1, 2010.
  * Geweke (1991) J. Geweke.  Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments.  Staff Report 148, Federal Reserve Bank of Minneapolis, 1991.  URL https://ideas.repec.org/p/fip/fedmsr/148.html.
  * Griffis et al. (2016) J. C. Griffis, J. B. Allendorfer, and J. P. Szaflarski.  Voxel-based Gaussian naïve Bayes classification of ischemic stroke lesions in individual T1-weighted MRI scans.  _Journal of neuroscience methods_ , 257:97–108, 2016.
  * Guhaniyogi (2020) R. Guhaniyogi.  Bayesian methods for tensor regression.  _Wiley StatsRef: Statistics Reference Online_ , pages 1–18, 2020.
  * Guhaniyogi and Spencer (2021) R. Guhaniyogi and D. Spencer.  Bayesian tensor response regression with an application to brain activation studies.  _Bayesian Analysis_ , 16(4):1221–1249, 2021.
  * Guhaniyogi et al. (2017) R. Guhaniyogi, S. Qamar, and D. B. Dunson.  Bayesian tensor regression.  _The Journal of Machine Learning Research_ , 18(1):2733–2763, 2017.
  * Hahn and Carvalho (2015) P. R. Hahn and C. M. Carvalho.  Decoupling shrinkage and selection in Bayesian linear models: a posterior summary perspective.  _Journal of the American Statistical Association_ , 110(509):435–448, 2015.
  * Kolda and Bader (2009) T. G. Kolda and B. W. Bader.  Tensor decompositions and applications.  _SIAM review_ , 51(3):455–500, 2009.
  * Kundu et al. (2019) S. Kundu, B. K. Mallick, and V. Baladandayuthapan.  Efficient Bayesian regularization for graphical model selection.  _Bayesian Analysis_ , 14(2):449, 2019.
  * Kundu et al. (2023) S. Kundu, A. Reinhardt, S. Song, J. Han, M. L. Meadows, B. Crosson, and V. Krishnamurthy.  Bayesian longitudinal tensor response regression for modeling neuroplasticity.  _Human Brain Mapping_ , 2023.
  * Lee et al. (2014) K.-J. Lee, G. L. Jones, B. S. Caffo, and S. S. Bassett.  Spatial Bayesian variable selection models on functional magnetic resonance imaging time-series data.  _Bayesian Analysis (Online)_ , 9(3):699, 2014.
  * Lock (2018) E. F. Lock.  Tensor-on-tensor regression.  _Journal of Computational and Graphical Statistics_ , 27(3):638–647, 2018.
  * Ma and Kundu (2022) X. Ma and S. Kundu.  High-dimensional measurement error models for lipschitz loss, 2022.  Preprint at https://arxiv.org/abs/2210.15008.
  * Morales et al. (2013) D. A. Morales, Y. Vives-Gilabert, B. Gómez-Ansón, E. Bengoetxea, P. Larrañaga, C. Bielza, J. Pagonabarraga, J. Kulisevsky, I. Corcuera-Solano, and M. Delfino.  Predicting dementia development in Parkinson’s disease using Bayesian network classifiers.  _Psychiatry Research: NeuroImaging_ , 213(2):92–98, 2013.
  * Nicholas G. Polson and Windle (2013) J. G. S. Nicholas G. Polson and J. Windle.  Bayesian inference for logistic models using Pólya–Gamma latent variables.  _Journal of the American Statistical Association_ , 108(504):1339–1349, 2013.
  * Pan et al. (2018) Y. Pan, Q. Mai, and X. Zhang.  Covariate-adjusted tensor classification in high dimensions.  _Journal of the American statistical association_ , 2018.
  * Peng et al. (2016) P. Peng, K. Lekadir, A. Gooya, L. Shao, S. E. Petersen, and A. F. Frangi.  A review of heart chamber segmentation for structural and functional analysis using cardiac magnetic resonance imaging.  _Magnetic Resonance Materials in Physics, Biology and Medicine_ , 29:155–195, 2016.
  * Plant et al. (2010) C. Plant, S. J. Teipel, A. Oswald, C. Böhm, T. Meindl, J. Mourao-Miranda, A. W. Bokde, H. Hampel, and M. Ewers.  Automated detection of brain atrophy patterns based on MRI for the prediction of Alzheimer’s disease.  _Neuroimage_ , 50(1):162–174, 2010.  10.1016/j.neuroimage.2009.11.046.
  * Polson and Scott (2011) N. G. Polson and S. L. Scott.  Data augmentation for support vector machines.  _Bayesian Analysis_ , 6:1–24, 2011.  URL http://ba.stat.cmu.edu/journal/2011/vol06/issue01/polson.pdf.
  * Rathore et al. (2017) S. Rathore, M. Habes, M. A. Iftikhar, A. Shacklett, and C. Davatzikos.  A review on neuroimaging-based classification studies and associated feature extraction methods for Alzheimer’s disease and its prodromal stages.  _NeuroImage_ , 155:530–548, 2017.
  * Sen and Parhi (2021) B. Sen and K. K. Parhi.  Predicting biological gender and intelligence from fMRI via dynamic functional connectivity.  _IEEE Transactions on Biomedical Engineering_ , 68(3):815–825, 2021.  10.1109/TBME.2020.3011363.
  * Smith and Fahrmeir (2007) M. Smith and L. Fahrmeir.  Spatial Bayesian variable selection with application to functional magnetic resonance imaging.  _Journal of the American Statistical Association_ , 102(478):417–431, 2007.
  * Sun et al. (2018) W. Sun, C. Chang, Y. Zhao, and Q. Long.  Knowledge-guided Bayesian support vector machine for high-dimensional data with application to analysis of genomics data.  In _2018 IEEE International Conference on Big Data (Big Data)_ , pages 1484–1493. IEEE, 2018.
  * Tokdar and Ghosh (2007) S. T. Tokdar and J. K. Ghosh.  Posterior consistency of logistic Gaussian process priors in density estimation.  _Journal of statistical planning and inference_ , 137(1):34–42, 2007.
  * Tustison et al. (2010) N. J. Tustison, B. B. Avants, P. A. Cook, Y. Zheng, A. Egan, P. A. Yushkevich, and J. C. Gee.  N4ITK: improved N3 bias correction.  _IEEE transactions on medical imaging_ , 29(6):1310–1320, 2010.
  * Tustison et al. (2014) N. J. Tustison, P. A. Cook, A. Klein, G. Song, S. R. Das, J. T. Duda, B. M. Kandel, N. van Strien, J. R. Stone, J. C. Gee, et al.  Large-scale evaluation of ANTs and FreeSurfer cortical thickness measurements.  _Neuroimage_ , 99:166–179, 2014.
  * Tustison et al. (2019) N. J. Tustison, A. J. Holbrook, B. B. Avants, J. M. Roberts, P. A. Cook, Z. M. Reagh, J. T. Duda, J. R. Stone, D. L. Gillen, M. A. Yassa, et al.  Longitudinal mapping of cortical thickness measurements: An Alzheimer’s Disease Neuroimaging Initiative-based evaluation study.  _Journal of Alzheimer’s Disease_ , 71(1):165–183, 2019.
  * Weston et al. (2016) P. S. Weston, J. M. Nicholas, M. Lehmann, N. S. Ryan, Y. Liang, K. Macpherson, M. Modat, M. N. Rossor, J. M. Schott, S. Ourselin, et al.  Presymptomatic cortical thinning in familial Alzheimer disease: A longitudinal MRI study.  _Neurology_ , 87(19):2050–2057, 2016.
  * Wu et al. (2022) Y. Wu, S. Kundu, J. S. Stevens, N. Fani, and A. Srivastava.  Elastic shape analysis of brain structures for predictive modeling of PTSD.  _Frontiers in Neuroscience_ , 16:954055, 2022.  10.3389/fnins.2022.954055.
  * Zhou et al. (2013) H. Zhou, L. Li, and H. Zhu.  Tensor regression with applications in neuroimaging data analysis.  _Journal of the American Statistical Association_ , 108(502):540–552, 2013.