Selection Criterion of Working Correlation Structure for Spatially Correlated Data

Abstract To obtain regression parameter estimates in generalized estimation equation modeling, whether in longitudinal or spatially correlated data, it is necessary to specify the structure of the working correlation matrix. The regression parameter estimates can be affected by the choice of this matrix. Within spatial statistics, the correlation matrix also influences how spatial variability is modeled. Therefore, this study proposes a new method for selecting a working matrix, based on conditioning the variance-covariance matrix naive. The method performance is evaluated by an extensive simulation study, using the marginal distributions of normal, Poisson, and gamma for spatially correlated data. The correlation structure specification is based on semivariogram models, using the Wendland, Matérn, and spherical model families. The results reveal that regarding the hit rates of the true spatial correlation structure of simulated data, the proposed criterion resulted in better performance than competing criteria: quasi-likelihood under the independence model criterion QIC, correlation information criterion CIC, and the Rotnizky–Jewell criterion RJC. The application of an appropriate spatial correlation structure selection was shown using the first-semester average rainfall data of 2021 in the state of Pernambuco, Brazil.


Introduction
The Generalized Estimating Equations (GEE) approach has been sparsely explored in the context of spatially correlated data analysis. This approach first requires specifying a correlation matrix, commonly called the working correlation matrix. This matrix is unknown and represented by R(ϕ), where ϕ is an unknown parametric correlation vector that completely characterizes matrix R. GEE's strength is to generate consistent estimates for regression parameters, even if the working matrix is not specified correctly. However, specifying the working matrix close to the true correlation matrix that generated the data can improve the efficiency of the parameter estimators, especially for small samples.
In this aspect, the literature has several criteria to help choose the best working correlation matrix. Among the most popular are the Rotnizky-Jewell criterion (RJC), quasi-likelihood under the independence model criterion (QIC), and Correlation Information Criterion (CIC). The RJC was initially proposed by Rotnitzky and Jewell (1990), who under the hypothesis of the true correlation matrix evaluated the asymptotic behavior of the Wald test statistic. Based on these ideas, Hin, Carey, and Wang (2007) compared the RJC criterion with the QIC when selecting the working correlation structure, using simulations with the Gaussian or binomial response. The last criterion (QIC) was presented by Pan (2001). For its development, the author used the well-known Akaike's information criterion for model selection, proposing a modification to accommodate model selection under GEE. Despite being created for selecting models in GEE this criterion has been widely used for choosing the working correlation matrix. Therefore, Hin and Wang (2009) suggested using half of the second term of the QIC, producing the CIC. Justifying it as the first term of the QIC, as it uses only the quasi-likelihood, and does not influence the information in the working correlation matrix.
These criteria, as well as some less-popular ones, were established in the context of longitudinal data analysis. Furthermore, little is known about the relative performance of each criterion regarding the selection of an optimal working correlation structure to improve the efficiency of the regression estimators. When data is georeferenced and spatially correlated with a correlation structure generated from some correlation or a semivariance function, nothing is known about the performance of these criteria.
When the specification of the response variable covariance matrix is incorrect, Zhao, Prentice, and Self (1992) provides important results about the loss of efficiency of the regression parameters estimates under GEE. Sutradhar and Das (1999) examined the loss of regression estimator efficiency due to incorrect specification of the correlation structure, based on the original work of Liang and Zeger (1986) with covariates at cluster levels. Sutradhar and Das (2000) studied the problem by performing a computational evaluation, based on an independent correlation matrix. Thus, the authors compared the estimators under independence and an AR(1) correlation structure. They concluded that estimator efficiency is strongly linked to the correctness of the working correlation matrix.
In the context of spatially correlated data analysis, the GEE approach starts by defining a vector Z = (Z(s 1 ) = z 1 , . . . , Z(s n ) = z n ) of observed response variables in various locations s 1 , . . . , s n , with s i ∈ S ⊂ R 2 . Then, from the proposal of Liang and Zeger (1986), the spatial covariance matrix is defined by V = A 1/2 R(ϕ)A 1/2 , provided that a consistent estimator is assumed for the spatial correlation parameter ϕ and A = diag(c(σ 2 1 , . . . , σ 2 n )) is a matrix of data variances and, where diag(.) is a function that takes the n-element vector and outputs an n × n diagonal matrix with the vector being the diagonal of the corresponding matrix. Without loss of generality, we assume c = 1. Thus, the quasi-likelihood estimates for the parameters of the marginal regression model are from a modified Fisher score algorithm. Hence, given a consistent estimator for the parameters of the correlation structure, a consistent estimator for the marginal regression model parameters can be consulted in the works of Prentice (1988), Albert and McShane (1995), Venezuela, Botter, and Sandoval (2007), and Nava et al. (2017).
Therefore, whatever the working correlation matrix, we havê V rob −V N ≥ 0, that is, the difference between the robust and naive estimators V rob and V N , respectively, results in a nonnegative definite matrix. Therefore, a poorly specified work matrix can result in a significant loss of efficiency. However, this problem can be controlled by conditioning a naive variancecovariance matrix.
The conditioning or matrix condition number is related to the concept of numerical efficiency, especially in problems involving matrix inversion. For this study, we are motivated by the fact that the robust estimator depends on the inversion of terms equivalent to estimatorV N . Here, ifV N is poorly conditioned, thenV rob tends to be poorly conditioned. Nonetheless, regarding the numerical efficiency of the regression model estimator, the iterative process depends on the matrix inversion repeatedly until convergence. Therefore, bad conditioning of the matrixV N can harm the numerical efficiency of the parameter estimator and even dilate it using rounding errors and numerical approximation.
Thus, in this study, we propose a new criterion to assist in the selection of the spatial working correlation matrix in GEE, based on the condition number of the variance-covariance matrix naive. This proposal is directed to the analysis of spatially correlated data. The criterion is constructed in a manner similar to the RJC, with the difference that instead of the trace of the estimated variance-covariance matrices, we introduce the concepts of norm and matrix conditioning. To assess the performance of the proposed criterion, we conducted an extensive simulation study and compared the results with three well-known criteria in the literature. The originality of this work lies in the fact that no criterion in the literature can help choose the spatial working correlation matrix, that is, a correlation matrix generated from some semivariance function under a geostatistical model, especially with the concept of matrix conditioning.
The remainder of this article is organized as follows: Section 2 summarizes the quasi-likelihood models for spatially correlated data from the results of the Wedderburn (1974), Liang and Zeger (1986), and Zeger and Liang (1986). Additionally, this section exposes the selection technique of the correlation structure proposed by Rotnitzky and Jewell (1990), Pan (2001), and Hin and Wang (2009) and describes the method of estimation of the correlation parameters. The description of the proposed criterion is provided in Section 3. Section 4 describes the simulation study conducted to assess the performance of the criteria and reports the results obtained. In Section 5, we illustrate the applicability of the proposed criterion and its competitors with a real dataset. Finally, Section 6 presents the final considerations.

