Introduction to EgoCor

This is an introduction to the package EgoCor. EgoCor offers a user friendly interface to fit and display a range semi-variogram model graphics and tables of parameters to facilitate the decision making about which exponential parameters fit either raw data or residuals best. The modeling function vario.mod() is based on the functions of the R package gstat. A further function providing the measure of uncertainty proposed by Dyck and Sauzet (2023). has been implemented in the package. With the R package EgoCor modelling the spatial correlation structure of health outcome with a measure of uncertainly is made available to non specialists.

Functions

We use this dataset to illustrate the following functions:

coords.plot(): for graphical description of locations;
distance.info(): for descriptive information about distances between observations;
vario.mod(): to fit empirical semi-variograms and based on that exponential models with graphical presentation;
vario.reg.prep(): to extract the spatial correlation structure of residuals of a regression model;
par.uncertainty(): to obtain bootstrap standard errors for the parameters of the exponential semi-variogram model.

While detailed information can be retrieved with help(function.name), a presentation and demonstration of all functions is provided in the following sections.

1. coords.plot()

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.

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.

coords.plot(birth)

2. distance.info()

This function provides further information about the distribution of pairwise Euclidean distances. It displays the following descriptive statistics:

a histogram of Euclidean distances;
minimum, 1st quartile, median, mean, 3rd quartile and maximum.

distance.info(birth)

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

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. vario.mod()

This function enables the simultaneous output of multiple exponential semi-variogram models fitted for a range of maximal distances and a number of bins. 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.

In the default setting shinyresults = TRUE, an interactive shiny application is opened automatically that displays the results. For the purpose of this vignette though we set shinyresults = FALSE and windowplots = TRUE.

The chosen maximal distance value specifies the subset of data pairs that are incorporated in the semi-variogram estimation. Only data pairs with an Euclidean distance ≤ 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:

mod = vario.mod(birth, max.dist = c(1000,800,600), shinyresults = FALSE, windowplots = TRUE)

#> 
#> Parameter estimates:
#>   max.dist nbins nbins.used    nugget partial.sill     shape prac.range
#> 1     1000    13         13  83189.30    152193.65  70.60961   180.7374
#> 2      800    13         13  88230.81    145523.85  62.85904   158.5179
#> 3      600    13         13 152494.69     91823.68 155.88956   314.4496
#>         RSV rel.bias
#> 1 0.6465789 1.027968
#> 2 0.6225495 1.020856
#> 3 0.3758362 1.066990
#> 
#> Parameter Legend:
#> Each row contains the estimated parameters of the exponential semi-variogram model with the stated maximal distance.
#>    - Index = model number,
#>    - max.dist = maximal distance considered in the empirical semi-variogram estimation,
#>    - nbins = number of bins specified for the empirical semi-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,
#>    - partial.sill = the estimated partial sill, 
#>    - shape = the estimated shape parameter,
#>    - prac.range = the practical range of the exponential semi-variogram model,
#>    - RSV = the relative structured variability, 
#>    - rel.bias = the relative bias between sample variance and estimated variance according to the model.
#>

You can also get the estimated parameters later by

mod$infotable
#>   max.dist nbins nbins.used    nugget partial.sill     shape prac.range
#> 1     1000    13         13  83189.30    152193.65  70.60961   180.7374
#> 2      800    13         13  88230.81    145523.85  62.85904   158.5179
#> 3      600    13         13 152494.69     91823.68 155.88956   314.4496
#>         RSV rel.bias
#> 1 0.6465789 1.027968
#> 2 0.6225495 1.020856
#> 3 0.3758362 1.066990

The maximal distance of 800 m seems to provide a good fit for this dataset. We can now refine the analysis by varying the number of bins or lags. The nbins parameter specifies the number of lags of the empirical semi-variogram to be estimated. A high number of lags might lead to a small within-lag-sample-size and thus to an unstable estimate. Simultaneously, a too small number of lags might lead to a model that does not detect a spatial correlation structure at all.

mod_2 = vario.mod(birth, max.dist = 800, nbins = c(11,12,13), 
                  shinyresults = FALSE, windowplots = TRUE)

#> 
#> Parameter estimates:
#>   max.dist nbins nbins.used    nugget partial.sill    shape prac.range
#> 1      800    11         11 109262.65     126784.8 83.78199   198.9144
#> 2      800    12         12  74966.36     159845.8 64.32218   167.9552
#> 3      800    13         13  88230.81     145523.9 62.85904   158.5179
#>         RSV rel.bias
#> 1 0.5371156 1.030869
#> 2 0.6807390 1.025475
#> 3 0.6225495 1.020856
#> 
#> Parameter Legend:
#> Each row contains the estimated parameters of the exponential semi-variogram model with the stated maximal distance.
#>    - Index = model number,
#>    - max.dist = maximal distance considered in the empirical semi-variogram estimation,
#>    - nbins = number of bins specified for the empirical semi-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,
#>    - partial.sill = the estimated partial sill, 
#>    - shape = the estimated shape parameter,
#>    - prac.range = the practical range of the exponential semi-variogram model,
#>    - RSV = the relative structured variability, 
#>    - rel.bias = the relative bias between sample variance and estimated variance according to the model.
#>

