Package 'georob'

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$:

$\omega$

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

$\phi$

is 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

$\zeta$

is the angle between $\boldsymbol{p}_2$ and the direction of the line, say $y^\prime$, defined by the intersection between the $x$-$y$-plane and the plane orthogonal to $\boldsymbol{p}_1$ running through the origin ($\zeta$ is measured from $y^\prime$ positive 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:

1. 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.$

2. 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 $\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 by

• introducing 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:

1. 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).

2. Extracting estimated model components (see residuals.georob, rstandard.georob, ranef.georob).

3. Robustly estimating sample variograms and for fitting variogram model functions to them (see sample.variogram and fit.variogram.model).

4. 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).

5. Assessing the goodness-of-fit and predictive power of the model by K-fold cross-validation (see cv.georob and validate.predictions).

6. Computing customary and robust external drift point and block Kriging predictions (see predict.georob, control.predict.georob).

7. Unbiased back-transformation of both point and block Kriging predictions of log-transformed data to the original scale of the measurements (see lgnpp).

8. Computing unconditional and conditional Gaussian simulations from a fitted spatial linear model (see condsim).

Author(s)

Andreas Papritz [email protected].

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.

---

Compact Storage of Symmetric and Triangular Matrices

Description

The utility function compress stores symmetric or triangular matrices compactly by retaining only the diagonal and either the lower or upper off-diagonal elements. The function expand restores such compressed matrices again to a square form.

Usage

compress(m)

  expand(object)

Arguments

m	either a single symmetric, lower or upper triangular matrix or a list of such matrices. The type of m (or of its component matrices) must be defined by the attribute struc with possible values "sym" (symmetric), "lt" (lower triangular) or "ut" (upper triangular).
object	a single compressed matrix or a list of such matrices generated by compress, see Value. The type of object (or of its components) must be defined by the attribute struc with possible values "sym" (symmetric), "lt" (lower triangular) or "ut" (upper triangular).

Value

If m is a single square matrix then compress generates a compressed matrix, which is a list with two components:

diag	a vector with the diagonal elements of m.
tri	a vector with off-diagonal elements.

If m is a list of square matrices then the result is also a list of compressed matrices.

expand creates a square matrix if object is a list with components diag and tri and a list of square matrices if object is a list of such lists. If m or objects are lists that contain other components than square or compressed matrices then these additional components are returned unchanged.

Author(s)

Andreas Papritz [email protected].

See Also

georob for (robust) fitting of spatial linear models.

Examples

data(meuse)

r.logzn.rob <- georob(log(zinc) ~ sqrt(dist) + ffreq, data = meuse,
    locations = ~ x + y, variogram.model = "RMexp",
    param = c(variance = 0.15, nugget = 0.05, scale = 200),
    tuning.psi = 1)

cov2cor(expand(r.logzn.rob[["cov"]][["cov.betahat"]]))

---

Control Parameters for georob

Description

This page documents parameters used to control georob. It describes the arguments of the functions control.georob, param.transf, fwd.transf, dfwd.transf, bwd.transf, control.rq, control.nleqslv, control.nlminb and control.optim, which all serve to control the behaviour of georob.

Usage

control.georob(ml.method = c("REML", "ML"), reparam = TRUE,
    maximizer = c("nlminb", "optim"), initial.param = TRUE,
    initial.fixef = c("lmrob", "rq", "lm"), bhat = NULL,
    min.rweight = 0.25,
    param.tf = param.transf(), fwd.tf = fwd.transf(),
    deriv.fwd.tf = dfwd.transf(), bwd.tf = bwd.transf(),
    psi.func = c("logistic", "t.dist", "huber"),
    irwls.maxiter = 50,
    irwls.ftol = 1.e-5, force.gradient = FALSE,
    min.condnum = 1.e-12, zero.dist = sqrt(.Machine[["double.eps"]]),
    error.family.estimation = c("gaussian", "long.tailed"),
    error.family.cov.effects = c("gaussian", "long.tailed"),
    error.family.cov.residuals = c("gaussian", "long.tailed"),
    cov.bhat = TRUE, full.cov.bhat = FALSE, cov.betahat = TRUE,
    cov.delta.bhat = TRUE, full.cov.delta.bhat = TRUE,
    cov.delta.bhat.betahat = TRUE,
    cov.ehat = TRUE, full.cov.ehat = FALSE,
    cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
    hessian = TRUE,
    rq = control.rq(), lmrob = lmrob.control(),
    nleqslv = control.nleqslv(),
    optim = control.optim(), nlminb = control.nlminb(),
    pcmp = control.pcmp(), ...)

param.transf(variance = "log", snugget = "log", nugget = "log", scale = "log",
    alpha = c(
      RMaskey = "log", RMdewijsian = "logit2", RMfbm = "logit2", RMgencauchy = "logit2",
      RMgenfbm = "logit2", RMlgd = "identity", RMqexp = "logit1", RMstable = "logit2"
    ),
    beta = c(RMdagum = "logit1", RMgencauchy = "log", RMlgd = "log"),
    delta = "logit1", gamma = c(RMcauchy = "log", RMdagum = "logit1"),
    kappa = "logit3", lambda = "log", mu = "log", nu = "log",
    f1 = "log", f2  ="log", omega = "identity", phi = "identity", zeta = "identity")

fwd.transf(...)

dfwd.transf(...)

bwd.transf(...)

control.rq(tau = 0.5, rq.method = c("br", "fnb", "pfn"),
    rq.alpha = 0.1, ci = FALSE, iid = TRUE,
    interp = TRUE, tcrit = TRUE, rq.beta = 0.99995, eps = 1e-06,
    Mm.factor = 0.8, max.bad.fixup = 3, ...)

control.nleqslv(method = c("Broyden", "Newton"),
    global = c("dbldog", "pwldog", "qline", "gline", "none"),
    xscalm = c("fixed", "auto"), control = list(ftol = 1e-04), ...)

control.optim(method = c("BFGS", "Nelder-Mead", "CG",
        "L-BFGS-B", "SANN", "Brent"), lower = -Inf, upper = Inf,
    control = list(reltol = 1e-05), ...)

control.nlminb(control = list(rel.tol = 1.e-5), lower = -Inf,
    upper = Inf, ...)

Arguments

