R/georob.exported.functions.R

Defines functions profilelogLik default.aniso default.fit.aniso default.fit.param param.bounds param.names expand compress control.pcmp control.nlminb control.optim control.nleqslv control.rq bwd.transf dfwd.transf fwd.transf param.transf control.georob pmm georob

Documented in bwd.transf compress control.georob control.nleqslv control.nlminb control.optim control.pcmp control.rq default.aniso default.fit.aniso default.fit.param dfwd.transf expand fwd.transf georob param.bounds param.names param.transf pmm profilelogLik

georob <- 
  function( 
    formula, data, subset, weights, na.action,
    model = TRUE, x = FALSE, y = FALSE, 
    contrasts = NULL, offset, 
    locations,
    variogram.model = c( "RMexp", "RMaskey", "RMbessel", "RMcauchy", 
      "RMcircular", "RMcubic", "RMdagum", "RMdampedcos", "RMdewijsian", "RMfbm",
      "RMgauss", "RMgencauchy", "RMgenfbm", "RMgengneiting", "RMgneiting", "RMlgd",
      "RMmatern", "RMpenta", "RMqexp", "RMspheric", "RMstable",
      "RMwave", "RMwhittle"
    ), 
    param, fit.param = default.fit.param()[names(param)],
	  aniso = default.aniso(), fit.aniso = default.fit.aniso(),
    variogram.object = NULL,
    tuning.psi = 2., control = control.georob(), verbose = 0,
    ...
  )
{
  
  ## wrapper function to georob.fit with "standard" interface for
  ## statistical modelling
  
  ## ToDos:
  
  ## - Variante fuer Klasse SpatialPointsDataFrame{sp}
  ## - Ueberpruefung der Konsistenz des Inputs hier statt in georob.fit
  ## - ausgewaehlte Elemente von control in Resultate speichern
  
  
  ## History:
  
  ## georob:
  
  ## - Implementierung von Schaetzung von Startwerten der Kovarianzparameter mittels REML
  ##   von um Ausreisser bereinigten Datensatz
  ## - Implementierung von Schaetzung von Nugget und Varianz von mikroskaliger Variation (snugget )
  ## - Implementierung der Schaetzungen fuer replizierte Messungen an gleicher Messtelle    
  ## - neue Varianten um Startwerte von betahat zu rechnen (Ersatz von use.lm durch 
  ##   Argument initial.fixef = c( "rq", "lmrob", "rlm" ), vgl. control.georob)
  ## - Kontrolle, ob Designmatrix vollen Spaltenrang hat
  ## - Modifikation fuer Fall, dass alle Variogrammparameter fixiert sind
  ## - IRWLS Berechung von betahat und bhat entweder von Werten im initial.object 
  ##   oder von Schaetzwerten aus vorangehender Iteration
  ## - Berechung der Kovarianzen zwischen betahat und (B - bhat)
  ## - Steuerung der Berechnung der verschiedenen Kovarianzen via control.georob
  ## - vollstaendige Implementierung von Standard Interfaces fuer Input und
  ##   Output fuer statistische Modellierung (analog lm, lmrob)
  ## - korrekte Behandlung von NAs und Implementierung von subset
  ## - Verzicht auf Berechnung von initialem Variogramm
  ## - Startwerte fuer bhat alle gleich Null gesetzt
  ## - neue Struktur von initial.objects
  ## - teilweise geaenderte Namen der Argumente (model -> variogram.model)
  
  ## f.glsrob803:
  
  ## - Umbenennung einiger Argument
  
  ## f.glsrob801:
  
  ## - teilweise Implementation von Standard Interfaces fuer Input und
  ##   Output fuer statistische Modellierung 
  ## - Berechnung des Objekts mit den Startwerten
  
  
  ## 2012-04-21 AP
  ## 2012-05-03 AP bounds for safe parameter values
  ## 2012-05-07 AP correction of error for constant trend
  ## 2012-05-28 AP handle missing names of coefficients after calling update
  ## 2013-04-23 AP new names for robustness weights
  ## 2013-05-23 AP correct handling of missing observations and to construct model.frame
  ## 2013-06-03 AP handling design matrices with rank < ncol(x)
  ## 2013-06-12 AP substituting [["x"]] for $x in all lists
  ## 2013-07-02 AP new transformation of rotation angles
  ## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
  ## 2013-09-06 AP exclusive use of nleqslv for solving estimating equations
  ## 2014-02-18 AP correcting error when fitting models with offset
  ## 2014-05-15 AP changes for version 3 of RandomFields
  ## 2014-05-22 AP correcting error when selecting initial.param
  ## 2014-06-02 AP partial matching of names of variogram parameters
  ## 2014-08-18 AP changes for parallelized computations
  ## 2014-08-18 AP changes for Gaussian ML estimation
  ## 2014-08-26 AP changes to return ml.method if fitted object
  ## 2015-06-30 AP function and arguments renamed
  ## 2015-07-17 AP change of control of computation of hessian for Gaussian (RE)ML, 
  ##               changes for singular design matrices, nlminb optimizer added
  ## 2015-08-19 AP control about error families for computing covariances added
  ## 2015-08-28 AP computation of hessian suppressed; correction of error when using georob.object;
  ##               control arguments hessian, initial.param, initial.fixef newly organized
  ## 2015-11-25 AP new way to control which variogram parameters are fitted
  ## 2016-07-15 AP allowing use of lmrob for computing initial fixed effects for rank-deficient model matrix
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-08-08 AP changes for nested variogram models
  ## 2017-05-07 AP correcting error for interface to rq.fit
  
  ## check validity of tuning.psi
  
  if( identical( control[["psi.func"]], "t.dist" ) && tuning.psi <= 1. ) 
    stop( "'tuning.psi' must be greater than 1 for t-dist psi-function" )
    
  ## check whether all mandatory arguments have been provided
  
  if( missing( formula ) || missing( locations ) || 
    ( is.null( variogram.object ) && missing( param )) ) stop( 
    "some mandatory arguments are missing" 
  )
    
  ## setup or check contents of variogram.object
  
  if( is.null( variogram.object ) ){
    
    ## match variogram model
    
    variogram.model <- match.arg( variogram.model )  
        
    ## match names of param, aniso, fit.param, fit.aniso
    ## !! CARE: code works only if fit.param has not been used before (lazy evaluation)
  
    if( !missing( param ) ){
      tmp <- names( param )
      tmp <- sapply(tmp, function(x, choices){
          match.arg(x, choices)
        },
        choices = names( default.fit.param() )
      )
      names( param ) <- tmp
    }
    
    if( !missing( fit.param ) ){
      tmp <- names( fit.param )
      tmp <- sapply(tmp, function(x, choices){
          match.arg(x, choices)
        },
        choices = names( default.fit.param() )
      )
      names( fit.param ) <- tmp
      fit.param <- fit.param[names( fit.param ) %in% names( param )]
    }
    
    if( !missing( aniso ) ){
      tmp <- names( aniso )
      tmp <- sapply(tmp, function(x, choices){
          match.arg(x, choices)
        },
        choices = names( default.aniso() )
      )
      names( aniso ) <- tmp
    }
    
    if( !missing( fit.aniso ) ){
      tmp <- names( fit.aniso )
      tmp <- sapply(tmp, function(x, choices){
          match.arg(x, choices)
        },
        choices = names( default.aniso() )
      )
      names( fit.aniso ) <- tmp
    }
    
    ## create variogram.object
    
    variogram.object <- list(
      list( 
        variogram.model = variogram.model,
        param = param, fit.param = fit.param,
        aniso = aniso, fit.aniso = fit.aniso
      )
    )
  
  } else {
    
    if( !missing( param ) ) cat( 
      "\n information on initial parameters in 'param' ignored because 'variogram.object' specified\n\n"
    )
    if( !missing( fit.param ) ) cat( 
      "\n information on initial parameters in 'fit.param' ignored because 'variogram.object' specified\n\n"
    )
    if( !missing( aniso ) ) cat( 
      "\n information on initial parameters in 'aniso' ignored because 'variogram.object' specified\n\n"
    )
    if( !missing( fit.aniso ) ) cat( 
      "\n information on initial parameters in 'fit.aniso' ignored because 'variogram.object' specified\n\n"
    )
    
    variogram.object <- lapply(
      variogram.object,
      function( y, vm ){
            
        ## match names of components
        
        tmp <- names( y )
        tmp <- sapply(tmp, function(x, choices){
            match.arg(x, choices)
          },
          choices = c( "variogram.model", "param", "fit.param", "aniso", "fit.aniso" )
        )
        names( y ) <- tmp
        
        ## check whether component variogram.model and param are present
        
        if( is.null( y[["variogram.model"]] ) ) stop(
          "'variogram.model' missing in some component of variogram.object"       
        )
        if( is.null( y[["param"]] ) ) stop(
          "'param' missing in some component of variogram.object"       
        )
        
        ## match variogram.model
        
        y[["variogram.model"]] <- match.arg( y[["variogram.model"]], vm )
        
                
        ## match names of param, aniso, fit.param, fit.aniso and set
        ## default values if missing
        
        tmp <- names( y[["param"]] )
        tmp <- sapply(tmp, function(x, choices){
            match.arg(x, choices)
          },
          choices = names( default.fit.param() )
        )
        names( y[["param"]] ) <- tmp
        
        if( is.null( y[["fit.param"]] ) ){
          y[["fit.param"]] <- default.fit.param()
        } else {
          tmp <- names( y[["fit.param"]] )
          tmp <- sapply(tmp, function(x, choices){
              match.arg(x, choices)
            },
            choices = names( default.fit.param() )
          )
          names( y[["fit.param"]] ) <- tmp
        }
        y[["fit.param"]] <- y[["fit.param"]][names(y[["fit.param"]] ) %in% names( y[["param"]] )]
        
        if( is.null( y[["aniso"]] ) ){
          y[["aniso"]] <- default.aniso()
        } else {
          tmp <- names( y[["aniso"]] )
          tmp <- sapply(tmp, function(x, choices){
              match.arg(x, choices)
            },
            choices = names( default.aniso() )
          )
          names( y[["aniso"]] ) <- tmp
        }
        
        if( is.null( y[["fit.aniso"]] ) ){
          y[["fit.aniso"]] <- default.fit.aniso()
        } else {
          tmp <- names( y[["fit.aniso"]] )
          tmp <- sapply(tmp, function(x, choices){
              match.arg(x, choices)
            },
            choices = names( default.aniso() )
          )
          names( y[["fit.aniso"]] ) <- tmp
        }
        y
      }, vm = variogram.model
    )
  }
  
  ## get model frame, response vector, weights, offset and design
  ## matrix (cf.  lm, lmrob)
  
  ret.x <- x
  ret.y <- y
  
#   ## vector with row number of included observations
#   
#   in.subset <- 1L:NROW( data )
#   if( !missing( subset ) ) in.subset <- in.subset[subset]
  
  ## build combined formula for fixed effects and locations
  
  extended.formula <- update( 
    formula,
    paste( 
      paste( as.character( formula )[c(2L, 1L, 3L)], collapse = " " ),
      as.character( locations )[2L], sep = " + "
    )
  )
  
  ## setting-up model frame
  
  cl <- match.call()
  mf <- match.call( expand.dots = FALSE )
  m <- match( 
    c( "formula", "data", "subset", "weights", "na.action", "offset"),
    names(mf), 0L 
  )
  mf <- mf[c(1L, m)]
  mf[["formula"]] <- extended.formula
  mf[["drop.unused.levels"]] <- TRUE
  mf[[1L]] <- as.name( "model.frame" )
  
  mf <- eval( mf, parent.frame() )
  
  ## setting-up terms objects
  
  mt     <- terms( formula )
  mt.loc <- terms( locations )
  
  ## eliminate intercept from mt.loc
  
  attr( mt.loc, "intercept" ) <- 0L
    
#   ## ... and assign fixed effects terms object as attribute to model.frame
#   
#   attr( mf, "terms" ) <- mt
  
  ## check whether 'empty' models have been entered
  
  if( is.empty.model( mt ) )
    stop( "an 'empty' fixed effects model has been specified" )
  if( is.empty.model( mt.loc ) )
    stop( "an 'empty' locations model has been specified" )
  
  ## check whether fixed effects model includes an intercept if an
  ## intrinsic variogram model is used
  
  if( identical( attr( mt, "intercept" ), 0L ) && 
    variogram.model %in% control[["irf.model"]] )
  stop(
    "the fixed effects model must include an intercept ",
    "if an unbounded variogram model is used"
  )
  
  ## extract fixed effects response variable, weights, offset and design matrix
  
  y <- model.response( mf, "numeric" )
  
  w <- as.vector( model.weights( mf ) )
  if( !is.null(w) )
    stop( "weights are not yet implemented for this estimator" )
  
  offset <- as.vector( model.offset(mf) )
  if( !is.null(offset) ) {
    if( length( offset ) != NROW(y) )
      stop( gettextf(
        "number of offsets is %d, should equal %d (number of observations)", 
        length(offset), NROW(y) ), domain = NA )
  }
  
  x <- model.matrix( mt, mf, contrasts )
  
  ## check if optionally provided bhat has correct length
  
  if( !is.null( control[["bhat"]] ) && length( y ) != length( control[["bhat"]] ) ) stop(
    "lengths of response vector and 'bhat' do not match"    
  )
  
  ## adjust initial.param if all variogram parameters are fixed
  
  if( all( sapply(
        variogram.object,
        function( x ) !any( c( x[["fit.param"]], x[["fit.aniso"]] ) )
      ))) control[["initial.param"]] <- FALSE
  
  ## adjust choice for initial.fixef to compute regression coefficients
  
  if( tuning.psi < control[["tuning.psi.nr"]] ){
    if( control[["initial.param"]] ) control[["initial.fixef"]] <- "lmrob"
  } else {
    control[["initial.param"]] <- FALSE
  }
  
  ## check whether design matrix has full column rank
  
  col.rank.XX <- list( deficient = FALSE, rank = NCOL( x ) )
  
  sv <- svd( crossprod( x ) )[["d"]]
  min.max.sv <- range( sv )
  condnum <- min.max.sv[1L] / min.max.sv[2L] 
  
  if( condnum <= control[["min.condnum"]] ){
    col.rank.XX[["deficient"]] <- TRUE
    col.rank.XX[["col.rank"]] <- sum( sv / min.max.sv[2L] > control[["min.condnum"]] )
    if( control[["initial.fixef"]] == "rq" ){
      cat( 
        "design matrix has not full column rank (condition number of X^T X: ", 
        signif( condnum, 2L ), ")\ninitial values of fixed effects coefficients are computed by 'lmrob' instead of 'rq'\n\n"
      )
      control[["initial.fixef"]] <- "lmrob"
      warning( 
        "design matrix has not full column rank (condition number of X^T X: ", 
        signif( condnum, 2L ), ")\ninitial values of fixed effects coefficients are computed by 'lmrob' instead of 'rq'"
      )
    }
  }
  
  ## subtract offset
  
  yy <- y
  if( !is.null( offset ) ) yy <- yy - offset
  
  ## compute initial guess of fixed effects parameters (betahat)
  
  r.initial.fit <- switch(
    control[["initial.fixef"]],
    rq = {
      
      Rho <- function( u, tau) u * (tau - (u < 0))
      tau <- control[["rq"]][["tau"]]
      process <- (tau < 0 || tau > 1)
      
      if(tau == 0) tau <- control[["req"]][["eps"]]
      if(tau == 1) tau <- 1. - control[["req"]][["eps"]]
      
      fit <- switch(
        control[["rq"]][["method"]],
        br = rq.fit( 
          x = x, y = yy, tau = tau, method = control[["rq"]][["method"]], 
          alpha = control[["rq"]][["alpha"]], ci = control[["rq"]][["ci"]],
          iid = control[["rq"]][["iid"]], interp = control[["rq"]][["interp"]], 
          tcrit = control[["rq"]][["tcrit"]]
        ),
        fnb = rq.fit( 
          x = x, y = yy, tau = tau, method = control[["rq"]][["method"]], 
          beta = control[["rq"]][["beta"]], eps = control[["rq"]][["eps"]]
        ),
        pfn = rq.fit( 
          x = x, y = yy, tau = tau, method = control[["rq"]][["method"]], 
          Mm.factor = control[["rq"]][["Mm.factor"]], 
          max.bad.fixup = control[["rq"]][["max.bad.fixup"]]
        )
      )
       
      if (process){
        rho <- list(x = fit$sol[1, ], y = fit$sol[3, ])
      } else {
        if (length(dim(fit$residuals))) 
        dimnames(fit$residuals) <- list(dimnames(x)[[1]], 
          NULL)
        rho <- sum(Rho(fit$residuals, tau))
      }
      if( control[["rq"]][["method"]] == "lasso" ){
        class(fit) <- c("lassorq", "rq")
      } else if( control[["rq"]][["method"]] == "scad"){
        class(fit) <- c("scadrq", "rq")
      } else {
        class(fit) <- ifelse(process, "rq.process", "rq")
      }
      fit[["na.action"]] <- attr( mf, "na.action" )
      fit[["formula"]] <- formula
      fit[["terms"]] <- mt
      fit[["xlevels"]] <- .getXlevels(mt, mf)
      fit[["call"]] <- cl
      fit[["tau"]] <- tau
      fit[["weights"]] <- w
      fit[["residuals"]] <- drop( fit[["residuals"]] )
      fit[["rho"]] <- rho
      fit[["method"]] <- control[["rq"]][["method"]]
      fit[["fitted.values"]] <- drop( fit[["fitted.values"]] )
      attr(fit, "na.message") <- attr( m, "na.message" )
      if( model ) fit[["model"]] <- mf
      fit
      
    },
    lmrob = {
      
      fit <- lmrob.fit( x, yy, control = control[["lmrob"]] )
      fit[["na.action"]] <- attr(mf, "na.action")
      fit[["offset"]] <- offset
      fit[["contrasts"]] <- attr(x, "contrasts")
      fit[["xlevels"]] <- .getXlevels(mt, mf)
      fit[["call"]] <- cl
      fit[["terms"]] <- mt
      if( control[["lmrob"]][["compute.rd"]] && !is.null(x) )
      fit[["MD"]] <- robMD( x, attr( mt, "intercept" ) )
      if( !is.null( offset ) ) fit[["fitted.values"]] + offset
      fit
      
    },
    lm = {
      
      fit <- if( is.null(w) ){
        lm.fit(x, y, offset = offset, singular.ok = TRUE )
      } else {
        lm.wfit(x, y, w, offset = offset, singular.ok = TRUE )
      }
      class(fit) <- c(if (is.matrix(y)) "mlm", "lm")
      fit[["na.action"]] <- attr(mf, "na.action")
      fit[["offset"]] <- offset
      fit[["contrasts"]] <- attr(x, "contrasts")
      fit[["xlevels"]] <- .getXlevels(mt, mf)
      fit[["call"]] <- cl
      fit[["terms"]] <- mt
      if (model) fit[["model"]] <- mf
      if (ret.x) fit[["x"]] <- x
      if (ret.y) fit[["y"]] <- y
      fit[["qr"]] <- NULL
      fit
      
    }
  )
  
  ## compute coordinates of locations and distance object
  
  locations.coords <- model.matrix( mt.loc, mf )
  
  if( 
    !( missing( aniso ) || missing( fit.aniso ) ) && 
    ( NCOL( locations.coords ) < 2L || NCOL( locations.coords ) > 3L )
  ) stop( 
    "anisotropic variogram models are implemented only for 2 or 3 dimensions" 
  )

  names( yy ) <- rownames( mf )
  
  ##  check whether argument "object."  has been provided in call (e.g. by
  ##  update ) and extract'locations' exists in workspace
  
  extras <- match.call( expand.dots = FALSE )$...
  georob.object <- extras[names(extras) %in% "object."]
  if( 
    length( georob.object ) && 
    exists( as.character( georob.object ), envir = parent.frame() ) #&&
    #     !control[["initial.param"]]
  ){
    if( verbose > 4. ) cat(
      "\n    georob: using some components of 'object.'\n"    
    )
    georob.object <- eval( 
      georob.object[[1L]], parent.frame() 
    )[c( "variogram.object", "locations.objects", "Valphaxi.objects" )]
    georob.object[["Valphaxi.objects"]] <- c( 
      list(error=FALSE), 
      georob.object[["Valphaxi.objects"]] 
    )
  } else {
    georob.object <- NULL
  }
  
  ## compute initial values of variogram and anisotropy parameters
  
  if( tuning.psi < control[["tuning.psi.nr"]] && control[["initial.param"]] ){
    
    if( verbose > 0. ) cat( "\ncomputing robust initial parameter estimates ...\n" )
    
    t.sel <- r.initial.fit[["rweights"]] > control[["min.rweight"]]
    
    if( verbose > 0. ) cat( 
      "\ndiscarding", sum( !t.sel ), "of", length( t.sel ), 
      "observations for computing initial estimates of variogram\nand anisotropy parameter by gaussian reml\n"
    )
    
    ## collect.items for initial object
    
    initial.objects <- list(
      x = as.matrix( x[t.sel, ] ),
      y = yy[t.sel],
      betahat = coef( r.initial.fit ),
      bhat = if( is.null( control[["bhat"]] ) ){
        rep( 0., length( yy ) )[t.sel]
      } else {
        control[["bhat"]][t.sel]
      },
      initial.fit = r.initial.fit,
      locations.objects = list(
        locations = locations,
        coordinates = locations.coords[t.sel, ]
      )
    )
    
    ## estimate model parameters with pruned data set
    
    t.georob <- georob.fit(
      initial.objects = initial.objects,
      variogram.object = variogram.object,
      param.tf = control[["param.tf"]],
      fwd.tf = control[["fwd.tf"]], 
      deriv.fwd.tf = control[["deriv.fwd.tf"]],
      bwd.tf = control[["bwd.tf"]], 
      georob.object = georob.object,
      safe.param = control[["safe.param"]],
      tuning.psi = control[["tuning.psi.nr"]],
      error.family.estimation = control[["error.family.estimation"]],
      error.family.cov.effects = control[["error.family.cov.effects"]],
      error.family.cov.residuals = control[["error.family.cov.residuals"]],
      cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
      cov.betahat = control[["cov.betahat"]], 
      cov.bhat.betahat = control[["cov.bhat.betahat"]],
      cov.delta.bhat = control[["cov.delta.bhat"]],
      full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
      cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
      cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
      cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
      aux.cov.pred.target = control[["aux.cov.pred.target"]],
      min.condnum = control[["min.condnum"]],
      col.rank.XX = col.rank.XX,
      psi.func = control[["psi.func"]],
      tuning.psi.nr = tuning.psi,
      ml.method = control[["ml.method"]],
      maximizer = control[["maximizer"]],
      reparam = control[["reparam"]],
      irwls.initial = control[["irwls.initial"]],
      irwls.maxiter = control[["irwls.maxiter"]],
      irwls.ftol = control[["irwls.ftol"]],
      force.gradient = control[["force.gradient"]],
      zero.dist = control[["zero.dist"]],
      control.nleqslv =  control[["nleqslv"]],
      control.optim = control[["optim"]],
      control.nlminb = control[["nlminb"]],
      hessian = FALSE,
      control.pcmp = control[["pcmp"]],
      verbose = verbose
    )
    
    variogram.object <- lapply(
      1:length(variogram.object),
      function( i, x, fvo ){
        x   <- x[[i]]
        fvo <- fvo[[i]]
        x[["param"]] <- fvo[["param"]][names(x[["param"]])]
        x[["aniso"]] <- fvo[["aniso"]][names(x[["aniso"]])]
        x
      }, x = variogram.object, fvo = t.georob[["variogram.object"]]
    )
    
  }
  
  ## collect.items for initial object
  
  initial.objects <- list(
    x = as.matrix( x ),
    y = yy,
    betahat = coef( r.initial.fit ),
    bhat = if( is.null( control[["bhat"]] ) ){
      rep( 0., length( yy ) )
    } else {
      control[["bhat"]]
    },
    initial.fit = r.initial.fit,
    locations.objects = list(
      locations = locations,
      coordinates = locations.coords
    )
  )
  
  ## estimate model parameters

  if( verbose > 0. ) cat( "\ncomputing final parameter estimates ...\n" )

  r.georob <- georob.fit(
    initial.objects = initial.objects,
    variogram.object = variogram.object,
    param.tf = control[["param.tf"]],
    fwd.tf = control[["fwd.tf"]], 
    deriv.fwd.tf = control[["deriv.fwd.tf"]],
    bwd.tf = control[["bwd.tf"]], 
    georob.object = georob.object,
    safe.param = control[["safe.param"]],
    tuning.psi = tuning.psi,
    error.family.estimation = control[["error.family.estimation"]],
    error.family.cov.effects = control[["error.family.cov.effects"]],
    error.family.cov.residuals = control[["error.family.cov.residuals"]],
    cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
    cov.betahat = control[["cov.betahat"]], 
    cov.bhat.betahat = control[["cov.bhat.betahat"]],
    cov.delta.bhat = control[["cov.delta.bhat"]],
    full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
    cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
    cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
    cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
    aux.cov.pred.target = control[["aux.cov.pred.target"]],
    min.condnum = control[["min.condnum"]],
    col.rank.XX = col.rank.XX,
    psi.func = control[["psi.func"]],
    tuning.psi.nr = control[["tuning.psi.nr"]],
    ml.method = control[["ml.method"]],
    maximizer = control[["maximizer"]],
    reparam = control[["reparam"]],
    irwls.initial = control[["irwls.initial"]],
    irwls.maxiter = control[["irwls.maxiter"]],
    irwls.ftol = control[["irwls.ftol"]],
    force.gradient = control[["force.gradient"]],
    zero.dist = control[["zero.dist"]],
    control.nleqslv =  control[["nleqslv"]],
    control.optim = control[["optim"]],
    control.nlminb = control[["nlminb"]],
    hessian = control[["hessian"]],
    control.pcmp = control[["pcmp"]],
    verbose = verbose
  )
  
  ## add offset to fitted values
  
  if( !is.null( offset ) )
    r.georob[["fitted.values"]] <- r.georob[["fitted.values"]] + offset
    
  ##
    
  r.georob[["control"]] <- control
  
  ## add remaining items to output
  
  if( control[["lmrob"]][["compute.rd"]] && !is.null( x ) )
    r.georob[["MD"]] <- robMD( x, attr(mt, "intercept") )

  if( model ) r.georob[["model"]] <- mf
  if( ret.x ) r.georob[["x"]] <- x
  if( ret.y ) r.georob[["y"]] <- y

  r.georob[["df.residual"]] <- length(yy) - col.rank.XX[["rank"]]
  r.georob[["na.action"]] <- attr(mf, "na.action")
  r.georob[["offset"]] <- offset
  r.georob[["contrasts"]] <- attr(x, "contrasts")
  r.georob[["xlevels"]] <- .getXlevels(mt, mf)
  r.georob[["rank"]] <- col.rank.XX[["rank"]]
  r.georob[["call"]] <- cl
  r.georob[["terms"]] <- mt
  
  ## set missing names of coefficients (bug of update)
  
  if( length( r.georob[["coefficients"]] ) == 1L && is.null( names( r.georob[["coefficients"]] ) ) ){
    names( r.georob[["coefficients"]] ) <- "(Intercept)"
  }

  
  class( r.georob ) <- c( "georob" )
  
  invisible( r.georob )
  
}


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

