georob-package: The georob Package
In georob: Robust Geostatistical Analysis of Spatial Data
Description
Details
Author(s)
References
See Also
Description
This is a summary of the features of georob,
a package in R for robust geostatistical analyses.
Details
georob is a package for robust analyses of geostatistical data.
Such data, say y(s_i), are
recorded at a set of locations,
s_i, i=1,2, …, n, in a
domain G in R^d, d in (1,2,3), along
with covariate information
x_j(s_i), j=1,2,
…, p.
Model
We use the following model for the data
Y(s) = S(s) + ε(s) = x(s)^T β + B(s) + ε(s),
where
S(s) = x(s)^T β + B(s) is the so-called signal,
x^Tβ
is the external drift,
B(s) is a stationary or
intrinsic spatial Gaussian random field with zero mean, and
ε(s) are
i.i.d errors from a possibly long-tailed distribution with scale parameter
τ (τ^2 is usually called nugget effect).
In vector form the model leads to
Y = X β + B + ε,
where X is the model matrix consisting
of the rows
x^T(s_i).
The (generalized) covariance matrix of the vector of unobserved
spatial Gaussian random effects
B
is denoted by
E[B B^T] = σ_n^2 I+σ^2 V_α,
where σ_n^2 is the variance of seemingly
uncorrelated micro-scale variation in
B(s)
that cannot be resolved with the chosen sampling design, and
σ^2 is the variance of the captured auto-correlated variation in
B(s). To estimate both
σ_n^2 and τ^2 (and not only their sum), one needs
replicated measurements for some of the
s_i.
We define
V_α
to be the matrix with elements
(V_α)_ij = γ_0 - γ(|A (s_i - s_j )|),
where the constant γ_0 is chosen large enough so that
V_α
is positive definite,
v() is a valid stationary or intrinsic variogram, and
A = A(α, f_1, f_2; ω, φ, ζ)
is a matrix that is used to model geometrically anisotropic auto-correlation.
In more detail,
A
maps an arbitrary point on an ellipsoidal surface with constant semivariance in
R^d,
centred on the origin, and having lengths of semi-principal axes,
P_1,
P_2,
P_3,
equal to
|P_1|=α,
|P_2|=f_1 α and
|P_3|=f_2 α,
0 < f_2 <= f_1 <= 1,
respectively, onto the surface of the unit ball centred on the origin.
The orientation of the ellipsoid is defined by the three angles
ω, φ and ζ:
- ω
is the azimuth of P_1
(= angle between north and the projection
of
P_1
onto the x-y-plane,
measured from north to south positive clockwise in degrees),
- φ
is 90 degrees minus the altitude of
P_1
(= angle between the zenith and
P_1,
measured from zenith to nadir positive clockwise in degrees), and
- ζ
is the angle between
P_2
and the direction of the line, say y',
defined by the intersection between the
x-y-plane and the plane orthogonal to
P_1 running through the origin
(ζ is measured from y' positive counterclockwise in degrees).
The transformation matrix is given by
A=diag(1/α, 1/(f_1\,α),1/(f_2\,α)) (C_1, C_2, C_3)
where
C_1^T=(sinφ sinω, -cosζ cosω + sinζ cosφ sinω, -cosω sinζ - cosζ cosφ sinω)
C_2^T=(-sinφ cosω, cosω sinζ cosφ + cosζ sinω, -cosζ cosω cosφ + sinζ sinω)
C_3^T=(cosφ, -sinφ sinζ, cosζ sinφ)
To model geometrically anisotropic variograms in
R^2
one has to set φ=90 and f_2 = 1,
and for f_1 = f_2 = 1
one obtains the model for isotropic auto-correlation
with range parameter α.
Note that for isotropic auto-correlation d may exceed 3.
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. σ^2, …, α, … and
f_1, …, ζ are denoted by the term “anisotropy
parameters”.
Estimation
The spatial random effects
B
and the model parameters
β and
θ^T =
(σ^2, σ^2_n, τ^2, α, f_1, f_2,
ω, φ, ζ, …)
are estimated by robust Restricted Maximum Likelihood (REML). 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
the kriging predictions of B
and the Gaussian Maximum Likelihood estimates of
β and
θ can be obtained
simultaneously by maximizing
- log(det( τ^2 I + Γ_θ) )
- ∑_i ( y_i - x^T(s_i)β - B(s_i) )^2 / τ^2
- B^T Γ_θ^-1 B
with respect to
B,
β,
θ.
The respective estimating equations can then by robustified by
replacing the standardized residuals, say
r/s,
by a bounded function,
ψ_c(r/s),
of them and by
introducting 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 B and
β for given
θ
and nonlinear root finding by the function
nleqslv of the R package nleqslv
to get θ.
The robustness of the procedure is controlled by the tuning parameter c
of the ψ_c-function. For c>=1000 the algorithm computes
Gaussian REML estimates and customary plug-in kriging predictions.
Instead of solving the Gaussian REML estimating equations, our software then
maximizes the Gaussian restricted loglikelihood using
optim.
georob uses variogram models implemented in the R package
RandomFields (see Variogram). Currently,
estimation of the parameters of the following models is implemented:
"bessel", "cauchy", "cauchytbm", "circular",
"cubic", "dagum", "dampedcosine", "DeWijsian",
"exponential", "fractalB", "gauss", "genB",
"gencauchy", "gengneiting", "gneiting", "lgd1",
"matern", "penta", "power", "qexponential",
"spherical", "stable", "wave", "whittle".
For most variogram parameters, closed-form expressions of
dγ/dθ_i
are used in the computations.
However, for the parameter ν of the models
"bessel", "matern" and "whittle"
dγ/dθ_i
is evaluated numerically by the function
numericDeriv, and this results in an
increase in computing time when ν is estimated.
Prediction
Robust plug-in external drift point kriging predictions
can be computed for an unsampled location
s_0
from the covariates
x(s_0),
the estimated parameters
hatβ,
hatθ
and the predicted random effects
hatB
by
hatY(s_0) = x^T(s_0) hatβ + γ_^Thatθ(s_0) Γ_hatθ^-1 hatB,
where
Γ_hatθ
is the estimated (generalized) covariance matrix of
B and
gamma_hatθ(s_0)
is the vector with the estimated (generalized) covariances between
B and
B(s_0).
Kriging variances can be computed as well, based on approximated covariances of
hatB and
hatβ
(see Künsch et al., 2011, and Appendices of
Nussbaum et al., 2012, for details).
The package georob provides in addition software for computing
robust external drift block kriging predictions. The required
integrals of the generalized covariance function are computed by
functions of the R package constrainedKriging.
Main functionality
For the time being, the functionality of georob is limited
to robust geostatistical analyses of single response variables.
No software is currently available for 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 georob, which also
fits such models efficiently by Gaussian REML).
Assessing the goodness-of-fit of the model by
K-fold cross-validation (see cv.georob).
Computing robust external drift point and block kriging
predictions (see 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).
Robustly estimating sample variograms and for fitting
variogram model functions to them (see
sample.variogram and fit.variogram.model).
Author(s)
Andreas Papritz andreas.papritz@env.ethz.ch
http://www.step.ethz.ch/people/scientific-staff/andreas-papritz
References
Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2012)
Organic Carbon Stocks of Swiss Forest Soils,
Institute of Terrestrial Ecosystems, ETH Zurich and
Swiss Federal Institute for Forest, Snow and Landscape Research
(WSL), pp. 51.
http://e-collection.library.ethz.ch/eserv/eth:6027/eth-6027-01.pdf
Kuensch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (in
preparation) Robust Geostatistics.
Kuensch, 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.
http://e-collection.library.ethz.ch/eserv/eth:7080/eth-7080-01.pdf
See Also
georob for (robust) fitting of spatial linear models;
georobObject for a description of the class georob;
plot.georob for display of REML variogram estimates;
georob.control for controlling the behaviour of georob;
cv.georob for assessing the goodness of a fit by georob;
predict.georob for computing robust kriging predictions; and finally
georobModelBuilding for stepwise building models of class georob;
georobMethods for further methods for the class georob,
sample.variogram and fit.variogram.model
for robust estimation and modelling of sample variograms.
georob documentation built on May 2, 2019, 6:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Details Author(s) References See Also
Description
This is a summary of the features of georob, a package in R for robust geostatistical analyses.
Details
georob is a package for robust analyses of geostatistical data. Such data, say y(s_i), are recorded at a set of locations, s_i, i=1,2, …, n, in a domain G in R^d, d in (1,2,3), along with covariate information x_j(s_i), j=1,2, …, p.
Model
We use the following model for the data
Y(s) = S(s) + ε(s) = x(s)^T β + B(s) + ε(s),
where S(s) = x(s)^T β + B(s) is the so-called signal, x^Tβ is the external drift, B(s) is a stationary or intrinsic spatial Gaussian random field with zero mean, and ε(s) are i.i.d errors from a possibly long-tailed distribution with scale parameter τ (τ^2 is usually called nugget effect). In vector form the model leads to
Y = X β + B + ε,
where X is the model matrix consisting of the rows x^T(s_i).
The (generalized) covariance matrix of the vector of unobserved spatial Gaussian random effects B is denoted by
E[B B^T] = σ_n^2 I+σ^2 V_α,
where σ_n^2 is the variance of seemingly uncorrelated micro-scale variation in B(s) that cannot be resolved with the chosen sampling design, and σ^2 is the variance of the captured auto-correlated variation in B(s). To estimate both σ_n^2 and τ^2 (and not only their sum), one needs replicated measurements for some of the s_i.
We define V_α to be the matrix with elements
(V_α)_ij = γ_0 - γ(|A (s_i - s_j )|),
where the constant γ_0 is chosen large enough so that V_α is positive definite, v() is a valid stationary or intrinsic variogram, and A = A(α, f_1, f_2; ω, φ, ζ) is a matrix that is used to model geometrically anisotropic auto-correlation. In more detail, A maps an arbitrary point on an ellipsoidal surface with constant semivariance in R^d, centred on the origin, and having lengths of semi-principal axes, P_1, P_2, P_3, equal to |P_1|=α, |P_2|=f_1 α and |P_3|=f_2 α, 0 < f_2 <= f_1 <= 1, respectively, onto the surface of the unit ball centred on the origin.
The orientation of the ellipsoid is defined by the three angles ω, φ and ζ:
- ω
is the azimuth of P_1 (= angle between north and the projection of P_1 onto the x-y-plane, measured from north to south positive clockwise in degrees),
- φ
is 90 degrees minus the altitude of P_1 (= angle between the zenith and P_1, measured from zenith to nadir positive clockwise in degrees), and
- ζ
is the angle between P_2 and the direction of the line, say y', defined by the intersection between the x-y-plane and the plane orthogonal to P_1 running through the origin (ζ is measured from y' positive counterclockwise in degrees).
The transformation matrix is given by
A=diag(1/α, 1/(f_1\,α),1/(f_2\,α)) (C_1, C_2, C_3)
where
C_1^T=(sinφ sinω, -cosζ cosω + sinζ cosφ sinω, -cosω sinζ - cosζ cosφ sinω)
C_2^T=(-sinφ cosω, cosω sinζ cosφ + cosζ sinω, -cosζ cosω cosφ + sinζ sinω)
C_3^T=(cosφ, -sinφ sinζ, cosζ sinφ)
To model geometrically anisotropic variograms in R^2 one has to set φ=90 and f_2 = 1, and for f_1 = f_2 = 1 one obtains the model for isotropic auto-correlation with range parameter α. Note that for isotropic auto-correlation d may exceed 3.
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. σ^2, …, α, … and f_1, …, ζ are denoted by the term “anisotropy parameters”.
Estimation
The spatial random effects B and the model parameters β and θ^T = (σ^2, σ^2_n, τ^2, α, f_1, f_2, ω, φ, ζ, …) are estimated by robust Restricted Maximum Likelihood (REML). 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 the kriging predictions of B and the Gaussian Maximum Likelihood estimates of β and θ can be obtained simultaneously by maximizing
- log(det( τ^2 I + Γ_θ) ) - ∑_i ( y_i - x^T(s_i)β - B(s_i) )^2 / τ^2 - B^T Γ_θ^-1 B
with respect to B, β, θ. The respective estimating equations can then by robustified by
replacing the standardized residuals, say r/s, by a bounded function, ψ_c(r/s), of them and by
introducting 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 B and
β for given
θ
and nonlinear root finding by the function
nleqslv of the R package nleqslv
to get θ.
The robustness of the procedure is controlled by the tuning parameter c
of the ψ_c-function. For c>=1000 the algorithm computes
Gaussian REML estimates and customary plug-in kriging predictions.
Instead of solving the Gaussian REML estimating equations, our software then
maximizes the Gaussian restricted loglikelihood using
optim.
georob uses variogram models implemented in the R package
RandomFields (see Variogram). Currently,
estimation of the parameters of the following models is implemented:
"bessel", "cauchy", "cauchytbm", "circular",
"cubic", "dagum", "dampedcosine", "DeWijsian",
"exponential", "fractalB", "gauss", "genB",
"gencauchy", "gengneiting", "gneiting", "lgd1",
"matern", "penta", "power", "qexponential",
"spherical", "stable", "wave", "whittle".
For most variogram parameters, closed-form expressions of
dγ/dθ_i
are used in the computations.
However, for the parameter ν of the models
"bessel", "matern" and "whittle"
dγ/dθ_i
is evaluated numerically by the function
numericDeriv, and this results in an
increase in computing time when ν is estimated.
Prediction
Robust plug-in external drift point kriging predictions can be computed for an unsampled location s_0 from the covariates x(s_0), the estimated parameters hatβ, hatθ and the predicted random effects hatB by
hatY(s_0) = x^T(s_0) hatβ + γ_^Thatθ(s_0) Γ_hatθ^-1 hatB,
where Γ_hatθ is the estimated (generalized) covariance matrix of B and gamma_hatθ(s_0) is the vector with the estimated (generalized) covariances between B and B(s_0). Kriging variances can be computed as well, based on approximated covariances of hatB and hatβ (see Künsch et al., 2011, and Appendices of Nussbaum et al., 2012, for details).
The package georob provides in addition software for computing robust external drift block kriging predictions. The required integrals of the generalized covariance function are computed by functions of the R package constrainedKriging.
Main functionality
For the time being, the functionality of georob is limited to robust geostatistical analyses of single response variables. No software is currently available for 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
georob, which also fits such models efficiently by Gaussian REML).Assessing the goodness-of-fit of the model by K-fold cross-validation (see
cv.georob).Computing robust external drift point and block kriging predictions (see
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).Robustly estimating sample variograms and for fitting variogram model functions to them (see
sample.variogramandfit.variogram.model).
Author(s)
Andreas Papritz andreas.papritz@env.ethz.ch
http://www.step.ethz.ch/people/scientific-staff/andreas-papritz
References
Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2012) Organic Carbon Stocks of Swiss Forest Soils, Institute of Terrestrial Ecosystems, ETH Zurich and Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), pp. 51. http://e-collection.library.ethz.ch/eserv/eth:6027/eth-6027-01.pdf
Kuensch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (in preparation) Robust Geostatistics.
Kuensch, 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. http://e-collection.library.ethz.ch/eserv/eth:7080/eth-7080-01.pdf
See Also
georob for (robust) fitting of spatial linear models;
georobObject for a description of the class georob;
plot.georob for display of REML variogram estimates;
georob.control for controlling the behaviour of georob;
cv.georob for assessing the goodness of a fit by georob;
predict.georob for computing robust kriging predictions; and finally
georobModelBuilding for stepwise building models of class georob;
georobMethods for further methods for the class georob,
sample.variogram and fit.variogram.model
for robust estimation and modelling of sample variograms.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.