R/aovSar.crd.R

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

#' @name aovSar.crd
#'
#' @title Using a SAR model to handle spatial dependence in a Completely Randomized Design
#' @description Fit a completely randomized design when the experimental units have some degree of
#' spatial dependence using a Spatial Lag Model (SAR).
#' @usage aovSar.crd(resp, treat, coord, seq.radius)
#'
#' @param resp Numeric or complex vector containing the values of the response variable.
#' @param treat Numeric or complex vector containing the treatment applied to each experimental unit.
#' @param coord Matrix of point coordinates.
#' @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,
#' if no sample is found in this interval the sequence upper limit is increased by 10\% and so on.
#'
#' @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
#' \code{aovSar.crd} returns an object of \code{\link[base]{class}} "SARanova".
#' The functions \code{\link[base]{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:
#'
#' \item{DF}{degrees of freedom of rho, treatments, residual and total.}
#' \item{SS}{sum of squares of rho, treatments and residual.}
#' \item{MS}{mean square 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 radius 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("crd_simulated")
#' resp <- crd_simulated$y
#' treat <- crd_simulated$trat
#' coord <- cbind(crd_simulated$coordX, crd_simulated$coordY)
#' cv <- aovSar.crd(resp, treat, coord)
#'
#' #Summary for class SARanova
#' summary(cv)
#'
#' #Anova for class SARanova
#' anova(cv)
#'
#' @importFrom gtools stars.pval
#' @importFrom stats resid
#' @importFrom spdep nb2listw dnearneigh nb2mat
#' @importFrom spatialreg lagsarlm summary.sarlm print.sarlm anova.sarlm
#' @export
aovSar.crd <- function(resp, treat, 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.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)

  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, 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))
  model.adj <- aov(Y_ajus ~ treat)
  aov.adj <- anova(model.adj)
  Sqt.nadj <- sum(aov.cl[,2])

  #Degres of freedom
  gltrat <- aov.adj[1][[1]][1]
  glerror <- aov.adj[1][[1]][2]
  gltot <- sum(gltrat,glerror)

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

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

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

  outpt <- list(DF = round(c(gltrat, glerror, gltot),0),
                SS = c(sqtrat, sqerror, sqtotcor),
                MS = c(mstrat, mserror),
                Fc = c(ftrat),
                residuals = resid(model.adj),
                p.value = c(pvalue), 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","SARcrd",class(aov.adj))
  return(outpt)
}


# Print method for this class
#' @export
#' @method print SARcrd
print.SARcrd <- 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])),
                    Residuals = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])))
  rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
  print(trm)
  rse <- sqrt(x$MS[2])
  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 SARcrd
summary.SARcrd <- function(object, ...) {
  cat("      Summary of CRD","\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 SARcrd
anova.SARcrd <- 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)
  anova.p1 <- data.frame("DF" = round(object$DF[1:3],0),
                         "SS" = round(object$SS[1:3],4),
                         "MS" = c(round(object$MS[1:2],4), ""),
                         "Fc" = c(round(object$Fc,4), "", ""),
                         "Pv" = c(format.pval(object$p.value), "", ""),
                         "St" = c(star, "", "")
  )
  colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value" ,"Pr(>F)", "")
  rownames(anova.p1) <- c("Treatment","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