Description Usage Arguments Details Value Author(s) References See Also Examples
View source: R/georob.exported.functions.R
The function georob
fits a linear model with spatially correlated
errors to geostatistical data that are possibly contaminated by
independent outliers. The regression coefficients and the parameters of
the variogram model are estimated by robust or Gaussian restricted
maximum likelihood (REML) or by Gaussian maximum likelihood (ML).
1 2 3 4 5 6 7 8 9 10 11 12  georob(formula, data, subset, weights, na.action, model = TRUE,
x = FALSE, y = FALSE, contrasts = NULL, offset, locations,
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,
tuning.psi = 2, control = control.georob(),
verbose = 0, ...)

formula 
a symbolic description of the regression model for the
external drift to be fit (mandatory argument). See

data 
an optional data frame, a

subset 
an optional vector specifying a subset of observations to be used in the fitting process. 
weights 
an optional vector of weights to be used in the fitting process, currently ignored. 
na.action 
a function which indicates what should happen when the
data contain 
model, x, y 
logicals. If 
contrasts 
an optional list. See the 
offset 
this optional argument can be used to specify an a
priori known component to be included in the linear predictor during
fitting. An 
locations 
a onesided formula defining the variables that are used as coordinates of the locations were the data was recorded (mandatory argument). 
variogram.model 
a character keyword defining the variogram model
to be fitted. Currently, most basic variogram models
provided by the package RandomFields can be fitted (see Details
and 
param 
a named numeric vector with initial values of the
variogram parameters (mandatory argument). The names of

fit.param 
a named logical vector (or a function such as

aniso 
a named numeric vector with initial values (or a function such as

fit.aniso 
a named logical vector (or a function such as

variogram.object 
an optional list that defines a possibly nested variogram model. Each component is itself a list with the following components:
Note that the arguments 
tuning.psi 
positive numeric. The tuning constant c of the ψ_cfunction of the robust REML algorithm. 
control 
a list specifying parameters that control the behaviour of

verbose 
positive integer controlling logging of diagnostic
messages to the console during model fitting. 
... 
further arguments passed to function (e.g. 
georob
fits a spatial linear model by robust or Gaussian RE(ML)
(Künsch et al., 2011, Künsch
et al., in preparation). georobIntro
describes the
employed model and briefly sketches the robust REML estimation and the
robust external drift Kriging method. Here, we describe further details
of georob
.
Currently, most basic variogram models provided by the
package RandomFields can be fitted by
georob
(see argument variogram.model
for a list of
implemented models). Some of these models have in addition to
variance
, snugget
, nugget
and scale
further
parameters. Initial values of these parameters (param
) and
fitting flags (fit.param
) must be passed to georob
by the
same names as used by the functions RM...
of the package
RandomFields (see RMmodel
). Use the
function param.names
to list additional parameters of a
given variogram.model.
The arguments fit.param
and fit.aniso
are used to control
what variogram and anisotropy parameters are estimated and which are
kept at the constant initial values. The functions
default.fit.param
and default.fit.aniso
set
reasonable default values for these arguments. Note, as an aside, that
the function default.aniso
sets (default) values of the
anisotropy parameters for an isotropic variogram.
The intrinsic variogram model RMfbm
is overparametrized when
both the variance
(plus possibly snugget
) and the
scale
are estimated. Therefore, to estimate the parameters of
this model, scale
must be kept fixed at an arbitrary value by
using fit.param["scale"] = FALSE
.
The subsection Model of georobIntro
describes
how such models are parametrized and gives definitions the various
elements of aniso
. Some additional remarks might be helpful:
The first semiprincipal axis points into the direction with
the farthest reaching autocorrelation, which is described by the range
parameter scale
(α).
The ranges in the direction of the second and third semiprincipal axes are given by f_1α and f_2 α, with 0 < f_2 <= f_1 <= 1.
The default values for aniso
(f_1=1, f_2=1)
define an isotropic variogram model.
Valid ranges for the angles characterizing the orientation of the semivariance ellipsoid are (in degrees): ω [0, 180], φ [0, 180], ζ [90, 90].
Simultaneous estimation of the variance of the microscale variation
(snugget
, sigma^2_n), appears seemingly
as spatially uncorrelated with a given sampling design, and of the
variance (nugget
, τ^2) of the independent errors
requires that for some locations
s_i replicated observations are
available. Locations less or equal than zero.dist
apart are
thereby considered as being coincident (see
control.georob
).
Parameters of variogram models can vary only within certain bounds (see
param.bounds
and RMmodel
for allowed ranges). georob
uses three mechanisms to constrain
parameter estimates to permissible ranges:
Parameter transformations: By default, all variance
(variance
, snugget
, nugget
), the range
scale
, the anisotropy parameters f1
and f2
and
many of the additional parameters are logtransformed before
solving the estimating equations or maximizing the restricted
loglikelihood and this warrants that the estimates are always
positive (see control.georob
for detailed explanations
how to control parameter transformations).
Checking permissible ranges: The additional parameters
of the variogram models such as the smoothness parameter ν
of the WhittleMat¬érn model are forced to stay in the permissible
ranges by signalling an error to nleqslv
, nlminb
or
optim
if the current trial values are invalid. These
functions then graciously update the trial values of the parameters
and carry their task on. However, it is clear that such a
procedure likely gets stuck at a point on the boundary of the
parameter space and is therefore just a workaround for avoiding
runtime errors due to invalid parameter values.
Exploiting the functionality of nlminb
and
optim
: If a spatial model is fitted nonrobustly, then the
arguments lower
, upper
(and method
of
optim
) can be used to constrain the parameters
(see control.optim
how to pass them to optim
).
For optim
one has to use the arguments method =
"LBFGSB"
, lower = l
, upper = u
, where
l and u are numeric vectors with the lower and upper
bounds of the transformed parameters in the order as they
appear in
c(variance, snugget, nugget, scale,
...)[fit.param], aniso[fit.aniso])
,
where ...
are
additional parameters of isotropic variogram models (use
param.names(variogram.model)
to display the names and the
order of the additional parameters for variogram.model
).
To solve the robustified estimating equations for B and β the following initial estimates are used:
hatB=0, if this turns out to be
infeasible, initial values can be passed to georob
by the
argument bhat
of control.georob
.
hatβ is
either estimated robustly by the function
lmrob
, rq
or
nonrobustly by lm
(see argument
initial.fixef
of control.georob
).
Finding the roots of the robustified estimating equations of the
variogram and anisotropy parameters is more sensitive to a good choice
of initial values than maximizing the Gaussian (restricted)
loglikelihood with respect to the same parameters. If the initial
values for param
and aniso
are not sufficiently close to
the roots of the system of nonlinear equations, then
nleqslv
may fail to find them.
Setting
initial.param = TRUE
allows one to find initial values that are
often sufficiently close to the roots so that
nleqslv
converges. This is achieved by:
Initial values of the regression parameters are computed by
lmrob
irrespective of the choice for
initial.fixef
(see control.georob
).
Observations with “robustness weights” of the
lmrob
fit, satisfying
ψ_c(r_i/s)/(r_i/s)
<=min.rweight, are discarded (see
control.georob
).
The model is fit to the pruned data set by Gaussian REML using
optim
.
The resulting estimates of the variogram parameters
(param
, aniso
) are used as initial estimates for the
subsequent robust fit of the model by nleqslv
.
Note that for step 3 above, initial values of param
and
aniso
must be provided to georob
.
Unlike robust REML, where robustified estimating equations are solved
for the variance parameters nugget
(τ^2),
variance
(σ^2), and possibly snugget
(σ_n^2), for Gaussian (RE)ML the
variances can be reparametrized to
the signal variance σ_B^2 = σ^2+σ_n^2,
the inverse relative nugget η = σ_B^2 / τ^2 and
the relative autocorrelated signal variance ξ = σ^2 / σ_B^2.
georob
maximizes then a (restricted) profile
loglikelihood that depends only on η, ξ,
α, ..., and σ_B^2 is estimated by an explicit
expression that depends on these parameters (e.g. Diggle and Ribeiro,
2006, p. 113). This is usually more efficient than maximizing the
(restricted) loglikelihood with respect to the original variance
parameters τ^2, σ_n^2 and
σ^2. georob
chooses the parametrization
automatically, but the user can control it by the argument
reparam
of the function control.georob
.
An object of class georob
representing a robust (or Gaussian) (RE)ML
fit of a spatial linear model. See
georobObject
for the components of the fit.
Andreas Papritz [email protected]
http://www.step.ethz.ch/people/scientificstaff/andreaspapritz.html
with contributions by Cornelia Schwierz.
Diggle, P. J. and Ribeiro, P. J. R. (2006) Modelbased Geostatistics. Springer.
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. http://ecollection.library.ethz.ch/eserv/eth:7080/eth708001.pdf
georobIntro
for a description of the model and a brief summary of the algorithms;
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 backtransformation of Kriging prediction
of logtransformed 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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57  ## Not run:
################
## meuse data ##
################
data(meuse)
## Gaussian REML fit
r.logzn.reml < 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)
summary(r.logzn.reml, correlation = TRUE)
## robust REML fit
r.logzn.rob < update(r.logzn.reml, tuning.psi = 1)
summary(r.logzn.rob, correlation = TRUE)
plot(r.logzn.reml, lag.dist.def = seq(0, 2000, by = 100))
lines(r.logzn.rob, col = "red")
###################
## wolfcamp data ##
###################
data(wolfcamp, package = "geoR")
d.wolfcamp < data.frame(x = wolfcamp[[1]][,1], y = wolfcamp[[1]][,2],
pressure = wolfcamp[[2]])
## fitting isotropic IRF(0) model
r.irf0.iso < georob(pressure ~ 1, data = d.wolfcamp, locations = ~ x + y,
variogram.model = "RMfbm",
param = c(variance = 10, nugget = 1500, scale = 1, alpha = 1.5),
fit.param = default.fit.param(scale = FALSE, alpha = TRUE),
tuning.psi = 1000)
summary(r.irf0.iso)
## fitting anisotropic IRF(0) model
r.irf0.aniso < georob(pressure ~ 1, data = d.wolfcamp, locations = ~ x + y,
variogram.model = "RMfbm",
param = c(variance = 5.9, nugget = 1450, scale = 1, alpha = 1),
fit.param = default.fit.param(scale = FALSE, alpha = TRUE),
aniso = default.aniso(f1 = 0.51, omega = 148.),
fit.aniso = default.fit.aniso(f1 = TRUE, omega = TRUE),
tuning.psi = 1000)
summary(r.irf0.aniso)
plot(r.irf0.iso, lag.dist.def = seq(0, 200, by = 7.5))
plot(r.irf0.aniso, lag.dist.def = seq(0, 200, by = 7.5),
xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180.),
add = TRUE, col = 2:5)
pchisq(2*(r.irf0.aniso[["loglik"]]  r.irf0.iso[["loglik"]]), 2, lower = FALSE)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.