Nothing
R/georob.S3methods.R
In georob: Robust Geostatistical Analysis of Spatial Data
#####################################
# #
# Methoden fuer Klasse "georob" #
# #
#####################################
# add1.georob
# deviance.georob
# drop1.georob
# extractAIC.georob
# logLik.georob
# step.default
# step.georob
# model.frame.georob
# model.matrix.georob
# nobs.georob
# print.georob
# ranef.georob
# residuals.georob
# resid.georob
# rstandard.georob
# rstudent.georob
# summary.georob,
# print.summary.georob
# vcov.georob
# waldtest.georob
# cv.georob (in separatem File)
# rstudent.cv.georob
# summary.cv.georob
# print.summary.cv.georob
# 2011-08-11 A. Papritz
## ##############################################################################
model.frame.georob <-
function(
formula, ...
)
{
## model.frame method for class georob
## 2012-12-19 A. Papritz
class( formula ) <- "lm"
model.frame( formula, ... )
}
## ##############################################################################
model.matrix.georob <-
function(
object, ...
)
{
## model.matrix method for class georob
## 2012-12-19 A. Papritz
class( object ) <- "lm"
model.matrix( object, ... )
}
## ##############################################################################
nobs.georob <-
function(
object, ...
)
{
## nobs method for class georob
## 2012-12-19 A. Papritz
class( object ) <- "lm"
nobs( object, ... )
}
## ##############################################################################
print.georob <-
function(
x, digits = max(3, getOption("digits") - 3), ...
)
{
## Print method for class "georob".
## Arguments:
## x an object generated by f.georob.initial.guess
## digits number of digits
## 2011-08-13 A. Papritz
## 2012-02-07 AP change for anisotropic variograms
## 2012-12-18 AP invisible(x)
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-02 AP new transformation of rotation angles
## code borrowed from print.lmrob for printing fixed effects coeffcients
cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
cat("\nFixed effects coefficients:\n")
print(
format( coef( x), digits = digits ), print.gap = 2,
quote = FALSE
)
## print variogram parameters
cat("\n")
cat( "Variogram: ", x[["variogram.model"]], "\n" )
param <- x[["param"]]
names( param ) <- ifelse(
x[["initial.objects"]][["fit.param"]],
names( param ),
paste( names( param ), "(fixed)", sep = "" )
)
print(
format( param, digits = digits ), print.gap = 2,
quote = FALSE
)
## print anisotropy parameters
if( !x[["aniso"]][["isotropic"]] ){
cat("\n")
cat( "Anisotropy parameters: ", "\n" )
aniso <- x[["aniso"]][["aniso"]]
names( aniso ) <- ifelse(
x[["initial.objects"]][["fit.aniso"]],
names( aniso ),
paste( names( aniso ), "(fixed)", sep = "" )
)
print(
format( aniso, digits = digits ), print.gap = 2,
quote = FALSE
)
}
invisible( x )
}
## ##############################################################################
ranef.georob <- random.effects.georob <-
function(
object,
standard = FALSE,
...
)
{
## Function extracts the random effects (bhat) from georob object
## (similar to ranef.lme{nlme})
## Arguments:
## object fitted georob object
## standard an optional logical value indicating whether the estimated random effects
## should be "standardized" (i.e. divided by the estimated standard error.
## Defaults to FALSE.
## ... further arguments passed to method (currently not used)
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2012-01-05 AP modified for compress storage of matrices
## 2012-10-18 AP changes for new definition of eta
## 2012-11-26 AP method for random.effects
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-05-06 AP changes for solving estimating equations for xi
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## temporarily redefine na.action component of object
object.na <- object[["na.action"]]
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
Valpha.objects <- expand( object[["Valpha.objects"]] )
covmat <- expand( object[["cov"]] )
bhat <- object[["bhat"]]
if( standard ){
if( is.null( covmat[["cov.bhat"]] ) ){
## compute standard errors of residuals
X <- model.matrix(
terms( object ),
model.frame( object )
)[!duplicated( object[["Tmat"]] ), , drop = FALSE]
r.cov <- compute.covariances(
Valpha.objects = Valpha.objects,
Aalpha = object[["Aalpha"]],
Palpha = object[["Palpha"]],
rweights = object[["rweights"]],
XX = X, TT = object[["Tmat"]], names.yy = rownames( model.frame( object ) ),
nugget = object[["param"]]["nugget"],
eta = sum( object[["param"]][c( "variance", "snugget")] ) / object[["param"]]["nugget"],
expectations = object[["expectations"]],
cov.bhat = TRUE, full.cov.bhat = FALSE,
cov.betahat = FALSE,
cov.bhat.betahat = FALSE,
cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
cov.delta.bhat.betahat = FALSE,
cov.ehat = FALSE, full.cov.ehat = FALSE,
cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = FALSE,
verbose = 0
)
if( r.cov[["error"]] ) stop(
"an error occurred when computing the variances of the random effects"
)
se <- sqrt( r.cov[["cov.bhat"]] )
} else {
## extract standard errors of residuals from georob object
if( is.matrix( covmat[["cov.bhat"]] ) ){
se <- sqrt( diag( covmat[["cov.bhat"]] ) )
} else {
se <- sqrt( covmat[["cov.bhat"]] )
}
}
bhat <- bhat / se
}
naresid( object.na, bhat )
}
## ##############################################################################
fixef.georob <- fixed.effects.georob <-
function(
object,
...
)
{
## Function extracts residuals from georob object (based on residuals.lm {stats})
## Arguments:
## object fitted georob object
## type character, type of resdiuals to computed
## ... further arguments passed to methods
## 2012-11-26 A. Papritz
## 2013-06-12 AP substituting [["x"]] for $x in all lists
object[["coef"]]
}
## ##############################################################################
residuals.georob <- resid.georob <-
function(
object,
type = c("working", "response", "deviance", "pearson", "partial" ),
terms = NULL,
level = 1,
...
)
{
## Function extracts residuals from georob object (based on residuals.lm {stats})
## Arguments:
## object fitted georob object
## type character, type of resdiuals to computed
## level integer scalar to select whether ehat (level == 1) or
## ehat + bhat (level == 0) is returned,
## only effective for type %in% c( "working", "response", "partial" )
## ... further arguments passed to methods
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2013-05-31 AP modified for computing partial residuals for single terms
## 2013-06-12 AP substituting [["x"]] for $x in all lists
type <- match.arg( type )
if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
## temporarily redefine na.action component of object
object.na <- object[["na.action"]]
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
r <- object[["residuals"]]
res <- switch(
type,
working = ,
response = r,
deviance = ,
pearson = if( is.null(object[["weights"]]) ) r else r * sqrt(object[["weights"]]),
partial = r
)
if( level == 0 && any( type %in% c( "working", "response", "partial" ) ) ){
res <- res + ranef( object, standard = FALSE )[object[["Tmat"]]]
}
if( type == "partial" )
res <- res + predict( object, type = "terms", terms = terms )[["fit"]]
drop( res )
naresid( object.na, res )
}
## ##############################################################################
rstandard.georob <-
function( model, level = 1, ... )
{
## Function extracts standardized residuals from georob object
## Arguments:
## model fitted georob object
## level integer scalar to select whether ehat (level == 1) or
## ehat + bhat (level == 0) is returned,
## ... further arguments (currently not used)
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2012-01-05 AP modified for compress storage of matrices
## 2012-10-18 AP changes for new definition of eta
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-05-06 AP changes for solving estimating equations for xi
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## temporarily redefine na.action component of model
model.na <- model[["na.action"]]
if( identical( class( model[["na.action"]] ), "exclude" ) ){
class( model[["na.action"]] ) <- "omit"
}
Valpha.objects <- expand( model[["Valpha.objects"]] )
covmat <- expand( model[["cov"]] )
if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
if(
( is.null( covmat[["cov.ehat"]] ) & level == 1 ) ||
( is.null( covmat[["cov.ehat.p.bhat"]] ) & level == 0 )
){
## compute standard errors of residuals
X <- model.matrix(
terms( model),
model.frame( model )
)[!duplicated( model[["Tmat"]] ), , drop = FALSE]
r.cov <- compute.covariances(
Valpha.objects = Valpha.objects,
Aalpha = model[["Aalpha"]],
Palpha = model[["Palpha"]],
rweights = model[["rweights"]],
XX = X, TT = model[["Tmat"]], names.yy = rownames( model.frame( model ) ),
nugget = model[["param"]]["nugget"],
eta = sum( model[["param"]][c( "variance", "snugget")] ) / model[["param"]]["nugget"],
expectations = model[["expectations"]],
cov.bhat = FALSE, full.cov.bhat = FALSE,
cov.betahat = FALSE,
cov.bhat.betahat = FALSE,
cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
cov.delta.bhat.betahat = FALSE,
cov.ehat = level == 1, full.cov.ehat = FALSE,
cov.ehat.p.bhat = level == 0, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = FALSE,
verbose = 0
)
if( r.cov[["error"]] ) stop(
"an error occurred when computing the variance of the residuals"
)
if( level == 1 ){
covmat[["cov.ehat"]] <-r.cov[["cov.ehat"]]
} else {
covmat[["cov.ehat.p.bhat"]] <-r.cov[["cov.ehat.p.bhat"]]
}
}
## extract standard errors of residuals from georob object
if( level == 1 ){
se <- covmat[["cov.ehat"]]
} else {
se <- covmat[["cov.ehat.p.bhat"]]
}
if( is.matrix( se ) ){
se <- sqrt( diag( se ) )
} else {
se <- sqrt( se )
}
## compute standardized residuals
naresid( model.na, residuals( model, level = level ) / se )
}
## ##############################################################################
rstudent.georob <-
function( model, ... )
{
## Function computes studentized residuals for fitted georob object
## Arguments:
## model fitted georob object
## data data frame that was used to fit model
## ... further arguments passed to cv.georob
## 2011-12-22 A. Papritz
## 2013-06-12 AP substituting [["x"]] for $x in all lists
if( !identical( class( model )[1], "georob" ) ) stop(
"model is not of class 'georob'"
)
r.cv <- cv( model, ... )
rstudent( model = r.cv )
}
## ##############################################################################
summary.georob <-
function (
object, correlation = FALSE,
signif = 0.95,
...
)
{
## ToDos:
## - Terms Objekt einfuegen
## - ausgewaehlte Angaben zur Fitmethode ausgeben
## - Wald-Test des Modells y ~ 1
## 2012-01-05 A. Papritz
## 2012-01-05 AP modified for compress storage of matrices
## 2012-01-31 AP corretion of error for computing CI for variogram parameters
## 2012-02-07 AP change for anisotropic variograms
## 2012-03-29 AP warning messages inserted
## 2012-05-23 ap correction of error for computing correlation matrix of variogram parameters
## 2012-11-04 AP handling compressed cov.betahat
## 2012-11-27 AP changes in parameter back-transformation
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-03 AP new transformation of rotation angles
covmat <- expand( object[["cov"]] )
ans <- object[c(
"call", "residuals", "bhat", "rweights", "converged", "convergence.code",
"iter", "loglik", "variogram.model", "gradient",
"tuning.psi", "df.residual"
)]
ans <- c( ans, object[["initial.objects"]]["fit.param"] )
if( !object[["aniso"]][["isotropic"]] ) ans[["fit.param"]] <- c(
ans[["fit.param"]], object[["initial.objects"]][["fit.aniso"]]
)
ans[["terms"]] <- NA
ans[["scale"]] <- sqrt(object[["param"]]["nugget"])
ans[["control"]][["method"]] <- "TODO: PRINT GLSROB CONTROL PARAMETERS HERE"
se <- sqrt(diag(covmat[["cov.betahat"]]))
est <- object[["coefficients"]]
tval <- est/se
ans[["coefficients"]] <- cbind(
est, se, tval, 2 * pt( abs(tval), object[["df.residual"]], lower.tail = FALSE )
)
dimnames( ans[["coefficients"]] ) <- list(
names(est), c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
)
if( correlation ){
ans[["correlation"]] <- covmat[["cov.betahat"]] / outer(se, se)
}
ans[["param"]] <- as.matrix( object[["param"]], ncol = 1 )
if( !object[["aniso"]][["isotropic"]] ) ans[["param"]] <- rbind(
ans[["param"]],
as.matrix( object[["aniso"]][["aniso"]], ncol = 1 )
)
colnames( ans[["param"]] ) <- "Estimate"
## compute confidence intervals of variogram parameters from observed
## Fisher information matrix (Gaussian REML only)
if( !is.null( object[["hessian"]] ) ){
## initialization
cor.tf.param <- cov.tf.param <- matrix(
NA, nrow = nrow( object[["hessian"]] ), ncol = nrow( object[["hessian"]] ),
dimnames = dimnames( object[["hessian"]] )
)
se <- rep( NA, nrow( object[["hessian"]] ) )
names( se ) <- rownames( object[["hessian"]])
ci <- matrix( NA, nrow = nrow( ans[["param"]] ), ncol = 2 )
colnames( ci ) <- c( "Lower", "Upper" )
rownames( ci ) <- rownames( ans[["param"]] )
## select parameters that are not on boundary of parameter space
sr <- !apply( object[["hessian"]], 1, function( x ) all( is.na( x ) ) )
if( sum( sr ) > 0 ){
t.chol <- try( chol( object[["hessian"]][sr, sr] ), silent = TRUE )
if( !identical( class( t.chol ), "try-error" ) ){
## compute covariance matrix of fitted transformed parameters
cov.tf.param[sr, sr] <- chol2inv( t.chol )
## correlation matrix and standard errors of fitted transformed
## parameters
if( sum( sr ) > 1 ){
cor.tf.param[sr, sr] <- cov2cor( cov.tf.param[sr, sr] )
} else {
cor.tf.param[sr, sr] <- 1.
}
se[sr] <- sqrt( diag( cov.tf.param )[sr] )
## compute confidence interval on original scale of parameters
sel.names <- names( object[["param"]][object[["initial.objects"]][["fit.param"]]] )
if( !object[["aniso"]][["isotropic"]] ) sel.names <- c(
sel.names,
names( object[["aniso"]][["aniso"]][object[["initial.objects"]][["fit.aniso"]]] )
)
sel.names <- sel.names[sr]
ci[sel.names, ] <- t(
sapply(
sel.names,
function( x, param, se, param.tf, trafo.fct, inv.trafo.fct ){
inv.trafo.fct[[param.tf[x]]](
trafo.fct[[param.tf[x]]]( param[x] ) +
c(-1, 1) * se[x] * qnorm( (1-signif)/2, lower.tail = FALSE )
)
},
param = c( object[["param"]], object[["aniso"]][["aniso"]] ),
se = se,
param.tf = object[["param.tf"]],
trafo.fct = object[["fwd.tf"]],
inv.trafo.fct = object[["bwd.tf"]]
)
)
if( !object[["aniso"]][["isotropic"]] ){
## map angles to halfcircle
if( !object[["aniso"]][["isotropic"]] ){
sel <- match( "omega", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < 0., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 180., ci[sel, ] - 180., ci[sel, ] )
}
sel <- match( "phi", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < 0., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 180., ci[sel, ] - 180., ci[sel, ] )
}
sel <- match( "zeta", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < -90., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 90., ci[sel, ] - 180., ci[sel, ] )
}
}
}
ans[["param"]] <- cbind( ans[["param"]], ci )
if( correlation ) ans[["cor.tf.param"]] <- cor.tf.param
} else {
warning(
"Hessian not positive definite:",
"\nconfidence intervals of variogram parameters cannot be computed"
)
}
}
}
ans[["se.residuals"]] <- if( !is.null( covmat[["cov.ehat"]] ) ){
if( is.matrix( covmat[["cov.ehat"]] ) ){
sqrt( diag( covmat[["cov.ehat"]] ) )
} else {
sqrt( covmat[["cov.ehat"]] )
}
} else NULL
class( ans ) <- c( "summary.georob" )
ans
}
## ##############################################################################
print.summary.georob <-
function (
x, digits = max(3, getOption("digits") - 3),
signif.stars = getOption("show.signif.stars"),
...
)
{
## ToDos:
## - Ausgabe df
## - Wald-Test des Modells y ~ 1
## - ausgewaehlte Angaben zur Fitmethode ausgeben
## 2012-01-05 A. Papritz
## 2012-01-31 AP small changes
## 2012-02-07 AP change for anisotropic variograms
## 2012-12-18 AP invisible(x)
## 2013-04-23 AP new names for robustness weights
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-01-23 AP prints maximized loglik even if all parameters fixed
cat("\nCall:")
cat( paste( deparse(x[["call"]]), sep = "\n", collapse = "\n"), "\n", sep = "" )
cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
if( is.na( x[["converged"]] ) ){
cat( "\nEstimation with fixed variogram parameters\n" )
} else {
if(!(x[["converged"]])) {
cat(
"\nAlgorithm did not converge, diagnostic code: ",
x[["convergence.code"]], "\n"
)
} else {
cat(
"\nConvergence in", x[["iter"]][1], "function and",
x[["iter"]][2], "Jacobian/gradient evaluations\n"
)
}
attr( x[["gradient"]], "eeq.emp" ) <- NULL
attr( x[["gradient"]], "eeq.exp" ) <- NULL
cat( "\nEstimating equations (gradient)\n")
print( x[["gradient"]], digits = digits, ... )
}
if( x[["tuning.psi"]] >=
georob.control()[["tuning.psi.nr"]] ) cat(
"\nMaximized restricted log-likelihood:",
x[["loglik"]], "\n"
)
df <- x[["df.residual"]]
bhat <- x[["bhat"]]
cat( "\nPredicted latent variable (z):\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(bhat), names = nam )
print( rq, digits = digits, ...)
}
else print( bhat, digits = digits, ...)
resid <- x[["residuals"]]
cat( "\nResiduals (epsilon):\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(resid), names = nam )
print( rq, digits = digits, ...)
}
else print( resid, digits = digits, ...)
if( !is.null( x[["se.residuals"]] ) ){
resid <- x[["residuals"]] / x[["se.residuals"]]
cat( "\nStandardized residuals:\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(resid), names = nam )
print( rq, digits = digits, ...)
}
else print( resid, digits = digits, ...)
}
cat( "\nVariogram: ", x[["variogram.model"]], "\n" )
rownames( x[["param"]] ) <- ifelse(
x[["fit.param"]],
rownames( x[["param"]] ),
paste( rownames( x[["param"]] ), "(fixed)", sep = "" )
)
## print( format( x[["param"]], digits = digits ), print.gap = 2, quote = FALSE )
printCoefmat(
x[["param"]], digits = digits, signif.stars = FALSE, ...
)
if( !is.null( x[["cor.tf.param"]] ) ){
correl <- x[["cor.tf.param"]]
p <- NCOL(correl)
if( p > 1 ){
cat("\nCorrelation of (transformed) variogram parameters:\n")
correl <- format(round(correl, 2), nsmall = 2,
digits = digits)
correl[!lower.tri(correl)] <- ""
print(correl[-1, -p, drop = FALSE], quote = FALSE)
}
}
cat( "\nFixed effects coefficients:\n" )
printCoefmat(
x[["coefficients"]], digits = digits, signif.stars = signif.stars, ...
)
cat(
"\nResidual standard error (sqrt(nugget)):",
format(signif(x[["scale"]], digits)), "\n"
)
correl <- x[["correlation"]]
if( !is.null(correl) ){
p <- NCOL(correl)
if( p > 1 ){
cat("\nCorrelation of fixed effects coefficients:\n")
correl <- format(round(correl, 2), nsmall = 2,
digits = digits)
correl[!lower.tri(correl)] <- ""
print(correl[-1, -p, drop = FALSE], quote = FALSE)
}
}
cat("\n")
summarizeRobWeights(x[["rweights"]], digits = digits, ... )
invisible( x )
}
## ##############################################################################
vcov.georob <-
function( object, ... )
{
## 2012-11-04 AP handling compressed cov.betahat
## 2013-06-12 AP substituting [["x"]] for $x in all lists
result <- expand( object[["cov"]][["cov.betahat"]] )
attr( result, "struc" ) <- NULL
result
}
## ##############################################################################
waldtest.georob <-
function(
object, ..., vcov = NULL, test = c("Chisq", "F"), name = NULL,
fixed = TRUE
)
{
## waldtest method for class georob
## 2012-02-08 AP change for anisotropic variograms
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-02-19 AP correcting option for verbose output
## 2014-05-15 AP changes for version 3 of RandomFields
test <- match.arg( test )
## refit object with fixed variogram parameters
if( fixed ) {
cat( "\nWald-Test with fixed variogram parameters of model 1\n\n" )
object <- update(
object,
param = object[["param"]],
aniso = object[["aniso"]][["aniso"]],
fit.param = c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )],
fit.aniso = c(
f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
),
verbose = 0
)
}
## Wald-Test
waldtest.default(
object = object, ..., vcov = vcov, test = test, name = name
)
}
## ##############################################################################
logLik.georob <-
function( object, REML = FALSE, ... )
{
## 2012-12-22 method for extracting (restricted) loglikelihood
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-02-27 AP computing 'pseudo' likelihood for robustly fitted models
val <- if( REML ){
val <- object[["loglik"]]
} else {
D <- deviance( object, ... )
-0.5 * (
D + attr( D, "log.det.covmat" ) + length( object[["residuals"]] ) * log( 2 * pi )
)
}
attr(val, "nall") <- length(object[["residuals"]])
attr(val, "nobs") <- object[["df.residual"]]
attr(val, "df") <- length(object[["coefficients"]]) +
sum( object[["initial.objects"]][["fit.param"]] ) +
sum( object[["initial.objects"]][["fit.aniso"]])
class(val) <- "logLik"
val
}
## ##############################################################################
deviance.georob <-
function(
object, warn = TRUE, ...
)
{
## deviance method for class georob
## 2012-12-22 A. Papritz
## 2013-05-23 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-02-27 AP computing 'pseudo' deviance for robustly fitted models
## 2014-03-12 AP option for suppressing warnings
## redefine na.action component of object
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
## determine factor to compute heteroscedastic nugget
if( object[["tuning.psi"]] < georob.control()[["tuning.psi.nr"]] ){
if( warn ) warning(
"deviance for robustly fitted model approximated by deviance of\n",
" equivalent model with heteroscedastic nugget"
)
w <- object[["rweights"]]
} else {
w <- 1.
}
## invert covariance matrix of observations
Valpha.objects <- expand( object[["Valpha.objects"]] )
G <- sum( object[["param"]][c("variance", "snugget")] ) * Valpha.objects[["Valpha"]]
G <- G[object[["Tmat"]], object[["Tmat"]]]
diag( G ) <- diag( G ) + object[["param"]]["nugget"] / w
iucf <- try(
backsolve(
chol( G ),
diag( length( object[["Tmat"]] ) ),
k = length( object[["Tmat"]] )
),
silent = TRUE
)
## compute deviance
if( identical( class( iucf ), "try-error" ) ) {
stop( "(generalized) covariance matrix of observations not positive definite" )
} else {
result <- sum( colSums( residuals( object, level = 0 ) * iucf )^2 )
attr( result, "log.det.covmat" ) <- -2 * sum( log( diag( iucf ) ) )
}
result
}
## ##############################################################################
extractAIC.georob <- function( fit, scale = 0, k = 2, ... )
{
## extractAIC method for class georob
## 2014-02-27 AP
edf <- sum( !is.na( fit[["coefficients"]] ) ) +
sum( fit[["initial.objects"]][["fit.param"]] ) +
sum( fit[["initial.objects"]][["fit.aniso"]] )
loglik <- logLik( fit, ... )
c(edf, -2 * loglik + k * edf)
}
## ##############################################################################
safe_pchisq <- function(q, df, ...)
## same code as safe_pchisq{stats}
{
df[df <= 0] <- NA
pchisq(q=q, df=df, ...)
}
## ##############################################################################
add1.georob <- function( object, scope, scale = 0, test=c("none", "Chisq" ),
k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
ncores = 1, ... )
{
## add1 method for class georob based on add1.default{stats}
## 2014-03-12 AP
## 2014-04-24 AP function returns only variogram parameters of best fitted object
## 2014-05-15 AP changes for version 3 of RandomFields
## auxiliary function for fitting model and extracting aic
f.aux <- function(
tt, scale, k, trace,
fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
verbose, ...
){
converged <- TRUE
if(trace > 1) {
cat("\ntrying +", tt, "\n", sep = "")
utils::flush.console()
}
if( fixed ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . +", tt)),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . +", tt)), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
} else {
if( use.fitted.param ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . +", tt)),
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . +", tt)), data = data,
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
}
}
else {
if( is.null( data ) ){
nfit <- update( object, as.formula(paste("~ . +", tt)),
verbose = verbose
)
} else {
nfit <- update( object, as.formula(paste("~ . +", tt)), data = data,
verbose = verbose
)
}
}
if( !nfit$converged ){
converged <- FALSE
if( verbose > 0 ) cat(
"lack of convergence when fitting model", paste("~ . +", tt),
"\nconvergence code:", nfit$convergence.code, "\n"
)
}
}
nnew <- nobs( nfit, use.fallback = TRUE )
if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
"number of rows in use has changed: remove missing values?"
)
df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
names( df.aic ) <- c("df", "AIC", "converged")
attr( df.aic, "nfit" ) <- list(
param = nfit[["param"]],
aniso = nfit[["aniso"]][["aniso"]],
initial.objects = nfit[["initial.objects"]],
call = nfit[["call"]]
)
df.aic
}
## check scope
if( missing(scope ) || is.null( scope ) ) stop("no terms in scope")
if( !is.character( scope ) ) scope <- add.scope( object, update.formula(object, scope) )
if( !length(scope)) stop( "no terms in scope for adding to object" )
## initialize result
ns <- length( scope )
ans <- matrix(
nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
)
ans[1L, ] <- extractAIC( object, scale, k = k, ... )
## loop over all components of scope
n0 <- nobs( object, use.fallback = TRUE )
env <- environment( formula(object) )
param <- object[["param"]]
aniso <- object[["aniso"]][["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
## prepare for parallel execution
if( ncores > 1 ){
ncores <- min( ncores, ns, detectCores() )
trace <- 0
verbose <- 0
}
if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
if( is.null( data ) ) stop(
"argument 'data' required for parallel execution on windows OS"
)
## create a SNOW cluster on windows OS
cl <- makePSOCKcluster( ncores, outfile = "")
## export required items to workers
junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
result <- parLapply(
cl,
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
)
stopCluster(cl)
} else {
## fork child processes on non-windows OS
result <- mclapply(
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
mc.cores = ncores, mc.allow.recursive = FALSE, ...
)
}
fits <- list(
list(
param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
initial.objects = object[["initial.objects"]],
call = object[["call"]]
)
)
fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
result <- t( simplify2array( result ) )
ans[2:NROW(ans), ] <- result[, 1:2]
if(!all( as.logical(result[, "converged"]) ) ) warning(
"lack of convergence when fitting models"
)
## store results
dfs <- ans[, 1L] - ans[1L, 1L]
dfs[1L] <- NA
aod <- data.frame( Df = dfs, AIC = ans[, 2L] )
## likelihood ratio test
test <- match.arg(test)
if(test == "Chisq") {
dev <- ans[, 2L] - k*ans[, 1L]
dev <- dev[1L] - dev; dev[1L] <- NA
nas <- !is.na(dev)
P <- dev
P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail=FALSE)
aod[, c("LRT", "Pr(>Chi)")] <- list(dev, P)
}
## format results
head <- c(
"Single term additions", "\nModel:", deparse(formula(object)),
if(scale > 0) paste("\nscale: ", format(scale), "\n")
)
class(aod) <- c("anova", "data.frame")
attr( aod, "heading") <- head
attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
aod
}
## ##############################################################################
drop1.georob <- function( object, scope, scale = 0, test=c( "none", "Chisq" ),
k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
ncores = 1, ... )
{
## drop1 method for class georob based on drop1.default{stats}
## 2014-03-12 AP
## 2014-04-24 AP function returns only variogram parameters of best fitted object
## 2014-05-15 AP changes for version 3 of RandomFields
## auxiliary function for fitting model and extracting aic
f.aux <- function(
tt, scale, k, trace,
fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
verbose, ...
){
converged <- TRUE
if(trace > 1) {
cat("\ntrying -", tt, "\n", sep = "")
utils::flush.console()
}
if( fixed ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . -", tt)),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . -", tt)), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
} else {
if( use.fitted.param ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . -", tt)),
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . -", tt)), data = data,
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
}
}
else {
if( is.null( data ) ){
nfit <- update( object, as.formula(paste("~ . -", tt)),
verbose = verbose
)
} else {
nfit <- update( object, as.formula(paste("~ . -", tt)), data = data,
verbose = verbose
)
}
}
if( !nfit$converged ){
converged <- FALSE
if( verbose > 0 ) cat(
"lack of convergence when fitting model", paste("~ . -", tt),
"\nconvergence code:", nfit$convergence.code, "\n"
)
}
}
nnew <- nobs( nfit, use.fallback = TRUE )
if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
"number of rows in use has changed: remove missing values?"
)
df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
names( df.aic ) <- c("df", "AIC", "converged")
attr( df.aic, "nfit" ) <- list(
param = nfit[["param"]],
aniso = nfit[["aniso"]][["aniso"]],
initial.objects = nfit[["initial.objects"]],
call = nfit[["call"]]
)
df.aic
}
## check scope
tl <- attr( terms(object), "term.labels" )
if( missing(scope) ) {
scope <- drop.scope( object )
} else {
if( !is.character(scope) ) scope <- attr(
terms(update.formula(object, scope)), "term.labels"
)
if( !all( match(scope, tl, 0L) > 0L) ) stop("scope is not a subset of term labels")
}
## initialize result
ns <- length( scope )
ans <- matrix(
nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
)
ans[1, ] <- extractAIC( object, scale, k = k, ... )
## loop over all components of scope
n0 <- nobs( object, use.fallback = TRUE )
env <- environment( formula(object) )
param <- object[["param"]]
aniso <- object[["aniso"]][["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
## prepare for parallel execution
if( ncores > 1 ){
ncores <- min( ncores, ns, detectCores() )
trace <- 0
verbose <- 0
}
if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
if( is.null( data ) ) stop(
"argument 'data' required for parallel execution on windows OS"
)
## create a SNOW cluster on windows OS
cl <- makePSOCKcluster( ncores, outfile = "")
## export required items to workers
junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
result <- parLapply(
cl,
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
)
stopCluster(cl)
} else {
## fork child processes on non-windows OS
result <- mclapply(
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
mc.cores = ncores, mc.allow.recursive = FALSE, ...
)
}
fits <- list(
list(
param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
initial.objects = object[["initial.objects"]],
call = object[["call"]]
)
)
fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
result <- t( simplify2array( result ) )
ans[2:NROW(ans), ] <- result[, 1:2]
if(!all( as.logical(result[, "converged"]) ) ) warning(
"lack of convergence when fitting models"
)
## store results
dfs <- ans[1L , 1L] - ans[, 1L]
dfs[1L] <- NA
aod <- data.frame( Df = dfs, AIC = ans[,2] )
## likelihood ratio test
test <- match.arg(test)
if(test == "Chisq") {
dev <- ans[, 2L] - k*ans[, 1L]
dev <- dev - dev[1L] ; dev[1L] <- NA
nas <- !is.na(dev)
P <- dev
P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail = FALSE)
aod[, c("LRT", "Pr(>Chi)")] <- list(dev, P)
}
## format results
head <- c("Single term deletions", "\nModel:", deparse(formula(object)),
if(scale > 0) paste("\nscale: ", format(scale), "\n"))
class(aod) <- c("anova", "data.frame")
attr( aod, "heading") <- head
attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
aod
}
## ############################################################################
step <- function( object, ... ) UseMethod( "step" )
## ############################################################################
step.default <- stats::step
## ############################################################################
step.georob <- function( object, scope, scale = 0,
direction = c( "both", "backward", "forward" ), trace = 1, keep = NULL, steps = 1000,
k = 2, data = NULL, fixed = TRUE, use.fitted.param =TRUE, verbose = 0,
ncores = 1, ... )
{
## step method for class georob
## 2014-03-12 AP
## code of step{stats} complemented by argument fixed to control whether
## variogram parameters should be kept fixed
mydeviance <- function(x, ...)
{
dev <- deviance(x, ...)
if(!is.null(dev)) dev else extractAIC(x, k=0, ...)[2L]
}
cut.string <- function(string)
{
if(length(string) > 1L)
string[-1L] <- paste0("\n", string[-1L])
string
}
re.arrange <- function(keep)
{
namr <- names(k1 <- keep[[1L]])
namc <- names(keep)
nc <- length(keep)
nr <- length(k1)
array(unlist(keep, recursive = FALSE), c(nr, nc), list(namr, namc))
}
step.results <- function(models, fit, object, usingCp=FALSE)
{
change <- sapply(models, "[[", "change")
rd <- sapply(models, "[[", "deviance")
dd <- c(NA, abs(diff(rd)))
rdf <- sapply(models, "[[", "df.resid")
ddf <- c(NA, diff(rdf))
AIC <- sapply(models, "[[", "AIC")
heading <- c("Stepwise Model Path \nAnalysis of Deviance Table",
"\nInitial Model:", deparse(formula(object)),
"\nFinal Model:", deparse(formula(fit)),
"\n")
aod <- data.frame(Step = I(change), Df = ddf, Deviance = dd,
"Resid. Df" = rdf, "Resid. Dev" = rd, AIC = AIC,
check.names = FALSE)
if(usingCp) {
cn <- colnames(aod)
cn[cn == "AIC"] <- "Cp"
colnames(aod) <- cn
}
attr(aod, "heading") <- heading
##stop gap attr(aod, "class") <- c("anova", "data.frame")
fit$anova <- aod
fit
}
Terms <- terms(object)
object$call$formula <- object$formula <- Terms
md <- missing(direction)
direction <- match.arg(direction)
backward <- direction == "both" | direction == "backward"
forward <- direction == "both" | direction == "forward"
if(missing(scope)) {
fdrop <- numeric()
fadd <- attr(Terms, "factors")
if(md) forward <- FALSE
}
else {
if(is.list(scope)) {
fdrop <- if(!is.null(fdrop <- scope$lower))
attr(terms(update.formula(object, fdrop)), "factors")
else numeric()
fadd <- if(!is.null(fadd <- scope$upper))
attr(terms(update.formula(object, fadd)), "factors")
}
else {
fadd <- if(!is.null(fadd <- scope))
attr(terms(update.formula(object, scope)), "factors")
fdrop <- numeric()
}
}
models <- vector("list", steps)
if(!is.null(keep)) keep.list <- vector("list", steps)
n <- nobs(object, use.fallback = TRUE) # might be NA
fit <- object
bAIC <- extractAIC(fit, scale, k = k, ...)
edf <- bAIC[1L]
bAIC <- bAIC[2L]
if(is.na(bAIC))
stop("AIC is not defined for this model, so 'step' cannot proceed")
if(bAIC == -Inf)
stop("AIC is -infinity for this model, so 'step' cannot proceed")
nm <- 1
## Terms <- fit$terms
if(trace) {
cat("Start: AIC=", format(round(bAIC, 2)), "\n",
cut.string(deparse(formula(fit))), "\n\n", sep = "")
utils::flush.console()
}
## FIXME think about df.residual() here
models[[nm]] <- list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = "", AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
usingCp <- FALSE
while(steps > 0) {
steps <- steps - 1
AIC <- bAIC
ffac <- attr(Terms, "factors")
scope <- factor.scope(ffac, list(add = fadd, drop = fdrop))
aod <- NULL
change <- NULL
if(backward && length(scope$drop)) {
aod <- drop1.georob(
fit, scope$drop, scale = scale,
k = k, trace = trace, data = data, fixed = fixed,
use.fitted.param = use.fitted.param, verbose = verbose,
ncores = ncores, ...
)
rn <- row.names(aod)
row.names(aod) <- c(rn[1L], paste("-", rn[-1L], sep=" "))
## drop zero df terms first: one at time since they
## may mask each other
if(any(aod$Df == 0, na.rm=TRUE)) {
zdf <- aod$Df == 0 & !is.na(aod$Df)
change <- rev(rownames(aod)[zdf])[1L]
}
}
if(is.null(change)) {
if(forward && length(scope$add)) {
aodf <- add1.georob(
fit, scope$add, scale = scale,
k = k, trace = trace, data = data, fixed = fixed,
use.fitted.param = use.fitted.param, verbose = verbose,
ncores = ncores, ...
)
rn <- row.names(aodf)
row.names(aodf) <- c(rn[1L], paste("+", rn[-1L], sep=" "))
aod <-
if(is.null(aod)) aodf
else rbind(aod, aodf[-1, , drop = FALSE])
}
attr(aod, "heading") <- NULL
## need to remove any terms with zero df from consideration
nzdf <- if(!is.null(aod$Df))
aod$Df != 0 | is.na(aod$Df)
aod <- aod[nzdf, ]
if(is.null(aod) || ncol(aod) == 0) break
nc <- match(c("Cp", "AIC"), names(aod))
nc <- nc[!is.na(nc)][1L]
o <- order(aod[, nc])
if(trace) print(aod[o, ])
if(o[1L] == 1) break
change <- rownames(aod)[o[1L]]
}
usingCp <- match("Cp", names(aod), 0L) > 0L
## may need to look for a `data' argument in parent
## was:
## fit <- update(fit, paste("~ .", change), evaluate = FALSE)
## fit <- eval.parent(fit)
nfit <- switch(
substr( change, 1, 1 ),
"-" = attr( aod, "nfit" ),
"+" = attr( aodf, "nfit" )
)
param <- nfit[["param"]]
aniso <- nfit[["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( fit[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
if( is.null( data ) ){
fit <- update(
fit, paste("~ .", change),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
fit <- update(
fit, paste("~ .", change), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
fit[["call"]] <- nfit[["call"]]
fit[["initial.objects"]] <- nfit[["initial.objects"]]
nnew <- nobs(fit, use.fallback = TRUE)
if(all(is.finite(c(n, nnew))) && nnew != n)
stop("number of rows in use has changed: remove missing values?")
Terms <- terms(fit)
bAIC <- extractAIC(fit, scale, k = k, ...)
edf <- bAIC[1L]
bAIC <- bAIC[2L]
if(trace) {
cat("\nStep: AIC=", format(round(bAIC, 2)), "\n",
cut.string(deparse(formula(fit))), "\n\n", sep = "")
utils::flush.console()
}
## add a tolerance as dropping 0-df terms might increase AIC slightly
if(bAIC >= AIC + 1e-7) break
nm <- nm + 1
## FIXME: think about using df.residual() here.
models[[nm]] <-
list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = change, AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
}
if(!is.null(keep)) fit$keep <- re.arrange(keep.list[seq(nm)])
step.results(models = models[seq(nm)], fit, object, usingCp)
}
Try the georob package in your browser
Any scripts or data that you put into this service are public.
georob documentation built on May 2, 2019, 6:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#####################################
# #
# Methoden fuer Klasse "georob" #
# #
#####################################
# add1.georob
# deviance.georob
# drop1.georob
# extractAIC.georob
# logLik.georob
# step.default
# step.georob
# model.frame.georob
# model.matrix.georob
# nobs.georob
# print.georob
# ranef.georob
# residuals.georob
# resid.georob
# rstandard.georob
# rstudent.georob
# summary.georob,
# print.summary.georob
# vcov.georob
# waldtest.georob
# cv.georob (in separatem File)
# rstudent.cv.georob
# summary.cv.georob
# print.summary.cv.georob
# 2011-08-11 A. Papritz
## ##############################################################################
model.frame.georob <-
function(
formula, ...
)
{
## model.frame method for class georob
## 2012-12-19 A. Papritz
class( formula ) <- "lm"
model.frame( formula, ... )
}
## ##############################################################################
model.matrix.georob <-
function(
object, ...
)
{
## model.matrix method for class georob
## 2012-12-19 A. Papritz
class( object ) <- "lm"
model.matrix( object, ... )
}
## ##############################################################################
nobs.georob <-
function(
object, ...
)
{
## nobs method for class georob
## 2012-12-19 A. Papritz
class( object ) <- "lm"
nobs( object, ... )
}
## ##############################################################################
print.georob <-
function(
x, digits = max(3, getOption("digits") - 3), ...
)
{
## Print method for class "georob".
## Arguments:
## x an object generated by f.georob.initial.guess
## digits number of digits
## 2011-08-13 A. Papritz
## 2012-02-07 AP change for anisotropic variograms
## 2012-12-18 AP invisible(x)
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-02 AP new transformation of rotation angles
## code borrowed from print.lmrob for printing fixed effects coeffcients
cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
cat("\nFixed effects coefficients:\n")
print(
format( coef( x), digits = digits ), print.gap = 2,
quote = FALSE
)
## print variogram parameters
cat("\n")
cat( "Variogram: ", x[["variogram.model"]], "\n" )
param <- x[["param"]]
names( param ) <- ifelse(
x[["initial.objects"]][["fit.param"]],
names( param ),
paste( names( param ), "(fixed)", sep = "" )
)
print(
format( param, digits = digits ), print.gap = 2,
quote = FALSE
)
## print anisotropy parameters
if( !x[["aniso"]][["isotropic"]] ){
cat("\n")
cat( "Anisotropy parameters: ", "\n" )
aniso <- x[["aniso"]][["aniso"]]
names( aniso ) <- ifelse(
x[["initial.objects"]][["fit.aniso"]],
names( aniso ),
paste( names( aniso ), "(fixed)", sep = "" )
)
print(
format( aniso, digits = digits ), print.gap = 2,
quote = FALSE
)
}
invisible( x )
}
## ##############################################################################
ranef.georob <- random.effects.georob <-
function(
object,
standard = FALSE,
...
)
{
## Function extracts the random effects (bhat) from georob object
## (similar to ranef.lme{nlme})
## Arguments:
## object fitted georob object
## standard an optional logical value indicating whether the estimated random effects
## should be "standardized" (i.e. divided by the estimated standard error.
## Defaults to FALSE.
## ... further arguments passed to method (currently not used)
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2012-01-05 AP modified for compress storage of matrices
## 2012-10-18 AP changes for new definition of eta
## 2012-11-26 AP method for random.effects
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-05-06 AP changes for solving estimating equations for xi
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## temporarily redefine na.action component of object
object.na <- object[["na.action"]]
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
Valpha.objects <- expand( object[["Valpha.objects"]] )
covmat <- expand( object[["cov"]] )
bhat <- object[["bhat"]]
if( standard ){
if( is.null( covmat[["cov.bhat"]] ) ){
## compute standard errors of residuals
X <- model.matrix(
terms( object ),
model.frame( object )
)[!duplicated( object[["Tmat"]] ), , drop = FALSE]
r.cov <- compute.covariances(
Valpha.objects = Valpha.objects,
Aalpha = object[["Aalpha"]],
Palpha = object[["Palpha"]],
rweights = object[["rweights"]],
XX = X, TT = object[["Tmat"]], names.yy = rownames( model.frame( object ) ),
nugget = object[["param"]]["nugget"],
eta = sum( object[["param"]][c( "variance", "snugget")] ) / object[["param"]]["nugget"],
expectations = object[["expectations"]],
cov.bhat = TRUE, full.cov.bhat = FALSE,
cov.betahat = FALSE,
cov.bhat.betahat = FALSE,
cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
cov.delta.bhat.betahat = FALSE,
cov.ehat = FALSE, full.cov.ehat = FALSE,
cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = FALSE,
verbose = 0
)
if( r.cov[["error"]] ) stop(
"an error occurred when computing the variances of the random effects"
)
se <- sqrt( r.cov[["cov.bhat"]] )
} else {
## extract standard errors of residuals from georob object
if( is.matrix( covmat[["cov.bhat"]] ) ){
se <- sqrt( diag( covmat[["cov.bhat"]] ) )
} else {
se <- sqrt( covmat[["cov.bhat"]] )
}
}
bhat <- bhat / se
}
naresid( object.na, bhat )
}
## ##############################################################################
fixef.georob <- fixed.effects.georob <-
function(
object,
...
)
{
## Function extracts residuals from georob object (based on residuals.lm {stats})
## Arguments:
## object fitted georob object
## type character, type of resdiuals to computed
## ... further arguments passed to methods
## 2012-11-26 A. Papritz
## 2013-06-12 AP substituting [["x"]] for $x in all lists
object[["coef"]]
}
## ##############################################################################
residuals.georob <- resid.georob <-
function(
object,
type = c("working", "response", "deviance", "pearson", "partial" ),
terms = NULL,
level = 1,
...
)
{
## Function extracts residuals from georob object (based on residuals.lm {stats})
## Arguments:
## object fitted georob object
## type character, type of resdiuals to computed
## level integer scalar to select whether ehat (level == 1) or
## ehat + bhat (level == 0) is returned,
## only effective for type %in% c( "working", "response", "partial" )
## ... further arguments passed to methods
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2013-05-31 AP modified for computing partial residuals for single terms
## 2013-06-12 AP substituting [["x"]] for $x in all lists
type <- match.arg( type )
if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
## temporarily redefine na.action component of object
object.na <- object[["na.action"]]
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
r <- object[["residuals"]]
res <- switch(
type,
working = ,
response = r,
deviance = ,
pearson = if( is.null(object[["weights"]]) ) r else r * sqrt(object[["weights"]]),
partial = r
)
if( level == 0 && any( type %in% c( "working", "response", "partial" ) ) ){
res <- res + ranef( object, standard = FALSE )[object[["Tmat"]]]
}
if( type == "partial" )
res <- res + predict( object, type = "terms", terms = terms )[["fit"]]
drop( res )
naresid( object.na, res )
}
## ##############################################################################
rstandard.georob <-
function( model, level = 1, ... )
{
## Function extracts standardized residuals from georob object
## Arguments:
## model fitted georob object
## level integer scalar to select whether ehat (level == 1) or
## ehat + bhat (level == 0) is returned,
## ... further arguments (currently not used)
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
## 2012-01-05 AP modified for compress storage of matrices
## 2012-10-18 AP changes for new definition of eta
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-05-06 AP changes for solving estimating equations for xi
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## temporarily redefine na.action component of model
model.na <- model[["na.action"]]
if( identical( class( model[["na.action"]] ), "exclude" ) ){
class( model[["na.action"]] ) <- "omit"
}
Valpha.objects <- expand( model[["Valpha.objects"]] )
covmat <- expand( model[["cov"]] )
if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
if(
( is.null( covmat[["cov.ehat"]] ) & level == 1 ) ||
( is.null( covmat[["cov.ehat.p.bhat"]] ) & level == 0 )
){
## compute standard errors of residuals
X <- model.matrix(
terms( model),
model.frame( model )
)[!duplicated( model[["Tmat"]] ), , drop = FALSE]
r.cov <- compute.covariances(
Valpha.objects = Valpha.objects,
Aalpha = model[["Aalpha"]],
Palpha = model[["Palpha"]],
rweights = model[["rweights"]],
XX = X, TT = model[["Tmat"]], names.yy = rownames( model.frame( model ) ),
nugget = model[["param"]]["nugget"],
eta = sum( model[["param"]][c( "variance", "snugget")] ) / model[["param"]]["nugget"],
expectations = model[["expectations"]],
cov.bhat = FALSE, full.cov.bhat = FALSE,
cov.betahat = FALSE,
cov.bhat.betahat = FALSE,
cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
cov.delta.bhat.betahat = FALSE,
cov.ehat = level == 1, full.cov.ehat = FALSE,
cov.ehat.p.bhat = level == 0, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = FALSE,
verbose = 0
)
if( r.cov[["error"]] ) stop(
"an error occurred when computing the variance of the residuals"
)
if( level == 1 ){
covmat[["cov.ehat"]] <-r.cov[["cov.ehat"]]
} else {
covmat[["cov.ehat.p.bhat"]] <-r.cov[["cov.ehat.p.bhat"]]
}
}
## extract standard errors of residuals from georob object
if( level == 1 ){
se <- covmat[["cov.ehat"]]
} else {
se <- covmat[["cov.ehat.p.bhat"]]
}
if( is.matrix( se ) ){
se <- sqrt( diag( se ) )
} else {
se <- sqrt( se )
}
## compute standardized residuals
naresid( model.na, residuals( model, level = level ) / se )
}
## ##############################################################################
rstudent.georob <-
function( model, ... )
{
## Function computes studentized residuals for fitted georob object
## Arguments:
## model fitted georob object
## data data frame that was used to fit model
## ... further arguments passed to cv.georob
## 2011-12-22 A. Papritz
## 2013-06-12 AP substituting [["x"]] for $x in all lists
if( !identical( class( model )[1], "georob" ) ) stop(
"model is not of class 'georob'"
)
r.cv <- cv( model, ... )
rstudent( model = r.cv )
}
## ##############################################################################
summary.georob <-
function (
object, correlation = FALSE,
signif = 0.95,
...
)
{
## ToDos:
## - Terms Objekt einfuegen
## - ausgewaehlte Angaben zur Fitmethode ausgeben
## - Wald-Test des Modells y ~ 1
## 2012-01-05 A. Papritz
## 2012-01-05 AP modified for compress storage of matrices
## 2012-01-31 AP corretion of error for computing CI for variogram parameters
## 2012-02-07 AP change for anisotropic variograms
## 2012-03-29 AP warning messages inserted
## 2012-05-23 ap correction of error for computing correlation matrix of variogram parameters
## 2012-11-04 AP handling compressed cov.betahat
## 2012-11-27 AP changes in parameter back-transformation
## 2013-04-23 AP new names for robustness weights
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-03 AP new transformation of rotation angles
covmat <- expand( object[["cov"]] )
ans <- object[c(
"call", "residuals", "bhat", "rweights", "converged", "convergence.code",
"iter", "loglik", "variogram.model", "gradient",
"tuning.psi", "df.residual"
)]
ans <- c( ans, object[["initial.objects"]]["fit.param"] )
if( !object[["aniso"]][["isotropic"]] ) ans[["fit.param"]] <- c(
ans[["fit.param"]], object[["initial.objects"]][["fit.aniso"]]
)
ans[["terms"]] <- NA
ans[["scale"]] <- sqrt(object[["param"]]["nugget"])
ans[["control"]][["method"]] <- "TODO: PRINT GLSROB CONTROL PARAMETERS HERE"
se <- sqrt(diag(covmat[["cov.betahat"]]))
est <- object[["coefficients"]]
tval <- est/se
ans[["coefficients"]] <- cbind(
est, se, tval, 2 * pt( abs(tval), object[["df.residual"]], lower.tail = FALSE )
)
dimnames( ans[["coefficients"]] ) <- list(
names(est), c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
)
if( correlation ){
ans[["correlation"]] <- covmat[["cov.betahat"]] / outer(se, se)
}
ans[["param"]] <- as.matrix( object[["param"]], ncol = 1 )
if( !object[["aniso"]][["isotropic"]] ) ans[["param"]] <- rbind(
ans[["param"]],
as.matrix( object[["aniso"]][["aniso"]], ncol = 1 )
)
colnames( ans[["param"]] ) <- "Estimate"
## compute confidence intervals of variogram parameters from observed
## Fisher information matrix (Gaussian REML only)
if( !is.null( object[["hessian"]] ) ){
## initialization
cor.tf.param <- cov.tf.param <- matrix(
NA, nrow = nrow( object[["hessian"]] ), ncol = nrow( object[["hessian"]] ),
dimnames = dimnames( object[["hessian"]] )
)
se <- rep( NA, nrow( object[["hessian"]] ) )
names( se ) <- rownames( object[["hessian"]])
ci <- matrix( NA, nrow = nrow( ans[["param"]] ), ncol = 2 )
colnames( ci ) <- c( "Lower", "Upper" )
rownames( ci ) <- rownames( ans[["param"]] )
## select parameters that are not on boundary of parameter space
sr <- !apply( object[["hessian"]], 1, function( x ) all( is.na( x ) ) )
if( sum( sr ) > 0 ){
t.chol <- try( chol( object[["hessian"]][sr, sr] ), silent = TRUE )
if( !identical( class( t.chol ), "try-error" ) ){
## compute covariance matrix of fitted transformed parameters
cov.tf.param[sr, sr] <- chol2inv( t.chol )
## correlation matrix and standard errors of fitted transformed
## parameters
if( sum( sr ) > 1 ){
cor.tf.param[sr, sr] <- cov2cor( cov.tf.param[sr, sr] )
} else {
cor.tf.param[sr, sr] <- 1.
}
se[sr] <- sqrt( diag( cov.tf.param )[sr] )
## compute confidence interval on original scale of parameters
sel.names <- names( object[["param"]][object[["initial.objects"]][["fit.param"]]] )
if( !object[["aniso"]][["isotropic"]] ) sel.names <- c(
sel.names,
names( object[["aniso"]][["aniso"]][object[["initial.objects"]][["fit.aniso"]]] )
)
sel.names <- sel.names[sr]
ci[sel.names, ] <- t(
sapply(
sel.names,
function( x, param, se, param.tf, trafo.fct, inv.trafo.fct ){
inv.trafo.fct[[param.tf[x]]](
trafo.fct[[param.tf[x]]]( param[x] ) +
c(-1, 1) * se[x] * qnorm( (1-signif)/2, lower.tail = FALSE )
)
},
param = c( object[["param"]], object[["aniso"]][["aniso"]] ),
se = se,
param.tf = object[["param.tf"]],
trafo.fct = object[["fwd.tf"]],
inv.trafo.fct = object[["bwd.tf"]]
)
)
if( !object[["aniso"]][["isotropic"]] ){
## map angles to halfcircle
if( !object[["aniso"]][["isotropic"]] ){
sel <- match( "omega", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < 0., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 180., ci[sel, ] - 180., ci[sel, ] )
}
sel <- match( "phi", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < 0., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 180., ci[sel, ] - 180., ci[sel, ] )
}
sel <- match( "zeta", rownames(ci) )
if( !is.na( sel ) ){
ci[sel, ] <- ifelse( ci[sel, ] < -90., ci[sel, ] + 180., ci[sel, ] )
ci[sel, ] <- ifelse( ci[sel, ] > 90., ci[sel, ] - 180., ci[sel, ] )
}
}
}
ans[["param"]] <- cbind( ans[["param"]], ci )
if( correlation ) ans[["cor.tf.param"]] <- cor.tf.param
} else {
warning(
"Hessian not positive definite:",
"\nconfidence intervals of variogram parameters cannot be computed"
)
}
}
}
ans[["se.residuals"]] <- if( !is.null( covmat[["cov.ehat"]] ) ){
if( is.matrix( covmat[["cov.ehat"]] ) ){
sqrt( diag( covmat[["cov.ehat"]] ) )
} else {
sqrt( covmat[["cov.ehat"]] )
}
} else NULL
class( ans ) <- c( "summary.georob" )
ans
}
## ##############################################################################
print.summary.georob <-
function (
x, digits = max(3, getOption("digits") - 3),
signif.stars = getOption("show.signif.stars"),
...
)
{
## ToDos:
## - Ausgabe df
## - Wald-Test des Modells y ~ 1
## - ausgewaehlte Angaben zur Fitmethode ausgeben
## 2012-01-05 A. Papritz
## 2012-01-31 AP small changes
## 2012-02-07 AP change for anisotropic variograms
## 2012-12-18 AP invisible(x)
## 2013-04-23 AP new names for robustness weights
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-01-23 AP prints maximized loglik even if all parameters fixed
cat("\nCall:")
cat( paste( deparse(x[["call"]]), sep = "\n", collapse = "\n"), "\n", sep = "" )
cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
if( is.na( x[["converged"]] ) ){
cat( "\nEstimation with fixed variogram parameters\n" )
} else {
if(!(x[["converged"]])) {
cat(
"\nAlgorithm did not converge, diagnostic code: ",
x[["convergence.code"]], "\n"
)
} else {
cat(
"\nConvergence in", x[["iter"]][1], "function and",
x[["iter"]][2], "Jacobian/gradient evaluations\n"
)
}
attr( x[["gradient"]], "eeq.emp" ) <- NULL
attr( x[["gradient"]], "eeq.exp" ) <- NULL
cat( "\nEstimating equations (gradient)\n")
print( x[["gradient"]], digits = digits, ... )
}
if( x[["tuning.psi"]] >=
georob.control()[["tuning.psi.nr"]] ) cat(
"\nMaximized restricted log-likelihood:",
x[["loglik"]], "\n"
)
df <- x[["df.residual"]]
bhat <- x[["bhat"]]
cat( "\nPredicted latent variable (z):\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(bhat), names = nam )
print( rq, digits = digits, ...)
}
else print( bhat, digits = digits, ...)
resid <- x[["residuals"]]
cat( "\nResiduals (epsilon):\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(resid), names = nam )
print( rq, digits = digits, ...)
}
else print( resid, digits = digits, ...)
if( !is.null( x[["se.residuals"]] ) ){
resid <- x[["residuals"]] / x[["se.residuals"]]
cat( "\nStandardized residuals:\n")
if(df > 5){
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- structure( quantile(resid), names = nam )
print( rq, digits = digits, ...)
}
else print( resid, digits = digits, ...)
}
cat( "\nVariogram: ", x[["variogram.model"]], "\n" )
rownames( x[["param"]] ) <- ifelse(
x[["fit.param"]],
rownames( x[["param"]] ),
paste( rownames( x[["param"]] ), "(fixed)", sep = "" )
)
## print( format( x[["param"]], digits = digits ), print.gap = 2, quote = FALSE )
printCoefmat(
x[["param"]], digits = digits, signif.stars = FALSE, ...
)
if( !is.null( x[["cor.tf.param"]] ) ){
correl <- x[["cor.tf.param"]]
p <- NCOL(correl)
if( p > 1 ){
cat("\nCorrelation of (transformed) variogram parameters:\n")
correl <- format(round(correl, 2), nsmall = 2,
digits = digits)
correl[!lower.tri(correl)] <- ""
print(correl[-1, -p, drop = FALSE], quote = FALSE)
}
}
cat( "\nFixed effects coefficients:\n" )
printCoefmat(
x[["coefficients"]], digits = digits, signif.stars = signif.stars, ...
)
cat(
"\nResidual standard error (sqrt(nugget)):",
format(signif(x[["scale"]], digits)), "\n"
)
correl <- x[["correlation"]]
if( !is.null(correl) ){
p <- NCOL(correl)
if( p > 1 ){
cat("\nCorrelation of fixed effects coefficients:\n")
correl <- format(round(correl, 2), nsmall = 2,
digits = digits)
correl[!lower.tri(correl)] <- ""
print(correl[-1, -p, drop = FALSE], quote = FALSE)
}
}
cat("\n")
summarizeRobWeights(x[["rweights"]], digits = digits, ... )
invisible( x )
}
## ##############################################################################
vcov.georob <-
function( object, ... )
{
## 2012-11-04 AP handling compressed cov.betahat
## 2013-06-12 AP substituting [["x"]] for $x in all lists
result <- expand( object[["cov"]][["cov.betahat"]] )
attr( result, "struc" ) <- NULL
result
}
## ##############################################################################
waldtest.georob <-
function(
object, ..., vcov = NULL, test = c("Chisq", "F"), name = NULL,
fixed = TRUE
)
{
## waldtest method for class georob
## 2012-02-08 AP change for anisotropic variograms
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-02-19 AP correcting option for verbose output
## 2014-05-15 AP changes for version 3 of RandomFields
test <- match.arg( test )
## refit object with fixed variogram parameters
if( fixed ) {
cat( "\nWald-Test with fixed variogram parameters of model 1\n\n" )
object <- update(
object,
param = object[["param"]],
aniso = object[["aniso"]][["aniso"]],
fit.param = c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )],
fit.aniso = c(
f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
),
verbose = 0
)
}
## Wald-Test
waldtest.default(
object = object, ..., vcov = vcov, test = test, name = name
)
}
## ##############################################################################
logLik.georob <-
function( object, REML = FALSE, ... )
{
## 2012-12-22 method for extracting (restricted) loglikelihood
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-02-27 AP computing 'pseudo' likelihood for robustly fitted models
val <- if( REML ){
val <- object[["loglik"]]
} else {
D <- deviance( object, ... )
-0.5 * (
D + attr( D, "log.det.covmat" ) + length( object[["residuals"]] ) * log( 2 * pi )
)
}
attr(val, "nall") <- length(object[["residuals"]])
attr(val, "nobs") <- object[["df.residual"]]
attr(val, "df") <- length(object[["coefficients"]]) +
sum( object[["initial.objects"]][["fit.param"]] ) +
sum( object[["initial.objects"]][["fit.aniso"]])
class(val) <- "logLik"
val
}
## ##############################################################################
deviance.georob <-
function(
object, warn = TRUE, ...
)
{
## deviance method for class georob
## 2012-12-22 A. Papritz
## 2013-05-23 AP correct handling of missing observations
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-02-27 AP computing 'pseudo' deviance for robustly fitted models
## 2014-03-12 AP option for suppressing warnings
## redefine na.action component of object
if( identical( class( object[["na.action"]] ), "exclude" ) ){
class( object[["na.action"]] ) <- "omit"
}
## determine factor to compute heteroscedastic nugget
if( object[["tuning.psi"]] < georob.control()[["tuning.psi.nr"]] ){
if( warn ) warning(
"deviance for robustly fitted model approximated by deviance of\n",
" equivalent model with heteroscedastic nugget"
)
w <- object[["rweights"]]
} else {
w <- 1.
}
## invert covariance matrix of observations
Valpha.objects <- expand( object[["Valpha.objects"]] )
G <- sum( object[["param"]][c("variance", "snugget")] ) * Valpha.objects[["Valpha"]]
G <- G[object[["Tmat"]], object[["Tmat"]]]
diag( G ) <- diag( G ) + object[["param"]]["nugget"] / w
iucf <- try(
backsolve(
chol( G ),
diag( length( object[["Tmat"]] ) ),
k = length( object[["Tmat"]] )
),
silent = TRUE
)
## compute deviance
if( identical( class( iucf ), "try-error" ) ) {
stop( "(generalized) covariance matrix of observations not positive definite" )
} else {
result <- sum( colSums( residuals( object, level = 0 ) * iucf )^2 )
attr( result, "log.det.covmat" ) <- -2 * sum( log( diag( iucf ) ) )
}
result
}
## ##############################################################################
extractAIC.georob <- function( fit, scale = 0, k = 2, ... )
{
## extractAIC method for class georob
## 2014-02-27 AP
edf <- sum( !is.na( fit[["coefficients"]] ) ) +
sum( fit[["initial.objects"]][["fit.param"]] ) +
sum( fit[["initial.objects"]][["fit.aniso"]] )
loglik <- logLik( fit, ... )
c(edf, -2 * loglik + k * edf)
}
## ##############################################################################
safe_pchisq <- function(q, df, ...)
## same code as safe_pchisq{stats}
{
df[df <= 0] <- NA
pchisq(q=q, df=df, ...)
}
## ##############################################################################
add1.georob <- function( object, scope, scale = 0, test=c("none", "Chisq" ),
k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
ncores = 1, ... )
{
## add1 method for class georob based on add1.default{stats}
## 2014-03-12 AP
## 2014-04-24 AP function returns only variogram parameters of best fitted object
## 2014-05-15 AP changes for version 3 of RandomFields
## auxiliary function for fitting model and extracting aic
f.aux <- function(
tt, scale, k, trace,
fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
verbose, ...
){
converged <- TRUE
if(trace > 1) {
cat("\ntrying +", tt, "\n", sep = "")
utils::flush.console()
}
if( fixed ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . +", tt)),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . +", tt)), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
} else {
if( use.fitted.param ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . +", tt)),
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . +", tt)), data = data,
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
}
}
else {
if( is.null( data ) ){
nfit <- update( object, as.formula(paste("~ . +", tt)),
verbose = verbose
)
} else {
nfit <- update( object, as.formula(paste("~ . +", tt)), data = data,
verbose = verbose
)
}
}
if( !nfit$converged ){
converged <- FALSE
if( verbose > 0 ) cat(
"lack of convergence when fitting model", paste("~ . +", tt),
"\nconvergence code:", nfit$convergence.code, "\n"
)
}
}
nnew <- nobs( nfit, use.fallback = TRUE )
if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
"number of rows in use has changed: remove missing values?"
)
df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
names( df.aic ) <- c("df", "AIC", "converged")
attr( df.aic, "nfit" ) <- list(
param = nfit[["param"]],
aniso = nfit[["aniso"]][["aniso"]],
initial.objects = nfit[["initial.objects"]],
call = nfit[["call"]]
)
df.aic
}
## check scope
if( missing(scope ) || is.null( scope ) ) stop("no terms in scope")
if( !is.character( scope ) ) scope <- add.scope( object, update.formula(object, scope) )
if( !length(scope)) stop( "no terms in scope for adding to object" )
## initialize result
ns <- length( scope )
ans <- matrix(
nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
)
ans[1L, ] <- extractAIC( object, scale, k = k, ... )
## loop over all components of scope
n0 <- nobs( object, use.fallback = TRUE )
env <- environment( formula(object) )
param <- object[["param"]]
aniso <- object[["aniso"]][["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
## prepare for parallel execution
if( ncores > 1 ){
ncores <- min( ncores, ns, detectCores() )
trace <- 0
verbose <- 0
}
if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
if( is.null( data ) ) stop(
"argument 'data' required for parallel execution on windows OS"
)
## create a SNOW cluster on windows OS
cl <- makePSOCKcluster( ncores, outfile = "")
## export required items to workers
junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
result <- parLapply(
cl,
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
)
stopCluster(cl)
} else {
## fork child processes on non-windows OS
result <- mclapply(
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
mc.cores = ncores, mc.allow.recursive = FALSE, ...
)
}
fits <- list(
list(
param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
initial.objects = object[["initial.objects"]],
call = object[["call"]]
)
)
fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
result <- t( simplify2array( result ) )
ans[2:NROW(ans), ] <- result[, 1:2]
if(!all( as.logical(result[, "converged"]) ) ) warning(
"lack of convergence when fitting models"
)
## store results
dfs <- ans[, 1L] - ans[1L, 1L]
dfs[1L] <- NA
aod <- data.frame( Df = dfs, AIC = ans[, 2L] )
## likelihood ratio test
test <- match.arg(test)
if(test == "Chisq") {
dev <- ans[, 2L] - k*ans[, 1L]
dev <- dev[1L] - dev; dev[1L] <- NA
nas <- !is.na(dev)
P <- dev
P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail=FALSE)
aod[, c("LRT", "Pr(>Chi)")] <- list(dev, P)
}
## format results
head <- c(
"Single term additions", "\nModel:", deparse(formula(object)),
if(scale > 0) paste("\nscale: ", format(scale), "\n")
)
class(aod) <- c("anova", "data.frame")
attr( aod, "heading") <- head
attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
aod
}
## ##############################################################################
drop1.georob <- function( object, scope, scale = 0, test=c( "none", "Chisq" ),
k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
ncores = 1, ... )
{
## drop1 method for class georob based on drop1.default{stats}
## 2014-03-12 AP
## 2014-04-24 AP function returns only variogram parameters of best fitted object
## 2014-05-15 AP changes for version 3 of RandomFields
## auxiliary function for fitting model and extracting aic
f.aux <- function(
tt, scale, k, trace,
fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
verbose, ...
){
converged <- TRUE
if(trace > 1) {
cat("\ntrying -", tt, "\n", sep = "")
utils::flush.console()
}
if( fixed ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . -", tt)),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . -", tt)), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
} else {
if( use.fitted.param ){
if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . -", tt)),
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
} else {
nfit <- update(
object, as.formula(paste("~ . -", tt)), data = data,
param = param, aniso = aniso, initial.param = "no", verbose = verbose
)
}
}
else {
if( is.null( data ) ){
nfit <- update( object, as.formula(paste("~ . -", tt)),
verbose = verbose
)
} else {
nfit <- update( object, as.formula(paste("~ . -", tt)), data = data,
verbose = verbose
)
}
}
if( !nfit$converged ){
converged <- FALSE
if( verbose > 0 ) cat(
"lack of convergence when fitting model", paste("~ . -", tt),
"\nconvergence code:", nfit$convergence.code, "\n"
)
}
}
nnew <- nobs( nfit, use.fallback = TRUE )
if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
"number of rows in use has changed: remove missing values?"
)
df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
names( df.aic ) <- c("df", "AIC", "converged")
attr( df.aic, "nfit" ) <- list(
param = nfit[["param"]],
aniso = nfit[["aniso"]][["aniso"]],
initial.objects = nfit[["initial.objects"]],
call = nfit[["call"]]
)
df.aic
}
## check scope
tl <- attr( terms(object), "term.labels" )
if( missing(scope) ) {
scope <- drop.scope( object )
} else {
if( !is.character(scope) ) scope <- attr(
terms(update.formula(object, scope)), "term.labels"
)
if( !all( match(scope, tl, 0L) > 0L) ) stop("scope is not a subset of term labels")
}
## initialize result
ns <- length( scope )
ans <- matrix(
nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
)
ans[1, ] <- extractAIC( object, scale, k = k, ... )
## loop over all components of scope
n0 <- nobs( object, use.fallback = TRUE )
env <- environment( formula(object) )
param <- object[["param"]]
aniso <- object[["aniso"]][["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
## prepare for parallel execution
if( ncores > 1 ){
ncores <- min( ncores, ns, detectCores() )
trace <- 0
verbose <- 0
}
if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
if( is.null( data ) ) stop(
"argument 'data' required for parallel execution on windows OS"
)
## create a SNOW cluster on windows OS
cl <- makePSOCKcluster( ncores, outfile = "")
## export required items to workers
junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
result <- parLapply(
cl,
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
)
stopCluster(cl)
} else {
## fork child processes on non-windows OS
result <- mclapply(
scope, f.aux,
scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
object = object, data = data, param = param, aniso = aniso,
fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
mc.cores = ncores, mc.allow.recursive = FALSE, ...
)
}
fits <- list(
list(
param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
initial.objects = object[["initial.objects"]],
call = object[["call"]]
)
)
fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
result <- t( simplify2array( result ) )
ans[2:NROW(ans), ] <- result[, 1:2]
if(!all( as.logical(result[, "converged"]) ) ) warning(
"lack of convergence when fitting models"
)
## store results
dfs <- ans[1L , 1L] - ans[, 1L]
dfs[1L] <- NA
aod <- data.frame( Df = dfs, AIC = ans[,2] )
## likelihood ratio test
test <- match.arg(test)
if(test == "Chisq") {
dev <- ans[, 2L] - k*ans[, 1L]
dev <- dev - dev[1L] ; dev[1L] <- NA
nas <- !is.na(dev)
P <- dev
P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail = FALSE)
aod[, c("LRT", "Pr(>Chi)")] <- list(dev, P)
}
## format results
head <- c("Single term deletions", "\nModel:", deparse(formula(object)),
if(scale > 0) paste("\nscale: ", format(scale), "\n"))
class(aod) <- c("anova", "data.frame")
attr( aod, "heading") <- head
attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
aod
}
## ############################################################################
step <- function( object, ... ) UseMethod( "step" )
## ############################################################################
step.default <- stats::step
## ############################################################################
step.georob <- function( object, scope, scale = 0,
direction = c( "both", "backward", "forward" ), trace = 1, keep = NULL, steps = 1000,
k = 2, data = NULL, fixed = TRUE, use.fitted.param =TRUE, verbose = 0,
ncores = 1, ... )
{
## step method for class georob
## 2014-03-12 AP
## code of step{stats} complemented by argument fixed to control whether
## variogram parameters should be kept fixed
mydeviance <- function(x, ...)
{
dev <- deviance(x, ...)
if(!is.null(dev)) dev else extractAIC(x, k=0, ...)[2L]
}
cut.string <- function(string)
{
if(length(string) > 1L)
string[-1L] <- paste0("\n", string[-1L])
string
}
re.arrange <- function(keep)
{
namr <- names(k1 <- keep[[1L]])
namc <- names(keep)
nc <- length(keep)
nr <- length(k1)
array(unlist(keep, recursive = FALSE), c(nr, nc), list(namr, namc))
}
step.results <- function(models, fit, object, usingCp=FALSE)
{
change <- sapply(models, "[[", "change")
rd <- sapply(models, "[[", "deviance")
dd <- c(NA, abs(diff(rd)))
rdf <- sapply(models, "[[", "df.resid")
ddf <- c(NA, diff(rdf))
AIC <- sapply(models, "[[", "AIC")
heading <- c("Stepwise Model Path \nAnalysis of Deviance Table",
"\nInitial Model:", deparse(formula(object)),
"\nFinal Model:", deparse(formula(fit)),
"\n")
aod <- data.frame(Step = I(change), Df = ddf, Deviance = dd,
"Resid. Df" = rdf, "Resid. Dev" = rd, AIC = AIC,
check.names = FALSE)
if(usingCp) {
cn <- colnames(aod)
cn[cn == "AIC"] <- "Cp"
colnames(aod) <- cn
}
attr(aod, "heading") <- heading
##stop gap attr(aod, "class") <- c("anova", "data.frame")
fit$anova <- aod
fit
}
Terms <- terms(object)
object$call$formula <- object$formula <- Terms
md <- missing(direction)
direction <- match.arg(direction)
backward <- direction == "both" | direction == "backward"
forward <- direction == "both" | direction == "forward"
if(missing(scope)) {
fdrop <- numeric()
fadd <- attr(Terms, "factors")
if(md) forward <- FALSE
}
else {
if(is.list(scope)) {
fdrop <- if(!is.null(fdrop <- scope$lower))
attr(terms(update.formula(object, fdrop)), "factors")
else numeric()
fadd <- if(!is.null(fadd <- scope$upper))
attr(terms(update.formula(object, fadd)), "factors")
}
else {
fadd <- if(!is.null(fadd <- scope))
attr(terms(update.formula(object, scope)), "factors")
fdrop <- numeric()
}
}
models <- vector("list", steps)
if(!is.null(keep)) keep.list <- vector("list", steps)
n <- nobs(object, use.fallback = TRUE) # might be NA
fit <- object
bAIC <- extractAIC(fit, scale, k = k, ...)
edf <- bAIC[1L]
bAIC <- bAIC[2L]
if(is.na(bAIC))
stop("AIC is not defined for this model, so 'step' cannot proceed")
if(bAIC == -Inf)
stop("AIC is -infinity for this model, so 'step' cannot proceed")
nm <- 1
## Terms <- fit$terms
if(trace) {
cat("Start: AIC=", format(round(bAIC, 2)), "\n",
cut.string(deparse(formula(fit))), "\n\n", sep = "")
utils::flush.console()
}
## FIXME think about df.residual() here
models[[nm]] <- list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = "", AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
usingCp <- FALSE
while(steps > 0) {
steps <- steps - 1
AIC <- bAIC
ffac <- attr(Terms, "factors")
scope <- factor.scope(ffac, list(add = fadd, drop = fdrop))
aod <- NULL
change <- NULL
if(backward && length(scope$drop)) {
aod <- drop1.georob(
fit, scope$drop, scale = scale,
k = k, trace = trace, data = data, fixed = fixed,
use.fitted.param = use.fitted.param, verbose = verbose,
ncores = ncores, ...
)
rn <- row.names(aod)
row.names(aod) <- c(rn[1L], paste("-", rn[-1L], sep=" "))
## drop zero df terms first: one at time since they
## may mask each other
if(any(aod$Df == 0, na.rm=TRUE)) {
zdf <- aod$Df == 0 & !is.na(aod$Df)
change <- rev(rownames(aod)[zdf])[1L]
}
}
if(is.null(change)) {
if(forward && length(scope$add)) {
aodf <- add1.georob(
fit, scope$add, scale = scale,
k = k, trace = trace, data = data, fixed = fixed,
use.fitted.param = use.fitted.param, verbose = verbose,
ncores = ncores, ...
)
rn <- row.names(aodf)
row.names(aodf) <- c(rn[1L], paste("+", rn[-1L], sep=" "))
aod <-
if(is.null(aod)) aodf
else rbind(aod, aodf[-1, , drop = FALSE])
}
attr(aod, "heading") <- NULL
## need to remove any terms with zero df from consideration
nzdf <- if(!is.null(aod$Df))
aod$Df != 0 | is.na(aod$Df)
aod <- aod[nzdf, ]
if(is.null(aod) || ncol(aod) == 0) break
nc <- match(c("Cp", "AIC"), names(aod))
nc <- nc[!is.na(nc)][1L]
o <- order(aod[, nc])
if(trace) print(aod[o, ])
if(o[1L] == 1) break
change <- rownames(aod)[o[1L]]
}
usingCp <- match("Cp", names(aod), 0L) > 0L
## may need to look for a `data' argument in parent
## was:
## fit <- update(fit, paste("~ .", change), evaluate = FALSE)
## fit <- eval.parent(fit)
nfit <- switch(
substr( change, 1, 1 ),
"-" = attr( aod, "nfit" ),
"+" = attr( aodf, "nfit" )
)
param <- nfit[["param"]]
aniso <- nfit[["aniso"]]
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
alpha = FALSE, beta = FALSE, delta = FALSE,
gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( fit[["param"]] )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
if( is.null( data ) ){
fit <- update(
fit, paste("~ .", change),
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
fit <- update(
fit, paste("~ .", change), data = data,
param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
}
fit[["call"]] <- nfit[["call"]]
fit[["initial.objects"]] <- nfit[["initial.objects"]]
nnew <- nobs(fit, use.fallback = TRUE)
if(all(is.finite(c(n, nnew))) && nnew != n)
stop("number of rows in use has changed: remove missing values?")
Terms <- terms(fit)
bAIC <- extractAIC(fit, scale, k = k, ...)
edf <- bAIC[1L]
bAIC <- bAIC[2L]
if(trace) {
cat("\nStep: AIC=", format(round(bAIC, 2)), "\n",
cut.string(deparse(formula(fit))), "\n\n", sep = "")
utils::flush.console()
}
## add a tolerance as dropping 0-df terms might increase AIC slightly
if(bAIC >= AIC + 1e-7) break
nm <- nm + 1
## FIXME: think about using df.residual() here.
models[[nm]] <-
list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = change, AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
}
if(!is.null(keep)) fit$keep <- re.arrange(keep.list[seq(nm)])
step.results(models = models[seq(nm)], fit, object, usingCp)
}
Try the georob package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.