On Statistical Properties of A Veracity Scoring Method for Spatial Data

TL;DR

This paper develops a veracity scoring method for spatial data without reference data, providing consistent parameter estimation and demonstrating advantages over traditional methods through simulations and real data application.

Contribution

It introduces a novel veracity scoring approach based on local summaries for spatial data lacking reference data, along with a consistent estimation method and theoretical analysis.

Findings

VS-based estimators are consistent under non-stationary noise.

VS estimators outperform OLS in asymptotic mean squared error.

Method successfully applied to real coal seam data.

Abstract

Measuring veracity or reliability of noisy data is of utmost importance, especially in the scenarios where the information are gathered through automated systems. In a recent paper, Chakraborty et. al. (2019) have introduced a veracity scoring technique for geostatistical data. The authors have used a high-quality `reference' data to measure the veracity of the varying-quality observations and incorporated the veracity scores in their analysis of mobile-sensor generated noisy weather data to generate efficient predictions of the ambient temperature process. In this paper, we consider the scenario when no reference data is available and hence, the veracity scores (referred as VS) are defined based on `local' summaries of the observations. We develop a VS-based estimation method for parameters of a spatial regression model. Under a non-stationary noise structure and fairly general assumptions on the underlying spatial process, we show that the VS-based estimators of the regression parameters are consistent. Moreover, we establish the advantage of the VS-based estimators as compared to the ordinary least squares (OLS) estimator by analyzing their asymptotic mean squared errors. We illustrate the merits of the VS-based technique through simulations and apply the methodology to a real data set on mass percentages of ash in coal seams in Pennsylvania.

Figures

Tables

Table 1. Table 1 : Comparison of orders of MSEs of OLS and VS under the example case: τ i 2 ≈ n c superscript subscript 𝜏 𝑖 2 superscript 𝑛 𝑐 \tau_{i}^{2}\approx n^{c} .

n	c	OLS-LB	VS-UB	RE-LB
100	0.1	0.016	0.010	1.585
	0.5	0.100	0.010	10.000
	0.8	0.398	0.010	39.811
500	0.1	0.004	0.002	1.862
	0.5	0.045	0.002	22.361
	0.8	0.289	0.002	144.270
1000	0.1	0.002	0.001	1.995
	0.5	0.032	0.001	31.623
	0.8	0.251	0.001	251.189

Table 2. Table 2 : Performance of regression parameter estimators under additive-multiplicative noise model.

		VS-based MSE		LS-based MSE		RE(Med-VS, OLS)=
${\sigma }_A$	$𝐧$	Med-VS	Mean-VS	OLS	EGLS	MSE(OLS)/MSE(Med-VS)
5	500	1.807	1.959	3.237	5.620	1.791
	1000	1.166	1.184	1.662	3.416	1.425
	5000	0.567	0.574	0.746	1.436	1.317
50	500	1.806	1.906	4.783	14.862	2.648
	1000	1.134	1.165	2.534	6.112	2.234
	5000	0.701	0.693	1.055	2.114	1.505
100	500	2.350	2.905	12.938	36.208	5.505
	1000	0.969	1.102	7.817	13.895	8.067
	5000	0.618	0.595	1.618	6.064	2.619
		VS-based MSE		LS-based MSE		R.E.(Med-VS, OLS)=
${\sigma }_M$	$𝐧$	Med-VS	Mean-VS	OLS	EGLS	MSE(OLS)/MSE(Med-VS)
0.447	500	1.889	1.945	2.771	4.481	1.467
	1000	1.005	1.033	1.559	2.942	1.552
	5000	0.572	0.578	0.699	1.352	1.220
0.707	500	2.219	2.253	4.098	13.248	1.847
	1000	1.063	1.117	2.685	6.031	2.525
	5000	0.647	0.641	0.919	2.527	1.420
0.953	500	2.221	2.422	6.182	20.996	2.784
	1000	1.201	1.206	3.578	9.048	2.978
	5000	0.669	0.671	1.476	4.746	2.206
		VS-based MSE		LS-based MSE		R.E.(Med-VS, OLS)=
$q_e$	$𝐧$	Med-VS	Mean-VS	OLS	EGLS	MSE(OLS)/MSE(Med-VS)
0.95	500	2.121	2.160	3.080	5.697	1.452
	1000	1.121	1.110	1.447	2.527	1.291
	5000	0.742	0.741	0.836	1.281	1.126
0.9	500	1.783	2.026	4.274	8.825	2.397
	1000	1.129	1.177	2.362	5.475	2.092
	5000	0.634	0.635	0.833	2.288	1.314
0.8	500	2.017	2.351	6.074	15.545	3.011
	1000	1.087	1.267	3.567	9.290	3.283
	5000	0.550	0.565	1.105	4.618	2.010

Table 3. Table 3 : Performance of covaraince parameter estimators under additive-multiplicative noise model.

		psill MSE		range MSE
${\sigma }_A$	$𝐧$	VS	WLS	VS	WLS
5	500	19.707	2058.760	2.471	3.099
	1000	1.255	1992.875	0.123	2.697
	5000	0.532	609.690	0.032	17.537
50	500	36.275	1334.724	4.549	0.580
	1000	1.800	1820.486	0.146	0.639
	5000	0.757	1163.299	0.046	11.401
100	500	45.653	682.913	9.477	0.502
	1000	1.675	913.053	0.186	0.506
	5000	0.933	1802.586	0.060	10.526
		psill MSE		range MSE
${\sigma }_M$	$𝐧$	VS	WLS	VS	WLS
0.447	500	7.712	1996.900	0.928	0.757
	1000	1.401	1863.389	0.164	5.779
	5000	0.580	577.577	0.047	6.312
0.707	500	14.324	2045.811	2.097	13.983
	1000	1.663	1566.315	0.170	4.272
	5000	0.820	792.043	0.038	8.928
0.953	500	28.337	1578.990	3.621	54.341
	1000	1.507	1537.173	0.201	16.551
	5000	1.134	1185.805	0.091	1.645
		psill MSE		range MSE
$q_e$	$𝐧$	VS	WLS	VS	WLS
0.95	500	11.786	1803.304	1.302	2.062
	1000	1.230	2385.300	0.146	2.884
	5000	0.665	531.284	0.047	12.393
0.9	500	5.216	2318.057	0.906	5.762
	1000	1.292	1535.030	0.098	5.707
	5000	1.063	817.137	0.026	2.803
0.8	500	36.674	1745.520	1.865	0.658
	1000	6.300	1308.637	0.116	13.993
	5000	4.584	1428.831	0.049	0.790

Table 4. Table 4 : Comparison of VS-based cov. estimation with GR

${\sigma }_{𝐀}$	$𝐧$	MSE(${\widehat{{\sigma }_{ϵ}}}^2$).VS	MSE(${\widehat{{\sigma }_{ϵ}}}^2$).GR	MSE(${\hat{\rho }}^2$).VS	MSE(${\hat{\rho }}^2$).GR
5	500	6.03	9.64e+22	1.38	3.95e+22
	1000	1.08	1.96e+23	0.11	1.77e+23
50	500	48.24	1.22e+23	7.14	1.38e+23
	1000	3.27	5.39e+22	0.68	1.75e+23
100	500	43.69	2.44e+22	9.44	9.19e+22
	1000	1.26	1.62e+21	0.18	2.40e+21

Table 5. Table 6 : Time comparison between VS-based and robust-REML. Machine configuration: DELL R7425 Dual Processor AMD Epyc 32 core 2.2 GHz machines with 512GB RAM each running 64Bit Ubuntu Linux Version 18.04.

Computation Time	VS-based	robust-REML
Min.	0.45	1.47
1st Qu.	0.54	1.99
Median	0.60	2.18
Mean	0.61	8.24
3rd Qu.	0.67	14.22
Max.	0.93	29.75

Equations

$$\begin{split}Y(\mbox{${\mathbf{s}}$})=\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$})^{\prime}\mbox{$\bm{\beta}$}+\epsilon(\mbox{${\mathbf{s}}$}),\end{split}$$

$$\begin{split}Z(\mbox{${\mathbf{s}}$}_{i})=\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i})^{\prime}\mbox{$\bm{\beta}$}+w(\mbox{${\mathbf{s}}$}_{i}),\end{split}$$

$$\begin{split}Z(\mbox{${\mathbf{s}}$}_{i})=\epsilon_{M_{i}}Y(\mbox{${\mathbf{s}}$}_{i})+\epsilon_{A_{i}}(\mbox{${\mathbf{s}}$}_{i}),\end{split}$$

$$\begin{split}V(\mbox{${\mathbf{s}}$}_{i})=\phi\left(\frac{\lvert Z(\mbox{${\mathbf{s}}$}_{i})-\mathcal{C}(\mathbf{Z}_{i})\rvert}{\alpha+D(\mathbf{Z}_{i})}\right),\end{split}$$

$$\begin{split}\hat{\mbox{$\bm{\beta}$}}_{\text{vs}}&=\underset{\mbox{$\bm{\beta}$}}{\text{argmin}}\;\sum_{i=1}^{n}V(\mbox{${\mathbf{s}}$}_{i})\left(Z(\mbox{${\mathbf{s}}$}_{i})-\mathbf{x}(\mbox{${\mathbf{s}}$}_{i})^{\prime}\mbox{$\bm{\beta}$}\right)^{2}\\ &=\left(\mbox{${\mathbf{X}}$}^{\prime}\mathbf{D}_{v}\mbox{${\mathbf{X}}$}\right)^{-1}\mbox{${\mathbf{X}}$}^{\prime}\mathbf{D}_{v}\mathbf{Z},\end{split}$$

$$\begin{split}\tilde{\epsilon}(\mbox{${\mathbf{s}}$}_{i})=V(\mbox{${\mathbf{s}}$}_{i})^{q}\hat{\epsilon}_{\text{vs}}(\mbox{${\mathbf{s}}$}_{i})+(1-V(\mbox{${\mathbf{s}}$}_{i})^{q})Q_{2}(\mathbf{Z}_{i}-\mathbf{X}_{i}\hat{\mbox{$\bm{\beta}$}}_{\text{vs}}),\end{split}$$

$$\text{Cov}\left(\epsilon(\mbox{${\mathbf{s}}$}_{i}),\epsilon(\mbox{${\mathbf{s}}$}_{j})\right)=\begin{cases}\sigma_{\epsilon}^{2}\;\rho_{\epsilon}\left(\mbox{${\mathbf{s}}$}_{i}-\mbox{${\mathbf{s}}$}_{j}\right)\;\;&\text{if}\;\;\mbox{${\mathbf{s}}$}_{i}\neq\mbox{${\mathbf{s}}$}_{j}\\ \sigma_{\epsilon}^{2}\;\;&\text{otherwise.}\end{cases}$$

$$\begin{array}{r}{ϵ}_{M_i}\sim \{\begin{array}{ll}\Delta (1)& \text{if}i\in G_n\\ F_M\left(\frac{\cdot -1}{{\sigma }_M}\right)& \text{o.w.}\end{array};{ϵ}_{A_i}\sim \{\begin{array}{ll}\Delta (0)& \text{if}i\in G_n\\ F_A\left(\frac{\cdot }{{\sigma }_A}\right)& \text{o.w.},\end{array}\end{array}$$

$$\begin{split}Z(\mbox{${\mathbf{s}}$}_{i})=\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i})^{\prime}\mbox{$\bm{\beta}$}+w(\mbox{${\mathbf{s}}$}_{i}),\end{split}$$

$${\sigma }_i^2$$

$$\begin{array}{r}{\lambda }_n^{-1}+n^{-1}{\lambda }_n^2\to 0.\end{array}$$

$$\begin{array}{r}{\alpha }_T(u,v)=\sup \left\{{\tilde{\alpha }}_T(A,B):d(A,B)\ge u;|A|\le v,|B|\le v\right\},\end{array}$$

$${\tilde{\alpha }}_T(A,B)=\sup \left\{|P(V_1\cap V_2)-P(V_1)P(V_2)|:V_1\in {\mathcal{F}}_T(A),V_2\in {\mathcal{F}}_T(B)\right\}.$$

$$V^{(m)}(\mbox{${\mathbf{s}}$}_{i})=\exp\left(-\frac{\lvert w(\mbox{${\mathbf{s}}$}_{i})\rvert}{\mathcal{I}_{n}(\mbox{${\mathbf{s}}$}_{i})}\right)+O_{p}(a_{n}),$$

$$\displaystyle C_{\epsilon}^{(l)}(q_{e})\leq\mathcal{I}_{n}(\mbox{${\mathbf{s}}$}_{i})=\xi_{n(i)}^{0.75}-\xi_{n(i)}^{0.25}\leq C_{\epsilon}^{(u)}(q_{e})$$

