# Parameter uncertainty estimation for exponential semi-variogram models: Two
generalized bootstrap methods with check- and quantile-based filtering

Julia Dycka and Odile Sauzeta,b CONTACT Julia Dyck. Email: j.dyck@uni-
bielefeld.de aDepartment of Business Administration and Economics, Bielefeld
University, Bielefeld, Germany; bSchool of Public Health, Bielefeld
University, Bielefeld, Germany.

###### Abstract

The estimation of parameter standard errors for semi-variogram models is
challenging, given the two-step process required to fit a parametric model to
spatially correlated data.

Motivated by an application in the social-epidemiology, we focus on
exponential semi-variogram models fitted to data between 500 to 2000
observations and little control over the sampling design. Previously proposed
methods for the estimation of standard errors cannot be applied in this
context. Approximate closed form solutions are too costly using generalized
least squares in terms of memory capacities. The generalized bootstrap
proposed by Olea and Pardo-Igúzquiza is nonetheless applicable with weighted
instead of generalized least squares. However, the standard error estimates
are hugely biased and imprecise.

Therefore, we propose a filtering method added to the generalized bootstrap.
The new development is presented and evaluated with a simulation study which
shows that the generalized bootstrap with check-based filtering leads to
massively improved results compared to the quantile-based filter method and
previously developed approaches. We provide a case study using birthweight
data.

###### keywords:

Exponential semi-variogram model; generalized bootstrap; standard error

GLS: generalized least squares; LS: least squares; WLS: weighted least squares

##  1 Introduction

###  1.1 Contextual backround

Semi-variogram modelling is a method used in spatial statistics to model the
relationship between pairwise correlation of observations and the distance
between these observations. It is commonly used in fields like agriculture or
ecology but has also seen applications in social epidemiology to study
neighbourhood effects on health. In the latter, the estimation of exponential
parametric models is of particular interest because it provides a measure of
which part of the total variance is spatially structured and which local
spatial scale is relevant. However the assessment of parameter uncertainty
presents some difficulties while being necessary for proper inference [1].

Contrary to some other fields of application, data in social epidemiology
cannot be sampled according to a scheme designed to facilitate the estimation
of parameters. On the contrary, analyses are usually based on survey data in
which there are few observations at very small distances [2] \- providing very
few data for the estimation of the first points of the empirical semi-
variogram - , but with possibly large sample sizes of over 1000 - making, for
instance, inverting matrices for the estimation of standard errors difficult.
Exemplary applications of semi-variogram modelling can be found in [2] and
[3].

###  1.2 Technical challenges

Three characteristics of semi-variogram estimation are a source of technical
difficulties for the construction of parameter estimators and, thereby, for
the construction of parameter uncertainty estimators. First, the estimation of
parameters associated with correlated data has to take the correlation
structure into account. Consequently, the simplifying assumption of an
underlying data generating process with independent and identically
distributed elements is not appropriate [4, p.3]. Second, the estimation of
the semi-variogram is not based on the spatial observations themselves but on
their pairwise Eucledian distances [5]. Third, the fitting of a semi-variogram
model requires a two-step process that includes, first, the estimation of a
discrete semi-variogram and, second, based on that the fitting of a mostly
continuous model. Multiple steps of transformation and estimation intuitively
lead to multiple potential sources of uncertainty that influence the precision
of the semi-variogram model parameters.

###  1.3 Existing approaches

Methods for the parameter standard error estimation have been proposed in the
past.

An analytic solution was developed by [6] for Gaussian data generating
processes and parameter estimates obtained using the generalized least squares
(GLS) estimation approach. While it leads to reliable standard error estimates
for datasets with up to 100 observations, results are often not obtained for
larger sample sizes of 1000 or higher due to the at least quadratically
increasing amount of memory capacities required for the computation. Moreover,
the theoretical statistical assumptions made in order to construct the
analytic solution are often not fulfilled, such that solutions, if feasible,
have to be interpreted to be approximate rather than exact.

An empirical approach using the so called generalized bootstrap is provided by
[7]. The reduced amount of theoretical assumptions and the use of the rather
simple bootstrap structure which can be implemented quite efficiently is
advantageous. However, as for the previous method, the generalized bootstrap
is constructed with focus on data with less than 100 observations, high
control over the sampling design and locations, and parameter estimates
obtained with the GLS approach as is typically the case in the geological
context. Larger datasets and the often conducted parameter estimation based on
bootstrap samples forced us to substitute the GLS parameter estimation by the
WLS estimation procedure in order to make it computationally feasible. Given
these changes, the standard error estimates are highly biased upwards and very
imprecise.

A more recent attempt uses an approximate maximum likelihood estimation
algorithm with sparse distance matrices for the parameter estimation to fasten
the computation procedure of the standard error estimates via generalized
bootstrap [8]. The approach focuses on applications to very large datasets
with 10000 or more data points and assumes approximate Gaussianity for the
maximum likelihood estimation. Datasets with a sample size of about 1000
appear to be too small to account for the errors due to sparse computation, so
that the semi-variogram parameter estimators themselves are hugely biased and
therefore not an appropriate base for uncertainty estimation.

###  1.4 Goal

Estimating standard errors based on one of the existing methods does not
provide satisfactory results in the data scenarios of our interest. Focusing
on typical social-epidemiological data to which exponential semi-variogram
models are to be fitted, we need a standard error estimation tool that

  1. 1.

works under as few distributional assumptions as possible,

  2. 2.

is based on the commonly used and computationally feasible WLS model fitting
approach,

  3. 3.

is applicable to data with sample sizes of about 500 to 2000 and little or no
control over the sampling design.

Given the limitation of existing methods in this context we propose a
generalized bootstrap with filtering based on the original generalized
bootstrap by Olea and Pardo-Igúzquiza [7] using WLS rather than GLS.

To compensate for the higher variance of the standard error estimator when
using WLS for the parameter estimation, the generalized bootstrap procedure is
extended by a filter. The filter evaluates each bootstrap estimate and
excludes results of a bootstrap estimation if the semi-variogram model fitting
algorithm with WLS is suspected to have not converged. Two filter options are
presented and evaluated: the check-based filtering and the quantile-based
filtering.

##  2 Methods

###  2.1 Statistical model

We assume an isotropic and second-order stationary spatial process
$\\{Z(\mathbf{s})\mid\mathbf{s}\in D\\}$ with a collection of random variables
$Z(\mathbf{s})$ and a set of locations $D$. In this case, the semi-variogram
of a spatial process characterizes the semi-variance structure of the process
depending on the Euclidean distance $h=||\mathbf{s}-\mathbf{s}^{\prime}||$
between two subjects [9, p.135–140]. The semi-variogram is defined as

$\displaystyle\gamma(h)=\frac{1}{2}Var(Z(\mathbf{s})-Z(\mathbf{s}^{\prime}))$
(1)

for two subjects at locations $\mathbf{s}$ and $\mathbf{s}^{\prime}$ with
distance $h:=||\mathbf{s}-\mathbf{s}^{\prime}||$.

Under the stated conditions the semi-variogram and the covariance function of
the spatial process have the relationship:

$\displaystyle C(h)=c_{0}+\sigma_{0}^{2}-\gamma(h)$ (2)

where $c_{0}+\sigma_{0}^{2}=Var[Z(\mathbf{s})]$ is the variance of the outcome
to be analysed. It is assumed to be constant over the studied surface. The
summand $c_{0}$ is called nugget effect and defines the value of $\gamma(h)$
for $h$ approaching zero. The so called partial sill $\sigma_{0}^{2}$
describes the asymptotic covariance of two observations on the same location.

Given data, Matheron’s estimator provides an unbiased empirical estimate of
the semi-variogram [10]. It is defined as

$\displaystyle\hat{\gamma}(h)=\frac{1}{2\cdot|N(h)|}\sum_{(\mathbf{s_{i}},\mathbf{s_{j}})\in
N(h)}\\{Z(\mathbf{s_{i}})-Z(\mathbf{s_{j}})\\}^{2},$ (3)

