R/aovSar.rcbd.R

In lrcastro/spANOVA: Spatial Analysis of Field Trials Experiments using Geostatistics and Spatial Autoregressive Model

#' @name aovSar.rcbd
#'
#' @title Using a SAR model to handle spatial dependence in a Randomized Complete Block Design
#' @description Fit a randomized complete block design when the experimental units have some degree of
#' spatial dependence using a Spatial Lag Model (SAR).
#' @usage aovSar.rcbd(resp, treat, block, coord, seq.radius)
#'
#' @param resp Numeric or complex vector containing the values of response variable.
#' @param treat Numeric or complex vector containing the treatment applied to each experimental unit.
#' @param block Numeric or complex vector specifying the blocks.
#' @param coord Matrix of point coordinates or a SpatialPoints Object.
#' @param seq.radius Complex vector containing a radii sequence used to set the neighborhood pattern.
#' The default sequence has ten numbers from 0 to half of the maximum distance between the samples.
#'
#' @return \code{aovSar.rcbd} returns an object of \code{\link[base]{class}} "SARanova".
#' The functions summary and anova are used to obtain and print a summary and analysis of variance
#' table of the results.
#' An object of class "SARanova" is a list containing the following components:
#'
#' @details
#' Three assumptions are made about the error in the analysis of variance (ANOVA):
#'
#' 1. the errors are normally distributed and, on average, zero;
#'
#' 2. the errors all have the same variance (they are homoscedastic), and
#'
#' 3. the errors are unrelated to each other (they are independent across observations).
#'
#'When these assumptions are not satisfied, data transformations in the response variable are
#'often used to circumvent this problem. For example, in absence of normality, the Box-Cox
#'transformation can be used.
#'
#'However, in many experiments, especially field trials, there is a type of correlation generated by
#'the sample locations known as spatial correlation, and this condition violates the independence assumption.
#'In this setting, this function provides an alternative for using ANOVA when the errors are spatially
#'correlated, by using a data transformation discussed in Long (1996)
#'
#' \deqn{Y_{adj} = Y - (\hat{\rho}WY - \hat{\rho}\beta_0),}
#'
#'where \eqn{\hat{\rho}} denotes the autoregressive spatial parameter of the SAR model estimated by
#'lagsarlm, \eqn{\beta_0} is the overall mean and \eqn{W} is a spatial neighborhood matrix which neighbors are defined as the
#'samples located within a radius, this radius is specified as a sequence in \code{seq.radius}. For each radius
#'in \code{seq.radius} the model is computed as well its AIC, then the radius chosen is the one
#'that minimizes AIC.
#'
#' The aim of this transformation is converting autocorrelated observations into non-correlated observations
#' in order to apply the analysis of variance and obtain suitable inferences.
#'
#' @return
#' \item{DF}{degrees of freedom of rho, treatments, residual and total.}
#' \item{SS}{sum of squares of rho, treatments and residual.}
#' \item{Fc}{F statistic calculated for treatment.}
#' \item{residuals}{residuals of the adjusted model.}
#' \item{p.value}{p-value associated to F statistic for treatment.}
#' \item{rho}{the autoregressive parameter.}
#' \item{Par}{data.frame with the radii tested and its AIC.}
#' \item{y_orig}{vector of response.}
#' \item{y_ajus}{vector of adjusted response.}
#' \item{treat}{vector of treatment applied to each experimental unit.}
#' \item{modelAdj}{model of class \code{\link[stats]{aov}} using the adjusted response.}
#' \item{modelstd}{data frame containing the ANOVA table using non-adjusted response.}
#' \item{namey}{response variable name.}
#' \item{namex}{treatment variable name.}
#'
#' @references
#' Long, D.S., 1996. Spatial statistics for analysis of variance of agronomic field trials. In:
#' Arlinghaus, S.L. (Ed.), Practical Handbook of Spatial Statistics. CRC Press, Boca Raton, FL,
#' pp. 251–278.
#'
#' Rossoni, D. F.; Lima, R. R. . Autoregressive analysis of variance for experiments with spatial
#' dependence between plots: a simulation study. Revista Brasileira de Biometria, 2019
#'
#' Scolforo, Henrique Ferraço, et al. "Autoregressive spatial analysis and individual
#' tree modeling as strategies for the management of Eremanthus erythropappus." Journal of
#' forestry research 27.3 (2016): 595-603.
#'
#'
#' @examples
#' data("rcbd_simulated")
#'
#' # Fitting the model
#' model <- aovSar.rcbd(rcbd_simulated$y, rcbd_simulated$trat, rcbd_simulated$block,
#'                      cbind(rcbd_simulated$coordX, rcbd_simulated$coordY))
#'
#' # Summary for class SARanova
#' summary(model)
#'
#' # Anova for class SARanova
#' anova(model)
#'
#' @importFrom gtools stars.pval
#' @export
aovSar.rcbd <- function(resp, treat, block, coord, seq.radius) {

  # Defensive programming
  if(!(is.vector(resp) | is.numeric(resp))) {
    stop("'resp' must be a vector or numeric")
  }

  if(!(is.vector(treat) | is.numeric(treat))) {
    stop("'treat' must be a vector or numeric")
  }

  if(!(is.vector(block) | is.numeric(block))) {
    stop("'block' must be a vector or numeric")
  }

  if(!(is.matrix(coord) | class(coord)=="SpatialPoints")) {
    stop("'coord' must be a matrix or SpatialPoints object")
  }

  if(ncol(coord) < 2){
    stop("'coord' must have at least two columns")
  }

  if(missing(seq.radius)){
    max.dist <- max(dist(coord))
    seq.radius <- seq(0, 0.5*max.dist, l = 11)[-1]
  }

  params <- data.frame(radius = 0, rho = 0, AIC = 0)
  anova.list <- list()
  n <- length(resp)
  p.radius <- length(seq.radius)
  Y_ajus <- NULL
  treat <- factor(treat)
  block <- factor(block)

  for (i in 1:p.radius) {
    nb <- dnearneigh(coord, 0, seq.radius[i])
    w <- try(nb2mat(nb, style = "W"), silent = TRUE)
    test <- grepl("Error", w)

    # Se caso nao forem encontradas amostras dentro do raio especificado
    k <- 0.1 # incremento
    while(test[1] == TRUE){
      seq.radius <- seq(0, (0.5+k)*max.dist, l = 11)[-1]
      nb <- dnearneigh(coord, 0, seq.radius[i])
      w <- try(nb2mat(nb, style = "W"), silent = TRUE)
      test <- grepl("Error", w)
      k <- k + 0.1
    }

    listw <- nb2listw(nb, glist = NULL, style = "W")

    # SAR model
    SAR <- lagsarlm(resp ~ treat + block, listw = listw,
                    method = "eigen", tol.solve = 1e-15)
    ajuste <- summary(SAR)
    rho <- as.numeric(ajuste["rho"]$rho)
    params[i, ] <- c(raio = seq.radius[i], rho = rho, AIC = AIC(SAR))
  }

  # Adjusting the data and constructing the ANOVA table
  best.par <- which.min(params$AIC)
  beta <- mean(resp)
  nb <- dnearneigh(coord, 0, seq.radius[best.par])
  w <- nb2mat(nb, style = "W")
  Y_ajus <- resp - (params[best.par,"rho"] * w%*%resp - params[best.par,"rho"] * beta)
  aov.cl <- anova(aov(resp ~ treat + block))
  model.adj <- aov(Y_ajus ~ treat + block)
  aov.adj <- anova(model.adj)
  Sqt.nadj <- sum(aov.cl[,2])

  #Degres of freedom
  gltrat <- aov.adj[1][[1]][1]
  glblock <- aov.adj[1][[1]][2]
  glerror <- aov.adj[1][[1]][3]

  gltot <- sum(gltrat, glblock, glerror)

  #Sum of squares
  sqtrat <- aov.adj[2][[1]][1]
  sqblock <- aov.adj[2][[1]][2]
  sqerror <- aov.adj[2][[1]][3]
  sqtot <- sum(sqtrat, sqblock, sqerror)
  sqtotcor <- Sqt.nadj - sqtot

  #Mean Squares
  mstrat <- sqtrat/gltrat
  msblock <- sqblock/glblock
  mserror <- sqerror/glerror

  #F statistics
  ftrat <- mstrat/mserror
  fblock <- msblock/mserror
  pvalue.trat <- pf(ftrat, gltrat, glerror, lower.tail = FALSE)
  pvalue.block <- pf(fblock, glblock, glerror, lower.tail = FALSE)
  name.y <- paste(deparse(substitute(resp)))
  name.x <- paste(deparse(substitute(treat)))

  outpt <- list(DF = round(c(gltrat, glblock, glerror, gltot),0),
                SS = c(sqtrat, sqblock, sqerror, sqtotcor),
                MS = c(mstrat, sqblock, mserror),
                Fc = c(ftrat, fblock),
                residuals = resid(model.adj),
                p.value = c(pvalue.trat, pvalue.block),
                rho = params[best.par,"rho"],
                Par = params, y_orig = resp, y_ajus =  Y_ajus,
                treat = treat , modelAdj = model.adj, modelstd = aov.cl,
                namey = name.y, namex = name.x)
  class(outpt)<-c("SARanova","SARrcbd",class(aov.adj))
  return(outpt)
}


