Bayesian Inﬂuence Analysis of the Skew-Normal Spatial Autoregression Models

: In spatial data analysis, outliers or inﬂuential observations have a considerable inﬂuence on statistical inference. This paper develops Bayesian inﬂuence analysis, including the local inﬂuence approach and case inﬂuence measures in skew-normal spatial autoregression models (SSARMs). The Bayesian local inﬂuence method is proposed to evaluate the impact of small perturbations in data, the distribution of sampling and prior. To measure the extent of different perturbations in SSARMs, the Bayes factor, the φ -divergence and the posterior mean distance are established. A Bayesian case inﬂuence measure is presented to examine the inﬂuence points in SSARMs. The potential inﬂuence points in the models are identiﬁed by Cook’s posterior mean distance and Cook’s posterior mode distance φ -divergence. The Bayesian inﬂuence analysis formulation of spatial data is given. Simulation studies and examples verify the effectiveness of the presented methodologies.


Introduction
In the fields of econometrics, air quality monitoring and epidemic monitoring, spatial data are often encountered. Spatial data are a kind of data with location attributes, which usually have spatial dependence. If the spatial correlation of data is not considered, the conclusion will often be inaccurate or even wrong. Spatial autoregressive models (SARMs) can effectively deal with the data with spatial correlation [1]. Therefore, SARMs have attracted the extensive attention of many scholars over the past few years. For example, Piribauer and Cuaresma [2] proposed two Bayesian variable selection approaches and compared their performance in SARMs. Xie et al. [3] investigated the variable selection problem in SARMs; Du et al. [4] studied the generalized method of moments estimator for partially linear additive SARMs; Li et al. [5] considered a variable selection method for SARMs based on the minimum prediction error criterion; Jay et al. [6] introduced SARMs for the statistical inference of ecological data, which discussed model selection, spatial regression, the estimation of autocorrelation, the estimation of other connectivity parameters, spatial prediction and the spatial smoothing of practical ecological inference; Anik et al. [7] investigated a Lagrange multiplier test of spatial dependence for SARMs with latent variables; Song et al. [8] proposed a class of penalized robust regression estimators based on exponential squared loss for SARMs with independent identically distributed Mathematics 2022, 10, 1306 2 of 19 errors. The above studies of SSARMs are mainly based on the hypothesis that the dependent variable follows a normal distribution. However, the normal hypothesis is often excessively restrictive of the variability of spatial data in the actual statistical inference, and unreasonable results may be generated.
In the actual spatial data, there are few cases of strict symmetry. At this time, if the spatial data is studied based on the assumption of normality, it is difficult to capture the changes in the data. In recent years, the skew-normal distribution has attracted widespread attention. A large number of studies on skew-normal distribution have appeared under various statistical models. For example, Marcos et al. [9] developed nonlinear mixed-effects models for a mixture of scales with a skew-normal distributions family and gained an estimate of the parameters by the maximum likelihood method; Yin et al. [10] introduced a variable selection procedure for a finite mixture of regression models using the skew-normal distribution; Tatsuya et al. [11] investigated the decision-theoretic properties of Stein-type shrinkage estimators in multivariate skew-normal distribution under the condition of quadratic loss; Liu et al. [12] considered an autoregressive model based on the skewnormal distribution; Mahdi [13] adopted the EM algorithm to research the maximum likelihood estimation of the mixed model with a skewed-normal assumption. However, as far as we know, there are few studies on SARMs with the response variables following skew-normal distribution.
In recent years, the Bayesian methods of local influence and case influence analysis have become widely used statistical diagnostic methods. For Bayesian local influence analysis, it is extensively studied within different models based on several objective functions. For example, Zhu et al. [14] established a universal framework of Bayesian influence analysis to evaluate the impact of simultaneous perturbations of prior, data and sample distribution for a set of statistical models; Zhang et al. [15] proposed Bayesian local influence analysis to evaluate the impact of different perturbations of individual observations, prior distribution and nonignorable missing data mechanisms in general EEs. Ouyang et al. [16] applied the Bayesian local influence method to semiparametric structural equation models and the effects of minor perturbations were evaluated using different perturbation options; Dai et al. [17] studied two Bayesian local diagnostic procedures for heteroscedastic SARMs; Ju et al. [18] developed a new SDPDMs by assuming that the random effects and error terms obey a skew-normal distribution and also developed a Bayesian local influence analysis method for it. For Bayesian case influence diagnostics, they are widely researched in various statistical models based on conditional predictive coordinates and K-L distance. Cancho et al. [19] presented Bayesian case deletion influence diagnostics for nonlinear regression models with scale mixtures of skew-normal distributions based on the K-L divergence; Zhu et al. [20] proposed three Bayesian case influence measures to identify influential points for a class of statistical models with missing data, which included the Cook's posterior mean and mode distance and φ-divergence; Tang et al. [21] detected the influence observations of generalized partial linear mixed models for longitudinal data based on φ-divergence and Cook's posterior mean distance; Hao et al. [22] investigated Bayesian case influence analysis based on the K-L divergence of a generalized autoregressive conditional heteroscedasticity model; Duan et al. [23] developed a Bayesian case deletion influence measure based on the φ-divergence of a semiparametric reproductive dispersion mixed model and presented computationally feasible formulas. So far, there are few research results regarding the Bayesian influence analysis of skew-normal spatial autoregression models (SSARMs). The importance of SSARMs is mainly reflected in the fact that SSARMs not only consider the spatial correlation of economic individuals in different regions but also consider the skewness of spatial data. It is more effective to make statistical inferences by the autoregressive model of skewed normal space than the traditional regression model. However, the observed data may deviate greatly from the established SSARMs. The motivation of this paper is that there are outliers or strong influential cases in the actual spatial data, which have an effect on the statistical analysis and inference. Therefore, based on references [14,20], a methodology of Bayesian statistical diagnosis for SSARMs is developed in this paper. The contribution of this paper is that Bayesian influence analysis, including the local influence approach and case influence measures in SSARMs, is proposed. The Bayesian local influence method is proposed to evaluate the impact of small perturbations in the data, distribution of sampling and prior. To measure the extent of different perturbations in SSARMs, the Bayes factor, the φ-divergence and the posterior mean distance are established. A Bayesian case influence measure is presented to examine the influence points in SSARMs. The potential influence points in the models are identified by Cook's posterior mean distance and Cook's posterior mode distance φ-divergence. The Bayesian influence analysis formulation of spatial data is given.
The rest of this paper is structured as follows. The Markov chain Monte Carlo (MCMC) algorithm and skew-normal spatial autoregression models are introduced in Section 2. Three Bayesian local influence methods are used to evaluate the effects of minor perturbation on the data, priors and distribution of sampling in Section 3. Section 4 shows three Bayesian case influence methods, which are used to recognize the outliers or influential points. Section 5 illustrates the effectiveness of the proposed approach through simulation studies and examples. The concluding section is presented in Section 6.

