Bayesian Analysis of Censored Spatial Data Based on a Non-Gaussian Model

TL;DR

This paper introduces a Bayesian non-Gaussian spatial model using a skew Gaussian-log Gaussian distribution to effectively analyze censored spatial data with skewness and heavy tails, employing MCMC for inference.

Contribution

It extends the skew log Gaussian model to handle censored data and heavy tails in spatial analysis, providing a new Bayesian framework with data augmentation and MCMC methods.

Findings

Effective modeling of censored spatial data with skewness and heavy tails.

Successful application to simulated and real datasets.

Demonstrates improved inference over traditional Gaussian models.

Abstract

In this paper, we suggest using a skew Gaussian-log Gaussian model for the analysis of spatial censored data from a Bayesian point of view. This approach furnishes an extension of the skew log Gaussian model to accommodate to both skewness and heavy tails and also censored data. All of the characteristics mentioned are three pervasive features of spatial data. We utilize data augmentation method and Markov chain Monte Carlo (MCMC) algorithms to do posterior calculations. The methodology is illustrated using simulated data, as well as applying it to a real data set. Keywords: Censored data, data augmentation, non-Gaussian spatial models, outlier, unified skew Gaussian.

Figures

Tables

Table 1. Table 1 . Identifiability of the correlation parameters based on the balanced case of censoring.

${\theta }_w$	True value	0.2	1	1.5
	Estimated value	0.98	1.15	1.41
${\theta }_{\lambda }$	True value	0.1	0.7	1.1
	Estimated value	0.14	0.81	0.98

Table 2. Table 2 . Prior Sensitivity Analysis: Setup of the Experiment and Maximum Relative Change (MRC).

				MRC
	Hyperparameter	Benchmark	Alternative values	$17.5\%$	$67\%$
${\beta }_0$	$c_0$	${10}^4$	${10}^2$ , ${10}^6$	0.0611	0.0926
$\alpha$	$c_1$	${10}^5$	${10}^3$ , ${10}^7$	0.132	0.1915
${\sigma }^2$	$(c_2,c_3)$	$({10}^{-6},{10}^{-6})$	$({10}^{-4},{10}^{-8})$,$({10}^{-8},{10}^{-4})$	0.1709	0.1135
${\omega }^2$	$(c_4,c_5)$	$(0.1,9)$	$(0.7,1.5)$ , $(0.5,0.7)$	0.3147	0.652
${\theta }_w$	$c_8$	$0.7$	$0.4$ , $1.4$	0.393	0.4117
${\theta }_{\lambda }$	$c_9$	$0.7$	$0.35$ , $1.60$	0.421	0.4705
$\nu$	$(c_6,c_7)$	$(0.5,1.5)$	$(0.30,1.0)$ , $(0.75,0.8)$	0.8003	0.96

Table 3. Table 3 . The estimation of the parameters in four models: GAUS, SUG, GLG and SUGLG.

			${\beta }_0$	$\alpha$	${\sigma }^2$	${\omega }^2$	${\theta }_w$	${\theta }_{\lambda }$	$\nu$
		True value	0	3	1	0.1	0.5	0.5	1
Estimated values	GAUS	$17.5\%$	1.67	-	2.17	0.9	0.38	0.15	-
		$67\%$	2.03	-	1.86	0.91	0.46	0.4	-
	SUG	$17.5\%$	1.21	2.03	1.68	0.83	0.45	0.12	-
		$67\%$	1.42	1.35	1.4	0.71	0.58	0.39	-
	GLG	$17.5\%$	0.98	-	1.74	0.66	0.22	0.34	0.51
		$67\%$	1.5	-	1.38	0.59	0.55	0.59	1.23
	SUGLG	$17.5\%$	0.05	2.77	1.33	0.2	0.56	0.45	1.12
		$67\%$	0.8	2.13	1.15	0.23	0.46	0.53	0.92

Table 4. Table 4 . The root of mean square error (RMSE) of the predicted response values based on GAUS, SUG, GLG and SUGLG models.

		GAUS	SUG	GLG	SUGLG
RMSE	$17.5\%$	0.191	0.083	0.095	0.031
	$67\%$	0.211	0.187	0.166	0.102

Table 5. Table 5 . Shapiro-Wilk’s normality test.

Shapiro-Wilk	P-value
$0.8765$	$0.002343$

Table 6. Table 6 . Estimation (Est) and Standard deviation (SD) of the parameters in four models: GAUS, SUG, GLG and SUGLG with the DIC and the LPML criterion to comparison these models.

		${\beta }_0$	$\alpha$	${\sigma }^2$	${\omega }^2$	${\theta }_w$	${\theta }_{\lambda }$	$\nu$	DIC	LPML
GAUS	Est	0.006	-	2.13	0.291	0.036	0.521	-	57.3	19.4
	SD	1.71	-	1.03	0.84	2.17	1.97	-
SUG	Est	0.043	2.76	1.98	0.284	0.251	0.56	-	40.7	31.1
	SD	1.65	1.19	1.11	0.761	2.061	2.01	-
GLG	Est	0.044	-	1.9	0.276	0.232	0.473	0.32	39.6	30
	SD	1.13	-	1.86	0.874	1.81	1.68	0.24
SUGLG	Est	0.011	2.22	1.4	0.197	0.461	0.465	2.61	11.4	43.7
	SD	0.33	0.95	0.41	0.29	0.925	0.71	0.13

Table 7. Table 7 . Bayesian spatial prediction of rainfall values in five stations.

	Observed value	Predicted value
Darab	0	0.0070121
Firuzabad	0	0.009761
Lamerd	0	0.0099013
Larestan	0	0.0095678
Qir-o-Karzin	0	0.0095411

Table 8. Table 8 . Rainfall amounts of Fars province measured in inch, geographical locations, and elevation in meter.