Response-Correlated Quasi-Likelihood Model
Let Z = (Z 1 (s 1 ), . . . , Z n (s n )) be a vector of observed response variables at locations s 1 , . . . , s n belonging to the space S ⊂ R 2 . For simplicity, the notation Z(s i ) = z i , i = 1, . . . , n is used to indicate the response variable observed in the location s i . Then, using a generalized estimating equation approach to model the marginal structure of the mean, it is assumed that the marginal distribution of z i belongs to the exponential family, that is is defined as a monotonous, differentiable, and invertible function, such that . , x p ) and β is a vector of unknown parameters and order p×1, and X a planning matrix associated with the covariates. The covariance matrix of the response is specified as provided that R is an n × n square matrix depends on some qvector of correlation parameters ϕ between the true correlation matrix of z. Therefore, the score function and information matrix in GEE , respectively, are defined by is the matrix of the variance function, where diag(.) is a function that takes an n-element vector and outputs an n × n diagonal matrix with the vector being the diagonal of the matrix. An estimator for the regression parameters which is consistent and asymptotically normal, results from a modification of the Fisher score algorithm, given a consistent estimator of ϕ is available in the literature. For example, it can be found in the works of Prentice (1988); Venezuela, Botter, and Sandoval (2007), and Nava et al. (2017), which is given by Equation (2).
In the original article, Liang and Zeger (1986) advocates that if the working correlation matrix is correctly specified, then an estimator for the variance-covariance matrix, associated with the estimates of the regression parameters is consistent, being obtained byV However, in practice, the true correlation matrix is unknown, which, therefore, must be estimated. Here, a consistent estimator for the variance-covariance matrix is given bŷ where cov(z) = z −μ z −μ . The (3) and (4) estimators are named naive and robust, respectively.

