# EgoCor: an R package to fit exponential semi-variograms to model the local
spatial correlation structure of health outcomes

Julia Dyck
Bielefeld University &Jan-Ole Koslik
Bielefeld University &Odile Sauzet
Bielefeld University<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
EgoCor: an R package to fit exponential semi-variograms to model the local
spatial correlation structure of health outcomes R package EgoCor As an
alternative to using administrative areas for the evaluation of small-area
health inequalities, Sauzet et al suggested to take an ego-centred approach
and model the spatial correlation structure of health outcomes at individual
level. Existing tools for the analysis of spatial data in R may appear too
complex to non-specialists which may limit the use of the approach. We present
the R package EgoCor which offers a user-friendly interface displaying in one
function a range of graphics and tables of parameters to facilitate the
decision making about which exponential parameters fit best either raw data or
residuals. This function is based on the functions of the R package gstat.
Moreover, we implemented a function providing the measure of uncertainty
proposed by Dyck and Sauzet. With the R package EgoCor the modelling of
spatial correlation structure of health outcomes with a measure of uncertainty
is made available to non specialists. R package, semi-variogram, exponential
models, small-area health inequalities R package, semi-variogram, exponential
models, small-area health inequalities Odile Sauzet
School of Public Health _and_
Department of Business Administration and Economics
Bielefeld University
Bielefeld, Germany
E-mail:

Julia Dyck
Department of Business Administration and Economics
Bielefeld University
Bielefeld, Germany
E-mail:

Jan-Ole Koslik
Department of Business Administration and Economics
Bielefeld University
Bielefeld, Germany
E-mail:

##  1 Background

The last 20 years has seen an increase in interest in the study of
associations between neighbourhood characteristics and health. Most of the
quantitative studies are based on neighbourhood defined as disjoint
administrative spatial units and non-measured spatial effects are estimated
via multilevel models. However, there are limits to this approach, which have
been discussed in the literature (Chaix _et al._ , 2009; Van Ham and Manley,
2012). An alternative approach is to conceptualise neighbourhood as ego-
centred such that everyone has its own neighbourhood centred on the place of
residence. Such a neighbourhood must have a relevance for the health outcome
of interest.

Some studies mentioned the use of the spatial correlation structure of health
outcomes as an alternative to estimating the correlation within administrative
areas (Chaix _et al._ , 2005). This idea has been brought further by Sauzet et
al (Sauzet _et al._ , 2021) by proposing to use the parameters of an
exponential model for the semi-variogram of health outcomes to quantitatively
assess this correlation structure. Such models provide a measure of the
presence of unmeasured spatial effects on health. It has the advantage of
fitting empirical data well as well as modelling the following concept: if the
place of residence has an effect on health, then health outcomes of neighbours
are correlated and this correlation evanesce with increasing distance between
neighbours. The procedure of investigating the spatial correlation includes
first estimating an empirical semi-variogram based on the data available and
then fit an exponential parametric model to this semi-variogram. Here only
small distances are considered to model the local correlation structure, i.e.
the correlation of health outcome between immediate neighbours.

The empirical semi-variogram is not unique and is based on the choice of the
maximal distance between observations (all pairs of observations which are
further apart than the maximal distance are not used to estimate the semi-
variogram) and the number of bins (distance intervals for pairs of
observations for which one point of the semi-variogram will be estimated). If
the number of observations is limited and only a small proportion of the
variance is spatially structured, then the ability of fitting an exponential
model to the semi-variogram may be strongly dependent on the mentioned meta
parameters. Moreover, the fit of the estimated exponential model must be
evaluated visually with a comparative evaluation of a range of maximal
distances and numbers of bins.

There are a number of possibilities to fit semi-variogram models using R
packages. But the range of modelling possibilities can be a deterrent to its
wider use by health researchers with limited experience with the analysis of
spatial data. Moreover, it remains presently difficult to obtain measures of
uncertainty for parametric semi-variogram models and as far as we know the
only method that has been implemented in current R packages is the BRISC
method by Saha and Datta (2018) provided in the R package BRISC (Saha and
Datta, 2022) which is not suitable to be used with health data (Dyck and
Sauzet, 2023). Dyck and Sauzet (2023) have investigated how to modify an
existing bootstrap approach to make it reliable in the context of health data
survey (sparse local data defined as a small number of pairs at small
distances, and large overall sample size) based on the work of Olea and Pardo-
Igúzquiza (2011). A filtered bootstrap estimate for the standard error of
parameter estimates of an exponential semi-variogram model has been proposed
and evaluated.

