R/aem.R

In adespatial: Multivariate Multiscale Spatial Analysis

#' Construct asymmetric eigenvector maps (AEM)
#' 
#' This function constructs eigenvectors of a site-by-link matrix. Weights can 
#' be applied to the links.
#' 
#' @param aem.build.binary Object created by function \code{aem.build.binary}.
#' @param binary.mat Site (n rows) by link (k columns) matrix. The 1s in the 
#'   matrix represents the presence of a link influencing a site, directly or 
#'   indirectly, otherwise the values are 0s.
#' @param weight Vector of weights of length k, to be applied to the links.
#' @param rm.link0 Logical (\code{TRUE}, \code{FALSE}) determining if the links 
#'   directly connecting a site to the origin (site 0) should be removed. 
#'   Default value: \code{FALSE}. This parameter is only used when an object of 
#'   class \code{aem.build.binary} is provided to the function.
#' @param print.binary.mat Logical (\code{TRUE}, \code{FALSE}) determining if 
#'   the site-by-link matrix used in the analysis should be printed.
#'   
#' @details
#' 
#' If only an object of class \code{aem.build.binary} is given to this function, The
#' argument \code{binary.mat} is not considered. \code{binary.mat} is only 
#' considered when the argument \code{aem.build.binary} is missing.
#' 
#' If weights are applied to the links, the length of vector \code{weight} has 
#' to take into account wether the links connecting real sites to the origin 
#' (the fictitious site 0) have been kept or removed.
#' 
#' @return
#' 
#' \item{\code{value}}{A vector of singular values associated with the AEM.} 
#' \item{\code{vectors}}{A matrix of eigenvector. Each column is an AEM 
#' eigenfunction (or variable).} \item{\code{mod.binary.mat}}{A site-by-link 
#' matrix modified through the function.}
#' @references
#' 
#' Blanchet F.G., P. Legendre and Borcard D. (2008) Modelling directional 
#' spatial processes in ecological data. \emph{Ecological Modelling}, 215, 
#' 325-336.
#' 
#' @author F. Guillaume Blanchet
#'   
#' @note
#' 
#' It sometimes happens that AEM eigenfunctions have equal singular values. In 
#' that case, different sets of AEM eigenfunctions may be produced on different 
#' plateforms.
#' 
#' Eigenvectors associated to an eigenvalue that is smaller than \eqn{10^{-12}} 
#' are considered negligeable. They have been removed from the created AEM 
#' eigenfunctions.
#' 
#' @seealso  \code{\link{aem.build.binary}}, \code{\link{svd}}
#'   
#' @importFrom spdep cell2nb
#'   
#' @examples
#' # Construction of object of class nb (spdep)
#' if(require("spdep", quietly = TRUE)){
#' nb <- cell2nb(5,5,"queen")
#' 
#' # Create fictitious geographical coordinates 
#' xy <- cbind(1:25,expand.grid(1:5,1:5))
#' 
#' # Build binary site-by-link matrix
#' bin.mat <- aem.build.binary(nb,xy)
#' 
#' # Construct AEM eigenfunctions from an object of class aem.build.binary
#' res <- aem(aem.build.binary=bin.mat,rm.link0=FALSE)
#' res$values
#' 
#' # Illustrate 4 AEM eigenfunctions using bubble plots
#' opal <- palette()
#' palette(c("black","white"))
#' oldpar <- par(mfrow=c(2,2))
#' symbols(x=xy[,2:3], circles=abs(res$vectors[,1]), inches=FALSE, asp=1,
#'  fg=ifelse(sign(-res$vectors[,1])+1>0,1,0), 
#'  bg=ifelse(sign(res$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 1")
#' symbols(x=xy[,2:3], circles=abs(res$vectors[,2]), inches=FALSE, 
#' asp=1, fg=ifelse(sign(-res$vectors[,2])+1>0,1,0),
#'  bg=ifelse(sign(res$vectors[,2])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 2")
#' symbols(x=xy[,2:3], circles=abs(res$vectors[,3]), inches=FALSE, 
#' asp=1, fg=ifelse(sign(-res$vectors[,3])+1>0,1,0), 
#' bg=ifelse(sign(res$vectors[,3])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 3")
#' symbols(x=xy[,2:3], circles=abs(res$vectors[,4]), inches=FALSE, asp=1,
#'  fg=ifelse(sign(-res$vectors[,4])+1>0,1,0), 
#'  bg=ifelse(sign(res$vectors[,4])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 4")
#' 
#' # Construct AEM eigenfunctions using only a site-by-link matrix
#' res2 <- aem(binary.mat=bin.mat[[1]])
#' res2$values
#' 
#' # Illustrate 4 AEM eigenfunctions using bubble plots
#' par(mfrow=c(2,2))
#' symbols(x=xy[,2:3], circles=abs(res2$vectors[,1]), inches=FALSE, 
#' asp=1, fg=ifelse(sign(-res2$vectors[,1])+1>0,1,0), 
#' bg=ifelse(sign(res2$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 1")
#' symbols(x=xy[,2:3], circles=abs(res2$vectors[,2]), inches=FALSE,
#' asp=1, fg=ifelse(sign(-res2$vectors[,2])+1>0,1,0), 
#' bg=ifelse(sign(res2$vectors[,2])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 2")
#' symbols(x=xy[,2:3], circles=abs(res2$vectors[,3]), inches=FALSE,
#' asp=1, fg=ifelse(sign(-res2$vectors[,3])+1>0,1,0), 
#' bg=ifelse(sign(res2$vectors[,3])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 3")
#' symbols(x=xy[,2:3], circles=abs(res2$vectors[,4]), inches=FALSE,asp=1,
#'  fg=ifelse(sign(-res2$vectors[,4])+1>0,1,0), 
#'  bg=ifelse(sign(res2$vectors[,4])+1>0,1,0), xlab="x", ylab="y")
#' title("AEM 4")
#' 
#' palette(opal)
#' par(oldpar)
#' 
#' # Construct AEM eigenfunctions with a function of the distance
#' # as weights to put on the links
#' 
#' # Construction of object of class nb (spdep)
#' nb<-cell2nb(5,5,"queen")
#' 
#' # Create fictitious geographical coordinates
#' xy <- cbind(1:25,expand.grid(1:5,1:5))
#' 
#' # Build binary site-by-link matrix
#' bin.mat <- aem.build.binary(nb,xy)
#' 
#' # Construct a matrix of distances
#' long.lien.mat<-as.matrix(dist(xy))
#' 
#' # Extract the edges, remove the ones directly linked to site 0
#' lien.b<-bin.mat$edges[-1:-5,]
#' 
#' # Construct a vector giving the length of each edge
#' long.lien<-vector(length=nrow(lien.b))
#' 
#' for(i in 1:nrow(lien.b)){
#' 	long.lien[i]<-long.lien.mat[lien.b[i,1],lien.b[i,2]]
#' }
#' 
#' # Construct a vector of weights based on distance
#' weight.vec<-1-(long.lien/max(long.lien))^2
#' 
#' # Construct AEM eigenfunctions from an object of class aem.build.binary
#' res <- aem(aem.build.binary=bin.mat,weight=weight.vec,rm.link0=TRUE)
#' res
#' 
#' # Computing Moran's I for AEMs
#' 
#' # Building AEMs
#' xy <- cbind(1:25,expand.grid(1:5,1:5))
#' Wdist <- 1/as.matrix(dist(xy[,2:3]))
#' 
#' nb <- cell2nb(5,5,"queen")
#' bin.mat <- aem.build.binary(nb,xy)
#' linkBase <- bin.mat[[2]]
#' link <- linkBase[-which(linkBase[,1] == 0),]
#' weight <- numeric()
#' 
#' for(i in 1:nrow(link)){
#'    weight[i] <- Wdist[link[i,1],link[i,2]]
#' }
#' 
#' AEM <- aem(bin.mat, weight = weight, rm.link0 = TRUE)
#' 
#' # Constructing asymmetric matrix
#' matasym <- matrix(0,ncol=25, nrow=25)
#' 
#' for(i in 1:nrow(link)){
#'     matasym[link[i,1],link[i,2]]<- weight[i]
#' }
#' 
#' # Build a listw object from the asymmetric matrix
#' listwAsym <-  mat2listw(matasym, style = "B", zero.policy = TRUE)
#' 
#' # Calculate Moran's I for AEM
#' MoranIAEM <- moran.randtest(AEM$vectors, listwAsym)
#' 
#' }
#' @keywords spatial
#' @export