where $N(h):=\\{(\mathbf{s_{i}},\mathbf{s_{j}})\in
D|\;\;||\mathbf{s_{i}}-\mathbf{s_{j}}||=h\\}$. To guarantee a sufficient
cardinality of $N(h)$ for each realized $h$, multiple lag distances are
grouped into lag intervals of the form $[h-\varepsilon,h+\varepsilon]$ for a
discrete number of lags $h_{1},\dots,h_{K}$.

Whereas the empirical semi-variogram provides point estimates for a finite
number of lag intervals, an almost continuous model is obtained by a
parametric semi-variogram model. The exponential semi-variogram model is
defined as

$\displaystyle\gamma_{\text{exp}}(h;c_{0},\sigma^{2}_{0},\phi)=c_{0}+\sigma_{0}^{2}\Big{\\{}1-\exp\big{(}-\frac{h}{\phi}\big{)}\Big{\\}}\;\;$
$\displaystyle\text{for}\;\;h>0$ (4)

and $\gamma_{\text{exp}}(h;c_{0},\sigma^{2}_{0},\phi)=0$ for $h=0$ with nugget
effect $c_{0}>0$, partial sill $\sigma_{0}^{2}>0$ and range parameter $\phi>0$
[4, p.61].

###  2.2 Model fitting

The semi-variogram model fitting is commonly done via the method of least
squares (LS), WLS or GLS [9, p.163–166]. In all three cases the parameter
vector $\theta=(c_{0},\sigma_{0}^{2},\phi)^{t}$ that minimizes the loss
function

$\displaystyle
LOSS(\theta)=[\hat{\mathbf{\gamma}}-\mathbf{\gamma_{\text{exp}}(\theta)}]^{t}\;\mathbf{\Gamma}^{-1}\;[\hat{\mathbf{\gamma}}-\mathbf{\gamma_{\text{exp}}(\theta)}]$
(5)

is chosen. The vector
$\hat{\mathbf{\gamma}}=(\hat{\gamma(d_{1})},\dots,\hat{\gamma(d_{K})})^{t}$
contains the empirical semi-variogram evaluated at $K$ lag distances and
$\mathbf{\gamma_{\text{exp}}(\theta)}=(\gamma_{\text{exp}}(d_{1};\mathbf{\theta}),\dots,\gamma_{\text{exp}}(d_{K};\mathbf{\theta}))^{t}$
is the corresponding semi-variogram model evaluated at the same lag distances.
In practice, an estimate for the covariance matrix of the empirical semi-
variogram denoted by $\mathbf{\Gamma}\in\mathbb{R}^{K\times K}$ is needed. In
case of the LS method, $\mathbf{\Gamma}=\sigma^{2}\mathbf{I}$ is assumed
meaning that the variance for each lag’s estimate is constant the lags’ semi-
variances are uncorrelated.

In case of the WLS method, $\mathbf{\Gamma}$ is a diagonal matrix with varying
entries and thus assumes varying variances for each lag’s estimate. Still, a
covariance of zero between two lags is preset.

In case of the GLS method, varying variances of the lags’ semi-variogram
estimates are allowed, as well as covariances between two lags. This is
reflected by a $\mathbf{\Gamma}$ matrix with potentially non zero entries on
the diagonal and off-diagonal [11]. Hence, the latter can be expected to lead
to the best results. However, the computational effort regarding the memory
capacities increases at minimum quadratically with the sample size $N$ since
the number of pairwise distances that go into the estimation of
$\mathbf{\Gamma}$ increases proportional to $N^{2}$ making calculations for
large N infeasible.

Therefore, the WLS method will be used as next best and applicable solution
for the semi-variogram parameter estimation. We choose weights depending on
the number of data pairs that are considered for the estimation of the
empirical semi-variogram according to Matheron’s estimator. The covariance
matrix estimate has entries $\mathbf{\Gamma}_{kk}=1/|N(d_{k})|$ on the
diagonal and $\mathbf{\Gamma}_{kl}=0$ on the non-diagonal for $k,l=1,\dots,K$
and $k\not=l$. Thereby the minimizing the difference between empirical semi-
variogram and model at lag distances with more information is prioritized. The
general formula of the loss function to be minimized then reduces to

$\displaystyle
LOSS_{\text{WLS}}(\theta)=\sum_{k=1}^{K}|N(d_{k})|\;\cdot[\hat{\gamma(d_{k})}-\gamma_{\text{exp}}(d_{k};\theta)]^{2}$
(6)

which is implemented as default option in the `variofit` function from within
the Rpackage geoR [12].

###  2.3 Generalized bootstrap with filtering

The generalized bootstrap approach [7] serves as starting point, but requires
modification to be applicable in the stated data framework of our interest.
Each semi-variogram model is fitted using the WLS instead of the GLS method.
Moreover, substantial changes are made with respect to the purpose of the
bootstrap estimates. Whereas the original algorithm’s goal is the estimation
of the covariance matrix and confidence bands of the empirical semi-variogram
for each defined lag interval, the focus here lies on the estimation of model
parameter standard errors of a fitted exponential semi-variogram model. As a
completely new element, a filter within the bootstrapping process tackles the
exclusion of improbable bootstrap estimations due to non-convergence of the
semi-variogram model fitting algorithm.

Considering these aspects, the following algorithm was developed:

  1. 0.

Exponential semi-variogram model: An empirical semi-variogram is estimated
based on the original spatial dataset $\mathbf{z}=(z_{1},\dots,z_{N})^{t}$
plus the spatial coordinates $\mathbf{s}$ with
$\mathbf{s_{i}}\in\mathbb{R}^{2}$ for $i=1,\dots,N$. Subsequently, an
exponential semi-variogram model $\gamma_{\text{exp}}$ is fitted using the
weighted least squares method producing the parameter estimate
$\mathbf{\hat{\theta}}$ for the true parameter vector $\mathbf{\theta}$ with
$(\theta_{1},\theta_{2},\theta_{3})=(c_{0},\sigma^{2}_{0},\phi)$.

  2. 1.

Normal score transformation: The data vector $\mathbf{z}$ is mapped into a
Gaussian space by the empirical normal score transformation function $\varphi$
that maps the data from the attribute space into a Gaussian space.
Consequently, $\mathbf{y}=\varphi(\mathbf{z})$ is a realization vector of a
standard normal random variable [13, 14].

  3. 2.

Exponential semi-variogram model for transformed data: An empirical semi-
variogram is fitted to $\mathbf{y}$ combined with the original data
coordinates. Based on that, an exponential semi-variogram model
$\tilde{\gamma}_{\text{exp}}$ is fitted using the WLS method providing the
parameter estimate $\mathbf{\tilde{\theta}}$.

  4. 3.

Covariance estimation: Making use of the semivarigram model
$\tilde{\gamma}_{\text{exp}}$, the covariance between two data points $z_{i}$
and $z_{j}$ is calculated based on their Euclidean distance
$d_{ij}=d(s_{i},s_{j})$ as

$c_{ij}=\tilde{c_{0}}+\tilde{\sigma}_{0}^{2}-\tilde{\gamma}_{\text{exp}}(d_{ij}).$
(7)

The $N\times N$ covariance matrix with entries $c_{ij}$ is denoted as
$\mathbf{C}$.

  5. 4.

Decorrelation of the data: The decomposition of the covariance matrix
$\mathbf{C}$ into the product of the lower triangular matrix $\mathbf{L}$ and
its transposed is obtained by the Cholsky decomposition algorithm as

$\mathbf{C}=\mathbf{L}\mathbf{L}^{t}.$ (8)

This is possible because $\mathbf{C}$ is square symmetric and positive-
definite [15]. The vector

$\mathbf{x}=\mathbf{L}^{-1}\mathbf{y}$ (9)