The aim of the R package EgoCor is to offer a user-friendly interface
displaying in one function a range of graphics and tables of parameters to
facilitate the decision making about which exponential semi-variogram model
parameters fit the data best together with a measure of uncertainty not
available until now. Our package is based on the functions of the R package
gstat (R Core Team, 2020; Pebesma and Wesseling, 1998) and on the work of Dyck
and Sauzet (2023). A measure of uncertainty for the model parameters can be
obtained from raw data as well as from the residuals from regression models
directly for adjusted analyses.

##  2 Statistical Methods

###  2.1 Semi-variogram model

We provide some basic concepts on semi-variogram modelling. For more
background we refer to Schabenberger and Gotway (2017). In the sequel, we
assume that the spatial process is isotropic and second-order stationary.
Under those conditions (Schabenberger and Gotway, 2017), the correlation
$C(h)$ between health outcomes $Z(s)$ and $Z(s^{\prime})$ (or the residuals of
a regression model) at locations $s$ and $s^{\prime}$ depends only on the
(lag) distance $h=||s-s^{\prime}||$ between those observations.

The empirical semi-variogram $\gamma(h)$ defined for a positive lag $h$ as

$\gamma(h)=\frac{1}{2}Var[Z(s)-Z(s+h)]$

is estimated from the data and is provided by the function variogram from the
R package gstat (Pebesma and Wesseling, 1998).

The covariance function $C(h)$ under some regularity conditions is given by

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

where $\sigma^{2}_{0}+c_{0}=Var[Z(s)]$ is the variance of the health outcome
and the nugget effect $c_{0}$ is the value of the semi-variogram when the
distance between two observations tend to 0. The parameter $\sigma_{0}^{2}$ is
called the partial sill.

The Matheron’s estimator provides an unbiased estimator for the empirical
semi-variogram at distance h between two observations (Matheron, 1962)

$\hat{\gamma}(h)=\frac{1}{2|N(h)|}\sum_{(s,s^{\prime})\in
N(h))}\\{Z(s)-Z(s^{\prime})\\}$

where $N(h)$ is the set of all observations lagging at distance $h$.

To the empirical semi-variogram we fit an exponential model

$\displaystyle\hat{\gamma}_{exp}(h)=\begin{cases}\hat{c}_{0}+\hat{\sigma}_{0}^{2}\Big{(}1-\exp\big{(}-\frac{h}{\hat{\phi}}\big{)}\Big{)}\;\;&\text{for}\;\;h>0,\\\
0\;\;&\text{for}\;\;h=0.\end{cases}$

The practical range (distance above which the correlation between observations
is less than 5% of the total variance) for this model is given by

$H=\hat{\phi}\log\left(\frac{\hat{\sigma}_{0}^{2}}{0.05(\hat{c}_{0}+\hat{\sigma}_{0}^{2})}\right).$

The relative structured variability (RSV) calculated as partial sill divided
by total variance is a measure of the degree of spatial structure:

$RSV=\frac{\hat{\sigma_{0}}^{2}}{\hat{c_{0}}+\hat{\sigma_{0}}^{2}}.$

The relative bias between the estimated variance according to the model and
the sample variance of the health outcome is obtained as

$RB=\frac{\hat{c_{0}}+\hat{\sigma_{0}}^{2}}{\widehat{Var[Z(s)]}}.$

###  2.2 Filtered bootstrap standard error

The algorithm used to obtain standard errors for the parameters of the
exponential semi-variogram model is based on the generalized bootstrap method
explained in Olea and Pardo-Igúzquiza (2011) and Pardo-Igúzquiza and Olea
(2012) and was adapted in Dyck and Sauzet (2023) to the case of the
characteristics of population data: large sample sizes over all, low number of
pairs at small distances. The estimation of the exponential model parameter
uncertainties of the fitted semi-variogram model is obtained by weighted least
squares. A filter was set up within the bootstrapping process to remove all
bootstrap parameter estimates for which the estimation algorithm did not
converge (due to the small number of pairs at small distances). We briefly
recall the steps of the filtered bootstrap algorithm:

  1. 0.