$$\begin{split}\hat{\mbox{$\bm{\beta}$}}_{\text{vs}}^{(m)}=\mbox{$\bm{\beta}$}+\left(\frac{1}{n}\mbox{${\mathbf{X}}$}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{${\mathbf{X}}$}\right)^{-1}\left(\frac{1}{n}\mbox{${\mathbf{X}}$}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}\right)+O_{p}\left(a_{n}\right),\end{split}$$

$$\begin{array}{r}E\left(\Vert {\mathbf{l}}_n^{\text{vs}}{\Vert }^2\right)\le C_3\cdot \left(n^{-1}{\left(C_1(q_e)\right)}^2+C_2{\lambda }_n^{-4}\right),\end{array}$$

$$\begin{array}{r}E\left(\Vert {\mathbf{l}}_n^{\text{ols}}{\Vert }^2\right)>n^{-1}\epsilon \left({\sigma }_{ϵ}^2\left(C^{(0)}-\epsilon \right){\left(C^{(1)}-\epsilon \right)}^{-1}+n^{-1}\sum _{i\in {\mathcal{M}}_n(\epsilon )\cap G_n^c}{\tau }_i^2\right)-C_1{\lambda }_n^{-4},\end{array}$$

$$\begin{array}{r}E\left(\Vert {\mathbf{l}}_n^{\text{ols}}{\Vert }^2\right)>n^{-1}\left(\epsilon {\sigma }_{ϵ}^2+n^{-1}\sum _{i\in G_n^c}{\tau }_i^2\right)-C_1{\lambda }_n^{-4},\end{array}$$

$$\begin{split}Y(\mbox{${\mathbf{s}}$}_{i})=\beta_{0}+\left(\beta_{x},\beta_{y}\right)^{\prime}\mbox{${\mathbf{s}}$}_{i}+\beta_{h}\;h(\mbox{${\mathbf{s}}$}_{i})+\epsilon(\mbox{${\mathbf{s}}$}_{i}),\quad\text{for}\;\;i\in\left\{1,\dots,n\right\},\end{split}$$

$$\begin{array}{rl}{ϵ}_{M_i}\underset{\text{indep.}}{\sim }2\times \text{Beta}({\alpha }_M,& {\alpha }_M),{ϵ}_{A_i}\underset{\text{indep.}}{\sim }N(0,{\sigma }_A^2).\end{array}$$

$$\begin{split}Y(\mbox{${\mathbf{s}}$})=\beta_{0}+\beta_{x}s_{x}+\epsilon(\mbox{${\mathbf{s}}$}),\end{split}$$

$${\alpha }_w\left(u,v\right)\le {\alpha }_{ϵ}\left(u,v\right),$$

$$\begin{array}{r}{M}_{n(i)}^p={\xi }_{n(i)}^p+\frac{p-{\hat{F}}_{n(i)}({\xi }_{n(i)}^p)}{{\bar{f}}_{n(i)}({\xi }_{n(i)}^p)}+R_{n(i)},\end{array}$$

$$\begin{split}Q_{2}\left(\mathbf{Z}_{i}\right)&=Q_{2}\left(\left(\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i_{1}})^{\prime}\mbox{$\bm{\beta}$}+w(\mbox{${\mathbf{s}}$}_{i_{1}}),\dots,\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i_{n(i)}})^{\prime}\mbox{$\bm{\beta}$}+w(\mbox{${\mathbf{s}}$}_{i_{n(i)}})\right)^{\prime}\right)\\ &=\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i})^{\prime}\mbox{$\bm{\beta}$}+Q_{2}\left(\underbrace{\left(w(\mbox{${\mathbf{s}}$}_{i_{1}}),\dots,w(\mbox{${\mathbf{s}}$}_{i_{n(i)}})\right)^{\prime}}_{=\mathbf{w}_{i}}\right)+O(\delta_{n}\lambda_{n}^{-1}),\end{split}$$

$$\begin{split}\left(\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i})-\mbox{${\mathbf{x}}$}(\mbox{${\mathbf{s}}$}_{i_{j}})\right)^{\prime}\mbox{$\bm{\beta}$}&=\sum_{k=1}^{p}\beta_{k}\left(x_{k}(\mbox{${\mathbf{s}}$}_{i})-x_{k}(\mbox{${\mathbf{s}}$}_{i_{j}})\right)\\ &=\sum_{k=1}^{p}\beta_{k}\left(\nabla x_{k}(\mbox{${\mathbf{s}}$}^{*}_{i_{j}})\right)^{\prime}\left(\mbox{${\mathbf{s}}$}_{i}-\mbox{${\mathbf{s}}$}_{i_{j}}\right)=O\left(\delta_{n}\lambda_{n}^{-1}\right),\end{split}$$

$$\displaystyle V^{(m)}(\mbox{${\mathbf{s}}$}_{i})$$

$$\displaystyle=\exp\left(-\frac{\lvert w(\mbox{${\mathbf{s}}$}_{i})+O_{p}(a_{n})\rvert}{\mathcal{I}_{n}(\mbox{${\mathbf{s}}$}_{i})+O_{p}(a_{n})}\right)$$

$$\displaystyle=\exp\left(-\frac{\lvert w(\mbox{${\mathbf{s}}$}_{i})\rvert}{\mathcal{I}_{n}(\mbox{${\mathbf{s}}$}_{i})}\right)+O_{p}(a_{n}),$$

$$p\ge \frac{1}{n(i)}\sum _{j\in G_{i,n}}{F}_{i_j}({\xi }_{n(i)}^p)\ge q_e{F}_{ϵ}({\xi }_{n(i)}^p/{\sigma }_{ϵ}),$$

$$\hat{\mbox{$\bm{\beta}$}}^{(m)}_{\text{vs}}-\mbox{$\bm{\beta}$}\;=\;\left(\frac{1}{n}\mbox{${\mathbf{X}}$}^{\prime}\mathbf{D}_{v}^{(m)}\mbox{${\mathbf{X}}$}\right)^{-1}\left(\frac{1}{n}\mbox{${\mathbf{X}}$}^{\prime}\mathbf{D}_{v}^{(m)}\mathbf{w}\right),$$

Full text

1 Introduction

With the advancement of sensor-related technologies in the field of data science and big-data analytics, ‘veracity’ is becoming one of the most important V’s, along with velocity, variety and volume. In a recent article, Lukoianovaand and Rubin (2014) mentioned how the data source, data collection technologies and data processing methodologies could induce bias, ambiguity and noise in real-world big-data applications. Practitioners often not only want to detect the noisy observations but also want to assign a reliability score to each of the observations, which efficiently indicates the ‘relative’ magnitude of the associated noise with it. For example, studies conducted by Evans (1997) show that dynamic display of veracity of the contents can help the map-users (e.g. in google maps) to take advantage of the reliability information in making decisions. Though there are some works on veracity analysis in media and communications (for example, see Conroy et al. 2015; Rendon et al. 2018 etc.) reliability assessment in the analysis of spatial data is not common in literature. Standard geostatistical analysis of the corrupted observations without taking the noisy nature of the data into account can result in erroneous inference and prediction. In a recent study, Chakraborty et al. (2019) introduced a reliability metric, namely veracity score (VS), to assess the credibility of the varying-quality crowdsourced observations in a geostatistical setting. In their analysis of mobile-sensor-generated noisy weather data, the authors have used a ‘high-quality’ reference data – coming from the ground weather stations – to define the veracity scores (VS) of the noisy observations. However, often in practice, such high-quality reference data is not available. In this paper, we consider the VS when there is no reference data available.

Although Chakraborty et al. (2019) have proposed veracity scoring based on ‘local’ summaries from the noisy data, statistical properties of the VS and corresponding VS-based estimators remain unexplored. In this paper, we investigate statistical properties of the VS-based methodology under fairly general conditions on the underlying spatial process. We assume an additive-multiplicative noise structure for the corrupted observations. In contrast to the more commonly used i.i.d. additive measurement errors (see e.g., Diggle et al. 2010), the underlying noise structure here is non-stationary due to the presence of the multiplicative component. Assuming a suitable spatial asymptotic framework (e.g. mixed-increasing domain, Hall and Patil 1994), we establish the consistency of the VS-based regression parameter estimator. The key result used here is a simplified asymptotic representation of the VS based on an extension of the Ghosh-Bahadur representation of sample quantiles of irregularly spaced spatial data. In addition, we show that the MSE of the VS-based estimator does not depend on the noise variances associated with the ‘bad’ observations asymptotically. This explains the empirical observation on robustness of the VS-based estimator against the noise inherent in crowdsourced spatial data as reported in Chakraborty et al. (2019). We also considered asymptotic behavior of the MSE of the OLS estimator. We show that the MSE of the OLS estimator is bounded below by a quantity that involves the sum of the associated noise variances. In situations where the combined noise level of the ‘bad’ observations is high, accuracy of the OLS estimator can be substantially inferior to the VS-based estimator.

Next, we evaluate finite sample properties of VS-based estimation though extensive simulations. As the proportion of noisy observations increases from $5\%$ to $20\%$, the relative efficiency (defined as ratio of MSEs) of the VS-based estimators of the regression parameters increases from $145\%$ to more than $300\%$ for a sample of size $500$ (see Table 2). Similar advantages of the VS-based estimators are noted at varying levels of the additive and multiplicative noise variances. We also carried out a comparative numerical analysis with the robust-REML method, a competing robust geostatistical technique proposed by Künsch et al. (2011) (also, see Papritz 2018a, Papritz 2018b). The VS-based technique outperforms the robust-REML method both in terms of scalability and statistical accuracy (see Table 4). Finally, we applied the VS-based analysis on a real data set containing measurements of mass percentages of ash in coal seams. The VS-based methodology identified all of the outliers previously detected in Cressie (1993) and Papritz (2018b)’s analysis of the coal ash data. More importantly, our analysis indicates the presence of spatial correlation (see Table 5) while the earlier analyses completely failed to capture the spatial dependence. Incorporation of the spatial dependence in the VS-based analysis led to comparable or much smaller point-wise prediction errors (indicated by the sizes of the circles in Figure 4(b)b).

The organization of the rest of the paper is as follows. In Section 2, we provide the background, introduce some notation, and briefly review the VS-based methodology. Section 3 defines the theoretical framework, states the assumptions and the results regarding asymptotic properties of the VS-based as well as the OLS estimators. In Section 4, we summarize the simulation design and simulation results for evaluating the VS-based methodology. Section 5 provides the details of the real data example. Finally, in Section 6, we summarize our findings and future directions to this work. Proofs of the main results of this paper are given in Appendix A. A supplementary material of this paper is available which contains the proofs of the auxiliary results and additional simulation results.

2 Veracity Score (VS) Methods

In this section, we briefly review the VS-based method in geostatistics. Before that, here we provide the background and a brief recap of some of the preliminary notation introduced in Chakraborty et al. (2019).