is a set of independent and identically distributed, hence uncorrelated
values.

  6. 5.

Classical bootstrap: Sampling $N$ times with replacement from $\mathbf{x}$
generates a bootstrap sample $\mathbf{x^{*}}$.

  7. 6.

Recorrelation: The bootstrap sample $\mathbf{x^{*}}$ reinherits the
correlation structure by applying the inverse operation of step 4, that is

$\mathbf{y^{*}}=\mathbf{Lx^{*}}.$ (10)

  8. 7.

Normal score back transformation: The back transformation of $\mathbf{y^{*}}$
to the attribute space through the inverse normal score function in step 1 is
done by

$\mathbf{z^{*}}=\varphi^{-1}(\mathbf{y^{*}}).$ (11)

  9. 8.

Analysis of the bootstrap sample: An exponential semi-variogram model is
estimated based on $\mathbf{z^{*}}$ combined with the original coordinates
providing an estimate $\mathbf{\theta^{*}}$.

From here on, the algorithm depends on the choice of filter method. If the
check-based filter is used, steps 9c, 10c and 11 are executed.

  1. (9c)

Check-based filtering: A test is applied to the bootstrap estimate
$\theta^{*}$ to indicate whether the exponential semi-variogram fitting
algorithm did converge. If the algorithm is assumed to have failed, the
bootstrap estimate is discarded.

  2. (10c)

Repetition: The steps 5 to 9 are repeated until a set of $B$ bootstrap
estimates ${\\{\theta^{*}_{b}\\}}_{b=1,\dots,B}$ has aggregated.

If the quantile-based filter is used, steps 9q, 10q and 11 are conducted.

  1. (9q)

Repetition: Given the quantile threshold $\alpha\in(0,1)$, the steps 5 to 8
are repeated $\tilde{B}=\frac{1}{\alpha}\cdot B$ times providing a set of
$\theta^{*}_{b}$ for $b=1,\dots,\tilde{B}$.

  2. (10q)

Quantile-based filtering: The first $B$ values of the increasingly ordered set
of bootstrap estimates $\theta^{*}_{b}$ is saved. This equals the set of
bootstrap estimates between the minimum value and the $\alpha$-quantile.

  3. (11)

Parameter standard error estimation: Based on the filtered set of repeatedly
estimated parameters $\mathbf{\theta}^{*}_{b}$ for $b=1,\dots B$ the parameter
standard error estimates are defined as the empirical standard deviations

$\displaystyle\widehat{se(\theta_{j})}=sd(\theta^{*}_{j})=\sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\Big{\\{}\theta_{bj}^{*}-\overline{\theta_{j}^{*}}\Big{\\}}^{2}},$
(12)

for $j=1,...,3$ referring to the three parameters $c_{0},\sigma^{2}_{0},\phi$.

####  2.3.1 Check-based filtering

The check-based filter (steps 9c, 10c) is applied in each repetition. The
bootstrap estimate $\theta^{*}$ is used to compute the corresponding overall
variance according to the model as $c_{0}^{*}+\sigma_{0}^{2*}$. If

$\displaystyle c_{0}^{*}+\sigma_{0}^{2*}>\tau\widehat{Var(Z)},$ (13)

i.e. if the variance indicated by the model exceeds the sample variance times
some threshold factor $\tau$, the bootstrap estimate $\theta^{*}$ is
discarded. Otherwise, it is retained. The size of the bootstrap sample $B$ is
specified beforehand and the procedure of estimating and checking one
bootstrap estimate is repeated until $B$ bootstrap estimates are accepted.

####  2.3.2 Quantile-based filtering

The quantile-based filter (steps 9q, 10q) is applied after the computation of
all bootstrap estimates. Given a threshold percentile $\alpha$ all bootstrap
estimates $\theta^{*}_{b}$ beyond the $\alpha\cdot 100\%$-quantile
$q_{\alpha}(\theta^{*}_{j})$ for $j=1,2,3$ are discarded while all bootstrap
estimates up to the $\alpha\cdot 100\%$-quantile are saved.

To make sure the set of bootstrap estimates used for the computation of the
parameter standard error contains $B$ elements, the number of bootstrap
samples to be generated in the first place is set to
$\tilde{B}=B\cdot\frac{1}{\alpha}$.

##  3 Simulation Study

The simulation study is conducted in order to test and evaluate the
performance of the presented filter methods in various sampling frameworks and
provide recommendations for meta- and tuning parameter specification. The
implementation is constructed using the statistical software environment R
[16]. The R package RandomFields [17] is used to simulate spatial data. The
packages geoR [12], SpatialTools [18] and additional individual functions [14]
are used for the implementation of the generalized bootstrap algorithm with
each filter method.

###  3.1 Data generating process

Data is generated based on a zero-mean Gaussian random process with an
exponential semi-variogram model of the form

$\displaystyle\gamma(h)=60+40\cdot\Big{\\{}1-\exp\Big{(}-\frac{h}{200}\Big{)}\Big{\\}}\;\;$
$\displaystyle\text{for}\;\;h>0.$ (14)

The population is located on a two-dimensional square grid
$D=\\{0,\dots,10000\\}^{2}$. One complete realization of locations is
generated by drawing 250000 data points using the `runif` command. The
attribute values $z$ at each realized location are computed with the
`RFsimulate` command.

###  3.2 Sampling schemes

Samples are produced by sampling from the simulated population without
replacement using the `sample` function. Sample scenarios vary with respect to
the sample size $N$ ($500,1000,2000$) and the sampling density (low, middle,
high). In case of a low density, the sampled locations are uniformly
distributed over the whole grid surface. Consequently, many point pairs with
similarly valued pairwise distances occur in the middle of the pairwise
distances range. In comparison, a high sampling density provides more small
valued pairwise distances (Figure 1).

Figure 1: Example for sampling densities and resulting distribution of
pairwise distances. Three subsets are sampled from the simulated population
with $N=1000$ and either a low, middle or high sampling density. The upper row
shows the sample locations on the map $D$ according to the sampling density
schemes. The empirical distribution of the pairwise distances is provided
below.

###  3.3 Metaparameters for the empirical semi-variogram

The maximal distance specifies the maximum value up to which pairwise
distances are considered within the semi-variogram estimation process. Three
values $1.1\cdot h^{*},1.5\cdot h^{*}$ and $2\cdot h^{*}$ are used. The true
practical range according to the given semi-variogram model is $h^{*}=599.15$.

The number of lag intervals is chosen to be $10$. According to that, the range
of pairwise distances from 0 to the maximal distance is divided into 10
equidistant lag intervals. For each interval, a point estimate is calculated
based on all point pairs within the lag bin for the empirical semi-variogram.

###  3.4 Tuning parameters for the filter methods

For the check-based filter approach, the factor $\tau$ scaling the threshold
of the test has to be specified. Within the frame of the simulation study, the
values $1.1,1.5,2.0,2.5,3.0$ are applied.

In regards to the quantile-based filter method, the percentage of the quantile
threshold $\alpha$ has to be chosen. Values ranging from $0.75$ to $1.0$ with
a steplength of $0.05$ are put to the test. The value 1.0 refers to the case
with no filtering.

###  3.5 Conduction of the simulation study and elimination of improbable runs

Each scenario is simulated $n_{sim}=3000$ times. Very large parameter
estimates for at least one of the estimated semi-variogram parameters are an
indicator that the algorithm did not converge in the estimation process. Only
estimated models fulfilling

$\displaystyle\hat{c_{0}}<1000\;\;\text{and}\;\;\hat{\sigma_{0}^{2}}<1000\;\;\text{and}\;\;\hat{\phi}<1000$
(15)

are approved and included into the evaluation of the results. This leaves
$\tilde{n}_{sim}<n_{sim}$ simulation runs to be analysed. The proportion
$\tilde{n}_{sim}/n_{sim}$ will be referred to as the convergence rate of the
parameter estimates in the next section.