pmm <- 
  function( 
    A, B, control = control.pcmp() 
  )
{
  
  ## function for parallelized matrix multiplication inspired by function
  ## parMM{snow}
  
  ## 2014-06-25 A. Papritz
  ## 2015-03-13 AP small changes in f.aux
  ## 2016-07-20 AP arguments renamed, new control of recursive parallelization
  ## 2018-01-05 AP improved memory management in parallel computations
  
  ## auxiliary function 
  
  f.aux <- function( i ){
    
    ## s, e, A, B are taken from parent environment
    
    A %*% B[ , s[i]:e[i], drop = FALSE]
    
  }
  
  ## determine number of cores
  
  ncores <- control[["pmm.ncores"]]
  if( !control[["allow.recursive"]] ) ncores <- 1L
  
  ## determine columns indices of matrix blocks
  
  k <- control[["f"]] * ncores
  n <- NCOL(B)
  dn <- floor( n / k )
  s <- ( (0L:(k-1L)) * dn ) + 1L
  e <- (1L:k) * dn
  e[k] <- n
  
  ##
  
  ## set default value for control of forking if missing (required for backward compatibility)
  
  if( is.null( control[["fork"]] ) ){
    control[["fork"]] <- !identical( .Platform[["OS.type"]], "windows" )
  }
    
  if( ncores > 1L ){
    
    if( !control[["fork"]] ){
      
      ## create a SNOW cluster on windows OS
      
      if( !sfIsRunning() ){
        options( error = f.stop.cluster )
        junk <- sfInit( parallel = TRUE, cpus = ncores )
      }
      
      junk <- sfExport( "s", "e", "A", "B" )
      
      res <- sfLapply( 1L:k, f.aux )
      
      if( control[["sfstop"]] ){
        junk <- sfStop()
        options( error = NULL )
      }
      
    } else {
      
      res <- mclapply( 1L:k, f.aux, mc.cores = ncores )
      
    }
    
    junk <- gc()
    
    matrix( unlist(res), nrow = NROW( A ) )
    
  } else A %*% B    
  
}

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

