R/georob.S3methods.R

Defines functions step.georob step drop1.georob add1.georob safe_pchisq extractAIC.georob deviance.georob logLik.georob waldtest.georob vcov.georob print.summary.georob summary.georob rstandard.georob resid.georob fixef.georob ranef.georob print.georob nobs.georob model.matrix.georob model.frame.georob print.coef.georob coef.georob

Documented in add1.georob coef.georob deviance.georob drop1.georob extractAIC.georob fixef.georob logLik.georob model.frame.georob model.matrix.georob nobs.georob print.coef.georob print.georob print.summary.georob ranef.georob resid.georob rstandard.georob safe_pchisq step step.georob summary.georob vcov.georob waldtest.georob

#####################################
#                                   #
#   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
  
  ## 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
  
  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"
  )]
  
  ans[["terms"]] <- NA
  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"]] ) ){
    
    ## initialization
    
    cor.tf.param <- cov.tf.param <- matrix( 
      NA, nrow = NROW( object[["hessian"]] ), ncol = NROW( object[["hessian"]] ),
      dimnames = dimnames( object[["hessian"]] )
    )
    
#     se <- rep( NA_real_, NROW( object[["hessian"]] ) )
#     names( se ) <- rownames( object[["hessian"]])
#     
#     ci <- matrix( NA_real_, nrow = NROW( object[["hessian"]] ), ncol = 2L )
#     colnames( ci ) <- c( "Lower", "Upper" )
#     rownames( ci ) <- rownames( object[["hessian"]] )
    
    ## select parameters that are not on boundary of parameter space
    
    sr  <- !apply( object[["hessian"]], 1L, function( x ) all( is.na( x ) ) )
    
    if( sum( sr ) > 0L ){
      
      t.chol <- try( chol( object[["hessian"]][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
        )
        
        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

  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
  
  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
  
  ## 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
  
  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
    
  ## 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
  }
    
  test <- match.arg(test)
  
  ## 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
    
  ## 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
  }
    
  test <- match.arg(test)
  
  ## 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

  ## 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)
  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
    
  ## 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, 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, 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)
}

Try the georob package in your browser

Any scripts or data that you put into this service are public.

georob documentation built on Aug. 28, 2018, 1:04 a.m.