station	precipitation	longitude	latitude	elevation
Abadeh	1.09345684	52.40	31.11	2030
Arsanjan	0.48489401	53.16	29.56	1703
Bavanat	0.21459461	53.40	30.28	2231
Darab	0	54.17	28.47	1098
Eqlid	1.29816493	52.38	30.54	2300
Estahban	0.41391564	54.02	29.05	1609
Farashband	0.37737203	52.06	28.48	782
Fasa	0.47246634	53.41	28.56	1288
Firuzabad	0	52.33	28.53	1362
Gerash	0.29725655	54.15	27.69	403
Jahrom	0.21737661	53.32	28.29	1082
Kavar	0.27805494	52.65	29.16	651
Kazerun	0.31700190	51.39	29.36	860
Kherameh	0.26663226	53.29	29.63	875
Khonj	0.55668826	53.40	27.98	511
Khorrambid	0.77678033	53.09	30.35	2251
Lamerd	0	53.12	27.22	405
Larestan	0	54.17	27.42	792
Mamasani	0.83764997	51.32	30.04	972
Marvdasht	0.51173630	52.54	29.56	1605
Mohr	0.23188202	52.88	27.55	659
Neyriz	0.58328311	54.20	29.12	1632
Pasargad	0.74427868	53.21	30.19	1614
Qir-o-Karzin	0	53.03	28.28	746
Rostam	1.49020491	51.51	30.25	864
Sarvestan	0.42793923	53.21	29.27	719
Sepidan	1.61606073	52.00	30.14	2201
Shiraz	0.68990878	52.36	29.32	1484
Zarghan	0.22756008	52.43	29.47	1596
Zarrindasht	0.34310578	54.25	28.21	1029

Equations

$$\mathcal{D}=\left\{Y\left({\text{s}}_{k_j}\right)=y_{k_j};j\in J\right\}\cup \left\{Y\left({\text{s}}_{k_i}\right)\in A_i;i\in I\right\},$$

$$\mathrm{Y}\left(s\right)=\mu \left(s\right)+\frac{\mathrm{W}\left(s\right)}{\sqrt{\lambda \left(s\right)}}+\tau \rho \left(s\right),\tau \ge 0,$$

$$\mathbf{W}\sim {\mathrm{SUN}}_{n,n}\left(0,\left({\alpha }^2+{\sigma }^2\right)C_{{\theta }_w},\alpha C_{{\theta }_w},0,C_{{\theta }_w}\right).$$

$$\displaystyle{\bf Y}{\mathop{=}\limits^{d}}{\rm X}\mbox{\boldmath$\beta$\unboldmath}+\alpha{\Lambda^{-\frac{1}{2}}}{\bf U}+\sigma{\Lambda^{-\frac{1}{2}}}{\bf V}+\tau\mbox{\boldmath$\rho$\unboldmath}.$$

$${\left(\ln \left({\lambda }_1\right),\dots ,\ln \left({\lambda }_n\right)\right)}^{\prime }\sim {\mathrm{N}}_n\left(-\frac{\nu }{2}\mathrm{J},\nu C_{{\theta }_{\lambda }}\right),$$

$$\displaystyle p\left({\bf y}\left|\mbox{\boldmath$\lambda$\unboldmath},\eta\right.\right)\sim{{\rm SUN}_{n,n}}\left({\rm X}\mbox{\boldmath$\beta$\unboldmath},{\Sigma_{y}},\alpha{\Lambda^{-\frac{1}{2}}}{C_{\theta_{w}}},0,{C_{\theta_{w}}}\right).$$

$$L\left(\eta |\mathcal{D}\right)=\underset{{\mathbb{R}}^{{+}^n}}{\int }\underset{{\mathbb{R}}^{{}^n}}{\int }\underset{{\prod }_{i\in I}}{\int }$$

$$\displaystyle{{\bf Y}{{}_{J}}}{\left|\mbox{\boldmath$\lambda$\unboldmath},\eta\right.}\sim{{\rm SUN}_{m,m}}\left({X_{J}}\mbox{\boldmath$\beta$\unboldmath},{\Sigma_{y_{J}}},\alpha{{\Lambda_{J}}^{-\frac{1}{2}}}{C_{{JJ}_{\theta_{w}}}},0,{C_{{JJ}_{\theta_{w}}}}\right)$$

$$\displaystyle{\mu}_{y_{I}\left|{y_{J}},\mbox{\boldmath$\lambda$\unboldmath},\eta\right.}=X_{I}\mbox{\boldmath$\beta$\unboldmath}+{\rm A}\left({{y_{J}}}-X_{J}\mbox{\boldmath$\beta$\unboldmath}\right),$$

$$\displaystyle{\Sigma}_{y_{I}\left|{y_{J}},\mbox{\boldmath$\lambda$\unboldmath},\eta\right.}=\left({\sigma^{2}}+{\alpha^{2}}\right){\Lambda_{I}^{-\frac{1}{2}}}{C_{II_{\theta_{w}}}}{\Lambda_{I}^{-\frac{1}{2}}}+{\tau^{2}}{I_{n-m}}-{\rm A}\left({\Lambda_{J}^{-\frac{1}{2}}}{C_{JI_{\theta_{w}}}}{\Lambda_{I}^{-\frac{1}{2}}}\right),$$

$$\displaystyle{\Gamma}_{y_{I}\left|{y_{J}},\mbox{\boldmath$\lambda$\unboldmath},\eta\right.}={\Lambda_{I}^{-\frac{1}{2}}}{C_{IJ_{\theta_{w}}}}-A{\Lambda_{J}^{-\frac{1}{2}}}{C_{JJ_{\theta_{w}}}},$$

$$\displaystyle{\Delta}_{y_{I}\left|{y_{J}},\mbox{\boldmath$\lambda$\unboldmath},\eta\right.}={C_{JJ_{\theta_{w}}}}-{\alpha^{2}}{C_{JJ_{\theta_{w}}}}{\Lambda_{J}^{-\frac{1}{2}}}{\Sigma_{y_{J}}^{-1}}{\Lambda_{J}^{-\frac{1}{2}}}{C_{JJ_{\theta_{w}}}},$$

$$\mathrm{A}=\left({\sigma }^2+{\alpha }^2\right){\Lambda }_I^{-\frac{1}{2}}{C}_{I{J}_{{\theta }_w}}{\Lambda }_J^{-\frac{1}{2}}{\Sigma }_{y_J}^{-1},$$

$$C_{I{J}_{{\theta }_w}}={\left[C_{{\theta }_w}\left(\Vert s_{I_i}-s_{J_j}\Vert \right)\right]}_{\left(n-m\right)\times m}$$