ml.method	a character keyword defining whether Gaussian maximum likelihood (ML) or restricted maximum likelihood (REML default) estimates will be computed (ignored if tuning.psi <= tuning.psi.nr).
reparam	a logical scalar. If TRUE (default) the re-parametrized variance parameters $\sigma_B^2$, $\eta$ and $\xi$ are estimated by Gaussian (RE)ML, otherwise the original parameters $\tau^2$, $\sigma_{\mathrm{n}}^2$ and $\sigma^2$ (cf. subsection Estimating variance parameters by Gaussian (RE)ML, section Details of georob).
maximizer	a character keyword defining whether the Gaussian (restricted) log-likelihood is maximized by nlminb (default) or optim.
initial.param	a logical scalar, controlling whether initial values of variogram parameters are computed for solving the robustified estimating equations of the variogram and anisotropy parameters. If initial.param = TRUE (default) robust initial values of parameters are computed by discarding outlying observations based on the “robustness weights” of the initial fit of the regression model by lmrob and fitting the spatial linear model by Gaussian REML to the pruned data set. For initial.param = FALSE no initial parameter values are computed and the estimating equations are solved with the initial values passed by param and aniso to georob (see Details of georob).
initial.fixef	a character keyword defining whether the function lmrob or rq is used to compute robust initial estimates of the regression parameters $\boldsymbol{\beta}$ (default "lmrob"). If the fixed effects model matrix has not full columns rank, then lm is used to compute initial values of the regression coefficients.
bhat	a numeric vector with initial values for the spatial random effects $\widehat{\boldsymbol{B}}$, with $\widehat{\boldsymbol{B}}=\boldsymbol{0}$ if bhat is equal to NULL (default).
min.rweight	a positive numeric with the “robustness weight” of the initial lmrob fit that observations must exceed to be used for computing robust initial estimates of variogram parameters by setting initial.param = TRUE (see georob; default 0.25).
param.tf	a function such as param.transf, which returns a named list of character strings that define the transformations to be applied to the variogram parameters for model fitting, see Details.
fwd.tf	a function such as fwd.transf, which returns a named list of invertible functions to be used to transform variogram parameters, see Details.
deriv.fwd.tf	a function such as dfwd.transf, which returns a named list of functions corresponding to the first derivatives of fwd.tf, see Details.
bwd.tf	a function such as bwd.transf, which returns the named list of inverse functions corresponding to fwd.tf, see Details.
psi.func	a character keyword defining what $\psi_c$-function should be used for robust model fitting. Possible values are "logistic" (a scaled and shifted logistic CDF, default), "t.dist" (re-descending $\psi_c$-function associated with Student $t$-distribution with $c$ degrees of freedom) and "huber" (Huber's $\psi_c$-function).
irwls.maxiter	a positive integer equal to the maximum number of IRWLS iterations to solve the estimating equations of $\boldsymbol{B}$ and $\boldsymbol{\beta}$ (default 50).
irwls.ftol	a positive numeric with the convergence criterion for IRWLS. Convergence is assumed if the objective function change of a IRWLS iteration does not exceed ftol.
force.gradient	a logical scalar controlling whether the estimating equations or the gradient of the Gaussian restricted log-likelihood are evaluated even if all variogram parameters are fixed (default FALSE).
min.condnum	a positive numeric with the minimum acceptable ratio of smallest to largest singular value of the model matrix $\boldsymbol{X}$ (default 1.e-12).
zero.dist	a positive numeric equal to the maximum distance, separating two sampling locations that are still considered as being coincident.
error.family.estimation	a character keyword, defining the probability distribution for $\varepsilon$ (default: "gaussian") that is used to approximate the covariance of $\widehat{\boldsymbol{B}}$ when solving the estimating equations, see Details.
error.family.cov.effects	a character keyword, defining the probability distribution for $\varepsilon$ (default: "gaussian") that is used to approximate the covariances of $\widehat{\boldsymbol{\beta}}$, $\widehat{\boldsymbol{B}}$ and $\boldsymbol{B}-\widehat{\boldsymbol{B}}$, see Details.
error.family.cov.residuals	a character keyword, defining the probability distribution for $\varepsilon$ (default: "long.tailed") that is used to approximate the covariances of $\widehat{\boldsymbol{\varepsilon}}=\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}} - \widehat{\boldsymbol{B}}$ and $\widehat{\boldsymbol{\varepsilon}} + \widehat{\boldsymbol{B}} =\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}}$, see Details.
cov.bhat	a logical scalar controlling whether the covariances of $\widehat{\boldsymbol{B}}$ are returned by georob (default FALSE).
full.cov.bhat	a logical scalar controlling whether the full covariance matrix (TRUE) or only the variance vector of $\widehat{\boldsymbol{B}}$ is returned (default FALSE).
cov.betahat	a logical scalar controlling whether the covariance matrix of $\widehat{\boldsymbol{\beta}}$ is returned (default TRUE).
cov.delta.bhat	a logical scalar controlling whether the covariances of $\boldsymbol{B}- \widehat{\boldsymbol{B}}$ are returned (default TRUE).
full.cov.delta.bhat	a logical scalar controlling whether the full covariance matrix (TRUE) or only the variance vector of $\boldsymbol{B}- \widehat{\boldsymbol{B}}$ is returned (default TRUE).
cov.delta.bhat.betahat	a logical scalar controlling whether the covariance matrix of $\boldsymbol{B}- \widehat{\boldsymbol{B}}$ and $\widehat{\boldsymbol{\beta}}$ is returned (default TRUE).
cov.ehat	a logical scalar controlling whether the covariances of $\widehat{\boldsymbol{\varepsilon}}=\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}} - \widehat{\boldsymbol{B}}$ are returned (default TRUE).
full.cov.ehat	a logical scalar controlling whether the full covariance matrix (TRUE) or only the variance vector of $\widehat{\boldsymbol{\varepsilon}}=\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}} - \widehat{\boldsymbol{B}}$ is returned (default FALSE).
cov.ehat.p.bhat	a logical scalar controlling whether the covariances of $\widehat{\boldsymbol{\varepsilon}} + \widehat{\boldsymbol{B}} =\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}}$ are returned (default FALSE).
full.cov.ehat.p.bhat	a logical scalar controlling whether the full covariance matrix (TRUE) or only the variance vector of $\widehat{\boldsymbol{\varepsilon}} + \widehat{\boldsymbol{B}} =\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}}$ is returned (default FALSE).
hessian	a logical scalar controlling whether for Gaussian (RE)ML the Hessian should be computed at the MLEs.
rq	a list of arguments passed to rq or a function such as control.rq that generates such a list (see rq for allowed arguments).
lmrob	a list of arguments passed to the control argument of lmrob or a function such as lmrob.control that generates such a list (see lmrob.control for allowed arguments).
nleqslv	a list of arguments passed to nleqslv or a function such as control.nleqslv that generates such a list (see nleqslv for allowed arguments).
nlminb	a list of arguments passed to nlminb or a function such as control.nlminb that generates such a list (see nlminb for allowed arguments).
optim	a list of arguments passed to optim or a function such as control.optim that generates such a list (see optim for allowed arguments).
pcmp	a list of arguments, passed e.g. to pmm or a function such as control.pcmp that generates such a list (see control.pcmp for allowed arguments).
...	for fwd.transf, dfwd.transf and bwd.transf a named vector of functions, extending the definition of transformations for variogram parameters (see Details).
variance, snugget, nugget, scale, alpha, beta, delta, gamma, kappa, lambda, mu, nu	character strings with names of transformation functions of the variogram parameters.
f1, f2, omega, phi, zeta	character strings with names of transformation functions of the anisotropy variogram parameters.
tau, rq.method, rq.alpha, ci, iid, interp, tcrit	arguments passed as ... to rq. Note that only "br", "fnb" and "pfn" methods of rq() are currently supported.
rq.beta, eps, Mm.factor, max.bad.fixup	arguments passed as ... to rq.
method, global, xscalm, control, lower, upper, reltol, rel.tol	arguments passed to related arguments of nleqslv, nlminb and optim, respectively.

Details

Parameter transformations

The arguments param.tf, fwd.tf, deriv.fwd.tf, bwd.tf define the transformations of the variogram parameters for RE(ML) estimation. Implemented are currently "log", "logit1", "logit2", "logit3" (various variants of logit-transformation, see code of function fwd.transf) and "identity" (= no) transformations. These are the possible values that the many arguments of the function param.transf accept (as quoted character strings) and these are the names of the list components returned by fwd.transf, dfwd.transf and bwd.transf. Additional transformations can be implemented by:

1. Extending the function definitions by arguments like

   fwd.tf = fwd.transf(my.fun = function(x) your transformation),
   deriv.fwd.tf = dfwd.transf(my.fun = function(x) your derivative),
   bwd.tf = bwd.transf(my.fun = function(x) your back-transformation),

2. Assigning to a given argument of param.transf the name of the new function, e.g.
   variance = "my.fun".

