{'fraction_non_alphanumeric': 0.035909787606744036, 'fraction_numerical': 0.011422523903364718, 'mean_word_length': 4.625744200369534, 'pattern_counts': {'":': 0, '<': 5, '<?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': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Spatial Homogeneity Pursuit for Survival Data with an Application to the SEER Respiration Cancer

In this work, we propose a new Bayesian spatial homogeneity pursuit method for survival data under the proportional hazards model to detect spatially clustered patterns in baseline hazard and regression coefficients. Specially, regression coefficients and baseline hazard are assumed to have spatial homogeneity pattern over space. To capture such homogeneity, we develop a geographically weighted Chinese restaurant process prior to simultaneously estimate coefficients and baseline hazards and their uncertainty measures. An efficient Markov chain Monte Carlo (MCMC) algorithm is designed for our proposed methods. Performance is evaluated using simulated data, and further applied to a real data analysis of respiration cancer in the state of Louisiana.


Introduction
Clinical data on individuals are often collected from different geographical regions and then aggregated and analyzed in public health studies. Analysis of such data conducted on a higher level often assume that covariate effects are constant over the entire spatial domain. This is a rather strong assumption, as all intrinsic heterogeneities in data are ignored. For example, if one was to study the hazard for patients with lung cancer, it is expected that the true hazard is not the same in areas where there is little air pollution and severe air pollution, even for patients with similar characteristics. From Tobler's first law of geography (Tobler, 1970), "Everything is related to everything else, but near things are more related than distant things," it is reasonable to consider similarities between nearby locations in survival data due to environmental circumstances in geographically close regions.
Existing approaches that account for such patterns in survival data can be put into two major categories. The first one is to incorporate spatial random effects in survival models, such that spatial variations are accounted for by different intercepts for different regions, while parameters for covariates are held constant. Banerjee and Dey (2005) used a frailty model to study infant mortality with two specifications for the spatial random effects -the areal conditionally autoregressive (CAR) prior, and a geostatistical spatial process. Zhou et al. (2008) proposed a joint model that uses the CAR prior to account for spatial variations in modeling both the age at diagnosis, and the subsequent conditional survival time.
The CAR prior has also been used in the accelerated failure time (AFT) model to analyze prostate cancer survival in Louisiana (Zhang and Lawson, 2011). In addition to normal spatial random effects, or log-normal frailties, gamma frailties have been used with the proportional hazards model to investigate spatial leukemia survival (Henderson et al., 2012) Another important approach, instead of assuming all covariate effects are constant, allows parameters to be spatially varying.  proposed an AFT model that allows the parameters, instead of spatial frailties, to follow certain spatial prior distributions. Similar to geographically weighted regression, where observations are assigned weights depending on their distance to the target location of interest. Hu and Huffer (2019) weighed observations with exponential kernel and proposed modified Nelson-Aalen and Kaplan-Meier estimators for spatial survival functions. Moving to the study of covariate effects, Xue et al. (2019) used the similar weighting scheme and proposed a modified partial likelihood to estimate varying coefficients in a spatial proportional hazards model.
Despite their flexibility, the aforementioned spatially varying coefficients models can be unnecessarily large. Imposing certain constraints on nearby regions so that they have the same parameter values provides an efficient way of reducing the model size without sacrificing too much of its flexibility. While similar endeavor have been made to cluster spatial survival responses (Huang et al., 2007;Bhatt and Tiwari, 2014), the clustering of covariate effects and baseline hazards have yet to be studied for survival data.
Two challenges are to be tackled for clustering of coefficients and baseline hazards for spatial survival models. First, the spatial structure needs to be appropriately incorporated into the clustering process. Contiguousness constraints should be added so that truly similar neighbors are driven to the same cluster. The constraints, however, should not be overly emphasized, as two distant regions may still share similar geographical and demographical characteristics and thus parameters. Existing methods, such as in Lee et al. (2017Lee et al. ( , 2019 and Li and Sang (2019), do not allow for globally discontiguous clusters, which is a serious limitation. Second, the true number of clusters is unknown, and needs to be estimated. With the probabilistic Bayesian framework, simultaneous estimation of the number of clusters and the clustering configuration for each region is achieved by complicated search algorithms (e.g., reversible jump MCMC, Green, 1995) in variable dimensional parameter spaces. Such algorithms assign a prior to the number of clusters that needs to be updated in every MCMC iteration, which made them difficult to implement or automate, and suffer from mixing issues as well as lack of scalability. Nonparametric Bayesian approaches, such as the Chinese restaurant process (CRP; Pitman, 1995), provide another approach to allow for uncertainties in the number of clusters. Its extension, the distance dependent CRP (ddCRP; Blei and Frazier, 2011), considers spatial information, and makes a flexible class of distributions over partitions that allows for dependencies between their elements. The CRP framework, however, has been shown to be inconsistent in its estimation of number of clusters (Miller and Harrison, 2013). Lu et al. (2018) proposed the powered CRP that suppresses the small tail clusters. Similar to the traditional CRP, however, it does not consider distance information, and therefore is not well-suited when spatial homogeneity is to be detected.
To address these challenges, in this work, we consider a spatial proportional hazards model, and propose a geographically weighted Chinese restaurant process (gwCRP) to capture the spatial homogeneity of both the regression coefficients and baseline hazards over subareas under piecewise constant hazards models framework (Friedman et al., 1982). Our main contributions in this paper are in three folds. First, we develop a new nonparametric Bayesian method for spatial clustering which combines the ideas of geographical weights and Dirichlet mixture models to leverage geographical information. Compared with existing methods, our proposed approach is able to capture both locally spatially contiguous clusters and globally discontiguous clusters. Moreover, an efficient Markov chain Monte Carlo (MCMC) algorithm is proposed for our proposed model without reversible jumps to simultaneously estimate the number of clusters and clustering configuration. Second, we introduce a model diagnostic technique, Logarithm of the Pseudo-Marginal Likelihood (LPML;Gelfand and Dey, 1994), to tune the range parameter in gwCRP. After tuning of range parameter based on LPML, the gwCRP has nearly consistent estimator of the number of clusters.
Finally, we apply the method to the analysis of the Surveillance, Epidemiology, and End Results (SEER) Program data in the state of Louisiana among different counties, which provide important information to study spatial survival rates.
The remainder of the paper is organized as follows. In Section 2, we develop a homogeneity pursuit of survival data in the piecewise constant proportional hazard framework with gwCRP prior. In Section 3, a collapsed Gibbs sampler algorithm and post MCMC inference are discussed. The extensive simulation studies are carried out in Section 4. For illustra-tion, our proposed methodology is applied to respiratory cancer survival data in Section 5.
Finally, we conclude this paper with a brief discussion in Section 6. For ease of exposition, additional technical results are given in supplementary material.

Spatial Piecewise Constant Hazards Models
Let T i denote the survival time for patient at location s i , δ i = 1 representing the event and δ i = 0 indicating censored, and X (s i ) denote the vector of covariates corresponding to T i for i = 1, 2, ..., n, and = 1, 2, ..., n i , with n i denoting the number of the patients at s i . All data that have been observed is denoted by D = {(T i , δ i , X (s i )), i = 1, 2, ..., n, = 1, 2, ..., n i }.
The semi-parametric proportional hazards model of Cox (1972) has long been used for timeto-event data. Under the Cox model, without the need of estimating the baseline hazard function, the covariate effects are estimated through maximizing the logarithm of the partial likelihood function given by where R(t) = {m : t m ≥ t} is the set of indices for subjects at risk at time t, and β is a p × 1-dimensional vector of regression coefficients. It can be seen from (1) that only the vector of regression coefficients is of interest, while the nonparametric baseline hazard function is canceled out. However, if one is interested not only in covariates effects, but also in prediction, the baseline hazard is needed. To circumvent this problem, we consider a proportional hazards model with piecewise constant baseline hazard function. Assume we partition [0, ∞) into J intervals (0 = a 0 < a 1 < · · · < a J = ∞), then the hazard function is given by with piecewise constant baseline hazard function λ 0 (t) = λ j for a j−1 ≤ t < a j , j = 1, . . . , J.
For the piecewise constant hazard function mentioned in (2), the baseline hazards λ 1 , . . . , λ J and regression coefficients β are constants over different regions. Due to observed environmental factors, spatially varying patterns in baseline hazards and regression coefficients of hazard function need to be considered. The piecewise constant hazard function with spatially varying pattern is therefore given by where λ 0(s i ) (t) = λ j (s i ) for a j−1 ≤ t < a j , j = 1, . . . , J. Under this model, λ(s i ) = (λ 1 (s i ), . . . , λ J (s i )) and β(s i ) represent the location-specific baseline hazards and regression coefficients.
Based on the location-specific survival model, the likelihood function for observed survival data D is obtained as where is the survival function for patient at location s i .
After some algebra, the logarithm of likelihood function for observed survival data D is obtained as the collection of parameters, then the maximized likelihood estimate (MLE) θ(s i ) can be obtained by solving the score function, which is the derivative of the logarithm of likelihood function in (6), and the estimated variance-covariance matrix of MLE is Σ i = (−H) −1 , respectively, where (−H) denotes the negative Hessian matrix. Based on the MLEs and estimated variance-covariance matrices, we have the following approximation of the likelihood.
is approximated as where MVN stands for the multivariate normal distribution.
For ease of exposition, the derivations of Σ i and the proof of Proposition 1 are given in Appendix A and Appendix B.
Based on the normal approximation given in Proposition 1, a natural way which follows Gelfand et al. (2003) for spatially varying pattern of baseline hazards and regression is to give a Gaussian process prior to θ(s i ), i = 1, . . . , n. The Gaussian process for θ(s i ), i = 1, . . . , n is defined as where θ = (θ(s 1 ) , . . . , θ(s n ) ) , µ is a p + J dimensional vector, H(φ) is a n × n spatial correlation matrix depending on the distance matrix with parameter φ, Σ is a (p + J) × (p + J) covariance matrix, and ⊗ denotes Kronecker product.
where |s i − s j | is the distance between s i and s j , and φ is the range parameter for spatial correlation. For the Gaussian process prior, the parameters of closer locations have stronger correlations.
For many spatial survival data, some regions will share same covariate effects or baseline hazards with their nearby regions. In addition, some regions will share similar parameters regardless of their geographical distances, due to the similarities of regions' demographical information such as income distribution (Ma et al., 2019;, food environment index, air pollution (Zhao et al., 2020), and etc.. A spatially varying pattern for θ(s i ), i = 1, . . . , n is not always valid. Based on the homogeneity pattern, we focus on the clustering of spatially-varying parameters. In our setting, we assume that the n parameter vectors can be clustered into k groups, i.e., θ(s i ) = θ z i ∈ {θ 1 , . . . , θ k }.

Geographically Weighted Chinese Restaurant Process
A latent clustering structure can be introduced to accommodate the spatial heterogeneity on parameters of sub-areas. Under the frequentist framework, the clustering problem could be solved in a two-stage approach: first obtain the estimate of number of clusters, k, then detect the optimal clustering assignment among all possible clusterings of n elements into k clusters.
However, in this approach, the performance of the estimation of cluster assignments highly relies on the estimated number of clusters, it may ignore uncertainty in the first stage and cause redundant cluster assignments. Bayesian nonparametric method is a natural remedy to simultaneously estimate the number of clusters and cluster assignments. The Chinese restaurant process (CRP; Pitman, 1995;Neal, 2000) offers choices to allow uncertainty in the number of clusters by assigning a prior distribution on (z 1 , z 2 , . . . , z n ). In CRP, z i , i = 2, . . . , n are defined through the following conditional distribution (also called a Pólya urn scheme, Blackwell et al., 1973).
|c|, at an existing cluster labeled c, α, at a new cluster.
Here |c| refers to the size of cluster labeled c, and α is the concentration parameter of the underlying Dirichlet process. Based on the Pólya urn scheme shown in (9), the costumer will have no preference for sitting with different costumers. For the spatial survival data, nearby regions will share similar environmental effects such as P.M. 2.5, water quality, etc..
These similar effects will lead the nearby sub-regions will share similar parameters. In order to consider similar effects caused by geographical distance. We modify the tradition CRP to geographically weighted CRP (gwCRP) in order to have the costumer will have higher probability sitting with their familiar costumers which are geographically nearby. We have the conditional distribution of θ(s i ) given θ(s 1 ), . . . , θ(s i−1 ) based on following definition.
where K * denote the number of clusters excluding the i-th observation, Based on the Definition 1, we have similar Pólya urn scheme called gwCRP for conditional distribution in (10) with CRP.
Proposition 2. A Pólya urn scheme of gwCRP is defined as at an existing cluster labeled c, α, at a new cluster. (11) For the geographical weights, we adapt the Stochastic Neighborhood Autoregressive (SNCAR) model (White and Ghosh, 2009) for areal data. SNCAR is an extension of the ordinary Conditional Autoregressive (Banerjee et al., 2014) model. Unlike the general adjacency matrix, whose diagonal elements are all 0 and off diagonal element a ij = 1 if areas A i and A j share a common boundary, the SNCAR model allows the off-diagonal elements to depend on unknown parameters. Compared the existing geographically weighted regression literatures, our weights are obtained by graph distance between different areas. Following Müller et al. (1987); Bhattacharyya and Bickel (2014) graph distance between two vertices v i and v j is defined as: where |V (e)| represents the cardinality of edges in e. In this way, we can calculate the graph distance among the counties in the data set. Based on the graph distance calculated by (12), we calculate the geographical weights by: where d v i v j is the graph distance between areas i and j. For the weighting function in (13), we give the largest weight (w ij ≡ 1) for the areas sharing the same boundary, which follows the first law of geography (Tobler, 1970). For the simplicity, we refer to gwCRP introduced above as gwCRP(α, h).
Remark 1. Based on the Pólya urn scheme defined in (11) and geographical weighting scheme defined in (13), we find that (i) when h = 0, the gwCRP reduces to traditional CRP, which leads to over-clustering problem in estimating of the number of clusters; (ii) when h → ∞, a new costumer just only choose the table representing spatialy contiguous regions. This will also lead to the same over-clustering problem as CRP.

Bayesian Inference
In this section, we will introduce the MCMC sampling algorithm, post MCMC inference method, and Bayesian model selection criterion.
The marginalization over k can avoid complicated reversible jump MCMC algorithms or even allocation samplers. The full conditionals of z 1 , . . . , z n are given in Proposition 3.
Proposition 3. The full conditional distributions P (z i = c | z −i , θ, θ) of z 1 , . . . , z n is given The proof of Proposition 3 is shown in C.

Inference of MCMC Results
Another important task for Bayesian nonparametric method is the inference of MCMC results. The posterior mean or median of clustering configurations z is not suitable. Dahl's method (Dahl, 2006)   with B(i, j) representing the empirical probability that observations i and j are in the same cluster among these M posterior samples. Furthermore, the iteration with the least squares distance to the Euclidean mean matrix B is obtained via and the clustering configuration z and estimated parameters from the C LS -th posterior sample are the posterior estimates selected by Dahl's method. The advantages of Dahl's method is that it utilizes information of all posterior clusterings, and selects the "average" clustering rather than forming one from an external algorithm.

Model Selection Criterion
In gwCRP, the decaying effect parameter h for geographical weights needs to be tuned, and we use the Logarithm of the Pseudo-Marginal Likelihood (LPML; Ibrahim et al., 2001) based on conditional predictive ordinate (CPO; Gelfand et al., 1992;Geisser, 1993;Gelfand and Dey, 1994) to select h. The LPML is defined as where CPO i is the i-th conditional predictive ordinate. Following Chen et al. (2000), a Monte Carlo estimate of the CPO, within the Bayesian framework, can be obtained as where B is the total number of Monte Carlo iterations, θ b z i is the b-th posterior sample from Algorithm 1, and f (·) is the likelihood function define in (6). An estimate of the LPML can subsequently be calculated as: A model with a larger LPML value is preferred.

Simulation
In this section, we will illustrate our proposed methods based on the simulation results.

Simulation Setting and Evaluation Metrics
In this section, we present simulation studies under four different designs to illustrate the performance of our proposed gwCRP method and compare with traditional CRP, in terms of both clustering configuration and estimation of regression coefficients and piecewise constant baseline hazards under proportional hazards model. Survival datasets that resemble the SEER respiratory cancer data for Louisiana are generated. The censoring rate is around 30%.
We design four different geographical clustering patterns in Louisiana state, which are shown in Figure 1. Designs I and III have three true clusters, and Designs II and IV have two true clusters. In addition, Designs II and III both have one cluster consisting of two disjoint areas since, in practice, it is still possible for two distant counties to belong to the same cluster.
Design IV has two clusters both consisting of disjoint areas. For the the matchability of clustering configurations, we use Rand Index (Rand, 1971) to measure the agreement between z (C LS ) and the true clustering configuration. The Rand where a denotes the number of pairs of objects in S that are in the same cluster in S 1 and the same cluster in S 2 , b denotes the number of pairs of objects in S that are in different clusters in S 1 and different clusters in S 2 , c denotes the number of pairs of objects in S that are in the same cluster in S 1 and different clusters in S 2 , d denotes the number of pairs of objects in S that are in different clusters in S 1 and the same cluster in S 2 . The Rand Index lies between 0 and 1, with a higher value indicating better agreement between two partitions.
Specially, when two partitions agree perfectly, the Rand Index is 1. In our simulation, the Rand Index is obtained by using R-package fossil (Vavrek, 2011).
In addition to clustering performance, we further evaluate the estimation performance of covariates coefficients and baseline hazards, which is assessed by average of mean squared error (AMSE) defined as follows. Let z = (z 1 , . . . , z n ) be the true clustering label vector, θ r (s i ) be the true parameter value of cluster j, κ r = 64 i=1 1(z i = r) be the number of counties in cluster r, r = 1, . . . , k, k r=1 κ r = n, and for simulated data set t, let θ (t) (s i ) be Dahl's method estimate at location s i . Then AMSE is calculated as which calculates mean squared errors for each cluster first, and then average across clusters. Figure 2 shows the histogram of k estimates and boxplots of Rand Index under different h and the optimal selected by LPML for four simulation designs. We see that when h = 0, the proposed gwCRP method is identical to the traditional CRP method, and in this case, CRP always tends to over-cluster and often yields smaller Rand Index than the results under h > 0.

Simulation Results
Another important trend is that, as h increases, the estimated number of clusters decreases first and then increases, and the Rand Index increases first and then decreases as h becomes too large. As we discussed in Remark 1, this is because when h increase from 0, the spatial patterns in the data is captured by the proposed gwCRP method. However, as h → ∞, the geographical weights w ij for spatial-discontiguous counties become 0, which means only adjacent counties can be classified into the same cluster, therefore leading to over-clustering phenomenon again. It is also discovered that the clustering perfomance under optimal h selected by LPML is very well, with the probability of selecting true number of clusters always greater than 0.75, and Rand Index larger than or similar to the highest results attained by some fixed value of h. h for different designs. For simplicity of summary results, here the AMSE of β is the average of AMSE of β 1 , β 2 , β 3 since they have similar scales, and the value of log λ is the average of AMSE of log λ 1 , log λ 2 , log λ 3 , respectively. Similar to the clustering performance, traditional CRP has the largest AMSE and AMSE decrease as h increase from 0 to moderate values, and increase again as h increase to relatively large values. The results of optimal h selected by LPML also has the best performance in estimation.
In a brief conclusion based on our simulation studies, gwCRP models have better performance than CRP both for clustering and parameter estimation. Our proposed model selection criterion, LPML, can nearly select the best performance h value for both clustering and parameter estimation.   Figure 3. We first fit the Cox model of patients for each county using the covariates selected. The regression coefficients are visualized in Figure 4. From Figures 3 and 4, it is seen that some counties have similar characteristics, no limited to only adjacent counties, indicating possibilities of globally discontiguous clusters.

Data Analysis
We set three pieces for the baseline hazard with cutting points 2.01 and 6.01, then run h from 0 to 10 with grid 0.1, and for each h, 2000 MCMC iterations are run and drop the first 500 as burn-in. The optimal h selected by LPML is 3.9, the corresponding estimate of number of clusters is 3, while the traditional CRP classifies the counties into 5 clusters. The plots of clustering patterns of CRP and gwCRP Optimal (h = 3.9) are shown in Figure 5, from which it is seen that the gwCRP captures the globally discontiguous clusters very well.
The estimates of regression covariates coefficients and baseline hazards obtained by gwCRP Optimal (h = 3.9) are given in Table 3, from which we see that, though the baseline hazards are similar, the regression covariates coefficients are quite different across different clusters. From the clustering results shown in Figure 5, we see that our proposed method successfully detects both spatially contiguous cluster and discontinuous cluster simultaneously. Most counties in Cluster 1 are spatially contiguous and most counties in Cluster 3 are spatially discontinuous. From the estimation results shown in Table 3, we see that (i) for the counties in Cluster 1, all four covariates have moderate hazard effects compared with other counties; (ii) for the counties in Cluster 2, diagnostic age has least hazard effects, but male, later cancer stage and distant historical stage will have higher hazards effects than other counties; (iii) for the counties in Cluster 3, the diagnostic age has most hazards effects than other counties.

Discussion
In this paper, we proposed a geographically weighted Chinese restaurant process to capture spatial homogeneity of regression coefficients and baseline hazards based on piecewise constant hazard model. An efficient MCMC algorithm is proposed for our methods without complicated reversible jump algorithm. Extensive simulation results are carried out to show that our proposed method has better clustering performance than the traditional CRP in spatial homogeneity pursuit for survival data. Simulation studies also show that our proposed methods have promising results in coefficients and baseline hazard estimation. An application to analysis of SERR data provides an interesting illustration of our proposed methods.
Furthermore, three topics beyond the scope of this paper are worth further investigation.
In this paper, our proposed algorithm is based on two step estimation under piecewise constant proportional hazard model assumption. Proposing an efficient sampling algorithm without Laplace approximation is an important future work. Furthermore, we fixed the number of pieces of baseline hazards in both simulation studies and real data analysis.
Imposing adaptive number of pieces model in baseline hazards is devoted for future research.
Finally, variable selection approaches based on hierarchical CRP (Griffiths et al., 2004) is also worth being investigated.

A Derivations of Σ i
After taking the second-order partial derivatives of the log-likelihood function, elements of the Hessian matrix are given by Let η(s i ) = log λ(s i ). Then, after this reparameterization, the (negative of) Hessian matrix is given by The estimated variance of the MLEs is Σ i = (−H) −1 .

B Proof of Proposition 1
For the likelihood L(s i ) at particular location s i , MLE θ(s i ) and its variance estimates Σ i .
We can used the second order of the Taylor And then we have L(s i ) ≈ MVN( θ(s i )|θ(s i ), Σ i ). (24)

C Full Conditionals of Collapsed Gibbs Sampling
Full contitonal of θ k :