Models Introduction and MCMC Algorithm
Y = (y 1 , . . . , y n ) T is a n × 1 dimensional dependent variables vector; X is a n × p dimensional explanatory variables matrix; W is a n × n dimensional spatial weights matrix; ρ represents the strength of the spatial dependence, which can be used to measure the effect of geographic correlation on the dependent variable; β is a p × 1 dimensional regression coefficients vector, which can indicate the degree of influence of the corresponding independent variable on the dependent variable; and ε = (ε 1 , . . . , ε n ) T is a n × 1 dimensional vector of disturbances. The traditional assumption on ε is ε ∼ N n (0, σ 2 ε I n ). We assume the distribution of ε follows skew-normal, that is, ε ind.
∼ SN n 0, σ 2 ε I n , δ 2 ε I n . The advantages we assume are (1) the skewness and heavy tails can be simultaneously explained; (2) the assumption of normal distribution can be relaxed; and (3) the accurate representation of the structure present in spatial data can be provided. Therefore, SSARMs have the following expression: Similar to Arellano-Valle et al. [24], skewed-normal distribution is defined as follows.
Definition 1. If a n-dimensional random vector z obeys a n-variate skew-normal distribution, its position vector τ ∈ R n , the scale matrix ∑ is a n × n positive definite matrix and ∆ is a n × m skewness matrix, the density function is expressed as which is denoted as z ∼ SN n,m (τ, ∑ , ∆). z ∼ SN n (τ, ∑ , ∆) when n = m, which can be given a stochastic representation as z represents 'distributed as', with z 0 and z 1 being independent.
According to the Definition 1, Equation (1) is represented as follows: where δ ε represents the skewness parameters, R ε ind.
∼ N n (0, I n )I{R ε > 0}, and I represents the indicator function.