Note that the values given for the arguments of param.transf must match the names of the functions returned by fwd.transf, dfwd.transf and bwd.transf.

Approximation of covariances of fixed and random effects and residuals

The robustified estimating equations of robust REML depend on the covariances of $\widehat{\boldsymbol{B}}$. These covariances (and the covariances of $\boldsymbol{B}-\widehat{\boldsymbol{B}}$, $\widehat{\boldsymbol{\beta}}$, $\widehat{\boldsymbol{\varepsilon}}$, $\widehat{\boldsymbol{\varepsilon}} + \widehat{\boldsymbol{B}}$) are approximated by expressions that in turn depend on the variances of $\varepsilon$, $\psi(\varepsilon/\tau)$ and the expectation of $\psi'(\varepsilon/\tau) (= \partial / \partial \varepsilon \, \psi(\varepsilon/\tau))$. The arguments error.family.estimation, error.family.cov.effects and
error.family.cov.residuals control what parametric distribution for $\varepsilon$ is used to compute the variance of $\varepsilon$, $\psi(\varepsilon/\tau)$ and the expectation of $\psi'(\varepsilon/\tau)$ when

• solving the estimating equations (error.family.estimation),

• computing the covariances of $\widehat{\boldsymbol{\beta}}$, $\widehat{\boldsymbol{B}}$ and $\boldsymbol{B}-\widehat{\boldsymbol{B}}$ (error.family.cov.effects) and

• computing the covariances of $\widehat{\boldsymbol{\varepsilon}}=\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}} - \widehat{\boldsymbol{B}}$ and $\widehat{\boldsymbol{\varepsilon}} + \widehat{\boldsymbol{B}} =\boldsymbol{Y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}}$
  (error.family.cov.residuals).

Possible options are: "gaussian" or "long.tailed". In the latter case the probability density function of $\varepsilon$ is assumed to be proportional to $1/\tau \, \exp(-\rho_c(\varepsilon/\tau))$, where $\psi_c(x)=\rho_c^\prime(x)$.

Value

control.georob, control.rq, control.nleqslv, control.optim and control.nlminb all create lists with control parameters passed to georob, rq, nleqslv, optim or nlminb, see arguments above and the help pages of the respective functions for information about the components of these lists. Note that the list returned by control.georob contains some components (irwls.initial, tuning.psi.nr, cov.bhat.betahat, aux.cov.pred.target) that cannot be changed by the user.

param.transf generates a list with character strings that define what transformations are used for estimating the variogram parameters, and fwd.transf, bwd.transf and dfwd.transf return lists of functions with forward and backward transformations and the first derivatives of the forward transformations, see section Parameter transformations above.

Author(s)

Andreas Papritz [email protected].

See Also

georobPackage for a description of the model and a brief summary of the algorithms;

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;

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.

Examples

data(meuse)

r.logzn.rob <- georob(log(zinc) ~ sqrt(dist), data = meuse, locations = ~ x + y,
    variogram.model = "RMexp",
    param = c(variance = 0.15, nugget = 0.05, scale = 200),
    tuning.psi = 1, control = control.georob(cov.bhat = TRUE,
    cov.ehat.p.bhat = TRUE, initial.fixef = "rq"), verbose = 2)

qqnorm(rstandard(r.logzn.rob, level = 0)); abline(0, 1)
qqnorm(ranef(r.logzn.rob, standard = TRUE)); abline(0, 1)

---

Generic Cross-validation

Description

Generic function for cross-validating models.

Usage

cv(object, ...)

Arguments

object	any model object.
...	additional arguments as required by the methods.

Value

will depend on the method function used; see the respective documentation.

Author(s)

Andreas Papritz [email protected].

See Also

georob for (robust) fitting of spatial linear models;

cv.georob for assessing the goodness of a model fitted by georob.

---

Cross-Validating a Spatial Linear Model Fitted by georob

Description

This function assesses the goodness-of-fit of a spatial linear model by K-fold cross-validation. In more detail, the model is re-fitted K times by robust (or Gaussian) (RE)ML, excluding each time 1/Kth of the data. The re-fitted models are used to compute robust (or customary) external Kriging predictions for the omitted observations. If the response variable is log-transformed then the Kriging predictions can be optionally transformed back to the original scale of the measurements. S3methods for evaluating and plotting diagnostic summaries of the cross-validation errors are described for the function validate.predictions.

Usage

## S3 method for class 'georob'
cv(object, formula = NULL, subset = NULL,
    method = c("block", "random"), nset = 10L, seed = NULL,
    sets = NULL, duplicates.in.same.set = TRUE, re.estimate = TRUE,
    param = object[["variogram.object"]][[1]][["param"]],
    fit.param = object[["variogram.object"]][[1]][["fit.param"]],
    aniso = object[["variogram.object"]][[1]][["aniso"]],
    fit.aniso = object[["variogram.object"]][[1]][["fit.aniso"]],
    variogram.object = NULL,
    use.fitted.param = TRUE, return.fit = FALSE,
    reduced.output = TRUE, lgn = FALSE,
    mfl.action = c("offset", "stop"),
    ncores = min(nset, parallel::detectCores()), verbose = 0, ...)

Arguments