Exponential semi-variogram model: An empirical semi-variogram is fitted to the
original spatial dataset to which an exponential model is fitted. This
provides 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}=(z_{1},\dots,z_{N})^{t}$ is mapped into a Gaussian space by the
empirical normal score transformation function $\varphi$. Consequently,
$\mathbf{y}=\varphi(\mathbf{z})$ is a realization vector of a standard normal
random variable (Deutsch and Journel, 1998).

  3. 2.

Exponential semi-variogram model for transformed data: An empirical semi-
variogram and exponential semi-variogram model $\tilde{\gamma}_{\text{exp}}$
are fitted to the transformed data $\mathbf{y}$ combined with the raw data’s
geo-coding providing the parameter estimate $\mathbf{\tilde{\theta}}$.

  4. 3.

Covariance estimation: Making use of the exponential semi-variogram model
characterized by $\mathbf{\tilde{\theta}}$, the covariance between two data
points $z_{i}$ and $z_{j}$ is calculated based on the Euclidean distance
$d_{ij}$ between these two points:

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

The 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 a lower and triangular matrix $\mathbf{L}$
and its transpose is obtained by the Cholesky decomposition algorithm as
$\mathbf{C}=\mathbf{L}\mathbf{L}^{t}.$

This decomposition is used to remove the correlation structure within the
sample $\mathbf{y}$. The resulting vector
$\mathbf{x}=\mathbf{L}^{-1}\mathbf{y}$ contains independent and identically
distributed, hence uncorrelated values (Solow, 1985).

  6. 5.

Classical bootstrap: Sampling with replacement from $\mathbf{x}$ leads to a
bootstrap sample $\mathbf{x^{*}}$ of the same size as the original spatial
dataset.

  7. 6.

Recorrelation: The resample $\mathbf{x^{*}}$ reinherits the correlation
structure by applying the inverse operation of step 4, that is
$\mathbf{y^{*}}=\mathbf{Lx^{*}}.$

  8. 7.

Normal score back transformation: The back transformation of $\mathbf{y^{*}}$
to the attribute space through the inverse normal score function is done by
$\mathbf{z^{*}}=\varphi^{-1}(\mathbf{y^{*}}).$

  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^{*}}$.

  10. 9.

Filtering: A check-filter based test is applied to the bootstrap estimate
$\mathbf{\theta^{*}}$ to indicate whether the exponential semi-variogram
fitting algorithm did converge. If
$c_{0}^{*}+\sigma_{0}^{2*}>\tau\widehat{Var(\mathbf{z})},$ i.e. if the
variance indicated by the model exceeds the estimated sample variance times
the threshold factor $\tau$, the bootstrap estimate is discarded. Otherwise it
is saved. Within the EgoCor package the threshold is set to $\tau=3$ by
default as in a simulation study it was found to provide the best results with
respect to the standard error estimates (Dyck and Sauzet, 2023).

  11. 10.

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

  12. 11.

Parameter standard error estimation: Based on the set of repeatedly estimated
parameters ${\\{\theta^{*}_{b}\\}}_{b=1,\dots,B}$ the parameters’ standard
error estimates are obtained as

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

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

##  3 The EgoCor package

In this section we describe functions provided by the package. We then apply
those functions to the simulated dataset birth.

###  3.1 Dataset

The simulated dataset birth is provided with the package EgoCor. The dataset
is based on the spatial distribution of real birthweight data (Spallek _et
al._ , 2017). It contains eight variables for 903 births:

  * •

x: x-coordinate in meters for a fictive Cartesian coordinate system,

  * •

y: y-coordinate in meters,

  * •

birthweight: birthweight in grams,

  * •

primiparous: first pregnancy (1) or subsequent pregnancy (0),

  * •

datediff: number of days to due date,

  * •

bmi: BMI of the mother at first medical appointment,

  * •

weight: weight of the mother at first medical appointment,

  * •

inc: income quintile.

###  3.2 Functions

We use the birth dataset to illustrate the following functions:

  * •

coords.plot(): for graphical description of locations,

  * •

distance.info(): for descriptive information about distances between
observations,

  * •