control.georob <- 
  function(
    ml.method = c( "REML", "ML" ),
    reparam = TRUE,
    maximizer = c( "nlminb", "optim" ),
    initial.param = TRUE,
    initial.fixef = c("lmrob", "rq", "lm"),
    bhat = NULL,
    min.rweight = 0.25,
    param.tf = param.transf(),
    fwd.tf = fwd.transf(), 
    deriv.fwd.tf = dfwd.transf(), 
    bwd.tf = bwd.transf(),
    psi.func = c( "logistic", "t.dist", "huber" ),
    irwls.maxiter = 50, irwls.ftol = 1.e-5,
    force.gradient = FALSE,
    min.condnum = 1.e-12,
    zero.dist = sqrt( .Machine[["double.eps"]] ),
    error.family.estimation    = c( "gaussian", "long.tailed" ),
    error.family.cov.effects   = c( "gaussian", "long.tailed" ),
    error.family.cov.residuals = c( "gaussian", "long.tailed" ),
    cov.bhat = TRUE, full.cov.bhat = FALSE,
    cov.betahat = TRUE, 
    cov.delta.bhat = TRUE, full.cov.delta.bhat = TRUE,
    cov.delta.bhat.betahat = TRUE,
    cov.ehat = TRUE, full.cov.ehat = FALSE,
    cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
    hessian = TRUE,
    rq = control.rq(),
    lmrob = lmrob.control(),
    nleqslv = control.nleqslv(),
    optim = control.optim(),
    nlminb = control.nlminb(),
    pcmp = control.pcmp(),
    ...
  )
{
  
  ## auxiliary function to set meaningful default values for georob

  ## Arguments: 
  
  ## 2012-04-21 A. Papritz   
  ## 2012-05-03 AP bounds for safe parameter values
  ## 2012-05-04 AP modifications for lognormal block kriging
  ## 2013-04-23 AP new names for robustness weights
  ## 2013-06-12 AP changes in stored items of Valphaxi object
  ## 2013-06-12 AP substituting [["x"]] for $x in all lists
  ## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
  ## 2014-05-15 AP changes for version 3 of RandomFields
  ## 2014-08-18 AP changes for parallelized computations
  ## 2014-08-18 AP changes for Gaussian ML estimation
  ## 2015-03-10 AP extended variogram parameter transformations
  ## 2015-06-30 AP function and arguments renamed
  ## 2015-07-17 AP nlminb optimizer added
  ## 2015-08-19 AP control about error families for computing covariances added
  ## 2015-08-28 AP control arguments hessian, initial.param, initial.fixef newly organized
  ## 2016-07-15 AP some arguments hard coded
  ## 2016-08-24 AP new default value for error.family.cov.residuals
  
  irwls.initial <- TRUE
  tuning.psi.nr <- 1000.
  cov.bhat.betahat <- FALSE
  aux.cov.pred.target <- FALSE
  safe.param <- 1.e12
  sepstr <- ".__...__."
  
  if( 
    !( all( unlist( param.tf ) %in% names( fwd.tf ) ) &&
       all( unlist( param.tf ) %in% names( deriv.fwd.tf ) ) &&
       all( unlist( param.tf ) %in% names( bwd.tf ) )
    )
  ) stop( 
    "undefined transformation of variogram parameters; extend respective function definitions" 
  )
  
  if( !irwls.initial && irwls.ftol >= 1.e-6 ) warning( 
    "'irwls.initial == FALSE' and large 'ftol' may create problems for root finding"
  )
  
  list(
    ml.method = match.arg( ml.method ), reparam = reparam,
    maximizer = match.arg( maximizer ),
    initial.param = initial.param,
    initial.fixef = match.arg( initial.fixef ),
    bhat = bhat,
    min.rweight = min.rweight,
    param.tf = param.tf, fwd.tf = fwd.tf, deriv.fwd.tf = deriv.fwd.tf, bwd.tf = bwd.tf,
    safe.param = safe.param, sepstr = sepstr,
    psi.func = match.arg( psi.func ), 
    tuning.psi.nr = tuning.psi.nr,
    irwls.initial = irwls.initial, irwls.maxiter = irwls.maxiter, irwls.ftol = irwls.ftol,
    force.gradient = force.gradient,
    min.condnum = min.condnum,
    zero.dist = zero.dist,
    error.family.estimation = match.arg( error.family.estimation ),
    error.family.cov.effects = match.arg( error.family.cov.effects ),
    error.family.cov.residuals = match.arg( error.family.cov.residuals ),    
    cov.bhat = cov.bhat, full.cov.bhat = full.cov.bhat,
    cov.betahat = cov.betahat, 
    cov.bhat.betahat = cov.bhat.betahat,
    cov.delta.bhat = cov.delta.bhat, full.cov.delta.bhat = full.cov.delta.bhat,
    cov.delta.bhat.betahat = cov.delta.bhat.betahat,
    cov.ehat = cov.ehat, full.cov.ehat = full.cov.ehat,
    cov.ehat.p.bhat = cov.ehat.p.bhat, full.cov.ehat.p.bhat = full.cov.ehat.p.bhat,
    aux.cov.pred.target = aux.cov.pred.target,
    hessian = hessian,
    irf.models = c( "RMdewijsian", "RMfbm", "RMgenfbm" ),
    rq = rq, lmrob = lmrob, nleqslv = nleqslv, 
    optim = optim, 
    nlminb = nlminb,
    pcmp = pcmp
  )
  
}

