R/check_spatialcorr.R

In rscherrer/sagreicolor:

#' Test spatial autocorrelation of the residuals
#'
#' This function tests for positive spatial autocorrelation of the (uni- or multivariate) residuals, using either Moran's I or Mantel's test.
#'
#' @param specdata A data frame containing at least columns for the dependent variables, as well as a column "island" and a column "habitat".
#' @param vars A character or integer vector. The names, or indices, of the dependent variables in \code{specdata}.
#' @param plotit Logical. Whether to plot residuals versus latitude and longitude.
#' @param test A character, what test to perform. \code{"moran"} (default) performs Moran's I test on each dependent variable. \code{"mantel"} performs Mantel's test to compare distance matrices in residual space and geographic space, meaning that it can perform a single test on multivariate residuals if several dependent variables are input.
#' @param method A character, the method used for P-value correction. See \code{?p.adjust}.
#' @param model A character string, either \code{lm} or \code{gls}, defining what type of linear model to extract residuals from. If GLS is chosen, AIC comparison of variance structures is performed (one residual variance per island, per habitat or per island per habitat), and the Pearson residuals of the best model are extracted (note that Pearson residuals are identical to raw residuals for a regular linear model, whereas they differ for GLS models with particular variance structures -- they are adjusted such that heteroscedasticity is reduced).
#' @return Either the output of the \code{mantel.test} function from ape, or a table summarizing the (multiple) Moran's I test(s). If \code{plotit = T}, also returns a scatter plot of the residuals versus latitude and longitude.
#' @note Beware that Mantel's test is a permutation test (999 permutations), and will take longer tha Moran's I test to run. Also, note that latitudes and longitudes are used to calculate geographical distances at the surface of the Earth using the Haversine distance function implemented in the package \code{geosphere}.
#' @author Raphael Scherrer
#' @export

# Function to check (uni or multivariate) spatial autocorrelation of the residuals
# Mantel allows one test for multivariate dataset
# Moran performs several tests, one for each
check_spatialcorr <- function(specdata, vars, plotit = T, test = "moran", method = "bonferroni", model = "gls") {

  library(ape)
  library(geosphere)
  library(nlme)

  # Extract residuals
  residuals <- sapply(vars, function(curr.variable) {

    specdata$X <- specdata[,curr.variable]

    # Fit a linear model
    if(model == "lm") {

      mod <- lm(X ~ island*habitat, data = specdata)

    } else if (model == "gls") {

      # Or a generalized least square model
      # Design several variance structures first
      mod1 <- gls(X ~ island*habitat, data = specdata, weights = varIdent(form = ~1 | island * habitat))
      mod2 <- gls(X ~ island*habitat, data = specdata, weights = varIdent(form = ~1 | island))
      mod3 <- gls(X ~ island*habitat, data = specdata, weights = varIdent(form = ~1 | habitat))
      mod4 <- gls(X ~ island*habitat, data = specdata)

      # AIC comparison and best model selection
      aicValues <- AIC(mod1, mod2, mod3, mod4)$AIC
      bestIdx <- which(aicValues == min(aicValues))
      mod <- get(c("mod1", "mod2", "mod3", "mod4")[bestIdx])

    }

    # Extract residuals
    # Note: Pearson residuals are the same as raw residuals for a regular linear model. They do differ for a GLS with more complex variance structure.
    residuals <- resid(mod, type = "pearson")

  })

  # Plots
  if(plotit) {
    ylab <- "Residuals"
    xlab1 <- "Latitude"
    xlab2 <- "Longitude"
    par(mfrow = c(2, 2))
    apply(residuals, 2, function(residuals) {
      plot(residuals ~ specdata$latitude, ylab = ylab, xlab = xlab1)
      plot(residuals ~ specdata$longitude, ylab = ylab, xlab = xlab2)
    })
    par(mfrow = c(1,1))
  }


  # Distance matrix in residual space
  dR <- as.matrix(dist(residuals))

  # Distance matrix in geography
  dG <- as.matrix(distm(specdata[, c("longitude", "latitude")], fun = distHaversine))

  # Weigth matrix is the inverse of the distance matrix (needed for Moran's test)
  w <- 1/dG

  # Replace infinite values (same locality) with zeros
  # These will not be taken into account in Moran's I calculations
  w[w == Inf] <- 0

  # Perform Mantel's test...
  if(test == "mantel") {

    message("Performing Mantel's test...")

    res <- ape::mantel.test(dR, dG, graph = F, alternative = "greater")
    res <- with(res, c(Z = z.stat, p.value = p))

  } else if(test == "moran") {

    message("Performing Moran's I test...")

    # Or perform (multiple) Moran's I tests
    res <- apply(residuals, 2, function(residuals) {
      ape::Moran.I(residuals, weight = w, alternative = "greater")
    })

    res <- as.data.frame(do.call("rbind", res))

    res$padj <- p.adjust(res$p.value, method)

  }

  return(res)

}