object	an object of class of "georob", see georobObject.
formula	an optional formula for the regression model passed by update to georob.
subset	an optional vector specifying a subset of observations to be used in the fitting process.
method	a character keyword, controlling whether subsets are formed by partitioning data set into contiguous spatial blocks by kmeans (default) or randomly. Ignored if sets is non-NULL.
nset	a positive integer defining the number K of subsets into which the data set is partitioned (default: nset = 10). Ignored if sets is non-NULL.
seed	an optional integer seed to initialize random number generation, see set.seed. Ignored if sets is non-NULL.
sets	an optional vector of the same length as the response vector of the fitted model and with positive integers taking values in $(1,2,\ldots,K)$, defining in this way the $K$ subsets into which the data set is split. If sets = NULL (default) the partition is randomly generated by kmeans or runif (using possibly seed).
duplicates.in.same.set	a logical scalar controlling whether replicated observations at a given location are assigned to the same subset when partitioning the data (default TRUE).
re.estimate	a logical scalar controlling whether the model is re-fitted to the reduced data sets before computing the Kriging predictions (TRUE, default) or whether the model passed in object is used to compute the predictions for the omitted observations, see Details.
param	a named numeric vector or a matrix or data frame with initial values of variogram parameters passed by update to georob. If param is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[1]][["param"]]) columns with initial values of variogram parameters for the nset cross-validation sets, and colnames(param) must match names(object[["variogram.object"]][[1]][["param"]]).
fit.param	a named logical vector or a matrix or data frame defining which variogram parameters should be adjusted by update. If fit.param is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[1]][["fit.param"]]) columns with variogram parameter fitting flags for the nset cross-validation sets, and colnames(param) must match names(object[["variogram.object"]][[1]][["fit.param"]]).
aniso	a named numeric vector or a matrix or data frame with initial values of anisotropy parameters passed by update to georob. If aniso is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[1]][["aniso"]]) columns with initial values of anisotropy parameters for the nset cross-validation sets, and colnames(aniso) must match names(object[["variogram.object"]][[1]][["aniso"]]).
fit.aniso	a named logical vector or a matrix or data frame defining which anisotropy parameters should be adjusted by update. If fit.aniso is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[1]][["fit.aniso"]]) columns with anisotropy parameter fitting flags for the nset cross-validation sets, and colnames(param) must match names(object[["variogram.object"]][[1]][["fit.aniso"]]).
variogram.object	an optional list that gives initial values for fitting a possibly nested variogram model for the cross-validation sets. Each component is a list with the following components: param: an optional named numeric vector or a matrix or data frame with initial values of variogram parameters passed by update to georob. If param is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[i]][["param"]]) columns with initial values of variogram parameters for the nset cross-validation sets (i is the ith variogram structure), and colnames(param) must match names(object[["variogram.object"]][[i]][["param"]]). fit.param: an optional named logical vector or a matrix or data frame defining which variogram parameters should be adjusted by update. If fit.param is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[i]][["fit.param"]]) columns with variogram parameter fitting flags for the nset cross-validation sets (i is the ith variogram structure), and colnames(param) must match names(object[["variogram.object"]][[i]][["fit.param"]]). aniso: an optional named numeric vector or a matrix or data frame with initial values of anisotropy parameters passed by update to georob. If aniso is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[i]][["aniso"]]) columns with initial values of anisotropy parameters for the nset cross-validation sets (i is the ith variogram structure), and colnames(aniso) must match names(object[["variogram.object"]][[i]][["aniso"]]). fit.aniso: an optional named logical vector or a matrix or data frame defining which anisotropy parameters should be adjusted by update. If fit.aniso is a matrix (or a data frame) then it must have nset rows and length(object[["variogram.object"]][[i]][["fit.aniso"]]) columns with anisotropy parameter fitting flags for the nset cross-validation sets(i is the ith variogram structure), and colnames(param) must match names(object[["variogram.object"]][[i]][["fit.aniso"]]).
use.fitted.param	a logical scalar controlling whether fitted values of param (and aniso are used as initial values when variogram parameters are fitted for the cross-validation sets (default TRUE).
return.fit	a logical scalar controlling whether information about the fit should be returned when re-estimating the model with the reduced data sets (default FALSE).
reduced.output	a logical scalar controlling whether the complete fitted model objects, fitted to the reduced data sets, are returned (FALSE) or only some components (TRUE, default, see Value). Ignored if return.fit = FALSE.
lgn	a logical scalar controlling whether Kriging predictions of a log-transformed response should be transformed back to the original scale of the measurements (default FALSE).
mfl.action	a character keyword controlling what is done when some levels of factor(s) are not present in any of the subsets used to fit the model. The function either stops ("stop") or treats the respective factors as model offset ("offset", default).
ncores	a positive integer controlling how many cores are used for parallelized computations, see Details.
verbose	a positive integer controlling logging of diagnostic messages to the console during model fitting. Passed by update to georob.
...	additional arguments passed by update to georob, see Details.

Details

Note that the data frame passed as data argument to georob must exist in the user workspace when calling cv.georob.

cv.georob uses the packages parallel and snowfall for parallelized computations. By default, the function uses $K$ CPUs but not more than are physically available (as returned by detectCores).

cv.georob uses the function update to re-estimated the model with the reduced data sets. Therefore, any argument accepted by georob except data can be changed when re-fitting the model. Some of them (e.g. formula, subset, etc.) are explicit arguments of cv.georob, but also the remaining ones can be passed by ... to the function.

Practitioners in geostatistics commonly cross-validate a fitted model without re-estimating the model parameters with the reduced data sets. This is clearly an unsound practice (see Hastie et al., 2009, sec. 7.10). Therefore, the argument re.estimate should always be set to TRUE. The alternative is provided only for historic reasons.

Value

The method cv.georob returns an object of class cv.georob, which is a list with the two components pred and fit.

pred is a data frame with the coordinates and the cross-validation prediction results with the following variables:

subset	an integer vector defining to which of the $K$ subsets an observation was assigned.
data	the values of the (possibly log-transformed) response.
pred	the Kriging predictions.
se	the Kriging standard errors.

If lgn = TRUE then pred has the additional variables:

lgn.data	the untransformed response.
lgn.pred	the unbiased back-transformed predictions of a log-transformed response.
lgn.se	the Kriging standard errors of the back-transformed predictions of a log-transformed response.

The second component fit contains either the full outputs of georob, fitted for the $K$ reduced data sets (reduced.output = FALSE), or $K$ lists with the components tuning.psi, converged,
convergence.code, gradient, variogram.object, coefficients along with the standard errors of $\widehat{\boldsymbol{\beta}}$, see georobObject.

Author(s)

Andreas Papritz [email protected].

References

Hastie, T., Tibshirani, R. and Friedman, J. (2009) The Elements of Statistical Learning; Data Mining, Inference and Prediction, Springer, New York, doi:10.1007/978-0-387-84858-7

See Also

georobPackage for a description of the model and a brief summary of the algorithms;

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;

georobMethods for further methods for the class georob;

predict.georob for computing robust Kriging predictions;

validate.predictions for validating 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.

Examples

## define number of cores for parallel computations
if(interactive()) ncpu <- 10L else ncpu <- 1L

data(meuse)

r.logzn <- georob(log(zinc) ~ sqrt(dist), data = meuse, locations = ~ x + y,
  variogram.model = "RMexp",
  param = c(variance = 0.15, nugget = 0.05, scale = 200),
  tuning.psi = 1000)

if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  r.logzn.cv.1 <- cv(r.logzn, seed = 1, lgn = TRUE, ncores = 1, verbose = 1)

  r.logzn.cv.2 <- cv(r.logzn, formula = .~. + ffreq, seed = 1, lgn = TRUE,
      ncores = ncpu)

  plot(r.logzn.cv.1, type = "bs")
  plot(r.logzn.cv.2, type = "bs", add = TRUE, col = "red")

  legend("topright", lty = 1, col = c("black", "red"), bty = "n",
      legend = c("log(Zn) ~ sqrt(dist)", "log(Zn) ~ sqrt(dist) + ffreq"))
}

---

Setting Default Values of Variogram Parameters

Description

Auxiliary functions to set sensible default values for anisotropy parameters and for controlling what variogram and anisotropy parameters should be estimated.

Usage

default.aniso(f1 = 1., f2 = 1., omega = 90., phi = 90., zeta = 0.)

default.fit.param(
  variance = TRUE, snugget = FALSE, nugget = TRUE, scale = TRUE,
  alpha = FALSE, beta = FALSE, delta = FALSE, gamma = FALSE,
  kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE)

default.fit.aniso(f1 = FALSE, f2 = FALSE, omega = FALSE,
  phi = FALSE, zeta = FALSE)

Arguments

variance	variance (sill $\sigma^2$) of the auto-correlated component of the Gaussian random field $B(\boldsymbol{s})$.
snugget	variance (spatial nugget $\sigma^2_{\mathrm{n}}$) of the seemingly spatially uncorrelated component of $B(\boldsymbol{s})$ (micro-scale spatial variation; default value snugget = 0).
nugget	variance (nugget $\tau^2$) of the independent error $\varepsilon(\boldsymbol{s})$.
scale	range parameter ($\alpha$) of the variogram.
alpha, beta, delta, gamma, kappa, lambda, mu, nu	names of additional variogram parameters such as the smoothness parameter $\nu$ of the Whittle-Matérn model (see gencorr and param.names).
f1	positive ratio $f_1$ of lengths of second and first semi-principal axes of an ellipsoidal surface with constant semi-variance in $\mathrm{I}\!\mathrm{R}^3$ (default f1 = 1), see subsection Model of georobPackage.
f2	positive ratio $f_2$ of lengths of third and first semi-principal axes of the semi-variance ellipsoid (default f2 = 1), see subsection Model of georobPackage.
omega	azimuth in degrees of the first semi-principal axis of the semi-variance ellipsoid (default omega = 90), see subsection Model of georobPackage.
phi	90 degrees minus altitude of the first semi-principal axis of the semi-variance ellipsoid (default phi = 90), see subsection Model of georobPackage.
zeta	angle in degrees between the second semi-principal axis and the direction of the line defined by the intersection between the $x$-$y$-plane and the plane orthogonal to the first semi-principal axis of the semi-variance ellipsoid through the origin (default zeta = 0), see subsection Model of georobPackage.