###  3.6 Performance measures

The performance statistics are presented as averages grouped by the sample
size, the sampling density, the maximal distance factor and the tuning
parameter specification for each filter method. The empirical standard
deviation over the parameter estimates of all simulation runs
$\eta_{\theta_{j}}=sd(\hat{\theta_{j}})$ serves as approximation to the true
standard error. Its variability is obtained by the Monte Carlo standard error
estimator
$se_{MC}(\eta_{\theta_{j}})=sd(\hat{\theta_{j}})/\sqrt{2(\tilde{n}_{sim}-1)}$.
The performance will be evaluated in terms of the estimates’ standard
deviation, bias and mean squared error.

The grouped average estimate $\widehat{\eta_{\theta_{j}}}$ according to each
grouping parameter is provided along with its grouped average of the
corresponding standard deviations $sd(\widehat{\eta_{\theta_{j}}})$. For
instance, the grouped average estimate given the grouping factor $N=500$ is
obtained by averaging over all estimates for sampling scenarios with $N=500$
and arbitrary realizations for the sampling density, maximal distance factor
and tuning parameter. The standard deviations of the same scenarios are
averaged to obtain the grouped average standard deviation.

The bias is defined as the difference

$\displaystyle
bias(\widehat{\eta_{\theta_{j}}})=\widehat{\eta_{\theta_{j}}}-\eta_{\theta_{j}}.$
(16)

The overall performance is measured in terms of the mean squared error (MSE)

$\displaystyle
MSE(\widehat{\eta_{\theta_{j}}})=sd(\widehat{\eta_{\theta_{j}}})^{2}+bias(\widehat{\eta_{\theta_{j}}})^{2}$
(17)

taking into account the performance in regards to both, the location of the
estimated standard error and the variability described by the squared
empirical standard deviation of the estimated standard error estimates.

##  4 Results

In the following, the simulation study results are described. Compact result
tables can be found in the appendix.

###  4.1 Convergence rate for the estimation of the parameters

The convergence rate $\tilde{n}_{sim}/n_{sim}$ is lowest for a sample size of
500 with $87.5\%$. It increases with increasing sample size leading to
$97.0\%$ for $N=1000$ and $99.9\%$ for $N=2000$. An increase in the sampling
density leads to a slight increase in the convergence rate ranging from
$92.9\%$ for a low sampling density to $95.6\%$ in case of middle and $95.9\%$
in case of a high sampling density. Similarly, an increase in the maximal
distance factor provides a slight increase in the convergence rate from
$91.5\%$ for a maximal distance factor of $1.1$ over $96.1\%$ for a maximal
distance factor of $1.5$ to $96.8\%$ for a maximal distance factor of $2$
(Table A.1).

###  4.2 Performance Evaluation: Check Filter Method

Because the performance of the standard error estimates for the nugget
parameter and partial sill are in a similar range, the effects of the sampling
framework on them will be presented together first. The standard error
estimates of the shape parameter differ in that their performance statistics
are higher. They will be summarized afterwards. Tables containing the results
of the simulation study results can be found in the appendix (Tables A.2, A.3,
A.4).

####  4.2.1 Effects on standard error estimates of nugget effect and partial
sill

##### Effect of sample size

The standard deviation of $\widehat{\eta_{c_{0}}}$ varies between 0.88 and
5.14 and decreases for increasing sample size. The bias is positively valued
in all cases and varies between 1.006 for a sample size of 1000 and 0.14 for a
sample size of 2000. A sample size of 500 leads to a value in between. The MSE
as overall performance measure varies between 0.80 and 28.83 and decreases
with increasing sample size.

The standard deviations of $\widehat{\eta_{\sigma^{2}_{0}}}$ are located
between 1.04 and 4.87. The bias varies between 0.52 and 1.83 and the MSE
between 1.36 and 26.22 decreasing for larger sample size.

##### Effect of sampling density

The standard deviations of $\widehat{\eta_{c_{0}}}$ are located between 2.70
and 3.84. The middle sampling density provides the lowest standard deviation
of 2.70, closely followed by the value of 2.81 when given a high sampling
density. The bias varies between the minimum of 0.42 for a high sampling
density and the maximum of 0.77 for a middle sampling density. The MSE ranges
between 7.86 for a middle sampling density, closely followed by 8.09 for a
high sampling density, and 15.14 for a low sampling density.

The standard deviation of $\widehat{\eta_{\sigma^{2}_{0}}}$ takes values
between 2.70 and 3.78 decreasing with increasing sampling density. The bias
varies between 1.11 given a high sampling density and 1.60 given a middle
sampling density. The MSE lies within the range of 8.54 and 15.82 decreasing
with increasing sampling density.

##### Effect of maximal distance

The standard deviation of $\widehat{\eta_{c_{0}}}$ varies between 3.01 and
3.26. The bias ranges from 0.15 to 0.97. Together they lead to MSE values
between 9.07 and 11.54 decreasing for larger maximal distance factor.

The standard deviation of $\widehat{\eta_{\sigma^{2}_{0}}}$ lies within the
range of 2.87 and 3.40. The bias ranges from 0.44 to 2.22. The MSE is located
between 8.45 and 16.50 and decreases when the maximal distance becomes larger.

##### Effect of tuning parameter

Observed standard deviation, bias and MSE of $\widehat{\eta_{c_{0}}}$ are
almost equal for all choices of $\tau$. For example, the MSE ranges between
10.03 and 10.21.

The standard deviation and bias of $\widehat{\eta_{\sigma^{2}_{0}}}$ show
slight variation depending on $\tau$. The MSE as overall performance measure
lies between 10.10 and 13.58. It is smallest for $\tau=1.5$, closely followed
by $\tau=1.1$ and becomes larger for $\tau>1.5$.

##### Summary of effects

The impact of varying the sampling scenarios and tuning parameter
specification on the performance of the nugget effect and partial sill can be
summarized as follows. Increasing the sample size improves the overall
performance. A sampling design with a middle to high density instead of a low
density induces an improvement. Increasing the maximal distance leads to
slightly better performance for $\widehat{\eta_{c_{0}}}$ and more clearly for
$\widehat{\eta_{\sigma^{2}_{0}}}$. The choice for the tuning parameter $\tau$
has no clear effect on the performance of $\widehat{\eta_{c_{0}}}$. The
performance of $\widehat{\eta_{\sigma^{2}_{0}}}$ is slightly improved when
smaller $\tau=1.1,1.5$ are used.

####  4.2.2 Effects on standard error estimates of shape parameter

While the overall performance measure MSE reaches at maximum a value of 26.83
for the nugget effect and partial sill standard errors, the MSE results for
$\widehat{\eta_{\phi}}$ varies within $1.95\text{\times}{10}^{3}$ at best and
$1.88\text{\times}{10}^{5}$ at worst regarding all grouped scenario averages.
Because of that, we confine ourselves to a summary of the effects.

##### Summary of effects

Regarding the sampling conditions, clear effect trends can be observed.
According to all three performance statistics, relatively best performance is
obtained, when the highest sample size of $2000$, the middle sampling density,
the low maximal distance factor of $1.1$ and the minimal tuning parameter
$\tau=1.1$ are given.

###  4.3 Performance evaluation: Quantile filter method

While the performance measures of the standard error estimates for the nugget
parameter lie in a similar range as in case of the check-filter method, the
performance of the standard errors obtained for the partial sill and shape
parameter are very high valued. Therefore, only a short summary of the
observed effects is given for each parameter’s standard error estimation.
Tables containing the simulation study results can be found in the appendix
(Tables A.5, A.6, A.7).

####  4.3.1 Effects on standard error estimates of nugget parameter

The performance statistics for the standard error estimates
$\widehat{\eta_{c_{0}}}$ lie in the following value ranges. Standard deviation
varies between $0.77$ and $4.68$, bias between $-1.54$ and $0.68$, MSE between
$1.01$ and $23.51$ across all calculated grouped average statistics.

