{'fraction_non_alphanumeric': 0.027586409076251064, 'fraction_numerical': 0.01187300284365747, 'mean_word_length': 4.673985362608117, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian tensor regression using the Tucker decomposition for sparse spatial modeling

Modeling with multidimensional arrays, or tensors, often presents a problem due to high dimensionality. In addition, these structures typically exhibit inherent sparsity, requiring the use of regularization methods to properly characterize an association between a tensor covariate and a scalar response. We propose a Bayesian method to efficiently model a scalar response with a tensor covariate using the Tucker tensor decomposition in order to retain the spatial relationship within a tensor coefficient, while reducing the number of parameters varying within the model and applying regularization methods. Simulated data are analyzed to compare the model to recently proposed methods. A neuroimaging analysis using data from the Alzheimer's Data Neuroimaging Initiative is included to illustrate the benefits of the model structure in making inference.


Introduction
Image analysis is an important application area with the development of computer memory that allows for storing large datasets on local computing machines. Machine learning applications of image analysis are now present in many fields of research, such as medical imaging, character translation, and self-driving cars.
In many of these cases, a fast modeling structure that provides inferential and/or prediction capabilities is all that is required. Such models have been shown to be very effective in settings with billions or trillions of observed data points. However, in some cases, low signal-to-noise ratios can present a challenge to such models, and additional model constraints and assumptions are needed to improve model power.
The field of neuroimaging is particularly rich in high-dimensional data with low signal-to-noise ratios.
Often, the goal in the analysis of these data is to make popoulation-level inferences about the function of the brain and the pathology of specific neurological disorders using high-resolution structural brain images taken from living subjects. Due to the clinical importance of inference from these models, assumptions and constraints must be carefully applied.
Several methods are already in use for these types of data within the neuroimaging community. One of the most commonly-used methods is referred to as the general linear model (GLM), which is not to be confused with the generalized linear model that is commonly used in statistics. The GLM performs a massive univariate analysis in which the response is regressed independently on each voxel within a tensor covariate, in addition to any additional vector covariates [Friston et al., 1995, Penny et al., 2011. This model has the advantage of being relatively computationally inexpensive. However, the GLM also assumes that the associations between different voxels in the tensor and the response are independent and no assumptions are made about the underlying sparsity of the tensor coefficient, which negatively impacts model power. In practice, different multiple testing corrections are used to preserve spatial relationships among neighboring voxels, though work by Eklund et al. [2016] suggests that these inflate the false discovery rate. Methods that control the false discovery rate are appealing for their interpretability [Benjamini andHochberg, 1995, Lindquist andMejia, 2015], but they fail to take spatial relationships within the tensor coefficient into account.
A different class of approaches takes advantage of the tensor structure of the data by decomposing the tensor covariates using one of two tensor decompositions and imposing regularization constraints on the tensor factors. Work done by Zhou et al. [2013] uses the parallel factorization/canonical polyadic (PARAFAC/CP) tensor decomposition in a classical tensor regression approach, which assumes that dimension margins are independent components for construction of the tensor coefficients. This was expanded in work by Li et al. [2018] in the use of the more flexible Tucker decomposition. Guhaniyogi et al. [2017] created a novel Bayesian prior structure on the PARAFAC/CP tensor decomposition elements, which improved on the uncertainty quantification from the model inference.
We propose a Bayesian tensor regression model using the Tucker decomposition (BTRT) that assumes sparsity and spatial dependence within a tensor-valued coefficient. This is accomplished through the use of the Tucker tensor decomposition [Tucker, 1966], which can be thought of as an analog to an independent components analysis in which the components are dimension margins. We show the efficiency and accuracy of the model over a number of competitors through simulation studies. The BTRT model is then used to perform a neuroimaging analysis of data from the Alzheimer's Disease Neuroimaging Initiative. The BTRT model and tensor regression competitors are implemented in the R software package bayestensorreg [Spencer, 2022].
We formally define terminology and use of tensor notation, and explain the BTRT and competitor models in section 2. Inferential accuracy is examined through the use of simulated datasets in section 3. An analysis of data from the Alzheimer's Disease Neuroimaging Initiative in section 4 illustrates improvements to inferential power using real data. We end with a discussion of our findings and directions for future research in section 5.

Methodology
This section begins with the introduction of tensor notation and the formal definition of two tensor decompositions. Next, the BTRT model and its prior structure is described. Model identifiability is then discussed, followed by a guide on how to choose the model rank. Finally, competitor models are outlined, which validate our proposed method as a useful and reliable tool in sparse tensor regression scenarios.

Tensor terminology and notation
A tensor of order D is a multi-dimensional array data structure A ∈ R p1,...,p D . To illustrate, a vector is a tensor of order 1, a matrix is a tensor of order 2, a cubic array is a tensor of order 3, and so forth.
The vectorization of a tensor A ∈ R p1×···×p D of order ≥ 2 results in a tensor of order 1 of length D j=1 p j , that is vec(B) ∈ R p1...p D . The inner product of two tensors A and A is the crossproduct of the vectorized elements of the tensors, that is A, A = vec(A) vec(A ).
The kth-mode matricization of a tensor, represented as A (k) is a matrix representation of a tensor of order 2 or higher such that the kth index becomes the first index and all other tensor indices are combined in order into a second index. That is, Tensors can be approximated via different tensor decompositions. One such decomposition is called the canonical decomposition/parallel factorization, also known as CANDECOMP/PARAFAC, or CP [Tucker, 1966]. This decomposition represents the tensor B as in which β j,r ∈ R pj is one of R principal components for the jth dimension of B. Here, R is known as the rank of the CP decomposition. The • operator denotes the outer product. At each value of r ∈ {1, . . . , R}, the series of outer products β 1,r • · · · • β D,r results in a D-dimensional tensor summand B r ∈ R p1,...,p D .
As the value of R increases, the number of summands in the decomposition increases and the resolution improves. All tensors can be represented exactly by the CP decomposition given sufficiently high rank R.
In practice, low-rank decomposition approximations have been found to be adequate to estimate sparse tensor coefficients with spatial dependence. Using the CP decomposition reduces the parameter space for estimating B from D j=1 p j to R D j=1 p j , while incorporating this sparse spatial structure. However, it is possible that not all dimension margins in a tensor decomposition require all R components in order to faithfully represent structure within a tensor. In such cases, the CP decomposition may be extended to the Tucker decomposition, which can be written as where (g 1,...,1 , g 2,...,1 , . . . , g R1,...,R D ) = G ∈ R R1,...,R D is the core tensor composed of elements, which assigns weights of importance to each of the D j=1 R j tensor summands that compose the tensor B [Tucker, 1966]. This representation is more flexible than the CP decomposition, as it allows for different decomposition ranks for different tensor dimension margins, and all ranks of different margin factors interact with every rank in other dimension margins. As an example, consider an order 2 tensor. A rank 2 CP decomposition uses two tensor summands to approximate the tensor, while a rank 2, 2 Tucker decomposition uses a weighted sum of four tensor summands at the cost of only four additional model parameters from the core tensor. The parameter space for estimating B using the Tucker decomposition is With these basic conventions, a linear model can be built using the Tucker decomposition that is more flexible with only a modest increase in the number of parameters in the parameter space.

Tensor Regression Model
Assuming a scalar response y i , vector covariates η i ∈ R q , and tensor covariates X i ∈ R p1×···×p D for subjects i = 1, . . . , n, the observed linear model can be represented as in which B ∈ R p1×···×p D is a tensor coefficient, γ ∈ R q is a vector coefficient, and i is an error term, which follows any distribution centered at zero. Note that, since we are focusing on the methods to address tensor regression, we are restricting our attention to scenarios in which q is relatively small. However, the model framework here does allow for a seamless extension to high-dimensional modeling of η i .
In the potentially large predictor space within the tensor-valued covariate, it is reasonable in certain applications to assume spatial dependence in the association between X i and y i . In order to impose this structure in the model while simultaneously reducing the parameter space, we use the Tucker tensor decomposition, as outlined in (2). We use Bayesian methods to define assumptions about the sparsity of the tensor coefficient using state-of-the-art regularization priors. As the other terms in the model are not the focus of our proposal, typical Bayesian modeling techniques are applied to the γ and i terms, as outlined below.

Bayesian Tensor Regression Tucker prior structure
We focus our model on applications in which a degree of sparsity can be reasonably assumed in the tensorvalued coefficient. Classical regularization methods used to meet such a goal include various penalized regression algorithms like the LASSO [Tibshirani, 1996] or RIDGE regression [Hoerl and Kennard, 1970].
However, these methods lack the ability to provide a measure of uncertainty quantification on the estimates for the parameters, which is important for determining the practical significance of the estimates for such parameters. Shifting to a Bayesian modeling structure can serve to fill this gap, leading to improved inference. Li et al. [2018] proposes a classical model for tensor regression using the Tucker decomposition in which a penalty is applied to either the core tensor G or both the core tensor G and each dimension component β j,rj in the decomposition. We propose shrinkage priors on both the core tensor and all dimension components in order to more strongly induce parameter regularization.
Following the previous work done by Guhaniyogi et al. [2017], an adapted generalized double-Pareto prior is applied to the dimension components within the Tucker tensor decomposition. That is, where W j,rj is a diagonal matrix with elements ω j,rj , for = 1, . . . , p j . Integrating over the element-specific scale parameters reduces the prior on β j,rj , to a double-exponential (Laplace) distribution centered at 0 with a scale parameter of λj,r j √ τ , which has heavier tails than a Gaussian distribution while also allocating higher densities around zero. This prior is a Bayesian analog to the adaptive LASSO, inheriting its desirable oracle properties of the maximum a posteriori estimator [Armagan et al., 2013].
In order to adequately select the proper rank for each dimension and to reduce noise in the tensor coefficient estimates, the generalized double-Pareto prior is also imposed on the elements of the core tensor G: This combination ensures that only tensor summands within the Tucker decomposition that explain additional variance are included in the estimate of B.
For the purposes of our analyses, a multivariate normal prior with mean µ γ and covariance Σ γ is placed on the elements of γ. This is done to maintain conjugacy while maintaining control over the expected effects of the vector-valued coefficients. The errors i s in this model are assumed to be independent identically distributed following a normal distribution with a mean of zero, and a variance of σ 2 y . An inverse gamma prior is placed on σ 2 y with hyperparameters a σ and b σ in order to maintain conjugacy.

Identifiability
The decomposition of the tensor coefficient affects the identifiability of the parameters in the model. That is, the value of any voxel within the tensor coefficient is estimated aŝ where v is the voxel location within the tensor (v = (v 1 , . . . , v D )). Any two of these summands can be multiplied by c and 1 c , respectively, and the estimate for that voxel within the tensor coefficientb v would remain unchanged. Indeed, the values of g r1,...,r D and β j,rj ,vj are not identifiable, butb v remains identifiable because of the shrinkage priors on the tensor decomposition components. In testing the BTRT model in very low signal-to-noise settings, models in which at least one margin rank R j is equal to one can result in divergent MCMC chains because core tensor values multiply with only the rank 1 tensor decomposition factor. In our experiments, increasing all R j to be ≥ 2 has alleviated these issues.
Something to consider when using the model in Equation (3) for a data analysis is that an identifiability problem exists between B and γ if there are voxels within X i , ∀i = 1, . . . , n that are collinear with any of the values in η i , ∀i = 1, . . . , n. This is a problem that may be undetected in two-step frequentist models with regularization, as either elements in B or γ may simply be assigned values of zero.

Selection of Rank
A key consideration in the use of the BTRT model structure is the selection of each dimension rank. Increases in a margin's rank can be made with an attempt to increase the spatial resolution on the inference of the tensor coefficient. For larger tensor covariates, this requires higher ranks if the nonzero coefficients are not arranged hypercubically. Choice of unequal ranks may be prudent if some of the tensor dimensions are much larger or smaller than others. Consider an example in which the tensor covariate has dimensions 100 × 100 × 4. The ranks for the first two dimensions may need to be considerably larger than the rank for the third dimension, as there are only four possible margin locations in the final dimension. One such realistic application would be to combine magnetic resonance images that use different sequences (e.g. T1 weighted, T2 weighted, effective T2, etc.) into a single tensor, using a low rank on the dimension that represents different sequences. This is a clear advantage over the CP/PARAFAC decomposition methods, as each dimension is not forced to have the same number of ranks as all of the others.
Once a set of reasonable ranks is decided for each tensor dimension, the different combinations of the ranks should be used to fit models and compare results using some model selection criterion like the deviance information criterion (DIC) [Gelman et al., 1995].

Competitor Models
The effectiveness of the proposed Bayesian sparse tensor regression models is shown by making direct comparisons to models commonly used in the field of neuroscience. The first is the general linear model (GLM), in which the response (y 1 , . . . , y n ) = y is regressed on each voxel v within the tensor covariate X ∈ R p1×...×p D ×n independently of the vector-valued coefficients. The linear model for the GLM is written in equation (4).
Note that the same responses (y 1 , . . . , y n ) = y are used to fit separate models for each voxel within the tensor covariate.
This model is an industry standard for its ease of implementation and high computational efficiency. However, this model often suffers from a high false discovery rate in multiple-subject analyses [Eklund et al., 2016].
We adjust for the high false discovery rate in our implementation by using the Benjamini-Hochberg multiple testing correction [Benjamini and Hochberg, 1995], which fixes the false discovery rate (FDR). In the following trials, the FDR is fixed at 0.05, in agreement with standard practice in neuroimaging.
Direct tensor regression competitors are also used to validate the performance of the proposed model structure. Specifically, the frequentist tensor regression using the PARAFAC/CP tensor decomposition (FTR CP) [Zhou et al., 2013] and the Tucker tensor decomposition (FTRT) [Li et al., 2018] are used to compare point estimates to classical methods. These methods apply can the LASSO to the margin factors in their respective tensor decompositions in Equations (1) and (2). The FTRT model also allows for the application of the LASSO to the elements of the core tensor. In our tests, the FTR CP model performs best with the LASSO applied to the margin factors and the FTRT model performs best when the LASSO is only applied to the elements within the core tensor, so these are the configurations used in our comparisons.
We also compare to the Bayesian tensor regression using the PARAFAC/CP tensor decomposition (BTR CP) [Guhaniyogi et al., 2017] provides a more direct comparison in terms of point estimates and uncertainty quantification.

Simulated Data Analysis
On order to demonstrate the efficacy of the BTRT model, data were simulated from the linear model in (3) under the following conditions. B ∈ R 50×50 , where nonzero-valued elements take the value of 1 in the middle of the regions, decreasing to lower nonzero values using the specifyregion function within the neuRosim package in R [Welvaert et al., 2011]. , where the order of the coefficient tensor is set to D = 2. The prior on the observation error variance is set to be relatively noninformative in the context of the simulation, with a mean of 10 and a variance of 100. The priors for λ and ϕ are set to have modes between 1.5 and 2, approaching 2 as the tensor dimension increases, which places the prior expected value for v r1,...,r D and the elements within W j,rj to be between 2 3 and 1 2 . The priors for τ and z have a mean of 1 when all R 1 , . . . , R D = 1, increasing sublinearly with both the minimum value of R j for j = 1, . . . , D, and the tensor dimension D. The prior variance for τ and z increases linearly in rank when D = 2, and superlinearly in rank when D > 2. This prior specification allows for a slightly higher prior variance for β j,rj and g r1,...,r D as rank and dimension increase in order to allow for moderately faster exploration of the parameter space as the tensor dimension and rank increase.
The prior for γ is set to be relatively noninformative.
For all 16 models with R 1 , R 2 ∈ {1, 2, 3, 4}, 11,000 Markov Chain Monte Carlo (MCMC) simulations from the posterior distribution were drawn. After discarding the first 1000 draws, the remaining draws were used to estimate B and γ.
Continuous distribution priors do not produce parameter estimates equal to zero. In order to perform variable selection, the sequential 2-means variable selection method proposed by Li and Pati [2017]  Point estimates for the tensor coefficient from the Bayesian tensor regression (BTR) models, frequentist tensor regression (FTR) models, and the GLM are shown in Figure 1. The ranks of the models shown were selected using the DIC for the Bayesian models or the Akaike's Information Criterion (AIC) for the frequentist models. The estimate shown for the GLM tensor coefficient assigns values of zero to voxels within the tensor that are not statistically significant after using the Benjamini-Hochberg multiple testing correction with a Type I error rate of 0.05. Included in this figure are the values for the root mean squared error, defined as whereB v is the point estimate for tensor coefficient element v obtained from a model, B v is the true value taken by the coefficient tensor element v, and V is the total number of elements in the coefficient tensor. These results show that all models except for the BTR CP model choose higher model ranks in order to capture the shapes and magnitudes of the three distinct nonzero regions in the tensor coefficient. Higher ranks are necessary in order to capture multiple nonzero regions, but the BTR CP model shows increased posterior variance as model rank increases, which biases the rank selection for that model toward lower values. The BTRT estimate is both visually and numerically most similar to the true values in this simulated data analysis. Notably, all of the tensor regression models perform significantly better than the two-step GLM model, which is often used in neuroimaging studies. Appendix A contains Table 2, which gives the RMSE values for all ranks of all competing models, and Figure 5, which shows the point estimates for all ranks of the BTRT model. The posterior densities for {γ 1 , γ 2 , γ 3 } under different models with selected ranks can be seen in Figure   2. These results show that the posterior densities for each γ j are centered at or around the true values. The

Model Convergence
In Bayesian modeling settings, it is important to ensure that the MCMC converges around an area within the neighborhood of the global posterior mode. Given that the posterior inference matches the inference shown within the frequentist models, the model does converge to a global mean in the simulated data settings when the assumption of normal error is satisfied, and the prior distributions are specified to be relatively uninformative. The log-likelihood values using the posterior draws, shown in Figure 6 in Appendix A, also show rapid convergence to a mode and posterior stability. The median effective sample sizes, as defined in Gong and Flegal [2014], across all voxels in B range from about 5400 to about 8100 out of 10000 draws from the posterior distribution, indicating sufficient numbers of effective samples from the posterior distribution were taken to make reliable inference.

Hyperparameter Sensitivity
Modeling in such a high-dimensional space with several hierarchical levels in the prior structure requires careful selection of hyperparameter values. However, the model is still expected to be somewhat robust to differences in prior specification. In order to test this expectation, a large grid of hyperparameter values for a σ , b σ , a λ , a τ , a z , and a ϕ was created. First, the values for the hyperparameters used to model the simulated data are taken and multiplied by 0.1, 1, or 10. One hundred configurations were randomly selected from the 729 possible configurations, and the model with ranks R 1 = R 2 = 3 was run with the same simulated data as above for 11,000 iterations with each combination. In each case, 1,000 iterations in the MCMC were discarded as a burn-in.
Of the 100 configurations, 93 converged to have an RMSE for the tensor coefficient B between 0.03279 and 0.03336. The remaining 7 configurations all had RMSEs within 10 −8 of 0.1336523, which is the RMSE that results from predicting B with a tensor of all zeros. The priors for these seven configurations all had the smallest values for the prior variance for the elements of B (approximately 2 × 10 −3 ) and/or the highest prior means (approximately 10 10 ) and variances (over 10 20 ) for σ 2 y . Indeed, individual inspection of these models with larger RMSE values show that the estimates for γ are within the same range as the estimates from the models that had a lower RMSE for the tensor coefficient, but they all estimate that the tensor coefficient is zero. This highlights the need to carefully balance the amount of shrinkage that can be applied to elements within B before no image coefficients are recovered. In addition, the prior specification for σ 2 y must be reasonable in order to achieve interpretable results.

Analysis of ADNI data
Structural differences in the brain can be measured using tensor-based morphometry (TBM), which provides a measure for the amount of deformation that needs to be applied to each voxel in order to match a given template [Brun et al., 2010]. Larger values indicate that the neighborhood around a given voxel is larger relative to the standard to which subject structural scans were registered, and smaller values indicate atrophy. These measures have been used to deduce phenotypes for certain groups, such as those diagnosed with HIV/AIDS [Lepore et al., 2006]. In our analysis, we apply the BTRT method and competing methods to measure associations between a subject's performance on the mini mental state exam (MMSE) and TBM data in the "TBM Jacobian Maps MDT-SC" collection from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database , which includes data for 817 subjects (342 female) from the first phase of the ADNI (ADNI-1). The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of the ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). All data were accessed April 21, 2022. Of these subjects, 188 were diagnosed with early Alzheimer's, 400 subjects exhibited mild cognitive impairment, and 229 were healthy controls. Phenotype data were also included for each subject in the analysis, matched to patient ID using the ADNIMERGE package in R [the ADNI team, 2022].
The MMSE is a diagnostic tool used to classify adults based on levels of cognitive impairment [Pangman et al., 2000]. The exam itself poses a series of questions to test the subject's ability to perform everyday tasks, i.e. reading an analog clock or reproducing a drawn shape. The maximum score that one can achieve on the test is 30, and the scores take integer values. These scores can then be used in conjunction with other information to diagnose an individual. MMSE values for the subjects in this analysis ranged from 18 to 30.
We include age, gender, years of education, and the number of Apolipoprotein E4 (APOE4) genetic variants as additional covariates for each subject. Age and gender are included to account for any confounding of those factors on the MMSE score. The number of years of education show a significant association with the MMSE scores in an exploratory data analysis. APOE4 has been identified as a genetic risk factor for Alzheimer's disease [Strittmatter and Roses, 1996], and this variable can take the values 0, 1, or 2. We treat the APOE4 covariate as continuous, as exploratory data analysis suggests a linear association with the MMSE.
We compare the performance of the different models in this analysis with the root mean squared predictive error (RMSPE) and the Pearson correlation between the predictions and the actual values in the dataset.
The Pearson correlation is included because it was used as a performance metric in the DREAM challenge held by ADNI [Zhu et al., 2016]. The BTRT model is run for 1,500 iterations, after which 500 iterations are discarded as a burn-in. The BTR CP model was successfully tested using simulated data under several different scenarios. However, MCMC chains for the TBM data fail to maintain numerical stability despite starting from several different initial conditions, and the BTR CP model was unable to produce more than about 30 samples from the posterior distribution for any model rank. Thus, the BTR CP results are excluded from comparison here. For the BTRT model, the point predictions for the response are calculated as the medians of the posterior predictive distributions, which are found for each subject aŝ for each sample s from the posterior distribution. The frequentist tensor regression models are run until the relative log-likelihood change between steps is less than 0.1%, and they predict the response value for each subject using their point estimates for B and γ.
Final point estimates for the tensor coefficient are found in the Bayesian models by using the sequential 2-means post-hoc variable selection method [Li and Pati, 2017], which is outlined in Appendix B. The frequentist tensor regression models are not corrected for multiple testing, as they use the LASSO to select the coefficients that are significantly different from zero. The GLM model estimate is found by applying the Benjamini-Hochsberg multiple testing correction to the image coefficient p-values and setting voxels that are not significant to have values of zero. The deviance information criterion (DIC) is used to select which BTRT model should be used to fit the data, and the Akaike information criterion (AIC) is used to select the rank of the FTR models. We choose to use the AIC for the frequentist model selection because it is analogous to the DIC for the Bayesian models. The plots for the final estimates can be seen in Figure 3.
The estimates shown here have had their values thresholded to be between ±1 × 10 −5 in order to show how models perform in terms of indicating areas of nonzero values, as they produce estimates that differ in orders of magnitude. We used the automated anatomical labelling atlas [Tzourio-Mazoyer et al., 2002] to highlight the amygdala in the background image, as it has been implicated in previous studies of Alzheimer's disease [Poulin et al., 2011]. We mask the image coefficient estimates to be within the brain to match our inferential goal of finding associations within brain structures, rather than measures of deformation outside the brain that would correspond to associations with total brain volume.
Estimates for the image coefficient greatly vary among the different models, which may suggest a very low signal-to-noise ratio within these data to conclude that there are nonzero coefficient values. Nonetheless, the point estimate from the BTRT model clearly highlights positive associations between the relative size of the amygdala region and the MMSE scores. The GLM finds no locations that are significantly different from zero after performing a multiple testing correction to limit the false discovery rate and using a Type while the FTR CP model may be too sensitive.
Posterior distributions for the non-tensor coefficients (γ) are shown in Figure 4, along with the point estimates for the parameters from the GLM and the FTR models. These densities show an agreement between the competitor models and the BTRT model in the estimation of these coefficients in the BTRT model.   Table 1. We can see from this table that the FTRT, GLM, and no-image competitor models all perform very similarly, and the selected FTR CP and BTRT models perform dramatically better in both measures. However, the BTRT model outperforms the FTR CP model while also providing more interpretable estimates for the image coefficient.

Discussion
We propose a Bayesian tensor regression model using the Tucker tensor decomposition that accounts for spatial dependence through the structure of the prior without explicitly defining a Gaussian process or an  [Spencer, 2022]. We compare the performance of the Bayesian tensor regression model using the Tucker decomposition to other recent tensor regression models using simulated data to show improved inference and model fit. We also apply the BTRT model and competitors in the analysis of neuroimaging data from the ADNI project.
The BTRT model shows improvement over competing models within simulation studies. The BTRT models outperform the frequentist Tucker tensor regression models by more closely matching the true tensor coefficients visually and numerically via the root mean squared error metric. The selected BTRT model dramatically outperforms the BTR CP model in terms of image coefficient estimates due to lower variance within the MCMC chains, which remains relatively consistent among different model ranks. The BTRT model also produces more specific posterior distributions than the BTR CP model for the non-image-valued parameters centered around the true generated values, indicating improved inferential performance. All of the tensor regression models show a marked improvement over a general linear model, which is often used in neuroimaging literature.
Analysis of the TBM data from the ADNI illustrates the value of the BTRT model in providing meaningful inference while performing well on model reliability metrics, suggesting that the model assumptions are able to detect sparse, spatially-dependent nonzero image coefficients better than FTR models. While the BTRT model performs impressively in terms of model fit on these two dimensional cases, additional work needs to be done in order to make this model a feasible option for fitting whole-brain volumes. Computational costs can be driven down using expectation-maximization or variational inference methods, and implementation in a more efficient programming language is also expected to improve efficiency. A Other performance metrics for simulated data analysis   B Sequential 2-means algorithm 11 Set elements inθ with then z smallest absolute values to 0 Figure 6: The log-likelihoods for the Bayesian models applied to the simulated data. Labels at the end of the lines give the model rank. These plots suggest rapid convergence of the MCMC in all models, though the variance in the log-likelihood from iteration-to-iteration is much lower in the Tucker decomposition models than in the CP/PARAFAC decomposition models.