$$\displaystyle\pi\left(\eta\right)=\pi\left(\mbox{\boldmath$\beta$\unboldmath}\right)\pi\left(\alpha\right)\pi\left({\sigma^{2}}\right)\pi\left({\omega^{2}}\right)\pi\left(\nu\right)\pi\left({\theta_{w}}\right)\pi\left({\theta_{\lambda}}\right).$$

$$\displaystyle p\left({y_{0}}\left|{\bf y}\right.\right)=\int_{\prod\nolimits_{i\in I}{{A_{i}}}}{\int_{R_{+}^{n}}{\int_{R_{+}^{p}}{p\left({{y_{0}}\left|{{\bf y},\mbox{\boldmath$\psi$\unboldmath},{\mbox{\boldmath$\psi$\unboldmath}_{0}}}\right.}\right)}}}p\left({{\mbox{\boldmath$\psi$\unboldmath}_{0}}\left|{{\bf y},\mbox{\boldmath$\psi$\unboldmath}}\right.}\right)p\left({\mbox{\boldmath$\psi$\unboldmath}\left|{\bf y}\right.}\right)d{\mbox{\boldmath$\psi$\unboldmath}_{0}}d\mbox{\boldmath$\psi$\unboldmath}d{y_{I}},$$

$$\left\{12,21,28,40\right\}\times \left\{10,20,30,40\right\},$$

$$\mathbf{Y}={\beta }_0+\alpha \mathbf{U}+\sigma \mathbf{V}+\tau \rho ,$$

$$\mathrm{Y}\left(s\right)={\beta }_0+\frac{\mathrm{W}\left(s\right)}{\sqrt{\lambda \left(s\right)}}+\tau \rho \left(s\right).$$

$$\displaystyle{\bf Y}{{}_{I}}\left|{\bf u},{v},\mbox{\boldmath$\lambda$\unboldmath},\eta,\mathscr{D}\right.\sim{{\rm N}_{n-m}}\left({\mu^{*}},{\Sigma^{*}}\right){\prod\limits_{i\in I}{{I_{\left\{{{y_{i}}\in{A_{i}}}\right\}}}}},$$

$${\mu }^{\ast }$$

$${\Sigma }^{\ast }$$

$$\displaystyle p\left({{\bf U}\left|{{\bf y},\mbox{\boldmath$\lambda$\unboldmath},\eta}\right.}\right)$$

$${\mu }_u^{\ast }$$

$$H_u$$

$$\displaystyle\pi\left({\mbox{\boldmath$\beta$\unboldmath}\left|{{y_{I}},{\bf u},\mbox{\boldmath$\lambda$\unboldmath},{\eta_{-\mbox{\boldmath$\beta$\unboldmath}}},\mathscr{D}}\right.}\right)$$

$$\displaystyle{\mu_{\mbox{\boldmath$\beta$\unboldmath}}^{*}}={{\bf H}_{\mbox{\boldmath$\beta$\unboldmath}}^{-1}}{X^{\prime}}{\Sigma_{y}^{-1}}\left({\bf y}-\alpha{\Lambda^{-\frac{1}{2}}}{\bf u}\right),$$

$$\displaystyle{{\bf H}_{\mbox{\boldmath$\beta$\unboldmath}}}={X^{\prime}}{\Sigma_{y}^{-1}}+\frac{1}{c_{0}}{{\rm I}_{k}}.$$

$$\displaystyle\pi\left({\alpha\left|{{y_{I}},{\bf u},\mbox{\boldmath$\lambda$\unboldmath},{\eta_{-\alpha}},\mathscr{D}}\right.}\right)$$

$${\mu }_{\alpha }^{\ast }=h_{\alpha }^{-1}{\mathbf{u}}^{\prime }{\Lambda }^{-\frac{1}{2}}{\Sigma }_y^{-1}Z,$$

$$h_{\alpha }={\mathbf{u}}^{\prime }{\Lambda }^{-\frac{1}{2}}{\Sigma }_y^{-1}{\Lambda }^{-\frac{1}{2}}\mathbf{u}+\frac{1}{c_1}.$$

$$\displaystyle\pi\left({{\sigma^{2}}\left|{{y_{I}},{\bf u},\mbox{\boldmath$\lambda$\unboldmath},{\eta_{-{\sigma^{2}}}},\mathscr{D}}\right.}\right)$$

$$\displaystyle\pi\left({{\omega^{2}}\left|{{y_{I}},{\bf u},\mbox{\boldmath$\lambda$\unboldmath},{\eta_{-{\omega^{2}}}},\mathscr{D}}\right.}\right)$$

$$\displaystyle\pi\left({\nu\left|{{y_{I}},{\bf u},\mbox{\boldmath$\lambda$\unboldmath},{\eta_{-\nu}},\mathscr{D}}\right.}\right)$$

Full text

MSC 2010: 62H11, 62F15.

11footnotetext:  Department of Statistics, Eghlid Higher Education Institute, Eghid, Iran.

Email: vahidtadayon24@gmail.com

1. Introduction

The customary approach to spatial data modeling is to accept that the random field of interest is Gaussian. Nevertheless, this assumption always is threatened by the existence of non-Gaussianity features (such as: heavier tails or skewness) in data sets. Moreover, wrong Gaussian assumptions affect the accuracy of spatial predictions and cause bias in the resulting parameter estimates as well. So, the acceptance of the Gaussianity could be overly restrictive to interpret the quantity of interest. Therefore, there is a demand for a more flexible class of sampling models to perform efficient inference and spatial prediction in the resulting models.

For modeling only skewed spatial data without heavy tails, the most commonly adopted strategy is to use the previous model after data transformation (Tadayon and Khaledi, 2015). Even so, an appropriate transformation if there exist, may not easy to obtain (see, e.g., De Oliveira et al., 1997; Tadayon and Rasekh, 2018; Tadayon and Torabi, 2018). Kim and Mallick (2004) developed the skew Gaussian random field based on the skew normal distribution which has been introduced by Azzalini and Capitanio (1999). Zhang and El-Shaarawi (2010) recently introduced a new class of stationary process with skew normal marginal distributions. Using generalized skew Gaussian spatial field is also one of the other strategies which has been suggested by Prates et al. (2012). More details on skew spatial can be found in Rimstad and Omre (2014) and Karimi and Mohammadzadeh (2011).

