Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
57,679 characters · 16 sections · 46 citation commands
Bayesian Clustered Coefficients Regression with Auxiliary Covariates Assistant Random Effects
Analysis of spatial data referenced over different locations has received widespread attention in many fields such as environmental science hu2018stat,yang2019bayesian, social science bradley2018computationally, and biostatistics xu2019latent. There are two major approaches to analyzing spatial data. The first approach, where spatial variations in the outcome are accounted for by an additive spatial random effect term at each location, has been studied both for the linear model cressie1992statistics and generalized linear model diggle1998model settings. The second approach, a regression model with spatially varying coefficients, is developed to capture spatial variations within the covariate effects themselves. For the spatially varying coefficients model, there are two major approaches for estimation of the regression coefficients. One is geographically weighted regression brunsdon1996geographically, the basic idea of which is to assign different weights to the observations based on a certain measure of distance between them and the target location. This work also has various extensions in generalized linear regression nakaya2005geographically and analysis of survival data hu2018modified,xue2019geographically, as well as under the Bayesian paradigm ma2019bayesian. It has been further extended to multiscale models that allow kernels of different variables to be on different scales, which greatly enhances its model flexibility fotheringham2017multiscale. Another major approach is to give a Gaussian process prior to the spatially varying coefficients gelfand2003spatial, which provides a natural and flexible way to view the coefficient surface as a realization from a spatial process. This work is universally applied into different models such as Poisson regression model reich2010bayesian and survival models hu2020comparison.
Most recently, heterogeneous covariate effects in many different fields, such as real estate applications, spatial econometrics, and environmental science are receiving increasing attention. For example, a country's big cities and small cities could be put into separate clusters and analyzed, as each subgroup share more similarities in development patterns, and thus similar covariate effects in regression models can be expected. Such clustering information is of great interest to regional economics researchers. Existing frequentist approaches include those based on scan statistics kulldorff1995spatial, two-step spatial hypothesis testing lee2017cluster,leespatial2019, and a penalized method based on minimum spanning tree li2019spatial. Under the Bayesian paradigm, an integrated framework ma2019bayesianspatial to detect clusters in the covariate effects as well as producing the parameter estimates for spatially dependent data are proposed, in which clustering is done via the Dirichlet process mixture model neal2000markov,ishwaran2002exact. The DPM, however, has turned out to produce extremely small clusters and make the estimation for number of clusters inconsistent miller2013simple. The mixture of finite mixture (MFM) model proposed by miller2018mixture provides a remedy to over-clustering problem for Bayesian nonparametric methods. Both DPM and MFM allow for uncertainty in the number of clusters instead of relying on a given number, which needs to be tuned based on certain criteria.
While spatial random effects have been used to account for the influence of geographical proximity on the similarity of outcomes for two neighboring observations, in existing approaches, the correlation between spatial random effect terms only depends on distance, and all other factors are ignored. For improvement in describing spatial correlation, some works have been done to bring in auxiliary information to help the estimation of spatial regression model white2009stochastic,lee2014bayesian,gao2019bayesian. These works focus on the estimation of the edge based on some covariates or spatial structure. The regression relationship between the $p$-dimensional covariance matrix and auxiliary information has been explored by zou2017covariance and liu2020semiparametric, which can help reveal the true correlation structure of spatial data. Given the importance of spatial random effects in a spatial regression model, their covariance structure need to be appropriately specified.
In this work, we propose a Bayesian clustered linear regression model with MFM, which when compared to the DPM, consistently estimates the number of clusters. The estimation of parameters remains precise. To our best knowledge, we firstly introduce the MFM in clustered coefficients regression model. In addition, a weighted average correlation structure of auxiliary covariates information is incorporated in linear mixed effects regression model based on Dirichlet prior.
The remainder of the paper is organized as follows. In Section (ref), we present the spatial clustered linear regression model with MFM. In Section (ref) , a MCMC sampling algorithm based on nimble de2017programming and post MCMC inference are discussed. Extensive simulation studies are performed in Section (ref). For illustration, our proposed methodology is applied to Georgia housing cost data in Section (ref). We conclude this paper with a brief discussion in Section (ref).
The basic geostatistical model gelfand2016spatial for observations made over a spatial domain can be written as:
where $\bm{Y} = (Y(s_1),\ldots, Y(s_n))$ denotes the $n$-dimensional vector of spatial responses observed at locations $s_1,\ldots,s_n$, $\bm{X} =
$ denotes the~$n\times p$ matrix of covariates, $\bm{\beta}$ is the vector of coefficients, $\bm{w} = (w(s_1),\ldots, w(s_n))^\top$ is a vector of spatial random effects, and $\bm{\epsilon}\sim MVN(\bm{0}, \tau_y^{-1}\bm{I})$ is the ``nugget effect'' with $\tau_y$ being the precision of the response carlin2014hierarchical. The above spatial regression model can be formulated alternatively as
where $\bm{w}$ denotes the spatial structure with a covariance matrix $\bm{\Sigma}_{\bm{w}}$, and N and MVN denote the univariate and multivariate normal distributions, respectively. The covariance matrix $\bm{\Sigma_w}$ is often defined to be $\sigma_w^2\bm{H}$, with $\bm{H}$ constructed using the great circle distance (GCD) matrix among different locations via the following three popular weighting schemes:
where $\phi$ is a bandwidth parameter that controls the spatial correlation.
Instead of using spatial random effects, i.e., location-wise intercepts to account for spatial variations in $\bm{Y}$, the spatially varying coefficients model gelfand2003spatial attributes such effects to variations in the parameters over the spatial domain, i.e., the parameter $\bm{\beta}$ itself varies. Such a model is formulated as:
where $\bm{X}(s_i)$ is vector of covariates at location of the $i$th subject $s_i$, and $\widetilde{\bm{\beta}}(s_i)$ is assumed to be generated from a $p$-variate spatial process model. With observations $(Y(s_i),\bm{X}(s_i))$ for $i=1,\ldots, n$, the model can be written as
where $\bm{Y}=(Y(s_1), \ldots, Y(s_n))^\top$, $\bm{X}^\top$ is an $n\times (np)$ block diagonal matrix whose $i$-th diagonal entry is $\bm{X}^\top(s_i)$, $\widetilde{\bm{\beta}} = (\widetilde{\bm{\beta}}(s_1)^\top, \ldots, \widetilde{\bm{\beta}}(s_n)^\top)^\top$, and $\bm{\epsilon}\sim \mbox{MVN}(0,\tau_y^{-1}\bm{I})$ with $\bm{I}$ being the identity matrix. To characterize the $p$-variate spatial process that generates $\widetilde{\bm{\beta}}(s_i)$, we rewrite the model as
where $\bm{\mu}_{\bm{\beta}}$ is a $p \times 1$ vector, $\bm{H}(\phi)$ is an $n\times n$-dimensional matrix measuring spatial correlations between the $n$ observed locations, $\bm{T}$ is a $p \times p$ covariance matrix associated with an observation vector at any spatial location, and $\otimes$ denotes the Kronecker product.
In the above formulation, each location has its own $\bm{\beta}$ vector. However, such a model could be too flexible, as there are certain regions that have very similar $\bm{\beta}$ values. From a modeling perspective, clustering such regions and having them all share one parameter vector encourages a parsimonious model without compromising the model's explanatory power. Another potential drawback of formulation in ((ref)) is that all variations are accounted for by $\widetilde{\bm{\beta}}$, and spatial random effects are ignored. However, the random effects term is rather important, and influenced not only by distance but also other factors such as demographics, transportation, etc. Therefore, a model with clustered coefficients, and also random effects terms that help account for the intricate connections between regions is desired.
Based on the heterogeneity pattern, we focus on the clustering of spatially-varying coefficients. A latent clustering structure can be introduced to accommodate the spatial heterogeneity on parameters of sub-areas. Let us denote the cluster belongings of the $i$th observation as $z_i$ for $i=1,\ldots, n$. The Dirichlet process ferguson1973bayesian offers a nonparametric approach for capturing heterogeneity effects in the data. The DP prior for the cluster belonging of the $i$th observation can be written as
where $k \rightarrow \infty$, $\pi_h$ denotes the random probability weight, $\delta_h$ is the Dirac $\delta$ with point mass at $h$, and $\alpha$ is the concentration parameter. The first equation in (ref) can be expressed equivalently as a multinomial distribution:
and we will use the $\sum_{h=1}^k \pi_h \delta_{h}$ notation for simplicity for the rest of this paper. The joint distribution for $z_1,\cdots,z_n$ can also be written as a conditional distribution, known as the Chinese restaurant process pitman1995exchangeable, neal2000markov. The distribution of $z_i$ is marginally represented by the stick-breaking construction of sethuraman1994constructive as
miller2013simple showed that the posterior distribution on the number of clusters does not converge to the true number of components, and extraneous clusters are often produced by the CRP. Later, miller2018mixture proposed a modification of the CRP called a mixture of finite mixtures (MFM) model to circumvent this issue, which can be formulated as
where $p(\cdot)$ is a proper probability mass function (p.m.f) on $\{1, 2, \ldots\}$. A default choice of $p(\cdot)$ is a $\mbox{Poisson}(1)$ distribution truncated to be positive miller2018mixture, which is assumed through the rest of the paper. Analogous to the stick-breaking representation in (ref), the MFM also has a similar construction. If we choose $k-1 \sim \mbox{Poisson}(\lambda)$ and $\gamma=1$ in (ref), the mixture weights $\pi_1,\cdots,\pi_k$ can be constructed as:
Gibbs samplers are easily constructed in stick-breaking framework ishwaran2001gibbs. For ease of exposition, we refer to the formulation in (ref) as $\text{MFM}(\gamma,\lambda)$.
In the regression model (ref) where spatial random effects are present, their covariance structure often depend on the geographical distance between pairs of locations, which in general indicates that closer locations have stronger correlation, and is in accordance with Tobler's first law of geography that “everything is related to everything else, but near things are more related than distant things”. In most economics problems, however, spatial proximity might not be the sole indicator for similarity, as there can be geographically distant locations that share similar demographical characteristics. For example, while the GCD between New York City and Albany is only 135 miles, which is far less than the 2569 miles between New York City and San Francisco (calculated using ggmap and geosphere packages in R), the population density of Albany is only 4525.3 per square mile, which is far smaller than those of New York City and San Francisco, which are, respectively, 27,709.4 and 19,104.4 per square mile worldpopulation2020. To incorporate such similarities into the covariance structure for random effects, motivated by covariance regression zou2017covariance and Bayesian model averaging raftery1997bayesian, we propose the following auxiliary covariates assistant covariance (ACAC) matrix for random effects in a mixture regression model:
where $\sigma^2$ is a constant accounting for the overall magnitude of variance, $\bm{w}=(w(s_1),\cdots,w(s_n))^\top$, and $\bm{W}(\bm{Z}_j),j=1,\cdots,J$ is the similarity matrix of the $j$th auxiliary covariate. Entries of the similarity matrix have values between 0 and 1, and are usually decreasing with respect to the absolute difference between values of the auxiliary covariates. The three aforementioned weighting schemes in (ref) can be used to define $\bm{W}(\bm{Z}_j)$. For example, an exponential decay similarity matrix $\bm{W}(\bm{Z}_j)$ can be constructed so that its $(\ell, \ell^{'})$-th element is
where $\kappa_j>0$ is the range parameter for the exponential kernel, and $|\cdot|$ denotes the Euclidean distance. In order to solve the identifiability issue, we set a constraint for $\alpha_0,\alpha_1,\cdots,\alpha_J$:
Proof for this proposition is directed to Supplemental Section S.1.
Based on the constraint in (ref), a Dirichlet prior is assigned to $\alpha_0,\cdots,$ $\alpha_J$. The prior distribution of $\alpha_0,\alpha_1,\cdots,\alpha_J$ is given as
where $\text{Dirichlet}(\nu)$ is the Dirichlet distribution with parameter $\nu$.
Combining the MFM and ACAC matrix, we have our final spatial MFM clustered regression model with ACAC hierarchically as follows, for $i = 1, \cdots, n$:
where $\bm{w} = (w(s_1), \cdots, w(s_n))^\top$, $\bm{\beta}_{z_i} = (\beta_{z_i1}, \cdots, \beta_{z_ip})^\top$ with $p$ being the dimension of the covariates $\bm{X}(s_i)$, $\bm{W}(\cdot)$ is the similarity matrix of the corresponding auxiliary covariate defined in ((ref)), and MFM the clustering method introduced in Section (ref). The choice of hyperparameters will be discussed in Section (ref).
Let $\bm{\theta}= \{(\tau_y, \mu_{z_i\ell}, \tau_{z_i\ell}, \sigma^2, \lambda, \kappa_j, \bm{\alpha}): i = 1,\cdots,n; \ell = 1,\cdots,p; j=1,\cdots,J\}$ denote the set of unknown parameters in the proposed model, and we assume that they are independent a priori. Therefore, we assign commonly used priors for these parameters: $\tau_y \sim \text{Gamma}(1,1)$, $\mu_{z_i\ell} \sim \mbox{N}(0, 1)$, $\tau_{z_i\ell} \sim \text{Gamma}(1,1)$, $\sigma^2 \sim \mbox{InverseGamma(1, 1)}$, $\lambda \sim \text{log-normal}(0, 1)$, $1/\kappa_j \sim \text{Gamma}(1, 1)$, and $\bm{\alpha} \sim \text{Dirichlet}(1,$ $1, \cdots)$. With the prior distributions specified above, the posterior distribution of these unknown parameters based on the data $D = \{Y(s_i), \bm{X}(s_i),$ $\bm{Z}_j\}$ is given by
The analytical form of the posterior distribution of $\bm{\theta}$ is unavailable. Therefore, we employ the Markov chain Monte Carlo (MCMC) sampling algorithm to sample from the posterior distribution, and then obtain the posterior estimates for the unknown parameters. Computation is facilitated by the nimble package in R, which uses syntax similar to WinBUGS and JAGS, but generates C++ code for faster computation. With the nimble package, sampling algorithms for the parameters are default samplers. For parameters $\bm{\alpha}$, a random-walk Dirichlet sampler is used. For $\bm{\beta}_{z_i}$, a random-walk block sampler is used. The same sampler is used for parameter $\bm{w}$. The conjugate sampler is used for $\tau_y$ and a categorical sampler is used for $z_i$, while for the other parameters, a random-walk sampler is applied.
In the proposed spatial regression model, since the covariance matrix of the spatial random effects can be constructed in different ways, including the unity, exponential, and Gaussian weighting schemes in (ref), a model selection criterion needs to be used for deciding which form of the covariance matrix is the most suitable for the data. A commonly used Bayesian model selection criterion, logarithm of the pseudo-marginal likelihood ibrahim2013bayesian can be used for this purpose. The LPML can be obtained through the conditional predictive ordinate (CPO) values, which are the Bayesian estimates for the probability of observing $Y_i$ in the future after other observations are made. Let $Y^*_{(-i)} = \{Y_j: j = 1, \cdots, i-1, i+1, \cdots, n\}$ denote the observations with the $i$th subject response removed. The CPO for the $i$th subject is defined as:
where
with $c(Y^*_{(-i)})$ being the normalizing constant. As discussed in chen2012monte, $\mbox{CPO}_i$ is also called the cross-validated predictive density, and Equation (ref) is essentially integrating the predictive distribution of $y(s_i)$ given the rest of the observations. It is a useful quantity for model checking, as it describes how much the observation at location $s_i$ supports the model. An equivalent expression for $\mbox{CPO}_i$ is:
where $\bm{\theta}$ denote the parameters of the model. Let $\{\bm{\theta}_t,~t=1,\ldots,T\}$ denote a Gibbs sample of $\bm{\theta}$ from $p(\bm{\theta\mid \bm{y}})$, using Equation (ref), a Monte Carlo estimate of the CPO can be obtained as:
where $T$ is the total number of Monte Carlo iterations. Based on $\widehat{\mbox{CPO}_i}$, the LPML can be estimated as:
A larger LPML value indicates better model fit.
Similar to ma2019bayesianspatial, we use the Rand index rand1971objective to evaluate the clustering performance, i.e., whether the final inferred clusters align well with the truth. Consider two partitions of $\{1,2,\ldots,n\}$, denoted as $\mathcal{C}_1 = \{A_1,\ldots,A_r\}$ and $\mathcal{C}_2 = \{ B_1,\ldots, B_s\}$. Out of all $n \choose 2$ pairs of observations, denote:
With the above specifications, the RI is calculated as
It can be seen that the RI ranges from 0 to 1, with a larger value suggesting better concordance between two clustering partitions. Computation of the RI is done using the R package fossil Rpkg:fossil.
We study the estimation performance as well as the clustering performance in this section. Two designs of true cluster configuration of Georgia counties are considered. The first case is similar to in ma2019bayesianspatial, where there are, respectively, 51, 49, and 59 counties in each cluster. The second is less balanced with 26, 44, and 89 counties in each cluster. The two partition schemes used in designing the simulation study are visualized in Figure (ref).
We consider the following data generation model:
where $\bm{w}$ is the distance- and auxiliary covariates-dependent vector of spatial random effects such that $$\bm{w}\sim \mbox{MVN}(\bm{0}, 0.25(0.81\bm{I} + 0.04\exp(-\mbox{GCD} / 4) + 0.05\bm{W}(\bm{Z}_1) + 0.1\bm{W}(\bm{Z}_2)),$$ with $\bm{Z}_1$ and $\bm{Z}_2$ being the two auxiliary covariates.
The two similarity matrices $\bm{W}(\bm{Z}_1)$ and $\bm{W}(\bm{Z}_2)$ are constructed using (ref) with $\kappa_1 =5$ and $\kappa_2 = 3$. The true parameters for the similarity matrices are set to relatively small values compared to $\alpha_0$ following zou2017covariance. For both partition schemes shown in Figure (ref), the true parameter vector for cluster 1 is set to $(4,1,-2)$, for cluster 2 $(1, 1, 0)$, and for cluster 3 $(1, -2, -1)$. For each partition shown in Figure (ref), a total of 100 datasets are generated.
In addition to the proposed model, to verify that identifying clusters do help with better estimation of the underlying coefficients, two additional models are fitted. The first alternative model is a Bayesian regression model either without clusters or spatial random effects, but includes the auxiliary covariates as main effects. It can be written hierarchically as
where $\widetilde{\bm{X}}(s_i)$ in this case becomes $\widetilde{\bm{X}}(s_i) = \left(X_1(s_i), X_2(s_i), X_3(s_i), Z_1(s_i), Z_2(s_i)\right)^\top$, and $\bm{\beta}\in\mathbb{R}^5$. We set $\sigma_\beta^2=100$ to induce a non-informative prior for $\bm{\beta}$.
The second alternative model is the Bayesian mixed model with spatial random effects but without clustering, which can be written as
Again, the parameter $\sigma_\beta^2$ is set to 100 to make a noninformative prior for $\bm{\beta}$.
The proposed approach and the two alternative models are evaluated in terms of parameter estimation. For estimation of the vector of coefficients, $\bm{\beta}$, we employ the mean absolute bias (MAB), mean standard deviation (MSD), mean of mean squared error (MMSE), and mean coverage rate (MCR) for assessment:
where $\widehat{\beta}_{\ell m r}$ denotes the posterior estimate for the $m$th coefficient of county $\ell$ in the $r$th replicate, $\overline{\widehat{\beta}}_{\ell m} = \frac{1}{100}\sum_{r=1}^{100} \widehat{\beta}_{\ell m r}$ , $\beta_{\ell m}$ is the true underlying parameter value, $\mbox{HPD}_{{\widehat{\beta}_{\ell m r }}}$ is the 95% highest posterior density interval for $\beta_{\ell m}$ in the $r$th replicate, and $1(\cdot)$ denotes the indicator function. Also, note that for the first alternative model, as we are primarily interested in estimation of the three true main effects, we omit the performance measures for the coefficients for the two auxiliary variables.
For each replicate, we set the chain length to 25,000 with thinning interval 2. The first 9,500 of retained samples are discarded as burn-in, and we use the remaining 3,000 iterations for posterior inference. The final cluster belonging inferred for each county is taken as the first mode of the posterior samples for $z_i$, $i=1,\ldots, 159.$
First we check the estimation performance using the four performance measures defined above. For ease of reference, we name the three competitive models as Alternative 1, Alternative 2, and proposed. It can be seen from the first row that with the incorporation of different clusters, each cluster of locations are allowed to have their own parameter vector. This additional flexibility of the proposed model enables less biased parameter estimation. As the proposed model includes clustering process, the chains for each parameter may jump between several underlying clusters, which causes their MSD to be larger than those for Alternatives 1 and 2, which restrict that all locations have the same set of parameters. However, with improved MAB, parameter estimates produced by the proposed model still have smaller MMSE than the other two models. Finally, as Alternatives 1 and 2 do not allow for clusters of coefficients, their parameter estimates are essentially close to the average of parameters over the 159 locations, which leads to their very low MCR.
The clustering performance of the proposed approach is presented in Figure (ref). Comparing across the two panels, it can be seen that under Design 1 there are more replicates where $K$ is correctly inferred, while under Design 2 there are more under-clustering replicates, which is due to its class imbalance. The average Rand index (ARI)'s turned out to be 0.703 and 0.752 for the two cases, respectively.
Finally, to verify that LPML is capable of reflecting the degree of fitness of the model to the data, for each simulation replicate, the LPML values of the three models are calculated. A boxplot of the 100 LPML values for each model under Designs 1 and 2 is given in Figure (ref). As discussed before, larger LPML values indicate better model fit. As clearly seen in the figure, the proposed model has overall much larger LPML values than the two alternatives, indicating that LPML is indeed capable of identifying a more suitable model in the scope of the research problem considered here.
The Georgia monthly housing dataset can be accessed at https://github.com/ys-xue/Bayesian-clustered-coefficients-regression -ACAC in .csv format. The original data source is www.healthanalytics.gatech.edu, which contains visualizations of data concerning multiple dimensions of Georgia. For each of the 159 counties, the median monthly housing cost for occupied housing units is observed. In addition, several independent variables are available: the unemployment percentage for adults between 18 and 64 years of age ($X_1$), the average per individual real and personal property taxes ($X_2$), the median home market value in thousand dollars ($X_3$), the White race population percentage ($Z_1$), the median age ($Z_2$), and population size in thousands ($Z_3$).
In our analysis, the first three economy-related covariates, $X_1$, $X_2$ and $X_3$, are used in the spatial regression part, while the remaining three demographic covariates are used in constructing the covariance matrix of spatial random effects. The final model is written as, for $i=1,\cdots, 159$,
where the $(k,k')$-th element of $W(\bm{Z}_j)$ for $(j = 1, 2, 3)$ is $\exp(- \kappa_j |Z_j(s_k) - Z_j(s_k')|)$, respectively, while for $W(\bm{Z}_4)$, the entry is $\text{exp}(-\kappa_4 \cdot \mbox{GCD})$. The priors of the unknown parameters are assigned as mentioned in Section (ref). Similar to in the simulation study, after burning in the first 9,500 of 12,500 iterations, 3,000 MCMC samples are collected the parameters. Similar to in the simulation studies, the final cluster configuration is obtained as the first mode from the posterior samples in the chains corresponding to $z_1,\ldots,z_{159}$.
We firstly apply the LPML to select the most suitable covariance structure of spatial effects $\bm{\Sigma}_{\bm{w}}$ for the model. The LPML values of the proposed auxiliary covariates assistant covariance matrix, the unity scheme, the exponential scheme and the Gaussian scheme are shown in Table (ref). Comparison of the LPML values leads to the conclusion that the proposed auxiliary covariates assisted covariance matrix provides the most suitable approximation for the covariance structure of the spatial random effects for this dataset, as it has the largest LPML value among the candidate covariance structures. Therefore the auxiliary covariates assisted covariance matrix is used in all subsequent analyses. The two alternative models we considered in the simulation studies are also examined, and their LPML values are also included in Table (ref). Among the candidate models considered, the proposed model that employs the ACAC has the largest LPML model, indicating that it is the most suitable choice to capture the heterogeneity in the Georgia housing cost data.
Three clusters of the coefficients in the spatial regression part $\left\{\beta_\ell(\bm{s})\right\}_{\ell=1}^p$ are identified through the MFM approach, whose posterior estimates are shown in Table (ref) and the cluster belongings of the 159 counties are visualized in Figure (ref). In addition, the traceplot for the number of clusters, $k$, is included in the supplemental material to verify convergence of the results. Convergence is further verified with Dahl's method dahl2006model in Section S2 of the supplemental material. Cluster 1 includes 10 counties and cluster 3 consists of 5 counties, while the rest 144 counties all fall within cluster 2. Taking a closer look, cluster 1 consists of Fulton, Douglas, Paulding, Henry, Newton, Barrow, Chattahoochee Lee, Effingham and Liberty, which are all relatively economically developed counties in terms of per capita income (among the top 50 according to 2015 United States Census Data and the 2006-2010 American Community Survey 5-Year Estimates) except Liberty. Cluster 3 consists of Fannin, Union, Towns, Rabun, and Clay. Both clusters include neighboring counties and non-adjacent counties, which again echos the finding in our simulation study that the proposed method takes into consideration both spatial adjacency and the inherent similarity between covariates that influence the spatial random effects.
From Table (ref), for counties belonging to cluster 2, the percentage of unemployment and the median house market price can help explain the change of median monthly housing cost. However, for the other counties, neither the factors we selected has impact on the dependent variable. Also, the intercept term for cluster 2 is noticeably negative, indicating a difference in the overall level of housing cost between counties in cluster 2 and those in the other two clusters.
Table (ref) shows the posterior estimates of the overall variance term of spatial random effects, $\sigma^2$, and coefficients for the similarity matrices. By comparing the posterior estimates of $\alpha_j (j = 0, \cdots, 4)$ in the auxiliary covariates assisted covariance matrix of the spatial effects, we can see that the similarity matrices defined by the size of population and the percentage of White race population have greater impact on the covariance matrix of spatial effects.
In this paper, we propose a Bayesian clustered coefficients regression model with auxiliary covariates assistant random effects. Our proposed model has two practical merits. First, our model simultaneously estimates the number of clusters and clustering configurations of regression coefficients. Second, auxiliary covariates information are included in our random effects model. The usage of proposed method is illustrated in simulation studies, where it shows accurate estimation and clustering performance. For Georgia housing cost data, our method dominates the other benchmark methods in terms of LPML.
In addition, three topics beyond the scope of this paper are worth further investigation. First, in our real data application, auxiliary covariates are selected based on their natures, which is not always available or clearly categorized in all possible applications. Proposing a quantitative criterion for auxiliary covariates determination is an interesting future work. Furthermore, different clusters may have different sparsity patterns of the covariates. Incorporating different sparsity structure of regression coefficients into the model will enable selection and identification of most important covariates. Finally, considering geographical information for clustering detection hu2020bayesian,zhao2020bayesian,geng2020bayesian is also devoted to future research.