All models’ graphical representations look similar. According to the relative bias model 3 provides a slightly better fit. We choose model 3 with max.dist = 800 and nbins = 13 as final model.

4. vario.reg.prep()

To use vario.mod() to model the spatial correlation structure of residuals, the studentized residuals from a (hierachical) linear regression model can be extracted by vario.reg.prep(). We want to see if adjusting for some predictors of birthweight explain some or all of the spatial correlation structure seen:

res <- lm(formula = birthweight ~ datediff + primiparous + bmi, data = birth)
v.prep = vario.reg.prep(res, data = birth)
models = vario.mod(v.prep, max.dist = c(800,600), shinyresults = FALSE, windowplots = TRUE)

#> 
#> Parameter estimates:
#>   max.dist nbins nbins.used    nugget partial.sill    shape prac.range
#> 1      800    13         13 0.0000000    1.0055208 25.34798   75.93575
#> 2      600    13         13 0.5996697    0.3825593 34.50250   70.82643
#>         RSV  rel.bias
#> 1 1.0000000 1.0020982
#> 2 0.3894807 0.9788857
#> 
#> Parameter Legend:
#> Each row contains the estimated parameters of the exponential semi-variogram model with the stated maximal distance.
#>    - Index = model number,
#>    - max.dist = maximal distance considered in the empirical semi-variogram estimation,
#>    - nbins = number of bins specified for the empirical semi-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,
#>    - partial.sill = the estimated partial sill, 
#>    - shape = the estimated shape parameter,
#>    - prac.range = the practical range of the exponential semi-variogram model,
#>    - RSV = the relative structured variability, 
#>    - rel.bias = the relative bias between sample variance and estimated variance according to the model.
#>

The results point towards a reduced spatial correlation structure with a maximal distance reduced to half the one obtained before adjustment and much less regularity in the empirical semi-variogram.

5. par.uncertainty()

This function provides the filtered bootstrap standard errors for all three exponential model parameters. As the bootstrap can take some time, it is not called by the function vario.mod() directly so that the bootstrap is not executed automatically for all fitted models. This is left to the choice of the user by selecting one specific fitted model, thus saving execution time. The main output of the function is a table with parameter estimates and standard errors for all three parameters of the estimated model. The user can either provide a model estimated by vario.mod() and a model number that specifies which one to use or provide the parameter estimates, the data, the maximum distance and number of bins seperately. The first option is more convenient, therefore we chose it here. The second option enables standard error estimation for exponential semi-variogram models for which no result of class "vario.mod.output" is available. This might be useful if a fitted model which was not obtained with EgoCor (eg. estimated directly with gstat functions or other software making use of the same estimation algorithm) shall be investigated with respect to parameter standard errors. It is recommended to use more bootstrap samples than done here.

unc = par.uncertainty(mod_2, mod.nr = 2, B = 100)
#> Two approaches regarding the input arguments:
#>  1. Provide the arguments
#>     - vario.mod.output (output object from vario.mod function),
#>     - mod.nr (number of the model in the vario.mod.output$infotable). 
#>  2. Provide the arguments
#>     - par.est (vector with estimated nugget, partial sill and shape parameters),
#>     - data (used to estimate the semi-variogram model parameters),
#>     - max.dist (semi-variogram parameter, numeric of length 1),
#>     - nbins (semi-variogram parameter, numeric of length 1).
#>  In both cases a threshold factor can be set. If not specified a default value of 1.2 is used.
#> 
#> 
#>  Bootstrap started.
#>  This can take a few minutes depending on the number
#>  of bootstrap samples B to be generated.
#> 
#> Initial estimation finished
#>  86 of 100 Models converged.
#> Reestimating:
#> 
 98 of 100 models converged.
 100 of 100 models converged.
#>  Call:  par.uncertainty(vario.mod.output = mod_2 , mod.nr = 2 , par.est = NULL , data= NULL , max.dist= NULL , nbins = NULL , B = 100 , threshold.factor = 3 ) 
#> 
#>                   Estimate Std. Error
#> nugget effect  74966.36163 156852.228
#> partial sill  159845.81359 469168.308
#> shape             64.32218    512.258

Dyck, Julia, and Odile Sauzet. 2023. “Parameter Uncertainty Estimation for Exponential Semivariogram Models: Two Generalized Bootstrap Methods with Check- and Quantile-Based Filtering.” Preprint. https://arxiv.org/abs/2202.05752.