vario.reg.prep(): to model the spatial correlation structure of residuals of a
regression model,

  * •

vario.mod(): to fit exponential models to semi-variograms with graphical
presentation,

  * •

par.uncertainty(): to obtain bootstrap standard errors for the parameters of
the exponential semi-variogram model.

###  3.3 Data exploration

The data format required by the EgoCor functions is either a data frame or a
matrix.

The first three columns of the data frame or matrix should be ordered the
following way:

1st column: x-coordinate in meters for a Cartesian coordinate system; 2nd
column: y-coordinate in meters for a Cartesian coordinate system; 3rd column:
outcome of interest. Other columns will be ignored. A message appears
following the output of the function vario.mod() recalling the required order
for the variables.

    coords.plot(birth)

Figure 1:  Coordinates plot for outcome birthweight

The function coords.plot() provides a simple visualization of the locations on
a two dimensional map and indicates whether the outcome is observed (by a
black circle) or missing (by a red x) at a specific location. The purpose of
this function is to look at the spatial distribution of observations and if
there might be a spatial pattern in the distribution of missing values in the
outcome of interest or in covariates.

Figure 1 displays the location of the observations. For the outcome
birthweight there are no missing values.

To illustrate the display of missing values we first created a new matrix with
the coordinates and the variable inc and then inserted 30 random missing
values.

    x <- sample(1:903, 30)
    datainc <- birth[, c(1:2, 8)]
    datainc[x, 3] <- NA
    coords.plot(datainc)

Figure 2: Coordinates plot for outcome inc with 30 random missing values

In Figure 2 the missing values are marked with red crosses. As expected no
spatial pattern is visible.

Further information about the distribution of pairwise Euclidean distances is
provided by the function distance.info(). It calculates

  * •

the distance matrix containing all pairwise Euclidean distances,

  * •

the set of all pairwise Euclidean distances where duplicate values due to
symmetry are deleted.

Moreover, distance.info() displays the following descriptive statistics:

  * •

a histogram of the Euclidean distances,

  * •

minimum, 1st quartile, median, mean, 3rd quartile and maximum of the Euclidean
distances.

    distance.info(birth)

    ##  Summary of distance set:
    ##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    ##       0    4244    6970    7506   10221   25963

Figure 3: Histogram of pairwise distances for the birth data

The output for the birth data is as follows and illustrated with the histogram
shown in Figure 3: From all the 815 409 pairwise distances, 30 570 are of less
than 2 000 meters and will be used for modelling of the local spatial
correlation structure.

###  3.4 Semi-variogram model fitting

The function vario.mod() enables the simultaneous output of multiple
exponential semi-variogram models fitted for a range of maximal distances and
bin numbers. Thereby, the focus lies on the ability of the function to provide
multiple estimation results depending on various specifications for the meta
parameters max.dist and nbins.

It is advised to try out different values for both parameters and choose the
model with the best fit. Commonly, the fit is evaluated by visual checks. An
additional check can be performed by comparing the sample variance with the
estimated variance according to the semi-variogram model
$\hat{\sigma^{2}}=\hat{c_{0}}+\hat{\sigma_{0}^{2}}$ (Sauzet _et al._ , 2021).

The chosen maximal distance value specifies the subset of data pairs that are
actually used for the semi-variogram estimation. Only data pairs with an
Euclidean distance $\leq$ max.dist are taken into account. For a first
exploration, it might be useful to try a range of maximal distances to locate
where the range might be situated by

    vario.mod(birth, max.dist = c(2000, 1500, 1000, 500), nbins=13, pdf = T,
    pdf.directory = getwd(), pdf.name = "Birthweight")

Figure 4: Shiny ouput from vario.mod(birth, max.dist = c(2000,1500,1000,500),
nbins=13)

The above code will save a PDF file showing all fitted semi-variograms and
will produce the shiny (Chang _et al._ , 2022) output shown in Figure 4.

Each row of the printed output table contains the estimated parameters of the
exponential semi-variogram model with one of the stated maximal distances.
More precisely, the table columns contain:

  * •

index: model number,

  * •

max.dist: maximal distance used in the estimation of the empirical variogram,

  * •

nbins: number of bins specified for the empirical variogram estimation,

  * •

