R/likGRF.R

##
## Functions for likelihood based inference in the geoR package
## ------------------------------------------------------------
##
## Includes:
##    - the main function for ML estimation - likfit()
##    - methods for class likGRF
##    - some other functions

".check.geoRparameters.values" <- function(list, messages = TRUE)
{
  if(!is.null(list$nugget))
    if(list$nugget < 0) stop("value for nugget must be non-negative")
  if(!is.null(list$psiA))
    if(list$psiA < 0) stop("value for psiA must be non-negative")
  if(!is.null(list$psiR))
    if(list$psiR < 1) stop("value for psiA must be >= 1")
  if(!is.null(list$kappa)){
    if(is.null(list$cov.model)) stop("cov.model must be provided when checking values of kappa")
    if(list$kappa[1] <= 0) stop("parameter kappa[1] must be greater than 0")
    if(any(list$cov.model == c("powered.exponential", "matern", "gneiting.matern", "gencauchy","cauchy"))){
      if(any(list$cov.model == c("gneiting.matern", "gencauchy"))){
        if(length(list$kappa) != 2)
          stop(paste("kappa must be of length 2 for the", list$cov.model, "correlation function"))
        if(list$cov.model == "gencauchy" && (list$kappa[2] <=0 | list$kappa[2] >2))
          stop("for the gencauchy model the kappa must be within (0,2]")          
      }
      else{
        if(list$cov.model == "powered.exponential" && (list$kappa <=0 | list$kappa >2))
          stop("for the powered.exponential model the kappa must be within (0,2]")          
        if(list$cov.model == "stable" && (list$kappa <=0 | list$kappa >2))
          stop("for the stable model the kappa must be within (0,2]")          
      }
    }
    else{
      if(messages && !is.null(list$kappa))
        cat(paste("kappa not used for the",list$cov.model, "correlation function\n"))
    }
  }
  return(invisible())
}

