R/georob_S3methods.R

In georob: Robust Geostatistical Analysis of Spatial Data

#####################################
#                                   #
#   Methoden fuer Klasse "georob"   #
#                                   #
#####################################

# add1.georob
# coef.georob
# print.coef.georob
# deviance.georob
# drop1.georob
# extractAIC.georob
# logLik.georob
# step
# step.default
# step.georob
# model.frame.georob
# model.matrix.georob
# nobs.georob
# print.georob
# proflik
# proflik.likfit
# proflik.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

## ##############################################################################

coef.georob <- function(
  object, what = c("trend", "variogram"), ...
)
{
  ## coef method for class georob

  ## 2017-10-20 A. Papritz

  f.aux <- function(x){
    tmp <- x[["param"]]
    if( !x[["isotropic"]] ){
      tmp <- c(tmp, x[["aniso"]])
    }
    attr( tmp, "variogram.model" ) <- x[["variogram.model"]]
    tmp
  }

  what <- match.arg( what )

  res <- switch(
    what,
    trend = object[["coefficients"]],
    variogram = lapply(object[["variogram.object"]], f.aux)
  )

  if( is.list(res) ){
    if( identical( length( res ), 1L ) ){
      res <- res[[1]]
    } else {
      names( res ) <- sapply( object[["variogram.object"]], function(x) x[["variogram.model"]] )
    }
  }

  class( res ) <- "coef.georob"

  res
}

## ##############################################################################

print.coef.georob <- function(
  x, ...
)
{
  ## print method for class coef.georob

  ## 2017-10-20 A. Papritz

  xx <- unclass( x )

  if( is.list( xx ) ){
    lapply(
      xx,
      function(x){
        cat( "Variogram: ", attr(x, "variogram.model"), "\n" )
        attr( x, "variogram.model" ) <- NULL
        print( x )
        cat( "\n" )
      }
    )

  } else {
    attr( xx, "variogram.model" ) <- NULL
    print( xx )
  }

  invisible(x)

}



## ##############################################################################

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
  ## 2016-08-05 AP changes for nested variogram models

  ## 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 = 2L,
    quote = FALSE
  )

  ## print variogram parameters

  cat("\n")

  lapply(
    x[["variogram.object"]],
    function(x){

      cat( "Variogram: ", x[["variogram.model"]], "\n" )

      param <- x[["param"]]
      names( param ) <- ifelse(
        x[["fit.param"]],
        names( param ),
        paste( names( param ), "(fixed)", sep = "" )
      )
      print(
        format( param, digits = digits,
          width = 16L, justify = "right" ), print.gap = 2L,
        quote = FALSE
      )

      cat("\n")

      if( !x[["isotropic"]] ){

        aniso <- x[["aniso"]]
        names( aniso ) <- ifelse(
          x[["fit.aniso"]],
          names( aniso ),
          paste( names( aniso ), "(fixed)", sep = "" )
        )

        print(
          format( aniso, digits = digits,
            width = 16L, justify = "right" ), print.gap = 2L,
          quote = FALSE
        )
        cat("\n")
      }
    }
  )

  invisible( x )

}

## ##############################################################################

ranef.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
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-08-05 AP changes for nested variogram models
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## check arguments

  stopifnot(identical(length(standard), 1L) && is.logical(standard))

  ## temporarily redefine na.action component of object

  object.na <- object[["na.action"]]
  if( identical( class( object[["na.action"]] ), "exclude" ) ){
    class( object[["na.action"]] ) <- "omit"
  }

  Valphaxi.objects <- expand( object[["Valphaxi.objects"]] )
  zhat.objects     <- expand( object[["zhat.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 <- covariances.fixed.random.effects(
        Valphaxi.objects = Valphaxi.objects[c("Valphaxi", "Valphaxi.inverse")],
        Aalphaxi = zhat.objects[["Aalphaxi"]],
        Palphaxi = zhat.objects[["Palphaxi"]],
        Valphaxi.inverse.Palphaxi = zhat.objects[["Valphaxi.inverse.Palphaxi"]],
        rweights = object[["rweights"]],
        XX = X, TT = object[["Tmat"]], TtT = as.vector( table( 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"]], family = "gaussian",
        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,
        control.pcmp = control.pcmp(),
        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 )

}

random.effects.georob <- ranef.georob

## ##############################################################################

fixef.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
  ## 2016-07-21 AP correcting wrong component name

  object[["coefficients"]]

}

fixed.effects.georob <- fixef.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 == 0L && !identical( type, "pearson" ) ){
    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 )
}

residuals.georob <- resid.georob

## ##############################################################################

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
  ## 2104-08-19 AP correcting error when computing covariances of residuals
  ## 2016-08-05 AP changes for nested variogram models
