R/proximities.R

In thibautjombart/adephylo: Exploratory Analyses for the Phylogenetic Comparative Method

###########
# proxTips
###########


#' Compute some phylogenetic proximities between tips
#' 
#' The function \code{proxTips} computes a given proximity between a set of
#' tips of a phylogeny. A vector of tips is supplied: proximities between all
#' possible pairs of these tips are computed.  The proximities are computed
#' from the shortest path between the tips. \cr
#' 
#' Proximities are computed as the inverse (to the power \code{a}) of a
#' phylogenetic distance (computed by \code{\link{distTips}}. Denoting
#' \eqn{D=[d_{ij}]} a matrix of phylogenetic distances, the proximity matrix
#' \eqn{M=[m_{ij}]} is computed as: \deqn{m_{ij} = \frac{1}{d_{ij}^a} \forall i
#' \neq j}{ m_{ij} = (1/d_{ij})^a for all i different from j} and \deqn{m_{ii}
#' = 0}
#' 
#' Several distances can be used, defaulting to the sum of branch lengths (see
#' argument \code{method}).  Proximities are not true similarity measures,
#' since the proximity of a tip with itself is always set to zero.\cr
#' 
#' The obtained matrix of phylogenetic proximities (M) defines a bilinear
#' symmetric form when M is symmetric (default):\cr \deqn{f(x,y) = x^{T}My}
#' 
#' In general, M is not a metric because it is not positive-definite.  Such a
#' matrice can be used to measure phylogenetic autocorrelation (using Moran's
#' index): \deqn{I(x) = \frac{x^TMx}{var(x)}}{I(x) = (x^{T}Mx)/(var(x)) }
#' 
#' or to compute lag vectors (Mx) used in autoregressive models, like: \deqn{x
#' = Mx + ... + e} where '...' is the non-autoregressive part of the model, and
#' 'e' are residuals.
#' 
#' \code{Abouheif} proximity refers to the phylogenetic proximity underlying
#' the test of Abouheif (see references). Let P be the set of all the nodes in
#' the path going from \code{node1} to \code{node2}. Let DDP be the number of
#' direct descendants from each node in P. Then, the so-called 'Abouheif'
#' distance is the inverse of the product of all terms in DDP.
#' \code{oriAbouheif} returns a matrix with non-null diagonal elements, as
#' formulated in Pavoine \emph{et al.} (2008). This matrix is bistochastic (all
#' marginal sums equal 1), but this bilinear symmetric form does not give rise
#' to a Moran's index, since it requires a null diagonal. \code{Abouheif}
#' contains Abouheif's proximities but has a null diagonal, giving rise to a
#' Moran's index.\cr
#' 
#' \code{sumDD} refers to a phylogenetic proximity quite similar to that of
#' Abouheif. We consider the same sets P and DDP. But instead of taking the
#' inverse of the product of all terms in DDP, this proximity computes the
#' inverse of the sum of all terms in DDP. This matrix was denoted 'M' in
#' Pavoine \emph{et al.} (2008), who reported that it is related to May's index
#' (May, 1990).
#' 
#' @param x a tree of class \code{\link[ape:read.tree]{phylo}},
#' \code{\link[phylobase]{phylo4-class}} or \code{\link[phylobase]{phylo4d-class}}.
#' @param tips A vector of integers identifying tips by their numbers, or a
#' vector of characters identifying tips by their names. Distances will be
#' computed between all possible pairs of tips.
#' @param method a character string (full or abbreviated without ambiguity)
#' specifying the method used to compute proximities; possible values are:\cr -
#' \code{patristic}: (inversed sum of) branch length \cr - \code{nNodes}:
#' (inversed) number of nodes on the path between the nodes \cr -
#' \code{oriAbouheif}: original Abouheif's proximity, with diagonal (see
#' details) \cr - \code{Abouheif}: Abouheif's proximity without diagonal (see
#' details) \cr - \code{sumDD}: (inversed) sum of direct descendants of all
#' nodes on the path (see details) \cr
#' @param f a function to change a distance into a proximity.
#' @param normalize a character string specifying whether the matrix must be
#' normalized by row (\code{row}), column (\code{col}), or not (\code{none}).
#' Normalization amounts to dividing each row (or column) so that the marginal
#' sum is 1. Hence, default is matrix with each row summing to 1.
#' @param symmetric a logical stating whether M must be coerced to be symmetric
#' (TRUE, default) or not. This is achieved by taking (denoting N the matrix of
#' proximities before re-symmetrization): \deqn{M = 0.5 * (N + N^{T})} Note
#' that \eqn{x^{T}Ny = x^{T}My}, but the latter has the advantage of using a
#' bilinear symmetric form (more appropriate for optimization purposes).
#' @param useC a logical indicating whether computations of distances (before
#' transformation into proximities) should be performed using compiled C code
#' (TRUE, default), or using a pure R version (FALSE). C version is several
#' orders of magnitude faster, and R version is kept for backward
#' compatibility.
#' @return A matrix of phylogenetic proximities.
#' @author Thibaut Jombart \email{tjombart@@imperial.ac.uk}
#' @seealso \code{\link{distTips}} which computes several phylogenetic
#' distances between tips.
#' @references == About Moran's index with various proximities == \cr Pavoine,
#' S.; Ollier, S.; Pontier, D.; Chessel, D. (2008) Testing for phylogenetic
#' signal in life history variable: Abouheif's test revisited.
#' \emph{Theoretical Population Biology}: \bold{73}, 79-91.\cr
#' 
#' == About regression on phylogenetic lag vector == \cr Cheverud, J. M.; Dow,
#' M. M.; Leutenegger, W. (1985) The quantitative assessment of phylogentic
#' constaints in comparative analyses: sexual dimorphism in body weights among
#' primates. \emph{Evolution} \bold{39}, 1335-1351.\cr
#' 
#' Cheverud, J. M.; Dow, M. M. (1985) An autocorrelation analysis of genetic
#' variation due to lineal fission in social groups of Rhesus macaques.
#' \emph{American Journal of Phyisical Anthropology} \bold{67}, 113-121.\cr
#' 
#' == Abouheif's original paper ==\cr Abouheif, E. (1999) A method for testing
#' the assumption of phylogenetic independence in comparative data.
#' \emph{Evolutionary Ecology Research}, \bold{1}, 895-909.\cr
#' 
#' == May's index ==\cr May, R.M. (1990) Taxonomy as destiny. \emph{Nature}
#' \bold{347}, 129-130.
#' @keywords manip
#' @examples
#' 
#' if(require(ape) & require(phylobase)){
#' ## make a tree
#' x <- as(rtree(10),"phylo4")
#' plot(x, show.node=TRUE)
#' axisPhylo()
#' ## compute different distances
#' proxTips(x, 1:5)
#' proxTips(x, 1:5, "nNodes")
#' proxTips(x, 1:5, "Abouheif")
#' proxTips(x, , "sumDD")
#' 
#' ## see what one proximity looks like
#' M <- proxTips(x)
#' obj <- phylo4d(x,as.data.frame(M))
#' table.phylo4d(obj,symbol="sq")
#' }
#' 
#' @import phylobase
#' @export proxTips
proxTips <- function(x, tips="all",
                     method=c("patristic","nNodes","oriAbouheif","Abouheif","sumDD"),
                     f=function(x){1/x}, normalize=c("row","col","none"), symmetric=TRUE, useC=TRUE){

    ## if(!require(phylobase)) stop("phylobase package is not installed")

    ## handle arguments
    x <- as(x, "phylo4")
    method <- match.arg(method)
    normalize <- match.arg(normalize)
    N <- nTips(x)
    if(tips[1]=="all") { tips <- 1:N }
    tips <- getNode(x, tips)

    ## some checks
    if (is.character(checkval <- checkPhylo4(x))) stop(checkval)
    if(any(is.na(tips))) stop("wrong tips specified")

    ## compute distances
    distMethod <- method
    if(length(grep("Abouheif", distMethod)>1)){
        distMethod <- "Abouheif"
    }
    D <- distTips(x, tips=tips, method=distMethod, useC=useC)
    D <- as.matrix(D)

    ## compute proximities
    res <- f(D)
    diag(res) <- 0

    ## handle Abouheif with diagonal (Abouheif1)
    if(method=="oriAbouheif"){
        sumMarg <- apply(res,1,sum)
        diag(res) <- (1-sumMarg)
        normalize <- "none" # not needed (already bistochastic)
        symmetric <- FALSE # not needed (aleady symmetric)
    }

    ## standardization
    if(normalize=="row") {
        res <- prop.table(res, 1)
    }

    if(normalize=="col") {
        res <- prop.table(res, 2)
    }

    ## re-symmetrize
    if(symmetric){
        res <- 0.5 * (res + t(res))
    }

    ## set the output
    return(res)

} # end proxTips