georob-package: The georob Package

Description

This is a summary of the features and functionality of georob, a package in R for robust geostatistical analyses.

Arguments

Details

georob is a package for robust analyses of geostatistical data. Such data, say $y_i=y(\mbox{\boldmath$s$\unboldmath}_i)$, are recorded at a set of locations, $\mbox{\boldmath$s$\unboldmath}_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(\mbox{\boldmath$s$\unboldmath}_i)$, $j=1,2, \ldots, p$.

Model

We use the following model for the data $y_i=y(\mbox{\boldmath$s$\unboldmath}_{i})$: $$Y(\mbox{\boldmath$s$\unboldmath}_i) = Z(\mbox{\boldmath$s$\unboldmath}_i) + \varepsilon = \mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T} \mbox{\boldmath$\beta$\unboldmath} + B(\mbox{\boldmath$s$\unboldmath}_i) + \varepsilon_i,$$ where $Z(\mbox{\boldmath$s$\unboldmath}_i)=\mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T} \mbox{\boldmath$\beta$\unboldmath} + B(\mbox{\boldmath$s$\unboldmath}_i)$ is the so-called signal, $\mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T} \mbox{\boldmath$\beta$\unboldmath}$ is the external drift, $\{B(\mbox{\boldmath$s$\unboldmath})\}$ 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 $$\mbox{\boldmath $Y = X \beta + B + \varepsilon$\unboldmath},$$ where $\mbox{\boldmath $X$\unboldmath}$ is the model matrix with the rows $\mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T}$.

The (generalized) covariance matrix of the vector of spatial Gaussian random effects $\mbox{\boldmath$B$\unboldmath}$ is denoted by $$\mathrm{E}[ \mbox{\boldmath$B$\unboldmath}\, \mbox{\boldmath$B$\unboldmath}^\mathrm{T}] = \mbox{\boldmath$\Gamma$\unboldmath}_\theta = \sigma_{\mathrm{n}}^2\mbox{\boldmath$I$\unboldmath} + \sigma^2\mbox{\boldmath$V$\unboldmath}_\alpha = \sigma^2_Z \, \mbox{\boldmath$V$\unboldmath}_{\alpha,\xi} = \sigma^2_Z \, ((1-\xi) \, \mbox{\boldmath$I$\unboldmath} + \xi\, \mbox{\boldmath$V$\unboldmath}_\alpha ) ,$$ where $\sigma_{\mathrm{n}}^2$ is the variance of seemingly uncorrelated micro-scale variation in $B(\mbox{\boldmath$s$\unboldmath})$ that cannot be resolved with the chosen sampling design, $\sigma^2$ is the variance of the captured auto-correlated variation in $B(\mbox{\boldmath$s$\unboldmath})$, $\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 $\mbox{\boldmath$s$\unboldmath}_i$.

We define $\mbox{\boldmath$V$\unboldmath}_{\alpha}$ to be the matrix with elements $$(\mbox{\boldmath$V$\unboldmath}_{\alpha})_{ij} = \gamma_0 - \gamma(|\mbox{\boldmath$A$\unboldmath}\;( \mbox{\boldmath$s$\unboldmath}_i-\mbox{\boldmath$s$\unboldmath}_j)|),$$ where the constant $\gamma_0$ is chosen large enough so that $\mbox{\boldmath$V$\unboldmath}_{\alpha}$ is positive definite, $\gamma(\cdot)$ is a valid stationary or intrinsic variogram, and $\mbox{\boldmath$A$\unboldmath} = \mbox{\boldmath$A$\unboldmath}(\alpha, f_1, f_2; \omega, \phi, \zeta)$ is a matrix that is used to model geometrically anisotropic auto-correlation. In more detail, $\mbox{\boldmath$A$\unboldmath}$ 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, $\mbox{\boldmath$p$\unboldmath}_1$, $\mbox{\boldmath$p$\unboldmath}_2$, $\mbox{\boldmath$p$\unboldmath}_3$, equal to $|\mbox{\boldmath$p$\unboldmath}_1|=\alpha$, $|\mbox{\boldmath$p$\unboldmath}_2|=f_1\,\alpha$ and $|\mbox{\boldmath$p$\unboldmath}_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$:

$\omega$

is the azimuth of $\mbox{\boldmath$p$\unboldmath}_1$ (= angle between north and the projection of $\mbox{\boldmath$p$\unboldmath}_1$ onto the $x$-$y$-plane, measured from north to south positive clockwise in degrees),

$\phi$

is 90 degrees minus the altitude of $\mbox{\boldmath$p$\unboldmath}_1$ (= angle between the zenith and $\mbox{\boldmath$p$\unboldmath}_1$, measured from zenith to nadir positive clockwise in degrees), and

$\zeta$

is the angle between $\mbox{\boldmath$p$\unboldmath}_2$ and the direction of the line, say $y^\prime$, defined by the intersection between the $x$-$y$-plane and the plane orthogonal to $\mbox{\boldmath$p$\unboldmath}_1$ running through the origin ($\zeta$ is measured from $y^\prime$ positive counter-clockwise in degrees).

The transformation matrix is given by $$\mbox{\boldmath$A$\unboldmath}= \left(\begin{array}{ccc} 1/\alpha & 0 & 0\\ 0 & 1/(f_1\,\alpha) & 0\\ 0 & 0 & 1/(f_2\,\alpha) \\ \end{array}\right) ( \mbox{\boldmath$C$\unboldmath}_1, \mbox{\boldmath$C$\unboldmath}_2, \mbox{\boldmath$C$\unboldmath}_3, ) $$ where $$\mbox{\boldmath$C$\unboldmath}_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 ) $$ $$\mbox{\boldmath$C$\unboldmath}_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 ) $$ $$\mbox{\boldmath$C$\unboldmath}_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 $\mbox{\boldmath $Y$\unboldmath}$ is given by $$\mathrm{Cov}[\mbox{\boldmath $Y$\unboldmath},\mbox{\boldmath $Y$\unboldmath}^\mathrm{T}] = \tau^2 \mbox{\boldmath$I$\unboldmath} + \mbox{\boldmath$\Gamma$\unboldmath}_\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, \ldots$ and $f_1, \ldots, \zeta$ are denoted by the term “anisotropy parameters”. In the sequel $\mbox{\boldmath$\theta$\unboldmath}$ is used to denote all variogram and anisotropy parameters except the nugget effect $\tau^2$.

Estimation

The unobserved spatial random effects $\mbox{\boldmath$B$\unboldmath}$ at the data locations $\mbox{\boldmath$s$\unboldmath}_i$ and the model parameters $\mbox{\boldmath$\beta$\unboldmath}$, $\tau^2$ and $\mbox{\boldmath$\theta$\unboldmath}^\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 or robust 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<U+008E>rn model.

In brief, the robust REML method is based on the insight that for given $\mbox{\boldmath$\theta$\unboldmath}$ and $\tau^2$ the Kriging predictions (= BLUP) of $\mbox{\boldmath$B$\unboldmath}$ and the generalized least squares (GLS = ML) estimates of $\mbox{\boldmath$\beta$\unboldmath}$ can be obtained simultaneously by maximizing $$ - \sum_i \left( \frac{ y_i - \mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T} \mbox{\boldmath$\beta$\unboldmath} - B(\mbox{\boldmath$s$\unboldmath}_i) }{\tau} \right)^2 - \mbox{\boldmath$B$\unboldmath}^{\mathrm{T}} \mbox{\boldmath$\Gamma$\unboldmath}^{-1}_\theta \mbox{\boldmath$B$\unboldmath} $$ with respect to $\mbox{\boldmath$B$\unboldmath}$ and $\mbox{\boldmath$\beta$\unboldmath}$ e.g. Harville (1977). Hence, the BLUP of $\mbox{\boldmath$B$\unboldmath}$, ML estimates of $\mbox{\boldmath$\beta$\unboldmath}$, $\mbox{\boldmath$\theta$\unboldmath}$ and $\tau^2$ are obtained by maximizing $$ - \log(\det( \tau^2 \mbox{\boldmath$I$\unboldmath} + \mbox{\boldmath$\Gamma$\unboldmath}_\theta )) - \sum_i \left( \frac{ y_i - \mbox{\boldmath $x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_i)^\mathrm{T} \mbox{\boldmath$\beta$\unboldmath} - B(\mbox{\boldmath$s$\unboldmath}_i) }{\tau} \right)^2 - \mbox{\boldmath$B$\unboldmath}^{\mathrm{T}} \mbox{\boldmath$\Gamma$\unboldmath}^{-1}_\theta \mbox{\boldmath$B$\unboldmath} $$ jointly with respect to $\mbox{\boldmath$B$\unboldmath}$, $\mbox{\boldmath$\beta$\unboldmath}$, $\mbox{\boldmath$\theta$\unboldmath}$ and $\tau^2$ or by solving the respective estimating equations.