On the other hand, the observed data may contain some isolated or grouped outliers which have extreme values compared to their neighboring observation values. Since these observations could provide additional information and knowledge about the model, they play an important role in the spatial statistics. Measurement error is the mainspring of encountering outlying observations and belonging to a region with larger observational variance relative to the rest, is another reason. Since the non-Gaussian characteristics of the data may be derived from these observations, in order to identify and treat with these regions, Palacios and Steel (2006) proposed a flexible class of sampling models called Gaussian log-Gaussian (GLG). In this model which contains a scale mixed random variable, is assumed that the mixed variable is log-normally distributed. Therefore, observations with small mixed random variable belong to a region with larger observational variance. In addition, to identify individual outliers, Fonseca and Steel (2011) considered a similar mixing in the nugget effect component. However, this approach is suitable only for symmetric heavier tail distributions and will fail to handle the skewed data (Zareifard and Khaledi, 2013; Tadayon, 2018). An interesting extension, therefore, would be to develop models based on mixtures of skew normal, which would yield a highly flexible and computationally tractable parametric model that could accommodate both multimodality and extreme skewness. This extension which is compatible with Genton and Zhang (2012) remedies, is the unified skew Gaussian-log Gaussian (SUGLG) model introduced by Zareifard and Khaledi (2013). They applied the stochastic approximation expectation-maximization (SAEM) algorithm to maximize the likelihood function.

In addition, on account of some limitations of the data collection mechanism, the complete spatial data are not always available and so the quantity of interest may not be exactly observed in some sampling locations. In such setting, often the exact values can be recorded only if they fall within a specified range, which they have been called “censored data”. One of the conventional approaches to deal with these data is to impute the censored values with “artificial data”, such as, the detection limits. However, this leads biased estimates (see, e.g., Fridley and Dixon, 2007; Tadayon, 2018). Dubrule and Kostov (1986) worked on an interpolation method which motivated by the dual formulation of kriging. They also used geostatistical approaches requiring only the second-order specification of the random field which represent a type of censored data. Utilize likelihood-based approaches for the Gaussian random field is another previous work on the analysis of spatial models in the presence of censored data which well done by Militino and Ugarte (1999). De Oliveira (2005) introduced a Bayesian approach for inference and spatial prediction based on censored data while Fridley and Dixon (2007) incorporate the spatial correlation through an unobserved latent spatial process. Rathbun (2006) applied the Robbins-Monro (1951) stochastic approximation algorithm for estimating the parameters of a spatial regression model with left-censored observations. Eventually, Toscas (2010) proposed a modification on Bayesian approach of De Oliveira (2005) to correct the Bayesian bias in estimation and prediction of spatially correlated left-censored observations.

In all aforementioned works the inference was conducted based on the Gaussianity of the random field. However, in practice, we often observe that the exploratory data analysis shows some version of heavier tails than the normal distribution such as skewness or/and outlier region and consequently, it violates the normality assumption. The application of some transformation is the most commonly adopted strategy. However, the existence of an appropriate transformation is one of the most important problems in this setting and the difficulty in interpretation issues is another. Tadayon and Khaledi (2015) developed and implemented a Markov chain Monte Carlo (MCMC) sampling strategy for inference and fitted the skew Gaussian model based on censored data. In this article, we intend to fit the SUGLG model defined by Zareifard and Khaledi (2013) on a set of spatial data which could take both skewness and heavy tails simoltaneously and also contains some censored data.

In the following section we describe the form of spatial censored data which considered in this article. The third section introduces the mixture model and derives some of its properties. Section 4 considers Bayesian estimation of the model parameters and describes data augmentation method as well as Bayesian prediction. Using simulated data, the identifiability of the parameters, sensitivity analysis and the application of the suggested model is studied in Section 5. Section 6 illustrates the usage of the proposed methodology on a real spatial data set which contain values of precipitation. Conclusions are presented in Section 7.

2. Censored Spatial Data

Let ${\bf y}=\left(y_{1},\ldots,y_{n}\right)^{\prime}$ be a realization of ${\bf Y}=\left({Y}\left(s_{1}\right),\ldots,{Y}\left(s_{n}\right)\right)^{\prime}$, represents the data measured at the sampling locations $s_{1},s_{2},\ldots,s_{n}$ in $R$, where $R$ is the region under study. In what follows, we assume that the mechanism that produced the censoring is uninformative, meaning that the censoring process is independent of the spatial process $Y(\cdot)$, but may be deterministic or random (De Oliveira, 2005). To argue the credibility of this assumption, consider a case in which the censoring time is independent of the survival time (Militino and Ugarte, 1999). Moreover, notice that this assumption is similar to that of data being missing at random. Thus, the observed data consist of exact observations measured at some sampling locations and interval observations measured in the rest as the result of censoring (Tadayon and Khaledi, 2015). In this manner, to avoid some difficulties, we consider that the data of size n, consist of $m$ exact observations and $n-m$ interval observations. A reasonable methodology to define $m$ in practice is related to the nature of the data. Therefore, the data would be denoted by

[TABLE]

such that $\left\{y_{k_{j}}~{};j\in J\right\}$ are the observed values and for any $i\in I$, $A_{i}$ is the interval where $Y\left(\textbf{s}_{k_{i}}\right)$ is known to belong. The censored interval $A_{i}$ can depend on the sampling situation and it varies from a location to another location.

3. Statistical Model

Let ${\rm Y}\left(\cdot\right)=\left\{{{\rm Y}\left({{s}}\right)};s\in R\right\}$ be the random field of interest where $R\subseteq\mathbb{R}^{d^{*}}$ and ${d^{*}}\geq 1$. Our starting point is the model

[TABLE]