Value

Either a named numeric vector with initial values of anisotropy parameters (default.aniso) or named logical vectors, controlling what parameters should be estimated (default.fit.param,
default.fit.aniso).

Author(s)

Andreas Papritz [email protected].

See Also

georobPackage for a description of the model and a brief summary of the algorithms;

georob for (robust) fitting of spatial linear models.

Examples

default.aniso(f1 = 0.5, omega = 45)
default.fit.param(scale=FALSE, alpha = TRUE)
default.fit.aniso(f1 = TRUE, omega = TRUE)

---

Elevation Data

Description

Surface elevation data taken from Davis (1972).

Usage

data(elevation)

Format

A data frame with 52 observations on the following 3 variables:

x

a numeric vector with the easting coordinate in multiplies of 50 feet.

y

a numeric vector with the northing coordinate in multiplies of 50 feet..

height

a numeric vector with the elevation in multiples of 10 feet.

Note

The data were imported from the package geoR.

Source

Davis, J.C. (1973) Statistics and Data Analysis in Geology, Wiley, New York.

Examples

data(elevation)
summary(elevation)

---

Fitting Model Functions to Sample Variograms

Description

The function fit.variogram.model fits a variogram model to a sample variogram by (weighted) non-linear least squares. The function control.fit.variogram.model generates a list with steering parameters which control fit.variogram.model. There are print, summary and lines methods for summarizing and displaying fitted variogram models.

Usage

fit.variogram.model(sv,
    variogram.model = c("RMexp", "RMaskey", "RMbessel", "RMcauchy",
        "RMcircular", "RMcubic", "RMdagum", "RMdampedcos", "RMdewijsian",
        "RMfbm", "RMgauss", "RMgencauchy", "RMgenfbm", "RMgengneiting",
        "RMgneiting", "RMlgd", "RMmatern", "RMpenta", "RMqexp",
        "RMspheric", "RMstable", "RMwave", "RMwhittle"),
    param, fit.param = default.fit.param()[names(param)],
	  aniso = default.aniso(), fit.aniso = default.fit.aniso(),
    variogram.object = NULL,
    max.lag = max(sv[["lag.dist"]]), min.npairs = 30,
    weighting.method = c("cressie", "equal", "npairs"),
    control = control.fit.variogram.model(),
    verbose = 0)

control.fit.variogram.model(maximizer = c("nlminb", "optim"),
    param.tf = param.transf(), fwd.tf = fwd.transf(),
    deriv.fwd.tf = dfwd.transf(), bwd.tf = bwd.transf(),
    hessian = TRUE, optim = control.optim(), nlminb = control.nlminb())

## S3 method for class 'fitted.variogram'
print(x, digits = max(3, getOption("digits") - 3), ...)

## S3 method for class 'fitted.variogram'
summary(object, correlation = FALSE, signif = 0.95, ...)

## S3 method for class 'fitted.variogram'
lines(x, what = c("variogram", "covariance", "correlation"),
    from = 1.e-6, to, n = 501, xy.angle = 90, xz.angle = 90,
    col = 1:length(xy.angle), pch = 1:length(xz.angle), lty = "solid", ...)

Arguments

sv	an object of class sample.variogram, see sample.variogram.
variogram.model	a character keyword defining the variogram model to be fitted. Currently, most basic variogram models provided formerly by the now archived package RandomFields can be fitted (see Details of georob and gencorr).
param	a named numeric vector with initial values of the variogram parameters. The following parameter names are allowed (see Details of georob and georobPackage for information about the parametrization of variogram models): variance: variance (sill $\sigma^2$) of the auto-correlated component of the Gaussian random field $B(\boldsymbol{s})$. snugget: variance (spatial nugget $\sigma^2_{\mathrm{n}}$) of the seemingly spatially uncorrelated component of $B(\boldsymbol{s})$ (micro-scale spatial variation; default value snugget = 0). nugget: variance (nugget $\tau^2$) of the independent errors $\varepsilon(\boldsymbol{s})$. scale: range parameter ($\alpha$) of the variogram. names of additional variogram parameters such as the smoothness parameter $\nu$ of the Whittle-Matérn model (see gencorr and param.names).
fit.param	a named logical vector (or a function such as default.fit.param that creates this vector) with the same names as used for param, defining which parameters are adjusted (TRUE) and which are kept fixed at their initial values (FALSE) when fitting the model.
aniso	a named numeric vector with initial values (or a function such as default.aniso that creates this vector) for fitting geometrically anisotropic variogram models. The names of aniso are matched against the following names (see Details and georobPackage for information about the parametrization of variogram models): f1: ratio $f_1$ of lengths of second and first semi-principal axes of an ellipsoidal surface with constant semi-variance in $\mathrm{I}\!\mathrm{R}^3$ (default f1 = 1). f2: ratio $f_2$ of lengths of third and first semi-principal axes of the semi-variance ellipsoid (default f2 = 1). omega: azimuth in degrees of the first semi-principal axis of the semi-variance ellipsoid (default omega = 90). phi: 90 degrees minus altitude of the first semi-principal axis of the semi-variance ellipsoid (default phi = 90). zeta: angle in degrees between the second semi-principal axis and the direction of the line defined by the intersection between the $x$-$y$-plane and the plane orthogonal to the first semi-principal axis of the semi-variance ellipsoid through the origin (default zeta = 0).
fit.aniso	a named logical vector (or a function such as default.fit.aniso that creates this vector) with the same names as used for aniso, defining which parameters are adjusted (TRUE) and which are kept fixed at their initial values (FALSE) when fitting the model.
variogram.object	an optional list that defines a possibly nested variogram model. Each component is itself a list with the following components: variogram.model: a character keyword defining the variogram model, see respective argument above. param: a named numeric vector with initial values of the variogram parameters, see respective argument above. fit.param: a named logical vector defining which parameters are adjusted, see respective argument above. aniso: a named numeric vector with initial values for fitting geometrically anisotropic variogram models, see respective argument above. fit.param: a named logical vector defining which anisotropy parameters are adjusted, see respective argument above. Note that the arguments variogram.model, param, fit.param, aniso and fit.aniso are ignored when variogram.object is passed to fit.variogram.model.
max.lag	a positive numeric defining the maximum lag distance to be used for fitting or plotting variogram models (default all lag classes).
min.npairs	a positive integer defining the minimum number of data pairs required so that a lag class is used for fitting a variogram model (default 30).
weighting.method	a character keyword defining the weights for non-linear least squares. Possible values are: "equal": no weighting , "npairs": weighting by number of data pairs in a lag class, "cressie": “Cressie's weights” (default, see Cressie, 1993, sec. 2.6.2).
verbose	a positive integer controlling logging of diagnostic messages to the console during model fitting.
control	a list with the components maximizer, param.tf, fwd.tf, bwd.tf, hessian, optim and nlminb or a function such as control.fit.variogram.model that generates such a list. See control.georob for information on maximizer, param.tf, fwd.tf, bwd.tf, hessian, optim and nlminb.
maximizer	a character keyword defining the optimizer for nonlinear least squares. Possible values are nlminb (default) or optim.
hessian	a logical scalar controlling whether the Hessian should be computed at the nonlinear least squares estimates.
param.tf	a function such as param.transf, which returns a named vector of character strings that define the transformations to be applied to the variogram parameters for model fitting, see control.georob.
fwd.tf	a function such as fwd.transf, which returns a named list of invertible functions to be used to transform variogram parameters, see control.georob.
deriv.fwd.tf	a function such as dfwd.transf, which returns a named list of functions corresponding to the first derivatives of fwd.tf, see control.georob.
bwd.tf	a function such as bwd.transf, which returns the named list of inverse functions corresponding to fwd.tf, see see control.georob.
nlminb	a list of arguments passed to nlminb or a function such as control.nlminb that generates such a list (see nlminb for allowed arguments).
optim	a list of arguments passed to optim or a function such as control.optim that generates such a list (see optim for allowed arguments).
object, x	an object of class fitted.variogram.
digits	a positive integer indicating the number of decimal digits to print.
correlation	a logical scalar controlling whether the correlation matrix of the fitted variogram parameters is computed (default FALSE).
signif	a numeric with the confidence level for computing confidence intervals for variogram parameters (default 0.95).
what	a character keyword with the quantity that should be displayed (default "variogram").
from	a numeric with the minimal lag distance used in plotting variogram models.
to	a numeric with the maximum lag distance used in plotting variogram models (default: largest lag distance of current plot).
n	a positive integer specifying the number of equally spaced lag distances for which semi-variances are evaluated in plotting variogram models (default 501).
xy.angle	a numeric vector with azimuth angles (in degrees, clockwise positive from north) in $x$-$y$-plane for which semi-variances should be plotted.
xz.angle	a numeric vector with angles in $x$-$z$-plane (in degrees, clockwise positive from zenith to south) for which semi-variances should be plotted.
col	a vector with colours of curves to distinguish curves relating to different azimuth angles in $x$-$y$-plane.
pch	a vector with the plotting symbols added to lines to distinguish curves relating to different angles in $x$-$z$-plane.
lty	a vector with the line types for plotting variogram models.
...	additional arguments passed to methods.

