| georob-package {georob} | R Documentation |
The georob Package
Description
This is a summary of the features and functionality of georob, a package in R for customary and robust geostatistical analyses.
Details
georob is a package for customary and robust analyses of
geostatistical data.
Such data, say y_i=y(\boldsymbol{s}_i), are
recorded at a set of locations,
\boldsymbol{s}_i, i=1,2, \ldots, n, in a
domain G \in \mathrm{I}\!\mathrm{R}^d, d \in (1,2,3), along
with covariate information
x_j(\boldsymbol{s}_i), j=1,2,
\ldots, p.
Model
We use the following model for the data
y_i=y(\boldsymbol{s}_{i}):
Y(\boldsymbol{s}_i) =
Z(\boldsymbol{s}_i) + \varepsilon =
\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}
\boldsymbol{\beta} +
B(\boldsymbol{s}_i) +
\varepsilon_i,
where
Z(\boldsymbol{s}_i)=\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}
\boldsymbol{\beta} +
B(\boldsymbol{s}_i) is the so-called signal,
\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}
\boldsymbol{\beta}
is the external drift,
\{B(\boldsymbol{s})\} is an unobserved stationary or
intrinsic spatial Gaussian random field with zero mean, and
\varepsilon_i is an
i.i.d error from a possibly long-tailed distribution with scale parameter
\tau (\tau^2 is usually called nugget effect).
In vector form the model is written as
\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{B} + \boldsymbol{\varepsilon},
where \boldsymbol{X} is the model matrix with the
rows
\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}.
The (generalized) covariance matrix of the vector of
spatial Gaussian random effects
\boldsymbol{B}
is denoted by
\mathrm{E}[
\boldsymbol{B}\,
\boldsymbol{B}^\mathrm{T}] =
\boldsymbol{\Gamma}_\theta =
\sigma_{\mathrm{n}}^2\boldsymbol{I} +
\sigma^2\boldsymbol{V}_\alpha =
\sigma_B^2 \, \boldsymbol{V}_{\alpha,\xi} =
\sigma_B^2 \, ((1-\xi) \, \boldsymbol{I} +
\xi\, \boldsymbol{V}_\alpha )
,
where
\sigma_{\mathrm{n}}^2
is the variance of seemingly uncorrelated micro-scale variation in
B(\boldsymbol{s})
that cannot be resolved with the chosen sampling design,
\boldsymbol{I} is the identity matrix,
\sigma^2 is the variance of the captured auto-correlated variation in
B(\boldsymbol{s}),
\sigma_B^2=\sigma_{\mathrm{n}}^2+\sigma^2
is the signal variance, and
\xi=\sigma^2/\sigma_B^2.
To estimate both
\sigma_{\mathrm{n}}^2 and \tau^2 (and not only their sum), one needs
replicated measurements for some of the
\boldsymbol{s}_i.
We define
\boldsymbol{V}_{\alpha}
to be the (generalized) correlation matrix with elements
(\boldsymbol{V}_{\alpha})_{ij} =
\gamma_0 - \gamma(|\boldsymbol{A}\;(
\boldsymbol{s}_i-\boldsymbol{s}_j)|),
where the constant \gamma_0 is chosen large enough so that
\boldsymbol{V}_{\alpha}
is positive definite,
\gamma(\cdot) is a valid stationary or intrinsic variogram, and
\boldsymbol{A} =
\boldsymbol{A}(\alpha, f_1, f_2; \omega, \phi, \zeta)
is a matrix that is used to model geometrically anisotropic auto-correlation.
In more detail,
\boldsymbol{A}
maps an arbitrary point on an ellipsoidal surface with constant semi-variance in
\mathrm{I}\!\mathrm{R}^3,
centred on the origin, and having lengths of semi-principal axes,
\boldsymbol{p}_1,
\boldsymbol{p}_2,
\boldsymbol{p}_3,
equal to
|\boldsymbol{p}_1|=\alpha,
|\boldsymbol{p}_2|=f_1\,\alpha and
|\boldsymbol{p}_3|=f_2\,\alpha,
0 < f_2 \leq f_1 \leq 1,
respectively, onto the surface of the unit ball centred on the origin.
The orientation of the ellipsoid is defined by the three angles
\omega, \phi and \zeta:
\omegais the azimuth of
\boldsymbol{p}_1(= angle between north and the projection of\boldsymbol{p}_1onto thex-y-plane, measured from north to south positive clockwise in degrees),\phiis 90 degrees minus the latitude of
\boldsymbol{p}_1(= angle between the zenith and\boldsymbol{p}_1, measured from zenith to nadir positive clockwise in degrees), and\zetais the angle between
\boldsymbol{p}_2and the direction of the line, sayy^\prime, defined by the intersection between thex-y-plane and the plane orthogonal to\boldsymbol{p}_1running through the origin (\zetais measured fromy^\primepositive counter-clockwise in degrees).
The transformation matrix is given by
\boldsymbol{A}=
\left(\begin{array}{ccc}
1/\alpha & 0 & 0\\
0 & 1/(f_1\,\alpha) & 0\\
0 & 0 & 1/(f_2\,\alpha) \\
\end{array}\right)
(
\boldsymbol{C}_1,
\boldsymbol{C}_2,
\boldsymbol{C}_3,
)
where
\boldsymbol{C}_1^\mathrm{T} = ( \sin\omega \sin\phi, -\cos\omega \cos\zeta - \sin\omega \cos\phi \sin\zeta, \cos\omega \sin\zeta - \sin\omega \cos\phi \cos\zeta )
\boldsymbol{C}_2^\mathrm{T} = ( \cos\omega \sin\phi, \sin\omega \cos\zeta - \cos\omega \cos\phi \sin\zeta, -\sin\omega \sin\zeta - \cos\omega \cos\phi\cos\zeta )
\boldsymbol{C}_3^\mathrm{T} = (\cos\phi, \sin\phi \sin\zeta, \sin\phi \cos\zeta )
To model geometrically anisotropic variograms in
\mathrm{I}\!\mathrm{R}^2
one has to set \phi=90 and f_2 = 1,
and for f_1 = f_2 = 1
one obtains the model for isotropic auto-correlation
with range parameter \alpha.
Note that for isotropic auto-correlation the software processes data for
which d may exceed 3.
Two remarks are in order:
Clearly, the (generalized) covariance matrix of the observations
\boldsymbol{Y}is given by\mathrm{Cov}[\boldsymbol{Y},\boldsymbol{Y}^\mathrm{T}] = \tau^2 \boldsymbol{I} + \boldsymbol{\Gamma}_\theta.Depending on the context, the term “variogram parameters” denotes sometimes all parameters of a geometrically anisotropic variogram model, but in places only the parameters of an isotropic variogram model, i.e.
\sigma^2, \ldots, \alpha, \ldotsandf_1, \ldots, \zetaare denoted by the term “anisotropy parameters”. In the sequel\boldsymbol{\theta}is used to denote all variogram and anisotropy parameters except the nugget effect\tau^2.
Estimation
The unobserved spatial random effects
\boldsymbol{B} at the data locations
\boldsymbol{s}_i
and the model parameters
\boldsymbol{\beta}, \tau^2 and
\boldsymbol{\theta}^\mathrm{T} =
(\sigma^2, \sigma_{\mathrm{n}}^2, \alpha, \ldots, f_{1}, f_{2},
\omega, \phi, \zeta)
are unknown and are estimated in georob either by Gaussian
(Harville, 1977) or robust (Künsch et al., 2011)
restricted maximum likelihood (REML) or
Gaussian maximum likelihood (ML). Here ...
denote further parameters of the variogram such as the smoothness parameter
of the Whittle-Matérn model.
In brief, the robust REML method is based on the insight that for
given \boldsymbol{\theta} and \tau^2 the
Kriging predictions (= BLUP) of
\boldsymbol{B} and the generalized least
squares (GLS = ML) estimates of
\boldsymbol{\beta} can be obtained
simultaneously by maximizing
- \sum_i
\left(
\frac{
y_i -
\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}
\boldsymbol{\beta} -
B(\boldsymbol{s}_i)
}{\tau}
\right)^2 -
\boldsymbol{B}^{\mathrm{T}}
\boldsymbol{\Gamma}^{-1}_\theta
\boldsymbol{B}
with respect to
\boldsymbol{B} and
\boldsymbol{\beta}, e.g.
Harville (1977).
Hence, the BLUP of \boldsymbol{B},
ML estimates of \boldsymbol{\beta},
\boldsymbol{\theta} and \tau^2
are obtained by maximizing
- \log(\det(
\tau^2 \boldsymbol{I} +
\boldsymbol{\Gamma}_\theta
)) -
\sum_i
\left(
\frac{
y_i -
\boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T}
\boldsymbol{\beta} -
B(\boldsymbol{s}_i)
}{\tau}
\right)^2 -
\boldsymbol{B}^{\mathrm{T}}
\boldsymbol{\Gamma}^{-1}_\theta
\boldsymbol{B}
jointly with respect to
\boldsymbol{B},
\boldsymbol{\beta},
\boldsymbol{\theta} and \tau^2
or by solving the respective estimating equations.
The estimating equations can then by robustified by
replacing the standardized errors, say
\varepsilon_i/\tau = ( y_i - \boldsymbol{x}(\boldsymbol{s}_i)^\mathrm{T} \boldsymbol{\beta} - B(\boldsymbol{s}_i) ) / \tau, by a bounded or re-descending\psi-function,\psi_c(\varepsilon_i/\tau), of them (e.g. Maronna et al, 2006, chap. 2) and byintroducing suitable bias correction terms for Fisher consistency at the Gaussian model,
see Künsch et al. (2011) for details.
The robustified estimating equations
are solved numerically by a combination of iterated re-weighted least squares
(IRWLS) to estimate \boldsymbol{B} and
\boldsymbol{\beta} for given
\boldsymbol{\theta} and \tau^2
and nonlinear root finding by the function
nleqslv of the R package nleqslv
to get \boldsymbol{\theta} and \tau^2.
The robustness of the procedure is controlled by the tuning parameter c
of the \psi_c-function. For c \ge 1000 the algorithm computes
Gaussian (RE)ML estimates and customary plug-in Kriging predictions.
Instead of solving the Gaussian (RE)ML estimating equations, our software then
maximizes the Gaussian (restricted) log-likelihood using nlminb or
optim.
georob uses variogram models that were provided formerly by the
now archived R package RandomFields and are now implemented in
the function gencorr of georob.
Currently, estimation of the parameters of the following models is
implemented:
"RMaskey", "RMbessel", "RMcauchy",
"RMcircular", "RMcubic", "RMdagum",
"RMdampedcos", "RMdewijsian", "RMexp" (default),
"RMfbm", "RMgauss",
"RMgencauchy",
"RMgenfbm", "RMgengneiting", "RMgneiting",
"RMlgd",
"RMmatern", "RMpenta", "RMqexp",
"RMspheric", "RMstable", "RMwave",
"RMwhittle".
For most variogram parameters, closed-form expressions of \partial
\gamma/ \partial \theta_i and \partial \gamma/
\partial \tau^2 are used in the computations.
However, for the parameter \nu of the models "RMbessel",
"RMmatern" and "RMwhittle" \partial \gamma/ \partial
\nu is evaluated numerically by the function
numericDeriv, and this results in an increase in
computing time when \nu is estimated.
Prediction
Customary and robust plug-in external drift point Kriging predictions
can be computed for an non-sampled location
\boldsymbol{s}_0
from the covariates
\boldsymbol{x}(\boldsymbol{s}_0),
the estimated parameters
\widehat{\boldsymbol{\beta}},
\widehat{\boldsymbol{\theta}}
and the predicted random effects
\widehat{\boldsymbol{B}}
by
\widehat{Y}(\boldsymbol{s}_0) = \widehat{Z}(\boldsymbol{s}_0) =
\boldsymbol{x}(\boldsymbol{s}_0)^\mathrm{T}
\widehat{\boldsymbol{\beta}} +
\boldsymbol{\gamma}^\mathrm{T}_{\widehat{\theta}}(\boldsymbol{s}_0)
\boldsymbol{\Gamma}^{-1}_{\widehat{\theta}}
\widehat{\boldsymbol{B}},
where
\boldsymbol{\Gamma}_{\widehat{\theta}}
is the estimated (generalized) covariance matrix of
\boldsymbol{B} and
\boldsymbol{\gamma}_{\widehat{\theta}}(\boldsymbol{s}_0)
is the vector with the estimated (generalized) covariances between
\boldsymbol{B} and
B(\boldsymbol{s}_0).
Kriging variances can be computed as well, based on approximated covariances of
\widehat{\boldsymbol{B}} and
\widehat{\boldsymbol{\beta}}
(see Künsch et al., 2011, and appendices of
Nussbaum et al., 2014, for details).
The package georob provides in addition software for computing customary and robust external drift block Kriging predictions. The required integrals of the generalized covariance function are computed by functions of the R package constrainedKriging.
Functionality
For the time being, the functionality of georob is limited to geostatistical analyses of single response variables. No software is currently available for customary and robust multivariate geostatistical analyses. georob offers functions for:
Robustly fitting a spatial linear model to data that are possibly contaminated by independent errors from a long-tailed distribution by robust REML (see functions
georob— which also fits such models efficiently by Gaussian (RE)ML —profilelogLikandcontrol.georob).Extracting estimated model components (see
residuals.georob,rstandard.georob,ranef.georob).Robustly estimating sample variograms and for fitting variogram model functions to them (see
sample.variogramandfit.variogram.model).Model building by forward and backward selection of covariates for the external drift (see
waldtest.georob,step.georob,add1.georob,drop1.georob,extractAIC.georob,
logLik.georob,deviance.georob). For a robust fit, the log-likelihood is not defined. The function then computes the (restricted) log-likelihood of an equivalent Gaussian model with heteroscedastic nugget (seedeviance.georobfor details).Assessing the goodness-of-fit and predictive power of the model by K-fold cross-validation (see
cv.georobandvalidate.predictions).Computing customary and robust external drift point and block Kriging predictions (see
predict.georob,control.predict.georob).Unbiased back-transformation of both point and block Kriging predictions of log-transformed data to the original scale of the measurements (see
lgnpp).Computing unconditional and conditional Gaussian simulations from a fitted spatial linear model (see
condsim).
Author(s)
Andreas Papritz papritz@retired.ethz.ch.
References
Harville, D. A. (1977) Maximum likelihood approaches to variance component estimation and to related problems, Journal of the American Statistical Association, 72, 320–340, doi:10.1080/01621459.1977.10480998.
Künsch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (in preparation) Robust Geostatistics.
Künsch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (2011) Robust estimation of the external drift and the variogram of spatial data. Proceedings of the ISI 58th World Statistics Congress of the International Statistical Institute. doi:10.3929/ethz-a-009900710
Maronna, R. A., Martin, R. D. and Yohai, V. J. (2006) Robust Statistics Theory and Methods, Wiley, Hoboken, doi:10.1002/0470010940.
Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2014) Estimating soil organic carbon stocks of Swiss forest soils by robust external-drift kriging. Geoscientific Model Development, 7, 1197–1210. doi:10.5194/gmd-7-1197-2014.
See Also
georob for (robust) fitting of spatial linear models;
georobObject for a description of the class georob;
profilelogLik for computing profiles of Gaussian likelihoods;
plot.georob for display of RE(ML) variogram estimates;
control.georob for controlling the behaviour of georob;
georobModelBuilding for stepwise building models of class georob;
cv.georob for assessing the goodness of a fit by georob;
georobMethods for further methods for the class georob;
predict.georob for computing robust Kriging predictions;
lgnpp for unbiased back-transformation of Kriging prediction
of log-transformed data;
georobSimulation for simulating realizations of a Gaussian process
from model fitted by georob; and finally
sample.variogram and fit.variogram.model
for robust estimation and modelling of sample variograms.