R/PEM-utils.R

In guenardg/MPSEM: Modelling Phylogenetic Signals using Eigenvector Maps

## **************************************************************************
##
##    (c) 2010-2025 Guillaume Guénard
##        Department de sciences biologiques,
##        Université de Montréal
##        Montreal, QC, Canada
##
##    ** Phylogenetic Eigenvector Maps (PEM) - Utility Functions **
##
##    This file is part of MPSEM
##
##    MPSEM is free software: you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation, either version 3 of the License, or
##    (at your option) any later version.
##
##    MPSEM is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
##    GNU General Public License for more details.
##
##    You should have received a copy of the GNU General Public License
##    along with MPSEM. If not, see <https://www.gnu.org/licenses/>.
##
##    R source code file
##
## **************************************************************************
##
#' Utility Function for Phylogenetic Eigenvector Maps
#' 
#' @description A set of functions to help in calculating Phylogenetic
#' Eigenvector Maps (PEM).
#' 
#' @name PEM-utils
#' 
#' @param x A \code{\link{graph-class}} object.
#' @param e A numeric vector of indices or a character vector of names
#' representing one or more edges.
#' @param v A numeric vector of indices or a character vector of names
#' representing one or more vertices. The default for function
#' \code{aggregateOnVertex} is all the vertices of \code{x}.
#' @param data A \code{\link{matrix}} or \code{\link{data.frame}} containing
#' edge data to be aggregated on the graph's vertices.
#' @param fun A function to be applied on the incoming edges of the vertices.
#' @param ... Supplementary arguments to be internally passed to \code{fun}.
#' @param default A vector of values to be assigned to vertices that have no
#' incoming edges. The default is a vector of zeros with the same length as the
#' number of columns is \code{data}.
#' 
#' @details Functions \code{isUnderEdge} and \code{isUnderVertex} are used to
#' set individual graph locations, in terms of edges or vertices where shifts in
#' the values of selection (a) and/or evolution rate (psi) parameters are
#' occurring. The binary contrasts matrices generated by function
#' \code{isUnderEdge} and \code{isUnderVertex} are used as descriptors in models
#' representing these subordinate parameters using higher order ones.
#' 
#' Function \code{aggregateOnVertex} allows one to aggregate variables
#' associated with the edges on their destination vertices using a user-defined
#' function. 
#' 
#' @return Functions \code{isUnderEdge} and \code{isUnderVertex} return a matrix
#' of binary contrasts whose rows represent the edges of the graph and columns
#' the edges or vertices given as indices or names to function
#' \code{isUnderEdge} (argument \code{e}) or \code{isUnderVertex} (argument
#' \code{v}), respectively. Function \code{aggregateOnVertex} returns a matrix
#' whose rows represent the vertices specified by its argument \code{v} and
#' columns represent the variables given using its argument \code{data}.
#' 
#' @author \packageAuthor{MPSEM} --
#' Maintainer: \packageMaintainer{MPSEM}
#' 
#' @references
#' Guénard, G., Legendre, P., and Peres-Neto, P. 2013. Phylogenetic eigenvector
#' maps: a framework to model and predict species traits. Methods in Ecology 
#' and Evolution. 4: 1120--1131
#' 
#' Makarenkov, V., Legendre, L. & Desdevise, Y. 2004. Modelling phylogenetic
#' relationships using reticulated networks. Zoologica Scripta 33: 89--96
#' 
#' Blanchet, F. G., Legendre, P. & Borcard, D. 2008. Modelling directional
#' spatial processes in ecological data. Ecological Modelling 215: 325--336
#' 
#' @seealso \code{\link{PEM-functions}}, \code{\link{PEM-class}}
#' 
#' @examples
#' ## Synthetic example
#' data.frame(
#'   species = as.logical(c(0,0,1,0,0,0,0,0,0,0,1,1,1)),
#'   type = c(2,2,3,1,2,2,2,2,2,2,3,3,3),
#'   x = c(1,3,4,0,1.67,4,1,1.33,2.33,3.33,4.33,4,5),
#'   y = c(1,1,1,0,0.5,0,-1,0,0,-0.5,-1,-0.5,-0.5),
#'   row.names = sprintf("V%d",1:13)
#' ) %>%
#' st_as_sf(
#'   coords=c("x","y"),
#'   crs = NA
#' ) %>%
#' graph %>%
#' add.edge(
#'   from = c(1,2,1,5,4,4,5,9,4,8,9,4,7,7,6,6,9,10,10),
#'   to = c(2,3,5,2,1,5,9,2,8,9,6,7,8,9,3,13,10,12,11),
#'   data = data.frame(
#'     distance = c(4.2,4.7,3.9,3.0,3.6,2.7,4.4,3.4,3.6,3.3,
#'                  4.8,3.2,3.5,4.4,2.5,3.4,4.3,3.1,2.2),
#'     row.names = sprintf("E%d",1:19)
#'   )
#' ) -> gr_ex
#' 
#' ## Plot the graph:
#' plot(gr_ex, cex.min=3, cex.lab=0.6)
#' 
#' ## Show the edges of the graph:
#' edge(gr_ex)
#' 
#' ## Identify the edges that are under or at the edges E7 or E17 using a binary
#' ## contrast matrix:
#' isUnderEdge(gr_ex, c("E7","E17"))
#' 
#' ## Identify the edges that are under vertices V5 or V9 using a binary
#' ## contrast matrix:
#' tmp <- isUnderVertex(gr_ex, c("V5","V9"))
#' tmp
#' 
#' ## Aggregate the result of isUnderVertex() using the following function
#' ## enables one to determine which of the vertices are found under V5 and V9:
#' aggregateOnVertex(gr_ex, tmp, function(x) ifelse(any(as.logical(x)),1,0))
#' 
#' 
NULL
#' 
#' @describeIn PEM-utils
#' 
#' Edges Under or At an Edge
#' 
#' Identify which of the edges of a graph are located at or under one or more
#' of the graph's edges using a binary contrast matrix.
#' 
#' @export
isUnderEdge <- function(x, e) {
  
  ifl <- InflMat(x)
  
  if(is.numeric(e)) {
    nms <- edgenames(ifl)[e]
  } else if(is.character(e)) {
    nms <- e
    e <- match(e, edgenames(ifl))
  } else
    stop("Argument 'e' must be either numeric or character.")
  
  matrix(
    0,
    nedge(ifl),
    length(e),
    dimnames = list(
      edgenames(ifl),
      nms
    )
  ) -> out
  
  for(i in 1L:length(e))
    out[edge(ifl)[[2L]] %in% which(as.logical(ifl[,e[i]])),i] <- 1
  
  out
}
#' 
#' @describeIn PEM-utils
#' 
#' Edges Under a Vertex
#' 
#' Identify which of the edges of a graph are located under one or more of the
#' graph's vertices using a binary contrast matrix.
#' 
#' @export
isUnderVertex <- function(x, v) {
  
  ifl <- InflMat(x)
  
  if(is.numeric(v)) {
    nms <- rownames(x)[v]
  } else if(is.character(v)) {
    nms <- v
    v <- match(v, rownames(x))
  } else
    stop("Argument 'e' must be either numeric or character.")
  
  matrix(
    0,
    nedge(ifl),
    length(v),
    dimnames = list(
      edgenames(ifl),
      nms
    )
  ) -> out
  
  for(i in 1L:length(v)) {
    
    de <- which(edge(ifl)[[1L]] == v[i])
    
    if(length(de) > 1L) {
      tmp <- which(apply(ifl[,de], 1L, any))
    } else
      tmp <- which(as.logical(ifl[,de]))
    
    out[edge(ifl)[[2L]] %in% tmp,i] <- 1
  }
  
  out
}
#' 
#' @describeIn PEM-utils
#' 
#' Aggregate on Vertex
#' 
#' Aggregate incoming edge values on a vertex or a set of vertices.
#' 
#' @export
aggregateOnVertex <- function(x, data, fun, ..., v = 1L:NROW(x),
                              default = rep(0, NCOL(data))) {
  
  if(!inherits(x,"graph"))
    stop("Argument 'x' must be a graph-class object.")
  
  if(!is.matrix(data))
    data <- as.matrix(data)
  
  if(nrow(data) != nedge(x))
    stop("Graph 'x' has ", nedge(x), " edges, but table 'data' has ",
         nrow(data), " rows")
  
  if(ncol(data) != length(default))
    stop("Argument 'data' is a ", ncol(data), "-column table, but argument ",
         "'default' is length ",length(default))
  
  if(is.character(v)) {
    
    if(any(is.na(tmp <- match(v, rownames(x)))))
      stop("Unknown vertex: ", paste(v[is.na(tmp)], collapse=", "))
    
    v <- tmp
  }
  
  if(length(tmp <- which((v < 0) | (v > nrow(x)))))
    stop("Vertex out of range: ", paste(v[tmp], collapse=", "))
  
  matrix(NA, length(v), ncol(data),
         dimnames = list(rownames(x)[v],colnames(data))) -> out
  
  for(i in 1L:length(v)) {
    
    wh <- which(edge(x)[[2L]] == i)
    
    if(!length(wh)) {
      out[i,] <- default
    } else if(length(wh) > 1L) {
      out[i,] <- apply(data[wh,,drop=FALSE], 2L, fun, ...)
    } else
      out[i,] <- data[wh,]
  }
  
  out
}
#'