##### Summary of effects

The impact of varying the sampling scenarios and tuning parameter
specification can be summarized as follows. Increasing the sample size
provides a clear improvement of the overall performance according to the MSE.
A sampling design with a middle to high density instead of a low density
provides an improvement. Changing the maximal distance and tuning parameter
$\alpha$ have almost no impact on the performance.

####  4.3.2 Effects on standard error estimates of partial sill

The observed statistics for the standard error estimates
$\widehat{\eta_{\sigma_{0}^{2}}}$ are very high valued. The MSE describing the
overall performance ranges between $1.09\text{\times}{10}^{5}$ and
$2.19\text{\times}{10}^{8}$ considering all grouped sample scenarios.

##### Summary of effects

Regarding the sampling conditions, clear effect trends can be observed.
According to all three performance statistics, relatively best performance is
obtained when the highest sample size of $2000$, the high sampling density,
the highest maximal distance factor of $2.0$ and the minimal tuning parameter
$\alpha=0.75$ are given.

####  4.3.3 Effects on standard error estimates of shape parameter

The standard error estimates $\widehat{\eta_{\phi}}$ estimated within the
simulation study lead to very high numbers in the performance statistics. The
MSE takes values from $1.15\text{\times}{10}^{8}$ and
$2.81\text{\times}{10}^{11}$ over all grouped sample scenarios.

##### Summary of effects

Effect trends are observable for each sampling design and tuning parameter.
All three performance statistics indicate relatively best performance when the
highest sample size of $2000$, the middle or high sampling density, a maximal
distance factor of $1.5$ and the minimal tuning parameter $\alpha=0.75$ are
used.

##  5 Discussion

Two extensions of the the generalized bootstrap algorithm were implemented:
the generalized bootstrap with check-based filtering and the generalized
bootstrap with quantile-based filtering. In the following, we discuss the
observed effects on the estimator performance and the drivers behind them.
Finally, we provide recommendations for application in praxis.

###  5.1 Check-based filtering

The standard error estimators $\widehat{\eta_{c_{0}}}$ and
$\widehat{\eta_{\sigma^{2}_{0}}}$ show similar, good performance. The
estimator $\widehat{\eta_{\phi}}$ provides larger bias and standard deviation
values. The differing performances of $\widehat{\eta_{c_{0}}}$ and
$\widehat{\eta_{\sigma^{2}_{0}}}$ on the one hand and $\widehat{\eta_{\phi}}$
on the other hand are likely driven by the form of distribution of the
parameter estimators $\hat{c_{0}},\hat{\sigma^{2}_{0}}$ and $\hat{\phi}$
(Figure 2). In case of the nugget effect and partial sill, the empirical
distributions of the parameter estimates appear approximately symmetric with
$sd(\hat{c_{0}})=14.68$ and $sd(\hat{\sigma^{2}_{0}})=14.55$. The distribution
of the shape parameter is more skewed with some extremely large estimates and
an $sd(\hat{\phi})=115.24$. A fitted normal distribution assigns a significant
amount of probability to values below zero for $\hat{\phi}$. The shape
parameter’s standard error estimates lead to an empirical distribution that
assigns much probability to values below zero, too. Moreover, it is heavy
tailed on the right side. The skewness and variance in the empirical
distribution is suspected to be the reason for the higher variability and
overestimation in case of $\widehat{\eta_{\phi}}$.

Figure 2: Histograms of parameter and standard error estimates. Histograms of
the model parameter estimates $\hat{c}_{0},\hat{\sigma}^{2}_{0}$, and
$\hat{\phi}$ (upper row) and the standard error estimates
$\widehat{\eta_{c_{0}}},\widehat{\eta_{\sigma^{2}_{0}}}$, and
$\widehat{\eta_{\phi}}$ (lower row) are presented each from left to right.
Each histogram is based on the sum of simulation results obtained in all
simulation scenarios ($\tilde{N}_{sim}=76789$). The presented standard error
estimates are calculated with the generalized bootstrap method with check-
filtering and $\tau=1.5$. In addition, the mean, median and normal density
with empirical mean and standard deviation of the data are marked to visualize
the degree of symmetry.

The overall performance of the check-based standard error estimators is
improved under specific framework conditions. A larger sample size increases
the number of data pairs that are incorporated into the semi-variogram
estimation procedure. A higher sampling density increases the amount of
shortly distanced data pairs that lie under the specified maximal distance
threshold and thereby increase the number of data pairs incorporated into the
semi-variogram estimation procedure. Especially, the performance of the
variance parameter uncertainties $\widehat{\eta_{c_{0}}}$ and
$\widehat{\eta_{\sigma^{2}_{0}}}$ become better when more data pairs are at
hand. Moreover, for $\widehat{\eta_{c_{0}}}$ and
$\widehat{\eta_{\sigma^{2}_{0}}}$ the performance of the estimator is improved
with increasing maximal distance factor as it further increases the amount of
data incorporated into the semi-variogram estimation procedure. In contrast, a
maximal distance specification that is close to the true maximal distance -
within the simulation study the lowest factor of $1.1$ \- optimizes the
performance of the shape parameter’s standard error estimator
$\widehat{\eta_{\phi}}$. Larger maximal distance factors have a huge negative
impact on the $\widehat{\eta_{\phi}}$ estimates. Regarding the tuning
parameter $\tau$, best estimation performances are obtained for $\tau=1.1,1.5$
according to the simulation study. All in all, large sample size, large
sampling density, accurate maximal distance, and a small tuning parameter
$\tau$ are wished when the performance of all three parameters is considered
at the same time.

The check-based filtering provides a built in comparison of the empirical
overall variance according to the data with the overall varaince indicated by
each fitted bootstrap semi-variogram model. Thereby, it restricts the
bootstrap semi-variogram model to not diverge too much from the data and
empirical semi-variogram. Consequently, parameter estimations are more stable
than without filtering.

###  5.2 Quantile-based filtering

The estimator $\widehat{\eta_{c_{0}}}$ has shown a small standard deviation
and is slightly underestimated. Its performance is similar for all framework
scenarios and tuning parameter specifications. In contrast,
$\widehat{\eta_{\sigma_{0}^{2}}}$ and $\widehat{\eta_{\phi}}$ have large
standard deviations indicating that the corresponding estimators are very
unstable. Moreover, $\widehat{\eta_{\sigma_{0}^{2}}}$ and
$\widehat{\eta_{\phi}}$ are hugely overestimated in all simulation scenarios.

The framework conditions are best for the largest sample size, middle to high
sampling density, and the middle to large maximal distance factor
specification. The effect of sample size and sampling density follow the same
scheme as discussed in case of the check-filtering. The effect of the maximal
distance factor seems unintuitive and is supposedly unreliable because of the
instability of the estimation method.

Regarding the tuning parameter specifications, $\alpha=0.75$ provides best
results when the performance of all three parameters is considered at the same
time. An even stricter $\alpha<0.75$ might lead to further improvement for the
standard deviation and bias measurements for $\widehat{\eta_{\sigma_{0}^{2}}}$
and $\widehat{\eta_{\phi}}$, especially. However, the choice of $\alpha$ seems
arbitrary. It is not bound to the model and data framework conditions in any
way. Further work on the quantile-based filtering approach expected to be
expedient.

###  5.3 Recommendations for action

We recommend the generalized bootstrap method with check-based filtering and a
tuning parameter specification of $\tau\in[1.1,1.5]$ for the uncertainty
estimation of exponential semi-variogam model parameters. Sample size and
density should be preferably large in order to provide enough point pairs
within the model fitting procedure for accurate estimation of nugget effect
and partial sill uncertainties. In addition, the true maximal distance should
be guessed as accurate as possible for the empirical semi-variogram estimation
to improve the shape parameter’s uncertainty estimation.

##  6 Case Study