## ======================================================================
param.transf <-
  function( 
    variance = "log", snugget = "log", nugget = "log", scale = "log", 
    alpha = c( 
      RMaskey = "log", RMdewijsian = "logit2", RMfbm = "logit2", RMgencauchy = "logit2", 
      RMgenfbm = "logit2", RMlgd = "identity", RMqexp = "logit1", RMstable = "logit2"
    ), 
    beta = c( RMdagum = "logit1", RMgencauchy = "log", RMlgd = "log" ), 
    delta = "logit1", 
    gamma = c( RMcauchy = "log", RMdagum = "logit1" ), 
    kappa = "logit3", lambda = "log", 
    mu = "log", 
    nu = "log",
    f1 = "log", f2  ="log", omega = "identity", phi = "identity", zeta = "identity"
  )
{
  
  ## function sets meaningful defaults for transformation of variogram
  ## parameters
  
  ## 2013-07-02 A. Papritz
  ## 2014-05-15 AP changes for version 3 of RandomFields
  ## 2015-03-10 AP extended transformation
  ## 2015-04-07 AP changes for fitting anisotropic variograms
  
  list( 
    variance = variance, snugget = snugget, nugget = nugget, scale = scale,
    alpha = alpha, 
    beta = beta, 
    delta = delta, 
    gamma = gamma, 
    kappa = kappa, 
    lambda = lambda, 
    mu = mu, 
    nu = nu, 
    f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta
  )
  
}

