“gmGeostats” is a package for multivariate geostatistics, focusing in the usage of data from multivariate restricted sampling spaces. Such data include positive data, compositional data, distributional data and the like. Most of the times, the geostatistical analysis of such data includes three steps:
The package is loaded, as usual with
library(gmGeostats)
#> Welcome to 'gmGeostats', a package for multivariate geostatistical analysis.
#> Note: use 'fit_lmc' instead of fit.lmcand it fundamentally depends on packages “compositions”, “gstat” and
“sp”. Other dependencies are more instrumental and less fundamental.
NOTE: if you separately need “compositions” or “gstat”, load these
packages first, then load “gmGeostats”. This will ensure that the
overloaded functions work properly for all three packages.
Alternatively, use fully qualified names
(e.g. pkg::foo()).
This vignette very briefly presents the steps to follow in several analysis and modelling routes, illustrated with the case of compositional data. No explanations, theory or discussion is included.
The data can be visually inspected with scatterplots, in raw
representation (plot(), pairs()), in ternary
diagrams (compositions::plot.acomp()), and in scatterplots
of logratio transformed data (use the transformations
pwlr(), alr(), clr() or
ilr() from package “compositions”, then
pairs()). Function pairs() can be given panel
functions such as e.g. vp.lrdensityplot(),
vp.kde2dplot() or vp.lrboxplot() resp. for
creating histograms of pairwise logratios, 2D kernel density maps on the
scatterplots or boxplots of the pairwise logratios. Package
“compositions” provides the class “acomp” to directly deal with the
right representation in the several methods.
Principal component analysis is also a strong help. For this, you
need an isometry (not just an isomorphism). Transformation
clr() is the best for this, and is actually automatically
used when you do princomp(acomp(YOURDATA)). “gmGeostats”
provides generalised diagonalisation methods to account for the spatial
dependence, see ?genDiag for details.
Create your spatial objects by connecting the spatial coordinates to
the multivariate observations via the functions
sp:SpatialPointsDataFrame() or better the “gmGeostats”
functions make.gmMultivariateGaussianSpatialModel() for
multivariate data and
make.gmCompositionalGaussianSpatialModel() for
compositional data. The functions
make.gm******SpatialModel() produces objects of spatial
data container class “gmSpatialModel”, that are necessary for the rest
of the analysis and modelling.
Swath plots are available with command swath(). If you
give it an “acomp” object you will obtain a matrix of logratio swath
plots. Otherwise you will get an set of swath plots, one for each
variable. Function pairsmap() works similarly, but produces
bubble maps (you can control size and color of the symbols).
Empirical variograms can be obtained with function
variogram() out of the spatial data container. You can also
use the function logratioVariogram() for compositional
data. Both accept anisotropy. Their output can be plotted with
plot(), which has specific methods for compositional and
non-compositional data. In the case of anisotropic data, you can also
use a method of image() to visualise the variogram maps,
see ?image.logratioVariogramAnisotropy for details.
Finally you can also check for the strength of the spatial dependence
with the test noSpatCorr.test(). This is a permutations
test, which null hypothesis is that the data do not exhibit spatial
autocorrelation.
Modelling the empirical variograms obtained in the last step can be
done with the function fit_lmc(). This requires specifying
a variogram model, which parameters will be fitted by that function.
Variogram specifications are available with any of the following
functions: gstat::vgm() and
gmGeostats::gmCgram() for multivariate data,
compositions::CompLinModCoReg() and
gmGeostats::LMCAnisCompo() for compositional data.
CompLinModCoReg() is the only one not accepting anisotropy.
You can mix and merge empirical and theoretical models from different
packages, as fit_lmc() will take care to convert between
them for appropriate consistency.
Plotting of LMCs against their empirical variograms can be done with
function variogramModelPlot().
Neighbourhood descriptions are created with function
KrigingNeighbourhood(). Kriging neighbourhoods and LMC
variogram models and can be attached to the “gmSpatialModel” objects at
the moment of creation via make.gm*() functions, using
arguments ng and model (this last one strictly
requiring you to also specify the formula argument).
Validation of the model or of the neighbourhood can then be obtained
with the function validate(). This requires an
object (the complete “gmSpatialModel”) and a
strategy. Validation strategies are small S3-objects
describing what will exactly be done in the validation. They can be
quickly defined by means of configuration functions as
LeaveOneOut() or NfoldCrossValidation(). The
call to validate() will provide some output that can be
evaluated with functions such as xvErrorMeasures(),
accuracy() or prediction().
This way of working is common to the package. You always build a
model (with a make.*() function), define a method parameter
object (created with a specific, verbose, helper function), and feed
both to a common umbrella function describing what do you want to do:
validate() or predict(), the second one also
requires an argument newdata as is standard in R. The
output can then be postprocessed by specific functions.
The method parameter for cokriging is actually the neighbourhood. So,
you can give predict() an object resulting from
KrigingNeighbourhood(), otherwise it will take the standard
one stored in the “gmSpatialModel”, or produce a global neighbourhood
cokriging if no neighbourhood description is found.
Output of predict for “gstat” objects (the current
default) can be re-formed to compositional shape by means of the
function gsi.gstatCokriging2compo(); there is an
gsi.gstatCokriging2rmult() as well for multivariate data.
Maps can be obtained with function image_cokriged(), that
produces a choropleth map with legend, and returns the color scale
(i.e. some breaks and a palette of colors) to be used, e.g. for plotting
the initial data on top of the maps using the same color scale. NOTE:
this funtion does NOT freeze the plot! Most probably you will need to
call par(mfrow=c(1,1)) to create a clean slate device for
the next plot.
Gaussian cosimulation requires joint multivariate normality. The
package provides the flow anamorphosis algorithm for this goal. This is
obtained in two steps. First, you create the transformation by calling
function ana() with your data and storing it. Then, you
apply the stored transformation to the data, and obtain the normalised
scores. These scores can then be treated with all methods and techniques
of the preceding sections “Exploratory analysis” and “Linear model of
coregionalisation (LMC)”, in particular a call to
make.gmMultivariateGaussianSpatialModel() will create the
passing “gmSpatialModel”.
Cosimulation method parameters are created by calls to one of the
functions SequentialSimulation(),
TurningBands() or CholeskyDecomposition(). All
these functions include an argument nsim controlling the
number of realisations desired. These are then obtained by a call to
predict(), giving it the “gmSpatialModel”, the
newdata and the method parameters, in this order.
Multivariate cosimulation output can be seen analogue to a
3-dimensional array, with one dimension running along the simulated
locations, one dimension along the realisations and one dimension along
the variables. This structure is captured in “gmGeostats” with an object
of class “DataFrameStack” (extending “data.frame” and mimicking arrays).
Point-wise, simulation-wise and variable-wise transformations on this
array can be computed with function gmApply(), a wrapper on
base::apply() allowing for an easier management of the
dimensions of the simulation stack. Maps can also be produced by
function image_cokriged().
Multipoint cosimulation is available with the same strategy than
Gaussian based simulation. One needs first to define a “gmSpatialModel”
containing the conditioning data (the original data) and the stochastic
model (the training image). Second, a simulation grid must be created,
and provided to predict() as the newdata
argument. And third, we must provide some method parameters defining the
simulation algorithm to use: currently, only direct sampling is
available, and its parameters can be constructed calling
DSpars().
Training images are currently objects of class “SpatialPixelsDataFrame” or “SpatialGridDataFrame”, from package “sp”. The conditioning data will be migrated to the simulation grid internally; the grid topologies for simulation and training image will be checked for consistency.