The generalized bootstrap with check-based filtering is demonstrated on data
obtained from the BaBi birth cohort study [19]. Within the study period from
2013 to 2016 information were gathered about 995 newborn babies and their
families including health status attributes and household locations. Following
the same steps as in the work of Sauzet et al. we use geo-masked coordinates
of the household locations to model the spatial semi-variance of birthweights
[2]. Excluding subjects with incomplete information leaves 783 data points for
the modelling procedure. The geo-masked locations of the households and a
histogram of the pairwise distances are shown in Figure 3.

Figure 3: Plot of geo-masked coordinates (left) and a histogram of all
pairwise Euclidean distances resulting from the coordinates (right).

Besides the choice of the exponential semi-variogram model form, the maximal
distance and number of lag intervals have to be chosen. We recommend fitting
models with varying meta parameter specifications and choose the model with
the best fit. The fit is evaluated by visual inspection and comparing the
value of the loss function given the estimated parameters. Moreover, we
observe the estimated variance according to the semi-variogram model
$\hat{\sigma}^{2}_{model}=\hat{c_{0}}+\hat{\sigma_{0}^{2}}$ compared to the
sample variance $\hat{\sigma}_{sample}^{2}$ by calculating the relative bias
as the ratio

$\displaystyle RB=\frac{\hat{\sigma}^{2}_{model}}{\hat{\sigma}_{sample}^{2}}.$
(18)

First, we try out a range of maximal distances
$(400,500,600,700,800,900,1000)$ and one specific number of lags $(13)$. By
looking at visualizations of the models, we tend to a maximal distance of
$500$. The estimated variance is $2\%$ higher than the sample variance
according to the RB. The loss function takes the minimal value of
$9.08\text{\times}{10}^{11}$ in case of a maximal distance of $500$ compared
to all other implemented options. Second, we investigate a range of lag
interval numbers $(10,11,12,13,14,15)$. Visually, $10$ or $11$ lag intervals
provide the best fit. The comparison between estimated variance is $3\%$
higher provided 10 lags and $2\%$ higher when 11 lags are prespecified. The
loss function takes the two smallest values $5.93\text{\times}{10}^{11}$ given
10 lags and $5.86\text{\times}{10}^{11}$ given 11 lags.

Both Model 1 with 10 lags and Model 2 with 11 lags seem to have a similarly
good fit according to visual exploration. To select the final model, we
estimate the parameter standard errors using the generalized bootstrap with
check-based filtering and tuning parameter $\tau=1.2$. While the parameter
estimates for both models are similar, Model 2 provides smaller standard error
estimates and is selected as final model (Table 6).

Figure 4: Semi-variogram models fitted to birthweight data with a prespecified
maximal distance of 500 and 10 lag intervals (left) and 11 lag intervals
(right).

Semi-variogram model parameter and standard error estimates for birthweight
data. Model 1 refers to the exponential semi-variogram model with a
prespecified maximal distance of 500 and 10 lag intervals. Model 2 refers to
the same model form with a prespecified maximal distance of 500 and 11 lag
intervals.  Model 1 Model 2 Parameter estimate Standard error estimate
Parameter estimate Standard error estimate $\mathbf{c_{0}}$ 102571 50862
125728 48691 $\bm{\sigma_{0}^{2}}$ 127758 48914 103442 47939 $\bm{\phi}$ 63.30
47.67 64.76 43.81

##  7 Summary, Conclusion and Limits

###  7.1 Summary

The goal was to provide an estimator for standard errors of the parameters of
an exponential semi-variogram model. The estimated standard errors of each
semi-variogram model parameter provide a basis for inference and evaluation of
the fitted model that exceeds subjective evaluation by eye.

In this article we have presented two new approaches to estimate the standard
errors of the parameters of an exponential semi-variogram model: The
generalized bootstrap method with either check- or quantile-based filtering.
Both methods were presented and evaluated using a simulation study considering
multiple framework scenarios varying with regards to sample size, sample
density and maximal distance specifications.

###  7.2 Conclusion

Based on the simulation study results, we recommend the generalized bootstrap
method with check-based filtering and tuning parameter $\tau\in[1.1,1.5]$ for
application. For the variance parameters, namely the nugget effect and the
partial sill, it provides reliable standard error estimators. In comparison,
the shape parameter’s standard error estimator is more variable and biased
upwards. Nevertheless, it can serve as rough assessment of the true shape
parameter’s standard error.

###  7.3 Limits

So far, we focused on the exponential semi-variogram with nugget effect since
it is of great interest in the study of neighbourhood correlation. However,
the generalized bootstrap method with filtering is conceptualized for every
kind of parametric semi-variogram model, such as the Gaussian, Whittle or
Spherical model. A validation of the approach via a simulation study for these
parametric model forms is yet to be validated.

The data generating process in the course of the simulation study is a
Gaussian random field. The methodological approach of the generalized
bootstrap is empirical rather than analytical. Because of that, similar
estimator performance is expected for data generated from a non-Gaussian
process. Nevertheless, the filter methods have to be evaluated in case of non-
Gaussian random fields, yet.

## Acknowledgement(s)

We thank all contributors and authors of the BaBi cohort study [19] for the
permission to use the dataset.

## Disclosure statement(s)

The authors report there are no competing interests to declare.

## Data availability statement

Data from the BaBi study was used for demonstration purposes. Because of
ethical restrictions, the data are available only upon request.

## References

  * [1] Bard Y. Nonlinear parameter estimation. New York: Academic Press; 1974. p. 170.
  * [2] Sauzet O, Breiding JH, Zolitschka KA, et al. An ego-centred approach for the evaluation of spatial effects on health in urban areas based on parametric semi-variogram models: concept and validation. BMC medical research methodology. 2021;21(1):1–12.
  * [3] Breckenkamp J, Razum O, Spallek J, et al. A method to define the relevant ego-centred spatial scale for the assessment of neighbourhood effects: the example of cardiovascular risk factors. BMC public health. 2021;21(1):1–7.
  * [4] Cressie NAC. Statistics for spatial data. Rev. ed. ed. New York: Wiley; 1993. Available from: https://onlinelibrary.wiley.com/doi/book/10.1002/9781119115151.
  * [5] Webster R, Oliver MA. Geostatistics for environmental scientists. 2nd ed. Chichester and Hoboken, NJ: Wiley; 2007. p. 119. Statistics in practice.
  * [6] Pardo-Igúzquiza E, Dowd P. Variance-covariance matrix of the experimental variogram: Assessing variogram uncertainty. Mathematical Geology. 2001;33(4):397–419.
  * [7] Olea RA, Pardo-Igúzquiza E. Generalized bootstrap method for assessment of uncertainty in semivariogram inference. Mathematical Geosciences. 2011;43(2):203–228.
  * [8] Saha A, Datta A. Brisc: bootstrap for rapid inference on spatial covariances. Stat. 2018;7(1):e184.
  * [9] Schabenberger O, Gotway CA. Statistical methods for spatial data analysis. Boca Raton [u.a.]: Chapman and Hall/CRC; 2005. Texts in statistical science.
  * [10] Matheron G. Principles of geostatistics. Economic Geology. 1963 12;58(8):1246–1266. Available from: https://doi.org/10.2113/gsecongeo.58.8.1246.
  * [11] Cressie N. Fitting variogram models by weighted least squares. Journal of the international Association for mathematical Geology. 1985;17(5):563–586.
  * [12] Ribeiro PJ, Diggle PJ. Package geor; 2020. Available from: https://CRAN.R-project.org/package$=$geoR.
  * [13] Deutsch CV, Journel AG. Gslib : geostatistical software library and user’s guide. 2nd ed. New York (N.Y.) : Oxford university press; 1998. Chapter V.2 Gaussian-Related Algorithms; p. 141 – 149. Available from: http://lib.ugent.be/catalog/rug01:000429582.
  * [14] Shortridge A. nscore.r ; 2008. Available from: https://msu.edu/ ashton/research/code/nscore.R.
  * [15] Solow AR. Bootstrapping correlated data. Journal of the International Association for Mathematical Geology. 1985 Oct;17(7):769–775. Available from: https://doi.org/10.1007/BF01031616.
  * [16] R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2020. Available from: https://www.R-project.org/.
  * [17] Schlather M, Malinowski A, Oesting M, et al. RandomFields: Simulation and analysis of random fields; 2020. R package version 3.3.8; Available from: https://cran.r-project.org/package=RandomFields.
  * [18] French J. Spatialtools: Tools for spatial data analysis; 2018. R package version 1.0.4; Available from: https://CRAN.R-project.org/package=SpatialTools.
  * [19] Spallek J, Grosser A, Höller-Holtrichter C, et al. Early childhood health in bielefeld, germany (babi study): study protocol of a social-epidemiological birth cohort. BMJ Open. 2017;7(8). Available from: https://bmjopen.bmj.com/content/7/8/e018398.

