A Bayesian Nonparametric model for textural pattern heterogeneity

Different classes of features have been proposed and interrogated for their utility in clinical applications. Image texture analysis, the focus of this article, endeavors to describe tumour heterogeneity through patterns of gray-level spatial dependence. The predominate process for obtaining features for texture analysis involves mapping the image domain into a gray-level co-occurrence matrix (GLCM). Avoiding multivariate analysis, the GLCM object is then further summarized to produce sets of discrete, but highly redundant features for subsequent analysis via multiple regression and machine learning algorithms.

A GLCM is a symmetric matrix defined over an image domain with cells comprised of counts that describe the distribution of co-occurring gray-scale valued pixels (or voxels) at a given offset and all directions (Haralick et al., 1973). More specifically, the $\left(l,h\right)^{th}$ entry of a GLCM represents the frequency with which a pixel with gray-level $l$ is present in a spatial location either horizontally, vertically or diagonally to the adjacent pixel with gray-level $h$. Figure 1 provides a simple illustration of the mapping from an image in gray-scale to the associated GLCM. More specifically, figure 1 (a) depicts a toy image with only four gray levels (0 through 3). The associated GLCM for this image is shown in figure 1 (b), where each $(l,h)$ cell counts the number of times the gray level $l$ occurs with gray level $h$ at an adjacent neighborhood of $l$. For example, to determine the value in cell $(2,0)$ in figure 1 (b), we count the co-occurrence of gray level 2 and gray level 0 (blue arrows) in figure 1 (a).

If we assume $K$ gray-level values (or bins), then the GLCM contains at most $K(K+1)/2$ unique cell counts if all directions are considered. This representation facilitates image standardization, which is useful in diagnostic oncology settings wherein targeted lesions may present with irregular shapes and sizes. The complex tumor images are thus mapped to a lower-dimensional standardized lattice, consisting of $K(K+1)/2$ structured, discrete random variables, thus simplifying model formulation and analysis.

Current practice further reduces the multivariate functional structure inherent to GLCMs to sets of summary statistics (e.g., GLCM-derived textural features) for subsequent analysis via multi-variable regression and machine learning algorithms. These reductive mappings are potentially limiting, as they may result in information loss. Moreover, many textural features are sensitive to rotations of the image and/or GLCM (Zhang et al., 2017). This lack of invariance impacts the robustness and reproducibility of feature-based approaches. Furthermore, important patterns for discerning the extent to which a malignant mass has proliferated throughout a region of interest are potentially masked with analysis based on the summary statistics.

To the best of our knowledge, statistical methods able to take into account the information provided by the complete GLCM objects, instead of GLCM-derived features, are lacking in cancer radiomics. Li et al. (2019) have recently proposed to analyze a GLCM as a multivariate object by means of a hierarchical Bayesian spatial model. The technique, based on Gaussian Markov random field priors, improved classification accuracy when compared with radiomics feature-based classifiers.

As a limitation, however, Li et al. (2019) employed a Gaussian transformation of the observable GLCM data prior to analysis, and thus lacked a formal model for the structured count data. The Gaussian transformation may fail to account for the typical skewness and zero-inflation of the distribution of the GLCM counts, thereby resulting in information loss. Also, the assumption of a Gaussian Markov random field intrinsically implies stationarity of the spatial process which is a restrictive assumption. Finally, Li et al. (2019) investigate the classification of tissue samples only in a supervised learning context, by means of binary classifiers trained by known pathology status. The objective of image-analysis is often discerning distinct patterns intrinsic to diverse tissue environments in settings wherein the specific clinical subtypes may be unknown. These limitations motivated the development of the proposed model presented herein, which extends GLCM-based textural pattern identification and analysis beyond supervised learning.