Let \left\{Y(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in\mathcal{D}\subset\mbox{{\mathbb{R}^{2}}}\right\} be the spatial process of our interest. We assume a spatial regression model, i.e. \left\{Y(\mbox{{\mathbf{s}}})\right\} has a decomposition of the form

[TABLE]

where \mbox{{\mathbf{x}}}(\cdot)=\left(x_{1}(\cdot),...,x_{p}(\cdot)\right)^{\prime} is the $p$-dimensional deterministic vector process of known covariates, $\bm{\beta}$ is the unknown regression parameter vector, and \left\{\epsilon(\mbox{{\mathbf{s}}})\right\} is a spatially correlated residual process. We assume that $\left\{\epsilon(\cdot)\right\}$ is an intrinsically stationary process with an admissible parametric variogram function 2\gamma(\mbox{{\mathbf{h}}};\mbox{\bm{\theta}})=\text{Var}\left\{\epsilon(\mbox{{\mathbf{s}}})-\epsilon(\mbox{{\mathbf{s}}}+\mbox{{\mathbf{h}}})\right\} where $\bm{\theta}$ is the covariance parameter of interest. Under second-order stationarity, it can be shown that \gamma(\mbox{{\mathbf{h}}})=C(\bm{0})-C(\mbox{{\mathbf{h}}}) where C(\mbox{{\mathbf{h}}})=\text{Cov}\left(\epsilon(\mbox{{\mathbf{s}}}),\epsilon(\mbox{{\mathbf{s}}}+\mbox{{\mathbf{h}}})\right) is the covariance function of the $\epsilon$-process. For more details about variograms and covariograms see Chapter 2 in Cressie (1993).

Now, instead of observing \left\{Y(\mbox{{\mathbf{s}}}_{1}),\dots,Y(\mbox{{\mathbf{s}}}_{n})\right\} at locations (can be irregularly space) \mathcal{S}_{n}=\left\{\mbox{{\mathbf{s}}}_{1},\dots,\mbox{{\mathbf{s}}}_{n}\right\}, we observe corrupted observations \left\{Z(\mbox{{\mathbf{s}}}_{1}),\dots,Z(\mbox{{\mathbf{s}}}_{n})\right\}. We write the corrupted observations as,

[TABLE]

where \left\{w(\mbox{{\mathbf{s}}})\right\} is a mean-zero spatially correlated process, possibly nonstationary, which contains noise in addition to the small-scale local variations induced by \epsilon(\mbox{{\mathbf{s}}}). For example, as mentioned by Chakraborty et al. (2019), let the corruption is coming through an additive-multiplicative noise model, i.e.

[TABLE]

where $\epsilon_{M_{i}}$ and $\epsilon_{A_{i}}$ are the randomm variables corresponding to the multiplicative and additive noise parts. Then the $w$-process can be written as w(\mbox{{\mathbf{s}}}_{i})=\epsilon_{M_{i}}(\mu(\mbox{{\mathbf{s}}}_{i})-1)+\epsilon_{M_{i}}\epsilon(\mbox{{\mathbf{s}}}_{i})+\epsilon_{A_{i}}. In case there is no multiplicative component, w(\mbox{{\mathbf{s}}}_{i})=\epsilon(\mbox{{\mathbf{s}}}_{i})+\epsilon_{A_{i}}. Observe that existance of the multiplicative component make the noisy residuals \left\{w(\mbox{{\mathbf{s}}})\right\} dependent on the location ${\mathbf{s}}$ and hence non-stationary in nature.

2.1 Formulation of VS

Denote a square $\delta$-neighborhood ($\delta>0$) around location ${\mathbf{s}}$ as \mathcal{B}_{\delta}(\mbox{{\mathbf{s}}}), i.e. \mathcal{B}_{\delta}(\mbox{{\mathbf{s}}})=(\mbox{{\mathbf{s}}}-\delta,\mbox{{\mathbf{s}}}+\delta], where the subtractions are component-wise. Then the VS of the observation Z(\mbox{{\mathbf{s}}}_{i}) is defined as (Chakraborty et al., 2019),

[TABLE]

where \phi:\mbox{{\mathbb{R}}}^{+}\cup\left\{0\right\}\to\mbox{{\mathbb{R}}}^{+}\cup\left\{0\right\} is some non-increasing bounded above function, referred as veracity function; $\alpha\geq 0$ is a regularity parameter, called the baseline deviation; \mathbf{Z}_{i}=\left(Z(\mbox{{\mathbf{s}}}_{i_{1}}),\dots,Z(\mbox{{\mathbf{s}}}_{i_{n(i)}})\right)^{\prime} where \{\mbox{{\mathbf{s}}}_{i_{1}},\dots,\mbox{{\mathbf{s}}}_{i_{n(i)}}\} is the set of observation locations in the small $\delta$-neighborhood \mathcal{B}_{\delta}(\mbox{{\mathbf{s}}}_{i}); and finally, $\mathcal{C}(\mathbf{x})$ is a measure of central tendency and $D(\mathbf{x})$ is a measure of dispersion of the values in the vector $\mathbf{x}$. Clearly, the definition of VS given in Equation 2.4 measures the amount of deviation of the observation Z(\mbox{{\mathbf{s}}}_{i}) from the ‘local’ summary $\mathcal{C}(\mathbf{Z}_{i})$ relative to the ‘local’ variation $D(\mathbf{Z}_{i})$. The veracity function $\phi(\cdot)$ being a non-increasing function, larger the value of the scaled deviation \frac{\lvert Z(\mbox{{\mathbf{s}}}_{i})-\mathcal{C}(\mathbf{Z}_{i})\rvert}{\alpha+D(\mathbf{Z}_{i})} lower is the VS and for lower values of the scaled deviation, VS will be higher.

For the measure of center and dispersion we have considered the following two choices that are widely used by practitioners:

$\mathcal{C}(\mathbf{Z}_{i})=Q_{2}(\mathbf{Z}_{i})$ and $D(\mathbf{Z}_{i})=\text{IQR}(\mathbf{Z}_{i})=Q_{3}(\mathbf{Z}_{i})-Q_{1}(\mathbf{Z}_{i})$, where Q_{j}(\mbox{{\mathbf{x}}}) denotes the $j$-th quartile of the observations in ${\mathbf{x}}$, $j\in\left\{1,2,3,4\right\}$; 2. 2.

\mathcal{C}(\mathbf{Z}_{i})=\bar{Z}_{i_{.}}=\frac{1}{n(i)}\sum_{j=1}^{n(i)}Z(\mbox{{\mathbf{s}}}_{i_{j}}), and D(\mathbf{Z}_{i})=\text{sd}(\mathbf{Z}_{i})=\sqrt{\frac{1}{n(i)-1}\sum_{j}\left(Z(\mbox{{\mathbf{s}}}_{i_{j}})-\bar{Z}_{i_{.}}\right)^{2}}.

Clearly, the first set of choices are based on sample quantiles and are expected to be more robust and efficient as compared to the sample mean and standard deviation in the analysis of noisy data (for details, see Sen 1968). In this paper, we refer the VS with choice (a) as the measure of location and scale as ‘Medain-VS’ and denote it by V^{(m)}(\mbox{{\mathbf{s}}}_{i}). Similarly, for choice (b) we denote the VS as V^{(a)}(\mbox{{\mathbf{s}}}_{i}) ($(a)$ for average) and refer it as ‘Mean-VS’. For the veracity function we have adopted the choice used by Chakraborty et al. (2019), i.e. $\phi(x)=\exp(-x)$.

2.2 VS-based estimation in spatial regression

Under the setting discussed above, as the covariance parameter is unknown, the standard practice is to estimate the regression parameter using ordinary least squares (OLS): \hat{\mbox{\bm{\beta}}}_{\text{ols}}=\left(\mathbf{X}^{\prime}\mathbf{X}\right)^{-1}\mathbf{X}^{\prime}\mathbf{Z}, where \mathbf{X}=\left(\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}_{1}),\dots,\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}_{n})\right)^{\prime}. Next, the residuals from the OLS fit are used in weighted least squares-based variogram estimation (for details, see page 90, Chapter 2, Cressie 1993) to obtain the covariance parameter estimator \hat{\mbox{\bm{\theta}}}_{\text{wls}}. Once the covariance parameter is estimated, one can try to improve the mean parameter estimates by using estimated generalized least squares (EGLS) estimator, given by \hat{\mbox{\bm{\beta}}}_{\text{egls}}=\left(X^{\prime}\Sigma^{-1}(\hat{\mbox{\bm{\theta}}}_{\text{wls}})X\right)^{-1}X^{\prime}\Sigma^{-1}(\hat{\mbox{\bm{\theta}}}_{\text{wls}})\mathbf{Z}, where \Sigma(\hat{\mbox{\bm{\theta}}}_{\text{wls}}) is the estimated covariance matrix. But, this last stage of updating the OLS fit introduces additional variability due to the estimation of the covariance parameters and hence not necessarily preferable than the OLS estimator. Finally, the estimated mean and covariance structure is used to predict the process at a new location \mbox{{\mathbf{s}}}_{0} using best linear unbiased predictor or kriging predictor (for details, see Chapter 3, Cressie 1993).

The standard approach of estimation and kriging is not robust in nature and hence, may produce erroneous inference and prediction in case there are outliers in the data. To estimate the mean parameter $\bm{\beta}$ robustly, instead of using ordinary squared error loss, in VS-based estimation, weighted squared error loss is used with VS as the corresponding weights (Chakraborty et al., 2019) as shown below:

[TABLE]

where \mathbf{D}_{v}=\text{diag}(V(\mbox{{\mathbf{s}}}_{1}),\dots ,V(\mbox{{\mathbf{s}}}_{n})). Clearly, as the VS is expected to be smaller for ‘bad’ observations, the VS-based estimator given in Equation 2.5 is expected to be affected less by the high noise associated with those. If the VS can capture the noisy observations perfectly, i.e. V(\mbox{{\mathbf{s}}}_{i})\approx 0 if Z(\mbox{{\mathbf{s}}}_{i}) is contaminated and V(\mbox{{\mathbf{s}}}_{i})=1 otherwise, then \hat{\mbox{\bm{\beta}}}_{\text{vs}} is approximately equal to the OLS estimator computed only from the ‘noiseless’ observations in the data.

Depending on the versions of VS used to compute the VS-based estimator, we consider two versions of the VS-based estimator in this article: ‘Mean-VS’, denoted by \hat{\mbox{\bm{\beta}}}^{(a)}_{\text{vs}}, where the sample mean and standard deviations are used to assess veracity and ‘Median-VS’, denoted by \hat{\mbox{\bm{\beta}}}^{(m)}_{\text{vs}}, with sample median and IQR in the definition of VS. Though the focus of this work is mainly on the properties of the ‘Median-VS’ regression estimator, in the following sections we put some remarks about the properties and asymptotic behaviors of the ‘Mean-VS’ version as well.

The de-trended observations, i.e. \hat{\epsilon}_{\text{vs}}(\mbox{{\mathbf{s}}}_{i})=Z(\mbox{{\mathbf{s}}}_{i})-\mathbf{x}(\mbox{{\mathbf{s}}}_{i})^{\prime}\hat{\mbox{\bm{\beta}}}_{\text{vs}} for $i\in\left\{1,2,...,n\right\}$ – unlike the case of dealing with high-quality observations – are not free of noise and hence, it is not reasonable to use these directly in the covariance parameter estimation. Chakraborty et al. (2019) proposed a VS-based smoothing of the residuals:

[TABLE]

where \mathbf{X}_{i}=\left(\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}_{i_{1}}),...,\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}_{i_{n(i)}})\right)^{\prime} is the $n(i)\times p$ matrix of the covariates corresponding to the observations in \mathcal{B}_{\bm{\delta}}(\mbox{{\mathbf{s}}}_{i}); and $q\geq 0$ is another regularity parameter which determines the smoothness of the residuals defined in Equation 2.6. Next, these VS-based smoothed residual vector \tilde{\mbox{\bm{\epsilon}}}=\left(\tilde{\epsilon}(\mbox{{\mathbf{s}}}_{1}),\dots,\tilde{\epsilon}(\mbox{{\mathbf{s}}}_{n})\right)^{\prime} is used in least squares based variogram model fitting (for example, see Cressie and Douglas 1980) to get the VS-based covariance parameter estimator \hat{\mbox{\bm{\theta}}}_{\text{vs}}. Finally, the prediction of the $\epsilon$-process at a new location \mbox{{\mathbf{s}}}_{0}, denoted by \tilde{\epsilon}(\mbox{{\mathbf{s}}}_{0}), is obtained using ordinary kriging with the estimated covariance parameters and the VS-based smoothed residuals as discussed in Chakraborty et al. (2019).

In the next section, we have focused on the asymptotic behavior of the VS and VS-based regression parameter estimator. Though investigation of asymptotic properties VS-based covariance estimation is beyond the scope of this paper, we have evaluated the same using simulations in Section 4.

3 Asymptotic properties of the VS-based regression estimator

Before going to the main results of this section we need to introduce some notation, specify the theoretical framework and state the regularity conditions. From this section onward, by $C(\cdot),C,C_{1},C_{2},\dots$ we denote constants with respect to the sample size $n$. For $d\geq 2$, we denote the volume of a set A\subset\mbox{{\mathbb{R}^{d}}} as $\lvert A\rvert$, i.e., the Lebesgue measure of $A$ if it has nonzero volume and the cardinality of $A$ if $A$ is finite. For two sequences of positive reals $\left\{a_{n}\right\}_{n=1}^{\infty}$ and $\left\{b_{n}\right\}_{n=1}^{\infty}$ we say one is of the order of another (denoted as $a_{n}\sim b_{n}$) if $a_{n}/b_{n}\to C>0$.

3.1 Model specification