#
  ## temporarily redefine na.action component of model

  model.na <- model[["na.action"]]
  if( identical( class( model[["na.action"]] ), "exclude" ) ){
    class( model[["na.action"]] ) <- "omit"
  }

  Valphaxi.objects <- expand( model[["Valphaxi.objects"]] )
  zhat.objects     <- expand( model[["zhat.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 == 1L ) ||
    ( is.null( covmat[["cov.ehat.p.bhat"]] ) & level == 0L )
  ){

    ## compute standard errors of residuals

    X <- model.matrix(
      terms( model),
      model.frame( model )
    )[!duplicated( model[["Tmat"]] ), , drop = FALSE]

    r.cov <- covariances.fixed.random.effects(
      Valphaxi.objects = Valphaxi.objects,
      Aalphaxi = zhat.objects[["Aalphaxi"]],
      Palphaxi = zhat.objects[["Palphaxi"]],
      Valphaxi.inverse.Palphaxi = zhat.objects[["Valphaxi.inverse.Palphaxi"]],
      rweights = model[["rweights"]],
      XX = X, TT = model[["Tmat"]], TtT = as.vector( table( model[["Tmat"]] ) ),
      names.yy = rownames( model.frame( model ) ),
      nugget = model[["variogram.object"]][[1L]][["param"]]["nugget"],
      eta = f.reparam.fwd( model[["variogram.object"]] )[[1L]][["param"]]["nugget"],
      expectations = model[["expectations"]], family = "long.tailed",
      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 == 1L, full.cov.ehat = FALSE,
      cov.ehat.p.bhat = level == 0L, full.cov.ehat.p.bhat = FALSE,
      aux.cov.pred.target = FALSE,
      control.pcmp = control.pcmp(),
      verbose = 0.
    )

    if( r.cov[["error"]] ) stop(
      "an error occurred when computing the variance of the residuals"
    )

    if( level == 1L ){
      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 == 1L ){
    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
  ## 2014-08-26 AP changes to print ml.method
  ## 2015-04-07 AP changes for fitting anisotropic variograms
  ## 2016-08-05 AP changes for nested variogram models
  ## 2019-05-24 AP correction of error when printing confidence interval of variogram parameters
  ## 2019-10-22 AP terms component taken from georob object
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()
  ## 2020-03-27 AP computing Hessian (observed Fisher information) of untransformed parameters

  ## check arguments

  stopifnot(identical(length(correlation), 1L) && is.logical(correlation))

  stopifnot(identical(length(signif), 1L) && is.numeric(signif) && signif > 0 && signif < 1)


  d2r <- pi / 180.

  covmat       <- expand( object[["cov"]] )

  ans <- object[c(
    "call", "residuals", "bhat", "rweights", "converged", "convergence.code",
    "iter", "loglik", "variogram.object", "gradient",
    "tuning.psi", "df.residual", "control", "terms"
  )]

  ans[["scale"]] <- sqrt(object[["variogram.object"]][[1L]][["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"]] <- cov2cor( covmat[["cov.betahat"]] )
  }

  ans[["param.aniso"]] <- lapply(
    ans[["variogram.object"]],
    function(object){
      res <- as.matrix( object[["param"]], ncol = 1L )
      nme.res <- names( object[["param"]] )
      fit.res <- object[["fit.param"]]
      if( !object[["isotropic"]] ){
        res <- rbind( res, as.matrix( object[["aniso"]], ncol = 1L ) )
        nme.res <- c( nme.res, names( object[["aniso"]] ) )
        fit.res <- c( fit.res, object[["fit.aniso"]] )
      }
      colnames( res ) <- "Estimate"
      nme.res[!fit.res] <- paste( nme.res[!fit.res], "(fixed)", sep="" )
      rownames( res ) <- nme.res
      res
    }
  )

  ## compute confidence intervals of variogram parameters from observed
  ## Fisher information matrix (Gaussian REML only)

  if( !is.null( object[["hessian.tfpa"]] ) ){

    ## initialization

    cor.tf.param <- cov.tf.param <- matrix(
      NA, nrow = NROW( object[["hessian.tfpa"]] ), ncol = NROW( object[["hessian.tfpa"]] ),
      dimnames = dimnames( object[["hessian.tfpa"]] )
    )

#     se <- rep( NA_real_, NROW( object[["hessian.tfpa"]] ) )
#     names( se ) <- rownames( object[["hessian.tfpa"]])
#
#     ci <- matrix( NA_real_, nrow = NROW( object[["hessian.tfpa"]] ), ncol = 2L )
#     colnames( ci ) <- c( "Lower", "Upper" )
#     rownames( ci ) <- rownames( object[["hessian.tfpa"]] )

    ## select parameters that are not on boundary of parameter space

    sr  <- !apply( object[["hessian.tfpa"]], 1L, function( x ) all( is.na( x ) ) )

    if( sum( sr ) > 0L ){

      t.chol <- try( chol( object[["hessian.tfpa"]][sr, sr, drop = FALSE] ), silent = TRUE )

      if( !identical( class( t.chol ), "try-error" ) ){

        ## compute covariance matrix of fitted transformed parameters

        cov.tf.param <- chol2inv( t.chol )
        dimnames( cov.tf.param ) <- dimnames( t.chol )

        ## correlation matrix and standard errors of fitted transformed
        ## parameters

        if( correlation ){
          ans[["cor.tf.param"]] <- cov2cor( cov.tf.param )
          colnames( ans[["cor.tf.param"]] ) <- rownames( ans[["cor.tf.param"]] ) <-
            gsub( ".__...__.", ".", colnames( ans[["cor.tf.param"]] ), fixed = TRUE )
        }

        se <- sqrt( diag( cov.tf.param ) )

        tmp <- f.aux.tf.param.fwd(
          ans[["variogram.object"]], object[["control"]][["param.tf"]],
          object[["control"]][["fwd.tf"]]
        )

        param <- tmp[["transformed.param.aniso"]][tmp[["fit.param.aniso"]]]
        param <- param[names(param) %in% names(se)]

        param.tf <- tmp[["tf.param.aniso"]][names(param)]

        se <- se[match( names(se), names(param))]

        ## confidence intervals

        ci <- t( sapply(
            1L:length(se),
            function( i, m, se, signif ){
              m[i] + c(-1., 1.) * se[i] * qnorm( (1.-signif)/2., lower.tail = FALSE )
            }, m = param, se = se, signif = signif
          ))
        colnames( ci ) <- c("Lower", "Upper")
        rownames( ci ) <- names( se )

        ci <- apply(
          ci, 2L,
          function( x, nme, bwd.tf, param.tf ){
            sapply( nme,
              function( x, bwd.tf, param.tf, param ) bwd.tf[[param.tf[x]]]( param[x] ),
              bwd.tf = bwd.tf, param.tf = param.tf, param = x
            )
          }, nme = names(se), bwd.tf = object[["control"]][["bwd.tf"]],
          param.tf = param.tf
        )

        if(is.vector(ci)) ci <- matrix(ci, nrow = 1)

        colnames( ci ) <- c("Lower", "Upper")
        rownames( ci ) <- names( se )

        ## convert to list

        tmp <- strsplit( rownames(ci), ".__...__.", fixed = TRUE )
        name.tmp <- rownames(ci) <- sapply( tmp, function(x) x[1L] )
        cmp <- sapply( tmp, function(x) x[2L] )

        ci <- tapply(
          1L:NROW(ci), factor( cmp ),
          function( i, ci ){
            ci[i, , drop = FALSE]
          }, ci = ci, simplify = FALSE
        )

        ## merge into ans[["param.aniso"]]

        ans[["param.aniso"]] <- lapply(
          1L:length(ans[["param.aniso"]]),
          function( i, pa, ci ){
            pa <- pa[[i]]
            if( i <= length(ci) ){
              ci <- ci[[i]]
            } else {
              ci <- matrix( rep( NA_real_, 2L * NROW( pa ) ), ncol = 2L )
              rownames( ci ) <- rownames(pa)
            }

            pa <- cbind( pa, Lower = rep( NA_real_, NROW(pa) ),
              Upper = rep( NA_real_, NROW(pa) ) )
            i <- match( rownames( ci ), rownames( pa ) )
            pa[i, 2L:3L] <- ci
            pa
          }, pa = ans[["param.aniso"]], ci = ci
        )

      } 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
  ## 2014-08-26 AP changes to print ml.method
  ## 2016-08-05 AP changes for nested variogram models

  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"]][1L], "function and",
        x[["iter"]][2L], "Jacobian/gradient evaluations\n"
      )
    }

    attr( x[["gradient"]], "eeq.emp" ) <- NULL
    attr( x[["gradient"]], "eeq.exp" ) <- NULL

    cat( "\nEstimating equations (gradient)\n")

    f.aux.print.gradient( x[["gradient"]],
      reparam = x[["control"]][["reparam"]] &&
        x[["tuning.psi"]] >= x[["control"]][["tuning.psi.nr"]]
    )
    #     print( x[["gradient"]], digits = digits, ... )

  }

  if( x[["tuning.psi"]] >= x[["control"]][["tuning.psi.nr"]] ){
    switch(
      x[["control"]][["ml.method"]],
      ML = cat( "\nMaximized log-likelihood:", x[["loglik"]], "\n" ),
      REML = cat( "\nMaximized restricted log-likelihood:", x[["loglik"]], "\n" )
    )
  }

  df <- x[["df.residual"]]

  bhat <- x[["bhat"]]
  cat( "\nPredicted latent variable (B):\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 <- residuals( x )
  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 <- residuals( x, type = "working" ) / 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(
    if( x[["tuning.psi"]] >= x[["control"]][["tuning.psi.nr"]] ){
      switch(
        x[["control"]][["ml.method"]],
        ML = "\n\nGaussian ML estimates\n",
        REML ="\n\nGaussian REML estimates\n"
      )
    } else {
      "\n\nRobust REML estimates\n"
    }
  )

  lapply(
    1L:length(x[["param.aniso"]]),
    function(i, x, vo){
      x <- x[[i]]
      tmp <- x[ ,1L]
      tmp <- tmp[!is.na(tmp) & tmp > 0. ]
      t.digits <- -floor( log10( min( tmp ) ) )
      cat( "\nVariogram: ", vo[[i]][["variogram.model"]], "\n" )
      printCoefmat(
        x, digits = max( digits, t.digits) , signif.stars = FALSE, ...
      )
      #       print(format(
#         signif( x, digits = 7L ), width = 16L, scientific = TRUE, justify = "left"), quote = FALSE, ...
#       )
    }, x = x[["param.aniso"]], vo = x[["variogram.object"]]
  )

  if( !is.null( x[["cor.tf.param"]] ) ){

    correl <- x[["cor.tf.param"]]
    p <- NCOL(correl)
    if( p > 1L ){
      cat("\nCorrelation of (transformed) variogram parameters:\n")
      correl <- format(round(correl, 2L), nsmall = 2L,
        digits = digits)
      correl[!lower.tri(correl)] <- ""
      print(correl[-1L, -p, drop = FALSE], quote = FALSE)
    }

  }

  cat( "\n\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 > 1L ){
      cat("\nCorrelation of fixed effects coefficients:\n")
      correl <- format(round(correl, 2L), nsmall = 2L,
        digits = digits)
      correl[!lower.tri(correl)] <- ""
      print(correl[-1L, -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("F", "Chisq"), name = NULL
  )
{

  ## 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
  ## 2015-08-28 AP computation of hessian suppressed
  ## 2015-09-01 AP keep variogram parameters always fixed
  ## 2016-07-15 AP optimization
  ## 2016-08-08 AP changes for nested variogram models
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## check arguments

  stopifnot(is.null(vcov) || is.function(vcov))
  stopifnot(is.null(name) || is.function(name))

  test <- match.arg( test )

  cl <-  object[["call"]]

  ## fix all variogram parameters

  cl <- f.call.set_allfitxxx_to_false( cl )

  ## set hessian equal to FALSE in control argument of georob call and update
  ## call

  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )

  ## set verbose = 0 in call

  cl <- f.call.set_x_to_value( cl, "verbose", 0. )

  object[["call"]] <- cl

#   ## set verbose = 0 in call
#
#   object[["call"]] <- update( object, verbose = 3., evaluate = FALSE )

  ## Wald-Test

  waldtest.default(
    object = object, ..., vcov = vcov, test = test, name = name
  )

}

## ##############################################################################

logLik.georob <-
  function( object, warn = TRUE, 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
  ## 2015-03-16 AP computing 'pseudo' likelihood for REML case
  ## 2016-08-08 AP changes for nested variogram models
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## check arguments

  stopifnot(identical(length(warn), 1L) && is.logical(warn))
  stopifnot(identical(length(REML), 1L) && is.logical(REML))

  if(
    object[["tuning.psi"]] >= object[["control"]][["tuning.psi.nr"]] &&
    (
      (identical( object[["control"]][["ml.method"]], "REML" ) && REML) ||
      (identical( object[["control"]][["ml.method"]], "ML" ) && !REML)
    )
  ){
    val <- object[["loglik"]]
  } else {
    D <- deviance( object, warn = warn, REML = REML, ... )
    if( REML ){
      val <- -0.5 * (
        D + attr( D, "log.det.covmat" ) + attr( D, "log.det.xticovmatx" ) +
        (length( object[["residuals"]] ) - length( object[["coefficients"]]) ) * log( 2. * pi )
      )
    } else {
      val <- -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"]]

  tmp <- unlist(lapply(
    object[["variogram.object"]],
    function(x) c( x[["fit.param"]], x[["fit.aniso"]] )
  ))

  attr(val, "df") <- length(object[["coefficients"]]) + sum( tmp )

  class(val) <- "logLik"
  val

}

## ##############################################################################

deviance.georob <-
  function( object, warn = TRUE, REML = FALSE, ... )
{
  ## 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
  ## 2015-03-16 AP attribute for computing maximized restricted loglikelihood
  ## 2016-08-08 AP changes for nested variogram models
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## check arguments

  stopifnot(identical(length(warn), 1L) && is.logical(warn))
  stopifnot(identical(length(REML), 1L) && is.logical(REML))

  ## 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"]] < object[["control"]][["tuning.psi.nr"]] ){
    if( warn ) warning(
      "deviance for robustly fitted model approximated by deviance of\n",
      "  equivalent Gaussian model with heteroscedastic nugget"
    )
    w <- object[["rweights"]]
  } else {
    w <- 1.
  }

  ## invert covariance matrix of observations

  Valphaxi.objects <- expand( object[["Valphaxi.objects"]] )

  var.signal <- attr( f.reparam.fwd( object[["variogram.object"]] ), "var.signal" )

  G <- var.signal * Valphaxi.objects[["Valphaxi"]]
  G <- G[object[["Tmat"]], object[["Tmat"]]]
  diag( G ) <- diag( G ) + object[["variogram.object"]][[1L]][["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 = 0L ) * iucf )^2 )
    attr( result, "log.det.covmat" ) <- -2. * sum( log( diag( iucf ) ) )
    if( REML ){
      if( !is.null( object[["x"]] ) ){
        XX <- object[["x"]]
      } else {
        if( is.null( object[["model"]]  ) ) stop(
          "'model' component missing in 'object', update 'object' with argument 'model=TRUE"
        )
        XX <- model.matrix( object[["terms"]], object[["model"]] )
      }

      tmp <- t( iucf ) %*% XX
      tmp <- determinant( crossprod( tmp ) )
      attr( result, "log.det.xticovmatx" ) <- tmp[["modulus"]] * tmp[["sign"]]
    }
  }
  result
}


## ##############################################################################

extractAIC.georob <- function( fit, scale = 0, k = 2, ... )
{
  ## extractAIC method for class georob

  ## 2014-02-27 AP
  ## 2016-08-08 AP changes for nested variogram models
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## check arguments

  stopifnot(identical(length(k), 1L) && is.numeric(k) && k > 0)

  tmp <- unlist(lapply(
      fit[["variogram.object"]],
      function(x) c( x[["fit.param"]], x[["fit.aniso"]] )
    ))


  edf <- sum( !is.na( fit[["coefficients"]] ) ) + sum(tmp)
  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, 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
  ## 2015-07-23 AP warning and action for attempt to process Gaussian REML object
  ## 2015-07-23 AP changes for dropping terms from model with fixed variogram parameters
  ## 2015-08-28 AP computation of hessian suppressed
  ## 2016-05-22 AP changes for better computational efficiency
  ## 2016-05-27 AP diagnostics about convergence
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-08-09 AP changes for nested variogram models
  ## 2018-01-05 AP improved memory management in parallel computations
  ##               improved error handling during parallel computations
  ## 2018-08-27 AP elimination of data argument
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## match arguments

  test <- match.arg(test)

  ## check arguments

  stopifnot(identical(length(trace), 1L)            && is.logical(trace))
  stopifnot(identical(length(fixed), 1L)            && is.logical(fixed))
  stopifnot(identical(length(use.fitted.param), 1L) && is.logical(use.fitted.param))

  stopifnot(identical(length(k), 1L)       && is.numeric(k)       && k > 0)
  stopifnot(identical(length(verbose), 1L) && is.numeric(verbose) && verbose >= 0)
  stopifnot(identical(length(ncores), 1L)  && is.numeric(ncores)  && ncores >= 1)


  ## auxiliary function for fitting model and extracting aic by add1 and drop1

  f.aux.add1.drop1 <- function( tt, change, scale, trace, ... ){

    ## objects object, k, n0, verbose are taken from parent environment

    if(trace > 1.) {
      cat( paste( "\ntrying ", change, sep="" ), tt, "\n", sep = "")
      utils::flush.console()
    }

    nfit <- update(
      object, as.formula(paste("~ .", change, tt)),
      verbose = verbose, object. = object
    )

    converged <- TRUE

    if( !is.na(!nfit[["converged"]]) && !nfit[["converged"]] ){
      converged <- FALSE
      warning( "there were errors: call function with argument 'verbose' > 1" )
      if( verbose > 0. ) cat(
        "lack of convergence when fitting model", paste("~ .", change, 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, REML = FALSE, ... ), as.numeric(converged))
    names( df.aic ) <- c("df", "AIC", "converged")
    #   attr( df.aic, "nfit" ) <- list(
    #     variogram.object = nfit[["variogram.object"]],
    #     initial.objects = nfit[["initial.objects"]],
    #     call = nfit[["call"]]
    #   )
    df.aic
  }

  ## evaluate terms and 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" )

  ## manipulate call

  cl <- object[["call"]]

  ## set verbose

  cl <- f.call.set_x_to_value( cl, "verbose", verbose )

  ## update object call to avoid computation of hessian

  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )

  ## update object call to avoid computation of covariance matrices

  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 )

  object[["call"]] <- cl

  ## check if object is result of GAUSSIAN ML fit, manipulate its call and
  ## re-fit it by ML if required

  if(
    object[["tuning.psi"]] >= object[["control"]][["tuning.psi.nr"]] &&
    identical( object[["control"]][["ml.method"]], "REML" )
  ){
    warning( "'object' was estmated by Gaussian REML which cannot be used for ",
      "adding/dropping single terms\n=> 'object' is re-fitted by Gaussian ML ",
      "before adding/deleting terms"
    )

    cl <- object[["call"]]

    cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "ml.method", "ML" )

    object[["call"]] <- cl
    object <- update( object )

  }

  ## manipulate object call if variogram parameters are fixed or estimation
  ## should start from estimated values in object

  if( fixed || use.fitted.param ){

    cl <- f.call.set_allxxx_to_fitted_values( object )

    if( fixed ) cl <- f.call.set_allfitxxx_to_false( cl )

    object[["call"]] <- cl

    object <- update( object )

  }

  ## initialize result

  ns <- length( scope )
  ans <- matrix(
    nrow = ns + 1L, ncol = 3L,
    dimnames = list(c("<none>", scope), c("df", "AIC", "converged"))
  )
  ans[1L,  ] <- c(
    extractAIC( object, scale, k = k, REML = FALSE, ... ),
    object[["converged"]]
  )

  ## loop over all components of scope

  n0 <- nobs( object, use.fallback = TRUE )
  env <- environment( formula(object) )

  ## prepare for parallel execution

  if( ncores > 1L ){
    ncores <- min( ncores, ns, detectCores() )
    trace <- 0
    verbose <- 0.
  }

  ## 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 SNOW cluster on windows OS

    clstr <- makeCluster( ncores, type = "SOCK" )
    save( clstr, file = "SOCKcluster.RData" )
    options( error = f.stop.cluster )

    ## export required items to workers

    junk <- clusterEvalQ( clstr, require( georob, quietly = TRUE ) )
    junk <- clusterExport(
      clstr, c( "object", "k", "n0", "verbose" ), envir =  environment()
    )

    result <- parLapply(
      clstr,
      scope,
      f.aux.add1.drop1,
      change = "+", scale = scale, trace = trace,
      ...
    )

    f.stop.cluster( clstr )

  } else {

    ## fork child processes on non-windows OS

    result <- mclapply(
      scope, f.aux.add1.drop1,
      change = "+", scale = scale, trace = trace,
      mc.cores = ncores,
      mc.allow.recursive = object[["control"]][["pcmp"]][["allow.recursive"]],
      ...
    )

  }

  #   fits <- list(
  #     list(
  #       variogram.object = object[["variogram.object"]],
  #       initial.objects = object[["initial.objects"]],
  #       call = object[["call"]]
  #     )
  #   )
  #   fits[2:NROW(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
  result <- t( simplify2array( result ) )
  ans[2L:NROW(ans), ] <- result
  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], Converged = as.logical( ans[, 3L] ) )

  ## likelihood ratio 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)
    aod <- aod[, c( "Df", "AIC", "LRT", "Pr(>Chi)", "Converged" )]
  }

  ## 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, 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
  ## 2015-07-23 AP warning and action for attempt to process Gaussian REML object
  ## 2015-07-23 AP changes for dropping terms from model with fixed variogram parameters
  ## 2015-08-28 AP computation of hessian suppressed
  ## 2016-05-22 AP changes for better computational efficiency
  ## 2016-05-27 AP diagnostics about convergence
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-08-09 AP changes for nested variogram models
  ## 2018-01-05 AP improved memory management in parallel computations
  ##               improved error handling during parallel computations
  ## 2018-08-27 AP elimination of data argument
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## match arguments

  test <- match.arg(test)

  ## check arguments

  stopifnot(identical(length(trace), 1L)            && is.logical(trace))
  stopifnot(identical(length(fixed), 1L)            && is.logical(fixed))
  stopifnot(identical(length(use.fitted.param), 1L) && is.logical(use.fitted.param))

  stopifnot(identical(length(k), 1L)       && is.numeric(k)       && k > 0)
  stopifnot(identical(length(verbose), 1L) && is.numeric(verbose) && verbose >= 0)
  stopifnot(identical(length(ncores), 1L)  && is.numeric(ncores)  && ncores >= 1)

  ## auxiliary function for fitting model and extracting aic by add1 and drop1

  f.aux.add1.drop1 <- function( tt, change, scale, trace, ... ){

    ## objects object, k, n0, verbose are taken from parent environment

    if(trace > 1.) {
      cat( paste( "\ntrying ", change, sep="" ), tt, "\n", sep = "")
      utils::flush.console()
    }

    nfit <- update(
      object, as.formula(paste("~ .", change, tt)),
      verbose = verbose, object. = object
    )

    converged <- TRUE

    if( !is.na(!nfit[["converged"]]) && !nfit[["converged"]] ){
      converged <- FALSE
      warning( "there were errors: call function with argument 'verbose' > 1" )
      if( verbose > 0. ) cat(
        "lack of convergence when fitting model", paste("~ .", change, 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, REML = FALSE, ... ), as.numeric(converged))
    names( df.aic ) <- c("df", "AIC", "converged")
    #   attr( df.aic, "nfit" ) <- list(
    #     variogram.object = nfit[["variogram.object"]],
    #     initial.objects = nfit[["initial.objects"]],
    #     call = nfit[["call"]]
    #   )
    df.aic
  }

  ## manipulate call

  cl <- object[["call"]]

  ## set verbose

  cl <- f.call.set_x_to_value( cl, "verbose", verbose )

  ## update object call to avoid computation of hessian

  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )

  ## update object call to avoid computation of covariance matrices

  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 )

  object[["call"]] <- cl

  ## check if object is result of GAUSSIAN ML fit, manipulate its call and
  ## re-fit it by ML if required

  if(
    object[["tuning.psi"]] >= object[["control"]][["tuning.psi.nr"]] &&
    identical( object[["control"]][["ml.method"]], "REML" )
  ){
    warning( "'object' was estmated by Gaussian REML which cannot be used for ",
      "adding/dropping single terms\n=> 'object' is re-fitted by Gaussian ML ",
      "before adding/deleting terms"
    )

    cl <- object[["call"]]

    cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "ml.method", "ML" )

    object[["call"]] <- cl
    object <- update( object )

  }

  ## manipulate object call if variogram parameters are fixed or estimation
  ## should start from estimated values in object

  if( fixed || use.fitted.param ){

    cl <- f.call.set_allxxx_to_fitted_values( object )

    if( fixed ) cl <- f.call.set_allfitxxx_to_false( cl )

    object[["call"]] <- cl

    object <- update( object )

  }

  ## evaluate terms and 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 = 3L,
    dimnames = list(c("<none>", scope), c("df", "AIC", "converged"))
  )
  ans[1L,  ] <- c(
    extractAIC( object, scale, k = k, REML = FALSE, ... ),
    object[["converged"]]
  )

  ## loop over all components of scope

  n0 <- nobs( object, use.fallback = TRUE )
  env <- environment( formula(object) )

  ## prepare for parallel execution

  if( ncores > 1L ){
    ncores <- min( ncores, ns, detectCores() )
    trace <- 0
    verbose <- 0.
  }

  ## 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 SNOW cluster on windows OS

    clstr <- makeCluster( ncores, type = "SOCK" )
    save( clstr, file = "SOCKcluster.RData" )
    options( error = f.stop.cluster )

    ## export required items to workers

    junk <- clusterEvalQ( clstr, require( georob, quietly = TRUE ) )
    junk <- clusterExport(
      clstr, c( "object", "k", "n0", "verbose" ), envir =  environment()
    )

    result <- parLapply(
      clstr,
      scope,
      f.aux.add1.drop1,
      change = "-", scale = scale, trace = trace,
      ...
    )

    f.stop.cluster( clstr )

  } else {

    ## fork child processes on non-windows OS

    result <- mclapply(
      scope, f.aux.add1.drop1,
      change = "-", scale = scale, trace = trace,
      mc.cores = ncores,
      mc.allow.recursive = object[["control"]][["pcmp"]][["allow.recursive"]],
      ...
    )

  }

  #   fits <- list(
  #     list(
  #       variogram.object = object[["variogram.object"]],
  #       initial.objects = object[["initial.objects"]],
  #       call = object[["call"]]
  #     )
  #   )
  #   fits[2:NROW(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
  result <- t( simplify2array( result ) )
  ans[2L:NROW(ans), ] <- result
  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[, 2L], Converged = as.logical( ans[, 3L] ) )

  ## likelihood ratio 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)
    aod <- aod[, c( "Df", "AIC", "LRT", "Pr(>Chi)", "Converged" )]
  }

  ## 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, fixed.add1.drop1 = TRUE, fixed.step = fixed.add1.drop1,
  use.fitted.param = TRUE, verbose = 0, ncores = 1, ... )
{

  ## step method for class georob

  ## 2014-06-05 AP
  ## 2015-07-23 AP warning and action for attempt to process Gaussian REML object
  ## 2015-07-23 AP changes for dropping terms from model with fixed variogram parameters
  ## 2016-05-22 AP correction of error when variogram parameters are fixed
  ## 2016-05-22 AP changes for better computational efficiency
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-08-09 AP changes for nested variogram models
  ## 2018-01-05 AP improved memory management in parallel computations
  ## 2018-08-27 AP elimination of data argument
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## match arguments

  direction <- match.arg(direction)

  ## check arguments

  stopifnot(identical(length(use.fitted.param), 1L) && is.logical(use.fitted.param))
  stopifnot(identical(length(fixed.add1.drop1), 1L) && is.logical(fixed.add1.drop1))
  stopifnot(identical(length(fixed.step), 1L)       && is.logical(fixed.step))

  stopifnot(identical(length(trace), 1L)   && is.numeric(trace)   && trace >= 0)
  stopifnot(identical(length(k), 1L)       && is.numeric(k)       && k > 0)
  stopifnot(identical(length(verbose), 1L) && is.numeric(verbose) && verbose >= 0)
  stopifnot(identical(length(ncores), 1L)  && is.numeric(ncores)  && ncores >= 1)
  stopifnot(identical(length(steps), 1L)   && is.numeric(steps)   && steps > 0)

  ## code of step{stats} complemented by argument fixed to control whether
  ## variogram parameters should be kept fixed

  ## auxiliary functions

  mydeviance <- function(x, ...){
    dev <- deviance(x, REML = FALSE, ...)
    if(!is.null(dev)) dev else extractAIC(x, k=0., REML = FALSE, ...)[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
  }

  ## evaluate terms and scope

  Terms <- terms(object)
  object$call$formula <- object$formula <- Terms
  md <- missing(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

  ## store number of fitted variogram parameters

  n.fitted.param.aniso <- sum(
    unlist(
      lapply(
        object[["variogram.object"]],
        function(x) c( x[["fit.param"]], x[["fit.aniso"]] )
      )))


  ## check consistency of arguments

  if( !fixed.add1.drop1 && fixed.step ){
    cat( "\n'fixed.step' set equal to FALSE because fixed.add1.drop1 == FALSE\n\n" )
    fixed.step <- FALSE
  }

  cl <- object[["call"]]

  ## set verbose

  cl <- f.call.set_x_to_value( cl, "verbose", verbose )

  ## update object call to avoid computation of hessian

  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )

  ## update object call to avoid computation of covariance matrices

  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 )

  object[["call"]] <- cl

  ## check if object is result of GAUSSIAN ML fit and manipulate its call to
  ## refit it by ML if required

  if(
    object[["tuning.psi"]] >= object[["control"]][["tuning.psi.nr"]] &&
    identical( object[["control"]][["ml.method"]], "REML" )
  ){
    warning( "'object' was estmated by Gaussian REML which cannot be used for ",
      "adding/dropping single terms\n=> 'object' is re-fitted by Gaussian ML ",
      "before adding/deleting terms"
    )

    cl <- object[["call"]]

    cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "ml.method", "ML" )

    object[["call"]] <- cl
    object <- update( object )

  }

  ## manipulate object call if variogram parameters are fixed or estimation
  ## should start from estimated values in object

  if( fixed.step || use.fitted.param ){

    cl <- f.call.set_allxxx_to_fitted_values( object )

    if( fixed.step ) cl <- f.call.set_allfitxxx_to_false( cl )

    object[["call"]] <- cl

    object <- update( object )

  }

  ## now start the model search

  fit <- object

  bAIC <- extractAIC(fit, scale, k = k, REML = FALSE, ...)
  edf <- bAIC[1L]
  bAIC <- bAIC[2L]
  if( fixed.add1.drop1 && !fixed.step ) bAIC <- bAIC - k * n.fitted.param.aniso
  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 <- 1L
  ## 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 > 0L) {
    steps <- steps - 1L
    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(
        fit, scope$drop, scale = scale,
        k = k, trace = trace >= 1, fixed = fixed.add1.drop1,
        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(
          fit, scope$add, scale = scale,
          k = k, trace = trace >= 1, fixed = fixed.add1.drop1,
          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[-1L, , 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) == 0L) 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" )
    #     )

    fit <- update(
      fit, paste("~ .", change), verbose = verbose, object. = fit
    )

    cl <- object[["call"]]
    cl["formula"] <- fit[["call"]]["formula"]
    fit[["call"]] <- cl

    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, REML = FALSE, ...)
    edf <- bAIC[1L]
    bAIC <- bAIC[2L]
    if( fixed.add1.drop1 && !fixed.step ) bAIC <- bAIC - k * n.fitted.param.aniso
    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 + 1.e-7) break
#     nm <- nm + 1L
    ## 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)
}