Bayesian Covariance Estimation for Multi-group Matrix-variate Data

Matrix-variate datasets consisting of a sample of $n$ matrices $Y_{1},...,Y_{n}$ each with dimensionality $p_{1}\times p_{2}$ are increasingly common in modern applications. Examples of such datasets include repeated measurements of a multivariate response, two-dimensional images, and spatio-temporal observations. Often, a matrix-variate dataset may be partitioned into several distinct groups or subpopulations for which group-level inference is of particular interest. For example, a multi-group dataset may be subdivided by socio-economic population or geographic region.

In analyzing multi-group matrix-variate data of this type, describing heterogeneity across groups is often of particular interest. For instance, in remote-sensing studies, scientists may be interested in understanding variation across classes of land cover from repeated measurements of spectral information, as in Johnson et al. (2012). The information collected at each site may be represented as a $p_{1}\times p_{2}$ matrix where the rows represent $p_{1}$ wavelengths and the columns represent $p_{2}$ dates when the images were taken. Accurate multi-group covariance estimation is necessary for such a task. Moreover, in many such applications, the population covariance of matrix-variate data may be near separable in that the covariances across the $p_{2}$ columns are near each other or the covariances across the $p_{1}$ rows are near each other. Incorporation of this structural information in an estimation procedure can improve estimation accuracy. More generally, accurate covariance estimation of multi-group matrix-variate data is a pertinent task in many statistical methodologies including classification, principal component analysis, and multivariate regression analysis, among others. For example, classification of a new observation based on a labeled training dataset with quadratic discriminant analysis (QDA) requires group-level estimates of population means and covariance matrices. As a result, adequate performance of the classification relies on, among other things, accurate group-level covariance estimates.

One approach to analyzing multi-group matrix-variate data is to vectorize each data matrix and utilize methods designed for generic multivariate data separately for each group. In this way, direct covariance estimates which only make use of within-group samples may be obtained from matrix-variate data by vectorizing each observation and computing the standard sample covariance matrices. While a group’s sample covariance may be unbiased, the estimate may have large variance unless the group-specific sample size $n_{j}$ is appreciably larger than the dimension $p=p_{1}p_{2}$. This is a limiting requirement as modern datasets often consist of many features, that is, often $p\approx n_{j}$, or even $p\gg n_{j}$. As a result, more accurate covariance estimates may be obtained via indirect or model-based methods which incorporate auxiliary information. Such methods may introduce bias, but can correspond to estimates with lower variance than unbiased methods.

To improve the accuracy of covariance estimates for multi-group data, researchers may estimate each group’s covariance with the pooled sample covariance matrix. This implicitly imposes an assumption of homogeneity of covariances across the populations and greatly reduces the number of unknown parameters to be estimated. Such an estimator may be biased, but can have lower error than the population-specific sample covariance. For example, linear discriminant analysis (LDA), which assumes homoscedasticity across groups, has been shown to outperform QDA when sample sizes are small, even in the presence of heterogeneous population covariances (see, for example, Marks and Dunn (1974)). A more robust approach estimates each group’s covariance as a weighted sum of a pooled estimate and the group-specific sample covariance. Such an approach is often referred to as partial-pooling, or, in the Bayesian framework, hierarchical modeling. For a nice introduction to partial pooling and an empirical Bayesian implementation, see Greene and Rayens (1989). More work on the topic is found in Friedman (1989); Rayens and Greene (1991); Brown et al. (1999). Relatedly, there has been work which assumes pooled elements of common covariance decompositions (e.g. pooled eigenvectors across groups (Flury, 1987)) and proposals to shrink elements of decompositions to pooled values (Daniels, 2006; Hoff, 2009b).

Alternatively, as opposed to pooling information across groups, accuracy may be improved by imposing structural assumptions on the covariances separately for each group. Some common structural assumptions include diagonality (Daniels and Kass, 1999), bandability (Wu and Pourahmadi, 2009), and sparsity (Friedman et al., 2008), among others. For matrix-variate data, a separable or Kronecker structure covariance assumption (Dawid, 1981) may be more appropriate. A Kronecker structured covariance represents each $p\times p$ population covariance as the Kronecker product of two smaller covariance matrices of dimension $p_{1}\times p_{1}$ and $p_{2}\times p_{2}$ which respectively represent the across-row and across-column covariances. Again, while a separable covariance estimator may be biased, it may have smaller error than an unstructured covariance estimator. In practice, however, the population covariance may not be well represented by a Kronecker structure. To allow for robustness to misspecification, a researcher may proceed in a Bayesian manner and adaptively shrink to the Kronecker structure as in Hoff et al. (2022). Such an estimator can be consistent, but may have stability issues similar to that of an unstructured covariance matrix if the population covariance is not well represented by a Kronecker structure. In this instance, more accurate multi-group covariance estimates may be obtained by pooling across groups rather than shrinking within each group to a separable structure.