Details

The parametrization of geometrically anisotropic variograms is described in detail in georobPackage, and the section Details of georob describes how the parameter estimates are constrained to permissible ranges. The same mechanisms are used in fit.variogram.model.

The method summary computes confidence intervals of the estimated variogram and anisotropy parameters from the Hessian matrix of the residual sums of squares, based on the asymptotic normal distribution of least squares estimates. Note that the Hessian matrix with respect to the transformed variogram and anisotropy parameters is used for this. Hence the inverse Hessian matrix is the covariance matrix of the transformed parameters, confidence intervals are first computed for the transformed parameters and the limits of these intervals are transformed back to the original scale of the parameters. Optionally, summary reports the correlation matrix of the transformed parameters, also computed from the Hessian matrix.

Value

The function fit.variogram.model generates an object of class fitted.variogram which is a list with the following components:

sse	the value of the object function (weighted residual sum of squares) evaluated at the solution.
variogram.object	the estimated parameters of a possibly nested variograms model. This is a list that contains for each variogram model structure the following components: variogram.model: the name of the fitted parametric variogram model. param: a named numeric vector with the (estimated) variogram parameters. fit.param: a named logical vector with the flags defining what variogram parameters were estimated. isotropic: logical indicating whether an isotropic variogram was fitted. aniso: a named numeric vector with the (estimated) anisotropy parameters. fit.aniso: a named logical vector with the flags defining what anisotropy parameters were estimated. sincos: a list with sin and cos of the angles $\omega$, $\phi$ and $\zeta$ that define the orientation of the anisotropy ellipsoid. rotmat: the matrix $(\boldsymbol{C}_1, \boldsymbol{C}_2, \boldsymbol{C}_3)$ (see georobPackage). sclmat: a vector with the elements 1, $1/f_1$, $1/f_2$ (see georobPackage).
param.tf	a character vector listing the transformations of the variogram parameters used for model fitting.
fwd.tf	a list of functions for variogram parameter transformations.
bwd.tf	a list of functions for inverse variogram parameter transformations.
converged	a logical scalar indicating whether numerical maximization by nlminb or optim converged.
convergence.code	a diagnostic integer issued by nlminb or optim (component convergence) about convergence.
iter	a named integer vector of length two with the number of function and gradient evaluations by nlminb or optim.
call	the matched call.
residuals	a numeric vector with the residuals, that is the sample semi-variance minus the fitted values.
fitted	a numeric vector with the modelled semi-variances.
weights	a numeric vector with the weights used for fitting.
hessian.tfpa	a symmetric matrix with the Hessian at the solution with respect to the transformed variogram and anisotropy parameters (missing if hessian is false). This Hessian is used by summary.fitted.variogram to compute confidence intervals for the estimated parameters.
hessian.ntfpa	a symmetric matrix with the Hessian at the solution with respect to the non-transformed variogram and anisotropy parameters (missing if hessian is false).

The function control.fit.variogram.model returns a list with parameters to steer
fit.variogram.model, see arguments of the function above for its components.

The method print.fitted.variogram invisibly returns the fitted variogram model unchanged.

The method summary.fitted.variogram returns an object of class summary.fitted.variogram which is a list containing a subset of the components of the fitted variogram object (call, residuals, weights, converged, convergence.code, iter, sse, variogram.object), the matrix param.aniso with the estimated values of the variogram parameters along with the bounds of the confidence intervals and optionally the correlation matrix cor.tf.param of the estimated transformed parameters. There is a print method for objects of class summary.fitted.variogram which returns invisibly the object unchanged.

The method lines.fitted.variogram is called for its side effects and returns the value NULL invisibly.

Author(s)

Andreas Papritz [email protected].

References

Cressie, N. A. C. (1993) Statistics for Spatial Data, Wiley, New York, doi:10.1002/9781119115151.

See Also

georobPackage for a description of the model and a brief summary of the algorithms; 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.

Examples

data(wolfcamp)

## fitting an isotropic IRF(0) model
r.sv.iso <- sample.variogram(pressure~1, data = wolfcamp,
    locations = ~x + y, lag.dist.def = seq(0, 200, by = 15))
plot(r.sv.iso, type = "l")

if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  r.irf0.iso <- fit.variogram.model(r.sv.iso, variogram.model = "RMfbm",
      param = c(variance = 100, nugget = 1000, scale = 1., alpha = 1.),
      fit.param = default.fit.param(scale = FALSE, alpha = TRUE))
  summary(r.irf0.iso, correlation = TRUE)
  lines(r.irf0.iso, line.col = "red")
}

## fitting an anisotropic IRF(0) model
r.sv.aniso <- sample.variogram(pressure~1, data = wolfcamp,
    locations = ~x + y, lag.dist.def = seq(0, 200, by = 15),
    xy.angle.def = c(0., 22.5, 67.5, 112.5, 157.5, 180.))
plot(r.sv.aniso, type = "l")

if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  r.irf0.aniso <- fit.variogram.model(r.sv.aniso, variogram.model = "RMfbm",
      param = c(variance = 100, nugget = 1000, scale = 1., alpha = 1.5),
      fit.param = default.fit.param(scale = FALSE, alpha = TRUE),
      aniso = default.aniso(f1 = 0.4, omega = 135.),
      fit.aniso = default.fit.aniso(f1 = TRUE, omega = TRUE),
      control = control.fit.variogram.model(
        maximizer = "optim",
        optim = control.optim(
          method = "BFGS", hessian = TRUE, control = list(maxit = 5000)
        )
      ))
  summary(r.irf0.aniso, correlation = TRUE)

  lines(r.irf0.aniso, xy.angle = seq(0, 135, by = 45))
}

---

Variogram Models