As discussed before, the process of our interest is denoted as \left\{Y(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in\mathcal{D}\subset\mbox{{\mathbb{R}^{2}}}\right\} and instead of observing realizations from the ‘true’ process, we observe corrupted observations \left\{Z(\mbox{{\mathbf{s}}}_{1}),\dots,Z(\mbox{{\mathbf{s}}}_{n})\right\} where, the observations locations \left\{\mbox{{\mathbf{s}}}_{1},\dots,\mbox{{\mathbf{s}}}_{n}\right\} are possibly irregularly spaced. The results of this paper are specific to the spatial regression model defined in Equation 2.1. For the zero mean second-order stationary residual process \left\{\epsilon(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in\mathcal{D}\right\}, let us denote the marginal distribution function of \epsilon(\mbox{{\mathbf{s}}}) as $F_{\epsilon}(\frac{x}{\sigma_{\epsilon}})$ and, the spatial dependence structure can be formulated as,

[TABLE]

Here, $F_{\epsilon}(\cdot)$ is a distribution function of a zero mean unit-variance random variable and $\rho_{\epsilon}(\cdot)$ is a non-negative definite function on ${\mathbb{R}^{2}}$ with $\rho_{\epsilon}(\mathbf{0})=1$. To preserve the identifiability, we assume that the supports of $\left\{\epsilon_{M_{i}}:i\in\left\{1,\dots,n\right\}\right\}$ are contained in $[0,\infty)$.

We further assume that among \left\{Z(\mbox{{\mathbf{s}}}_{1}),\dots,Z(\mbox{{\mathbf{s}}}_{n})\right\} only a ‘small’ portion is ‘bad-quality data’ and the corrupted observations are coming through the additive-multiplicative noise model given in Equation 2.3. We formulate this assumption as follows.

[TABLE]

where, $\Delta(x)$ denotes a degenerate distribution with point mass at $-\infty<x<\infty$; $G_{n}\subset\left\{1,\dots,n\right\}$ is a subset of indices; and $\sigma_{M},\sigma_{A}$ are positive constants. With this model, if $i\in G_{n}$, we have no noise associated with the observation, i.e., Z(\mbox{{\mathbf{s}}}_{i})=Y(\mbox{{\mathbf{s}}}_{i}). If $i\notin G_{n}$, then Z(\mbox{{\mathbf{s}}}_{i})=\epsilon_{M_{i}}Y(\mbox{{\mathbf{s}}}_{i})+\epsilon_{A_{i}}, where $\epsilon_{M_{i}}$ and $\epsilon_{A_{i}}$ have positive variances. Also, we have taken $\left\{\epsilon_{M_{i}}\right\}_{i=1}^{n}$, $\left\{\epsilon_{A_{i}}\right\}_{i=1}^{n}$ and \left\{\epsilon(\mbox{{\mathbf{s}}}_{i})\right\}_{i=1}^{n} are independent of each other. We further assume that the proportion of “good” observations is a constant (w.r.t. $n$) denoted by $q_{e}$, i.e., $\lvert G_{n}\rvert/n\approx q_{e}$, and $1-q_{e}$ is the proportion of noisy observations in the data. Observe that, the variance corresponding to $\epsilon_{M_{i}}$, $\sigma_{M_{i}}=\sigma_{M}$ if $i\notin G_{n}$ and $\sigma_{M_{i}}=0$ if $i\in G_{n}$. Similar observation can be made for the additive noise variances (denoted by $\sigma_{A_{i}}$) as well.

Under the specified model we can rewrite our observations as,

[TABLE]

where E(w(\mbox{{\mathbf{s}}}_{i}))=0, \mathrm{Var}(w(\mbox{{\mathbf{s}}}_{i}))=\sigma_{i}^{2} and,

[TABLE]

where, \tau_{i}^{2}\equiv\tau(\mbox{{\mathbf{s}}}_{i})^{2}=(\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}_{i})^{\prime}\mbox{\bm{\beta}})^{2}\sigma_{M}^{2}+\sigma_{\epsilon}^{2}\sigma_{M}^{2}+\sigma_{A}^{2} is the additional noise variance associated with the corrupted observation at location \mbox{{\mathbf{s}}}_{i}. Clearly, if $i\in G_{n}$, $\tau_{i}^{2}=0$. We denote the marginal distribution function of w(\mbox{{\mathbf{s}}}_{i}) as $F_{i}(x)$ for $i\in\left\{1,\dots,n\right\}$, for any $n\geq 2$. Clearly, if $i\in G_{n}$, w(\mbox{{\mathbf{s}}}_{i})=\epsilon(\mbox{{\mathbf{s}}}_{i}) and $F_{i}(x)=F_{\epsilon}(x/\sigma_{\epsilon})$ and if $i\in G_{n}^{c}$, we denote the marginal distribution function of w(\mbox{{\mathbf{s}}}_{i}) as $F_{2}(\frac{x}{\sigma_{i}})$, where $F_{2}(\cdot)$ is again a distribution function of a centered and scaled random variable. Without loss of any generality, we take the baseline deviation parameter $\alpha$ in the definition of VS to be equal to [math]. The results of this paper can be straightforwardly extended for any other finite constant (w.r.t. $n$) $\alpha>0$.

3.2 Spatial framework and notations

We take the sampling region to be $\mathcal{D}\equiv\mathcal{D}_{n}=\lambda_{n}[0,1]^{2}$, i.e. a 2-dimensional square region with area $\lambda_{n}^{2}$ where $\left\{\lambda_{n}\right\}$ is a sequence of positive reals such that $\lambda_{n}\to\infty$ as $n\to\infty$ and

[TABLE]

Under condition (3.3), the spatial asymptotic framework we consider is similar to the mixed increasing-domain asymptotic framework used in Hall and Patil (1994), Lahiri et al. (1999), etc.The first component of (3.3) states the domain has to increase with the sample size, and the second component allows the possibility of infilling sampling points as well. The results stated here are not particular to the shape and position of the rectangular region $\mathcal{D}_{n}$. Our results can be extended to more general shapes using the asymptotic framework where the region is obtained by inflating a prototype region contained in $(-1/2,1/2]^{2}$ with center at the origin by a scaling constant as discussed in Lahiri et al. (2002).

Now, we introduce some notation required to state the regularity conditions, the results as well as the proofs of this paper. For any A\subset\mbox{{\mathbb{R}^{2}}} and a random field \left\{T(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in\mbox{{\mathbb{R}^{2}}}\right\} let us denote \mathcal{F}_{T}(A)=\sigma\langle T(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in A\rangle, the $\sigma$-field generated by the random variables \left\{T(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in A\right\}. For any $d\geq 2$ and for any \mathbf{x}\in\mbox{{\mathbb{R}^{d}}} the $L_{1}$ and $L_{2}$ norms are denoted by $\lvert\mathbf{x}\rvert$ and $\lVert\mathbf{x}\rVert$ respectively. For any two sets $A$ and $B$ in ${\mathbb{R}^{2}}$ we denote the the distance between them as d(A,B)=\inf\left\{\lvert\mbox{{\mathbf{s}}}-\mbox{{\mathbf{s}}}^{\prime}\rvert:\mbox{{\mathbf{s}}}\in A,\mbox{{\mathbf{s}}}^{\prime}\in B\right\}. For any random field \left\{T(\mbox{{\mathbf{s}}}):\mbox{{\mathbf{s}}}\in\mbox{{\mathbb{R}^{2}}}\right\} we define the strong-mixing coefficient as

[TABLE]

where the supremum is taken over the set of all 2-dimensional rectangles $A$, $B$ in ${\mathbb{R}^{2}}$, $u>0$, $v>0$; and

[TABLE]

For a distribution function $F(\cdot)$, we define the inverse of it as $F^{-1}(p)=\text{inf}\left\{x:F(x)\geq p\right\}$, for $p\in[0,1]$.

To definition of VS in Section 2 uses a square neighborhood \mathcal{B}_{\delta}(\mbox{{\mathbf{s}}}_{i}) around the spatial point \mbox{{\mathbf{s}}}_{i} with one side of $2\delta$ units. From this section we change the notation as: $\delta\equiv\delta_{n}$ and $\mathcal{B}_{\delta}(\cdot)\equiv\mathcal{B}_{\delta_{n}}(\cdot)$ to ensure that the size of the ‘local’ neighborhood also varies with the sample size $n$. Let $M_{n(i)}^{p}$ be the smallest $p$-th sample quantile of \left\{w(\mbox{{\mathbf{s}}}_{i_{1}}),\dots,w(\mbox{{\mathbf{s}}}_{i_{n(i)}})\right\}, the residuals associated with the observations in the square-neighborhood \mathcal{B}_{\delta_{n}}(\mbox{{\mathbf{s}}}_{i}). We denote \hat{F}_{n(i)}(x)=\left\{\text{No. of w(s_{i_{j}})\leq x}\right\}/n(i), the empirical distribution function of the realizations from the $w$-process in the neighborhood around \mbox{{\mathbf{s}}}_{i}. Clearly, $M_{n(i)}^{p}=\hat{F}_{n(i)}^{-1}(p)$. By $\bar{F}_{n(i)}(x)$ we denote the distribution function defined as $\bar{F}_{n(i)}(x)=\left(\sum_{j=1}^{n(i)}F_{i_{j}}(x)\right)/n(i)$, where $F_{i_{j}}$ is the distribution function of w(\mbox{{\mathbf{s}}}_{i_{j}}). The smallest $p$-th quantile of the distribution function $\bar{F}_{n(i)}$ is denoted by $\xi_{n(i)}^{p}$, i.e. $\xi_{n(i)}^{p}=\bar{F}_{n(i)}^{-1}(p)$.

3.3 Consistency of the VS-based regression parameter estimator

Here we focus on consistency of the \hat{\mbox{\bm{\beta}}}_{\text{vs}}^{(m)}, the Median-VS version of VS-based estimator. The conditions we need for that are following.

• (C.1)

Number of observations in any unit square in $\mathcal{D}_{n}$ is in between $[C_{1}\frac{n}{\lambda_{n}^{2}},C_{2}\frac{n}{\lambda_{n}^{2}}]$ for some positive constants $C_{1}$ and $C_{2}$ such that $C_{1}<C_{2}$.

• (C.2)

Number of non-noisy observations (Z(\mbox{{\mathbf{s}}}_{i}) such that $i\in G_{n}$) in any unit square of the sampling region is bounded below by $C_{1}q_{e}n\lambda_{n}^{-2}$.

• (C.3)

The covariate process \mathbf{x}(\mbox{{\mathbf{s}}}) is such that, \mathbf{x}(\mbox{{\mathbf{s}}})=\mathbf{x}_{0}(\lambda_{n}^{-1}\mbox{{\mathbf{s}}}) where \mathbf{x}_{0}:[0,1]^{2}\to\mbox{{\mathbb{R}^{p}}} is a differentiable function with bounded partial derivatives in $(0,1)^{2}$ and \frac{1}{n}\sum_{i=1}^{n}\lVert\mathbf{x}(\mbox{{\mathbf{s}}}_{i})\rVert^{2}=O(1) for any set of \left\{\mbox{{\mathbf{s}}}_{1},\dots,\mbox{{\mathbf{s}}}_{n}\right\}\subset\mathcal{D}_{n} and for any $n\geq 2$.

• (C.4)

\left\{\epsilon(\mbox{{\mathbf{s}}}):\;\mbox{{\mathbf{s}}}\in\mbox{{\mathbb{R}^{2}}}\right\} is a second-order stationary random field such that, \int_{\mbox{{\mathbb{R}^{2}}}}\;\lvert\rho_{\epsilon}\left(\mbox{{\mathbf{h}}}\right)\rvert d\mbox{{\mathbf{h}}}<\infty; and for some $\kappa>0$, $E\lvert\epsilon(\mathbf{0})\rvert^{2+\kappa}<\infty$ and $\alpha_{\epsilon}(u,v)\leq Cu^{-\nu_{1}}v^{\nu_{2}}$ with $\nu_{1}>4/\kappa$ and $\nu_{2}>0$.

• (C.5)

$\delta_{n}^{-1}+\lambda_{n}^{-1}\delta_{n}+n^{-1/2}\lambda_{n}\delta_{n}^{-1}\to 0.$

• (C.6)

$\underset{n}{\text{sup}}\underset{i\in\left\{1,\dots,n\right\}}{\text{sup}}\lvert\xi^{p}_{n(i)}\rvert<\infty$ for all $p\in\left\{0.25,0.75\right\}$.

• (C.7)

For all $i\in\left\{1,\dots,n\right\}$, $j\in\left\{1,\dots,n(i)\right\}$ and $p\in\left\{0.25,0.5,0.75\right\}$, $F_{i_{j}}(\cdot)$ is absolutely continuous in an neighborhood of $\xi^{p}_{n(i)}$, i.e. $f_{i_{j}}(x)=\left(d/dx\right)F_{i_{j}}(x)$ exists in that neighborhood.

• (C.8)

For all $i\in\left\{1,\dots,n\right\},\;\;$ $0<\underset{j\in\left\{1,\dots,n(i)\right\}}{\text{inf}}f_{i_{j}}(\xi^{p}_{n(i)})\leq\underset{j\in\left\{1,\dots,n(i)\right\}}{\text{sup}}f_{i_{j}}(\xi^{p}_{n(i)})<\infty$.

• (C.9)

$F_{i}^{-1}(0.5)=0$ for all $i\in\left\{1,\dots,n\right\}$.

In a later remark we have discussed how the conditions would change if we have considered the Mean-VS estimator, \hat{\mbox{\bm{\beta}}}_{\text{vs}}^{(a)}. We now comment briefly on the conditions. (C.1) assumes that the sampling design has to be such that the number of observations in any unit square varies as the same rate as of $n/\lambda_{n}^{2}$, the ‘average’ number of observations per unit square. In (C.2) we restrict the number of noisy observations in any unit square by $C_{1}(1-q_{e})n\lambda_{n}^{-2}$ which combining with (C.1) specifies that at most $1-q_{e}$ proportion of the observations in an unit block can be noisy, i.e. corresponding noise variances (denoted by $\tau_{i}^{2}$) are non-zero. (C.3) puts regularity conditions on the covariates; it states that the covariate process is an ‘inflated’ version of a ‘smooth’ process \mbox{{\mathbf{x}}}_{0}(\mbox{{\mathbf{s}}}) in the prototype sampling region $[0,1]^{2}$, making \mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}}) a slowly varying process over the space. The other condition in (C.3) is standard in regression theory which puts a bound on the magnitude of the covariates as compared to the sample size. (C.4) states the required moment and mixing conditions on the residual process \left\{\epsilon(\mbox{{\mathbf{s}}})\right\}. Standard Gaussian process with ‘nicely’ decaying covariance functions, e.g. Expoenential, Matérn, Gaussian, etc. satisfies this condition. (C.5) makes sure that the size of the $\delta_{n}$-neighborhoods used to compute the VS increases with $n$ resulting in increasing number of observations for VS computation. But, at the same time the second term of condition (C.5) states that the block has to be small enough as compared to the whole sampling region retaining the ‘local’ feature of VS. The third term in (C.5) states that the rate at which the samples are infilling in the $\delta_{n}$-neighborhood has to go to $\infty$ as $n\to\infty$. Conditions (C.6) - (C.8) are standard assumptions needed for quantile consistency of non-i.i.d. random variables (for details, see Sen 1968 and Ghosh 1971). (C.9) is needed to make the w(\mbox{{\mathbf{s}}}_{i})’s marginal population median is equal to [math].

