R/msr.4thcorner.R

In sdray/adespatial: Multivariate Multiscale Spatial Analysis

#' Moran spectral randomization for fourth-corner analysis
#'
#' This function allows to test fourth-corner statistics using constrained null
#' models (for traits and/or environmental variables). If the argument
#' \code{phyloORorthobasis} is specified, random traits are
#' phylogenetically-constrained to preserve the global autocorrelation (Moran's
#' I) and the phylogenetic structures at multiple scales. If not, standard
#' permutations are used. If the argument \code{listwORorthobasis} is specified,
#' random environmental variables are spatially-constrained to preserve the
#' global autocorrelation (Moran's I) and the spatial structures at multiple
#' scales. If not, standard permutations are used. Multiscale property is
#' defined by the power spectrum (i.e. decomposition of the variance of the
#' original variables) on a basis of orthonormal eigenvectors (Moran's
#' Eigenvector Maps, MEM).
#'
#'
#' @param x An object generated by the \code{fourthcorner} 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 phyloORorthobasis an object of the class \code{phylo} (phylogeny)
#'   created by the functions of the \pkg{ape} package or an object of class
#'   \code{orthobasis} generated by functions of \pkg{adephylo}
#'   (\code{me.phylo})
#' @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{4thcorner} randomized replicates.
#' @author Stephane Dray \email{stephane.dray@@univ-lyon1.fr}
#' @seealso \code{\link{msr.default}}, \code{\link[adephylo]{me.phylo}}
#' @references Braga, J., Thuiller, W., ter Braak, C.J.F. and Dray, S. (2018)
#'   Integrating spatial and phylogenetic information in the fourth-corner
#'   analysis to test trait-environment relationships. Ecology, 99:2667-2674.
#'
#' @keywords spatial
#' @examples
#'
#' if(require("ade4", quietly = TRUE) & require("adephylo", quietly = TRUE)
#' & require("spdep", quietly = TRUE) & require("ape", quietly = TRUE)){
#' data(mafragh, package = "ade4")
#' fr1 <- fourthcorner(mafragh$env, mafragh$flo, mafragh$traits$tabQuantitative, nrepet = 49)
#' phy <- read.tree(text = mafragh$tre)
#' lw <- nb2listw(mafragh$nb)
#' fr1.msr <- msr(fr1, listwORorthobasis = lw, phyloORorthobasis = phy)
#'
#' fr1
#' fr1.msr
#' }
#'
#' @importFrom adephylo me.phylo
#' @importFrom ade4 as.krandtest combine.4thcorner
#' @export
msr.4thcorner <- function(x, listwORorthobasis, phyloORorthobasis, nrepet = x$npermut, method = c("pair", "triplet", "singleton"), ...){
    
    appel <- as.list(x$call)
    L <- eval.parent(appel$tabL)
    R <- eval.parent(appel$tabR)
    if(any(!apply(R, 2, is.numeric)))
        stop("Not implemented: table R contains non-numeric variables")
    Q <- eval.parent(appel$tabQ)
    if(any(!apply(Q, 2, is.numeric)))
        stop("Not implemented: table Q contains non-numeric variables")
    
    if(missing(listwORorthobasis)){
        ## just permute
        Rmsr <- replicate(nrepet, R[sample(nrow(R)),, drop = FALSE], simplify = FALSE)
    } else {
        Rmsr <- msr(R, listwORorthobasis, nrepet = nrepet, method = method, simplify = FALSE, ...)
    }    
    
    if(missing(phyloORorthobasis)){
        ## just permute
        Qmsr <- replicate(nrepet, Q[sample(nrow(Q)),, drop = FALSE], simplify = FALSE)
    } else {
        if (inherits(phyloORorthobasis, "phylo"))
            phyloORorthobasis <- me.phylo(phyloORorthobasis)
        Qmsr <- msr(Q, phyloORorthobasis, nrepet = nrepet, method = method, simplify = FALSE, ...)    
    }
    
    func.4th <- as.character(appel[[1]])  
    if(func.4th == "fourthcorner")
        elems <- c("tabD", "tabD2", "tabG")
    #    else if(func.4th == "fourthcorner2") TODO
    #        elems <- c("tabG", "trRLQ")
    else
        stop(paste("Not yet implemented for", func.4th))
    
    test.Rrand <- lapply(Rmsr, function(x) do.call(func.4th, list(tabR = as.data.frame(x),
        tabL = L, tabQ = Q , modeltype = 2, nrepet = 0)))
    
    test.Qrand <- lapply(Qmsr, function(x) do.call(func.4th, list(tabR = R,
        tabL = L, tabQ = as.data.frame(x) , modeltype = 2, nrepet = 0)))
    
    res.R <- vector(mode = "list", length = length(elems))
    for(i in 1:length(elems))
        res.R[[i]] <- matrix(NA, nrepet, length(test.Rrand[[1]][[elems[i]]]$obs))
    names(res.R) <- elems
    res.Q <- res.R 
    for(i in 1:nrepet){
        for(j in elems){
            res.R[[j]][i,] <- test.Rrand[[i]][[j]]$obs   
            res.Q[[j]][i,] <- test.Qrand[[i]][[j]]$obs  
        }
    }
    
    fc.R <- fc.Q <- x
    
    for(j in elems){
        if(inherits(x[[j]], "krandtest")){
            output <- ifelse(inherits(x[[j]], "lightkrandtest"), "light", "full")
            fc.R[[j]] <- as.krandtest(sim = res.R[[j]], obs = x[[j]]$obs, 
                alter = x[[j]]$alter, call = match.call(), names = x[[j]]$names,
                p.adjust.method = x[[j]]$adj.method, output = output)
            fc.Q[[j]] <- as.krandtest(sim = res.Q[[j]], obs = x[[j]]$obs, 
                alter = x[[j]]$alter, call = match.call(), names = x[[j]]$names,
                p.adjust.method = x[[j]]$adj.method, output = output)
            fc.Q[[j]]$statnames <- fc.R[[j]]$statnames <- x[[j]]$statnames
        }
    }
    
    res <- combine.4thcorner(fc.R, fc.Q)
    res$call <- match.call()
    res$npermut <- nrepet
    res$model <- "msr"
    return(res)
    
}