Description

The function gencorr computes intrinsic or stationary isotropic generalized correlations (= negative semi-variances computed with sill (variance) parameter equal to 1) for a set of basic variogram models formerly made available by the function RFfctn of the now archived R package RandomFields.

Usage

gencorr(x, variogram.model, param)

Arguments

x	a numeric vector with scaled lag distances, i.e. lag distances divided by the range parameter param["scale"].
variogram.model	a character keyword defining the variogram model. Currently, the following models are 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", see Details.
param	a named numeric vector with values of the additional parameters of the variogram models such as the smoothness parameter of the Whittle-Matérn model, see param.names for the names of these parameters. Note that some variogram models (e.g. "RMcircular") do not have any additional parameters.

Details

The name and parametrization of the variogram models originate from the function RFfctn of RandomFields. The equations and further information are taken (with minor modifications) from the help pages of the respective functions of the archived R package RandomFields, version 3.3.14 (Schlather et al., 2022). Note that the variance (sill, param["variance"]) and the range parameters (param["scale"]) are assumed to be equal to 1 in the following formulae, and $x$ is the lag distance. The variogram functions are stationary and are valid for any number of dimensions if not mentioned otherwise.

The following intrinsic or stationary isotropic variogram functions $\gamma(x)$ are implemented in gencorr:

• RMaskey

  $\gamma(x)= 1 - (1-x)^\alpha 1_{[0,1]}(x)$

  $1_{[0,1]}(x)$ is the indicator function equal to 1 for $x \in [0,1]$ and 0 otherwise. This variogram function is valid for dimension $d$ if $\alpha \ge (d+1)/2$. For $\alpha=1$ we get the well-known triangle (or tent) model, which is only valid on the real line.

• RMbessel

  $\gamma(x) = 1 - 2^\nu \Gamma(\nu+1) x^{-\nu} J_\nu(x)$

  where $\nu \ge \frac{d-2}2$, $\Gamma$ denotes the gamma function and $J_\nu$ is a Bessel function of first kind. This models a hole effect (see Chilès and Delfiner, 1999, p. 92). An important case is $\nu=-0.5$ which gives the variogram function

  $\gamma(x)= 1 - \cos(x)$

  and which is only valid for $d=1$ (this equals RMdampedcos for $\lambda = 0$). A second important case is $\nu=0.5$ with variogram function

  $\gamma(x) = \left(1 - \frac{\sin(x)}{x}\right) 1_{x>0}$

  which is valid for $d \le 3$. This coincides with RMwave.

• RMcauchy

  $\gamma(x) = 1 - (1 + x^2)^{-\gamma}$

  where $\gamma > 0$. The parameter $\gamma$ determines the asymptotic power law. The smaller $\gamma$, the longer the long-range dependence. The generalized Cauchy Family (RMgencauchy) includes this family for the choice $\alpha = 2$ and $\beta = 2 \gamma$.

• RMcircular

  $\gamma(x) = 1 - \left(1 -\frac{2}{\pi} \left(x \sqrt{1-x^2} + \arcsin(x)\right)\right) 1_{[0,1]}(x)$

  This variogram function is valid only for dimensions $d \le 2$.

• RMcubic

  $\gamma(x) = 1 - (1-7 x^2 + 8.75 x^3 - 3.5 x^5 + 0.75 x^7) 1_{[0,1]}(x)$

  The model is only valid for dimensions $d \le 3$. It is a 2 times differentiable variogram function with compact support (see Chilès and Delfiner, 1999, p. 84).

• RMdagum

  $\gamma(x) = (1+x^{-\beta})^{-\gamma / \beta}$

  The parameters $\beta$ and $\gamma$ can be varied in the intervals $(0,1]$ and $(0,1)$, respectively. Like the generalized Cauchy model (RMgencauchy) the Dagum family can be used to model separately fractal dimension and Hurst effect (see Berg et al., 2008).

• RMdampedcos

  $\gamma(x) = 1 - \exp(-\lambda x) \cos(x)$

  The model is valid for any dimension $d$. However, depending on the dimension of the random field the following bound $\lambda \ge 1/{\tan(\pi/(2d))}$ has to be respected. This variogram function models a hole effect (see Chilès and Delfiner, 1999, p. 92). For $\lambda = 0$ we obtain

  $\gamma(x)= 1 - \cos(x)$

  which is only valid for $d=1$ and corresponds to RMbessel for $\nu=-0.5$.

• RMdewijsian

  $\gamma(x) = \log(1 + x^{\alpha})$

  where $\alpha \in (0,2]$. This is an intrinsic variogram function. Originally, the logarithmic model $\gamma(x) = \log(x)$ was named after de Wijs and reflects a principle of similarity (see Chilès and Delfiner, 1999, p. 90). But note that $\gamma(x) = \log(x)$ is not a valid variogram function.

• RMexp

  $\gamma(x) = 1 - e^{-x}$

  This model is a special case of the Whittle model (RMwhittle) if $\nu=0.5$ and of the stable family (RMstable) if $\nu = 1$. Moreover, it is the continuous-time analogue of the first order auto-regressive time series covariance structure.

• RMfbm

  $\gamma(x) = x^\alpha$

  where $\alpha \in (0,2)$. This is an intrinsically stationary variogram function. For $\alpha=1$ we get a variogram function corresponding to a standard Brownian Motion. For $\alpha \in (0,2)$ the quantity $H = \frac{\alpha}{2}$ is called Hurst index and determines the fractal dimension $D = d + 1 - H$ of the corresponding Gaussian sample paths where $d$ is the dimension of the random field (see Chilès and Delfiner, 1999, p. 89).

• RMgauss

  $\gamma(x) = 1 - e^{-x^2}$

  The Gaussian model has an infinitely differentiable variogram function. This smoothness is artificial. Furthermore, this often leads to singular matrices and therefore numerically instable procedures (see Stein, 1999, p. 29). The Gaussian model is included in the stable class (RMstable) for the choice $\alpha = 2$.

• RMgencauchy

  $\gamma(x) = 1 - (1 + x^\alpha)^{-\beta/\alpha}$

  where $\alpha \in (0,2]$ and $\beta > 0$. This model has a smoothness parameter $\alpha$ and a parameter $\beta$ which determines the asymptotic power law. More precisely, this model admits simulating random fields where fractal dimension D of the Gaussian sample path and Hurst coefficient H can be chosen independently (compare also with RMlgd): Here, we have $D = d + 1 - \alpha/2, \alpha \in (0,2]$ and $H = 1 - \beta/2, \beta > 0$. The smaller $\beta$, the longer the long-range dependence. The variogram function is very regular near the origin, because its Taylor expansion only contains even terms and reaches its sill slowly. Note that the Cauchy Family (RMcauchy) is included in this family for the choice $\alpha = 2$ and $\beta = 2 \gamma$.

• RMgenfbm

  $\gamma(x) = (1 + x^{\alpha})^{\delta/\alpha} - 1$

  where $\alpha \in (0,2)$ and $\delta \in (0,1)$. This is an intrinsic variogram function.