## ======================================================================
fwd.transf <- 
  function( 
    ...
  )
{
  
  ## definition of forward transformation of variogram parameters
  
  ## 2013-07-02 A. Papritz
  ## 2015-03-10 AP extended variogram parameter transformations
  ## 2015-04-07 AP changes for fitting anisotropic variograms
  
  list( 
    log = function(x) log(x), 
    logit1 = function(x) log( x / (1. - x) ), 
    logit2 = function(x) log( x / (2. - x) ), 
    logit3 = function(x) log( (x - 1.) / (3. - x) ), 
    identity = function(x) x, ...
  )
}

## ======================================================================
dfwd.transf<-
  function( 
    ...
  )
{
  
  ## definition of first derivative of forward transformation of variogram
  ## parameters
  
  ## 2013-07-02 A. Papritz
  ## 2015-03-10 AP extended variogram parameter transformations
  ## 2015-04-07 AP changes for fitting anisotropic variograms
  
  list( 
    log = function(x) 1./x, 
    logit1 = function(x) 1. / (x - x^2),
    logit2 = function(x) 2. / (2.*x - x^2),
    logit3 = function(x) 2. / (4.*x - 3. - x^2),
    identity = function(x) rep(1., length(x)), ... 
  )  
  
}

## ======================================================================
bwd.transf <-
function( 
  ...
)
{
  
  ## definition of backward transformation of variogram parameters
  
  ## 2013-07-02 A. Papritz
  ## 2015-03-10 AP extended variogram parameter transformations
  ## 2015-04-07 AP changes for fitting anisotropic variograms
  ## 2016-08-03 AP corrections for logitx if argument is +/-Inf
  
  list( 
    log = function(x) exp(x), 
    logit1 = function(x){
      if( !is.finite(x) ){
        if(sign(x) < 0.) 0. else 1.
      } else exp(x) / (1. + exp(x))
    },
    logit2 = function(x){
      if( !is.finite(x) ){
        if(sign(x) < 0.) 0. else 2.
      } else 2. * exp(x) / (1. + exp(x))
    },
    logit3 = function(x){
      if( !is.finite(x) ){
        if(sign(x) < 0.) 1. else 3.
      } else (3. * exp(x) + 1.) / (1. + exp(x))
    },
    identity = function(x) x, ...
  )
}