More generally, in covariance estimation based on matrix-variate data from multiple populations, it is rarely obvious whether shrinking each unstructured covariance separately towards a Kronecker structure or shrinking all unstructured covariances towards an unstructured pooled covariance will result in more accurate estimates. This is particularly difficult in the presence of small within-group sample sizes as popular classical statistical tests for both homogeneity of covariances (Box, 1949) and accuracy of a Kronecker structure assumption (Lu and Zimmerman, 2005) rely on approximations that require large sample sizes to achieve the desired precision.

To account for this uncertainty, we propose a hierarchical model that adaptively allows for both types of shrinkage. Specifically, in this article, we provide a model-based multi-group covariance estimation method for matrix-variate data that improves the overall accuracy of direct covariance estimates. We propose a hierarchical model for unstructured group-level covariances that adaptively shrinks each estimate either within-population towards a separable Kronecker structure, across-populations towards a shared pooled covariance, or towards a weighted additive combination of the two. The model features flexibility in the amount of each type of shrinkage. Furthermore, the proposed model has a latent-variable representation that results in straightforward Bayesian inference via a Metropolis-Hastings algorithm. The proposed model provides robustness to mis-specification of structural assumptions and improved stability if assumptions are wrong while maintaining coherence and interpretability.

The article proceeds as follows. In Section 2 we motivate our method and detail the proposed hierarchical model. We describe a Bayesian estimation algorithm in Section 3 and demonstrate properties of the proposed method via a simulation study in Section 4. The flexibility of the proposed method is shown in two examples in Section 5. In the first example, we demonstrate the usefulness of inference under the proposed model in speech recognition. In the second example, we perform inference on a chemical exposure data set where understanding heterogeneity across socio-demographic groups is of key interest. We conclude with a discussion in Section 6.

In this section we introduce the Shrinkage Within and Across Groups (SWAG) covariance model, a hierarchical model developed for simultaneous covariance estimation for multi-group, matrix-variate data. We are particularly motivated by improving the overall accuracy of group-specific estimates of population covariances when the true covariance structures are unknown and group-specific sample sizes are small relative to the number of features. The proposed model adaptively allows for flexible shrinkage either across groups, within a group to a Kronecker structure, or an additive combination of the two. The SWAG model is constructed from semi-conjugate priors to allow for straightforward Bayesian estimation and interpretable parameters. The section proceeds by introducing motivation in Sections 2.1 and 2.2, presenting the proposed hierarchical covariance model in Section 2.3, and elaborating on parameter interpretation in Section 2.4.

We demonstrate the performance of the proposed SWAG model by comparing the accuracy of various covariance estimators obtained under four different population covariance regimes. In particular, each regime features either homogeneous (Ho) or heterogeneous (He) covariances across groups, and the covariances in each regime are either Kronecker structured (K) or not Kronecker structured (N). Our goal in this section is to generally explore results when the true group-specific covariance matrices have off-diagonal values relatively far from zero within a given group. Specifically, in unstructured regimes, each group’s true covariance is an exchangeable correlation matrix of dimension $p\times p$ with a fixed correlation randomly generated between 0.35 and 0.9. In Kronecker structured regimes, each group’s true covariance is the Kronecker product of a $p_{2}\times p_{2}$ exchangeable correlation matrix and a $p_{1}\times p_{1}$ exchangeable correlation matrix. For a given regime and parameter size combination, the true covariance matrices do not vary within the simulation. Details of the within-group true covariances $\{\Sigma_{1},...,\Sigma_{J}\}$ for each regime are contained in Contrast Table 1. While we do not necessarily expect the truth in real-world scenarios to be one of these extreme cases, this study provides insight into the behavior of the SWAG model.

