Nothing
## ############################################################################
### georob
georob <-
function(
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,
...
)
{
## wrapper function to georob.fit with "standard" interface for
## statistical modelling
## 2012-04-21 AP
## 2012-05-03 AP bounds for safe parameter values
## 2012-05-07 AP correction of error for constant trend
## 2012-05-28 AP handle missing names of coefficients after calling update
## 2013-04-23 AP new names for robustness weights
## 2013-05-23 AP correct handling of missing observations and to construct model.frame
## 2013-06-03 AP handling design matrices with rank < ncol(x)
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-02 AP new transformation of rotation angles
## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
## 2013-09-06 AP exclusive use of nleqslv for solving estimating equations
## 2014-02-18 AP correcting error when fitting models with offset
## 2014-05-15 AP changes for version 3 of RandomFields
## 2014-05-22 AP correcting error when selecting initial.param
## 2014-06-02 AP partial matching of names of variogram parameters
## 2014-08-18 AP changes for parallelized computations
## 2014-08-18 AP changes for Gaussian ML estimation
## 2014-08-26 AP changes to return ml.method if fitted object
## 2015-06-30 AP function and arguments renamed
## 2015-07-17 AP change of control of computation of hessian for Gaussian (RE)ML,
## changes for singular design matrices, nlminb optimizer added
## 2015-08-19 AP control about error families for computing covariances added
## 2015-08-28 AP computation of hessian suppressed; correction of error when using georob.object;
## control arguments hessian, initial.param, initial.fixef newly organized
## 2015-11-25 AP new way to control which variogram parameters are fitted
## 2016-07-15 AP allowing use of lmrob for computing initial fixed effects for rank-deficient model matrix
## 2016-07-20 AP changes for parallel computations
## 2016-08-08 AP changes for nested variogram models
## 2017-05-07 AP correcting error for interface to rq.fit
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## 2020-02-19 AP correction of error in sanity checks of arguments and for if() and switch()
## 2024-02-21 AP added sfStop()
on.exit( if( sfIsRunning() ) sfStop() )
#### -- check arguments
## check whether all mandatory arguments have been provided
if( missing( formula ) || missing( locations ) ||
( is.null( variogram.object ) && missing( param )) ) stop(
"some mandatory arguments are missing"
)
## check validity of arguments
if(!missing(data)) stopifnot(is.data.frame(as.data.frame(data)))
if(!missing(subset)) stopifnot(is.null(subset) || is.logical(subset) || is.character(subset) || is.numeric(subset))
if(!missing(weights)) stopifnot(is.null(weights) || is.numeric(weights))
if(!missing(offset)) stopifnot(is.null(offset) || is.numeric(offset))
if(!missing(param)){
stopifnot(is.numeric(param))
stopifnot(is.logical(fit.param))
}
stopifnot(identical(length(model), 1L) && is.logical(model))
stopifnot(identical(length(x), 1L) && is.logical(x))
stopifnot(identical(length(y), 1L) && is.logical(y))
stopifnot(is.logical(fit.aniso))
stopifnot(identical(length(tuning.psi), 1L) && is.numeric(tuning.psi) && tuning.psi > 0)
stopifnot(identical(length(verbose), 1L) && is.numeric(verbose) && verbose >= 0)
stopifnot(is.numeric(aniso))
stopifnot(is.list(control))
stopifnot(is.null(contrasts) || is.list(contrasts))
stopifnot(is.null(variogram.object) || is.list(variogram.object))
if( identical( control[["psi.func"]], "t.dist" ) && tuning.psi <= 1. )
stop( "'tuning.psi' must be greater than 1 for t-dist psi-function" )
#### -- check/setup variogram object
## setup or check contents of variogram.object
if( is.null( variogram.object ) ){
## match variogram model
variogram.model <- match.arg( variogram.model )
## match names of param, aniso, fit.param, fit.aniso
## !! CARE: code works only if fit.param has not been used before (lazy evaluation)
if( !missing( param ) ){
tmp <- names( param )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.fit.param() )
)
names( param ) <- tmp
}
if( !missing( fit.param ) ){
tmp <- names( fit.param )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.fit.param() )
)
names( fit.param ) <- tmp
fit.param <- fit.param[names( fit.param ) %in% names( param )]
}
if( !missing( aniso ) ){
tmp <- names( aniso )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.aniso() )
)
names( aniso ) <- tmp
}
if( !missing( fit.aniso ) ){
tmp <- names( fit.aniso )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.aniso() )
)
names( fit.aniso ) <- tmp
}
## create variogram.object
variogram.object <- list(
list(
variogram.model = variogram.model,
param = param, fit.param = fit.param,
aniso = aniso, fit.aniso = fit.aniso
)
)
} else {
if( !missing( param ) ) cat(
"\n information on initial parameters in 'param' ignored because 'variogram.object' specified\n\n"
)
if( !missing( fit.param ) ) cat(
"\n information on initial parameters in 'fit.param' ignored because 'variogram.object' specified\n\n"
)
if( !missing( aniso ) ) cat(
"\n information on initial parameters in 'aniso' ignored because 'variogram.object' specified\n\n"
)
if( !missing( fit.aniso ) ) cat(
"\n information on initial parameters in 'fit.aniso' ignored because 'variogram.object' specified\n\n"
)
variogram.object <- lapply(
variogram.object,
function( y, vm ){
## match names of components
tmp <- names( y )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = c( "variogram.model", "param", "fit.param", "aniso", "fit.aniso" )
)
names( y ) <- tmp
## check whether component variogram.model and param are present
if( is.null( y[["variogram.model"]] ) ) stop(
"'variogram.model' missing in some component of variogram.object"
)
if( is.null( y[["param"]] ) ) stop(
"'param' missing in some component of variogram.object"
)
## match variogram.model
y[["variogram.model"]] <- match.arg( y[["variogram.model"]], vm )
## match names of param, aniso, fit.param, fit.aniso and set
## default values if missing
tmp <- names( y[["param"]] )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.fit.param() )
)
names( y[["param"]] ) <- tmp
if( is.null( y[["fit.param"]] ) ){
y[["fit.param"]] <- default.fit.param()
} else {
tmp <- names( y[["fit.param"]] )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.fit.param() )
)
names( y[["fit.param"]] ) <- tmp
}
y[["fit.param"]] <- y[["fit.param"]][names(y[["fit.param"]] ) %in% names( y[["param"]] )]
if( is.null( y[["aniso"]] ) ){
y[["aniso"]] <- default.aniso()
} else {
tmp <- names( y[["aniso"]] )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.aniso() )
)
names( y[["aniso"]] ) <- tmp
}
if( is.null( y[["fit.aniso"]] ) ){
y[["fit.aniso"]] <- default.fit.aniso()
} else {
tmp <- names( y[["fit.aniso"]] )
tmp <- sapply(tmp, function(x, choices){
match.arg(x, choices)
},
choices = names( default.aniso() )
)
names( y[["fit.aniso"]] ) <- tmp
}
y
}, vm = variogram.model
)
}
#### -- model frame etc.
## get model frame, response vector, weights, offset and design
## matrix (cf. lm, lmrob)
ret.x <- x
ret.y <- y
# ## vector with row number of included observations
#
# in.subset <- 1L:NROW( data )
# if( !missing( subset ) ) in.subset <- in.subset[subset]
## build combined formula for fixed effects and locations
extended.formula <- update(
formula,
paste(
paste( as.character( formula )[c(2L, 1L, 3L)], collapse = " " ),
as.character( locations )[2L], sep = " + "
)
)
## setting-up model frame
cl <- match.call()
mf <- match.call( expand.dots = FALSE )
m <- match(
c( "formula", "data", "subset", "weights", "na.action", "offset"),
names(mf), 0L
)
mf <- mf[c(1L, m)]
mf[["formula"]] <- extended.formula
mf[["drop.unused.levels"]] <- TRUE
mf[[1L]] <- as.name( "model.frame" )
mf <- eval( mf, parent.frame() )
## setting-up terms objects
mt <- terms( formula )
mt.loc <- terms( locations )
## eliminate intercept from mt.loc
attr( mt.loc, "intercept" ) <- 0L
# ## ... and assign fixed effects terms object as attribute to model.frame
#
# attr( mf, "terms" ) <- mt
## check whether 'empty' models have been entered
if( is.empty.model( mt ) )
stop( "an 'empty' fixed effects model has been specified" )
if( is.empty.model( mt.loc ) )
stop( "an 'empty' locations model has been specified" )
## check whether fixed effects model includes an intercept if an
## intrinsic variogram model is used
if( identical( attr( mt, "intercept" ), 0L ) &&
variogram.model %in% control[["irf.model"]] )
stop(
"the fixed effects model must include an intercept ",
"if an unbounded variogram model is used"
)
## extract fixed effects response variable, weights, offset and design matrix
y <- model.response( mf, "numeric" )
w <- as.vector( model.weights( mf ) )
if( !is.null(w) )
stop( "weights are not yet implemented for this estimator" )
offset <- as.vector( model.offset(mf) )
if( !is.null(offset) ) {
if( length( offset ) != NROW(y) )
stop( gettextf(
"number of offsets is %d, should equal %d (number of observations)",
length(offset), NROW(y) ), domain = NA )
}
x <- model.matrix( mt, mf, contrasts )
## check if optionally provided bhat has correct length
if( !is.null( control[["bhat"]] ) && length( y ) != length( control[["bhat"]] ) ) stop(
"lengths of response vector and 'bhat' do not match"
)
## adjust initial.param if all variogram parameters are fixed
if( all( sapply(
variogram.object,
function( x ) !any( c( x[["fit.param"]], x[["fit.aniso"]] ) )
))) control[["initial.param"]] <- FALSE
## adjust choice for initial.fixef to compute regression coefficients
if( tuning.psi < control[["tuning.psi.nr"]] ){
if( control[["initial.param"]] ) control[["initial.fixef"]] <- "lmrob"
} else {
control[["initial.param"]] <- FALSE
}
## check whether design matrix has full column rank
col.rank.XX <- list( deficient = FALSE, rank = NCOL( x ) )
sv <- svd( crossprod( x ) )[["d"]]
min.max.sv <- range( sv )
condnum <- min.max.sv[1L] / min.max.sv[2L]
if( condnum <= control[["min.condnum"]] ){
col.rank.XX[["deficient"]] <- TRUE
col.rank.XX[["col.rank"]] <- sum( sv / min.max.sv[2L] > control[["min.condnum"]] )
if( control[["initial.fixef"]] == "rq" ){
cat(
"design matrix has not full column rank (condition number of X^T X: ",
signif( condnum, 2L ), ")\ninitial values of fixed effects coefficients are computed by 'lmrob' instead of 'rq'\n\n"
)
control[["initial.fixef"]] <- "lmrob"
warning(
"design matrix has not full column rank (condition number of X^T X: ",
signif( condnum, 2L ), ")\ninitial values of fixed effects coefficients are computed by 'lmrob' instead of 'rq'"
)
}
}
## subtract offset
yy <- y
if( !is.null( offset ) ) yy <- yy - offset
#### -- initial values of fixed effects parameters
## compute initial guess of fixed effects parameters (betahat)
r.initial.fit <- switch(
control[["initial.fixef"]],
rq = {
Rho <- function( u, tau) u * (tau - (u < 0))
tau <- control[["rq"]][["tau"]]
process <- (tau < 0 || tau > 1)
if(tau == 0) tau <- control[["req"]][["eps"]]
if(tau == 1) tau <- 1. - control[["req"]][["eps"]]
fit <- switch(
control[["rq"]][["method"]],
br = rq.fit(
x = x, y = yy, tau = tau, method = control[["rq"]][["method"]],
alpha = control[["rq"]][["alpha"]], ci = control[["rq"]][["ci"]],
iid = control[["rq"]][["iid"]], interp = control[["rq"]][["interp"]],
tcrit = control[["rq"]][["tcrit"]]
),
fnb = rq.fit(
x = x, y = yy, tau = tau, method = control[["rq"]][["method"]],
beta = control[["rq"]][["beta"]], eps = control[["rq"]][["eps"]]
),
pfn = rq.fit(
x = x, y = yy, tau = tau, method = control[["rq"]][["method"]],
Mm.factor = control[["rq"]][["Mm.factor"]],
max.bad.fixup = control[["rq"]][["max.bad.fixup"]]
)
)
if (process){
rho <- list(x = fit$sol[1, ], y = fit$sol[3, ])
} else {
if (length(dim(fit$residuals)))
dimnames(fit$residuals) <- list(dimnames(x)[[1]],
NULL)
rho <- sum(Rho(fit$residuals, tau))
}
if( control[["rq"]][["method"]] == "lasso" ){
class(fit) <- c("lassorq", "rq")
} else if( control[["rq"]][["method"]] == "scad"){
class(fit) <- c("scadrq", "rq")
} else {
class(fit) <- ifelse(process, "rq.process", "rq")
}
fit[["na.action"]] <- attr( mf, "na.action" )
fit[["formula"]] <- formula
fit[["terms"]] <- mt
fit[["xlevels"]] <- .getXlevels(mt, mf)
fit[["call"]] <- cl
fit[["tau"]] <- tau
fit[["weights"]] <- w
fit[["residuals"]] <- drop( fit[["residuals"]] )
fit[["rho"]] <- rho
fit[["method"]] <- control[["rq"]][["method"]]
fit[["fitted.values"]] <- drop( fit[["fitted.values"]] )
attr(fit, "na.message") <- attr( m, "na.message" )
if( model ) fit[["model"]] <- mf
fit
},
lmrob = {
fit <- lmrob.fit( x, yy, control = control[["lmrob"]] )
fit[["na.action"]] <- attr(mf, "na.action")
fit[["offset"]] <- offset
fit[["contrasts"]] <- attr(x, "contrasts")
fit[["xlevels"]] <- .getXlevels(mt, mf)
fit[["call"]] <- cl
fit[["terms"]] <- mt
if( control[["lmrob"]][["compute.rd"]] && !is.null(x) )
fit[["MD"]] <- robMD( x, attr( mt, "intercept" ) )
if( !is.null( offset ) ) fit[["fitted.values"]] + offset
fit
},
lm = {
fit <- if( is.null(w) ){
lm.fit(x, y, offset = offset, singular.ok = TRUE )
} else {
lm.wfit(x, y, w, offset = offset, singular.ok = TRUE )
}
class(fit) <- c(if (is.matrix(y)) "mlm", "lm")
fit[["na.action"]] <- attr(mf, "na.action")
fit[["offset"]] <- offset
fit[["contrasts"]] <- attr(x, "contrasts")
fit[["xlevels"]] <- .getXlevels(mt, mf)
fit[["call"]] <- cl
fit[["terms"]] <- mt
if (model) fit[["model"]] <- mf
if (ret.x) fit[["x"]] <- x
if (ret.y) fit[["y"]] <- y
fit[["qr"]] <- NULL
fit
}
)
#### -- prepare coordinates
## compute coordinates of locations and distance object
locations.coords <- model.matrix( mt.loc, mf )
if(
!( missing( aniso ) || missing( fit.aniso ) ) &&
( NCOL( locations.coords ) < 2L || NCOL( locations.coords ) > 3L )
) stop(
"anisotropic variogram models are implemented only for 2 or 3 dimensions"
)
names( yy ) <- rownames( mf )
## check whether argument "object." has been provided in call (e.g. by
## update ) and extract'locations' exists in workspace
extras <- match.call( expand.dots = FALSE )$...
georob.object <- extras[names(extras) %in% "object."]
if(
length( georob.object ) &&
exists( as.character( georob.object ), envir = parent.frame() ) #&&
# !control[["initial.param"]]
){
if( verbose > 4. ) cat(
"\n georob: using some components of 'object.'\n"
)
georob.object <- eval(
georob.object[[1L]], parent.frame()
)[c( "variogram.object", "locations.objects", "Valphaxi.objects" )]
georob.object[["Valphaxi.objects"]] <- c(
list(error=FALSE),
georob.object[["Valphaxi.objects"]]
)
} else {
georob.object <- NULL
}
#### -- initial values of variogram parameters
## compute initial values of variogram and anisotropy parameters
if( tuning.psi < control[["tuning.psi.nr"]] && control[["initial.param"]] ){
if( verbose > 0. ) cat( "\ncomputing robust initial parameter estimates ...\n" )
t.sel <- r.initial.fit[["rweights"]] > control[["min.rweight"]]
if( verbose > 0. ) cat(
"\ndiscarding", sum( !t.sel ), "of", length( t.sel ),
"observations for computing initial estimates of variogram\nand anisotropy parameter by gaussian reml\n"
)
## collect.items for initial object
initial.objects <- list(
x = as.matrix( x[t.sel, ] ),
y = yy[t.sel],
betahat = coef( r.initial.fit ),
bhat = if( is.null( control[["bhat"]] ) ){
rep( 0., length( yy ) )[t.sel]
} else {
control[["bhat"]][t.sel]
},
initial.fit = r.initial.fit,
locations.objects = list(
locations = locations,
coordinates = locations.coords[t.sel, ]
)
)
## estimate model parameters with pruned data set
t.georob <- georob.fit(
initial.objects = initial.objects,
variogram.object = variogram.object,
param.tf = control[["param.tf"]],
fwd.tf = control[["fwd.tf"]],
deriv.fwd.tf = control[["deriv.fwd.tf"]],
bwd.tf = control[["bwd.tf"]],
georob.object = georob.object,
safe.param = control[["safe.param"]],
tuning.psi = control[["tuning.psi.nr"]],
error.family.estimation = control[["error.family.estimation"]],
error.family.cov.effects = control[["error.family.cov.effects"]],
error.family.cov.residuals = control[["error.family.cov.residuals"]],
cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
cov.betahat = control[["cov.betahat"]],
cov.bhat.betahat = control[["cov.bhat.betahat"]],
cov.delta.bhat = control[["cov.delta.bhat"]],
full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
aux.cov.pred.target = control[["aux.cov.pred.target"]],
min.condnum = control[["min.condnum"]],
col.rank.XX = col.rank.XX,
psi.func = control[["psi.func"]],
tuning.psi.nr = tuning.psi,
ml.method = control[["ml.method"]],
maximizer = control[["maximizer"]],
reparam = control[["reparam"]],
irwls.initial = control[["irwls.initial"]],
irwls.maxiter = control[["irwls.maxiter"]],
irwls.ftol = control[["irwls.ftol"]],
force.gradient = control[["force.gradient"]],
zero.dist = control[["zero.dist"]],
control.nleqslv = control[["nleqslv"]],
control.optim = control[["optim"]],
control.nlminb = control[["nlminb"]],
hessian = FALSE,
control.pcmp = control[["pcmp"]],
verbose = verbose
)
variogram.object <- lapply(
1:length(variogram.object),
function( i, x, fvo ){
x <- x[[i]]
fvo <- fvo[[i]]
x[["param"]] <- fvo[["param"]][names(x[["param"]])]
x[["aniso"]] <- fvo[["aniso"]][names(x[["aniso"]])]
x
}, x = variogram.object, fvo = t.georob[["variogram.object"]]
)
}
## collect.items for initial object
initial.objects <- list(
x = as.matrix( x ),
y = yy,
betahat = coef( r.initial.fit ),
bhat = if( is.null( control[["bhat"]] ) ){
rep( 0., length( yy ) )
} else {
control[["bhat"]]
},
initial.fit = r.initial.fit,
locations.objects = list(
locations = locations,
coordinates = locations.coords
)
)
#### -- final estimate of model parameters
## estimate model parameters
if( verbose > 0. ) cat( "\ncomputing final parameter estimates ...\n" )
r.georob <- georob.fit(
initial.objects = initial.objects,
variogram.object = variogram.object,
param.tf = control[["param.tf"]],
fwd.tf = control[["fwd.tf"]],
deriv.fwd.tf = control[["deriv.fwd.tf"]],
bwd.tf = control[["bwd.tf"]],
georob.object = georob.object,
safe.param = control[["safe.param"]],
tuning.psi = tuning.psi,
error.family.estimation = control[["error.family.estimation"]],
error.family.cov.effects = control[["error.family.cov.effects"]],
error.family.cov.residuals = control[["error.family.cov.residuals"]],
cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
cov.betahat = control[["cov.betahat"]],
cov.bhat.betahat = control[["cov.bhat.betahat"]],
cov.delta.bhat = control[["cov.delta.bhat"]],
full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
aux.cov.pred.target = control[["aux.cov.pred.target"]],
min.condnum = control[["min.condnum"]],
col.rank.XX = col.rank.XX,
psi.func = control[["psi.func"]],
tuning.psi.nr = control[["tuning.psi.nr"]],
ml.method = control[["ml.method"]],
maximizer = control[["maximizer"]],
reparam = control[["reparam"]],
irwls.initial = control[["irwls.initial"]],
irwls.maxiter = control[["irwls.maxiter"]],
irwls.ftol = control[["irwls.ftol"]],
force.gradient = control[["force.gradient"]],
zero.dist = control[["zero.dist"]],
control.nleqslv = control[["nleqslv"]],
control.optim = control[["optim"]],
control.nlminb = control[["nlminb"]],
hessian = control[["hessian"]],
control.pcmp = control[["pcmp"]],
verbose = verbose
)
#### -- prepare output
## add offset to fitted values
if( !is.null( offset ) )
r.georob[["fitted.values"]] <- r.georob[["fitted.values"]] + offset
##
r.georob[["control"]] <- control
## add remaining items to output
if( control[["lmrob"]][["compute.rd"]] && !is.null( x ) )
r.georob[["MD"]] <- robMD( x, attr(mt, "intercept") )
if( model ) r.georob[["model"]] <- mf
if( ret.x ) r.georob[["x"]] <- x
if( ret.y ) r.georob[["y"]] <- y
r.georob[["df.residual"]] <- length(yy) - col.rank.XX[["rank"]]
r.georob[["na.action"]] <- attr(mf, "na.action")
r.georob[["offset"]] <- offset
r.georob[["contrasts"]] <- attr(x, "contrasts")
r.georob[["xlevels"]] <- .getXlevels(mt, mf)
r.georob[["rank"]] <- col.rank.XX[["rank"]]
r.georob[["call"]] <- cl
r.georob[["terms"]] <- mt
## set missing names of coefficients (bug of update)
if( identical(length( r.georob[["coefficients"]] ), 1L) &&
is.null( names( r.georob[["coefficients"]] ) ) ){
names( r.georob[["coefficients"]] ) <- "(Intercept)"
}
class( r.georob ) <- c( "georob" )
invisible( r.georob )
}
## ##############################################################################
pmm <-
function(
A, B, control = control.pcmp()
)
{
## function for parallelized matrix multiplication inspired by function
## parMM{snow}
## 2014-06-25 A. Papritz
## 2015-03-13 AP small changes in f.aux
## 2016-07-20 AP arguments renamed, new control of recursive parallelization
## 2018-01-05 AP improved memory management in parallel computations
## 2023-12-20 AP added on.exit(options(old.opt)), deleted options(error = NULL)
## auxiliary function
f.aux <- function( i ){
## s, e, A, B are taken from parent environment
A %*% B[ , s[i]:e[i], drop = FALSE]
}
## determine number of cores
ncores <- control[["pmm.ncores"]]
## determine columns indices of matrix blocks
k <- control[["f"]] * ncores
n <- NCOL(B)
dn <- floor( n / k )
s <- ( (0L:(k-1L)) * dn ) + 1L
e <- (1L:k) * dn
e[k] <- n
## set default value for control of forking if missing (required for backward compatibility)
if( is.null( control[["fork"]] ) ){
control[["fork"]] <- !identical( .Platform[["OS.type"]], "windows" )
}
if( ncores > 1L ){
if( !control[["fork"]] ){
## create a SNOW cluster on windows OS
if( !sfIsRunning() ){
old.opt <- options( error = f.stop.cluster )
on.exit( options( old.opt ) )
junk <- sfInit( parallel = TRUE, cpus = ncores )
}
junk <- sfExport( "s", "e", "A", "B" )
res <- sfLapply( 1L:k, f.aux )
if( control[["sfstop"]] ){
junk <- sfStop()
}
} else {
res <- mclapply( 1L:k, f.aux, mc.cores = ncores )
}
junk <- gc()
matrix( unlist(res), nrow = NROW( A ) )
} else A %*% B
}
# ##############################################################################
control.georob <-
function(
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(),
...
)
{
## auxiliary function to set meaningful default values for georob
## 2012-04-21 A. Papritz
## 2012-05-03 AP bounds for safe parameter values
## 2012-05-04 AP modifications for lognormal block kriging
## 2013-04-23 AP new names for robustness weights
## 2013-06-12 AP changes in stored items of Valphaxi object
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
## 2014-05-15 AP changes for version 3 of RandomFields
## 2014-08-18 AP changes for parallelized computations
## 2014-08-18 AP changes for Gaussian ML estimation
## 2015-03-10 AP extended variogram parameter transformations
## 2015-06-30 AP function and arguments renamed
## 2015-07-17 AP nlminb optimizer added
## 2015-08-19 AP control about error families for computing covariances added
## 2015-08-28 AP control arguments hessian, initial.param, initial.fixef newly organized
## 2016-07-15 AP some arguments hard coded
## 2016-08-24 AP new default value for error.family.cov.residuals
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## set values of fixed control arguments
irwls.initial <- TRUE
tuning.psi.nr <- 1000.
cov.bhat.betahat <- FALSE
aux.cov.pred.target <- FALSE
safe.param <- 1.e12
sepstr <- ".__...__."
## match arguments
ml.method <- match.arg( ml.method )
maximizer <- match.arg( maximizer )
initial.fixef <- match.arg( initial.fixef )
psi.func <- match.arg( psi.func )
error.family.estimation <- match.arg( error.family.estimation )
error.family.cov.effects <- match.arg( error.family.cov.effects )
error.family.cov.residuals <- match.arg( error.family.cov.residuals )
## check arguments
stopifnot(identical(length(reparam), 1L) && is.logical(reparam))
stopifnot(identical(length(initial.param), 1L) && is.logical(initial.param))
stopifnot(identical(length(force.gradient), 1L) && is.logical(force.gradient))
stopifnot(identical(length(cov.bhat), 1L) && is.logical(cov.bhat))
stopifnot(identical(length(full.cov.bhat), 1L) && is.logical(full.cov.bhat))
stopifnot(identical(length(cov.betahat), 1L) && is.logical(cov.betahat))
stopifnot(identical(length(cov.delta.bhat), 1L) && is.logical(cov.delta.bhat))
stopifnot(identical(length(full.cov.delta.bhat), 1L) && is.logical(full.cov.delta.bhat))
stopifnot(identical(length(cov.delta.bhat.betahat), 1L) && is.logical(cov.delta.bhat.betahat))
stopifnot(identical(length(cov.ehat), 1L) && is.logical(cov.ehat))
stopifnot(identical(length(full.cov.ehat), 1L) && is.logical(full.cov.ehat))
stopifnot(identical(length(cov.ehat.p.bhat), 1L) && is.logical(cov.ehat.p.bhat))
stopifnot(identical(length(full.cov.ehat.p.bhat), 1L) && is.logical(full.cov.ehat.p.bhat))
stopifnot(identical(length(hessian), 1L) && is.logical(hessian))
stopifnot(identical(length(min.rweight), 1L) && is.numeric(min.rweight) && min.rweight >= 0)
stopifnot(identical(length(irwls.maxiter), 1L) && is.numeric(irwls.maxiter) && irwls.maxiter >= 1)
stopifnot(identical(length(irwls.ftol), 1L) && is.numeric(irwls.ftol) && irwls.ftol > 0)
stopifnot(identical(length(min.condnum), 1L) && is.numeric(min.condnum) && min.condnum >= 0)
stopifnot(identical(length(zero.dist), 1L) && is.numeric(zero.dist) && zero.dist > 0)
stopifnot(is.null(bhat) || is.numeric(bhat))
stopifnot(is.list(param.tf))
stopifnot(is.list(fwd.tf))
stopifnot(is.list(deriv.fwd.tf))
stopifnot(is.list(bwd.tf))
stopifnot(is.list(rq))
stopifnot(is.list(lmrob))
stopifnot(is.list(nleqslv))
stopifnot(is.list(optim))
stopifnot(is.list(nlminb))
stopifnot(is.list(pcmp))
if(
!( all( unlist( param.tf ) %in% names( fwd.tf ) ) &&
all( unlist( param.tf ) %in% names( deriv.fwd.tf ) ) &&
all( unlist( param.tf ) %in% names( bwd.tf ) )
)
) stop(
"undefined transformation of variogram parameters; extend respective function definitions"
)
if( !irwls.initial && irwls.ftol >= 1.e-6 ) warning(
"'irwls.initial == FALSE' and large 'ftol' may create problems for root finding"
)
list(
ml.method = ml.method, reparam = reparam,
maximizer = maximizer,
initial.param = initial.param,
initial.fixef = initial.fixef,
bhat = bhat,
min.rweight = min.rweight,
param.tf = param.tf, fwd.tf = fwd.tf, deriv.fwd.tf = deriv.fwd.tf, bwd.tf = bwd.tf,
safe.param = safe.param, sepstr = sepstr,
psi.func = psi.func,
tuning.psi.nr = tuning.psi.nr,
irwls.initial = irwls.initial, irwls.maxiter = irwls.maxiter, irwls.ftol = irwls.ftol,
force.gradient = force.gradient,
min.condnum = min.condnum,
zero.dist = zero.dist,
error.family.estimation = error.family.estimation,
error.family.cov.effects = error.family.cov.effects,
error.family.cov.residuals = error.family.cov.residuals,
cov.bhat = cov.bhat, full.cov.bhat = full.cov.bhat,
cov.betahat = cov.betahat,
cov.bhat.betahat = cov.bhat.betahat,
cov.delta.bhat = cov.delta.bhat, full.cov.delta.bhat = full.cov.delta.bhat,
cov.delta.bhat.betahat = cov.delta.bhat.betahat,
cov.ehat = cov.ehat, full.cov.ehat = full.cov.ehat,
cov.ehat.p.bhat = cov.ehat.p.bhat, full.cov.ehat.p.bhat = full.cov.ehat.p.bhat,
aux.cov.pred.target = aux.cov.pred.target,
hessian = hessian,
irf.models = c( "RMdewijsian", "RMfbm", "RMgenfbm" ),
rq = rq, lmrob = lmrob, nleqslv = nleqslv,
optim = optim,
nlminb = nlminb,
pcmp = pcmp
)
}
## ======================================================================
param.transf <-
function(
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"
)
{
## function sets meaningful defaults for transformation of variogram
## parameters
## 2013-07-02 A. Papritz
## 2014-05-15 AP changes for version 3 of RandomFields
## 2015-03-10 AP extended transformation
## 2015-04-07 AP changes for fitting anisotropic variograms
list(
variance = variance, snugget = snugget, nugget = nugget, scale = scale,
alpha = alpha,
beta = beta,
delta = delta,
gamma = gamma,
kappa = kappa,
lambda = lambda,
mu = mu,
nu = nu,
f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta
)
}
## ======================================================================
fwd.transf <-
function(
...
)
{
## definition of forward transformation of variogram parameters
## 2013-07-02 A. Papritz
## 2015-03-10 AP extended variogram parameter transformations
## 2015-04-07 AP changes for fitting anisotropic variograms
list(
log = function(x) log(x),
logit1 = function(x) log( x / (1. - x) ),
logit2 = function(x) log( x / (2. - x) ),
logit3 = function(x) log( (x - 1.) / (3. - x) ),
identity = function(x) x, ...
)
}
## ======================================================================
dfwd.transf<-
function(
...
)
{
## definition of first derivative of forward transformation of variogram
## parameters
## 2013-07-02 A. Papritz
## 2015-03-10 AP extended variogram parameter transformations
## 2015-04-07 AP changes for fitting anisotropic variograms
list(
log = function(x) 1./x,
logit1 = function(x) 1. / (x - x^2),
logit2 = function(x) 2. / (2.*x - x^2),
logit3 = function(x) 2. / (4.*x - 3. - x^2),
identity = function(x) rep(1., length(x)), ...
)
}
## ======================================================================
bwd.transf <-
function(
...
)
{
## definition of backward transformation of variogram parameters
## 2013-07-02 A. Papritz
## 2015-03-10 AP extended variogram parameter transformations
## 2015-04-07 AP changes for fitting anisotropic variograms
## 2016-08-03 AP corrections for logitx if argument is +/-Inf
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
list(
log = function(x) exp(x),
logit1 = function(x){
stopifnot(identical(length(x), 1L) && is.numeric(x))
if( !is.finite(x) ){
if(sign(x) < 0.) 0. else 1.
} else exp(x) / (1. + exp(x))
},
logit2 = function(x){
stopifnot(identical(length(x), 1L) && is.numeric(x))
if( !is.finite(x) ){
if(sign(x) < 0.) 0. else 2.
} else 2. * exp(x) / (1. + exp(x))
},
logit3 = function(x){
stopifnot(identical(length(x), 1L) && is.numeric(x))
if( !is.finite(x) ){
if(sign(x) < 0.) 1. else 3.
} else (3. * exp(x) + 1.) / (1. + exp(x))
},
identity = function(x) x, ...
)
}
## ======================================================================
control.rq <-
function(
## specific arguments for rq: tau = 0.5, method = "br"
tau = 0.5, rq.method = c("br", "fnb", "pfn"),
## specific arguments for rq.fit.br: tau = 0.5, alpha = 0.1, ci = FALSE, iid = TRUE, interp = TRUE, tcrit = TRUE
rq.alpha = 0.1, ci = FALSE, iid = TRUE, interp = TRUE, tcrit = TRUE,
## specific arguments for rq.fit.fnb: tau = 0.5, rhs = (1 - tau) * apply(x, 2, sum), beta = 0.99995, eps = 1e-06
rq.beta = 0.99995, eps = 1.e-06,
## specific arguments for rq.fit.pfn: tau = 0.5, Mm.factor = 0.8, max.bad.fixup = 3, eps = 1e-06
Mm.factor = 0.8, max.bad.fixup = 3,
...
)
{
## function sets meaningful defaults for selected arguments of function
## rq{quantreg}
## 2012-12-14 A. Papritz
## 2014-07-29 AP
## 2015-06-30 AP function and arguments renamed
## 2017-05-09 AP small changes
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## match arguments
rq.method = match.arg(rq.method)
## check values of arguments
stopifnot(identical(length(ci), 1L) && is.logical(ci))
stopifnot(identical(length(iid), 1L) && is.logical(iid))
stopifnot(identical(length(interp), 1L) && is.logical(interp))
stopifnot(identical(length(tcrit), 1L) && is.logical(tcrit))
stopifnot(identical(length(tau), 1L) && is.numeric(tau) && tau > 0 && tau < 1)
stopifnot(identical(length(rq.alpha), 1L) && is.numeric(rq.alpha) && rq.alpha > 0 && rq.alpha < 1)
stopifnot(identical(length(rq.beta), 1L) && is.numeric(rq.beta))
stopifnot(identical(length(eps), 1L) && is.numeric(eps))
stopifnot(identical(length(Mm.factor), 1L) && is.numeric(Mm.factor) && Mm.factor > 0)
stopifnot(identical(length(max.bad.fixup), 1L) && is.numeric(max.bad.fixup) && max.bad.fixup > 0)
list(
tau = tau, method = rq.method,
alpha = rq.alpha, ci = ci, iid = iid, interp = interp, tcrit = tcrit,
beta = rq.beta, eps = eps,
Mm.factor = Mm.factor, max.bad.fixup = max.bad.fixup
)
}
## ======================================================================
control.nleqslv <-
function(
method = c( "Broyden", "Newton"),
global = c( "dbldog", "pwldog", "qline", "gline", "none" ),
xscalm = c( "fixed", "auto" ),
control = list( ftol = 1.e-4 ),
...
)
{
## function sets meaningful defaults for selected arguments of function
## nleqslv{nleqslv}
## 2013-07-12 A. Papritz
## 2014-07-29 AP
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(is.list(control))
list(
method = match.arg( method ),
global = match.arg( global ),
xscalm = match.arg( xscalm ),
control = control
)
}
## ======================================================================
control.optim <-
function(
method = c( "BFGS", "Nelder-Mead", "CG", "L-BFGS-B", "SANN", "Brent" ),
lower = -Inf, upper = Inf,
control = list(reltol = 1.e-5),
...
)
{
## function sets meaningful defaults for selected arguments of function optim
## 2012-12-14 A. Papritz
## 2014-07-29 AP
## 2015-06-30 AP function and arguments renamed
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(is.numeric(lower))
stopifnot(is.numeric(upper))
stopifnot(is.list(control))
list(
method = match.arg( method ),
lower = lower, upper = upper,
control = control
)
}
## ======================================================================
control.nlminb <-
function(
control = list( rel.tol = 1.e-5 ),
lower = -Inf, upper = Inf,
...
)
{
## function sets meaningful defaults for selected arguments of function nlminb
## 2015-07-17 A. Papritz
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(is.numeric(lower))
stopifnot(is.numeric(upper))
stopifnot(is.list(control))
list(
control = control,
lower = lower, upper = upper
)
}
## ======================================================================
control.pcmp <-
function(
pmm.ncores = 1, gcr.ncores = 1, max.ncores = parallel::detectCores(),
f = 1, sfstop = FALSE, allow.recursive = FALSE,
fork = !identical( .Platform[["OS.type"]], "windows" ),
...
)
{
## function sets meaningful defaults for parallelized computations
## 2014-07-29 AP
## 2015-06-30 AP function and arguments renamed
## 2015-07-29 AP changes for elimination of parallelized computation of gradient or estimating equations
## 2016-07-20 AP renamed function, separate ncores arguments various parallelized computations
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## 2024-01-25 AP changed default for allow.recursive
## check arguments
stopifnot(identical(length(sfstop), 1L) && is.logical(sfstop))
stopifnot(identical(length(allow.recursive), 1L) && is.logical(allow.recursive))
stopifnot(identical(length(fork), 1L) && is.logical(fork))
stopifnot(identical(length(pmm.ncores), 1L) && is.numeric(pmm.ncores) && pmm.ncores >= 1)
stopifnot(identical(length(gcr.ncores), 1L) && is.numeric(gcr.ncores) && gcr.ncores >= 1)
stopifnot(identical(length(max.ncores), 1L) && is.numeric(max.ncores) && max.ncores >= 1)
stopifnot(identical(length(f), 1L) && is.numeric(f) && f >= 1)
pmm.ncores <- min( pmm.ncores, max.ncores )
gcr.ncores <- min( gcr.ncores, max.ncores )
list(
pmm.ncores = pmm.ncores, gcr.ncores = gcr.ncores, max.ncores = max.ncores,
f = f, sfstop = sfstop,
allow.recursive = allow.recursive,
fork = fork
)
}
## ======================================================================
compress <-
function( m )
{
## function stores a list of or a single lower, upper triangular or
## symmetric matrix compactly
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
aux <- function( x ){
struc <- attr( x, "struc" )
if( !is.null( struc ) ){
stopifnot(identical(length(struc), 1L))
switch(
struc,
sym = {
aux <- list( diag = diag( x ), tri = x[lower.tri(x)] )
attr( aux, "struc" ) <- "sym"
aux
},
lt = {
aux <- list( diag = diag( x ), tri = x[lower.tri(x)] )
attr( aux, "struc" ) <- "lt"
aux
}
,
ut = {
aux <- list( diag = diag( x ), tri = x[upper.tri(x)] )
attr( aux, "struc" ) <- "ut"
aux
}
)
} else {
x
}
}
if( is.list( m ) ){
lapply( m, aux )
} else {
aux ( m )
}
}
## ======================================================================
expand <-
function( object )
{
## function expands a list of or a compactly stored lower, upper
## triangular or symmetric matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
aux <- function( x ){
struc <- attr( x, "struc" )
if( !is.null( struc ) ){
stopifnot(identical(length(struc), 1L))
switch(
struc,
sym = {
n <- length( x[["diag"]] )
dn <- names( x[["diag"]] )
aux <- matrix( 0., n, n )
aux[lower.tri( aux )] <- x[["tri"]]
aux <- aux + t( aux )
diag( aux ) <- x[["diag"]]
dimnames( aux ) <- list( dn, dn )
attr( aux, "struc" ) <- "sym"
aux
},
lt = {
n <- length( x[["diag"]] )
dn <- names( x[["diag"]] )
aux <- matrix( 0., n, n )
aux[lower.tri( aux )] <- x[["tri"]]
diag( aux ) <- x[["diag"]]
dimnames( aux ) <- list( dn, dn )
attr( aux, "struc" ) <- "lt"
aux
}
,
ut = {
n <- length( x[["diag"]] )
dn <- names( x[["diag"]] )
aux <- matrix( 0., n, n )
aux[upper.tri( aux )] <- x[["tri"]]
diag( aux ) <- x[["diag"]]
dimnames( aux ) <- list( dn, dn )
attr( aux, "struc" ) <- "ut"
aux
}
)
} else {
x
}
}
ln <- names( object )
if( is.list( object ) ){
if( length( ln ) == 2L && all( ln == c( "diag", "tri" ) ) ){
aux( object )
} else {
lapply( object, aux )
}
} else {
object
}
}
## ======================================================================
param.names <-
function( model )
{
## function returns names of extra parameters of implemented variogram
## models (cf. Variogram{RandomFields})
## 2012-01-24 A. Papritz
## 2014-05-15 AP changes for version 3 of RandomFields
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(identical(length(model), 1L) && is.character(model))
switch(
model,
"RMaskey" = "alpha",
"RMbessel" = "nu",
"RMcauchy" = "gamma",
"RMcircular" = NULL,
"RMcubic" = NULL,
"RMdagum" = c( "beta", "gamma" ),
"RMdampedcos" = "lambda",
"RMdewijsian" = "alpha",
"RMexp" = NULL,
"RMfbm" = "alpha",
"RMgauss" = NULL,
"RMgencauchy" = c( "alpha", "beta" ),
"RMgenfbm" = c( "alpha", "delta" ),
"RMgengneiting" = c( "kappa", "mu" ),
"RMgneiting" = NULL,
"RMlgd" = c( "alpha", "beta" ),
"RMmatern" = "nu",
"RMpenta" = NULL,
"RMqexp" = "alpha",
"RMspheric" = NULL,
"RMstable" = "alpha",
"RMwave" = NULL,
"RMwhittle" = "nu",
stop( model, " variogram not implemented" )
)
}
## ##############################################################################
param.bounds <-
function( model, d )
{
## function returns range of parameters for which variogram models are
## valid (cf. Variogram{RandomFields})
## 2012-03-30 A. Papritz
## 2014-05-15 AP changes for version 3 of RandomFields
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## 2023-11-17 AP correction of error for RMlgd (valid only for d in (1, 2)
## 2023-11-17 AP correction of lower bound for RMdampedcos
## 2023-12-20 AP correction of upper bounds of alpha for RMfbm and RMgenfbm
## check arguments
stopifnot(identical(length(model), 1L) && is.character(model))
stopifnot(identical(length(d), 1L) && is.numeric(d) && d >= 1)
deps <- .Machine$double.eps * 2.
switch(
model,
"RMaskey" = list( alpha = c( 0.5 * (d + 1L), Inf ) ),
"RMbessel" = list( nu = c( 0.5 * (d - 2L), Inf ) ),
"RMcauchy" = list( gamma = c( deps, Inf ) ),
"RMcircular" = NULL,
"RMcubic" = NULL,
"RMdagum" = list( beta = c( deps, 1.), gamma = c( deps, 1.-deps) ),
"RMdampedcos" = list( lambda = c( 1./tan(pi/(2 * d)), Inf ) ),
"RMdewijsian" = list( alpha = c( deps, 2. ) ),
"RMexp" = NULL,
"RMfbm" = list( alpha = c( deps, 2. - deps) ),
"RMgauss" = NULL,
"RMgencauchy" = list( alpha = c(deps, 2.), beta = c(deps, Inf) ),
"RMgenfbm" = list( alpha = c(deps, 2. - deps), delta = c(deps, 1. - deps) ),
"RMgengneiting" = list( kappa = c(1, 3), mu = c( d/2, Inf ) ),
"RMgneiting" = NULL,
"RMlgd" = list(
alpha = c(
deps,
if( d <= 2L ) 0.5 * (3L-d) else stop("dimension > 2 not allowed for RMlgd model" )
),
beta = c(deps, Inf)
),
"RMmatern" = list( nu = c(deps, Inf) ),
"RMpenta" = NULL,
"RMqexp" = list( alpha = c(0., 1.) ),
"RMspheric" = NULL,
"RMstable" = list( alpha = c(deps, 2.) ),
"RMwave" = NULL,
"RMwhittle" = list( nu = c(deps, Inf) ),
stop( model, " variogram not implemented" )
)
}
## ##############################################################################
default.fit.param <-
function(
variance = TRUE, snugget = FALSE, nugget = TRUE, scale = TRUE,
alpha = FALSE, beta = FALSE, delta = FALSE, gamma = FALSE,
kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE )
{
## function sets default flags for fitting variogram parameters
## 2015-11-27 A. Papritz
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(identical(length(variance), 1L) && is.logical(variance))
stopifnot(identical(length(snugget), 1L) && is.logical(snugget))
stopifnot(identical(length(nugget), 1L) && is.logical(nugget))
stopifnot(identical(length(scale), 1L) && is.logical(scale))
stopifnot(identical(length(alpha), 1L) && is.logical(alpha))
stopifnot(identical(length(beta), 1L) && is.logical(beta))
stopifnot(identical(length(delta), 1L) && is.logical(delta))
stopifnot(identical(length(gamma), 1L) && is.logical(gamma))
stopifnot(identical(length(kappa), 1L) && is.logical(kappa))
stopifnot(identical(length(lambda), 1L) && is.logical(lambda))
stopifnot(identical(length(mu), 1L) && is.logical(mu))
stopifnot(identical(length(nu), 1L) && is.logical(nu))
c(
variance = variance, snugget = snugget, nugget = nugget, scale = scale,
alpha = alpha, beta = beta, delta = delta, gamma = gamma,
kappa = kappa, lambda = lambda, mu = mu, nu = nu
)
}
## ##############################################################################
default.fit.aniso <-
function( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
{
## function sets default flags for fitting anisotropy parameters
## 2015-11-27 A. Papritz
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(identical(length(f1), 1L) && is.logical(f1))
stopifnot(identical(length(f2), 1L) && is.logical(f2))
stopifnot(identical(length(omega), 1L) && is.logical(omega))
stopifnot(identical(length(phi), 1L) && is.logical(phi))
stopifnot(identical(length(zeta), 1L) && is.logical(zeta))
c( f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta )
}
## ##############################################################################
default.aniso <-
function(
f1 = 1., f2 = 1., omega = 90., phi = 90., zeta = 0. )
{
## function sets default values for anisotropy parameters
## 2015-11-26 A. Papritz
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## check arguments
stopifnot(identical(length(f1), 1L) && is.numeric(f1) && f1 > 0)
stopifnot(identical(length(f2), 1L) && is.numeric(f2) && f2 > 0)
stopifnot(identical(length(omega), 1L) && is.numeric(omega))
stopifnot(identical(length(phi), 1L) && is.numeric(phi))
stopifnot(identical(length(zeta), 1L) && is.numeric(zeta))
c( f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta )
}
## ##############################################################################
### profilelogLik
profilelogLik <- function( object, values, use.fitted = TRUE, verbose = 0,
ncores = min( parallel::detectCores(), NROW(values) ) ){
## function to compute (restricted) likelihood profile for a georob fit
## 2015-03-18 A. Papritz
## 2015-04-08 AP changes in returned results
## 2016-07-14 AP optimization
## 2016-07-20 AP changes for parallel computations
## 2016-07-28 AP returns gradient in results
## 2016-08-08 AP changes for nested variogra
## 2016-08-12 AP changes for nested variogram models
## 2017-12-22 AP improved memory management in parallel computations
## 2020-02-14 AP sanity checks of arguments and for if() and switch()
## 2023-12-20 AP checking class by inherits()
## 2023-12-20 AP added on.exit(options(old.opt)), replaced makeCluster by makePSOCKcluster
## 2023-12-20 AP replacement of identical(class(...), ...) by inherits(..., ...)
## 2024-02-01 AP saving SOCKcluster.RData to tempdir()
## 2024-02-21 AP conditionally prevent recursive parallelizations in pmm or f.aux.gcr
#### -- auxiliary function
## auxiliary function to fit model and return maximized (pseudo) log-likelihood
f.aux <- function( i ){
## values, object, data are taken from parent environment
## set fixed initial values
values <- values[i, ]
tmp <- strsplit( names(values), control.georob()[["sepstr"]], fixed = TRUE )
cl <- object[["call"]]
for( i in 1L:length(values) ){
cl <- f.call.set_onexxx_to_value( cl, tmp[[i]][1L], values[i], as.integer(tmp[[i]][2L]) )
}
object[["call"]] <- cl
fit <- update( object, data = data )
## extract fitted variogram parameters
param.aniso <- unlist(lapply(
1L:length(fit[["variogram.object"]]),
function( i, x ){
x <- x[[i]]
res <- c( x[["param"]][x[["fit.param"]]], x[["aniso"]][x[["fit.aniso"]]] )
if( length(res) ){
names(res) <- paste( names(res), i, sep = control.georob()[["sepstr"]] )
res
} else {
NULL
}
}, x = fit[["variogram.object"]]
))
## results
c(
loglik = logLik(
fit, warn = FALSE, REML = identical( object[["control"]][["ml.method"]], "REML" )
),
param.aniso,
coef( fit ),
gradient = fit[["gradient"]],
converged = fit[["converged"]]
)
}
#### -- check arguments
## check whether all mandatory arguments have been provided
if( missing(object) || missing(values) ) stop(
"some mandatory arguments are missing"
)
stopifnot(inherits(object, "georob"))
stopifnot(identical(length(use.fitted), 1L) && is.logical(use.fitted))
stopifnot(identical(length(verbose), 1L) && is.numeric(verbose) && verbose >= 0)
stopifnot(identical(length(ncores), 1L) && is.numeric(ncores) && ncores >= 1)
## warning for robust fits
if( object[["tuning.psi"]] < object[["control"]][["tuning.psi.nr"]] ){
warning(
"likelihood approximated for robustly fitted model by likelihood of\n",
" equivalent Gaussian model with heteroscedastic nugget"
)
}
if( !(is.matrix(values) || is.data.frame( values ) || is.list( values ) ) ) stop(
"'values' must be a dataframe or a matrix"
)
#### -- prepare data
## get data.frame with required variables (note that the data.frame passed
## as data argument to georob must exist in GlobalEnv)
data <- cbind(
get_all_vars(
formula( object ), data = eval( getCall(object)[["data"]] )
),
get_all_vars(
object[["locations.objects"]][["locations"]], eval( getCall(object)[["data"]] )
)
)
if( inherits( object[["na.action"]], "omit" ) ) data <- na.omit(data)
## select subset if appropriate
if( !is.null( getCall(object)[["subset"]] ) ){
data <- data[eval( getCall(object)[["subset"]] ), ]
}
## extract variogram.object
## check values
if( is.data.frame(values) ) values <- list( values )
if( !identical( length(values), length(object[["variogram.object"]]) ) ) stop(
"'values' must be a list of the same length as 'variogram.object'"
)
## check names of fixed variogram parameters and fix respective
## parameters in variogram.object
values <- lapply(
1L:length(values),
function( i, v, vo ){
v <- v[[i]]
vo <- vo[[i]]
fixed.param.aniso <- colnames( v )
## match fixed.param.aniso
fixed.param.aniso <- sapply(
fixed.param.aniso, match.arg,
choices = c( names(default.fit.param()), names(default.fit.aniso()) )
)
if(
any( !fixed.param.aniso %in% c( names( vo[["param"]] ), names( vo[["aniso"]] ) ) )
) stop( "column names of 'values' do not match names of variogram parameters" )
colnames(v) <- fixed.param.aniso
v
}, v = values, vo = object[["variogram.object"]]
)
#### -- manipulate call to compute only required items
## manipulate call so that fitted values are used as initial values
cl <- f.call.set_allxxx_to_fitted_values( object )
## fix initial values for parameters present in values
for( i in 1L:length(values) ){
for( nme in colnames(values[[i]]) ){
cl <- f.call.set_onefitxxx_to_value( cl, nme, FALSE, i )
}
}
## update object call to avoid computation of covariance matrices and
## hessian and to set reparam = FALSE if variance parameters are fitted
cl <- f.call.set_x_to_value( cl, "verbose", verbose )
reparam <- !any(
unique( sapply( values, colnames ) ) %in% c( "variance", "snugget", "nugget" )
)
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "reparam", reparam )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.bhat", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.betahat", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.delta.bhat", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.delta.bhat.betahat", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.ehat", FALSE )
cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.ehat.p.bhat", FALSE )
## conditionally prevent recursive parallelizations in pmm or f.aux.gcr
if( ncores > 1L && !object[["control"]][["pcmp"]][["allow.recursive"]] ){
cl <- f.call.prevent_recursive_parallelization( cl )
}
## update call in object
object[["call"]] <- cl
#### -- prepare sets of parameter values for which to compute likelihood
## rename variogram parameters in values
values <- lapply(
1L:length(values),
function( i, v, vo ){
v <- v[[i]]
names( v ) <- paste( names(v), i, sep = control.georob()[["sepstr"]] )
v
}, v = values
)
## expand values to data frame
if( length(values) > 1L ){
tmp <- values
values <- as.list( tmp[[1L]] )
for( i in 2L:length(tmp) ){
for( j in colnames( tmp[[i]] ) ){
values <- c( values, as.list( tmp[[i]][, j, drop=FALSE] ) )
}
}
values <- expand.grid( values )
} else {
values <- values[[1L]]
}
#### -- fit model for sets of parameter values for which to compute likelihood
## loop over all elements of values
values <- as.matrix( values )
## set default value for control of forking if missing (required for backward compatibility)
if( is.null( object[["control"]][["pcmp"]][["fork"]] ) ){
object[["control"]][["pcmp"]][["fork"]] <- !identical( .Platform[["OS.type"]], "windows" )
}
if( ncores > 1L && !object[["control"]][["pcmp"]][["fork"]] ){
## create a PSOCK cluster on windows OS
fname <- file.path( tempdir(), "SOCKcluster.RData" )
clstr <- makePSOCKcluster( ncores )
save( clstr, file = fname )
old.opt <- options( error = f.stop.cluster )
on.exit( options( old.opt ) )
## export required items to workers
junk <- clusterEvalQ( clstr, require( georob, quietly = TRUE ) )
junk <- clusterExport( clstr, c( "values", "object", "data" ), envir = environment() )
result <- parLapply(
clstr,
1L:NROW(values),
f.aux
)
f.stop.cluster( clstr, fname )
} else {
## fork child processes on non-windows OS
result <- mclapply(
1L:NROW(values),
f.aux,
mc.cores = ncores,
mc.allow.recursive = object[["control"]][["pcmp"]][["allow.recursive"]]
)
}
#### -- prepare output
## collect results
result <- as.data.frame( cbind( values, t( simplify2array( result ) ) ) )
if( identical( length( object[["variogram.object"]]), 1L ) ){
tmp <- colnames(result)
tmp <- gsub(
paste( control.georob()[["sepstr"]], "1", sep = "" ),
"", tmp
)
colnames(result) <- tmp
}
result
}
## ##############################################################################
### gencorr
gencorr <- function( x, variogram.model, param ){
## function computes generalized correlation (= negative semivariance) of
## some stationary and IRF autocorrelation models available in the
## RandomFields package. the function is used as a substitute for
## RFfctn{RandomFields}
## arguments:
## x vector with (possibly rotated) scaled lag distances
## variogram.model name of autocorrelation model
## param extra parameter of autocorrelation model
## function returns vector with negative semivariances
## 2023-11-17 Andreas Papritz
## 2024-02-01 AP moved to code file for exported functions
result <- switch(
variogram.model[1],
#### --- RMaskey (compact support)
RMaskey = {
A <- unname( param["alpha"] )
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- (1 - xx)^A - 1.
res
},
#### --- RMbessel
RMbessel = {
A <- unname( param["nu"] )
res <- rep( 0., length( x ) )
sel <- x > 0.
xx <- x[sel]
res[sel] <- ( 2.^A * gamma(1.+A) * xx^(-A) * besselJ( xx, A ) ) - 1.
res
},
#### --- RMcauchy
RMcauchy = {
A <- unname( param["gamma"] )
( 1. + x^2 )^(-A) - 1.
},
#### --- RMcircular (compact support)
RMcircular = {
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- -2./pi * ( xx * sqrt( 1. - xx^2 ) + asin( xx ) )
res
},
#### --- RMcubic (compact support)
RMcubic = {
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- -7*xx^2 + 8.75*xx^3 - 3.5*xx^5 + 0.75*xx^7
res
},
#### --- RMdagum
RMdagum = {
A <- unname( param["beta"] )
B <- unname( param["gamma"] )
-( 1. + x^(-A) )^(-B/A)
},
#### --- RMdampedcos
RMdampedcos = {
A <- unname( param["lambda"] )
exp( -A * x ) * cos( x ) - 1.
},
#### --- RMdewijsian (IRF0 model)
RMdewijsian = {
A <- unname( param["alpha"] )
-log( x^A + 1. )
},
#### --- RMexp
RMexp = {
exp( -x ) - 1.
},
#### --- RMfbm (IRF0 model)
RMfbm = {
A <- unname( param["alpha"] )
-x^A
},
#### --- RMgauss
RMgauss = {
exp( -x^2 ) - 1.
},
#### --- RMgencauchy
RMgencauchy = {
A <- unname( param["alpha"] )
B <- unname( param["beta"] )
( 1. + x^A )^(-B/A) - 1.
},
#### --- RMgenfbm (IRF0 model)
RMgenfbm = {
A <- unname( param["alpha"] )
B <- unname( param["delta"] )
-( ( x^A + 1)^B - 1.)
},
#### --- RMgengneiting (compact support)
RMgengneiting = {
A <- unname( param["kappa"] )
B <- unname( param["mu"] )
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- if( identical( as.integer(A), 1L ) ){
BB <- B + 2.5
( 1. + BB * xx ) * ( 1. - xx)^BB - 1.
} else if( identical( as.integer(A), 2L ) ){
BB <- B + 4.5
( 1. + BB * xx + (BB^2 - 1.) * xx^2 / 3. ) * (1. - xx)^BB - 1.
} else if( identical( as.integer(A), 3L ) ){
BB <- B + 6.5
(
1 + BB * xx + (2. * BB^2 - 3.) * xx^2 / 5. +
(BB^2 - 4.) * BB * xx^3 / 15.
) * (1. - xx)^BB - 1.
} else {
stop( "RMgengneiting model undefined for 'kappa' !%in% 1:3" )
}
res
},
#### --- RMgneiting (compact support)
RMgneiting = {
res <- rep( -1., length( x ) )
sel <- x < 1. / 0.301187465825
xx <- 0.301187465825 * x[sel]
res[sel] <- (1. + 8.* xx + 25. * xx^2 + 32. * xx^3 ) *
(1. - xx )^8 - 1.
res
},
#### --- RMlgd (cf Gneiting, T. / Schlather, M.
## Stochastic models that separate fractal dimension and the Hurst effect,
## 2004, SIAM Review , Vol. 46, No. 2, p. 269-282
RMlgd = {
A <- unname( param["alpha"] )
B <- unname( param["beta"] )
res <- rep( NA_real_, length( x ) )
sel <- x <= 1.
res[sel] <- -B / (A + B) * x[sel]^A
res[!sel] <- A / (A + B) * x[!sel]^(-B) - 1.
res
},
#### --- RMmatern
RMmatern = {
A <- unname( param["nu"] )
res <- rep( 0., length( x ) )
sel <- x > 0.
xx <- sqrt(2. * A) * x[sel]
res[sel] <- 2.^(1. - A) / gamma(A) * xx^A * besselK( xx, A ) - 1.
res
},
#### --- RMpenta (compact support)
RMpenta = {
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- -22./3. * xx^2 + 33. * xx^4 - 38.5 * xx^5 +
16.5 * xx^7 - 5.5 * xx^9 + 5./6. * xx^11
res
},
#### --- RMqexp
RMqexp = {
A <- unname( param["alpha"] )
( 2. * exp(-x) - A * exp(-2.*x) ) / (2. - A) - 1.
},
#### --- RMspheric (compact support)
RMspheric = {
res <- rep( -1., length( x ) )
sel <- x < 1.
xx <- x[sel]
res[sel] <- -1.5 * xx + 0.5 * xx^3
res
},
#### --- RMstable
RMstable = {
A <- unname( param["alpha"] )
exp( -x^A ) - 1.
},
#### --- RMwave
RMwave = {
res <- rep( 0., length( x ) )
sel <- x > 0.
xx <- x[sel]
res[sel] <- sin( xx ) / xx - 1.
res
},
#### --- RMwhittle
RMwhittle = {
A <- unname( param["nu"] )
res <- rep( 0., length( x ) )
sel <- x > 0.
xx <- x[sel]
res[sel] <- 2.^(1 - A) / gamma(A) * xx^A * besselK( xx, A ) - 1.
res
},
stop( variogram.model, " model undefined" )
)
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.