## ======================================================================
control.rq <-
  function(
    ## specific arguments for rq: tau = 0.5, method = "br"
    tau = 0.5, rq.method = c("br", "fnb", "pfn"),
    ## specific arguments for rq.fit.br: tau = 0.5, alpha = 0.1, ci = FALSE, iid = TRUE, interp = TRUE, tcrit = TRUE
    rq.alpha = 0.1, ci = FALSE, iid = TRUE, interp = TRUE, tcrit = TRUE,
    ## specific arguments for rq.fit.fnb: tau = 0.5, rhs = (1 - tau) * apply(x, 2, sum), beta = 0.99995, eps = 1e-06
    rq.beta = 0.99995, eps = 1.e-06,
    ## specific arguments for rq.fit.pfn: tau = 0.5, Mm.factor = 0.8, max.bad.fixup = 3, eps = 1e-06
    Mm.factor = 0.8, max.bad.fixup = 3,
    ...
  )
{
  
  ## function sets meaningful defaults for selected arguments of function
  ## rq{quantreg}
  
  ## 2012-12-14 A. Papritz
  ## 2014-07-29 AP
  ## 2015-06-30 AP function and arguments renamed
  ## 2017-05-09 AP small changes
  
  rq.method = match.arg(rq.method)
  
  list( 
    tau = tau,  method = rq.method,
    alpha = rq.alpha, ci = ci, iid = iid, interp = interp, tcrit = tcrit,
    beta = rq.beta, eps = eps,
    Mm.factor = Mm.factor, max.bad.fixup = max.bad.fixup
  )
}


## ======================================================================
control.nleqslv <-
  function( 
    method = c( "Broyden", "Newton"), 
    global = c( "dbldog", "pwldog", "qline", "gline", "none" ),
    xscalm = c( "fixed", "auto" ),
    control = list( ftol = 1.e-4 ), 
    ...
  )
{
  
  ## function sets meaningful defaults for selected arguments of function
  ## nleqslv{nleqslv} 
  
  ## 2013-07-12 A. Papritz
  ## 2014-07-29 AP
  ## 2015-06-30 AP function and arguments renamed
  
  list( 
    method = match.arg( method ),
    global = match.arg( global ),
    xscalm = match.arg( xscalm ),
    control = control
  )
}

## ======================================================================
control.optim <-
  function( 
    method = c( "BFGS", "Nelder-Mead", "CG", "L-BFGS-B", "SANN", "Brent" ),
    lower = -Inf, upper = Inf,
    control = list(reltol = 1.e-5),
    ...
  )
{
  
  ## function sets meaningful defaults for selected arguments of function optim
  ## 2012-12-14 A. Papritz
  ## 2014-07-29 AP
  ## 2015-06-30 AP function and arguments renamed
  
  list( 
    method = match.arg( method ),
    lower = lower, upper = upper,
    control = control
  )
}

## ======================================================================
control.nlminb <-
  function( 
    control = list( rel.tol = 1.e-5 ),
    lower = -Inf, upper = Inf,
    ...
  )
{
  
  ## function sets meaningful defaults for selected arguments of function nlminb
  ## 2015-07-17 A. Papritz
  
  list( 
    control = control,
    lower = lower, upper = upper
  )
}

## ======================================================================
control.pcmp <-
  function( 
    pmm.ncores = 1, gcr.ncores = 1, max.ncores = detectCores(), 
    f = 1, sfstop = FALSE, allow.recursive = TRUE, 
    fork = !identical( .Platform[["OS.type"]], "windows" ),
    ... 
  )
{
  
  ## function sets meaningful defaults for parallelized computations
  
  ## 2014-07-29 AP
  ## 2015-06-30 AP function and arguments renamed
  ## 2015-07-29 AP changes for elimination of parallelized computation of gradient or estimating equations
  ## 2016-07-20 AP renamed function, separate ncores arguments various parallelized computations
  
  pmm.ncores <- min( pmm.ncores, max.ncores )
  gcr.ncores <- min( gcr.ncores, max.ncores )
  
  list( 
    pmm.ncores = pmm.ncores, gcr.ncores = gcr.ncores, max.ncores = max.ncores, 
    f = f, sfstop = sfstop, 
    allow.recursive = allow.recursive,
    fork = fork
  )

}

## ======================================================================

compress <- 
  function( m )
{
  
  ## function stores a list of or a single lower, upper triangular or
  ## symmetric matrix compactly

  ## 2013-06-12 AP substituting [["x"]] for $x in all lists

  aux <- function( x ){
    struc <- attr( x, "struc" )
    if( !is.null( struc ) ){
      switch( 
        struc,
        sym = {
          aux <- list( diag = diag( x ), tri = x[lower.tri(x)] )
          attr( aux, "struc" ) <- "sym"
          aux
        },
        lt = {
          aux <- list( diag = diag( x ), tri = x[lower.tri(x)] )
          attr( aux, "struc" ) <- "lt"
          aux
        }
        ,
        ut = {
          aux <- list( diag = diag( x ), tri = x[upper.tri(x)] )
          attr( aux, "struc" ) <- "ut"
          aux
        }
      )
    } else {
      x
    }
  }
  
  if( is.list( m ) ){
    lapply( m, aux )
  } else {
    aux ( m )
  }
  
  
  
}