The estimating equations can then by robustified by

replacing the standardized residuals, say $\varepsilon_i/\tau$, by a bounded or re-descending $\psi$-function, $\psi_c(\varepsilon_i/\tau)$, of them (e.g. Marona et al, 2006, chap. 2) and by
introducing suitable bias correction terms for Fisher consistency at the Gaussian model,

see K<U+009F>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 $\mbox{\boldmath$B$\unboldmath}$ and $\mbox{\boldmath$\beta$\unboldmath}$ for given $\mbox{\boldmath$\theta$\unboldmath}$ and $\tau^2$ and nonlinear root finding by the function nleqslv of the R package nleqslv to get $\mbox{\boldmath$\theta$\unboldmath}$ 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 implemented in the R package RandomFields (see RMmodel). 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$ 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

Robust plug-in external drift point Kriging predictions can be computed for an non-sampled location $\mbox{\boldmath$s$\unboldmath}_0$ from the covariates $\mbox{\boldmath$x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_0)$, the estimated parameters $\widehat{\mbox{\boldmath$\beta$\unboldmath}}$, $\widehat{\mbox{\boldmath$\theta$\unboldmath}}$ and the predicted random effects $\widehat{\mbox{\boldmath$B$\unboldmath}}$ by $$ \widehat{Y}(\mbox{\boldmath$s$\unboldmath}_0) = \widehat{Z}(\mbox{\boldmath$s$\unboldmath}_0) = \mbox{\boldmath$x$\unboldmath}(\mbox{\boldmath$s$\unboldmath}_0)^\mathrm{T} \widehat{\mbox{\boldmath$\beta$\unboldmath}} + \mbox{\boldmath$\gamma$\unboldmath}^\mathrm{T}_{\widehat{\theta}}(\mbox{\boldmath$s$\unboldmath}_0) \mbox{\boldmath$\Gamma$\unboldmath}^{-1}_{\widehat{\theta}} \widehat{\mbox{\boldmath$B$\unboldmath}}, $$

where $\mbox{\boldmath$\Gamma$\unboldmath}_{\widehat{\theta}}$ is the estimated (generalized) covariance matrix of $\mbox{\boldmath$B$\unboldmath}$ and $\mbox{\boldmath$\gamma$\unboldmath}_{\widehat{\theta}}(\mbox{\boldmath$s$\unboldmath}_0)$ is the vector with the estimated (generalized) covariances between $\mbox{\boldmath$B$\unboldmath}$ and $B(\mbox{\boldmath$s$\unboldmath}_0)$. Kriging variances can be computed as well, based on approximated covariances of $\widehat{\mbox{\boldmath$B$\unboldmath}}$ and $\widehat{\mbox{\boldmath$\beta$\unboldmath}}$ (see K<U+009F>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.

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 functions georob --- which also fits such models efficiently by Gaussian (RE)ML --- profilelogLik and control.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.variogram and fit.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 (see deviance.georob for details).
Assessing the goodness-of-fit and predictive power of the model by K-fold cross-validation (see cv.georob and validate.predictions).
Computing 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).

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://dx.doi.org/10.3929/ethz-a-007555133

K<U+009F>nsch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (in preparation) Robust Geostatistics.

K<U+009F>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. http://e-collection.library.ethz.ch/eserv/eth:7080/eth-7080-01.pdf

Maronna, R. A., Martin, R. D. and Yohai, V. J. (2006) Robust Statistics Theory and Methods, John Wiley \& Sons.