# Print method for this class
#' @export
#' @method print SARrcbd
print.SARrcbd <- function(x, ...) {
  cat("Response: ", x$namey, "\n")
  cat("Terms:","\n")
  trm <- data.frame(treat = c(as.character(round(x$SS[1],3)),as.character(x$DF[1])),
                    block = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])),
                    Residuals = c(as.character(round(x$SS[3],3)),as.character(x$DF[3])))
  rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
  print(trm)
  rse <- sqrt(x$MS[3])
  cat("\n")
  cat("Residual standard error:",rse)
  cat("\n")
  cat("Spatial autoregressive parameter:", x$rho,"\n")
  cat("Samples considered neighbor within a",x$Par[which.min(x$Par[,3]),1],"units radius")
}


# Summary method for this class
#' @export
#' @method summary SARrcbd
summary.SARrcbd <- function(object, ...) {
  cat("      Summary of RCBD","\n","\n")
  cat("Parameters tested:","\n")
  print(object$Par)
  cat("\n")
  cat("Selected parameters:","\n")
  print(object$Par[which.min(object$Par[,3]),])
}

# Anova method for this class
#' @export
#' @method anova SARrcbd
anova.SARrcbd <- function(object, compare = FALSE, ...) {

  if(is.logical(compare) == FALSE){
    warning("'compare' must be logical. Assuming compare == FALSE")
    compare = FALSE
  }

  cat("Analysis of Variance With Spatially Correlated Errors","\n")
  cat("\n")
  cat("Response:", ifelse(length(object$namey)>1,"resp",object$namey), "\n")
  star <- stars.pval(object$p.value[1])
  star2 <- stars.pval(object$p.value[2])
  anova.p1 <- data.frame("DF" = round(object$DF[1:4],0),
                         "SS" = round(object$SS[1:4],4),
                         "MS" = c(round(object$MS[1:3],4), ""),
                         "Fc" = c(round(object$Fc[1],4), round(object$Fc[2],4), "", ""),
                         "Pv" = c(format.pval(object$p.value[1]), format.pval(object$p.value[2]), "", ""),
                         "St" = c(star, star2, "", "")
  )
  colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)", "")
  rownames(anova.p1) <- c("Treatment", "Block", "Residuals", "Corrected Total")
  print(anova.p1)
  cat("---","\n")
  cat("Signif. codes: ",attr(star, "legend"), "\n")

  if(compare){
    cat("\n", "\n")
    cat("---------------------------------------------------------------","\n")
    cat("Standard Analysis of Variance", "\n")
    cat("---------------------------------------------------------------")
    cat("\n")
    print(object$modelstd)
  }

  return(invisible(anova.p1))

}


#exportar a funcao stars.pval do pacote gtools