with the mean surface $\mu\left(s\right)$ and the scale parameter $\tau\in\mathbb{R}^{+}$. $\mu\left(s\right)$ is assumed to be a linear function of ${\bf f}^{\prime}\left(s\right)$, suchlike \mu\left(s\right)={\bf f}^{\prime}\left(s\right)\mbox{\boldmath\beta\unboldmath} where, ${\bf f}^{\prime}\left(s\right)$ is a vector of $k$ known functions of the spatial coordinate and \mbox{\boldmath\beta\unboldmath}\in{\mathbb{R}^{k}} is the regression coefficient vector. The term $\rho\left(s\right)$ denotes an uncorrelated Gaussian process with mean 0 and unitary variance, modeling the so-called nugget effect, which allows for measurement error and small-scale variation. This Further, the second-order stationary error process ${\rm W}\left(s\right)$, is a unified skew Gaussian random process corresponding to the model defined by Arellano-Valle and Azzalini (2006). On the other hand, for any vector from this random process like ${\bf W}=\left({\rm W}\left(s_{1}\right),\ldots,{\rm W}\left(s_{n}\right)\right)^{\prime}$, we consider that the distribution of vector $\bf W$ is ${\bf W}{\mathop{=}\limits^{d}}\alpha{\bf U}+\sigma{\bf V}$ in which, ${\bf V}$ and ${\bf U}$ have n-variate normal distribution ${{\rm N}_{n}}\left({\bf 0},{C_{\theta_{w}}}\right)$ and truncated normal distribution ${\rm TN}_{n}\left({\bf 0};{\bf 0},{C_{\theta_{w}}}\right)$ respectively. The $\left(i,j\right)$-th element of the $n\times n$ correlation matrix ${C_{\theta_{w}}}$ is ${C_{\theta_{w}}}\left(\|{s_{i}}-{s_{j}}\|\right)$. Thus, ${\bf W}$ has multivariate unified skew-normal (SUN) distribution111See Appendix C

[TABLE]

By setting \mbox{\boldmath\rho\unboldmath}=\left(\rho\left(s_{1}\right),\ldots,\rho\left(s_{n}\right)\right)^{\prime} and $\Lambda=diag\left(\lambda_{1},\ldots,\lambda_{1}\right)$, equation (1) can be written as

[TABLE]

In addition, like Palacios and Steel (2006) we assume that $\ln{\lambda\left(s\right)}$ is a Gaussian random field with finite-dimensional distributions:

[TABLE]

for ${\rm J}^{\prime}=\left(1,\ldots,1\right)_{n\times 1}$ and $\nu>0$, which implies a lognormal distribution for each $\lambda_{i}$. Clearly, by tending $\nu$ to 0, $\lambda_{i}$ tends to 1 and so equation (1) reduces to unified skew Gaussian model. For simplicity, we adopt a same correlation matrices for W and $\lambda$ but with different parameters $\theta_{w}$ and $\theta_{\lambda}$, respectively, and also we suppose that both W and $\lambda$ have an isotropic exponential correlation function, separately. However, the first assumption might be somewhat restrictive especially in the case of spatial data without temporal replications (as identifiability problem in these parameters, if exist, can affect the estimates of parameters), but we can hope that this problem is possibly solvable when there are many observations in space. Section 5 confirms this claim. Following Palacios and Steel (2006), we set ${\omega^{2}}=\frac{\tau^{2}}{\sigma^{2}}$ for simplicity and identifiability.

The model (1) is able to accommodate each of individual outliers and areas of the space with an inflated variance relative to the rest, because of observations with small values $\lambda_{i}$ fall out of the mean surface.

Let \eta^{\prime}=\left(\mbox{\boldmath\beta\unboldmath},\alpha,{\sigma^{2}},{\omega^{2}},\nu,{\theta_{w}},{\theta_{\lambda}}\right) and ${\bf y}=\left(y_{I},y_{J}\right)=\left(y_{1},\ldots,y_{n}\right)$ be a single realization of the considered random field with $y_{i}={\rm Y}\left(s_{i}\right)$, then the conditional joint distribution of $\bf y$ given \mbox{\boldmath\lambda\unboldmath}=\left(\lambda\left(s_{1}\right),\ldots,\lambda\left(s_{n}\right)\right)^{\prime} and $\eta$ is

[TABLE]

Furthermore, the likelihood function for $\eta$ given $\mathscr{D}$ is

[TABLE]

where ${\bf w}=\left(w\left(s_{1}\right),\ldots,w\left(s_{n}\right)\right)^{\prime}$, ${\Sigma_{y}}={\Sigma_{w}}+{\tau^{2}}{\rm{\bf I}}_{n}$ and ${\rm{\bf I}}_{n}$ is the identity matrix of dimension $n$. By setting ${C_{{JJ}_{\theta_{w}}}}=\left[{C_{\theta_{w}}}\left(\|{s_{J_{i}}}-{s_{J_{j}}}\|\right)\right]_{m\times m}$, we can easily see that

[TABLE]

and the conditional distributions of interval observations {\bf Y}{{}_{I}}\left|{y_{J}},\mbox{\boldmath\lambda\unboldmath},\eta\right. can be written as ${{\rm SUN}_{n-m,n-m}}{\prod\limits_{i\in I}{{I_{\left\{{{y_{i}}\in{A_{i}}}\right\}}}}}$ in which its parameters are

[TABLE]

where

[TABLE]

and the rest is [math]. It is necessary to note that $I_{\left\{{\rm A}\right\}}$ shows the indicator function of the set $\rm A$.

4. Bayesian Inference

To complete the Bayesian model specification, we need to select some prior distributions for all unknown parameters. One can then easily see that a convenient strategy of avoiding improper posterior distribution, is to utilize proper (but diffuse) priors. However, in this case, the insignificance of the prediction sensitivity to the hyperparameters should be established. For convenience but not necessary optimal, we assume elements of the parameter vector $\eta$ to be independent a priori, which means that

[TABLE]

The prior distributions adopted are as follows: \mbox{\boldmath\beta\unboldmath}\sim{{\rm N}_{k}}\left({\textbf{0}},{c_{0}}{{\rm{\textbf{I}}}_{k}}\right), $\alpha\sim{{\rm N}}\left({0},{c_{1}}\right)$, ${\sigma^{-2}}\sim{\rm Gamma}\left({c_{2}},{c_{3}}\right)$, ${\omega^{-2}}\sim{\rm GIG}\left(0,{c_{4}},{c_{5}}\right)$, $\nu\sim{\rm GIG}\left(0,{c_{6}},{c_{7}}\right)$, ${\theta_{w}}\sim{\rm Exp}\left(\frac{c_{8}}{{med}\left({d}\right)}\right)$ and ${\theta_{\lambda}}\sim{\rm Exp}\left(\frac{c_{9}}{{med}\left({d}\right)}\right)$, where the hyperparameters ${c_{1}},\ldots,{c_{9}}$ are chosen to reflect vague prior information and ${{med}\left({d}\right)}$ is the median of all distance between the data locations. Moreover, the closed form of generalized inverse-Gaussian (GIG) distribution and some justification for utilizing the above priors can be found in Palacios and Steel (2006).