Estimation of Spatial Correlation Parameters
Let ϕ = (φ 1 , φ 2 , φ 3 ) with φ 1 , φ 2 , and φ 3 denote the nugget effect, process contribution or variance, and scale parameter, respectively (with q = 3). The parameter φ 1 is seen as microscale variations, that is, with possible measurement errors. However, φ 2 measures the process variance, whereas φ 3 controls the decay of the correlation curve. Therefore, considering the variance-covariance matrix defined by Equation (1), the spatial correlation matrix R, of order n × n, is defined by Let h ij be the Euclidean distance between the spatial locations s i and s j , i, j = 1, . . . , n, then the (i, j)th element of matrix R, with i = j is given by where ρ : R 2 −→ R defines a spatial correlation function. For example, for the Wendland family of covariance functions, we have where B(·) denotes the Beta function. We consider κ ∈ Z + , which is called the smoothing parameter. The parameter μ characterizes the concavity of the curve and by relation (6), defines the positivity of the generated correlation matrix.
Further details about the covariance function given in Equation (5) and the condition given in Equation (6)

Preliminary Criteria for Correlation Matrix Selection
There are several criteria for selecting the best working correlation matrix in the context of GEE estimates available in the literature. However, the most popular measures are the Quasi-likelihood methods under the Independence model criterion (QIC), Correlation Information Criterion (CIC), and the Rotnizky-Jewell criterion (RJC), initially proposed by Pan (2001), Hin and Wang (2009), Rotnitzky and Jewell (1990), and Hin, Carey, and Wang (2007), respectively.
The QIC criterion is given by where Q(·) denotes the quasi-likelihood function,V I the variance-covariance matrix for the estimators of the parameters under independence, andV rob is the variance-covariance matrix estimated from the robust estimator (4). The CIC matrix is given in Equation (7), and is obtained from the middle of the second term of the QIC criterion.
The RJC is given by (3), and p is the number of covariates involved in the marginal regression model. All these measures can be used as selection criteria for the working correlation matrix.

Proposed Criterion for Correlation Matrix Selection
The RJC criterion considers that under the correct specification of the working matrix, the robust and naive estimators are close. Thus, the total variance is considered a measure of adjustment to assist in correlation selection. Therefore, there is a need for a criterion defined from an adjustment measure that carries all the information present in the variance-covariance matrices and the spatial covariance matrix.
Therefore, the criteria for selecting the correlation matrix depend on the robust variance estimator, which in turn depends on the matrix inversion process. Then, it is important to consider how rounding errors propagate, that is, to evaluate how small variations in the elements of V and V −1 influence the results of the selection statistic. Thus, we propose a criterion based on a modification in the expressions of the RJC criterion, introducing the norm and condition number of the matriceŝ V rob andV N , respectively. After this change, all the information brought from the spatial correlation matrix can be captured to estimate the variance-covariance matrices.
The introduction of the condition number of a matrix helps evaluate the sensitivity of the variance-covariance matrix. The poor conditioning of a matrix is related to matrix singularity and its definition is easily found in the main texts on numerical analysis as being where · : R n −→ R defines a matrix norm, here the Euclidean norm. Therefore, the magnitude of the matrix condition number B reveals the sensitivity of the problem solution to perturbations in the working matrix. See for instance Belsley, Kuh, and Welsch (2005) and Chapra and Canale (2016). Belsley, Kuh, and Welsch (2005) used the condition number to help in the diagnostics of collinearity between the regression variables. Chapra and Canale (2016) used the condition number in the general context of numerical analysis.
Furthermore, it stands out that the robust variance-covariance estimator assumes asymptotic properties; therefore, it is not always valid. However, if the marginal regression model and the working correlation matrix are correctly specified, the robust estimator reduces to the naive estimator. Thus, applying the norm on both sides of the iterative equation (2) to estimates of the regression parameters, we have is a solution of the GEE system, we can say thatβ Therefore, after some algebraic manipulation, we can show that Thus, we have an upper bound for the relative error of the parameter estimates of the marginal regression model under GEE. Therefore, the relative error is proportional to the inverse of the norm of the matrix estimated for variances-covariances, under the naive estimator. Here, the condition number is a proportionality factor, that is, it quantifies the sensitivity for estimates of the regression parameters. Consequently, if the norm of the estimated variance-covariance matrix is considerable, a small perturbation in the working matrix or in the response can mean a large relative error for estimates of the β's.
Thus, we replace the matrices obtained through the robust and naive estimators from their norm and condition number, respectively, in the expressions of the RJC criterion. Then, we propose a Condition-Number-Based Criterion (CNBC) given in Equation (8), which chooses the correlation structure that corresponds to the lowest value of where /p and p is the number of covariates involved in the marginal regression model.
The performance of the proposed CNBC criterion is evaluated and compared with the existing criteria using a simulation study in various scenarios of the Geostatistical framework. The main results are presented in the next section. We also evaluated the CNBC criterion in the context of longitudinal data and some results are given in the supplementary material.