nbins.used: number of bins used for the empirical semi-variogram estimation
(can differ from nbins in case of colocatted data points),

  * •

nugget: the estimated nugget effect $\hat{c_{0}}$,

  * •

partial.sill: the estimated partial sill $\hat{\sigma}_{0}^{2}$,

  * •

shape: the estimated shape parameter $\hat{\phi}$,

  * •

prac.range: the practical range of the exponential model,

  * •

RSV: the relative structured variability,

  * •

rel.bias: the relative bias between the sum of the estimated partill sill and
nugget and the sample variance (which theoretically are the same).

The maximal distance of 1000 meters seems to provide the best fit among the
tried and we can now refine the analysis by considering smaller maximal
distances

    vario.mod(birth, max.dist = c(1000, 800, 600), nbins=13)

Figure 5: Ouput from vario.mod(birth, max.dist = c(1000, 800, 600), nbins=13)

leading to the output presented in Figure 5.

Because a maximal distance of 800 meters provides the best fit for the
exponential model with a low relative bias and a good visual fit, we
investigate further the role of the number of bins for this maximal distance.

The nbins parameter specifies the number of lags of the empirical semi-
variogram to be estimated. On the one hand, a high number of lags might lead
to small within-lag-sample-size and thus to an unstable estimate. On the other
hand, a too small number of lags might lead to a model, that does not detect a
spatial correlation structure at all. To decide on one or multiple values for
nbins, taking a look at the histogram plot of the pairwise distances (see
Figure 3) obtained by distance.info() may help.

Trying out multiple nbins specifications by

    vario.mod(birth, max.dist = 800, nbins = c(11, 12, 13))

Figure 6: Ouput from vario.mod(birth, max.dist = 800, nbins = c(11,12,13))

we obtain the output presented in Figure 7. All models provide similar results
but model 3 with max.dist = 800 and nbins = 13 gives a slightly better fit and
could be selected as final model with respect to the health outcome
birthweight.

###  3.5 Modelling the spatial correlation structure of residuals

Instead of modelling the correlation structure of a health outcome, the
vario.mod() function can be used to model the spatial correlation structure of
residuals from a (hierarchical) linear regression. To do so, the studentized
residuals from a (hierarchical) linear regression model are extracted via the
vario.reg.prep() function.

We want to investigate if adjusting for some predictors of birthweight might
explain some or all of the observed spatial correlation structure. In the
first step, we fit the following regression model and investigate the output:

    res <- lm(birthweight ~ datediff + primiparous + bmi, data = birth)
    summary(res)

    ##
    ## Call:
    ## lm(formula = birthweight ~ datediff + primiparous + bmi, data = birth)
    ##
    ## Residuals:
    ##      Min       1Q   Median       3Q      Max
    ## -1109.92  -274.10   -14.14   260.87  1373.98
    ##
    ## Coefficients:
    ##             Estimate Std. Error t value Pr(>|t|)
    ## (Intercept) 3402.687     59.886  56.819  < 2e-16 ***
    ## datediff     -24.217      1.444 -16.773  < 2e-16 ***
    ## primiparous -108.669     28.424  -3.823 0.000141 ***
    ## bmi            6.551      2.506   2.614 0.009092 **
    ## ---
    ## Signif. codes:  0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
    ##
    ## Residual standard error: 415.9 on 899 degrees of freedom
    ## Multiple R-squared:  0.247,ΨAdjusted R-squared:  0.2445
    ## F-statistic: 98.32 on 3 and 899 DF,  p-value: < 2.2e-16

All predictors are significant. Using the vario.reg.prep() function we now
assign the studentized residuals of the regression model to the spatial
coordinates and investigate which maximal distance provides the best
exponential semi-variogram model fit for the studentized residuals.

We start at similar distances to the ones chosen for the raw data obtaining
the output shown in Figure 7. Then we opt out of the shiny display option to
save the output of vario.mod() in models for future use:

    v.prep <- vario.reg.prep(res, data = birth)
    models <- vario.mod(v.prep, max.dist = c(1000,800,600),
                        nbins = 13, shinyresults = FALSE)

Figure 7: Ouput from vario.mod(v.prep, max.dist = c(1000,800,600), nbins = 13)