Bayesian Local Influence Analysis
This section uses the Bayesian local influence analysis method to evaluate the impact of small perturbations in the data, prior distribution and sample distribution on the posterior distribution in the SSARMs.

Perturbation Model and Manifold
Following Zhu et al. [14], the following perturbation model with skew-normal spatial data is considered: and p(Y, θ X, ω)dYdθ = 1 , where Y = (y 1 , . . . , y n ) T is a n × 1 dimensional dependent variables vector of expression (1), ω p ∈ k p represents the perturbations to the prior, ω d ∈ k d represents the perturbations to the data and ω s ∈ k s represents the perturbations to sampling distribution. Denote k = k p + k d + k s . Suppose ω 0 = (ω 0 T p , ω 0 T d , ω 0 T s ) T ∈ k represents no perturbation. Similar to Zhu et al. [14], the perturbed model = p(Y, θ X, ω) : ω ∈ k is regarded as the k-dimensional manifold under some regular conditions. The tangent space T ω of at each ω ∈ is spanned by k functions ∂ ω q (ω), where ∂ ω q = ∂/∂ω q , (ω) = log p(Y, θ|X, ω) and ω q is the qth element of ω. So, where ∂ 2 ω j ω q (ω) = ∂ 2 (ω)/∂ω j ∂ω q , the expectation of the joint probability density function p(Y, θ|X, ω) is denoted as E ω . The diagonal element g jj (ω) is considered to measure the amount of perturbation ω j and the quantity r iq = g jq (ω)/ g jj (ω)g qq (ω) is considered to measure the association between perturbations ω j and ω q . If G ω 0 is a diagonal matrix, all elements of ω are mutually orthogonal and the relevant perturbation is referred to as a suitable perturbation. If G ω 0 is not a diagonal matrix, the relevant perturbation is referred to as an improper perturbation. However, we can always find an appropriate perturbation vector ω = ω 0 + G ω 0 1 2 ω − ω 0 such that G( ω) evaluated at ω 0 equals cI k , where c is a positive number.
The Bayesian local influence analysis method is designed to quantify the degree of dependence of the posterior distribution on these three key elements of Bayesian analysis, including the prior distribution, the sampling distribution and the data. To illustrate the effect of small perturbations of the prior distribution, the sampling distribution and the data on SSARMs, the following four forms of perturbation are considered. The perturbation of covariate is considered in Example 1. The perturbation of the dependent variable is considered in Example 2. The perturbation of the prior distribution is considered in Example 3. The perturbation of the sampling distribution is considered in Example 4. Example 1. The Bayesian perturbation model for the spatial data includes many perturbation schemes. The perturbation scheme to data points is proposed for identifying outliers and influential observations. Consider the perturbation covariate X : can be spanned by when It is easily shown that G ω 0 x is a diagonal matrix, the perturbation is referred to as a suitable perturbation. Otherwise, the perturbation is deemed to be inappropriate.