`aem` <-
    function (aem.build.binary, binary.mat, weight, rm.link0 = FALSE, print.binary.mat=FALSE) 
    {
        if (missing(weight)) {
            weight <- 1
        }
        if(!is.numeric(weight)){
            stop("weight needs to be numeric")
        }
        
        if (missing(aem.build.binary)) {
            if (is.matrix(binary.mat) == FALSE) {
                stop("binary.mat is not a matrix")
            }
            res.mat <- as.matrix(binary.mat)
        }
        else {
            res.mat <- aem.build.binary[[1]]
            if (rm.link0) {
                link0 <- which(aem.build.binary[[2]][, 1] == 0)
                res.mat <- res.mat[, -link0]
            }
        }
        nsite <- nrow(res.mat)
        res.mat <- t(t(res.mat) * weight)
        
        res.mat.ct <- apply(res.mat, 2, scale, scale = FALSE)
        val.vec <- svd(res.mat.ct, nu = (nsite - 1), nv = 0)
        val <- val.vec$d[1:(nsite - 1)]^2/(nsite - 1)
        
        lim<-10^{-12}
        val.sel<-which(val>=lim)
        vec.lim<-val.vec$u[,val.sel]
        val.lim<-val[val.sel]
        
        if(print.binary.mat){
            res<-list(val.lim,vec.lim,res.mat)
            names(res)<-c("values","vectors","mod.binary.mat")
        }else{
            res <- list(val, vec.lim)
            names(res) <- c("values", "vectors")
        }
        return(res)
    }