Now, we try to derive the posterior distribution of the parameters by combining the likelihood function (2) and the prior distribution (3). Since a direct Bayesian analysis of the model is computationally impractical because of multiple integrals in posterior distribution, we adopt the data augmentation method and use the latent variables U and $\lambda$ to produce some samples $\{({y_{I}^{\left(i\right)}},{{u}^{\left(i\right)}},{\lambda^{\left(i\right)}},{\eta^{\left(i\right)}})\}_{i=1}^{l}$ from the extended model p\left({y_{I}},{\bf U},{\mbox{\boldmath\lambda\unboldmath}},\eta,{\left|\mathscr{D}\right.}\right) via MCMC methods. Details are given in Appendix A.

It must be noted that in many applications, the prediction of values at unsampled locations is a usual task of interest. The spatial prediction of the response in an unobserved location $y_{0}$ would be based on the Bayesian predictive distribution:

[TABLE]

where {\mbox{\boldmath\psi\unboldmath}}\left(s\right)=\ln\mbox{\boldmath\lambda\unboldmath}\left(s\right). Evidently, the full predictive distribution can not be evaluated in closed form but can be approximated using Monte Carlo samples.

5. Simulation Studies

The main aim of this section is to assess the performance of the proposed methodology on detecting two levels of censoring and also represent prediction values at some new unobserved locations. But since the SUGLG model may have a problem to indentify the correlation parameters (see Section 3), so, we first use simulation to evaluate the identifiability of these parameters. Then we present a banal method to revalue the ability of the this model to correctly identify outlying observations as well. Throughout, results are based on every 100th draw from an MCMC chain of length 150,000 with a burn-in of 100,000 which provides enough output for convergence. The spatial sampling points were considered on coordinates of the data file 97data.dat available from GSLIB software (Deutsch and Journel, 1998) in which, 97 locations are taken on a pseudo-regular grid over a bidimensional region 50 by 50 miles.

First of all, we appoint a lattice locations

[TABLE]

as a hold-out data set to evaluate the performance of the model in spatial prediction and then data are generated on these 113 locations based on the sampling model

[TABLE]

with a constant mean surface $\beta_{0}=0$, $\alpha=3,\sigma=1,\tau=0.1$ and ${\theta_{\lambda}}={\theta_{w}}=0.5$. After that, a cluster of outliers at locations 29, 37, 59, 78 and 84 has been created by imposing two units to the simulated values. Then, by leaving the simulated value in the lattice locations (4), the model will be fitted on the rest (of size 97). For model validation, we now implement two levels of left-censored design, one 17.5$\%$ and other 67$\%$ as a balanced and extreme cases of censoring respectively, such that in the first of them, the 17th ordered observation will be substituted for 1st up to 17th smallest values and in the second, 1st up to 65th smallest observations will be replaced with 65th ordered observation. Fig. 1 shows a schematic description of the region that displays the sampling locations as well as the region with inflated variance and censored locations. As an exploratory analysis, Fig. 2 which shows the histogram with the nonparametric density estimator of generated data, confirm the existence of an outlier region in these data as well as skewness.

Since the inference may be challenging in identifying the correlation parameters $\theta_{w}$ and $\theta_{\lambda}$, here, we focus on this problem to see what extent information about these parameters can be recovered from data. This study will be done just based on the balanced case of censoring because we do not expect a high ability to identify the parameters in the extreme case. To check the aforementioned problem, three data sets were generated with $\theta_{w}=0.25,1,1.5$ and $\theta_{\lambda}=0.1,0.7,1.1$ (in which, the other values have been fixed) and the estimate of the model parameters are obtained. However there are quite moderate sample size, the results which are reported in Table 1, indicates that there is no significant problem on the identifying of the correlation parameters and so the data almost allow for acceptable inference on these parameters.

At the moment, we concentrate on the performance of the Bayesian estimation, in which we must do one of these actions: choose some benchmark values for the various hyperparameters used in the priors, using empirical Bayes or employing hierarchical Bayes. To facilitate the computations, we consider the first method. But, we must investigatethe robustness of the posterior results under changes in the prior hyperparameters. It must be noted that we do not expect stability of the posteriors per se, but rather our interest is often in prediction. For this purpose, we employ the relative change criterion, which is the absolute value of the induced change in the marginal posterior mean of each parameter divided by the standard deviation computed under the benchmark prior. Table 2 which lists the alternative prior hyperparameters as well as the maximum relative change recorded for each parameter, shows no significant changes.

Therefore, the sensitivity of posterior can be partitioned into three groups: First, the prior on $\left({\beta_{0}},\alpha,{\sigma^{2}}\right)$ which does not play a vital role, $\left({\omega^{2}},{\theta_{w}},{\theta_{\lambda}}\right)$ where the data is reasonably informative on these parameters as the second and the last is $\nu$ by the largest influence of prior changes. Since there is relatively little direct information in the data on $\nu$, choice of prior is quite important. However, the adoption of hierarchical or empirical Bayes can be helpful.

Table 3 summarizes the results of fitting four models:

• •

Gaussian (GAUS): ${\rm Y}\left(s\right)={\beta_{0}}+\tau\rho\left(s\right)$,

• •

Unified Skew Gaussian (SUG): ${\rm Y}\left(s\right)={\beta_{0}}+{\rm W}\left(s\right)+\tau\rho\left(s\right)$,

• •

Gaussian-Log Gaussian (GLG): ${\rm Y}\left(s\right)={\beta_{0}}+{{\lambda\left(s\right)}}^{-\frac{1}{2}}+\tau\rho\left(s\right)$,

• •

Unified Skew Gaussian-Log Gaussian (SUGLG):

[TABLE]