## ======================================================================

expand <- 
  function( object )
{
  
  ## function expands a list of or a compactly stored lower, upper
  ## triangular or symmetric matrices

  ## 2013-06-12 AP substituting [["x"]] for $x in all lists

  aux <- function( x ){
    struc <- attr( x, "struc" )
    if( !is.null( struc ) ){
      switch( 
        struc,
        sym = {
          n <- length( x[["diag"]] )
          dn <- names( x[["diag"]] )
          aux <- matrix( 0., n, n )
          aux[lower.tri( aux )] <- x[["tri"]]
          aux <- aux + t( aux )
          diag( aux ) <- x[["diag"]]
          dimnames( aux ) <- list( dn, dn )
          attr( aux, "struc" ) <- "sym"
          aux
        },
        lt = {
          n <- length( x[["diag"]] )
          dn <- names( x[["diag"]] )
          aux <- matrix( 0., n, n )
          aux[lower.tri( aux )] <- x[["tri"]]
          diag( aux ) <- x[["diag"]]
          dimnames( aux ) <- list( dn, dn )
          attr( aux, "struc" ) <- "lt"
          aux
        }
        ,
        ut = {
          n <- length( x[["diag"]] )
          dn <- names( x[["diag"]] )
          aux <- matrix( 0., n, n )
          aux[upper.tri( aux )] <- x[["tri"]]
          diag( aux ) <- x[["diag"]]
          dimnames( aux ) <- list( dn, dn )
          attr( aux, "struc" ) <- "ut"
          aux
        }
      )
    } else {
      x
    }
  }
  
  ln <- names( object )
  if( is.list( object ) ){
    if( length( ln ) == 2L && all( ln == c( "diag", "tri" ) ) ){
      aux( object )
    } else {
      lapply( object, aux )
    }
  } else {
    object
  }  
  
}

## ======================================================================

param.names <- 
  function( model )
{
  
  ## function returns names of extra parameters of implemented variogram
  ## models (cf. Variogram{RandomFields})
  
  ## 2012-01-24 A. Papritz
  ## 2014-05-15 AP changes for version 3 of RandomFields
  
  switch(
    model,
    "RMaskey"         = "alpha",
    "RMbessel"        = "nu",
    "RMcauchy"        = "gamma",
    "RMcircular"      = NULL,
    "RMcubic"         = NULL,
    "RMdagum"         = c( "beta", "gamma" ),
    "RMdampedcos"     = "lambda",
    "RMdewijsian"     = "alpha",
    "RMexp"           = NULL,
    "RMfbm"           = "alpha",
    "RMgauss"         = NULL,
    "RMgencauchy"     = c( "alpha", "beta" ),
    "RMgenfbm"        = c( "alpha", "delta" ),
    "RMgengneiting"   = c( "kappa", "mu" ),
    "RMgneiting"      = NULL,
    "RMlgd"           = c( "alpha", "beta" ),
    "RMmatern"        = "nu",
    "RMpenta"         = NULL,
    "RMqexp"          = "alpha",
    "RMspheric"       = NULL,
    "RMstable"        = "alpha",
    "RMwave"          = NULL,
    "RMwhittle"       = "nu",
    stop( model, " variogram not implemented" )
  )
}

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

param.bounds <- 
function( model, d )
{
  
  ## function returns range of parameters for which variogram models are
  ## valid (cf.  Variogram{RandomFields})
  
  ## 2012-03-30 A. Papritz
  ## 2014-05-15 AP changes for version 3 of RandomFields
  
  switch(
    model,
    "RMaskey"         = list( alpha = c( 0.5 * (d + 1L), Inf ) ),
    "RMbessel"        = list( nu = c( 0.5 * (d - 2L), Inf ) ),
    "RMcauchy"        = list( gamma = c( 1.e-18, Inf ) ),
    "RMcircular"      = NULL,
    "RMcubic"         = NULL,
    "RMdagum"         = list( beta = c( 1.e-18, 1.), gamma = c( 1.e-18, 1.-1.e-18) ),
    "RMdampedcos"     = list( lambda = c( if( d > 2L ) sqrt(3.) else 1., Inf ) ),
    "RMdewijsian"     = list( alpha = c( 1.e-18, 2. ) ),
    "RMexp"           = NULL,
    "RMfbm"           = list( alpha = c( 1.e-18, 2.) ),
    "RMgauss"         = NULL,
    "RMgencauchy"     = list( alpha = c(1.e-18, 2.), beta = c(1.e-18, Inf) ),
    "RMgenfbm"        = list( alpha = c(1.e-18, 2.), delta = c(1.e-18, 1.-1.e-18) ),
    "RMgengneiting"   = list( kappa = c(1, 3), mu = c( d/2, Inf ) ),
    "RMgneiting"      = NULL,
    "RMlgd"           = list( 
                        alpha = c( 
                          1.e-18, 
                          if( d <= 3L ) 0.5 * (3L-d) else stop("dimension > 3 not allowed for RMlgd model" ) 
                        ), 
                        beta = c(1.e-18, Inf)
                      ),
    "RMmatern"        = list( nu = c(1.e-18, Inf) ),
    "RMpenta"         = NULL,
    "RMqexp"          = list( alpha = c(0., 1.) ),
    "RMspheric"       = NULL,
    "RMstable"        = list( alpha = c(1.e-18, 2.) ),
    "RMwave"          = NULL,
    "RMwhittle"       = list( nu = c(1.e-18, Inf) ),
    stop( model, " variogram not implemented" )
  )
}

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

default.fit.param <- 
function(
  variance = TRUE, snugget = FALSE, nugget = TRUE, scale = TRUE, 
  alpha = FALSE, beta = FALSE, delta = FALSE, gamma = FALSE, 
  kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE )
{
  
  ## function sets default flags for fitting variogram parameters
  
  ## 2015-11-27 A. Papritz
  
  c( 
	variance = variance, snugget = snugget, nugget = nugget, scale = scale,
	alpha = alpha, beta = beta, delta = delta, gamma = gamma,
	kappa = kappa, lambda = lambda, mu = mu, nu = nu
  )
  
}

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