Here, we present a computationally feasible Bayesian nonparametric framework for inference of the GLCM as a multivariate count object. The modeling approach captures patterns of spatial dependence intrinsic to the GLCM lattices itself, thereby avoiding further processing steps to produce reductive and redundant GLCM summary statistics. A feature-based analysis requires careful variable selection and/or multiple comparisons, and often limits the reproducibility of the resultant classifiers. By way of contrast, our proposed model for the GLCM object appropriately accounts for over and under dispersion of the observed counts, and simultaneously adjusts for existing spatial autocorrelation at nearby cells. Model specification is facilitated by a rounded kernel mixture model (Canale and Dunson, 2011), such that the observed spatial distribution of GLCM cell counts are assumed to arise from a latent continuous random vector defined on the GLCM lattice. The rounded kernel framework allows a flexible and robust estimation of the data generating distribution for count data(Canale and Prünster, 2017), with respect to typical Poisson-Gamma mixtures (Karlis and Xekalaki, 2005; Guindani et al., 2014). Considering conjointly intra-lesion heterogeneity of the GLCMs objects with spatial correlation of the lattice counts, we further model the multivariate latent vectors using a spatial Dirichlet process (Gelfand et al., 2005). This is achieved by defining the base distribution with a conditionally auto-regressive (CAR) Gaussian structure. Hereafter, our method is referenced as the hierarchical rounded Gaussian spatial Dirichlet process (hierarchical rounded Gaussian SDP) model.

Flexible for actual applications, our modeling framework facilitates pattern identification with analyses of cohorts with heterogeneous lesion sizes. Our unsupervised Bayesian model is applied to identify textural patterns intrinsic to adrenal lesions with CT. Data was obtained in a retrospective study of patients who underwent concordant biopsy and imaging based on a dynamic sequence of CT scans with and without administration of contrast (Ng et al., 2018a, b, c). To elucidate whether textural subtypes identified exhibit pathological orientation, subtypes are evaluated for correspondence with each lesion’s diagnosed pathological status. Additionally, simulation is used to compare performance to a class of machine-learning approaches currently promoted with cancer radiomics applications. Our Bayesian nonparametric method demonstrates robustness to capturing the underlying spatial patterns of GLCM objects, whereas commonly employed feature-based algorithmic strategies fail to identify patterns in the multivariate structure intrinsic to GLCMs.

The manuscript is organized as follows. Section 2 describes the motivating adrenal lesion CT imaging study. Section 3 introduces the proposed hierarchical rounded Gaussian spatial Dirichlet process model. Sections 4 and 5 present our case and simulation studies.

Figure 2 depicts enhancement patterns on the Hounsfield unit domain of ROIs comprising adrenal lesions as well as their corresponding gray-level co-occurrence matrices for three representative subjects. The top image, reflects a high degree of co-occurrence at higher gray levels, which exploratory analyses suggest to be representative of cases of malignancies within our case study. Conversely, the bottom image tends to be indicative of benign lesions. The central image, for which high co-occurrences were observed at middle gray-levels, depicts an intermediate case. Importantly, the GLCM cells exhibit high correlation with their adjacent neighbors, as depicted by the fairly smooth changes observed when transversing away from the peak cell counts. This suggests that lesion heterogeneity may be differentiable by spatial pattern in the GLCM object by leveraging this dependence structure with an appropriate multivariate model. Here, we propose to describe the observed heterogeneity of the GLCM matrices by modeling the random effects $\bm{\theta_{t}}$ with a spatial Dirichlet process (Gelfand et al., 2005). In Bayesian modeling, nonparametric priors have been often employed for unsupervised learning, to achieve improved inference on subject-specific parameters by borrowing information within population subgroups. The spatial Dirichlet process assumes that each spatial random effect $\bm{\theta_{t}}$ is a sample from an a.s. discrete probability distribution $G$,

$\displaystyle\bm{\theta_{t}}\mid G\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}G,\qquad t=1,2,..,T.$

This is a realization of a Dirichlet Process,

$$G\sim DP(\upsilon G_{0}),$$

where $G_{0}$ denotes a base (or centering) parametric spatial model, such that $E(G(A))=G_{0}(A)$ for any measurable set $A\subset\mathcal{R}^{n}$, and $\nu>0$ is a precision parameter, characterizing the prior uncertainty of $G$ around the parametric model $G_{0}(\cdot)$. Since $G$ is almost surely discrete, there is a positive probability that some of the vectors $\bm{\theta_{t}}$’s take the same value for some of the subjects, thus inducing a probability-based clustering of the GLCM matrices in homogeneous subgroups.

Textural patterns in the GLCM object can be identified through model specifications that acknowledge spatial dependencies among adjacent cell counts. In particular, here we specify the base parametric model $G_{0}$ as a $n$-variate normal distribution by assuming a conditionally auto-regressive (CAR) structure (Besag, 1974; Banerjee et al., 2014)