##  Appendix A Simulation study results

###  A.1 Convergence behaviour of the parameter estimation

Percentage of averaged probable runs grouped by sample size, sample density
and maximal distance factor. The number of all executed runs $n_{sim}=3000$
equals 100%. Parameter estimates are obtained by the WLS approach (sec. 2.2)
grouping parameter % probable WLS runs $\varnothing$ by sample size 500 0.875
1000 0.970 2000 0.999 $\varnothing$ by sample density low 0.929 middle 0.956
high 0.959 $\varnothing$ by maximal distance factor 1.1 0.915 1.5 0.961 2.0
0.968

###  A.2 Performance tables - check filtering for nugget parameter

Check Filter Method - Performance measure results for the nugget effect
$c_{0}$ grouped by sample size, sample density, maximal distance factor and
tuning parameter $\tau$.  grouping parameter $\eta_{c_{0}}$ (Monte Carlo SE)
$\hat{\eta}_{c_{0}}$ (empirical SE) bias MSE $\varnothing$ grouped by sample
size 500 21.39 (0.2961) 22.05 (5.1373) 0.6589 26.8257 1000 11.22 (0.1474)
12.23 (3.3294) 1.0046 12.0937 2000 5.91 (0.0764) 6.05 (0.8837) 0.1428 0.8014
$\varnothing$ grouped by sample density low 15.51 (0.2118) 16.12 (3.8414)
0.6151 15.1350 middle 11.23 (0.1505) 12.00 (2.6954) 0.7728 7.8623 high 11.78
(0.1577) 12.20 (2.8135) 0.4184 8.0911 $\varnothing$ grouped by maximal
distance factor 1.1 13.00 (0.1800) 13.96 (3.2568) 0.9652 11.5385 1.5 12.74
(0.1703) 13.43 (3.0858) 0.6915 10.0005 2.0 12.78 (0.1696) 12.93 (3.0077)
0.1496 9.0687 $\varnothing$ grouped by tuning parameter $\tau$ 1.1 12.84
(0.1733) 13.45 (3.1353) 0.6145 10.2076 1.5 12.84 (0.1733) 13.44 (3.1176)
0.5971 10.0761 2.0 12.84 (0.1733) 13.44 (3.1122) 0.5987 10.0440 2.5 12.84
(0.1733) 13.44 (3.1101) 0.5997 10.0322 3.0 12.84 (0.1733) 13.44 (3.1088)
0.6004 10.0252

###  A.3 Performance tables - check filtering for partial sill

Check Filter Method - Performance measure results for the partial sill
$\sigma^{2}_{0}$ grouped by sample size, sample density, maximal distance
factor and tuning parameter $\tau$.  grouping parameter
$\eta_{\sigma^{2}_{0}}$ (Monte Carlo SE) $\hat{\eta}_{\sigma^{2}_{0}}$
(empirical SE) bias MSE $\varnothing$ grouped by sample size 500 20.99
(0.2906) 22.59 (4.8658) 1.5943 26.2181 1000 10.76 (0.1413) 12.60 (3.2890)
1.8381 14.1962 2000 5.83 (0.0753) 6.35 (1.0435) 0.5232 1.3626 $\varnothing$
grouped by sample density low 15.03 (0.2054) 16.28 (3.7775) 1.2460 15.8221
middle 11.05 (0.1481) 12.65 (2.7174) 1.5984 9.9393 high 11.50 (0.1538) 12.61
(2.7034) 1.1112 8.5432 $\varnothing$ grouped by maximal distance factor 1.1
12.57 (0.1742) 14.79 (3.4024) 2.2185 16.4979 1.5 12.40 (0.1657) 13.69 (2.9227)
1.2943 10.2175 2.0 12.61 (0.1674) 13.05 (2.8732) 0.4428 8.4516 $\varnothing$
grouped by tuning parameter $\tau$ 1.1 12.53 (0.1691) 13.44 (3.0439) 0.9158
10.1039 1.5 12.53 (0.1691) 13.52 (3.0082) 0.9974 10.0439 2.0 12.53 (0.1691)
13.80 (3.0102) 1.2717 10.6785 2.5 12.53 (0.1691) 14.09 (3.0781) 1.5621 11.9152
3.0 12.53 (0.1691) 14.37 (3.1902) 1.8457 13.5838

###  A.4 Performance tables - check filtering for shape parameter

Check Filter Method - Performance measure results for the shape parameter
$\phi$ grouped by sample size, sample density, maximal distance factor and
tuning parameter $\tau$.  grouping parameter $\eta_{\phi}$ (Monte Carlo SE)
$\hat{\eta}_{\phi}$ (empirical SE) bias MSE $\varnothing$ grouped by sample
size 500 154.10 (2.13) 314.33 (402.30) 160.24 $1.88\text{\times}{10}^{05}$
1000 114.11 (1.50) 185.40 (237.55) 71.29 $6.15\text{\times}{10}^{04}$ 2000
55.58 (0.72) 74.89 (39.71) 19.31 $1.95\text{\times}{10}^{03}$ $\varnothing$
grouped by sample density low 115.87 (1.57) 195.56 (230.45) 79.69
$5.95\text{\times}{10}^{04}$ middle 108.90 (1.46) 177.79 (201.92) 68.90
$4.55\text{\times}{10}^{04}$ high 99.02 (1.32) 201.27 (247.18) 102.25
$7.16\text{\times}{10}^{04}$ $\varnothing$ grouped by maximal distance factor
1.1 124.44 (1.70) 147.63 (129.06) 23.19 $1.72\text{\times}{10}^{04}$ 1.5
102.74 (1.37) 163.05 (171.17) 60.31 $3.29\text{\times}{10}^{04}$ 2.0 96.62
(1.28) 263.94 (379.32) 167.32 $1.72\text{\times}{10}^{05}$ $\varnothing$
grouped by tuning parameter $\tau$ 1.1 107.93 (1.45) 114.37 (48.97) 6.44
$2.44\text{\times}{10}^{03}$ 1.5 107.93 (1.45) 148.74 (106.64) 40.81
$1.30\text{\times}{10}^{04}$ 2.0 107.93 (1.45) 189.52 (206.62) 81.59
$4.93\text{\times}{10}^{04}$ 2.5 107.93 (1.45) 230.64 (325.77) 122.71
$1.21\text{\times}{10}^{05}$ 3.0 107.93 (1.45) 274.45 (444.59) 166.52
$2.25\text{\times}{10}^{05}$

###  A.5 Performance tables - quantile filtering for nugget parameter