In a simulation study under each regime, we compare the loss of various covariance estimators under a range of dimensionality and number of groups. Specifically, we consider $p_{1}\in\{2,4,8\}$, $p_{2}=3$, and $J\in\{4,10\}$. As we are particularly motivated by the small sample size case, we consider within-group sample sizes $n_{1}=...=n_{J}=p+1$. For each regime, we take the MLEs of the covariances under the simplest correctly specified model as the oracle estimator. In total, this includes the standard sample MLE for each group’s covariance $\boldsymbol{S}=\{{S}_{1},...,{S}_{J}\}$, the pooled sample MLE $\hat{S}_{p}$, the separable MLE for each group $\hat{\boldsymbol{K}}=\{\hat{K}_{1},...,\hat{K}_{J}\}$, and the pooled separable MLE $\hat{K}_{p}$. These oracle estimators for each group’s covariance under the four regimes considered are specified in Contrast Table 2. We note that the most accurate estimator of the He,U regime may vary depending on the problem dimension and the number of the groups due to high variance of the sample covariance estimator when sample size is small relative to number of features, as in this study.

The simulation study proceeds as follows. For each permutation of regime, dimension, and number of groups, we sample 50 data sets each from a mean-zero matrix normal distribution of the appropriate dimension with population covariance as given in Table 1. For each data set, we compute the reference covariance estimates $\boldsymbol{S}$, $\hat{S}_{p}$, $\hat{\boldsymbol{K}}$, and $\hat{K}_{p}$. Additionally, we run the proposed Metropolis-Hastings sampler for the SWAG model for 28,000 iterations removing the first 3,000 iterations as a burn-in period and saving every 10th iteration as a thinning mechanism. From the resulting 2,500 Monte Carlo samples, we obtain the Bayes estimate under an invariant loss, Stein’s loss, of each group’s covariance, $\hat{\boldsymbol{\Sigma}}=\{\hat{\Sigma}_{1},...,\hat{\Sigma}_{J}\}$ where $\hat{\Sigma}_{j}=E[\Sigma_{j}^{-1}|Y_{j}]^{-1}$. Then, we compute Stein’s loss averaged across the populations $\bar{L}(\boldsymbol{\Sigma},\hat{\boldsymbol{\Sigma}})=\frac{1}{J}\sum_{j=1}^{J}L_{S}(\Sigma_{j},\hat{\Sigma}_{j})$ where $L_{S}(\Sigma_{j},\hat{\Sigma}_{j})=tr\left({\Sigma}_{j}^{-1}\hat{\Sigma}_{j}\right)-\log\left|{\Sigma}_{j}^{-1}\hat{\Sigma}_{j}\right|-p_{1}p_{2}$. We report the average of the 50 $\bar{L}$ values to approximate frequentist risk for each scenario considered. In our analysis of results, we refer to $\bar{L}$ as the loss and do not discuss group-specific Stein losses.

In general, we expect inference with the ‘oracle’ estimator to outperform that of the SWAG model in each regime considered. However, as knowledge of true structural behavior is rare in practice, the oracle estimator of one regime obtained from a correctly specified model may perform arbitrarily poorly in a different regime. In contrast, due to the flexibility of the proposed SWAG model, we expect the SWAG estimator to perform nearly as well as each regime’s oracle estimator, and outperform the other estimators considered. Specifically, the SWAG model is correctly specified for all four cases as each regime corresponds to particular limiting choices of parameters in the SWAG model. Furthermore, given that we consider problem dimension size similar to each population’s sample size, we expect the SWAG estimator to outperform the sample covariance estimators in all cases.

The results of the simulation study are presented in Table 3. In the table, the oracle estimator for each regime has a grey background, and the smallest two average losses are in bold font. In summary, the SWAG estimator performs best or second best in all cases considered except those in the Ho, Kr regime. In nearly all cases where the SWAG estimator has the second smallest loss, it is beat by the respective regime’s oracle estimator. Therefore, given that the oracle is unknown in practice, we conclude the SWAG model is particularly effective in accurate population covariance estimation. While the overarching conclusion of the performance of the SWAG model remains the same across most cases considered, the dynamics differ across regime. To gain a better understanding of the variation around the average losses displayed in the table, Figure 2 displays the empirical densities of the 50 $\bar{L}$ values for each regime in the case where $J=4,p=12$.