According to this table, it is found that the SUGLG model has a better performance. Moreover, to compare four models in prediction, we compute the root of mean-square error of the predicted response values on the determined lattice locations. Table 4 suggets the desirability of using the SUGLG model.

Finally, by plotting the Bayesian estimate of $\lambda_{i}$’s for $i=1,\ldots,n$, it can be seen that how posterior values of $\lambda_{i}$’s might be used to identify outliers (because of the effect of small values of $\lambda_{i}$’s in the denominator, see equation (1)). Fig. 3 which shows mean of the posterior values $\lambda_{i}$’s for $i=1,\ldots,97$, confirms the existence of an outlier region in the clustered locations 29, 37, 59, 78 and 84.

6. Rainfall Data Analysis

The analyzed data set is comprised of rainfall amounts, measured in inch and accumulated on the first and coldest Wednesday in December 2012 at 30 stations located in Fars province of Iran. A schematic description of the region, the stations and rainfall amounts is shown in Fig. 4. In this figure it can be seen that there is a region with larger observation relative to the rest, along the north-eastern of the study region. In this data set (which are presented in Appendix B), the zero values are reported for five stations. Nevertheless, since the nonzero values are reported for the countrysides of these stations by residents and the stations operators, we consider their values as censored with censored intervals $\left[0,0.01\right)$, where the value 0.01 is the least value of precipitation in their countrysides. These locations are indicated by ”$\bigtriangledown$” as a left-censored observations in Fig. 4.

First of all, we performed the exploratory data analysis (EDA) such as the histogram, Q-Q plot, the Shapiro-Wilk’s test, the empirical semivariogram and the Haining (1991) method to detect outlier observations. Some of those results are presented in Fig. 5.

The histogram and Q-Q plot suggest that the data set has moderately right-skewed distribution which was confirmed by using the Shapiro-Wilk test where reported in Table 5. In addition, comparing the Haining bounds and the response values, indicate a region in the northings of studying space with larger observational variance relative to the rest, where we observed large amounts of precipitation. Finally, it is clear from the empirical semivariograms that there exists a strong spatial correlation as well as a nugget effect in the data set.

Before starting to compute estimations and spatial predictions, we mention that, since explanatory analysis of the data did not show any significant relation between response value and the spatial locations, the mean function is assumed to be constant, so we have just $\beta_{0}$ and thus $k=1$. This exploratory data analysis suggests that SUGLG model is a suitable option for doing the analysis. To perform Bayesian inference, benchmark values of the hyperparameter have chosen same as the given values in Table 2. The results is represented in Table 6. Although the DIC and the LMPL ($=\sum\nolimits_{i=1}^{n}{\log\left({CP{O_{i}}}\right)}$)222The conditional predictive ordinate (CPO): The CPO model comparison is a Bayesian cross-validation approach (e.g., Geisser, 1993) criterion present a versus results in comparison between two SUG and GLG models in some situations, but both of them confirm the better performance of the SUGLG model.

Morover, the prediction of rainfall amounts in the five mentioned stations, is presented in Table 7. One of the attractive feature of the Bayesian approach of course, is that it allows to make inference about the censored values, and in particular to quantify uncertainty about them. Fig. 6 displays the histograms of the posterior samples of the censored values, where vertical lines are placed at their respective censoring limits. It is clear that for some of these locations there is a lot of uncertainty about the true value, while for others the respective censoring limit is a poor estimate of the true value.

Finally, the contour map corresponding to the predictive mean surface under the SUGLG case, is shown in Fig. 7 which are computed from the predictive distribution over a regular grid of $20\times 20$. According to this figure, the predictions are highest in the area of inflated variance that contains the aforementioned observations. Fig. 8 presents the standard deviation surfaces for the SUGLG case as well.

7. Conclusions

In this article we have considered a unified skew Gaussian-log Gaussian model to analyze of non-Gaussian random fields based on censored data. Overall, the Bayesian framework, data augmentation method and MCMC algorithms have been extended for parameter estimation. However, using this algorithm entail more complication in the computations, recently, some studies have been performed to find replacing framework and it has led to some new methods such as Variational Bayes’ method (Ren et al. 2011). Although this method requires more complex theoretic calculations, it could increase the speed of calculations. So, assessing the performance of this method in Bayesian inference of Non-Gaussian random fields in the presence of censored values is an interesting area to investigate in further research.

Appendix A: The Full Conditional Distributions

In sequel, we try to present the full conditional distributions of all unknown quantities p\left({y_{I}},{\bf U},{\mbox{\boldmath\lambda\unboldmath}},\eta,{\left|\mathscr{D}\right.}\right) in the Gibbs sampler framework. Moreover, to facilitate the computations we put the notation ${\eta_{-\phi}}$ to show the vector $\eta$ without the element $\phi$ and set Z={\bf y}-X\mbox{\boldmath\beta\unboldmath}.

• $\bullet$

**Full conditional distribution of censored values

**By considering the partition of ${\bf u}=\left(u_{I},u_{J}\right)$ we have

[TABLE]

where,

[TABLE]

which is a $(n-m)$-variate truncated normal distribution and so it can easily generated by standard algorithm.

• $\bullet$

Full conditional distribution of latent variable ${\bf U}$

[TABLE]

which means that {\bf U}\left|{{\bf y},\mbox{\boldmath\lambda\unboldmath},\eta}\right.\sim{{\rm TN}_{n}}\left({\bf 0};{\mu_{u}^{*}},{H_{u}^{-1}}\right) by

[TABLE]

• $\bullet$

Full conditional distribution of regression coefficients

[TABLE]

So, \mbox{\boldmath\beta\unboldmath}\left|{\bf y},{\bf u},{\mbox{\boldmath\lambda\unboldmath}},{\eta_{-\mbox{\boldmath\beta\unboldmath}}}\right.\sim{{\rm N}_{k}}\left({\mu_{\mbox{\boldmath\beta\unboldmath}}^{*}},{{\bf H}_{\mbox{\boldmath\beta\unboldmath}}^{-1}}\right) with

[TABLE]

• $\bullet$

Full conditional distribution of parameter $\alpha$

[TABLE]

Thus, {\alpha}\left|{\bf y},{\bf u},{\mbox{\boldmath\lambda\unboldmath}},{\eta_{-\alpha}}\right.\sim{\rm N}\left({\mu_{\alpha}^{*}},{h_{\alpha}^{-1}}\right) where,