Example 2.
Bayesian local influence analysis methods have many more perturbation schemes for spatial data points. In Example 1, perturbations of X are considered. However, perturbations of the dependent variable Y also have an impact on the Bayesian analysis. Consider the perturbation to the response variables Y:Y ω y = (y 1 + ω 1 , . . . , y n + ω n ) T = Y + ω y , in which ω y = (ω 1 , . . . , ω n ) T and ω 0 y = (0, . . . , 0) T stands for no perturbation. A Riemann manifold is formed by a perturbation model = p(Y, θ X, ω y ) : ω y ∈ n . The tangent space T ω y of can be spanned by It is easily shown that G ω 0 When G ω 0 y is a diagonal matrix, the perturbation is referred to as a suitable perturbation. Otherwise, the perturbation is deemed to be inappropriate. Example 3. The Bayesian perturbation model for the prior distribution includes many existing schemes. Consider the perturbation to the prior distribution of β: p(β ω β ) ∼ N ω β β, M β , in which ω β is a positive number. Additionally, ω 0 β = 1 stands for no perturbation. A Riemann manifold is formed by a perturbation model = p(Y, θ X, ω β ) : ω β ∈ 1 . The tangent space T ω β of can be spanned by It is easily shown thatG ω 0 WhenG ω 0 β is a diagonal matrix, the perturbation is referred to as a suitable perturbation. Otherwise, the perturbation is deemed to be inappropriate. Example 4. The Bayesian perturbation model for the sampling distribution includes many perturbation schemes. Consider the perturbation to the distribution of the disturbances ε: p(ε ω ε ) ∼ SN n 0, σ 2 ε I n , ω ε δ ε I n , in which ω ε is a scalar and ω 0 ε = 1 means no perturbation. A Riemann manifold is formed by a perturbation model It is easily shown that G ω 0 is a diagonal matrix, the perturbation is referred to as a suitable perturbation. Otherwise, the perturbation is deemed to be inappropriate.

Local Influence Measures
Following Zhu et al. [14], let objective function f (ω) : → r be the objective function. The influence of perturbation can be measured by considering the Bayes factor, φ-divergence and posterior mean distance. If ω(t) is a geodesic on the finite dimensional manifold with ω(0) = ω 0 and ∂ t ω(t)| t=0 = h ∈ m , in practice, the geodesic can be taken as ω 0 Firstly, when ∇ f = 0, the first-order influence measure in the defined direction h ∈ m is where G = G ω 0 , W f represents a positive semi-definite matrix. Notably, the first-order influence measure can be rewritten with appropriate perturbations ω(ω), then further defined as Following Poon and Poon [25], the first-order adjusted influence measure FIC f ( ω 0 ),h can be defined as Let η 1 ≥ . . . η t > 0 be the non-zero eigenvalues of matrix and u 1 , . . . , u t be the corresponding orthogonal eigenvectors, where u i = (u i1 , . . . , u ik ) T for i = 1, . . . , ι. The local maximum influence degree ω is reflected by the maximum eigenvalue η 1 . The corresponding eigenvector u 1 is the direction of the most significant perturbations showing the perturbations of the three components of the SSARMs to obtain the largest local variation in f ω 0 and further identify influential observations, misspecified priors and insufficient sampling distributions. Regrettably, it is not sufficient to assess u 1 local impacts by inspection alone (Poon and Poon [25]). Through their arguments, local influence is evaluated using the overall contribution vector of all eigenvectors associated with all non-zero eigenvalues: where v j is a k × 1 primary perturbation vector, the jth element is 1, the rest of the features are 0 and b jj is the jth diagonal element of matrix . Similar to Zhu et al. [14], M(0) + 2SM(0) is considered as a benchmark, where M(0) is the mean of M(0) j : j = 1, . . . , k and SM (0) is the standard error of M(0) j : j = 1, . . . , k . That is to say, if M(0) + 2SM(0) ≤ b jj , the jth observation is identified as an influential point.
Example 5 (Bayes factor). In this example, the perturbation change is measured by the Bayes factor. It measures this change mainly by the difference between log p(Y ω 0 ) without perturbation and log{p(Y|ω)} with perturbation. Consider objective function as the Bayes factor: . . , S is generated via the MCMC algorithm, unifying theGibbs sampler with the MH algorithm.
Secondly, when ∇ f = 0, the second-order Taylor expansion of the objective function f (ω) Zhu et al. [14], the second-order influence measure in the direction h ∈ m is given as In particular, when the appropriate perturbation ω(ω) is added, the second-order influence measure can be further simplified as Therefore, the second-order adjusted influence measure can be defined as in The largest eigenvalue of s is an important indicator to evaluate the local influence behavior of the model and the corresponding eigenvector is the direction of the most significant local influence, which reveals how the three elements of SSARMs are perturbed and can be identified as the potential influence points, erroneous priors and improper modelling assumptions.
Example 6 (φ-divergence). In this example, the perturbation change is measured by the φdivergence. It measures this change mainly by the difference between two posterior probability density functions. Setting the objective function f (ω) between the two posterior probability density functions before and after the introduction of perturbation ω as φ-divergence can be defined as where R(θ Y, X, ω) = p(θ|Y,X,ω) p(θ|Y,X) andφ(·) is a convex function withφ(1) = 0. It is easy to prove that ∇ φ = 0 and where a ⊗2 = aa T . H φ is difficult to calculate directly because it involves higher dimensional integrals. To solve this problem, the approximate expression of H φ can be written as Example 7 (Posterior Mean Distance). In this example, the perturbation change is measured by Cook's posterior mean distance. It measures this change mainly by the difference between two posterior means. Consider the objective function f (ω) as Cook's posterior mean distance between two posterior means of known function g(θ) before and after the introduction of perturbation ω, which is defined as where Q g (ω) = g(θ)p(θ Y, X, ω)dθ is the posterior mean of g(θ) and For the above considered objective functions and perturbation schemes, the following steps are used to achieve Bayesian local influence analysis.
Step 4. For some given objective functions (0), the jth observation is detected as an influential observation.