In the regime where population covariances are homogeneous across population and Kronecker structured, the oracle estimator $\hat{K}_{p}$ has the smallest average loss in the cases explored, as expected. Interestingly, for a given $J$, the average loss corresponding to $\hat{K}_{p}$ and $\hat{\mathbf{K}}$ for $p=12$ is less than that for $p=6$. This is a consequence of the scaling used within Stein’s loss function as, in this case, the difference between the first two terms of the loss increases by less than the increase in $p$. In general, the pooled MLE and group-specific Kronecker MLEs perform well in this regime, which is not surprising given the variety of estimators considered and the overlap in their underlying assumptions. The SWAG estimator tends to have a similar $\bar{L}$ to these estimators, and a notably smaller average loss than the sample covariance. Furthermore, in most cases considered, the empirical density of the loss corresponding to these three estimators ($\hat{\boldsymbol{\Sigma}},\hat{S}_{p},$ and $\hat{\boldsymbol{K}}$) overlap, as seen in Figure 2.

Given population covariances homogeneous across population and not Kronecker structured, the pooled MLE and SWAG estimators have the two smallest average losses. All other estimators have comparatively large average loss. This is particularly apparent in Figure 2 where the densities of the losses corresponding to $\hat{\boldsymbol{\Sigma}}$ and $\hat{S}_{p}$ are smaller than and distinct from the densities corresponding to the other estimators considered. In high dimension cases in particular, incorrectly assuming a Kronecker covariance in this regime results in estimators nearly as inaccurate as the sample covariances.

In the heterogeneous and Kronecker structured regime, the oracle estimator and the SWAG estimator have the two smallest average losses. Overall in this regime, these two estimators behave similarly and are notable improvements over the other estimators considered, as seen clearly in Figure 2. Again, in this regime, the average loss corresponding to the oracle estimator $\hat{\mathbf{K}}$ for $p=12$ is less than that for $p=6$ as a result of the scaling used in Stein’s loss function.

In the fourth, and perhaps most realistic, regime considered, where population covariances are heterogeneous across population and not Kronecker structured, the SWAG estimator has the smallest average loss in all cases except when $J=10,p=6$, where it has the second smallest average loss. In this regime, the oracle estimator is outperformed by all other estimators considered. This is not surprising given the combination of the simplicity of the exchangeable population covariances in conjunction with large instability of the sample covariance due in part to the small sample sizes considered for each regime. Again, we see estimators based on an incorrect assumption of separability result in distinctly larger loss than the more flexible SWAG estimator (Figure 2).

In total, these results support the conclusion that the SWAG model outperforms standard alternatives given the true structure of the population covariances is unknown. The flexibility of the SWAG model is particularly useful given the large error under estimators based on incorrect assumptions. In particular, wrongly assuming a Kronecker structure can result in an error nearly as large as that obtained under the sample covariance. While the extreme regimes considered do not necessarily reflect the truth in real-world scenarios, this simulation study highlights the flexibility of the SWAG model.

We demonstrate the usefulness of the SWAG model for estimating covariance matrices in multi-group matrix-variate populations by analyzing two data sets. In the first example, we perform a speech recognition task on a publicly available spoken-word audio dataset. In the second example, we analyze chemical exposure data which features small group-specific sample sizes.

In general, a data matrix for a single group $Y$ of dimension $n\times p$ can be decomposed into two orthogonal matrices such that one may be used for mean estimation and the other, based on centered data, for covariance estimation. As such, while a mean estimation step could be included in the proposed SWAG Metropolis-Hastings algorithm, we will estimate the covariance matrices based on centered data matrices throughout the applications. To elaborate, first, define the $n$-dimensional centering matrix $\mathcal{C}=I_{n}-\mathcal{P}_{1}$ where $\mathcal{P}_{1}=\mathbbm{1}_{n}\mathbbm{1}_{n}^{T}/n$ is a rank-1 idempotent projection matrix and $\mathcal{C}$ is a rank-$(n-1)$ idempotent projection matrix (Christensen, 2011). Then, note that,

where $\bar{y}$ is the length $p$ vector of column means of $Y$ and $\mathcal{C}Y$ is the residual matrix. Then, for $Y\sim N_{n\times p}(\mathbbm{1}_{n}\mu^{T},\Sigma\otimes I_{n})$, $\bar{y}\sim N_{p}(\mu,\Sigma/n)$ and $nS\sim\text{Wishart}_{p}(\Sigma,n-1)$ where $S=(\mathcal{C}Y)^{T}\mathcal{C}Y/n$ is the sample covariance matrix. $\bar{y}$ and $\mathcal{C}Y$ are uncorrelated, so $\bar{y}$ and $S$ are uncorrelated (Mardia et al., 1979). In this way, $\Sigma$ may be estimated using centered data $\mathcal{C}Y$.