$\displaystyle G_{0}(\cdot\mid\sigma^{2},\rho)=MVN(\cdot\mid\bm{0},(\bm{D}-\rho\bm{W})^{-1}\sigma^{2}),$

where $\bm{W}$ is an adjacency matrix with $W_{kd}$ indicating whether lattice elements $k$ and $d$ are adjacent neighbors; $\bm{D}$ is a diagonal matrix with $k$th entry denoting the number of neighbors for lattice $k$; $\rho$ is a smoothing parameter controlling the degree of spatial dependence and $\sigma^{2}$ is a variance parameter (Banerjee et al., 2014). The covariance structure implied by the spatial Dirichlet process construction is of particular relevance. As a matter of fact, integrating out the unknown $G$, the marginal covariance matrix of the latent vector $\bm{y}(s)$ can be computed as

$\displaystyle cov(\bm{y(s)})$

$\displaystyle=$

$\displaystyle\gamma^{2}\sigma^{2}(\bm{D}-\rho\bm{W})^{-1}+\tau^{2}\bm{I},$

i.e. constant over the entire spatial domain. However, any realization $G$ of the spatial Dirichlet process is a.s. discrete, i.e. $G=\sum_{l=1}^{\infty}\,p_{l}\,\delta_{\bm{\theta^{*}_{l}(s)}}$, where $\delta_{\theta}$ denotes a point mass at $\theta$, $w_{l}$ is a sequence of weights, and $\bm{\theta_{l}^{*}(s)}$ are sequences of unique realizations from the base distribution $G_{0}$ (independently of the weights). Then, for any given $G$, the conditional covariance matrix of $\bm{y}(s)$ has elements

$\displaystyle cov(y(s_{i}),y(s_{j})|G)$

$\displaystyle=$

$\displaystyle\gamma^{2}[\sum_{l=1}^{\infty}p_{l}\theta^{*}_{l}(s_{i})\theta^{*}_{l}(s_{j})-(\sum_{l=1}^{\infty}p_{l}\theta^{*}_{l}(s_{i}))(\sum_{l=1}^{\infty}p_{l}\theta^{*}_{l}(s_{j}))].$

Therefore, the covariance changes over the spatial domain, and in general it is not constant in any two locations. Thus, the spatial Dirichlet process formulation provides a flexible model to describe the spatial dependence observed in the GLCMs, including “non-stationary” patterns. Furthermore, the effect of the scaling factor $\gamma$ on the covariance matrix reflects the general idea that regions of smaller size are characterized by smaller variance of the observed gray-scales over the spatial domain.

The hierarchical model is fully specified with the prior distributional assumptions for the remaining parameters. More specifically, for the regression coefficients, we can assume a multivariate Gaussian prior with mean $\bm{\beta_{0}}$ and covariance matrix $\bm{\Sigma_{\beta}}$. To promote maximal data learning, it is common to set $\bm{\beta_{0}}$ as $\bm{0}$ and $\bm{\Sigma_{\beta}}$ as a diagonal matrix $\sigma^{2}_{\beta}\,\bm{I}$ with large $\sigma^{2}_{\beta}$, to allow a vaguely non-informative specification spanning a large region of the parameter space. The parameters $\tau^{2}$ and $\sigma^{2}$ are assigned, respectively, (conditionally) conjugate inverse-gamma priors, $IG(\tau^{2}\mid a_{\tau},b_{\tau})$ and $IG(\sigma^{2}\mid a_{\sigma},b_{\sigma})$ (with means $b_{\tau}/(a_{\tau}-1)$ and $b_{\sigma}/(a_{\sigma}-1)$, respectively). The hyperparameters $a_{\tau}$, $b_{\tau}$, $a_{\sigma}$, $b_{\sigma}$ are typically fixed to represent vague prior information. A $Gamma(a_{\upsilon},b_{\upsilon})$ prior (with mean $a_{\upsilon}/b_{\upsilon}$) is imposed for $\upsilon$ as the corresponding mixture of gamma full conditional facilitates MCMC computation (Escobar and West, 1995). Lastly, we follow prevailing literature and consider a $Uniform(0,1)$ distribution for the $\rho$ parameter in the CAR model (Lee, 2013).

Figure 3 provides a graphical representation of the proposed hierarchical model for the GLCM counts. Conditional on the latent $\bm{y_{t}}$’s, the sampling distribution can be expressed as