Bayesian Case Influence Analysis
According to the case deletion approach provided by Cook and Weisberg [26,27], three Bayesian case influence measures are considered in this section. Let Y S and X S denote subsamples of Y and X, respectively, for which all observations are in S. Let Y [S] and X [S] represent subsamples of Y and X, respectively, and delete all observations in Y S and X S .
To order to measure the impact of deleting observations in {Y S , X S } on the posterior distribution of θ, the first type of Bayesian case influence measure is the φ-influence of Y [S] , defined as where φ(·) is a convex function with φ(1) = 0, and measures the difference between p(θ Y [S] , X [S] ) and p(θ|Y, X) posterior distributions. If D φ (S) is large, the observations in S are detected to be influential.
In order to measure the impact of deleting observations in S on the posterior mean of θ, the second type of Bayesian case influence measure that Cook's posterior mean distance, denoted by CM(S), is considered as: where W θ is a positive definite matrix and can be seen as the inverse of θ's posterior covariance matrix. θ = θp(θ Y)dθ and θ [S] = θp(θ Y [S] )dθ are the posterior means of θ based on data Y and Y [S] , respectively. If the value of CM(S) is large, it points out that deleting observations in S will have a significant impact on the posterior mean. Therefore, the observations in S are considered to be strong influence points. In order to measure the impact of deleting observations in S on the posterior mode of θ, the third type of Bayesian case influence measure is Cook's posterior mode distance, denoted by CP(S), which is considered as: where, G θ is a positive definite matrix andθ = argmax θ logp(θ Y) and ) represent the posterior mode of parameter θ concerning Y and Y [S] , respectively. If the value of CP(S) is large, it indicates that deleting observations in S will have a significant impact on the posterior mode. Therefore, the observations in S are considered strong influence points.
To obtain D φ (S), CM(S) and CP(S), it is necessary to calculate p(θ|Y) and p(θ Y [S] ) .
To compute CM(S), we need to evaluate θ and θ [S] . Where θ = E θ|Y (θ) and θ [S] In order to obtain CP(S), it is necessary to calculateθ andθ [S] . Typically, the posterior mode of θ is not a closed-form, so an iterative algorithm is needed to obtainθ andθ [S] . G θ in CP(S) can be analytically obtained by evaluating J(θ) = −∂ 2 θ log p(θ Y) = −∂ 2 θ log p(Y θ) − ∂ 2 θ log p(θ) atθ. W θ can be obtained by computing the value of J θ at θ, which can also be equal to the inverse of the posterior covariance matrix of θ. Thus, D φ (S), CM(S) and CP(S) are derived directly by averaging MCMC samples over the posterior distribution p(θ|Y) .

Simulation Studies and Real Examples
The proposed approaches in Sections 3 and 4 will be illustrated by three simulation studies and practical cases of air quality from China (2020) in this section.
To illustrate the previously proposed Bayesian theory of local influence analysis, the following two simulation studies were considered.