We now state the consistency results of this paper. The following two propositions are the key to prove the rest of the results of our paper. In the following proposition we provide an asymptotic approximation of the Median-VS.

Proposition 3.1**.**

Under conditions (C.1) to (C.9), for all $i\in\left\{1,2,\dots,n\right\}$,

[TABLE]

where \mathcal{I}_{n}(\mbox{{\mathbf{s}}}_{i})=\xi_{n(i)}^{0.75}-\xi_{n(i)}^{0.25}, i.e. the IQR of the distribution function $\bar{F}_{n(i)}(x)$ and $a_{n}$ is a sequence of positive reals $\downarrow 0$ given by, $a_{n}=n^{-1/2}\lambda_{n}\delta_{n}^{-1}+\delta_{n}\lambda_{n}^{-1}.$

The proof of this proposition is based on the extension of the Ghosh-Bahadur representation of sample quantiles for irregularly spaced spatial data (Lemma A.2 in Appendix A.1). Before we comment on the implication of Proposition 3.1, let us state the following proposition.

Proposition 3.2**.**

Under conditions (C.1) to (C.9), for all $i\in\left\{1,\dots,n\right\}$

[TABLE]

*where, $C_{\epsilon}^{(l)}(q_{e})=\sigma_{\epsilon}\left(F_{\epsilon}^{-1}(\text{max}\left\{1-(0.25/q_{e}),0\right\})-F_{\epsilon}^{-1}(\text{min}\left\{0.25/q_{e},1\right\})\right)$ and

$C_{\epsilon}^{(u)}(q_{e})=\sigma_{\epsilon}\left(F_{\epsilon}^{-1}(\text{min}\left\{0.75/q_{e},1\right\})-F_{\epsilon}^{-1}(\text{max}\left\{1-(0.75/q_{e}),0\right\})\right)$.*

From Proposition 3.2 we conclude that the scaling term, \mathcal{I}_{n}(\mbox{{\mathbf{s}}}_{i}), in the Median-VS approximation in Proposition 3.1 is bounded above by $C_{\epsilon}^{(u)}(q_{e})$, which is a finite constant if $q_{e}>0.75$, i.e. the proportion of noisy observations is less than $25\%$. Hence, under the assumption that $q_{e}>0.75$, from Proposition 3.1 we can see that if the noise variance associated with w(\mbox{{\mathbf{s}}}_{i}) – denoted by $\tau_{i}^{2}$ – is large, with higher probability \lvert w(\mbox{{\mathbf{s}}}_{i})\rvert will take larger values making \exp\left(-\lvert w(\mbox{{\mathbf{s}}}_{i})\rvert/\mathcal{I}_{n}(\mbox{{\mathbf{s}}}_{i})\right) (from here denoted by \tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{i})) – which is bounded above by \exp\left(-\lvert w(\mbox{{\mathbf{s}}}_{i})\rvert/C^{(u)}_{\epsilon}(q_{e})\right) – closer to [math]. Whereas, if w(\mbox{{\mathbf{s}}}_{i}) has less variance then with higher probability \lvert w(\mbox{{\mathbf{s}}}_{i})\rvert will take values in a neighborhood of [math] making \tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{i}) closer to $1$. Also, note that the lower bound $C^{(l)}_{\epsilon}(q_{e})>0$ if we have $q_{e}>0.5$. Hence, given that the poportion of good observations in the data is more than $75\%$, Proposition 3.1 along with Proposition 3.2 asymptotically justifies that Median-VS will give low scores to ‘bad’ observations and high scores to ‘good’ observations.

Next, we use the approximation in Proposition 3.1 to come up with a representation of the VS-based mean parameter estimator \hat{\mbox{\bm{\beta}}}^{(m)}_{\text{vs}} in the following theorem.

Theorem 3.1**.**

Under conditions (C.1) to (C.9),

[TABLE]

where \tilde{\mathbf{D}}_{v}=\text{diag}\left(\tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{1}),\dots,\tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{n})\right) and \mathbf{w}=\left(w(\mbox{{\mathbf{s}}}_{1}),\dots,w(\mbox{{\mathbf{s}}}_{n})\right)^{\prime}.

The proof of Theorem 3.1 is given in Appendix A.2.3. Now, we add some remarks regarding the representation of Median-VS-based estimator given in Equation 3.5.

Remark 3.1**.**

The representation in Equation 3.5 asymptotically approximates the estimation error in \hat{\mbox{\bm{\beta}}}^{(m)}_{\text{vs}} by the random variable given by \left(\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{D}_{v}\right)\mbox{{\mathbf{X}}}\right)^{-1}\left(\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\tilde{D}_{v}\mathbf{w}\right). Observe that, if w(\mbox{{\mathbf{s}}}_{i}) are marginally symmetric around [math], then E\left(\tilde{V}(\mbox{{\mathbf{s}}}_{i})w(\mbox{{\mathbf{s}}}_{i})\right)=0. Moreover, as proved in Section A.2.5, E\left(\tilde{V}(\mbox{{\mathbf{s}}}_{i})^{2}w(\mbox{{\mathbf{s}}}_{i})^{2}\right)\leq\exp(-2)(C_{\epsilon}^{(u)}(q_{e}))^{2}. Using this, in the next section we show that the MSE in VS-based estimation of the regression parameter is asymptotically not dependent on the high noise variances, i.e. $\tau_{i}^{2}$. On the contrary, the same is not true for the OLS estimator as, E\left(w(\mbox{{\mathbf{s}}}_{i})^{2}\right)=\sigma_{\epsilon}^{2}+\tau_{i}^{2}. This observation highlights the importance of the representation in Equation 3.5. We elaborately discuss the advantage of VS-based estimation in this regard in Section 3.4.

We use Theorem 3.1 to prove the consistency of \hat{\mbox{\bm{\beta}}}^{(m)}_{\text{vs}} in the following corollary. For that we need the a set of additional conditions as follows.

• (C.10)

\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}}\to\mbox{{\mathbf{C}}}_{X}\succ 0.

• (C.11)

$E\left(\tilde{D}_{v}\mathbf{w}\right)=\mathbf{0}$.

• (C.12)

The proportion of non-noisy observations, i.e. $q_{e}$ is strictly greater than $0.75$.

• (C.13)

For any real number $a>0$, \psi_{\epsilon}(a)=\int_{0}^{\infty}e^{-x}P\left(\lvert\epsilon(\mbox{{\mathbf{s}}})\rvert<ax\right)dx>0.

We now briefly comment on the above conditions. Condition (C.10) is standard assumption needed for the consistency of the least squares based mean parameter estimators. (C.11) is trivially valid if the w(\mbox{{\mathbf{s}}}_{i})’s are marginally symmetric around [math]. (C.12) is needed as to get a finite upper-bound on the scaling component of \tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{i}), i.e. I_{n}(\mbox{{\mathbf{s}}}_{i}), we need the proportion of noisy observations to be less than $0.25$. (C.13) is required to ensure the invertibility of the matrix \left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{{\mathbf{X}}}\right). This condition is trivially true when the marginal distribution of \epsilon(\mbox{{\mathbf{s}}}) has positive mass around a neighborhood of its center.

Corollary 3.1**.**

Under conditions (C.1) to (C.13), \hat{\mbox{\bm{\beta}}}_{\text{vs}}^{(m)}=\mbox{\bm{\beta}}+O_{p}\left(a_{n}\right).

Corollary 3.1 not only proves the consistency of the Median-VS-based estimator, but also provides an order of rate of convergence in probability for the Median-VS-based estimator. We state the proof of Corollary 3.1 in Appendix A.2.4. We now put some remarks on the consistency of the Mean-VS-based estimator.

Remark 3.2**.**

Similar consistency results of Mean-VS-based estimator \left(\hat{\mbox{\bm{\beta}}}_{\text{vs}}^{(a)}\right) can also be proved but that requires a set of different conditions than that are stated in (C.1) - (C.9). As the sample standard deviation is used in the definition of Mean-VS we need \frac{1}{n}\sum_{i=1}^{n}\lVert\mathbf{x}(\mbox{{\mathbf{s}}}_{i})\rVert^{4}=O(1) as well as all the error components including \epsilon(\mbox{{\mathbf{s}}}) must have finite fourth order moments. In addition, the mixing condition in (C.4) changes as we would require $\epsilon$-process to be fourth-order stationary and all the correlations of the form \text{cor}\left(\epsilon^{a}(\cdot),\epsilon^{b}(\cdot+\mbox{{\mathbf{h}}})\right), for all $a,b\in\left\{1,\dots,4\right\}$, have to be integrable over the space. Conditions (C.5) - (C.9) are not needed as those are specific to the existence of the population quantiles used in the results of the Median-VS-based estimator.

Though the consistency of the VS-based estimators has been proved in this section, that does not necessarily point out the advantages of the VS-based estimation over the standard methods, e.g. the OLS, when it comes to analyzing noisy spatial data. In the next subsection, we investigate the asymptotic mean squared errors of the VS-based as well as the OLS estimator. For the next section we only consider the Median-VS-based regression estimator (\hat{\mbox{\bm{\beta}}}^{(m)}_{\text{vs}}) and hence, for simplicity of nations, we denote it by \hat{\mbox{\bm{\beta}}}_{\text{vs}}.

3.4 Asymptotic efficiency of VS-based estimators

From Theorem 3.1, we can write \hat{\mbox{\bm{\beta}}}_{\text{vs}}-\mbox{\bm{\beta}}=\mathbf{l}_{n}^{\text{vs}}+O_{p}(a_{n}), where the leading error term in the VS-based estimation is given by, \mathbf{l}_{n}^{\text{vs}}=\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{{\mathbf{X}}}\right)^{-1}\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}\right). On the other hand, due to linearity of the OLS estimator, \hat{\mbox{\bm{\beta}}}_{\text{ols}}-\mbox{\bm{\beta}}=\mathbf{l}_{n}^{\text{ols}}, where it can be easily shown that \mathbf{l}_{n}^{\text{ols}}=\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}}\right)^{-1}\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}\mathbf{w}\right). Before stating the results of this section, we need to introduce a few more notations. It can be shown easily that E\left(\lVert\mathbf{l}_{n}^{\text{ols}}\rVert^{2}\right)=n^{-1}\text{tr}\left(\mbox{{\mathbf{H}}}_{n}\bm{\Sigma}_{w}\right) (details are given in Appendix A.2.6), where \mbox{{\mathbf{H}}}_{n}=n^{-1}\mathbf{X}(n^{-1}\mathbf{X}^{\prime}\mathbf{X})^{-2}\mathbf{X}^{\prime}, $\bm{\Sigma}_{w}=\mathrm{Var}\left(\mathbf{w}\right)$. For any symmetric matrix ${\mathbf{A}}$, by \lambda_{\text{min}}\left(\mbox{{\mathbf{A}}}\right), we denote the minimum eigenvalue of ${\mathbf{A}}$. Let us introduce the following notations: M_{X}=\underset{\mbox{{\mathbf{s}}}\in\mathcal{D}_{n}}{\text{sup}}\underset{1\leq j\leq p}{\text{max}}\;x^{2}_{j}(\mbox{{\mathbf{s}}}); C^{(0)}=2^{-1}\text{tr}\left(\mbox{{\mathbf{C}}}_{X}^{-1}\right); C^{(1)}=2p\left(\lambda_{\text{min}}\left(\mbox{{\mathbf{C}}}_{X}\right)\right)^{-2}M_{X}; and finally \mathcal{M}_{n}(\varepsilon)=\left\{1\leq i\leq n:\left(\mbox{{\mathbf{H}}}_{n}\right)_{ii}>n^{-1}\varepsilon\right\}, where by \left(\mbox{{\mathbf{A}}}\right)_{ij} we denote the $ij$-th element of the matrix ${\mathbf{A}}$. Now, we have the required notation to state the main results of this section.