$\displaystyle\noindent P(\bm{Z}|\bm{y})$

$\displaystyle\propto$

$\displaystyle\prod_{t=1}^{T}\prod_{i=1}^{n}I\left\{a_{z_{t}(s_{i})}\leq y_{t}(s_{i})<a_{z_{t}(s_{i})+1}\right\},$

where $I(\cdot)$ denotes the indicator function, for fixed cut-points $a_{0}$, $a_{1}$, …$a_{\max_{t}\,\{z_{t}\}}$. Then, the hierarchical model on the latent vector $\bm{y_{t}}$ can be summarized as follows:

	$\displaystyle\noindent\bm{y_{t}}\mid\bm{\beta},\bm{\theta_{t}},\tau^{2}$	$\displaystyle\sim MVN\left(\bm{y_{t}}\mid\bm{x_{t}}\bm{\beta}+\gamma_{t}\bm{\theta_{t}},\tau^{2}\bm{I}\right),\qquad t=1,2,..,T$
	$\displaystyle\bm{\theta_{t}}\mid G^{(n)}$	$\displaystyle\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}G^{(n)},\qquad t=1,2,..,T$
	$\displaystyle G^{(n)}\mid\upsilon,\sigma^{2},\rho$	$\displaystyle\sim DP(\upsilon G_{0}^{(n)})$
	$\displaystyle G_{0}^{(n)}(\cdot\mid\sigma^{2},\rho)$	$\displaystyle\equiv MVN(\cdot\mid\bm{0},\sigma^{2}(\bm{D}-\rho\bm{W})^{-1})$
	$\displaystyle\bm{\beta}$	$\displaystyle\sim MVN\left(\bm{\beta}\mid\bm{\beta_{0}},\bm{\Sigma_{\beta}}\right)$
	$\displaystyle\tau^{2}$	$\displaystyle\sim IG(a_{\tau},b_{\tau})$
	$\displaystyle\sigma^{2}$	$\displaystyle\sim IG(a_{\sigma},b_{\sigma})$
	$\displaystyle\upsilon$	$\displaystyle\sim Gamma(a_{\upsilon},b_{\upsilon})$
	$\displaystyle\rho$	$\displaystyle\sim Uniform(0,1).$

In this section, we outline the main steps of the MCMC algorithm used to obtain posterior samples from the proposed hierarchical rounded Gaussian SDP model, and we discuss how we can summarize the posterior inference about the heterogeneity of the GLCM objects across subjects.

To obtain appropriate synthetic data capable of representing the spatially-oriented patterns within GLCMs exhibited in our diagnostic cancer study of adrenal lesions, we generated GLCMs constructed with 16 gray-levels with normalized element-wise probability densities. Our simulation study was conducted using two data generating models. The first model assumed that density patterns within the GLCM objects were a sample from a bivariate normal distribution with

$$\mu=(2+c,14-c)^{{}^{\prime}}\quad\text{and}\quad\Sigma=s\begin{bmatrix}1&-0.7\\ -0.7&1\\ \end{bmatrix}.$$

Various simulation scenarios were obtained as different combinations of the shift parameter $c$ and the scale parameter $s$, as described further below. A final smoothed event rate surface, wherein the event rate is regarded as the probability of observing a co-occurrence count within a given cell, was obtained for each simulated GLCM by smoothing over the Gaussian-derived empirical rate surface calculated in proportion to the number of generated points in each grid and with respect to the total number of generated points of the entire grid surface. Additionally, to reflect inter-patient heterogeneity in image size we scaled the rate surface over the discrete GLCM space by a random integer that is sampled from 500 to 20000. The scale densities were then rounded to the nearest integer value to yield simulated counts. The simulation design yields a distribution over the GLCM space that generates reasonably small and large GLCM count surfaces, representing both over- and under-dispersed scenarios. Simulation scenarios were then obtained by generating GLCMs under different choices of the mean shift parameter $c$ and variance scale parameter $s$. In particular, for a given $s$ we generated GLCMs as a mixture of five components represented by $c\in\{5,5.5,6,6.5,7\}$, such that the center of the generating bivariate normal distribution shifts their location gradually from north-west toward south-east along the $135^{\circ}$ diagonal of the GLCM. We also evaluated performance by varying $s\in\{10,12,15,16,17\}$, with each value of $s$ representing the extent to which the noise covers the true signal in the underlying spatial patterns intrinsic to each tissue class. As we set the mean vectors of the bivariate normal distribution of the different components closer to each other, increasing the value of $s$ diminishes the extent to which the spatial patterns are differentiable by any method. For each $s$, we generated 100 subjects in total, with 20 patients in each of the five subtypes defined by $c$. To interrogate robustness to the Gaussian assumption, in the Supplementary material, Section B, we also present a simulation using the bivariate skew normal distribution of Azzalini and Capitanio (1999) as the latent generator. The results are consistent with those described below for the bivariate normal setting.