Simulation Study
To evaluate the performance of the proposed method, a simulation study is conducted in terms of the selection proportion of a true spatial correlation structure compared with the criteria QIC, CIC, and RJC. The idea is to evaluate the percentage in which the criteria identified the true correlation matrix that generated data variability. To generate data, we considered three types of distributions, Normal, Poisson, and Gamma, for spatially correlated data.
To generate the responses that have a spatial and correlated structure, we use the mvrnorm() function from the MASS version 18.144 package (Venables and Ripley 2002) of the software R Core Team (2020), which generates multivariate normal distribution values with the specification of a correlation structure. In the case of Poisson and Gamma responses, we use the NORmal To Anything (NORTA) algorithm initially proposed by Cario and Nelson (1997). Further details about this algorithm can be found in Niaki and Abbasi (2008), Yahav and Shmueli (2012), and Bonat (2017).
To define the regression model, we consider the identity linkage functions g(μ i ) = μ i , logarithmic g(μ i ) = log(μ i ), (log stands for the natural logarithmic), and inverse g(μ i ) = 1 μ i for Normal, Poisson, and Gamma distributions, respectively. The regression parameter (β) is considered fixed and equal to β 0 = 0.5. For the spatial correlation structure, we considered Wendland with κ = 1 and 3, Matérn with κ = 0.5, 3.5, and spherical. The parameter vector that characterizes the spatial structure was considered fixed for all scenarios, equal to ϕ = (φ 1 , φ 2 , φ 3 ) = (0.2, 0.6, 0.9) . Three types of grids are considered. First, we consider a grid of sizes 10 × 10 and 20 × 20, both with a spacing of 1 unit between the points. Subsequently, the grid of a set of actual data on the first-semester rainfall of 2021 in the state of Pernambuco, Brazil, was collected from 157 stations and measurements.
Replicas of 1000 Monte Carlo are generated for each scenario; and in each replica, a vector of responses z = (z 1 , . . . , z n ) is obtained. Therefore, the responses and model fitting of the regression model under GEE are performed in three situations 1. The true spatial correlation matrix is defined as Wκ: Wendland covariance with smoothing κ, which generates the responses. Then, it fits the model under GEE considering the true correlation matrix. Afterward, it fits the model specifying a correlation structure Mκ: Covariance of the Matérn family with smoothing κ, and with the spherical model;  Numbers in bold report a comparison of each criterion with our own performance. It means that the criteria presented a higher percentage of correct answers for the choice of the correlation structure than in relation to the percentages referring to the choice of the non-true correlation structure.
2. The true spatial correlation matrix is defined as Mκ: covariance of the Matérn family with smoothing κ, and generates the responses. Then, it fits the model under GEE considering the true correlation matrix. Then, it fits the model specifying a correlation structure Wκ: Wendland family covariance with smoothing κ and with the spherical model; 3. Define the true spatial correlation matrix as being from the spherical model; and generate the response variables for each of the distributions, normal, Poisson, and gamma. Then, it fits the model under GEE considering the true correlation matrix. Subsequently, it fits the model by specifying a correlation structure Wκ: Wendland family covariance with κ smoothing, and with Mκ: Matérn family covariance with κ smoothing.
Thus, for each case above, the values of the QIC, CIC, RJC, and CNBC criteria are calculated. Then, the percentage of correct selection is obtained. The results can be seen in Tables 1-3.
In Table 1 are the results of the simulations for the spatial correlation structure of Normal and Poisson correlated responses. For normally distributed and spatially correlated responses, we noticed that when the true correlation structure is Wendland, the QIC criterion did not present a good result in all scenarios. However, the RJC and CNBC criteria achieved excellent results, in which they obtained a percentage of 100% in all scenarios. Additionally, although the CIC criteria had reasonable results, their percentage of correct selection was less than 50%. When the answers are generated from the Matérn family, the CIC results remain slightly equivalent to the previous ones. However, QIC regained its performance and managed to be highly com-petitive with the RJC and CNBC criteria. Here, the QIC obtained a hit rate with a variation between 79.7% and 99.8%. When we consider a correlation structure from the spherical model, the QIC criterion again showed poor results. In contrast, the RJC and CNBC criteria are excellent for selecting a true spatial correlation structure.
In the case of Poisson correlated responses, the CIC criterion performed reasonably well. However, the QIC criterion was competitive in relation to RJC and CNBC for both the Wendland and Matérn structures. The performance of the QIC stands out, for the Matérn structure and smoothing parameter κ = 3.5, in which it achieved a percentage equal to 100%, surpassing the CNBC, which reached an equal rate at 98.7%. However, for the spherical model, QIC remained with lower results. In Table 2, with the correlated response of the Gamma, the CIC and QIC criteria maintained reasonable performance for both the Wendland and Matérn covariance structures and spherical models, although the QIC showed poor performance for the spherical model considering a grid of size 10 × 10. The RJC and CNBC criteria continued to show excellent results, with percentages above 83% and 93%, respectively.
When the grid used was based on the grid used to collect real data, the criteria maintained the same performance standards as for a simulated and regular grid, as shown in Table 3. The proposed CNBC criterion showed excellent results in all scenarios, with a percentage above 93%, reaching 100% in some situations. The RJC criteria achieved a percentage greater than 80%. Therefore, the proposed criteria CNBC and the RJC are the best to choose for the spatial correlation matrix. Furthermore, note that the CNBC criterion considers the effects of rounding errors involved in the matrix inversion process, which are present in the robust and naive estimators. However, note that the proposal assumes that the marginal average model is correctly specified. Therefore, it is interesting to perform future investigations to Numbers in bold report a comparison of each criterion with our own performance. It means that the criteria presented a higher percentage of correct answers for the choice of the correlation structure than in relation to the percentages referring to the choice of the nontrue correlation structure.
evaluate the criterion performance when selecting covariates in the linear predictor, under the assumption of the correct specification of the spatial correlation matrix. Regarding the QIC criterion, its limitations in selecting the true correlation structure had already been highlighted in the works of Hin, Carey, and Wang (2007), Hin and Wang (2009), Gosho, Hamada, andYoshimura (2011), Jaman et al. (2016), and Silva and Cirillo (2018). We conducted a simulation study to evaluate the performance of the CNBC criterion compared with QIC, CIC, and RJC in the context of longitudinal data for three correlation structures, exchangeable, autoregressive, and unstructured. The results show that the CNBC criterion performs similarly to the competitors and we present some of the simulation results in the supplementary material.