The results point towards a reduced spatial correlation structure with a well
fitting maximal distance reduced to only 600 meters and much less regularity
in the empirical semi-variogram (see Figure 7).

Analoguely to the unadjusted case (see section 3.4), we try out multiple nbins
values given the maximal distance set to 600 meters:

    vario.mod(v.prep, max.dist = 600, nbins = c(11, 12, 13))

Figure 8: Ouput from vario.mod(v.prep, max.dist = 600, nbins = c(11, 12, 13))

Based on the resulting graphics and table (see Figure 8) we come to the
conclusion that the models with max.dist = 600 and nbins = 12 (see Figure 8)
or max.dist = 600 and nbins = 13 (see Figure 7) provide visually very similar
and the best fits.

###  3.6 Filtered bootstrap standard errors

The function par.uncertainty() provides filtered bootstrap standard errors for
all three exponential model parameters. Standard errors are important to
conduct proper inference (Bard, 1974). Moreover, they can be helpful to get an
impression about the reliability of the estimated model and provide an
objective tool to compare two or more models which seem to provide an equally
good fit when evaluated visually as demonstrated in the last section 3.5.

Because the execution of the filtered bootstrap algorithm can take some time
(depending on the sample size and number of bootstrap repetitions), the
par.uncertainty() function is not automatically called within vario.mod() so
that the bootstrap is not executed in all models estimated by vario.mod().
This is left to the choice of the user by selecting the model number in the
option mod.nr and thereby saving execution time.

Due to the bootstrap component within the standard error calculation,
repeatedly estimating the standard errors for one fixed semi-variogram model
will slightly vary. If the filtered bootstrap results are wished to be
reproducible, a seed has to be set prior to the application of the
par.uncertainty() command.

We save the vario.mod.output object containing the two best semi-variogram
models of the residuals according to visual inspection:

    models <- vario.mod(v.prep, max.dist = 600, nbins = c(12, 13))

Based on that, we can estimate the parameter standard errors for model 1
(max.dist = 600, nbins = 12):

    unc1 <- par.uncertainty(models, mod.nr = 1, threshold.factor = 3)
    unc1$unc.table

    ##                 Estimate  Std. Error
    ## nugget effect  0.5575581   0.3177576
    ## partial sill   0.4275874   0.5015041
    ## shape         42.6818620 166.6720236

and model 2 (max.dist = 600, nbins = 13):

    unc1 <- par.uncertainty(models, mod.nr = 2, threshold.factor = 3)
    unc2$unc.table

    ##                 Estimate  Std. Error
    ## nugget effect  0.5996697   0.3172093
    ## partial sill   0.3825593   0.5561611
    ## shape         34.5024985 203.0448548

According to the above R output the standard error estimates of model 1 appear
to be slightly smaller than the standard error results of model 2 confirming
our impression that the models are very similar. The concrete standard error
values allow us to make an objective comparison and encourage us to select
model 1 as final model.

##  4 Conclusion

The R package EgoCor proposes a range of functions to explore which empirical
semi-variogram meta parameters (maximal distance and number of bins) provide
the best fitting exponential semi-variogram model. Practitioners can easily
use the package to apply the concept of correlation neighbourhood of Sauzet
_et al._ (2021) which is an ego-centred approach to assess the existence of
unmeasured neighbourhood effects. Moreover, the package provides an
implementation of the filtered bootstrap method for the estimation of standard
errors of Dyck and Sauzet (2023) based on the work of Olea and Pardo-Igúzquiza
(2011); Pardo-Igúzquiza and Olea (2012) which has been adapted to the
characteristics of survey data: large sample sizes and low number of
observations at close distances.

## Computational details

The results in this paper were obtained using R 4.2.1 with the EgoCor 1.1.0
package. R itself and all packages used are available from the Comprehensive R
Archive Network (CRAN) at https://CRAN.R-project.org/.