Simulation 1.
In this simulation study, the outliers were generated in the following two ways (i) x 16 was changed to x 16 + 2 and x 29 was changed to x 29 + 2, and (ii) y 7 was changed to y 7 + 8 and y 22 was changed to y 22 + 8. For case (i), for the covariate perturbation in Example 1, the Bayesian local influence measuresFIC B,e j ,SIC D φ ,e j and SIC M d ,e j were calculated by using an MCMC algorithm iteration 10,000 times, discarding the first 5000 times and using the last 5000 iteration values to calculate them. The results are shown in Figure 1a-c, respectively. Cases 16 and 29 turned out to be influential observations, as we expected. For case (ii), for the response variables perturbation in Example 2, the Bayesian local influence measures FIC B,e j , SIC D φ ,e j and SIC M d ,e j were calculated by using an MCMC algorithm iteration 10,000 times, discarding the first 5000 times and using the last 5000 iteration values to calculate them. The results are shown in Figure 2a-c, respectively. Cases 7 and 22 turned out to be influential observations, as expected.
where ω = (ω 1 , . . . , ω 30 , ω β , ω ε ) T and X ω = X + ω x 1 T p . In this case, ω 0 = (0, . . . , 0, 1, 1) T represents no perturbation. After some calculations, we obtained G ω 0 = diag(a 11 , a 22 , a 33 ), The Bayesian local influence measures FIC B,e j , SIC D φ ,e j and SIC M d ,e j were concluded using an MCMC algorithm iteration 10,000 times, discarding the first 5000 times and using the last 5000 iteration values to calculate them based on prior distribution N 10β 0 , M β of β and the distributions SN 30 0, σ 2 ε I n , 30δ ε R ε of ε. The results are shown in Figure 3a-c. The results showed that cases 16 and 29 were identified to be influential observations, and the prior distribution of β and distributions of ε were detected to be misspecified distributions with having a large effect on FIC B,e j , SIC D φ ,e j and SIC M d ,e j .
To demonstrate the Bayesian case influence diagnostics methods, we consider the following simulation research: Simulation 3. In this simulation research, the outliers were generated in the following two ways (i) x 16 was changed to x 16 + 2 and x 29 was changed to x 29 + 2, and (ii) y 7 was changed to y 7 + 8 and y 22 was changed to y 22 + 8. The Bayesian case influence measures D φ (S),CM(S) and CP(S) were calculated using an MCMC algorithm iteration 10,000 times, discarding the first 5000 times and using the last 5000 iteration values to calculate them. The results are shown in Figures 4a-c and 5a-c. In this simulation, samples 16 and 29 are created as outliers for x and the other points as normal. Therefore, results are correct when samples 16 and 29 are detected to be outliers and the other points are less influenced or have no influence. Similarly, samples 7 and 22 are creased as outliers for y and the other points as normal. Therefore, the results are correct when samples 7 and 22 are detected to be outliers and the other points are less influenced or have no influence. It can be seen from Figure 4 that cases 16 and 29 were detected to be influential; it can be seen from Figure 5 that cases 7 and 22 were diagnosed as influence points. As we expected, the results are the same as the previous diagnosis of Bayesian local influence measures.  To demonstrate the Bayesian case influence diagnostics methods, we consider t following simulation research: Simulation 3. In this simulation research, the outliers were generated in the following two wa (i) was changed to + 2 and was changed to + 2, and (ii) was changed + 8 and was changed to