Quantile Filter Method - Performance measure results for the nugget effect
$c_{0}$ grouped by sample size, sample density, maximal distance factor and
tuning parameter $\alpha$.  grouping parameter $\eta_{c_{0}}$ (Monte Carlo SE)
$\hat{\eta}_{c_{0}}$ (empirical SE) bias MSE $\varnothing$ grouped by sample
size 500 21.39 (0.2961) 20.12 (4.6813) -1.2636 23.5111 1000 11.22 (0.1474)
11.27 (3.1462) 0.0504 9.9013 2000 5.91 (0.0764) 5.27 (0.7742) -0.6412 1.0105
$\varnothing$ grouped by sample density low 15.51 (0.2118) 14.70 (3.5197)
-0.8106 13.0451 middle 11.23 (0.1505) 10.89 (2.4744) -0.3399 6.2384 high 11.78
(0.1577) 11.08 (2.6076) -0.7040 7.2952 $\varnothing$ grouped by maximal
distance factor 1.1 13.00 (0.1800) 12.74 (2.9417) -0.2563 8.7193 1.5 12.74
(0.1703) 12.19 (2.8412) -0.5463 8.3710 2.0 12.78 (0.1696) 11.73 (2.8188)
-1.0517 9.0517 $\varnothing$ grouped by tuning parameter $\alpha$ 0.75 12.84
(0.1733) 11.30 (2.8308) -1.5360 10.3725 0.80 12.84 (0.1733) 11.59 (2.8193)
-1.2495 9.5097 0.85 12.84 (0.1733) 11.91 (2.8261) -0.9343 8.8595 0.90 12.84
(0.1733) 12.27 (2.8525) -0.5650 8.4562 0.95 12.84 (0.1733) 12.74 (2.8969)
-0.0998 8.4022 1.00 12.84 (0.1733) 13.52 (2.9778) 0.6758 9.3242

###  A.6 Performance tables - quantile filtering for partial sill

Quantile Filter Method - Performance measure results for the partial sill
$\sigma^{2}_{0}$ grouped by sample size, sample density, maximal distance
factor and tuning parameter $\alpha$.  grouping parameter
$\eta_{\sigma^{2}_{0}}$ (Monte Carlo SE) $\hat{\eta}_{\sigma^{2}_{0}}$
(empirical SE) bias MSE $\varnothing$ grouped by sample size 500 20.99 (0.29)
5169.66 (5727.38) 5148.67 $5.93\text{\times}{10}^{07}$ 1000 10.76 (0.14)
1434.45 (2203.55) 1423.69 $6.88\text{\times}{10}^{06}$ 2000 5.83 (0.08) 118.78
(314.30) 112.95 $1.12\text{\times}{10}^{05}$ $\varnothing$ grouped by sample
density low 15.03 (0.21) 2860.49 (3151.04) 2845.45
$1.80\text{\times}{10}^{07}$ middle 11.05 (0.15) 2036.60 (2810.21) 2025.55
$1.20\text{\times}{10}^{07}$ high 11.50 (0.15) 1825.80 (2283.98) 1814.31
$8.51\text{\times}{10}^{06}$ $\varnothing$ grouped by maximal distance factor
1.1 12.57 (0.17) 3908.97 (4646.65) 3896.40 $3.68\text{\times}{10}^{07}$ 1.5
12.40 (0.17) 1770.32 (2215.36) 1757.92 $8.00\text{\times}{10}^{06}$ 2.0 12.61
(0.17) 1043.60 (1383.22) 1030.99 $2.98\text{\times}{10}^{06}$ $\varnothing$
grouped by tuning parameter $\alpha$ 0.75 12.53 (0.17) 86.33 (321.57) 73.80
$1.09\text{\times}{10}^{05}$ 0.80 12.53 (0.17) 226.05 (587.95) 213.53
$3.91\text{\times}{10}^{05}$ 0.85 12.53 (0.17) 505.10 (946.79) 492.58
$1.14\text{\times}{10}^{06}$ 0.90 12.53 (0.17) 993.00 (1407.84) 980.47
$2.94\text{\times}{10}^{06}$ 0.95 12.53 (0.17) 1941.15 (2020.74) 1928.62
$7.80\text{\times}{10}^{06}$ 1.00 12.53 (0.17) 9694.16 (11205.58) 9681.64
$2.19\text{\times}{10}^{08}$

###  A.7 Performance tables - quantile filtering for shape parameter

Quantile Filter Method - Performance measure results for the shape parameter
$\phi$ grouped by sample size, sample density, maximal distance factor and
tuning parameter $\alpha$.  grouping parameter $\eta_{\phi}$ (Monte Carlo SE)
$\hat{\eta}_{\phi}$ (empirical SE) bias MSE $\varnothing$ grouped by sample
size 500 154.10 (2.13) $1.96\text{\times}{10}^{05}$
($2.15\text{\times}{10}^{05}$) $1.96\text{\times}{10}^{05}$
$8.50\text{\times}{10}^{10}$ 1000 114.11 (1.50) $5.38\text{\times}{10}^{04}$
($8.21\text{\times}{10}^{04}$) $5.37\text{\times}{10}^{04}$
$9.61\text{\times}{10}^{09}$ 2000 55.58 (0.72) $3.66\text{\times}{10}^{03}$
($1.01\text{\times}{10}^{04}$) $3.60\text{\times}{10}^{03}$
$1.15\text{\times}{10}^{08}$ $\varnothing$ grouped by sample density low
115.87 (1.57) $1.06\text{\times}{10}^{05}$ ($1.11\text{\times}{10}^{05}$)
$1.06\text{\times}{10}^{05}$ $2.35\text{\times}{10}^{10}$ middle 108.90 (1.46)
$7.17\text{\times}{10}^{04}$ ($1.03\text{\times}{10}^{05}$)
$7.16\text{\times}{10}^{04}$ $1.57\text{\times}{10}^{10}$ high 99.02 (1.32)
$7.61\text{\times}{10}^{04}$ ($9.39\text{\times}{10}^{04}$)
$7.60\text{\times}{10}^{04}$ $1.46\text{\times}{10}^{10}$ $\varnothing$
grouped by maximal distance factor 1.1 124.44 (1.70)
$9.62\text{\times}{10}^{04}$ ($1.07\text{\times}{10}^{05}$)
$9.60\text{\times}{10}^{04}$ $2.06\text{\times}{10}^{10}$ 1.5 102.74 (1.37)
$7.74\text{\times}{10}^{04}$ ($9.74\text{\times}{10}^{04}$)
$7.72\text{\times}{10}^{04}$ $1.55\text{\times}{10}^{10}$ 2.0 96.62 (1.28)
$8.04\text{\times}{10}^{04}$ ($1.04\text{\times}{10}^{05}$)
$8.03\text{\times}{10}^{04}$ $1.72\text{\times}{10}^{10}$ $\varnothing$
grouped by tuning parameter $\alpha$ 0.75 107.93 (1.45)
$3.02\text{\times}{10}^{03}$ ($1.24\text{\times}{10}^{04}$)
$2.91\text{\times}{10}^{03}$ $1.63\text{\times}{10}^{08}$ 0.80 107.93 (1.45)
$1.11\text{\times}{10}^{04}$ ($2.46\text{\times}{10}^{04}$)
$1.10\text{\times}{10}^{04}$ $7.25\text{\times}{10}^{08}$ 0.85 107.93 (1.45)
$2.54\text{\times}{10}^{04}$ ($4.01\text{\times}{10}^{04}$)
$2.53\text{\times}{10}^{04}$ $2.24\text{\times}{10}^{09}$ 0.90 107.93 (1.45)
$4.63\text{\times}{10}^{04}$ ($5.69\text{\times}{10}^{04}$)
$4.62\text{\times}{10}^{04}$ $5.37\text{\times}{10}^{09}$ 0.95 107.93 (1.45)
$8.12\text{\times}{10}^{04}$ ($7.49\text{\times}{10}^{04}$)
$8.11\text{\times}{10}^{04}$ $1.22\text{\times}{10}^{10}$ 1.00 107.93 (1.45)
$3.41\text{\times}{10}^{05}$ ($4.06\text{\times}{10}^{05}$)
$3.41\text{\times}{10}^{05}$ $2.81\text{\times}{10}^{11}$