We fitted the simulated GLCMs using the proposed HRGSDP method, where the model hyper-parameters were set to obtain vague prior specification as in the case study. MCMC was implemented with 20,000 iterations, the first 10,000 of which were discarded as burn-in. Convergence was assessed by using the Geweke diagnostic [31] with implementation using the R package ‘coda’.

We then considered five alternative clustering approaches, which are typically employed for the analysis of imaging data in oncology studies. More specifically, we considered (i) hierarchical clustering based on a few GLCM-derived texture features, which are reported routinely in radiomics studies (such as contrast, correlation, homogeneity, energy and entropy) (Kumar et al., 2012). We refer this method as ‘FeaHC’. Additionally, for the same GLCM-derived texture features, (ii) spectral clustering, referenced as ‘FeaSC’, (iii) K-means clustering, referenced as ‘FeaKM’ (iv) consensus clustering algorithms, referenced as ‘FeaCO’ and (v) Gaussian mixture model, referenced as ‘FeaGMM’ were also implemented as typically reported in clinical imaging studies (Gensheimer et al., 2015; Wu et al., 2016; Parmar et al., 2015; Zhou et al., 2018).

In contrast to our unsupervised approach, all of the competing methods above require the pre-specification of the number of clusters in advance of analysis. In order to optimize their performance , we adopted the true number of clusters as the ground truth, i.e. we assumed five simulated subtypes. That is, clusters were obtained for each simulated dataset by cutting the resultant dendrogram with rank equal to the true number of clusters. Of course, in real data settings, such knowledge would not be available. By way of contrast, no truth-based specification was used for ascertaining subtypes using our Bayesian method, which provided an artificial advantage favoring competing approaches.

While our modeling framework was devised for the purpose of unsupervised clustering, comparing the relative quality of competing sets of unsupervised clusters is abstruse in the absence of an underlying mechanism for understanding realizations of the objects that are the basis for cluster analysis. Therefore, our simulation study was devised under the assumption that GLCM objects arise from disparate true classes of lesions. Each class was assigned a true probability distribution for generating individual GLCMs, as described in Section 5.1, providing a basis for comparison. Specifically, method performance was evaluated in relation to the matching matrix constructed with true memberships as rows and the predicted membership as columns such that the $(i,j)$th entry measures the counts of the $i$th class subject with membership assigned by the clustering approach $j$. The matching matrix was computed with our proposed method and the above mentioned approaches for each simulated cohort on 100 replicate simulations for each value of $s.$ To evaluate and compare the performance of our proposed HRGSDP method and the other competing approaches mentioned above, we used two performance metrics: a measure of the mis-assignment rate and the Pearson chi-square test statistic obtained from the resultant matching matrices.

With our proposed Bayesian nonparamteric approach, a membership assignment was obtained for each subject by posterior clustering of $\bm{\theta}$ (see Section 3.2.2). Then, the assigned membership along with its true class (defined by $c$) were interrogated to measure the extent to which GLCMs arising from truly identical pattern distributions tended to cluster together. The competing methods require the pre-specification of the total number of clusters to assign the GLCMs to cluster labels. However, the HRGSDP is unsupervised, as it provides a data-driven estimate of the cluster rank. Nevertheless, in the following, we measure the performance of the competing methods with respect to the true rank (5) used to simulate GLCM patterns, which would be unknown in practice.

Methods were compared on the basis of their mis-assignment rate, which was defined to characterize the proportion of GLCM objects that cluster with truly uncommon patterns. This is calculated by tabulating the distribution of model-derived predicted cluster assignments by their true class, and identifying the mode class for each cluster (class with maximum proportion). For each predicted cluster, the number of assignments to other classes (not the mode) is taken as an error. The errors are summed and normalized by the sample size. Then, the resultant mis-assignment rate was defined as the proportion of mis-assigned subjects overall.