• RMgengneiting This is a family of stationary models whose elements are specified by the two parameters $\kappa$ and $\mu$ with $\kappa$ being a non-negative integer and $\mu \ge \frac{d}{2}$ with $d$ denoting the dimension of the random field (the models can be used for any dimension). Let $\beta = \mu + 2\kappa +1/2$.

  For $\kappa = 0$ the model equals the Askey model (RMaskey) and is therefore not implemented.

  For $\kappa = 1$ the model is given by

  $\gamma(x) = 1 - \left(1+\beta x \right)(1-x)^{\beta} 1_{[0,1]}(x), \qquad \beta = \mu +2.5,$

  If $\kappa = 2$

  $\gamma(x) = 1 - \left(1 + \beta x + \frac{\beta^{2} - 1}{3} x^{2} \right)(1-x)^{\beta} 1_{[0,1]}(x), \qquad \beta = \mu+4.5,$

  and for $\kappa = 3$

  $\gamma(x) = 1 - \left( 1 + \beta x + \frac{(2\beta^{2}-3)}{5} x^{2}+ \frac{(\beta^2 - 4)\beta}{15} x^{3} \right)(1-x)^\beta 1_{[0,1]}(x), \beta = \mu+6.5,$

  A special case of this family is RMgneiting (with $s = 1$ there) for the choice $\kappa = 3, \mu = 3/2$.

• RMgneiting

  $\gamma(x) = 1 - (1 + 8 s x + 25 s^2 x^2 + 32 s^3 x^3)(1-s x)^8$

  if $0 \le x \le \frac{1}{s}$ and

  $\gamma(x)= 1$

  otherwise. Here, $s=0.301187465825$. This variogram function is valid only for dimensions less than or equal to 3. It is 6 times differentiable and has compact support. This model is an alternative to RMgauss as its graph is hardly distinguishable from the graph of the Gaussian model, but possesses neither the mathematical nor the numerical disadvantages of the Gaussian model. It is a special case of RMgengneiting for the choice $\kappa=3, \mu=1.5$.

• RMlgd

  $\gamma(x) = \frac{\beta}{\alpha + \beta} x^{\alpha} 1_{[0,1]}(x) + (1 - \frac{\alpha}{\alpha + \beta} x^{-\beta}) 1_{x>1}(x)$

  where $\beta >0$ and $0 < \alpha \le (3-d)/2$, with $d$ denoting the dimension of the random field. The model is only valid for dimension $d=1,2$. This model admits simulating random fields where fractal dimension $D$ of the Gaussian sample and Hurst coefficient $H$ can be chosen independently (compare also RMgencauchy): Here, the random field has fractal dimension $D = d+1 - \alpha/2$ and Hurst coefficient $H = 1-\beta/2$ for $0< \beta \le 1$.

• RMmatern

  $\gamma(x) = 1 - \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu}x)^\nu K_\nu(\sqrt{2\nu}x)$

  where $\nu > 0$ and $K_\nu$ is the modified Bessel function of second kind. This is one of 3 possible parametrizations (Whittle, Matérn, Handcock-Wallis) of the Whittle-Matérn model. The Whittle-Matérn model is the model of choice if the smoothness of a random field is to be parametrized: the sample paths of a Gaussian random field with this covariance structure are $m$ times differentiable if and only if $\nu > m$ (see Gneiting and Guttorp, 2010, p. 24). Furthermore, the fractal dimension $D$ of the Gaussian sample paths is determined by $\nu$: We have $D = d + 1 - \nu, \nu \in (0,1)$ and $D = d$ for $\nu > 1$ where $d$ is the dimension of the random field (see Stein, 1999, p. 32). If $\nu=0.5$ the Matérn model equals RMexp. For $\nu$ tending to $\infty$ a rescaled Gaussian model RMgauss appears as limit for the Handcock-Wallis parametrization.

• RMpenta

  $\gamma(x) = 1 - \left(1 - \frac{22}{3}x^{2} + 33 x^{4} - \frac{77}{2} x^{5} + \frac{33}{2} x^{7} - \frac{11}{2} x^{9} + \frac{5}{6}x^{11}\right) 1_{[0,1]}(x)$

  The model is only valid for dimensions $d \le 3$. It has a 4 times differentiable variogram function with compact support (cf. Chilès and Delfiner, 1999, p. 84).

• RMqexp

  $\gamma(x)= 1 - \frac{2 e^{-x} - \alpha e^{-2x}}{ 2 - \alpha }$

  where $\alpha \in [0,1]$.

• RMspheric

  $\gamma(x) = 1 - \left(1 - \frac{3}{2} x + \frac{1}{2} x^3\right) 1_{[0,1]}(x)$

  This variogram model is valid only for dimensions less than or equal to 3 and has compact support.

• RMstable

  $\gamma(x) = 1 - e^{-x^\alpha}$

  where $\alpha \in (0,2]$. The parameter $\alpha$ determines the fractal dimension $D$ of the Gaussian sample paths: $D = d + 1 - \alpha/2$ where $d$ is the dimension of the random field. For $\alpha < 2$ the Gaussian sample paths are not differentiable (see Gelfand et al., 2010, p. 25). The stable family includes the exponential model (RMexp) for $\alpha = 1$ and the Gaussian model ( RMgauss) for $\alpha = 2$.

• RMwave

  $\gamma(x) = \left(1 - \frac{\sin(x)}{x}\right) 1_{x>0}$

  The model is only valid for dimensions $d \le 3$. It is a special case of RMbessel for $\nu = 0.5$. This variogram models a hole effect (see Chilès and Delfiner, 1999, p. 92).

• RMwhittle

  $\gamma(x)=1 - \frac{2^{1- \nu}}{\Gamma(\nu)} x^{\nu}K_{\nu}(x)$

  where $\nu > 0$ and $K_\nu$ is the modified Bessel function of second kind. This is one of 3 possible parametrizations (Whittle, Matérn, Handcock-Wallis) of the Whittle-Matérn model, for further details, see information for entry RMmatern above.

Value

A numeric vector with generalized correlations (= negative semi-variances computed with variance parameter param["variance"] = 1).

Author(s)

Andreas Papritz [email protected]

References

Berg, C., Mateau, J., Porcu, E. (2008) The dagum family of isotropic correlation functions, Bernoulli, 14, 1134–1149, doi:10.3150/08-BEJ139.

Chilès, J.-P., Delfiner, P. (1999) Geostatistics Modeling Spatial Uncertainty, Wiley, New York, doi:10.1002/9780470316993.

Gneiting, T. (2002) Compactly supported correlation functions. Journal of Multivariate Analysis, 83, 493–508, doi:10.1006/jmva.2001.2056.

Gneiting, T., Schlather, M. (2004) Stochastic models which separate fractal dimension and Hurst effect. SIAM review 46, 269–282, doi:10.1137/S0036144501394387.

Gneiting, T., Guttorp, P. (2010) Continuous Parameter Stochastic Process Theory, In Gelfand, A. E., Diggle, P. J., Fuentes, M., Guttrop, P. (Eds.) Handbook of Spatial Statistics, CRC Press, Boca Raton, p. 17–28, doi:10.1201/9781420072884.

Schlather M., Malinowski A., Oesting M., Boecker D., Strokorb K., Engelke S., Martini J., Ballani F., Moreva O., Auel J., Menck P.J., Gross S., Ober U., Ribeiro P., Ripley B.D., Singleton R., Pfaff B., R Core Team (2022). RandomFields: Simulation and Analysis of Random Fields. R package version 3.3.14, https://cran.r-project.org/src/contrib/Archive/RandomFields/.

Stein, M. L. (1999) Interpolation of Spatial Data: Some Theory for Kriging, Springer, New York, doi:10.1007/978-1-4612-1494-6.

See Also

georobPackage for a description of the model and a brief summary of the algorithms;

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.

Examples

## scaled lag distances
x <- seq(0, 3, length.out = 100)

## generalized correlations stable model
y <- gencorr(x, variogram.model = "RMstable", param = c(alpha = 1.5))
plot(x, y)

## generalized correlations circular model
y <- gencorr(x, variogram.model = "RMcircular")
plot(x, y)

---