Application
As an example, precipitation data provided by the Agrometeorological Monitoring System (Agritempo), available at https:// www.agritempo.gov.br/agritempo/index.jsp, are considered to assess the applicability of the criterion proposed here. The response variable is the average precipitation for the first half of 2021 in the state of Pernambuco, Brazil, collected from 157 georeferenced stations. The marginal model with covariates is given by where β 0 , β 1 , β 2 , β 3 , and β 4 denote the regression parameter unknowns that are to be estimated. The explanatory variables x(s i ), y(s j ) and H, i, j = 1, . . . , n correspond to the points of Numbers in bold report a comparison of each criterion with our own performance. It means that the criteria presented a higher percentage of correct answers for the choice of the correlation structure than in relation to the percentages referring to the choice of the nontrue correlation structure.
latitude, longitude, and altitude, respectively, for spatial location (s i , s j ) belongs to S ⊂ R 2 . To model average precipitation data for the first half of 2021 in the state of Pernambuco, Brazil, we initially defined the spatial correlation structure. To do so, we first consider the parametric vector ϕ = (φ 1 , φ 2 , φ 3 ) , with entries defined in this order, the nugget effect, process variance, and scale. In this study, we consider the correlation functions of the Wendland and Matérn families, and spherical geostatistical models. The estimation of the parametric vector is performed by generalized least squares, with results given in Table 4.
The RJC values indicate the selection of any correlation structure considered, since it gives constant values. According to the CIC and CNBC criteria, the selected correlation structure should be Wendland with κ = 3. The QIC criterion indicates the selection of the spherical correlation structure. Therefore, the final chosen model is the Wendland correlation structure with κ = 3. The estimated parameters for the regression model for all the fitted correlation structures are shown in the Table 5, with the respective robust standard and naive errors.
We note that in all cases, the choice of the working correlation matrix has a direct impact on both the robust and naive standard error estimates. Note that the standard errors were different, and that the adjustment using the Wendland structure with κ = 1 presented smaller hard standard errors of the estimates. However, the smallest standard error of the naive type of estimate was obtained when using the Matérn structure with κ = 0.5.
The results show that the naive standard errors were greater than the robust standard errors. Assuming a significance of 5%, for the Wendland case with κ = 1, only the covariate β 2 (Longitude) is marginally significant according to the Wald hypothesis test. In the case of Maérn with κ = 0.5, the covariates are not significant. In the other cases, all the covariates are significant.
Checking the half-normal probability plot with a simulated envelope in Figure 1, we observe that just in the case when considering a Matérn correlation structure with smoothing 0.5 (middle left plot), most of the points are outside the simulated envelope. Therefore, in this case, there are indications that the adjustment is inadequate. Alternatively, when considering a correlation structure of the type Wendlandκ = 1, Wendlandκ = 3, Matérnκ = 0.5, Matérnκ = 3.5 and spherical model, most points are well adjusted within the confidence band.
Furthermore, we note that in cases where there are indications of good adjustments, the presence of a discrepant point is verified. Even so, in Figure 1, when considering the correlation structure of the spherical model (bottom left plot), we noticed the presence of a greater number of points close to the limits of the confidence band. Regarding the best correlation structure, the CNBC criterion indicated the Wendland modelκ = 3, which is plausible given the fit (upper right plot) and its standard errors of the estimates.
Thus, we verified that the methodology can accurately model the spatial behavior of precipitation in the state of Pernambuco, Brazil for the period studied. It is evident that the GEE technique, with the robust proposal of selection of the spatial correlation matrix, presents itself as a promising tool in the study of the spatial distribution of precipitation.