The Pearson chi-square test statistic was also used to measure the homogeneity of true class membership and the assigned class membership. Thus, larger values of the test statistic, which represent larger deviations from the null random frequency distribution, characterize better performance as cluster assignments derived from clustering techniques become more consistent with their true subtype memberships.

Figure 6 reports the mis-assignment rate distributions obtained in our simulation study as a function of the values of the dispersion parameter $s\in\{10,12,15,16,17\}$.

Among all currently employed feature-based approaches, the mean mis-assignment rates are similar, ranging from $32.4\%$ to $58.7\%$, with the Gaussian mixture model (FeaGMM) having the best performance with $s=15$ ($32.4\%$). The best performance for consensus clustering (FeaCO), hierarchical clustering (FeaHC), K-means (FeaKM) and spectral clustering (FeaSC) methods are $53.9\%$ with $s=17$, $42.8\%$ with $s=17$, $41.9\%$ with $s=16$ and $38.6\%$ with $s=16$, respectively. In all simulated scenarios, the Gaussian mixture model (FeaGMM) method yields the minimum mis-assignment rate among the feature-based approaches. As a contrast, the performance of the proposed method is globally superior at smaller values of $s$ (0.1% for $s=10$ and $5.8\%$ for $s=15$) ,i.e. for high signal-to-noise scenarios. As the random noise increases, the performance of the proposed method decreases but it is still better than or at least comparable to the performances of the competing methods, resulting in mean mis-assignment rates of $12.5\%$ for $s=16$ and $19.2\%$ for $s=17$. This behavior appears to be due to the combination of two factors. On the one hand, high noise results in imprecise cluster estimates. On the other hand, for large values of $s$, the noise dominates the signal conveyed by the true mean spatial patterns of dependencies that characterize the true GLCM subtypes. As our simulations demonstrate, different subtypes are best identified by large differences in the values of the mean location parameter $c$.

Figure 7 reports the distributions of the Pearson chi-square test statistic obtained as a function of the dispersion parameter $s\in\{10,12,15,16,17\}$, respectively.

We should note that, in our simulation scenario, the maximum possible value of the Pearson chi-square test statistic is 400. This is obtained whenever all clusters identified by a model consist entirely of GLCMs arising from the same true pattern distribution. We refer to this as perfect discrimination. The means of the statistic for all feature-based approaches are similar, ranging between 83 and 225. Among them, the Gaussian mixture model yielded the best clustering performance, reaching a maximum of 225 with $s=17$. The maximum values of the Pearson chi-square test statistic obtained for the consensus clustering (FeaCO), hierarchical clustering (FeaHC), K-means (FeaKM) and spectral clustering (FeaSC) methods were 124, 162, 160 and 177 at $s=17$, respectively.

Results for the HRGSDP model far exceeded those resulting from the above machine learning clustering approaches at all levels of signal-to-noise. As a matter of fact, results were perfect for $s<15$ with our proposed HRGSDP model resulting in a mean Pearson chi-square test statistic of 400 at $s=10$ and $s=12$. In the presence of higher levels of random noise, the statistic was 373 at $s=15$, 341 at $s=16$ and 309 at $s=17$, and thus globally superior.

The trends in the diagnostic performance that the HRGSDP method exhibited in our simulation study conform to intuition as it pertains to identification of spatial discriminatory patterns and random noise. Overall, the proposed HRGSDP method outperformed the machine-learning methods currently used in the literature for reasonable signal-to-noise ratios. The algorithmic-oriented techniques fail to take into account the multivariate lattice structure intrinsic to GLCMs, and thus are unable to properly capture the spatial dependencies and the distributional heterogeneity inherent to a sample of GLCM objects. In the presence of high signal-to-noise in the true GLCM generating distribution, as characterized by small values of $s$, the patterns conveyed in GLCMs were more separable from each other. This resulted in enhanced discrimination results for the HRGSDP method which leverages the spatial location information. As the variance scale parameter $s$ increases, the distinction of spatial patterns among different $c$ are diminished. The methods yielded similar results when adopting the skew bivariate normal distribution as the latent generator. Full results arising from the skew normal distribution are described in detail in Section B of the Supplementary materials.