Theorem 3.2**.**

Let assumptions (C.1) to (C.13) hold true. Then, for positive constants $C_{1}(q_{e})$, $C_{2}$ and $C_{3}$,

[TABLE]

where $C_{1}(q_{e})=\left(C_{\epsilon}^{(u)}(q_{e})/\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)\right)$ and $C_{2}$, $C_{3}$ are independent of the noise variance parameters ($\tau_{i}^{2}$ for $i\in G_{n}$, and $q_{e}$).

Theorem 3.3**.**

Let assumptions (C.1) to (C.4) and (C.10) hold true. Then, for any $0<\varepsilon<\text{min}\left(C^{(0)},C^{(1)}\right)$ and large enough $n$,

[TABLE]

with $\lvert\mathcal{M}_{n}(\varepsilon)\rvert\geq\frac{C^{(0)}-\varepsilon}{C^{(1)}-\varepsilon}n$. Here $C_{1}$ is a positive constant not dependent on the error variance parameters ($\tau_{i}^{2}$ for $i\in G_{n}$, and $q_{e}$).

Corollary 3.2**.**

Let assumptions (C.1) - (C.4) and (C.10) hold true. In addition, assume \text{inf}_{n}\;\underset{1\leq i\leq n}{\text{min}}\lVert\mbox{{\mathbf{X}}}[i,]\rVert^{2}>0. Then, there exists an $\varepsilon>0$ such that for large enough $n$,

[TABLE]

for some positive constant $C_{1}$ not dependent on the noise model parameters.

The proofs of Theorem 3.2, Theorem 3.3 and Corollary 3.2 are stated in Appendix A.2.5, A.2.6 and A.2.7 respectively. Next, we put some remarks on the above results.

Remark 3.3**.**

From Theorem 3.2 we see that the MSE of the leading error term in the Median-VS-based regression estimation can be bounded above by some constant $C_{3}$ times $n^{-1}\left(C_{1}(q_{e})\right)^{2}+C_{2}\lambda_{n}^{-4}$, which goes to [math] as $n\to\infty$, and is dependent on the noise model parameters only through the proportion of the good observations $q_{e}$. Clearly, the MSE of the Median-VS-based regression estimator is not dependent on the noise variances (denoted by $\tau_{i}^{2}$) associated with the ‘bad’ observations. Note that, $\psi_{\epsilon}(\cdot)$ is an increasing function. We have already mentioned that as $q_{e}$ increase, $C_{\epsilon}^{(u)}(q_{e})$ decreases and $C_{\epsilon}^{(l)}(q_{e})$ increases, i.e. $C_{1}(q_{e})$ is a decreasing function of $q_{e}$. Hence, if the proportion of good observations increases, the upper bound in Equation 3.6 decreases, restricting the range of $E\left(\lVert\mathbf{l}_{n}^{\text{vs}}\rVert^{2}\right)$ from above. The other term, $C_{2}\lambda_{n}^{-4}$ has hardly any effect on $E\left(\lVert\mathbf{l}_{n}^{\text{vs}}\rVert^{2}\right)$ as $\lambda_{n}^{-4}$ is of much smaller order as compared to $n^{-1}$.

Remark 3.4**.**

The lower-bound for the MSE of the OLS estimator, as we can see in Theorem 3.3 and Corollary 3.2, is dependent on the additional noise variances of the corrupted observations through the summation $n^{-1}\sum_{i\in G_{n}^{c}\cap\mathcal{M}_{n}(\varepsilon)}\tau_{i}^{2}$ and $n^{-1}\sum_{i\in G_{n}^{c}}\tau_{i}^{2}$ respectively. Clearly, if the noise variances corresponding to some of the observations are high, efficiency of the OLS estimator will be significantly affected. For example, consider the case where $\tau_{i}^{2}\approx n^{c}$ for some constant $c>0$. Then for a given $q_{e}$, the lower-bound in Equation 3.8, is of order $n^{c-1}$. In Table 1 we show the advantage of using VS-based estimator as opposed to the OLS one in terms of the asymptotic order of relative efficiency.

Here, we refer the order of the lower bound of the MSE of the OLS estimator by OLS-LB, the upper bound for the MSE of the VS-based estimator by VS-UB and, the lower bound for the relative efficiency of VS-based estimator to the OLS one (defines as MSE(VS)/MSE(OLS)) by RE-LB. Table 1 clearly portrays the advantage of the VS-based estimation as compared to the OLS one: if the order of the noise variances associated with even a ‘small’ portion of the observations is ‘high’ , due to the robustness of the VS-based estimator against the noisy observations, there is a significant gain in efficiency relative to the OLS.

4 Simulation Study

In this section, we summarize our efforts to evaluate the VS-based estimator numerically. We first describe the simulation designs and then, we state the results and draw inference from them.

4.1 Simulation Setup

We use the following spatial regression model to simulate the realizations from the the ‘true’ random field \left\{Y(\mbox{{\mathbf{s}}})\right\}:

[TABLE]

where, \mbox{\bm{\beta}}=\left(\beta_{0},\beta_{x},\beta_{y},\beta_{h}\right)^{\prime} is the vector of regression parameters; h(\mbox{{\mathbf{s}}}) is a deterministic function of location ${\mathbf{s}}$; and \left\{\epsilon(\mbox{{\mathbf{s}}})\right\} is a mean zero second-order stationary spatially correlated process. To define the function h(\mbox{{\mathbf{s}}}) over the sampling region, we use the deterministic function h(\mbox{{\mathbf{s}}})=H_{1}\cdot\sum_{j=1}^{H_{2}}w_{h}(j)f(\mbox{{\mathbf{s}}};\;\bm{\mu}_{j},\Sigma_{j})+H_{3}, where $f(\cdot;\bm{\mu},\Sigma)$ denotes the bivariate normal density with mean $\bm{\mu}$ and covariance matrix $\Sigma$ and $\left\{\left(\bm{\mu}_{j},\Sigma_{j}\right):1\leq j\leq H_{2}\right\}$ is a fixed set of vectors and matrices. The choice of this function is motivated from the elevation function used in the simulations of Chakraborty et al. (2019). The residual vector \left(\epsilon(\mbox{{\mathbf{s}}}_{1}),\dots,\epsilon(\mbox{{\mathbf{s}}}_{n})\right)^{\prime} are sampled from a mean zero Gaussian process with the Exponential covariance function given as C^{\text{exp}}(\mbox{{\mathbf{h}}};\bm{\theta})=\sigma_{\epsilon}^{2}\exp\left(\frac{\lVert\mbox{{\mathbf{h}}}\rVert}{\rho}\right)+\tau^{2}\mathbb{I}(d=0) where, $\mathbb{I}(\cdot)$ denotes the indicator function. The covariance parameter vector of interest is $\bm{\theta}=\left(\tau^{2},\sigma_{\epsilon}^{2},\rho\right)^{\prime}$, where $\tau^{2}$ is the nugget effect, $\sigma_{\epsilon}^{2},\rho$ are the partial sill and range parameters respectively (for details see Haskard 2007; page 37 Gelfand et al. 2010).

To include noise in the varying-quality observations \left(Z(\mbox{{\mathbf{s}}}_{1}),\dots,Z(\mbox{{\mathbf{s}}}_{n})\right)^{\prime} using the additive-multiplicative noise structure given in Equation 2.3, we use the following distributions for the additive and multiplicative components. Given a set of good observations $G_{n}$, for $i\in G_{n}$, the $\epsilon_{M_{i}}$ and $\epsilon_{A_{i}}$ are identically equal to $1$ and [math], respectively. For $i\in G_{n}^{c}$, we take

[TABLE]

Here, the variance associated with the multiplicative component of a noisy observation is given by $\sigma_{M}^{2}=\frac{1}{2\alpha_{M}+1}$. To select the set of ‘good’ observations, i.e. $G_{n}$, we have arbitrarily selected $\lfloor q_{e}n\rfloor$ many (fixed) indices among $\left\{1,\dots,n\right\}$ prior the simulations. Note that, the choice of multiplicative error distribution in 4.2 restricts its realizations to be in $[0,2]$, and also, ensures that the multiplicative errors are symmetric around 1.

Estimation of the regression parameters through VS-based and the OLS method do not require the covariance parameters to be estimated. Whereas, the GLS estimator uses the covariance information which is often unknown in real data analysis. Hence, the covariance parameters are estimated from the de-trended spatial observations and then the estimated covariance parameters are used to compute GLS estimator. Then the updated mean estimates are used to compute the covariance parameters and we continue until the parameter estimates converge. This estimator is called the estimated GLS (referred as EGLS, for details see page 23 Cressie 1993). In the following section, we compare the performance of VS-based (Mean-VS, Median-VS) and commonly used least squares based (OLS and EGLS) techniques for estimation of the regression parameters.

Though the theoretical analysis of the VS-based covariance parameter estimation is beyond the scope of this paper, in our simulation studies we evaluate the accuracy of the same to capture the dependence structure of the spatial process. The reasons are following. First, the boxplots of the regression parameter estimates in Section B.1 in the supplementary material suggest that the robust approach (robust regression with Huber’s loss, Huber and Ronchetti 2009) has similar performance as compared to the VS-based methods. But, extension of the robust approach for estimation of the covariance parameters is not straightforward. Whereas, the VS computed from the data can be incorporated in the covariance analysis to reduce the effects of ‘bad’ observations in the estimation. Hence, comparison of the VS-based technique of covariance estimation needs to be evaluated as compared to the other robust techniques of covariance parameter estimations. We consider the following two methods in that regard: the robust WLS based variogram model fitting by Cressie and Douglas (1980) (referred as WLS in Table 3) and the robust-REML technique proposed by Künsch et al. (2011) and implemented by Papritz (2018a). Second, we want to justify the use of VS-based covariance analysis in the real data example in Section 5.

For all the simulation results in this section we take \mbox{\bm{\beta}}=\left(70,5,-2,-0.05\right)^{\prime} and \mbox{\bm{\theta}}=\left(0,6,1\right)^{\prime}. To investigate the performances of the different types of estimators with varying noise level, we consider the following set of values of the noise parameters: (a) $\sigma_{A}\in\left\{5,50,100\right\}$, (b) $\alpha_{M}\in\left\{2,0.5,0.05\right\}$, and (c) $q_{e}\in\left\{0.95,0.85,0.75\right\}$. The set of values for each of the above parameters are written in the order of increasing noise variance and proportion. For example, if the first values of all the three sets are considered then the additional noise variance of a corrupted observation at location ${\mathbf{s}}$ is 0.2\left(\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}})^{\prime}\mbox{\bm{\beta}}\right)^{2}+26.2, and the proportion of such observations is $5\%$; whereas the same variance will be 0.91\left(\mbox{{\mathbf{x}}}(\mbox{{\mathbf{s}}})^{\prime}\mbox{\bm{\beta}}\right)^{2}+10005.46 if the last values of the sets are considered and the proportion of such noisy observations will be $25\%$. In the tables in the next section, when one of the noise model parameters are varying, the other parameters are fixed at the lowest noise level, i.e. if $\sigma_{A}$ is varying, then $\alpha_{M}$ is fixed at $2$ and $q_{e}$ is fixed at $0.95$.

4.2 Simulation Results

Table 2 shows the empirical mean square errors of the considered mean parameter estimators. In the final column of Table 2, the relative efficiency (RE) of the Median-VS-based regression estimator with respect to the OLS one is reported to emphasize the superiority of VS-based estimation in the analysis of noisy spatial data.