## References

  * Bard (1974) Bard Y (1974).  _Nonlinear parameter estimation_.  Academic Press, New York.  ISBN 0-12-078250-2.
  * Chaix _et al._ (2009) Chaix B, Merlo J, Evans D, Leal C, Havard S (2009).  “Neighbourhoods in eco-epidemiologic research: delimiting personal exposure areas. A response to Riva, Gauvin, Apparicio and Brodeur.”  _Social science & medicine_, 69(9), 1306–1310.
  * Chaix _et al._ (2005) Chaix B, Merlo J, Subramanian S, Lynch J, Chauvin P (2005).  “Comparison of a spatial perspective with the multilevel analytical approach in neighborhood studies: the case of mental and behavioral disorders due to psychoactive substance use in Malmö, Sweden, 2001.”  _American journal of epidemiology_ , 162(2), 171–182.
  * Chang _et al._ (2022) Chang W, Cheng J, Allaire J, Sievert C, Schloerke B, Xie Y, Allen J, McPherson J, Dipert A, Borges B (2022).  _shiny: Web Application Framework for R_.  R package version 1.7.4, URL https://CRAN.R-project.org/package=shiny.
  * Deutsch and Journel (1998) Deutsch CV, Journel AG (1998).  _GSLIB : geostatistical software library and user’s guide_ , chapter V.2 Gaussian-Related Algorithms, pp. 141 – 149.  2nd ed. edition. New York (N.Y.) : Oxford university press.  ISBN 0195100158.  URL http://lib.ugent.be/catalog/rug01:000429582.
  * Dyck and Sauzet (2023) Dyck J, Sauzet O (2023).  “Parameter uncertainty estimation for exponential semi-variogram models: Two generalized bootstrap methods with check-and quantile-based filtering.”  _arXiv preprint_.  URL https://arxiv.org/pdf/2202.05752.
  * Matheron (1962) Matheron G (1962).  _Traité de géostatistique appliquée. 1 (1962)_ , volume 1.  Editions Technip.
  * Olea and Pardo-Igúzquiza (2011) Olea RA, Pardo-Igúzquiza E (2011).  “Generalized Bootstrap Method for Assessment of Uncertainty in Semivariogram Inference.”  _Mathematical Geosciences_ , 43(2), 203–228.  10.1007/s11004-010-9269-6.
  * Pardo-Igúzquiza and Olea (2012) Pardo-Igúzquiza E, Olea RA (2012).  “VARBOOT: A spatial bootstrap program for semivariogram uncertainty assessment.”  _Computers & Geosciences_, 41, 188–198.
  * Pebesma and Wesseling (1998) Pebesma EJ, Wesseling CG (1998).  “Gstat: a program for geostatistical modelling, prediction and simulation.”  _Computers & Geosciences_, 24(1), 17–31.  10.1016/s0098-3004(97)00082-4.
  * R Core Team (2020) R Core Team (2020).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL https://www.R-project.org/.
  * Saha and Datta (2018) Saha A, Datta A (2018).  “BRISC: bootstrap for rapid inference on spatial covariances.”  _Stat_ , 7(1), e184.
  * Saha and Datta (2022) Saha A, Datta A (2022).  _BRISC: Fast Inference for Large Spatial Datasets using BRISC_.  R package version 1.0.5, URL https://CRAN.R-project.org/package=BRISC.
  * Sauzet _et al._ (2021) Sauzet O, Breiding JH, Zolitschka KA, Breckenkamp J, Razum O (2021).  “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_ , 21(1), 1–12.
  * Schabenberger and Gotway (2017) Schabenberger O, Gotway CA (2017).  _Statistical methods for spatial data analysis_.  CRC press.
  * Solow (1985) Solow AR (1985).  “Bootstrapping correlated data.”  _Journal of the International Association for Mathematical Geology_ , 17(7), 769–775.  ISSN 1573-8868.  10.1007/BF01031616.
  * Spallek _et al._ (2017) Spallek J, Grosser A, Höller-Holtrichter C, Doyle IM, Breckenkamp J, Razum O (2017).  “Early childhood health in Bielefeld, Germany (BaBi study): study protocol of a social-epidemiological birth cohort.”  _BMJ Open_ , 7(8).  ISSN 2044-6055.  10.1136/bmjopen-2017-018398.  https://bmjopen.bmj.com/content/7/8/e018398.full.pdf, URL https://bmjopen.bmj.com/content/7/8/e018398.
  * Van Ham and Manley (2012) Van Ham M, Manley D (2012).  “Neighbourhood Effects Research at a Crossroads. Ten Challenges for Future Research Introduction.”  _Environment and Planning A_ , 44(12), 2787–2793.