Conclusions
It is certain that selecting an appropriate working correlation structure in GEE plays an important role in obtaining efficient estimates of the parameters of the marginal regression model. It is not only important but also advantageous to propose and understand criteria that help select a spatial correlation matrix. Therefore, in this study, a new criterion, named Condition-Number-Based Criterion (CNBC), was developed based on the condition number of the variance-covariance matrices Naive, to select a spatial working correlation structure appropriate in GEE. We compared our proposal with the RJC, QIC, and CIC criteria, available in the literature, through an extensive simulation study. To generate spatially correlated responses, the NORTA algorithm was used, with regular grids of sizes 10 × 10, 20 × 20 and an irregular grid based on the precipitation dataset. The results show that the proposed criterion performs  better than the RJC criterion, in terms of the correct selection of the correlation matrix, in all cases when considering spatially correlated Poisson and Gamma responses. In the case of the normal distribution, the criteria obtained slightly equivalent results. Compared with the QIC and CIC criteria, CNBC performed better in all scenarios. We emphasize that the proposed criterion allows the information present in the spatial correlation matrix to contribute to the estimation of all parts of the regression model under GEE, contrary to the other QIC, CIC, and RJC criteria. Furthermore, in this work, we focus on the selection of the spatial correlation matrix under the assumption that the average model is correctly specified. The specification of the correlation structure is based on semivariogram models, using the Wendland, Matérn, and spherical model families. Thus, future investigation of the CNBC criterion includes the simultaneous selection of the spatial covariance matrix and the linear predictor in GEE.
Finally, we emphasize that to choose an adequate spatial correlation structure in GEE, no method is currently available in the literature other than the one proposed in this article. The QIC and CIC criteria, reviewed in this article, did not work well when the true correlation structure is the Wendland family. Although when the true correlation structure is the Matérn family, the CIC criterion, in some cases, selected the true correlation structure, then we conclude that it is not suitable for our context, given its effective selection rate of less than 50%. However, the QIC criterion worked well in the case of the Normal response, and Poisson spatially correlated with the Matérn family correlation structure. Additionally, even though the RJC criterion is equivalent to the CNBC for normal responses, the CNBC was superior in all scenarios for the Poisson and Gamma responses, both for the Wendland and Matérn family correlation structures and the spherical model. Therefore, the results showed that the proposed criterion performs much better for spatially correlated data than the criteria popularly known in the literature. In future research, it would be interesting to see more details about the performance of the CNBC criterion in the context of longitudinal data. Some results are addressed in the supplementary material.

Supplementary Materials
In the supplementary material we present some additional simulations results on the performance of the CNBC criterion when used in the context of longitudinal data.