It is evident from Table 2 that the VS-based estimation has better performance than the standard OLS and EGLS methods across all the cases. The increasing noise variance has very little effect on the performance of the VS-based estimation in terms of MSE. For instance, consider the case when $n=500$. If the standard deviation of the associated additive noise increases from $5$ to $50$, the MSE of Median-VS-based estimator has hardly changed, whereas the MSE of the EGLS and OLS estimators have increased by approximately $164.4\%$ and $47.8\%$, in the respective order. The consistency of the VS-based estimator is also prominent from Table 2 – in each of the scenarios, as sample size increases, the MSE is going close to [math]. The relative efficiency column justifies the superiority of the VS-based estimator to the least squares-based estimators. For example, when $\sigma_{A}=100$ and the sample size $n=1000$, i.e. very high noise variance is associated with $5\%$ of the observations – to achieve same accuracy as of the Median-VS-based estimators the OLS and the EGLS methods need approximately $8$ and $14.5$ times more observations. Also, given a sample size, the relative efficiency increases if the magnitude of the noise variances associated with the ‘bad’ observations increase (controlled by $\sigma_{M}$ and $\sigma_{A}$) or the proportion of the ‘bad’ observations increases. Note that, the advantage of the VS-based estimation for varying proportion of the non-noisy data (i.e. $q_{e}$) is not as ‘good’ as in the case of varying additive noise variance $\sigma_{A}$. The reason is that the effect of the increasing noise variances on the MSE of the VS-based estimators is much less as cmpared to the increase of proportion of noisy observations, which we have established theoretically in Section 3.4.

Observe that, in Table 2, the accuracy is hampered severely if instead of the OLS estimator EGLS is used to estimate the regression parameter. One of the possible reasons behind this is the inefficient estimation of the covariance parameters using robust weighted least squares technique.

Hence, in Table 3, we have investigated the performance of the VS-based covariance estimation and compared it with that of the standard approach: robust weighted least squares (WLS) based variogram estimator. Clearly, the WLS approach has failed completely to estimate the covariance parameters – partial sill (or psill, denoted by $\sigma_{\epsilon}^{2}$) and range (denoted by $\rho$). Due to the added noise variance in the data the optimization in the variogram fitting in WLS method is not stable. On the other hand, VS-based variogram estimation has reasonably good accuracy, given that the proportion of the noisy observations are not too large. Also, in most of the cases, if the noise model parameters are held fixed, the MSE of the VS-based estimators decreases, as the sample size $n$ increase, which justifies the consistency of the VS-based variogram model fitting technique. But, in the case where the proportion of noisy observations is large (e.g. $q_{e}=0.8$, i.e. $20\%$ of the observations are noisy) the accuracy of the VS-based estimation is hampered. Clearly, the accuracy of the VS-based covariance estimation is more sensitive to the proportion of noisy observations in the data, rather than the magnitude of the noise variances – the same fact that we have established theoretically for the VS-based mean estimation in Section 3.4.

Next, in Table 4, we briefly report the results of our comparative analysis of VS-based technique with another robust covariance estimation method in geostatistics – robust REML proposed by Künsch et al. (2011) and implemented in R-package georob by Papritz (2018a) (referred as GR in Table 4).

The results are surprisingly in favor of the VS-based methodology. Clearly, the robust REML technique has completely failed to provide reasonable covariance estimates. We can not provide the results for all other choices of $n$ and noise model parameters because, the REML estimation in georob has huge time complexity. For example, for $\sigma_{A}=5$ and $n=1000$, on an average georob takes $145$ sec., whereas, VS-based estiamtion takes $2.5$ sec. – i.e. georob takes on an average $58$ times more. Hence, both in terms of computational complexity and the accuracy of estimation, VS-based technique outperforms robust-REML method (as implemented in Papritz 2018a).

In the next section, we are going to apply VS-based estimation technique on the coal ash data (Gomez and Hazen 1970) and compare our results with that of Papritz (2018b), where the author has implemented robust REML estimation on the same data set.

5 Example: Coalash data

In this section we consider the coal ash data (Gomez and Hazen 1970), which reports measurement (in $\%$ mass) of ash present in coal seam in Pennsylvania. The data set is available in R-package gstat (Pebesma 2004) and has been previously used to demonstrate robust geostatistical techniques by Cressie (1993) and Papritz (2018b). Let us denote the target process as \left\{Y(\mbox{{\mathbf{s}}})\right\} where Y(\mbox{{\mathbf{s}}}) is the mass percentage (often referred as $\%$* mass*) of ash in the coal seam at location ${\mathbf{s}}$. The observations plotted in Figure 1(a) are denoted by \left\{Z(\mbox{{\mathbf{s}}}_{1}),\dots,Z(\mbox{{\mathbf{s}}}_{n})\right\} where $n=208$. Some of these observations are corrupted due to measurement errors or other types of unknown noises.

We first compute the VS (Median-VS) of the observations and plot it spatially in Figure 1(b).

Cressie (1993) identified a set of observations as outliers by investigating the deviation of the observations from their overall sample median. VS uses the similar idea – instead of using the overall median, it assigns a score to each of the observations according to its deviation from the ‘local’ median relative to the ‘local’ variation. The observations that were identified as outliers by Cressie (1993) and Papritz 2018b receive VS less than 0.18 (in a scale of $(0,1]$). As both in theoretical analysis and simulation studies we have established that VS-based estimation is hampered if proportion of noisy observations is significantly large, we check the proportion of observations which have VS less than $0.22$ – which implies the deviation from local median is at least $1.51\times$local-IQR. We find that approximately $12\%$ of the observations have ‘high’ noise associated with it.

Once the VS of the observations has been computed we incorporate those score to estimate the mean and covariance structure of the process robustly. Through some explanatory data analysis similar to Papritz (2018b), we consider the following spatial regression model for the target process \left\{Y(\mbox{{\mathbf{s}}})\right\},

[TABLE]

where \mbox{\bm{\beta}}=\left(\beta_{0},\beta_{x}\right)^{\prime} is the regression parameter, $s_{x}$ is longitude (i.e. \mbox{{\mathbf{s}}}=\left(s_{x},s_{y}\right)^{\prime}) and \left\{\epsilon(\mbox{{\mathbf{s}}})\right\} is a zero mean second order stationary process with Matérn covariance and the covariance parameter vector is given by \mbox{\bm{\theta}}=\left(\sigma_{\epsilon}^{2},\rho,\eta^{2},\nu\right)^{\prime} where $\sigma_{\epsilon}^{2}$ is the partiall sill, $\rho$ is the range, $\eta^{2}$ is the nugget and $\nu$ is the smoothness parameter (for details, see Haskard 2007; Gelfand et al. 2010).

VS-based regression yields an estimate of \hat{\mbox{\bm{\beta}}}_{\text{vs}}=\left(11.071,-0.188\right)^{\prime} for the regression parameter – which suggests that the ash proportion decreases as we go from west to east.

Next, the VS-based smoothing is implemented on the observed residuals from the VS-based regression and the effect of smoothing has been displayed through frequency plots in Figure 2. Finally, the VS-based smooth residuals are used to estimate the covariance parameters through VS-based variogram model fitting. We summarize the details of VS-based covariance estimation in Figure 3 and Table 5. Observe that, the range parameter estiamte ($0.486$) in Table 5 is significantly greater than [math], i.e. the VS-based analysis can capture the spatial correlation in the residual process, as opposed to the pure nugget model found in the analysis of Papritz (2018b).

Parameters Estiamtes

partial sill ($\sigma_{\epsilon}^{2}$) 0.219

range ($\rho$) 0.486

nugget ($\eta^{2}$) 0.021

smoothness ($\kappa$) 5.00

Table 5: Estimated Matérn parameters.

Next, we compare the VS-based analysis with the robust REML approach. To implement the robust REML methodology on the coal ash data the R-package georob (Papritz 2018a) is used. We use leave-one-out cross-validation technique to compare the two approaches: for each of the observations in the coal ash data, we consider it as the test data and try to predict (kriging) it using all other observations without the test one. As, trying to predict a ‘bad’ observation may lead to erroneous inference, from the test data set we exclude the observations marked as outliers by Papritz (2018b) (count $11$) in their analysis. But, the training data contains all the data points. Instead of using VS, we use the analysis done by Papritz (2018b) to remove the outliers from the test data so that the choice of test data in the cross-validated analysis is not biased towards the VS-based approach.

In Figure 4, we have summarized the results of comparative analysis between VS-based and robust-REML approach.

From the empirical distribution function plot in Figure 4(a) it is evident that the prediction errors corresponding to the robust REML predictor has higher mass in the tails than the VS-based predictions. Moreover, as we can see in Figure 4(b), the margins of prediction errors (ME), that is the half of the lengths of the prediction intervals, for the VS-based approach is either similar or smaller as compared to the robust REML for most of the test data points. To check the credibility of the VS-based prediction intervals we find that $98.5\%$ of the test data points were inside the $95\%$ VS-based prediction intervals. In terms of cross-validated mean squared prediction error (MSPE), the VS-based MSPE is $0.716$ whereas the robust-REML prediction has an MSPE of $0.743$, i.e. there is nearly $4\%$ gain in efficiency when VS-based analysis is carrired out in place of the robust REML approach. Moreover, in terms of the computational time VS-based approach has a huge advantage as compared to the robust REML, as shown in Table 6.

Hence, for the analysis of larger spatial data sets, VS-based analysis has advantage over the robust REML approach.

6 Conclusion

In this paper, we have investigated the large sample behavior of the veracity scoring technique proposed by Chakraborty et al. (2019) for irregularly spaced spatial data when there is no reference data available. Through the asymptotic approximation of the VS we have established that VS of an observation is expected to take ‘small’ value (close to [math]) if the noise variance associated with the observation is high and, it is going to take larger values if the associated noise variance is less. Moreover, under the mixed-increasing domain spatial asymptotic framework and for a non-stationary noise model, we have sheded light upon the circumstances under which the weighted least squares regression estimator with veracity scores as the corresponding weights – which are data-driven and hence, spatially correlated – is asymptotically consistent. We have also provided an order on the rate of convergence in probability for the VS-based regression estimator. Next, we have established that the asymptotic mean squared error of the VS-based estimation is highly resistant to the magnitudes of the noise variances associated with the corrupted observations. But, the accuracy of the estimation can be hampered if the proportion of noisy observations increases. In addition, we have also considered the asymptotic MSE of the ordinary least squares estimator and showed how it depends on both the proportion of noisy observations as well as the magnitudes of the noise variances. We have provided empirical justifications of the advantages of the VS-based analysis through rigorous simulations. Finally, we have implemented the VS-based approach on the coal ash data and compared the results with that of robust REML technique – which shows that the VS-based approach has advantage both in terms of prediction accuracy and computational time.

There are several future directions to this work. First, the theoretical properties of the VS-based covariance estimation and kriging under similar spatial asymptotic framework can be explored. Second, a version of VS for spatio-temporal data will be of high interest, especially for analyzing real-time noisy crowdsourced data. Third, the results of this paper can be extended for spatial processes on ${\mathbb{R}^{d}}$ with $d\geq 3$. That way the proerties of the VS and the VS-based estimation will be extended for spatio-temporal data as well.

Acknowledgements

This research is partially funded by National Science Foundation (NSF) grant DMS-1613192. The authors would like to thank Dr. Alyson Wilson for bringing the importance of this work to authors’ attention and also for her insightful remarks during the course of this work.

Appendix A Proofs

To prove the main results of this paper we need to state couple of lemmas in the next subsection.

A.1 Auxiliary lemmas

Let us define $\hat{G}_{n(i)}(x)=1-\hat{F}_{n(i)}(x)$ and $\bar{G}_{n(i)}(x)=1-\bar{F}_{n(i)}(x)$.

Lemma A.1**.**

The strong-mixing coefficient $\alpha_{w}\left(\cdot,\cdot\right)$ of \left\{w(\mbox{{\mathbf{s}}})\right\} is bounded above by the strong-mixing coefficient $\alpha_{\epsilon}\left(\cdot,\cdot\right)$ of \left\{\epsilon(\mbox{{\mathbf{s}}})\right\}, i.e.

[TABLE]

for any $u>0$ , $v>0$.

Lemma A.2**.**

For any $0<p<1$,

[TABLE]

where $\bar{f}_{n(i)}(x)=\left(\sum_{j=1}^{n(i)}f_{i_{j}}(x)\right)/n(i)$ and $\sqrt{n(i)}R_{n(i)}\to 0$ in probability.

Clearly, the representation in Equation A.1 is the Ghosh-Bahadur representation in our setting. We extended the proofs of Sen (1968) and Ghosh (1971) for heteroscedastic $\alpha$-mixing sequence of random variables. The proofs of these two lemmas have been provided in Section A.1 of the supplementary material.