[TABLE]

• $\bullet$

**Full conditional distribution of parameter ${\sigma^{2}}$

**By assume that {\bf b}={\bf y}-X\mbox{\boldmath\beta\unboldmath}-\alpha{\Lambda^{-\frac{1}{2}}}{\bf u}=Z-\alpha{\Lambda^{-\frac{1}{2}}}{\bf u} and ${\bf B}={\Lambda^{-\frac{1}{2}}}{C_{\theta_{w}}}{\Lambda^{-\frac{1}{2}}}+{\omega^{2}}{{\rm I}_{n}}$ we have

[TABLE]

Since this full conditional distribution does not have any standard form, we use the Metropolis-Hastings method.

• $\bullet$

**Full conditional distribution of parameter ${\omega^{2}}$

**By recalling $\bf b$ and $\bf B$ from the above item, we also have

[TABLE]

and hence {\omega^{2}}\left|{\bf y},{\bf u},\mbox{\boldmath\lambda\unboldmath},{\eta_{-\omega^{2}}}\right.\sim{\rm GIG}\left(-\frac{n}{2},\sqrt{{c_{4}^{2}}+\frac{{{\bf b}^{\prime}}{{\bf B}^{-1}}{{\bf b}}}{\sigma^{2}}},{c_{5}}\right).

• $\bullet$

Full conditional distribution of parameter $\nu$

[TABLE]

In this case, the Metropolis-Hastings method can be used for sampling from this full conditional as well. Recall that this full conditional also can be represent in terms of {\bf a}=\ln\mbox{\boldmath\lambda\unboldmath}+\frac{\nu}{2}{\rm J}=\mbox{\boldmath\psi\unboldmath}+\frac{\nu}{2}{\rm J} in the penultimate line which is used in the next item. Therefore,

[TABLE]

• $\bullet$

**Full conditional distributions of the correlation parameters ${\theta_{\lambda}}$ and $\theta_{w}$

**In terms of $\bf a$, $\bf b$ and $\bf B$, we have

[TABLE]

and

[TABLE]

These full conditional distributions are of nonstandard form, so a Metropolis-Hastings step would be used.

Appendix B: Rainfall Data

Appendix C: Details On the Used Distributions

An $n$-dimensional random vector $Y$ is said to have a multivariate closed skew-normal distribution, denoted by CSN_{n,m}\left(\mbox{\boldmath\mu\unboldmath},\Sigma,D,{\bf v},\Theta\right), if its density is \phi_{n}\left({\bf y};\mbox{\boldmath\mu\unboldmath},\Sigma\right)\Phi_{m}\left(D\left({\bf y}-\mbox{\boldmath\mu\unboldmath}\right);{\bf v},\Theta\right)/\Phi_{m}\left({\bf 0};{\bf v},\Theta+D\Sigma^{\prime}D\right), where \mbox{\boldmath\mu\unboldmath}\in{\mathbb{R}}^{n}, ${\bf v}\in{\mathbb{R}}^{m}$ and $\Sigma\in{\mathbb{R}}^{n\times n}$ and $\Theta\in{\mathbb{R}}^{m\times m}$ are both covariance matrices, $D\in{\mathbb{R}}^{m\times n}$, \phi_{n}\left({\bf y};\mbox{\boldmath\mu\unboldmath},\Sigma\right) and \Phi_{n}\left({\bf y};\mbox{\boldmath\mu\unboldmath},\Sigma\right) are the probability density function (pdf) and cumulative distribution function (cdf), respectively, of the $n$-dimensional normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$.

Arellano-Valle and Azzalini (2006) presented a family of skew-normal distributions which unifies a plethora of $SN$ distributions including $CSN$ with no over-parametrization problem. According to their approach, an $n$-dimensional random vector $Y$ has a multivariate unified skew-normal distribution, denoted by SUN_{n,m}\left(\mbox{\boldmath\mu\unboldmath},\Sigma,\Gamma,{\bf v},\Delta\right), if its density is

[TABLE]

where $\Delta\in{\mathbb{R}}^{m\times m}$ is a correlation matrix, $\Gamma\in{\mathbb{R}}^{n\times m}$. The $SUN$ and $CSN$ classes are equivalent when $\Delta=\Theta+D\Sigma D^{\prime}$ and $\Gamma=\Sigma D$.

Let $\left\{{G\left({\bf s}\right),s\in R\subseteq{\mathbb{R}^{d}}}\right\}$, $d\geq 1$, be a spatial, ergodic, stationary, zero-mean Gaussian process with stationary covariance function $c\left({\bf h}\right)$ and denoted the covariance matrix of the random vector ${\bf G}={\left({G\left({{s_{1}}}\right),\ldots,G\left({{s_{n}}}\right)}\right)^{\prime}}$ by $\Sigma$. Also consider $V\sim{N_{n}}\left({\bf 0},{C_{\theta}}\right)$ and $U\sim{TN_{n}}\left({\bf 0};{\bf 0},{C_{\theta}}\right)$ are independent in which, {TN_{n}}\left({\bf c};\mbox{\boldmath\mu\unboldmath},\Sigma\right) denotes the N_{n}({\mbox{\boldmath\mu\unboldmath}},\Sigma) distribution truncated below at a point ${\bf c}$. Then if ${{\bf W}_{\lambda}}\mathop{=}\limits^{d}\alpha{\Lambda^{-\frac{1}{2}}}{\bf U}+\sigma{\Lambda^{-\frac{1}{2}}}{\bf V}$ and $\bf T$ be a normal vector of dimension $n$, we have

[TABLE]

Now, we define a $SUGLG$ process $\{W_{\lambda}(s)\}$ based on the following hierarchical representation:

1. (1)

$\left[{{W_{\lambda}}\left({\bf s}\right)\left|{\lambda\left({\bf s}\right)}\right.}\right]\mathop{=}\limits^{d}\left[{G\left({\bf s}\right)\left|{{\bf T}>{\bf 0},\lambda\left({\bf s}\right)}\right.}\right]$, 2. (2)

${\lambda\left({\bf s}\right)}$ is a log-Gaussian stochastic process.