default.fit.aniso <- 
function( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
{
  
  ## function sets default flags for fitting anisotropy parameters
  
  ## 2015-11-27 A. Papritz
  
  c( f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta )
  
}

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

default.aniso <- 
function(
  f1 = 1., f2 = 1., omega = 90., phi = 90., zeta = 0. )
{
  
  ## function sets default values for anisotropy parameters
  
  ## 2015-11-26 A. Papritz
  
  c( f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta )
  
}

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

profilelogLik <- function( object, values, use.fitted = TRUE, verbose = 0, 
  ncores = min( detectCores(), NROW(values) ) ){
  
  ## function to compute (restricted) likelihood profile for a georob fit
  
  ## arguments:
  
  ## object             a fitted georob object
  ## values             a matrix or dataframe with the fixed variogram parameter 
  ## use.fitted         control whether fitted variogram parameters should be used as initial values
  ##                    values for which remaining parameters are estimated
  ## ncores             number of cores for parallelized computations
  
  ## 2015-03-18 A. Papritz
  ## 2015-04-08 AP changes in returned results
  ## 2016-07-14 AP optimization
  ## 2016-07-20 AP changes for parallel computations
  ## 2016-07-28 AP returns gradient in results
  ## 2016-08-08 AP changes for nested variogra
  ## 2016-08-12 AP changes for nested variogram models
  ## 2017-12-22 AP improved memory management in parallel computations
  
  ## auxiliary function to fit model and return maximized (pseudo) log-likelihood
  
  f.aux <- function( i ){
    
    ## values, object, data are taken from parent environment
    
    ## set fixed initial values 
    
    values <- values[i, ]
    tmp <- strsplit( names(values), control.georob()[["sepstr"]], fixed = TRUE )
    cl <- object[["call"]]
    for( i in 1L:length(values) ){
      cl <- f.call.set_onexxx_to_value( cl, tmp[[i]][1L], values[i], as.integer(tmp[[i]][2L]) )
    }
    object[["call"]] <- cl
    
    fit <- update( object, data = data )
    
    ## extract fitted variogram parameters
    
    param.aniso <- unlist(lapply( 
        1L:length(fit[["variogram.object"]]),
        function( i, x ){
          
          x <- x[[i]]
          res <- c( x[["param"]][x[["fit.param"]]], x[["aniso"]][x[["fit.aniso"]]] )
          
          if( length(res) ){
            names(res) <- paste( names(res), i, sep = control.georob()[["sepstr"]] )
            res
          } else {
            NULL
          }
        }, x = fit[["variogram.object"]]
      ))
    
    ## results
    
    c( 
      loglik = logLik( 
        fit, warn = FALSE, REML = identical( object[["control"]][["ml.method"]], "REML" ) 
      ), 
      param.aniso,      
      coef( fit ),
      gradient = fit[["gradient"]],
      converged = fit[["converged"]]
    )
    
  }
  
  ## check whether all mandatory arguments have been provided
  
  if( missing(object) || missing(values) ) stop(
    "some mandatory arguments are missing" 
  )
  
  ## warning for robust fits
  
  if( object[["tuning.psi"]] < object[["control"]][["tuning.psi.nr"]] ){
    warning( 
      "likelihood approximated for robustly fitted model by likelihood of\n",
      "  equivalent Gaussian model with heteroscedastic nugget"
    )
  }
  
  if( !(is.matrix(values) || is.data.frame( values ) || is.list( values ) ) ) stop( 
    "'values' must be a dataframe or a matrix" 
  )
  
  ## get data.frame with required variables (note that the data.frame passed
  ## as data argument to georob must exist in GlobalEnv)
  
  data <- cbind(
    get_all_vars( 
      formula( object ), data = eval( getCall(object)[["data"]] )
    ),
    get_all_vars( 
      object[["locations.objects"]][["locations"]], eval( getCall(object)[["data"]] )
    )
  )
  
  if( identical( class( object[["na.action"]] ), "omit" ) ) data <- na.omit(data)
  
  ## select subset if appropriate
  
  if( !is.null( getCall(object)[["subset"]] ) ){
    data <- data[eval( getCall(object)[["subset"]] ), ]
  }
  
  ## extract variogram.object
  
  ## check values
  
  if( is.data.frame(values) ) values <- list( values )
  
  if( !identical( length(values), length(object[["variogram.object"]]) ) ) stop(
    "'values' must be a list of the same length as 'variogram.object'"  
  )
  
  ## check names of fixed variogram parameters and fix respective
  ## parameters in variogram.object
  
  values <- lapply( 
    1L:length(values),
    function( i, v, vo ){
      
      v <- v[[i]]
      vo <- vo[[i]]
      
      fixed.param.aniso <- colnames( v )
      
      ## match fixed.param.aniso
      fixed.param.aniso <- sapply( 
        fixed.param.aniso, match.arg, 
        choices = c( names(default.fit.param()), names(default.fit.aniso()) )
      )
      if( 
        any( !fixed.param.aniso %in% c( names( vo[["param"]] ), names( vo[["aniso"]] ) ) ) 
      ) stop( "column names of 'values' do not match names of variogram parameters" )
      
      colnames(v) <- fixed.param.aniso
      v
      
    }, v = values, vo = object[["variogram.object"]]
  )
  
  ## manipulate call so that fitted values are used as initial values
  
  cl <- f.call.set_allxxx_to_fitted_values( object )
  
  ## fix initial values for parameters present in values
  
  for( i in 1L:length(values) ){
    for( nme in colnames(values[[i]]) ){
      cl <- f.call.set_onefitxxx_to_value( cl, nme, FALSE, i )     
    }
  }
  
  ## update object call to avoid computation of covariance matrices and
  ## hessian and to set reparam = FALSE if variance parameters are fitted
  
  cl <- f.call.set_x_to_value( cl,  "verbose", verbose )
  
  reparam <- !any( 
    unique( sapply( values, colnames ) )  %in% c( "variance", "snugget", "nugget" )
  )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "reparam", reparam )
  
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "hessian", FALSE )
  
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.bhat", FALSE )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.betahat", FALSE )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.delta.bhat", FALSE )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.delta.bhat.betahat", FALSE )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.ehat", FALSE )
  cl <- f.call.set_x_to_value_in_fun( cl, "control", "control.georob", "cov.ehat.p.bhat", FALSE )
  
  object[["call"]] <- cl
  
  ## rename variogram parameters in values
  
  values <- lapply( 
    1L:length(values),
    function( i, v, vo ){
      v <- v[[i]]
      names( v ) <- paste( names(v), i, sep = control.georob()[["sepstr"]] )
      v
    }, v = values
  )
  
  ## expand values to data frame
  
  if( length(values) > 1L ){
    tmp <- values
    values <- as.list( tmp[[1L]] )
    
    for( i in 2L:length(tmp) ){
      for( j in colnames( tmp[[i]] ) ){
        values <- c( values, as.list( tmp[[i]][, j, drop=FALSE] ) )
      }
    }
    values <- expand.grid( values )
  } else {
    values <- values[[1L]]
  }
  
  ## loop over all elements of values
  
  values <- as.matrix( values )
  
  ## set default value for control of forking if missing (required for backward compatibility)
  
  if( is.null( object[["control"]][["pcmp"]][["fork"]] ) ){
    object[["control"]][["pcmp"]][["fork"]] <- !identical( .Platform[["OS.type"]], "windows" )
  }
  
  if( ncores > 1L && !object[["control"]][["pcmp"]][["fork"]] ){
    
    ## create a 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( "values", "object", "data" ), envir = environment() )
    
    result <- parLapply(
      clstr, 
      1L:NROW(values),
      f.aux
    )
    
    f.stop.cluster( clstr )
    
  } else {
    
    ## fork child processes on non-windows OS
    
    result <- mclapply(
      1L:NROW(values),
      f.aux, 
      mc.cores = ncores,
      mc.allow.recursive = object[["control"]][["pcmp"]][["allow.recursive"]]
    )
    
  }
  
  ## collect results
  
  result <- as.data.frame( cbind( values, t( simplify2array( result ) ) ) )
  
  if( identical( length( object[["variogram.object"]]), 1L ) ){
    tmp <- colnames(result)
    tmp <- gsub( 
      paste( control.georob()[["sepstr"]], "1", sep = "" ),
      "", tmp
    )
    colnames(result) <- tmp    
  }
  result
  
}

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.