A.2 Proofs of the main results

We start by stating the proof of Proposition 3.1.

A.2.1 Proof of Proposition 3.1

Proof.

Note that,

[TABLE]

as, from condition (C.3),

[TABLE]

for some \mbox{{\mathbf{s}}}^{*}_{i_{j}}\in\mbox{{\mathbb{R}^{2}}} in between \mbox{{\mathbf{s}}}_{i} and \mbox{{\mathbf{s}}}_{i_{j}} component wise. Now, from Lemma A.2 and condition (C.9), we have $Q_{2}(\mathbf{w}_{i})=O_{p}\left((n(i))^{-1/2}\right)=O_{p}(n^{-1/2}\lambda_{n}\delta_{n}^{-1})$. So, from equation A.2, we have Z(\mbox{{\mathbf{s}}}_{i})-Q_{2}(\mathbf{Z}_{i})=w(\mbox{{\mathbf{s}}}_{i})+O_{p}(n^{-1/2}\lambda_{n}\delta_{n}^{-1}+\delta_{n}\lambda_{n}^{-1}). Similar argument using the first and third sample quartile we can show that \text{IQR}(\mathbf{Z}_{i})=\mathcal{I}_{n}(\mbox{{\mathbf{s}}}_{i})+O_{p}(n^{-1/2}\lambda_{n}\delta_{n}^{-1}+\delta_{n}\lambda_{n}^{-1}). Now, through simple probability calculations we have

[TABLE]

and hence proved. ∎

A.2.2 Proof of Proposition 3.2

Proof.

By definition of $\xi_{n(i)}^{p}$, $\bar{F}_{n(i)}(\xi_{n(i)}^{p})=p$. Now, observe that, $\frac{1}{n(i)}\sum_{j=1}^{n(i)}F_{i_{j}}(\xi_{n(i)}^{p})=\frac{1}{n(i)}\left\{\sum_{j\in G_{i,n}}F_{i_{j}}(\xi_{n(i)}^{p})+\sum_{j\in G_{i,n}^{c}}F_{i_{j}}(\xi_{n(i)}^{p})\right\}$, where $G_{i,n}\subset G_{n}$ is the set of indices of the observations in the square neighborhood around \mbox{{\mathbf{s}}}_{i} which are not affected by the additive-multiplicative noise. From condition (C.1) and (C.2), $\lvert G_{i,n}\rvert/n(i)\geq q_{e}$. Hence,

[TABLE]

which proves that $\xi_{n(i)}^{p}\leq\sigma_{\epsilon}\left(F_{\epsilon}^{-1}(\text{min}\left\{p/q_{e},1\right\})\right)$. In place of $p$, using the similar argument for the $(1-p)$-th quantile we can show that $\xi_{n(i)}^{p}\geq\sigma_{\epsilon}\left(F_{\epsilon}^{-1}(\text{max}\left\{1-\frac{1-p}{q_{e}},0\right\})\right)$. Taking $p=0.75$ and $0.25$ we have the proof. ∎

A.2.3 Proof of Theorem 3.1

Proof.

From definition of Median-VS-based regression estimator,

[TABLE]

where \mathbf{D}_{v}^{(m)}=\text{diag}\left(V^{(m)}(\mbox{{\mathbf{s}}}_{1}),\dots,V^{(m)}(\mbox{{\mathbf{s}}}_{n})\right). For simplicity of notation, as we are only considering Median-VS here, let us denote $\mathbf{D}_{v}^{(m)}$ as $\mathbf{D}_{v}$. Now, from Proposition 3.1, \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mathbf{D}_{v}\mathbf{w}=\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}+O_{p}(a_{n}) and \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mathbf{D}_{v}\mbox{{\mathbf{X}}}=\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{{\mathbf{X}}}+\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\left(\tilde{\mathbf{D}}_{v}-E\left(\tilde{\mathbf{D}}_{v}\right)\right)\mbox{{\mathbf{X}}}+O_{p}(a_{n}). Next we will show that \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\left(\tilde{\mathbf{D}}_{v}-E\left(\tilde{\mathbf{D}}_{v}\right)\right)\mbox{{\mathbf{X}}}=O_{p}(\lambda_{n}^{-1}).

Notice that \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\left(\tilde{\mathbf{D}}_{v}-E\left(\tilde{\mathbf{D}}_{v}\right)\right)\mbox{{\mathbf{X}}} is a $p\times p$ matrix and hence enough to show that a general element of of the matrix say $(k,l)$-th element is $O_{p}(\lambda_{n}^{-1})$. Now, from simple calculations,

[TABLE]

where, T_{kl}(\mbox{{\mathbf{s}}}_{i})=x_{k}\left(\mbox{{\mathbf{s}}}_{i}\right)x_{l}\left(\mbox{{\mathbf{s}}}_{i}\right)\left\{\tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{i})-E\left(\tilde{V}^{(m)}(\mbox{{\mathbf{s}}}_{i})\right)\right\}. As, E\left(T_{kl}(\mbox{{\mathbf{s}}}_{i})\right)=0 enough to show that \mathrm{Var}\left(\frac{1}{n}\sum_{i}T_{kl}(\mbox{{\mathbf{s}}}_{i})\right)=O(\lambda_{n}^{-2}). To do so we will adopt the similar strategy used in the proof of Lemma A.2. We will divide the whole sampling region $\mathcal{D}_{n}=\lambda_{n}[0,1]^{2}$ into disjoint unit squares $\left\{B_{1},\dots,B_{b_{n}}\right\}$ and $b_{n}$ is the minimal natural number such that $\mathcal{D}_{n}\subset\cup_{m=1}^{b_{n}}B_{m}$. Clearly $b_{n}=\lceil\lambda_{n}\rceil^{2}$ and we can rewrite the set of unit squares as $\left\{B_{m}:\;m\in\left\{1,\dots,b_{n}\right\}\right\}=\left\{B_{m_{1},m_{2}}:\;m_{1},m_{2}\in\left\{1,\dots,\sqrt{b_{n}}=\lceil\lambda_{n}\rceil\right\}\right\}$, i.e. by $m_{1}$ and $m_{2}$ we are indexing the squares row and column wise respectively. Then following similar arguments we can show that E\left(\frac{1}{n}\sum_{i}T_{kl}(\mbox{{\mathbf{s}}}_{i})\right)^{2}=O(\lambda_{n}^{-2}). Hence, \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mathbf{D}_{v}\mbox{{\mathbf{X}}}=\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{{\mathbf{X}}}+O_{p}(\lambda_{n}^{-1}+a_{n})=\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mbox{{\mathbf{X}}}+O_{p}(a_{n}). So,

[TABLE]

and hence proved. ∎

A.2.4 Proof of Corollary 3.1

Proof.

From condition (C.11), $\frac{1}{n}\mathbf{X}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}=O_{p}\left(\sqrt{\mathrm{Var}\left(\frac{1}{n}\mathbf{X}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}\right)}\right)$. Using similar argument as in the proof of Theorem 3.1, we can show that $\mathrm{Var}\left(\frac{1}{n}\mathbf{X}^{\prime}\tilde{\mathbf{D}}_{v}\mathbf{w}\right)=O(\lambda_{n}^{-2})$. Hence, \hat{\mbox{\bm{\beta}}}^{(a)}_{\text{VS}}=\mbox{\bm{\beta}}+\left(\frac{1}{n}\mathbf{X}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mathbf{X}\right)^{-1}O_{p}\left(\lambda_{n}^{-1}\right)+O_{p}(a_{n}). Now, for all $i\in\left\{1,\dots,n\right\}$,

[TABLE]

Clearly, as from condition (C.12) $q_{e}>0.75$, $C_{\epsilon}^{(l)}(q_{e})>0$ and under assumption (C.13) $\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)$ is strictly positive. Note that, $\frac{1}{n}\mathbf{X}^{\prime}\mathbf{X}\succeq\frac{1}{n}\mathbf{X}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mathbf{X}\succ\frac{\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)}{n}\mathbf{X}^{\prime}\mathbf{X}$ where, from condition (C.10), \frac{1}{n}\mathbf{X}^{\prime}\mathbf{X}\to\mbox{{\mathbf{C}}}_{X} for some positive definite matrix \mbox{{\mathbf{C}}}_{X}. Clearly, for large enough $n$, \left(\frac{1}{n}\mathbf{X}^{\prime}E\left(\tilde{\mathbf{D}}_{v}\right)\mathbf{X}\right)^{-1}\preceq\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)\mbox{{\mathbf{C}}}_{X} and hence we have the proof. ∎

A.2.5 Proof of Theorem 3.2

Proof.

Observe that,

[TABLE]

We can show that \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}}\succeq\mathbf{B}_{n}\succ\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)\frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}} where $\psi_{\epsilon}\left(C_{\epsilon}^{(l)}(q_{e})\right)>0$ and from condition (C.10) \frac{1}{n}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}}\to\mbox{{\mathbf{C}}}_{X} as $n\to\infty$. Hence, for large enough $n$,

[TABLE]

and so,

[TABLE]

and

[TABLE]

So, we have shown that,

[TABLE]

Using the similar trick as used in the proof of Lemma A.2, i.e. dividing the region into unit blocks and also using the fact that \lvert\left(\mbox{{\mathbf{A}}}_{n}\right)_{ij}\rvert=O(n^{-1}), we can show that

[TABLE]

which concludes the proof. ∎

A.2.6 Proof of Theorem 3.3

Proof.

Recall, \mathbf{l}_{n}^{\text{ols}}=\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}}\right)^{-1}\left(n^{-1}\mbox{{\mathbf{X}}}^{\prime}\mathbf{w}\right). So,

[TABLE]

where \mbox{{\mathbf{H}}}_{n}=n^{-1}\mathbf{X}(n^{-1}\mathbf{X}^{\prime}\mathbf{X})^{-2}\mathbf{X}^{\prime}, $\bm{\Sigma}_{w}=\mathrm{Var}\left(\mathbf{w}\right)$. Let us also denote n^{-1}\mbox{{\mathbf{X}}}^{\prime}\mbox{{\mathbf{X}}} by \mbox{{\mathbf{C}}}_{n} and by condition (C.10), \mbox{{\mathbf{C}}}_{n}\to\mbox{{\mathbf{C}}}_{X} as $n\to\infty$. Note that \text{tr}\left(\mbox{{\mathbf{H}}}_{n}\right)=\text{tr}\left(\mbox{{\mathbf{C}}}_{n}^{-1}\right)\to\text{tr}\left(\mbox{{\mathbf{C}}}_{X}^{-1}\right)=2C^{(0)}>0. So, for large enough $n$, \text{tr}\left(\mbox{{\mathbf{H}}}_{n}\right)\geq C^{(0)}>0. Also, for large enough $n$,

[TABLE]

For any $\varepsilon>0$ let us define,

[TABLE]

and let us denote the cardinality of this set by $m_{n}$ (dependent on $\varepsilon$). Now consider a $\varepsilon>0$ such that $\varepsilon<\text{min}\left(C^{(0)},C^{(1)}\right)$. Then,

[TABLE]

and by simple algebric calculation we have $m_{n}\geq\left(C^{(0)}-\varepsilon\right)\left(C^{(1)}-\varepsilon\right)^{-1}n$. Recall, E(w(\mbox{{\mathbf{s}}}_{i})^{2}) is equal to $\sigma_{\epsilon}^{2}$ if $i\in G_{n}$ and is equal to $\sigma_{\epsilon}^{2}+\tau_{i}^{2}$ otherwise. So,

[TABLE]

Finally, using the similar methodology as in proof of Theorem 3.2 we can show that \lvert\frac{1}{n}\sum_{i\neq j}(\mbox{{\mathbf{H}}}_{n})_{ij}E\left(w(\mbox{{\mathbf{s}}}_{i})w(\mbox{{\mathbf{s}}}_{j})\right)\rvert=O(\lambda_{n}^{-4}). Hence we have shown that, for any $0<\varepsilon<\text{min}\left(C^{0},C^{1}\right)$ and large enough $n$,

[TABLE]

∎

A.2.7 Proof of Corollary 3.2

Proof.

Observe that, for large enough $n$, for all $i\in\left\{1,\dots,n\right\}$,

[TABLE]

Clearly, from the assumption of the Corollary 3.2, $C_{2}>0$. Take any $\varepsilon>0$ such that $\varepsilon<C_{2}$. Then, for all $i\in\left\{1,\dots,n\right\}$, \left(\mbox{{\mathbf{H}}}_{n}\right)_{ii}>n^{-1}\varepsilon and hence,

[TABLE]

The rest of the proof is similar to that of Theorem 3.3. ∎