Real Example
A dataset concerning China's air quality was used to demonstrate our proposed approach in the subsection. Let X p be a 31 × 1 vector of the fine particulate matter (PM 2.5 ) in 31 provinces, let X s be a 31 × 1 vector of the sulfur dioxide (SO 2 ) in 31 provinces, let X c be a 31 × 1 vector of the carbon monoxide (CO) in 31 provinces and let Y = (y 1 , . . . , y 31 ) T be a 31 × 1 vector of the air quality index (AQI) in 31 provinces.
Firstly, the above data were exploited to illustrate the application of the Bayesian local influence analysis method. The following three perturbation options were considered:  Figure 8b show that the prior distribution of β, the distributions of ε were diagnosed to have a significant effect by SIC D φ ,e j . The results of Figure 8c showed that the prior distribution of β was diagnosed to have a significant effect by SIC M d ,e j . Secondly, the above data are used to illustrate the application of the Bayesian case influence analysis method. We calculated Bayesian case influence measures D φ (S), CM(S) and CP(S). The results are presented in Figures 9 and 10. Figure 9a-c demonstrate that cases 6 and 19 were diagnosed as influence points. Figure 10a-c demonstrate that cases 1 and 12 were diagnosed as influence points. As we expected, the results are the same as the previous diagnosis of Bayesian local influence measures. nosed as influence points by , , , and , in perturbation scheme (iii). The results of Figure 8b show that the prior distribution of , the distributions of ε were diagnosed to have a significant effect by , . The results of Figure 8c showed that the prior distribution of was diagnosed to have a significant effect by , . Figure 6. Plots of covariate perturbation in the real example.     Secondly, the above data are used to illustrate the application of the Bayesian case influence analysis method. We calculated Bayesian case influence measures ( ) , ( ) and ( ). The results are presented in Figures 9 and 10. Figure 9a-c demonstrate that cases 6 and 19 were diagnosed as influence points. Figure 10a-c demonstrate that cases 1 and 12 were diagnosed as influence points. As we expected, the results are the same as the previous diagnosis of Bayesian local influence measures.    Finally, by reviewing China's air quality data for April 2020, it was found that cases 6 and 19's influence points represent Harbin and Shenyang, respectively. In particular, Harbin's air quality is the second-lowest among the 31 provincial capitals in China, while Shenyang's air quality is the third-lowest among the 31 provincial capitals in China. This indicates that the air quality situation in these two cities is more prominent among the 31 provinces. Similarly, it was found that cases 1 and 12 of the diagnosed impact points represented Beijing and Lhasa, respectively. Among them, Beijing's air quality ranks 17th among 31 provincial capital cities. The air quality of Lhasa ranks fourth among the 31 Finally, by reviewing China's air quality data for April 2020, it was found that cases 6 and 19's influence points represent Harbin and Shenyang, respectively. In particular, Harbin's air quality is the second-lowest among the 31 provincial capitals in China, while Shenyang's air quality is the third-lowest among the 31 provincial capitals in China. This indicates that the air quality situation in these two cities is more prominent among the 31 provinces. Similarly, it was found that cases 1 and 12 of the diagnosed impact points represented Beijing and Lhasa, respectively. Among them, Beijing's air quality ranks 17th among 31 provincial capital cities. The air quality of Lhasa ranks fourth among the 31 provinces, indicating that the air quality of Lhasa is indeed outstanding among the 31 provincial capital cities. Although Beijing's air quality ranks in the middle of the 31 provincial capitals, there are many sandstorms in Beijing in April, which will lead to abnormal air quality in Beijing. It is further illustrated that the impact point diagnosed is consistent with the relevant air quality reality.

Discussion
Based on SSARMs with response variables obeying skew-normal distribution, this paper proposed a Bayesian local influence method to evaluate the impact of small perturbations in data, prior distribution and sample distribution. Perturbation models with separate or simultaneous perturbations of covariates, response variables, parameter prior and sample distributions were developed. A Bayesian perturbation manifold was constructed to describe the internal structure of perturbation models and quantify the perturbation degree of different perturbation models. The Bayesian local influence measures, including the Bayes factor, the φ-distance and the posterior mean distance, were used to evaluate the impact of various perturbations. In addition, three Bayesian case deletion influence measures, including φ-distance, Cook's posterior mean distance and Cook's posterior mode distance, were proposed to identify potential outliers in skew-normal spatial autoregression models. The effectiveness of our proposed approach was illustrated by three simulation studies and a real example. The results showed that: (1) our proposed Bayesian local influence approach can effectively identify the potential influence points, misspecified prior distribution and misspecified sampling distribution; (2) our proposed Bayesian case influence approach can be used to effectively detect the potential influence observations; (3) the outliers detected by the Bayesian local influence approach and Bayesian case influence approach are consistent, which further explains the rationality of the two methods.