"likfit" <-
  function (geodata, coords=geodata$coords, data=geodata$data,
            trend = "cte", ini.cov.pars,
            fix.nugget = FALSE, nugget = 0, 
            fix.kappa = TRUE, kappa = 0.5, 
            fix.lambda = TRUE, lambda = 1, 
            fix.psiA = TRUE, psiA = 0, fix.psiR = TRUE, psiR = 1, 
            cov.model, realisations, lik.method = "ML",
            components = TRUE, nospatial = TRUE, limits = pars.limits(), 
            print.pars = FALSE, messages, ...) 
{
  name.geodata <- deparse(substitute(geodata))
  ##
  ## Checking input
  ##
  call.fc <- match.call()
  ldots <- list(...)
  temp.list <- list()
  temp.list$print.pars <- print.pars
  if(missing(messages))
    messages.screen <- as.logical(ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages")))
  else messages.screen <- messages
  ##
  if(!missing(ini.cov.pars)){
    if(any(class(ini.cov.pars) == "eyefit")){
ini.cov.pars <- ini.cov.pars[[1]]
#      cov.model <- ini.cov.pars[[1]]$cov.model
#      kappa <- ini.cov.pars[[1]]$kappa
    }
    if(any(class(ini.cov.pars) == "variomodel")){
      cov.model <- ini.cov.pars$cov.model
      kappa <- ini.cov.pars$kappa
}
  }
  if(missing(cov.model)) cov.model <- "matern"
  cov.model <- match.arg(cov.model, choices = .geoR.cov.models)
  if(cov.model == "stable") cov.model <- "powered.exponential"
  if(any(cov.model == c("power", "gneiting.matern", "gencauchy")))
     stop(paste("parameter estimation for", cov.model, "is not yet implemented"))
  ##  if(any(cov.model == c("gneiting.matern", "gencauchy"))){
  ##    if(length(kappa != 2))
  ##      stop(paste(cov.model, "requires two values in the argument kappa"))
  ##    if(length(fix.kappa) == 1) fix.kappa <- rep(fix.kappa, 2) 
  ##    stop("parameter estimation for gneiting.matern model is not yet implemented")
  ##  }
  fixed.pars <- list(cov.model=cov.model)
  if(fix.nugget) fixed.pars$nugget <- nugget
  if(fix.kappa) fixed.pars$kappa <- kappa
  if(fix.psiA) fixed.pars$psiA <- psiA
  if(fix.psiR) fixed.pars$psiR <- psiR
  .check.geoRparameters.values(list=fixed.pars, messages = messages.screen)
  if(cov.model == "matern" & all(kappa == 0.5)) cov.model <- "exponential"
  temp.list$cov.model <- cov.model
  if(cov.model == "powered.exponential")
    if(limits$kappa["upper"] > 2) limits$kappa["upper"] <- 2
  if(cov.model == "gencauchy")
    if(limits$kappa2["upper"] > 2) limits$kappa2["upper"] <- 2
  ##
  ## Likelihood method
  ##
#####
##### temporary code back compatibility for argument "method"
  lik.MET <- c("ML", "ml", "RML", "REML", "rml", "reml")
  MET <- pmatch(names(ldots), "method") == 1
  if(!is.na(MET) && any(MET) && (ldots[[which(MET)]] %in% lik.MET)){
    warning("argument \"method\" has changed and is now used as an argument to be passed to optim(). Use \"lik.method\" to define the likelihood method")
    lik.method <- lik.MET[pmatch(ldots[[which(MET)]], lik.MET)]
    ldots[which(as.logical(pmatch(names(ldots), "method", nomatch=0)))] <- NULL
  }
#####
  method.lik <- lik.method
  if(method.lik %in% c("REML","reml","rml","RML"))  method.lik <- "RML"
  if(method.lik %in% c("ML", "ml")) method.lik <- "ML"
  if(method.lik == "ML" & cov.model == "power")
    stop("\n\"power\" model can only be used with method.lik=\"RML\".\nBe sure that what you want is not \"powered.exponential\"")
  temp.list$method.lik <- method.lik
  ##
  ## setting coordinates, data and covariate matrices
  ##
  coords <- as.matrix(coords)
  data <- as.vector(data)
  n <- length(data)
  if((nrow(coords) != n) | (2*n) != length(coords))
    stop("\nnumber of locations does not match with number of data")
  if(missing(geodata))
    xmat <- trend.spatial(trend=trend, geodata=list(coords = coords, data = data))
  else xmat <- unclass(trend.spatial(trend=trend, geodata=geodata))
  xmat.contrasts  <- attr(xmat,"contrasts")
  xmat <- unclass(xmat)
  if(nrow(xmat) != n)
    stop("trend matrix has dimension incompatible with the data")
  .solve.geoR(crossprod(xmat))
  beta.size <- temp.list$beta.size <- dim(xmat)[2]
  ##
  ## setting a factor to indicate different realisations
  ##
  if(missing(realisations))
    realisations <- as.factor(rep(1, n))
  else{
    if(!missing(geodata)){
        real.name <- deparse(substitute(realisations))
        if(all(isTRUE(as.logical(real.name))))
          if(is.null(geodata$realisations)) stop("element realisation not available in the geodata object")
          else realisations <- geodata$realisations
      else{
        if(!is.null(geodata[[real.name]]))
          realisations <- geodata[[real.name]]
      }
    }
    if(length(realisations) != n)
      stop("realisations must be a vector with the same length of the data")
    realisations <- as.factor(realisations)
  }
  temp.list$realisations <- realisations
  nrep <- temp.list$nrep <- length(levels(realisations))
  ind.rep <- split(1:n, realisations)
  vecdist <- function(x){as.vector(dist(x))}
  ##
  ## Initial values for parameters
  ##
  ## have to consider transformation, residuals from trend etc
#  var.data <- mean(tapply(data, realisations, var))
#  d.max <- max(by(ap$coords, ap$realisations, function(x) max(dist(x))))
#  if(missing(ini.cov.pars))
#    ini.cov.pars <- expand.grid(var.data/2, 3*var.data/4, var.data)
  if(any(class(ini.cov.pars) == "eyefit")){
    init <- nugget <- kappa <- NULL
    for(i in 1:length(ini.cov.pars)){
      init <- drop(rbind(init, ini.cov.pars[[i]]$cov.pars))
      nugget <- c(nugget, ini.cov.pars[[i]]$nugget)
      if(cov.model == "gneiting.matern")
        kappa <- drop(rbind(kappa, ini.cov.pars[[i]]$kappa))
      else
        kappa <- c(kappa, ini.cov.pars[[i]]$kappa)
    }
    ini.cov.pars <- init
  }
  if(any(class(ini.cov.pars) == "variomodel")){
    nugget <- ini.cov.pars$nugget
    kappa <- ini.cov.pars$kappa
    ini.cov.pars <- ini.cov.pars$cov.pars
  }
  if(is.matrix(ini.cov.pars) | is.data.frame(ini.cov.pars)){
    ini.cov.pars <- as.matrix(ini.cov.pars)
    if(nrow(ini.cov.pars) == 1)
      ini.cov.pars <- as.vector(ini.cov.pars)
    else{
      if((cov.model != "pure.nugget") & (ncol(ini.cov.pars) != 2))
        stop("\nini.cov.pars must be a matrix or data.frame with 2 components: \ninitial values for sigmasq and phi")
    }
  }
  if(is.vector(ini.cov.pars)){
    if((cov.model != "pure.nugget") & (length(ini.cov.pars) != 2))
      stop("\nini.cov.pars must be a vector with 2 components: \ninitial values for sigmasq and phi")
  }
  ##
  ## Checking for multiple initial values for preliminar search of   
  ## best initial value
  ##
  if(is.matrix(ini.cov.pars) | (length(nugget) > 1) | (length(kappa) > 1) | (length(lambda) > 1) | (length(psiR) > 1) | (length(psiA) > 1)){
    if(messages.screen) cat("likfit: searching for best initial value ...")
    ini.temp <- matrix(ini.cov.pars, ncol=2)
    grid.ini <- as.matrix(expand.grid(sigmasq=unique(ini.temp[,1]), phi=unique(ini.temp[,2]), tausq=unique(nugget),
                                      kappa=unique(kappa), lambda=unique(lambda), psiR=unique(psiR), psiA=unique(psiA)))
    assign(".likGRF.dists.vec",  lapply(split(as.data.frame(coords), realisations), vecdist), envir=.geoR.env)
    temp.f <- function(parms, coords, data, temp.list)
      return(loglik.GRF(geodata = geodata,
                        coords = coords, data = as.vector(data),
                        cov.model=temp.list$cov.model,
                        cov.pars=parms[1:2],
                        nugget=parms["tausq"], kappa=parms["kappa"],
                        lambda=parms["lambda"], psiR=parms["psiR"],
                        psiA=parms["psiA"], trend= trend,
                        method.lik=temp.list$method.lik,
                        compute.dists=FALSE,
                        realisations = realisations))
    grid.lik <- apply(grid.ini, 1, temp.f, coords = coords,
                      data = data, temp.list = temp.list)
    grid.ini <- grid.ini[(grid.lik != Inf) & (grid.lik != -Inf) & !is.na(grid.lik) & !is.nan(grid.lik),, drop=FALSE] 
    grid.lik <- grid.lik[(grid.lik != Inf) & (grid.lik != -Inf) & !is.na(grid.lik) & !is.nan(grid.lik)] 
    ini.temp <- grid.ini[which(grid.lik == max(grid.lik)),, drop=FALSE]
    if(all(ini.temp[,"phi"] == 0)) ini.temp <- ini.temp[1,, drop=FALSE]
    rownames(ini.temp) <- "initial.value"
    if(messages.screen){
      cat(" selected values:\n")
      print(rbind(format(ini.temp, digits=2), status=ifelse(c(FALSE, FALSE, fix.nugget, fix.kappa, fix.lambda, fix.psiR, fix.psiA), "fix", "est")))
      cat(paste("likelihood value:", max(grid.lik), "\n"))
    }
    dimnames(ini.temp) <- NULL
    ini.cov.pars <- ini.temp[1:2]
    nugget <- ini.temp[3]
    kappa <- ini.temp[4]
    lambda <- ini.temp[5]
    psiR <- ini.temp[6]
    psiA <- ini.temp[7]
    grid.ini <- NULL
    remove(".likGRF.dists.vec", envir=.geoR.env)
  }
  ##
  tausq <- nugget
  ##
  ## Box-Cox transformation for fixed lambda
  ##
  if(fix.lambda) {
    if(abs(lambda - 1) < 0.0001) {
      temp.list$log.jacobian <- 0
      temp.list$z <- as.vector(data)
    }
    else {
      if(any(data <= 0))
        stop("Transformation option not allowed when there are zeros or negative data")
      Jdata <- data^(lambda - 1)
      if(any(Jdata <= 0))
        temp.list$log.jacobian <- log(prod(Jdata))
      else temp.list$log.jacobian <- sum(log(Jdata))
      Jdata <- NULL
      if(abs(lambda) < 0.0001)
        temp.list$z <- log(data)
      else temp.list$z <- ((data^lambda) - 1)/lambda
    }
  }
  else{
    temp.list$z <- as.vector(data)
    temp.list$log.jacobian <- NULL
  }
  ##
  ## Coordinates transformation for fixed anisotropy parameters
  ##
  if(fix.psiR & fix.psiA){
    if(psiR != 1 | psiA != 0)
      coords <- coords.aniso(coords, aniso.pars=c(psiA, psiR))
      assign(".likGRF.dists.vec", lapply(split(as.data.frame(coords), realisations), vecdist), envir=.geoR.env)
    range.dist <- range(get(".likGRF.dists.vec", envir=.geoR.env))
    max.dist <- max(range.dist)
    min.dist <- min(range.dist)
  }
  ##
  ##
  ##
  ini <- ini.cov.pars[2]
  ##  fixed.pars <- NULL
  lower.optim <- c(limits$phi["lower"])
  upper.optim <- c(limits$phi["upper"])
  fixed.values <- list()
  if(fix.nugget) {
    ##    fixed.pars <- c(fixed.pars, 0)
    fixed.values$tausq <- nugget
  }
  else {
    ini <- c(ini, nugget/ini.cov.pars[1])
    lower.optim <- c(lower.optim, limits$tausq.rel["lower"])
    upper.optim <- c(upper.optim, limits$tausq.rel["upper"])
  }
  if(fix.kappa){
    ##    fixed.kappa <- c(fixed.pars, kappa)
    fixed.values$kappa <- kappa
  }
  else {
    ini <- c(ini, kappa)
    lower.optim <- c(lower.optim, limits$kappa["lower"])
    upper.optim <- c(upper.optim, limits$kappa["upper"])
  }
  if(fix.lambda){
    ##    fixed.pars <- c(fixed.pars, lambda)
    fixed.values$lambda <- lambda
  }
  else {
    ini <- c(ini, lambda)
    lower.optim <- c(lower.optim, limits$lambda["lower"])
    upper.optim <- c(upper.optim, limits$lambda["upper"])
  }
  if(fix.psiR){
    ##    fixed.pars <- c(fixed.pars, psiR)
    fixed.values$psiR <- psiR
  }
  else {
    ini <- c(ini, psiR)
    lower.optim <- c(lower.optim, limits$psiR["lower"])
    upper.optim <- c(upper.optim, limits$psiR["upper"])
  }
  if(fix.psiA){
    ##    fixed.pars <- c(fixed.pars, psiA)
    fixed.values$psiA <- psiA
  }
  else {
    ini <- c(ini, psiA)
    lower.optim <- c(lower.optim, limits$psiA["lower"])
    upper.optim <- c(upper.optim, limits$psiA["upper"])
  }
  ## This must be here, after the previous ones:
  if(fix.nugget & nugget > 0){
    ## Warning: Inverting order here, ini will be now: c(phi,sigmasg)
    ini <- c(ini, ini.cov.pars[1])
    lower.optim <- c(lower.optim, limits$sigmasq["lower"])
    upper.optim <- c(upper.optim, limits$sigmasq["upper"])
    ##    fixed.pars <- c(fixed.pars, ini.cov.pars[1])
    ##    fixed.values$sigmasq <- 0
  }
  ##
  names(ini) <- NULL
  if(length(ini) == 1) justone <- TRUE
  else justone <- FALSE
  ##
  ip <- list(f.tausq = fix.nugget, f.kappa = fix.kappa,
             f.lambda = fix.lambda,
             f.psiR = fix.psiR, f.psiA = fix.psiA)
  ##
  npars <- beta.size + 2 + sum(unlist(ip)==FALSE)
  temp.list$coords <- coords
  temp.list$xmat <- split(as.data.frame(unclass(xmat)), realisations)
  temp.list$xmat <- lapply(temp.list$xmat, as.matrix)
  temp.list$n <- as.vector(unlist(lapply(temp.list$xmat, nrow)))
  ##
  ## Constant term in the likelihood
  ##
  temp.list$loglik.cte <- rep(0, nrep)
  for(i in 1:nrep){
    if(method.lik == "ML"){
      if(ip$f.tausq & (tausq > 0))
        temp.list$loglik.cte[i] <-  (temp.list$n[i]/2)*(-log(2*pi))
      else
        temp.list$loglik.cte[i] <-  (temp.list$n[i]/2)*(-log(2*pi) +
                                                        log(temp.list$n[i]) -1)
    }
    if(method.lik == "RML"){
      xx.eigen <- eigen(crossprod(temp.list$xmat[[i]]),
                        symmetric = TRUE, only.values = TRUE)
      if(ip$f.tausq & (tausq > 0))
        temp.list$loglik.cte[i] <- - ((temp.list$n[i]-beta.size)/2)*(log(2*pi)) +
          0.5 * sum(log(xx.eigen$values))
      else
        temp.list$loglik.cte[i] <-  - ((temp.list$n[i]-beta.size)/2)*(log(2*pi)) +
          ((temp.list$n[i]-beta.size)/2)*(log(temp.list$n[i]-beta.size)) -
            ((temp.list$n[i]-beta.size)/2) + 0.5 * sum(log(xx.eigen$values))
    }
  }
  ##  
  if(messages.screen) {
    cat("---------------------------------------------------------------\n")
    cat("likfit: likelihood maximisation using the function ")
    if(is.R()){if(justone) cat("optimize.\n") else cat("optim.\n")} else cat("nlminb.\n")
    cat("likfit: Use control() to pass additional\n         arguments for the maximisation function.")
    cat("\n        For further details see documentation for ")
    if(is.R()){if(justone) cat("optimize.\n") else cat("optim.\n")} else cat("nlminb.\n")        
    cat("likfit: It is highly advisable to run this function several\n        times with different initial values for the parameters.\n")
    cat("likfit: WARNING: This step can be time demanding!\n")
    cat("---------------------------------------------------------------\n")
  }
  ##
  ## Numerical minimization of the -loglikelihood
  ##
  if(length(ini) == 1){
    if(upper.optim == Inf) upper.optim <- 50*max.dist
    lik.minim <- do.call("optimize", c(list(.negloglik.GRF,
                                            lower=lower.optim,
                                            upper=upper.optim,
                                            fp=fixed.values,
                                            ip=ip, temp.list = temp.list), ldots))
    lik.minim <- list(par = lik.minim$minimum,
                      value = lik.minim$objective,
                      convergence = 0,
                      message = "function optimize used")      
  }
  else{
    MET <- pmatch(names(ldots), names(formals(optim)))
    if(is.na(MET) || all(names(formals(optim))[MET] != "method"))
      ldots$method <- "L-BFGS-B"
    if(!is.null(names(ldots))){
      names(ldots)[which(as.logical(pmatch(names(ldots), "method", nomatch=0)))] <- "method"
    }
    if(!is.null(ldots$method) && ldots$method == "L-BFGS-B"){
      ldots$lower <- lower.optim
      ldots$upper <- upper.optim
    }
    lik.minim <- do.call("optim", c(list(par = ini, fn = .negloglik.GRF,
                                         fp=fixed.values, ip=ip, temp.list = temp.list), ldots))
    ##      lik.minim <- optim(par = ini, fn = .negloglik.GRF, method=optim.METHOD
    ##                         lower=lower.optim, upper=upper.optim,
    ##                         fp=fixed.values, ip=ip, temp.list = temp.list, ...)
  }
  ##
  if(messages.screen) cat("likfit: end of numerical maximisation.\n")
  par.est <- lik.minim$par
  if(any(par.est < 0)) par.est <- round(par.est, digits=12)
  phi <- par.est[1]
  ##
  ## Values of the maximised likelihood
  ##
  if(is.R())
    loglik.max <- - lik.minim$value
  else
    loglik.max <- - lik.minim$objective
  ##
  ## Assigning values for estimated parameters
  ##
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    psiA <- par.est[2]
  }
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    psiR <- par.est[2]
  }
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    psiR <- par.est[2]
    psiA <- par.est[3]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    lambda  <- par.est[2]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    lambda  <- par.est[2]
    psiA <- par.est[3]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    lambda  <- par.est[2]
    psiR <- par.est[3]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    lambda  <- par.est[2]
    psiR <- par.est[3]
    psiA <- par.est[4]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    kappa  <-  par.est[2]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    kappa  <-  par.est[2]
    psiA <- par.est[3]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    kappa  <-  par.est[2]
    psiR <- par.est[3]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    kappa  <-  par.est[2]
    psiR <- par.est[3]
    psiA <- par.est[4]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    kappa <-  par.est[2]
    lambda <- par.est[3]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    kappa <-  par.est[2]
    lambda <- par.est[3]
    psiA <- par.est[4]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    kappa <-  par.est[2]
    lambda <- par.est[3]
    psiR<- par.est[4]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    kappa <-  par.est[2]
    lambda <- par.est[3]
    psiR<- par.est[4]
    psiA<- par.est[5]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    psiA<- par.est[3]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    psiR<- par.est[3]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    psiR<- par.est[3]
    psiA<- par.est[4]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    lambda <- par.est[3]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    lambda <- par.est[3]
    psiA <- par.est[4]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    lambda <- par.est[3]
    psiR <- par.est[4]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    lambda <- par.est[3]
    psiR <- par.est[4]
    psiA <- par.est[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    psiA <- par.est[4]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    psiR <- par.est[4]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    psiR <- par.est[4]
    psiA <- par.est[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    lambda <- par.est[4]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    lambda <- par.est[4]
    psiA <- par.est[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    lambda <- par.est[4]
    psiR <- par.est[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- par.est[2]
    kappa <-  par.est[3]
    lambda <- par.est[4]
    psiR <- par.est[5]
    psiA <- par.est[6]
  }
  ##
  if(fix.nugget & nugget > 0){
    sigmasq <- par.est[length(par.est)]
    if(sigmasq > 1e-12) tausq <- nugget/sigmasq
    check.sigmasq <- TRUE
  }
  else check.sigmasq <- FALSE
  ##
  ##
  ## Transforming data according to the estimated lambda (Box-Cox) parameter
  ##
  if(!fix.lambda) {
    if(abs(lambda - 1) < 0.0001) {
      log.jacobian.max <- 0
    }
    else {
      if(any(data^(lambda - 1) <= 0))
        log.jacobian.max <- log(prod(data^(lambda - 1)))
      else log.jacobian.max <- sum(log(data^(lambda - 1)))
      temp.list$z <- ((data^lambda)-1)/lambda
    }
  }
  else{
    log.jacobian.max <- temp.list$log.jacobian
  }
  data.rep <- split(temp.list$z, realisations)
  coords.rep <- split(as.data.frame(coords), realisations)
  coords.rep <- lapply(coords.rep, as.matrix)
  ##
  ## Transforming coords for estimated anisotropy (if the case)
  ##
  if(fix.psiR & fix.psiA)
    remove(".likGRF.dists.vec", envir=.geoR.env)
  else{
    if(round(psiR, digits=6) != 1 | round(psiA, digits=6) != 0)
      coords <- coords.aniso(coords, aniso.pars=c(psiA, psiR))
    rangevecdist <- function(x){range(as.vector(dist(x)))}
    range.dist <- lapply(split(as.data.frame(coords), realisations), rangevecdist)
    range.dist <- range(as.vector(unlist(range.dist)))
    max.dist <- max(range.dist)
    min.dist <- min(range.dist)
  }      
#  gc(verbose=FALSE)
  ##
  ## Computing estimated beta and tausq/sigmasq (if the case)
  ##
  xivx <- matrix(0, ncol=beta.size, nrow=beta.size)
  xivy <- matrix(0, ncol=1, nrow=beta.size)
  yivy <- 0
  for(i in 1:nrep){
    ni <- temp.list$n[i]
    if((phi < 1e-12))
      V <- diag(x=(1+tausq), ni)
    else{
      if(check.sigmasq){
        if(sigmasq < 1e-12){
          if(!fix.nugget)
            V <- diag(x=(1+tausq), ni)
          else
            V <- diag(x=sqrt(tausq), ni)          
        }
        else
          V <- varcov.spatial(coords = coords.rep[[i]],
                              cov.model = cov.model,
                              kappa = kappa, nugget = tausq,
                              cov.pars = c(1, phi))$varcov
      }
      else
        V <- varcov.spatial(coords = coords.rep[[i]],
                            cov.model = cov.model,
                            kappa = kappa, nugget = tausq,
                            cov.pars = c(1, phi))$varcov
    }
    ivyx <- solve(V,cbind(data.rep[[i]],temp.list$xmat[[i]]))
    xivx <- xivx + crossprod(ivyx[,-1],temp.list$xmat[[i]])
    xivy <- xivy + crossprod(ivyx[,-1],data.rep[[i]])
    yivy <- yivy + crossprod(data.rep[[i]],ivyx[,1])
  }
  betahat <- .solve.geoR(xivx, xivy)
  res <- as.vector(temp.list$z - xmat %*% betahat)
  if(!fix.nugget | (nugget < 1e-12)){
    ssres <- as.vector(yivy - 2*crossprod(betahat,xivy) +
                       crossprod(betahat,xivx) %*% betahat)  
    if(method.lik == "ML")
      sigmasq <- ssres/n
    else
      sigmasq <- ssres/(n - beta.size)
  }
  if(fix.nugget){
    if(nugget > 0)
      tausq <- nugget
  }
  else tausq <- tausq * sigmasq
  betahat.var <- .solve.geoR(xivx)
  if(sigmasq > 1e-12) betahat.var <- sigmasq * betahat.var
#  if(!fix.nugget & phi < 1e-16){
#    tausq <- sigmasq + tausq
#    sigmasq <- 0
#  }
  ##
  ## Preparing output
  ##
  if((phi < 0.001*min.dist)){
    tausq <- tausq + sigmasq
    sigmasq <- 0
  }
  if((sigmasq < 1e-12)) phi <- 0
  ##
  n.model.pars <- beta.size + 7
  par.su <- data.frame(status=rep(-9,n.model.pars))
  ind.par.su <- c(rep(0, beta.size), ip$f.tausq, 0, 0, ip$f.kappa,
                  ip$f.psiR, ip$f.psiA,ip$f.lambda)
  par.su$status <- ifelse(ind.par.su,"fixed", "estimated")
  par.su$values <- round(c(betahat, tausq, sigmasq, phi, kappa, psiR, psiA, lambda), digits=4)
  if(beta.size == 1) beta.name <- "beta"
  else beta.name <- paste("beta", 0:(beta.size-1), sep="")
  row.names(par.su) <- c(beta.name, "tausq", "sigmasq", "phi", "kappa",
                             "psiR", "psiA", "lambda")
  par.su <- par.su[c((1:(n.model.pars-3)), n.model.pars-1, n.model.pars-2, n.model.pars),] 
  ##
  lik.results <- list(cov.model = cov.model,
                      nugget = tausq,
                      cov.pars=c(sigmasq, phi),
                      sigmasq = sigmasq,
                      phi = phi,
                      kappa = kappa,
                      beta = as.vector(betahat),
                      beta.var = betahat.var,
                      lambda = lambda,
                      aniso.pars = c(psiA = psiA, psiR = psiR),
                      tausq = tausq,
                      practicalRange = practicalRange(cov.model=cov.model,
                        phi = phi, kappa = kappa),
                      method.lik = method.lik, trend = trend,
                      loglik = loglik.max,
                      npars = npars,
                      AIC = -2 * (loglik.max - npars),
                      BIC = -2 * (loglik.max - 0.5 * log(n) * npars),
#                      residuals = res,
                      parameters.summary = par.su,
                      info.minimisation.function = lik.minim,
                      max.dist = max.dist,
                      trend = trend,
                      trend.matrix= xmat,
                      transform.info = list(fix.lambda = fix.lambda,
                        log.jacobian = log.jacobian.max))
  ##
  ## Likelihood results for the model without spatial correlation
  ##
  if(nospatial){
    if(fix.lambda){
      beta.ns <- .solve.geoR(crossprod(xmat), crossprod(xmat, temp.list$z))
      ss.ns <- sum((as.vector(temp.list$z - xmat %*% beta.ns))^2)
      if(method.lik == "ML"){
        nugget.ns <- ss.ns/n
        loglik.ns <- (n/2)*((-log(2*pi)) - log(nugget.ns) - 1) + temp.list$log.jacobian
      }
      if(method.lik == "RML"){
        nugget.ns <- ss.ns/(n-beta.size)
        loglik.ns <- ((n-beta.size)/2)*((-log(2*pi)) - log(nugget.ns) -1) +
          temp.list$log.jacobian
      }
      npars.ns <- beta.size + 1 + !fix.lambda
      lambda.ns <- lambda
    }
    else{
      if(is.R())
        lik.lambda.ns <- optim(par=1, fn = .negloglik.boxcox,
                               method = "L-BFGS-B",
                               lower = limits$lambda["lower"],
                               upper = limits$lambda["upper"],
                               data = data, xmat = xmat,
                               lik.method = method.lik)
      else
        lik.lambda.ns <- nlminb(par=1, fn = .negloglik.boxcox,
                                lower=limits$lambda["lower"],
                                upper=limits$lambda["upper"],
                                data = data, xmat = xmat,
                                lik.method = method.lik)
      lambda.ns <- lik.lambda.ns$par
      if(abs(lambda) < 0.0001) tdata.ns <- log(data)
      else tdata.ns <- ((data^lambda.ns)-1)/lambda.ns
      beta.ns <- .solve.geoR(crossprod(xmat),crossprod(xmat,tdata.ns))
      ss.ns <- sum((as.vector(tdata.ns - xmat %*% beta.ns))^2)
      if(is.R())
        value.min.ns <- lik.lambda.ns$value
      else
        value.min.ns <- lik.lambda.ns$objective
      if(method.lik == "ML"){
        loglik.ns <- (- value.min.ns)+ (n/2)*((-log(2*pi)) + log(n) - 1)
        nugget.ns <- ss.ns/n
      }
      if(method.lik == "RML"){
        nugget.ns <- ss.ns/(n-beta.size)
        loglik.ns <- (- value.min.ns)+ ((n-beta.size)/2)*((-log(2*pi)) +
                                                          log(n-beta.size) - 1)
      }      
      npars.ns <- beta.size + 1 + !fix.lambda
    }
    lik.results$nospatial <- list(beta.ns = beta.ns, variance.ns = nugget.ns,
                                  loglik.ns = loglik.ns, npars.ns = npars.ns,
                                  lambda.ns = lambda.ns, AIC.ns = -2 * (loglik.ns - npars.ns),
                                  BIC.ns = -2 * (loglik.ns - 0.5 * log(n) * npars.ns))
  }
  ##
  ## Assigning names to the components of the mean vector beta
  ##
  if(length(lik.results$beta.var) == 1)
    lik.results$beta.var <- as.vector(lik.results$beta.var)
  if(length(lik.results$beta) > 1){
    ##    if(inherits(trend, "formula") || (!is.null(class(trend)) && any(class(trend) == "trend.spatial")))
    if(inherits(trend, "formula") || (length(class(trend)) > 0 && any(class(trend) == "trend.spatial")))
      beta.names <- c("intercept", paste("covar", 1:(ncol(xmat)-1), sep = ""))
    else
      if(trend == "1st")
        beta.names <- c("intercept", "x", "y")
      else
        if(trend == "2nd")
          beta.names <- c("intercept", "x", "y", "x2", "xy", "y2")
    names(lik.results$beta) <- beta.names
  }
  ##
  ## Computing residuals and predicted values
  ## (isolated components of the model)
  ##
  if(components) {
    if(!fix.psiR & !fix.psiA)
      if(psiR != 1 | psiA != 0)
        coords <- coords.aniso(coords, aniso.pars=c(psiA, psiR))
    #coords.rep <- split(as.data.frame(coords), realisations)
    #res.rep <- split(res, realisations)
    trend.comp <- temp.list$z - res
    spatial.comp <- list()
    for(i in 1:nrep){
#      invcov <- varcov.spatial(coords = coords[ind.rep[[i]],], cov.model = cov.model, 
#                               kappa = kappa, nugget = tausq,
#                               cov.pars = c(sigmasq, phi), inv=TRUE)$inverse 
#      covmat.signal <- varcov.spatial(coords = coords[ind.rep[[i]],],
#                                      cov.model = cov.model, 
#                                      kappa = kappa, nugget = 0,
#                                      cov.pars = c(sigmasq, phi))$varcov
      spatial.comp[[i]] <- as.vector(varcov.spatial(coords = coords[ind.rep[[i]],],
                                                    cov.model = cov.model, 
                                                    kappa = kappa, nugget = 0,
                                                    cov.pars = c(sigmasq, phi))$varcov %*%
                                     varcov.spatial(coords = coords[ind.rep[[i]],],
                                                    cov.model = cov.model, 
                                                    kappa = kappa, nugget = tausq,
                                                    cov.pars = c(sigmasq, phi), inv=TRUE)$inverse %*%
                                     res[ind.rep[[i]]]) 
    }
    spatial.comp <- as.vector(unlist(spatial.comp))[as.vector(unlist(ind.rep))]
    predict.comp <- trend.comp + spatial.comp
    residual.comp <- as.vector(temp.list$z - predict.comp)
#    residual.std <- as.vector(invcov %*% residual.comp)
#    residual.trend.std <- as.vector(invcov %*% res)
    lik.results$model.components <-
      data.frame(trend = trend.comp, spatial = spatial.comp, residuals = residual.comp)
#    lik.results$s2.random <- (crossprod(res,invcov) %*% res)/(n - beta.size)
#    lik.results$s2 <- (crossprod(residual.comp,invcov) %*% residual.comp)/(n - beta.size)
  }
  ##
  lik.results$contrasts <- xmat.contrasts
  lik.results$call <- call.fc
  ##
  ## Assigning classes
  ##
  oldClass(lik.results) <- c("likGRF", "variomodel")
  ##
  ## Some warning messages about particular possible results
  ##
  if(messages.screen){
    if((lik.results$cov.pars[1] < (0.01 * (lik.results$nugget + lik.results$cov.pars[1])))& lik.results$cov.pars[2] > 0)
      cat("\nWARNING: estimated sill is less than 1 hundredth of the total variance. Consider re-examine the model excluding spatial dependence\n" )      
    if((lik.results$cov.pars[2] > (10 * max.dist)) & lik.results$cov.pars[1] > 0 )
      cat("\nWARNING: estimated range is more than 10 times bigger than the biggest distance between two points. Consider re-examine the model:\n 1) excluding spatial dependence if estimated sill is too low and/or \n 2) taking trends (covariates) into account\n" ) 
    if(((lik.results$cov.pars[2] < (0.1 * min.dist)) & (lik.results$cov.pars[1] > 0)) & lik.results$cov.pars[2] > 0)
      cat("\nWARNING: estimated range is less than 1 tenth of the minimum distance between two points. Consider re-examine the model excluding spatial dependence\n" ) 
  }
  ##
  attr(lik.results, "geodata") <- name.geodata
  return(lik.results)
}

".negloglik.GRF" <-
  function(pars, fp, ip, temp.list)
### pars : values for the parameters to be estimated
  ## sequence is c(phi, tausq, kappa, lambda, psiR, psiA, sigmasq)
### fixed pars: parameters considered fixed
### ind.pars : list indicating which are fixed and which are to be estimated
  ##
  ## Warning:
  ##  if fix.nugget = TRUE and nugget > 0 ,
  ## sigmasq should be passed and fp$nugget is the value of the nugget
  ## otherwise the RELATIVE nugget should be passed
{
  p <- temp.list$beta.size
  log.jacobian <- temp.list$log.jacobian
  ## Obligatory parameter:
  phi <- pars[1]
  ## Others
  if(ip$f.tausq){
    if(fp$tausq > 0){
      npars.min <- length(pars)
      sigmasq <- pars[npars.min]
    }
    else sigmasq <- 1
  }
  else sigmasq <- 1
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- pars[2]
  }
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- pars[2]
    psiA <- fp$psiA
  }
  if(ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- pars[2]
    psiA <- pars[3]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- pars[2]
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- pars[2]
    psiR <- fp$psiR
    psiA <- pars[3]
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- pars[2]
    psiR <- pars[3]
    psiA <- fp$psiA
  }
  if(ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- fp$kappa
    lambda <- pars[2]
    psiR <- pars[3]
    psiA <- pars[4]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- pars[3]
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- fp$lambda
    psiR <- pars[3]
    psiA <- fp$psiA
  }
  if(ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- fp$lambda
    psiR <- pars[3]
    psiA <- pars[4]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- pars[3]
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- pars[3]
    psiR <- fp$psiR
    psiA <- pars[4]
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- pars[3]
    psiR <- pars[4]
    psiA <- fp$psiA
  }
  if(ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- fp$tausq
    kappa <- pars[2]
    lambda <- pars[3]
    psiR <- pars[4]
    psiA <- pars[5]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- pars[3]
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- pars[3]
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- fp$lambda
    psiR <- pars[3]
    psiA <- pars[4]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- pars[3]
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- pars[3]
    psiR <- fp$psiR
    psiA <- pars[4]
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- pars[3]
    psiR <- pars[4]
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- fp$kappa
    lambda <- pars[3]
    psiR <- pars[4]
    psiA <- pars[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- fp$lambda
    psiR <- fp$psiR
    psiA <- pars[4]
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- fp$lambda
    psiR <- pars[4]
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & !ip$f.kappa & ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- fp$lambda
    psiR <- pars[4]
    psiA <- pars[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- pars[4]
    psiR <- fp$psiR
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- pars[4]
    psiR <- fp$psiR
    psiA <- pars[5]
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- pars[4]
    psiR <- pars[5]
    psiA <- fp$psiA
  }
  if(!ip$f.tausq & !ip$f.kappa & !ip$f.lambda & !ip$f.psiR & !ip$f.psiA){
    tausq <- pars[2]
    kappa <- pars[3]
    lambda <- pars[4]
    psiR <- pars[5]
    psiA <- pars[6]
  }
  ##
  if(temp.list$print.pars){
    running.pars <- c(phi = phi, tausq = tausq, kappa =kappa, psiA = psiA, psiR = psiR, lambda = lambda)
    if(ip$f.tausq && fp$tausq > 0)
      running.pars <- c(sigmasq=sigmasq, running.pars)
    print(running.pars)
  }
  ##
  ## Absurd values
  ##
  if(kappa < 1e-04 | (tausq+sigmasq) < (.Machine$double.eps^0.5) |
     any(c(phi, tausq, sigmasq, kappa) < 0))
    return(.Machine$double.xmax^0.5)
  ##
  ## Anisotropy
  ##
  if(!ip$f.psiR | !ip$f.psiA){
    coords.c <- coords.aniso(temp.list$coords, aniso.pars=c(psiA, psiR))
    vecdist <- function(x){as.vector(dist(x))}
    assign(".likGRF.dists.vec", lapply(split(as.data.frame(coords.c),
                                             temp.list$realisations), vecdist), envir=.geoR.env)
  }
  ##
  ## Box-Cox transformation
  ##
  if(!ip$f.lambda){
    if(abs(lambda - 1) < 0.0001) {
      log.jacobian <- 0
    }
    else {
      if(any(temp.list$z <= 0))
        stop("Transformation not allowed for zero or negative data")
      data <- temp.list$z^(lambda - 1)
      if(any(data <= 0)) log.jacobian <- log(prod(data))
      else log.jacobian <- sum(log(data))
      data <- NULL
    }
    if(abs(lambda) < 0.0001)
      data <- log(temp.list$z)
    else data <- ((temp.list$z^lambda) - 1)/lambda
  }
  else data <- temp.list$z
  data <- split(data, as.factor(temp.list$realisations))
  ##
  ## Computing likelihood
  ##
  sumnegloglik <- 0
  for(i in 1:temp.list$nrep){
    ## NOTE: Likelihood for Independent observations 
    ##       arbitrary criteria used here:
    ##       (phi < 1-e16) or (sigmasq < 1-e16)  ==> independence
    ##
    n <- temp.list$n[i]
    xmat <- temp.list$xmat[[i]]
    z <- data[[i]]
    if((phi < 1e-16) | (sigmasq < 1e-16)){
      if(ip$f.tausq)
        v <- list(varcov = diag(x=(tausq+sigmasq), n),
                  log.det.to.half = (n/2) * log(tausq+sigmasq))
      else
        v <- list(varcov = diag(x=(1+tausq), n),
                  log.det.to.half = (n/2) * log(1+tausq))
    }
    else
      v <- varcov.spatial(dists.lowertri = get(".likGRF.dists.vec", envir=.geoR.env)[[i]],
                          cov.model = temp.list$cov.model, kappa=kappa,
                          nugget = tausq, cov.pars=c(sigmasq, phi),
                          det = TRUE)
    if(!is.null(v$crash.parms)) return(.Machine$double.xmax^0.5)
    ivx <- solve(v$varcov,xmat)
    xivx <- crossprod(ivx,xmat)
    betahat <- try(.solve.geoR(xivx,crossprod(ivx,z)), silent=TRUE)
    if(inherits(betahat, "try-error")){
      t.ei <- eigen(xivx, symmetric = TRUE)
#      if(exists("trySilent"))
        betahat <- try(crossprod(t(t.ei$vec)/sqrt(t.ei$val)) %*% crossprod(ivx,z), silent=TRUE)
#      else{
#        error.now <- options()$show.error.message
#        options(show.error.messages = FALSE)
#        betahat <- try(crossprod(t(t.ei$vec)/sqrt(t.ei$val)) %*% crossprod(ivx,z))
#        if(is.null(error.now) || error.now) options(show.error.messages = TRUE)        
#      }
    }
    if(inherits(betahat, "try-error"))
      stop("Covariates have very different orders of magnitude. Try to multiply and/or divide them to bring them to similar orders of magnitude") 
    res <- z - xmat %*% betahat
    ssres <- drop(crossprod(res,solve(v$varcov,res)))
    if(temp.list$method.lik == "ML"){
      if(ip$f.tausq & (tausq > 0))
        negloglik <- v$log.det.to.half +  0.5 * ssres
      else
        negloglik <- (n/2) * log(ssres) +  v$log.det.to.half
    }
    if(temp.list$method.lik == "RML"){
      if(length(as.vector(xivx)) == 1) {
        choldet <- 0.5 * log(xivx)
      }
      else {
        chol.xivx <- chol(xivx)
        choldet <- sum(log(diag(chol.xivx)))
      }
      if(ip$f.tausq & (tausq > 0))
        negloglik <- v$log.det.to.half +  0.5 * ssres + choldet
      else
        negloglik <- ((n-p)/2) * log(ssres) +  v$log.det.to.half + choldet
    }  
    negloglik <- negloglik - temp.list$loglik.cte[i]
    sumnegloglik <- sumnegloglik + negloglik
  }
  sumnegloglik <- sumnegloglik - log.jacobian
  if(sumnegloglik > (.Machine$double.xmax^0.5) | sumnegloglik == Inf | sumnegloglik == -Inf)
    sumnegloglik <- .Machine$double.xmax^0.5
  if(temp.list$print.pars)
    cat(paste("log-likelihood = ", -sumnegloglik, "\n"))
  return(sumnegloglik) 
}

"likfit.limits" <- "pars.limits"

"print.likGRF" <-
  function (x, digits = max(3, getOption("digits") - 3), ...)
{
  est.pars <- as.vector(x$parameters.summary[x$parameters.summary[,1]
  == "estimated",2])
  names.est.pars <- dimnames(x$parameters.summary[x$parameters.summary[,1] == "estimated",])[[1]]
  names(est.pars) <- names.est.pars
  cat("likfit: estimated model parameters:\n")
  print.default(format(est.pars, digits = digits), ...)
  cat(paste("Practical Range with cor=0.05 for asymptotic range:", format(x$practicalRange, ...)))
  cat("\n")
  ##  print(round(est.pars, digits=digits))
  cat("\nlikfit: maximised log-likelihood = ")
  cat(format(x$loglik, digits = digits))
  ##  cat(round(x$loglik, digits=digits))
  cat("\n")
  return(invisible())
}  

"summary.likGRF" <-
  function(object, ...)
{
  names.pars <- dimnames(object$parameters.summary)[[1]]
  summ.lik <- list()
  if(object$method.lik == "ML")
    summ.lik$method.lik <- "maximum likelihood"
  if(object$method.lik == "RML")
    summ.lik$method.lik <- "restricted maximum likelihood"
  summ.lik$mean.component <- object$beta
  names(summ.lik$mean.component) <- names.pars[1:length(object$beta)]
  summ.lik$cov.model <- object$cov.model
  summ.lik$spatial.component <- object$parameters.summary[c("sigmasq", "phi"),]
  summ.lik$spatial.component.extra <- object$parameters.summary[c("kappa", "psiA", "psiR"),]
  summ.lik$nugget.component <- object$parameters.summary[c("tausq"),, drop=FALSE]
  summ.lik$transformation  <- object$parameters.summary[c("lambda"),, drop=FALSE]
  summ.lik$practicalRange <-  object$practicalRange
  summ.lik$likelihood <- list(log.L = object$loglik, n.params = as.integer(object$npars),
                               AIC = object$AIC, BIC = object$BIC)
  summ.lik$ns.likelihood <- list(log.L = object$nospatial$loglik.ns, n.params = object$nospatial$npars.ns, AIC=object$nospatial$AIC.ns, BIC=object$nospatial$BIC.ns)
  summ.lik$estimated.pars <- dimnames(object$parameters.summary[object$parameters.summary[,1] == "estimated",])[[1]]
  summ.lik$call <- object$call
  oldClass(summ.lik) <- "summary.likGRF"
  return(summ.lik)
}

"print.summary.likGRF" <-
  function (x, digits = max(3, getOption("digits") - 3), ...)
{
  if(length(class(x)) == 0 || all(class(x) != "summary.likGRF"))
    stop("object is not of the class \"summary.likGRF\"")
  cat("Summary of the parameter estimation\n")
  cat("-----------------------------------\n")
  cat(paste("Estimation method:", x$method.lik, "\n"))
  cat("\n")
  ##
  ## Estimates of the model components
  ## Model: Y(x) = X\beta + S(x) + e 
  ##
  cat("Parameters of the mean component (trend):")
  cat("\n")
  print.default(round(x$mean.component, digits = digits), ...)
  ##  print(round(x$mean.component, digits=digits))
  cat("\n")
  ##
  cat("Parameters of the spatial component:")
  cat("\n")
  cat(paste("   correlation function:", x$cov.model))
  cat(paste("\n      (estimated) variance parameter sigmasq (partial sill) = ", format(x$spatial.component[1,2], digits=digits)))
  cat(paste("\n      (estimated) cor. fct. parameter phi (range parameter)  = ", format(x$spatial.component[2,2], digits=digits)))
  if(any(x$cov.model == c("matern", "powered.exponential",
     "cauchy", "gencauchy", "gneiting.matern"))){
    kappa <- x$spatial.component.extra["kappa",2]
    if(x$spatial.component.extra["kappa",1] == "estimated")
      cat(paste("\n      (estimated) extra parameter kappa =", format(kappa, digits=digits)))
    else{
      cat(paste("\n      (fixed) extra parameter kappa = ", kappa))
      if(x$cov.model == "matern" & (format(kappa, digits=digits)  == 0.5))
      cat(" (exponential)")
    }
  }
  cat("\n")
  ##
  aniso <-  x$spatial.component.extra[c("psiA", "psiR"),]
  psiApsiR <- x$spatial.component.extra[c("psiA", "psiR"),2]
  cat("   anisotropy parameters:")
  if(aniso["psiA",1] == "estimated")
    cat(paste("\n      (estimated) anisotropy angle =",
              format(psiApsiR[1], digits=digits),
              " (",format((psiApsiR[1]*360)/(2*pi), dig=1), "degrees )"))
  else
    cat(paste("\n      (fixed) anisotropy angle =", psiApsiR[1],
              " (",(psiApsiR[1]*360)/(2*pi), "degrees )"))
  if(aniso["psiR",1] == "estimated")
    cat(paste("\n      (estimated) anisotropy ratio =",
              format(psiApsiR[2], digits=digits)))
  else
    cat(paste("\n      (fixed) anisotropy ratio =", psiApsiR[2]))
  cat("\n")
  cat("\n")  
  cat("Parameter of the error component:")
  if(x$nugget.component[,1] == "estimated")
    cat(paste("\n      (estimated) nugget = ", format(x$nugget.component[,2], digits=digits)))
  else
    cat(paste("\n      (fixed) nugget =", x$nugget.component[,2]))
  cat("\n")
  cat("\n")
  cat("Transformation parameter:")
  cat("\n")
  lambda <- x$transformation[,2]
  if(x$transformation[,1] == "estimated")
    cat(paste("      (estimated) Box-Cox parameter =", format(lambda, digits=digits)))
  else{
    cat(paste("      (fixed) Box-Cox parameter =", lambda))
    if(abs(lambda - 1) <  0.0001) cat(" (no transformation)")
    if(abs(lambda) < 0.0001) cat(" (log-transformation)")
  }
  cat("\n")
  cat("\n")
  cat("Practical Range with cor=0.05 for asymptotic range:",
      format(x$practicalRange, ...))
  cat("\n")
  cat("\n")
  cat("Maximised Likelihood:")
  cat("\n")
  print(format(x$likelihood, digits=digits), ...)
  cat("\n")
  cat("non spatial model:")
  cat("\n")
  print(format(x$ns.likelihood, digits=digits), ...)
  cat("\n")
  cat("Call:")
  cat("\n")
  print(x$call)
  cat("\n")
  invisible(x)
}

"loglik.GRF" <-
  function(geodata, coords=geodata$coords, data=geodata$data,
           obj.model = NULL,
           cov.model="exp", cov.pars,
           nugget=0, kappa=0.5, lambda=1, psiR=1, psiA=0,
           trend="cte", method.lik="ML",
           compute.dists = TRUE, realisations = NULL)
{
  if(!is.null(obj.model)){
    if(!is.null(obj.model$cov.model)) cov.model <- obj.model$cov.model
    if(!is.null(obj.model$cov.pars)) cov.pars <- obj.model$cov.pars
    if(!is.null(obj.model$nugget)) nugget <- obj.model$nugget
    if(!is.null(obj.model$kappa)) kappa <- obj.model$kappa
    if(!is.null(obj.model$lambda)) lambda <- obj.model$lambda
    if(!is.null(obj.model$psiR)) psiR <- obj.model$psiR
    if(!is.null(obj.model$psiA)) psiA <- obj.model$psiA      
   if(!is.null(obj.model$trend)) trend <- eval(obj.model$trend)
  ## a resolver: problema em passando  trend
  }
  sigmasq <- cov.pars[1]
  phi <- cov.pars[2]
  if(method.lik == "REML" | method.lik == "reml" | method.lik == "rml") 
    method.lik <- "RML"
  if(method.lik == "ML" | method.lik == "ml")
    method.lik <- "ML"
  if(is.null(realisations))
    realisations <- as.factor(rep(1, length(data)))
  else
    realisations <- as.factor(realisations)
  nrep <- length(levels(realisations))
  ##
  ## Absurd values
  ##
  if(kappa < 1e-04) return(-(.Machine$double.xmax^0.5))
  if((nugget+sigmasq) < 1e-16) return(-(.Machine$double.xmax^0.5))
  ##
  ## Trend matrix
  ##
  if(missing(geodata))
    xmat <- unclass(trend.spatial(trend=trend, geodata = list(coords = coords, data = data)))
  else
    xmat <- unclass(trend.spatial(trend=trend, geodata = geodata))
  if (nrow(xmat) != nrow(coords)) 
    stop("coords and trend have incompatible sizes")
  beta.size <- ncol(xmat)
  xmat <- split(as.data.frame(unclass(xmat)), realisations)
  xmat <- lapply(xmat, as.matrix)
  ##
  ## Anisotropy
  ##
  vecdist <- function(x){as.vector(dist(x))}
  if(psiR != 1 | psiA != 0){
    coords.c <- coords.aniso(coords, aniso.pars=c(psiA, psiR))
    .likGRF.dists.vec <- lapply(split(as.data.frame(coords.c),
                                      as.factor(realisations)), vecdist)
  }
  else if(compute.dists) .likGRF.dists.vec <- lapply(split(as.data.frame(coords),
                                                           as.factor(realisations)), vecdist)
  ##
  ## Box-Cox transformation
  ##
  z <- data
  if(abs(lambda - 1) < 0.0001)
    log.jacobian <- 0
  else {
    if(any(z <= 0))
      stop("Transformation not allowed for zero or negative data")
    data <- z^(lambda - 1)
    if(any(data <= 0)) log.jacobian <- log(prod(data))
    else log.jacobian <- sum(log(data))
    data <- NULL
    if(abs(lambda) < 0.0001)
      data <- log(z)
    else data <- ((z^lambda) - 1)/lambda
  }
  data <- split(data, as.factor(realisations))
  ##
  ## Computing likelihood
  ##
  sumnegloglik <- 0
  for(i in 1:nrep){
    ## NOTE: Likelihood for Independent observations 
    ##       arbitrary criteria used here:
    ##       (phi < 1-e16) or (sigmasq < 1-e16)  ==> independence
    ##
    n <- length(data[[1]])
    if((phi < 1e-16) | (sigmasq < 1e-16)){
      V <- list(varcov = diag(x=(nugget+sigmasq), n),
                log.det.to.half = (n/2) * log(nugget+sigmasq))
    }
    else{
      V <- varcov.spatial(dists.lowertri = .likGRF.dists.vec[[i]],
                          cov.model = cov.model, kappa=kappa,
                          nugget = nugget, cov.pars=c(sigmasq, phi),
                          det = TRUE)
    }
    if(!is.null(V$crash.parms)){
      cat("varcov.spatial: improper matrix for following the given parameters:")
      print(V$crash.parms)
      stop()
    }
    ivx <- solve(V$varcov,xmat[[i]])
    xivx <- crossprod(ivx,xmat[[i]])
    betahat <- .solve.geoR(xivx, crossprod(ivx,data[[i]]))
    res <- data[[i]] - xmat[[i]] %*% betahat
    ssres <- drop(crossprod(res, solve(V$varcov,res)))
    if(method.lik == "ML"){
      negloglik <- (n/2)*(log(2*pi)) + V$log.det.to.half +  0.5 * ssres
    }
    if(method.lik == "RML"){
      choldet <- sum(log(diag(chol(xivx))))
      negloglik <- V$log.det.to.half +  0.5 * ssres + choldet
      xx.eigen <- eigen(crossprod(xmat[[i]]), symmetric = TRUE, only.values = TRUE)
      negloglik <- negloglik + ((n-beta.size)/2)*(log(2*pi)) - 0.5 * sum(log(xx.eigen$values))
    }
    sumnegloglik <- sumnegloglik + negloglik
  }
  sumnegloglik <- sumnegloglik - log.jacobian
  if(sumnegloglik > (.Machine$double.xmax^0.5))
    sumnegloglik <- .Machine$double.xmax^0.5
  return(as.vector(-sumnegloglik))
}


##
## ----------------
## Other functions
## ----------------
##
##

"likfit.nospatial" <-
  function(temp.list, ...)
{
  results <- list()
  z <- temp.list$z
  n <- temp.list$n
  beta.size <- temp.list$beta.size
  xmat <- temp.list$xmat
  txmat <- temp.list$txmat
  ixx <- solve(crossprod(xmat))
  if(temp.list$fix.lambda == FALSE){
    if (temp.list$minimisation.function == "nlm"){
      assign(".temp.lower.lambda",-2, envir=.geoR.env)
      assign(".temp.upper.lambda", 2, envir=.geoR.env)
      results <- nlm(.proflik.lambda, 1, ...)
      if(exists(".temp.lambda", envir=.geoR.env)){
        results$lambda <- get(".temp.lambda", envir=.geoR.env)
        remove(".temp.lambda", envir=.geoR.env, inherits = TRUE)
      }
      else{
        results$lambda <- results$estimate
      }
      remove(".temp.lower.lambda", ".temp.upper.lambda", inherits = TRUE, pos=1)
    }
    if (temp.list$minimisation.function == "nlmP"){
      results <- .nlmP(.proflik.lambda, 1, lower=-2, upper=2,...)  
      results$lambda <- results$estimate
    }
    if (temp.list$minimisation.function == "optim"){
      results <- optim(1, .proflik.lambda, method="L-BFGS-B", lower=-2, upper=2,...)
      results$minimum <- results$value
      results$lambda <- results$par
    }
    if(results$lambda == 1) {
      temp.list$log.jacobian <- 0
    }
    else {
      if(any(z <= 0))
        stop("Transformation option not allowed when there are zeros or negative data")
      if(any(z^(results$lambda - 1) <= 0))
        temp.list$log.jacobian <- log(prod(z^(results$lambda - 1)))
      else temp.list$log.jacobian <- sum(log(z^(results$lambda - 1)))
      if(results$lambda == 0)
        z <- log(z)
      else z <- ((z^results$lambda) - 1)/results$lambda
    }
  }
  else{
    results$lambda <- temp.list$lambda
    results$code <- 1
    if (temp.list$minimisation.function == "optim") results$convergence <- 0
  }
  ssres <- (z %*% (diag(n) - xmat %*%
                   solve(crossprod(xmat)) %*% txmat) %*% z)
  if(temp.list$method == "ML"){
    results$tausqhat <- ssres/n
    if(temp.list$fix.lambda)
      results$minimum <- as.vector(((n/2) * log(2 * pi) +
                          (n/2) * log(results$tausqhat) +
                          (n/2)  -
                          temp.list$log.jacobian))
  }
  if(temp.list$method == "RML") {
    results$tausqhat  <- (ssres/(n-beta.size))
    if(temp.list$fix.lambda){
      results$minimum <- as.vector((((n - beta.size)/2) * log(2 * pi) +
                          ((n - beta.size)/2) * log(results$tausqhat) +
                          (n/2) -
                          temp.list$log.jacobian
                          ))
    }
  }
  if (temp.list$minimisation.function == "optim") results$value <- results$minimum    
  return(results)
}

"loglik.spatial" <- function(pars)
{
  tausq <- pars[1]
  sigmasq <- pars[2]
  sill.total <- tausq + sigmasq
  phi <- pars[3]
  lambda <- pars[4]
  z <- .temp.list$z
  n <- .temp.list$n
  if(!isTRUE(.temp.list$fix.lambda)) {
    if(lambda == 1) {
      .temp.list$log.jacobian <- 0
    }
    else {
      if(any(z < 0))
        stop("Transformation option not allowed when there are zeros or negative data"
             )
      if(any(z^(lambda - 1) <= 0))
        .temp.list$log.jacobian <- log(prod(z^(lambda -
                                               1)))
      else .temp.list$log.jacobian <- sum(log(z^(lambda - 1)))
      if(lambda == 0)
        z <- log(z)
      else z <- ((z^lambda) - 1)/lambda
    }
  }
  beta.size <- .temp.list$beta.size
  kappa <- .temp.list$kappa
#  covinf <- varcov.spatial(dists.lowertri = .temp.list$
#                           dists.lowertri, cov.model = .temp.list$cov.model,
#                           kappa = kappa, nugget = tausq,
#                           cov.pars = c(sigmasq, phi), scaled = FALSE,
#                           inv = TRUE, det = TRUE,
#                           only.inv.lower.diag = TRUE)
  COV <- varcov.spatial(dists.lowertri = .temp.list$
                        dists.lowertri, cov.model = .temp.list$cov.model,
                        kappa = kappa, nugget = tausq,
                        cov.pars = c(sigmasq, phi), scaled = FALSE,
                        det = TRUE)
#  xix <- as.double(rep(0, beta.size*beta.size))
#  xix <- .C("bilinearform_XAY",
#            as.double(covinf$lower.inverse),
#            as.double(covinf$diag.inverse),
#            as.double(as.vector(.temp.list$xmat)),
#            as.double(as.vector(.temp.list$xmat)),
#            as.integer(beta.size),
#            as.integer(beta.size),
#            as.integer(n),
#            res = xix, PACKAGE = "geoR")$res
#  attr(xix, "dim") <- c(beta.size, beta.size)
  xix <- crossprod(.temp.list$xmat, solve(COV$varcov,
                                           .temp.list$xmat))
  if(length(as.vector(xix)) == 1) {
    choldet <- 0.5 * log(xix)
  }
  else {
    choldet <- sum(log(diag(chol(xix))))
  }
  iy <- solve(COV$varcov, z)
  xiy <- crossprod(.temp.list$xmat, iy)
  beta.hat <- drop(solve(xix,xiy))
  yiy <- drop(crossprod(z, iy))
  ssresmat <- as.vector(yiy - 2*crossprod(beta.hat,xiy) +  beta.hat %*% xix %*% beta.hat)
  if(.temp.list$method == "ML") {
    loglik <- ( - (n/2) * log(2 * pi) -
               COV$log.det.to.half -
               0.5 * ssresmat + 
               .temp.list$log.jacobian)
  }
  if(.temp.list$method == "RML") {
    xx.eigen <- eigen(crossprod(.temp.list$xmat), symmetric = TRUE, only.values = TRUE)
    loglik <- ( - ((n - beta.size)/2) * log(2 * pi) +
               0.5 * sum(log(xx.eigen$values)) -
               COV$log.det.to.half -
               (0.5) * ssresmat -
               choldet + .temp.list$log.jacobian)
  }
  return(as.vector(loglik))
}

"logLik.likGRF" <- function(object, ...){
  val <- object$loglik
  attr(val, "df") <- object$npars
  attr(val, "method") <- object$method.lik
  class(val) <- "logLik"
  return(val)
}

"fitted.likGRF" <- 
  function(object, spatial = TRUE, ...)
{
  if(is.null(object$model.components)){
    object$call$components <- TRUE
    object$call$messages <- FALSE
    object$call$ini.cov.pars <- object$cov.pars
    object$call$nugget <- object$nugget
    object$call$lambda <- object$lambda
    cat("be patient ... this function currently require calling likfit again\n")
    object <- eval(object$call)
  }
  if(spatial)
    return(rowSums(object$model.components[,c("trend", "spatial")]))
  else
    return(object$model.components[,"trend"]) 
}

"resid.likGRF" <- "residuals.likGRF" <-
  function(object, spatial = FALSE, ...)
{
  if(is.null(object$model.components)){
    object$call$components <- TRUE
    object$call$messages <- FALSE
    object$call$ini.cov.pars <- object$cov.pars
    object$call$nugget <- object$nugget
    object$call$lambda <- object$lambda
    cat("be patient ... this function currently require calling likfit again\n")
    object <- eval(object$call)
  }
  if(spatial)
    return(rowSums(object$model.components[,c("spatial","residuals")]))
  else
    return(object$model.components[,"residuals"])
}

Try the geoR package in your browser

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

geoR documentation built on Aug. 9, 2022, 5:11 p.m.