R/msr.varipart.R

In sdray/adespatial: Multivariate Multiscale Spatial Analysis

#'Moran spectral randomization for variation partitioning
#'
#'The functions allows to evaluate the significance and estimate parts in
#'variation partitioning using Moran Spectral Randomization (MSR) as a
#'spatially-constrained null model to account for spatial autocorrelation in
#'table X. Hence, this function provides a variation partioning adujsted for
#'spurious correlation due to spatial autocorrelation in both the response and
#'one explanatory matrix.
#'
#'@param x An object generated by the \code{varipart} function.
#'@param listwORorthobasis an object of the class \code{listw} (spatial weights)
#'  created by the functions of the \pkg{spdep} package or an object of class
#'  \code{orthobasis}
#'@param nrepet an \code{integer} indicating the number of replicates
#'@param method an character specifying which algorithm should be used to
#'  produce spatial replicates (see \code{\link{msr.default}}).
#'@param \dots further arguments of the \code{\link{msr.default}} function.
#'@return An object of class \code{varipart} randomized replicates.
#'
#'@details The function corrects the biases due to spatial autocorrelation by
#'  using MSR procedure to produce environmental predictors that preserve the
#'  spatial autocorrelation and the correlation structures of the original
#'  environmental variables while being generated independently of species
#'  distribution.
#'
#'
#'@author(s) Stephane Dray \email{stephane.dray@@univ-lyon1.fr} and Sylvie
#'Clappe \email{sylvie.clappe@@univ-lyon1.fr}
#'
#'@seealso \code{\link{msr.default}}, \code{\link[ade4]{varipart}}
#'
#'@references
#'
#'Clappe, S., Dray S. and P.R. Peres-Neto (2018) Beyond
#'neutrality: disentangling the effects of species sorting and spurious
#'correlations in community analysis. Ecology 99:1737-1747.
#'
#'Wagner, H. H., and S. Dray (2015). Generating spatially constrained null models
#'for irregularly spaced data using Moran spectral randomization methods.
#'Methods in Ecology and Evolution 6:1169-1178.
#'
#' @examples
#' library(ade4)
#' library(spdep)
#' data(mafragh)
#' ## Performing standard variation partitioning
#' dudiY <- dudi.pca(mafragh$flo, scannf = FALSE, scale = FALSE)
#' mafragh.lw <- nb2listw(mafragh$nb)
#' me <- mem(mafragh.lw, MEM.autocor = "positive")
#' vprda <- varipart(dudiY, mafragh$env, me, type = "parametric")
#'
#' ## Adjust estimation and compute p-value by msr methods
#' vprda.msr <- msr(vprda, mafragh.lw, nrepet=99)
#' vprda.msr
#'@importFrom ade4 as.krandtest
#'@importFrom stats as.formula lm.wfit model.frame
#'@export
msr.varipart <-
    function(x,
        listwORorthobasis,
        nrepet = x$test$rep[1],
        method = c("pair", "triplet", "singleton"),
        ...) {

        appel <- as.list(x$call)
        
        Y <- eval.parent(appel$Y)
        X <- data.frame(eval.parent(appel$X))
        W <- data.frame(eval.parent(appel$W))
        
        if (!inherits(Y, "dudi")) {
            scale <- FALSE ## default in varipart (appel$scale is NULL if not specified in the call)
            if (!is.null(appel$scale))
                scale <-  eval.parent(appel$scale)
            response.generic <- as.matrix(scalewt(Y, scale = scale))
            lw <- rep(1/NROW(Y), NROW(Y))
            sqlw <- sqrt(lw)
            sqcw <- sqrt(rep(1, NCOL(Y)))
            wt <- outer(sqlw, sqcw)
            inertot <- sum((response.generic * wt)^2)
        } else {
            inertot <- sum(Y$eig)
            lw <- Y$lw
            sqlw <- sqrt(lw)
            sqcw <- sqrt(Y$cw)
            response.generic <- as.matrix(Y$tab)
            wt <- outer(sqlw, sqcw)
        }
        

        Xmsr <- msr(X,
            listwORorthobasis,
            nrepet = nrepet,
            method = method, simplify = FALSE)
       
        WXmsr <- lapply(Xmsr, cbind, W)
        
        # 
        # fast computation of R2/adjusted using MSR procedure
        
        R2msrtest.QR <- function(df) {
            response.generic <- response.generic * wt
            isim <- c()
            for (i in 1:nrepet) {
                df[[i]] <- data.frame(df[[i]])
                mf <- model.matrix(~., df[[i]])
                x.expl <- scalewt(mf[, -1, drop = FALSE], scale = FALSE, wt = lw) * sqrt(lw)
                Q <- qr(x.expl, tol = 1e-06)
                Yfit.X <- qr.fitted(Q, response.generic)
                isim[i] <- sum(Yfit.X^2) / inertot
            }
            return(isim)
        }
        
        R2msrtest.lmwfit <- function(df) {
            isim <- c()
            for (i in 1:nrepet) {
                df[[i]] <- data.frame(df[[i]])
                colnames(df[[i]])[1:ncol(X)] <- colnames(X)
                fmla <-
                    as.formula(paste("response.generic ~", paste(colnames(df[[i]]), collapse = "+")))
                mf <-
                    model.frame(fmla, data = cbind.data.frame(response.generic, df[[i]]))
                mt <- attr(mf, "terms")
                xm <- model.matrix(mt, mf)
   
                ## Fast function for computing sum of squares of the fitted table
                isim[i] <-
                    sum((lm.wfit(
                        y = response.generic,
                        x = xm,
                        w = lw
                    )$fitted.values * wt) ^ 2) / inertot
            }
            
            
            return(isim)
        }
        
        ab.ini <- sum(x$R2[c("a", "b")])
        bc.ini <- sum(x$R2[c("b", "c")])
        
        R2msrtest <- R2msrtest.lmwfit
        if (identical(all.equal(lw, rep(1/length(lw), length(lw))), TRUE))
            R2msrtest <- R2msrtest.QR
        
        msr.ab <- R2msrtest(Xmsr)
        msr.abc <- R2msrtest(WXmsr)
        
        #  abc.ini <- sum(x$R2[c("a", "b", "c")])
        #  test <- as.krandtest(
        #     obs = c(ab.ini, abc.ini),
        #     sim = cbind(msr.ab, msr.abc),
        #     names = c("ab", "abc"),
        #     call = match.call()
        # )
        # 

         test <- as.randtest(
             obs = ab.ini,
             sim = msr.ab,
             call = match.call()
         )
         
        
        ab.adj <- 1 - (1 - ab.ini) / (1 - mean(msr.ab))
        a.adj <- 1 - (1 - x$R2["a"]) / (1 - mean(msr.abc - bc.ini))
        b.adj <- ab.adj - a.adj
        c.adj <- sum(x$R2.adj[c("b", "c")]) - b.adj
        d.adj <- 1 - (a.adj + b.adj + c.adj)
        R2.adj.msr <- c(a.adj, b.adj, c.adj, d.adj)
        names(R2.adj.msr) <- names(x$R2)
        res <-
            list(
                test = test,
                R2 = x$R2,
                R2.adj.msr = R2.adj.msr
                )
        class(res) <- c("varipart", "list")
        return(res)
    }