From 574e4a9ffc7503feb10351a268c0523912f6257e Mon Sep 17 00:00:00 2001 From: paradis Date: Fri, 15 Jan 2010 12:41:31 +0000 Subject: [PATCH] add delta.plot + clarify doc in write.tree.Rd git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@104 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 5 +++- DESCRIPTION | 4 +-- R/delta.plot.R | 39 ++++++++++++++++++++++++++ man/delta.plot.Rd | 53 +++++++++++++++++++++++++++++++++++ man/mst.Rd | 5 ++-- man/write.tree.Rd | 2 +- src/delta_plot.c | 70 +++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 171 insertions(+), 7 deletions(-) create mode 100644 R/delta.plot.R create mode 100644 man/delta.plot.Rd create mode 100644 src/delta_plot.c diff --git a/ChangeLog b/ChangeLog index fa76411..0809491 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ - CHANGES IN APE VERSION 2.4-2 + CHANGES IN APE VERSION 2.5 NEW FEATURES @@ -6,6 +6,9 @@ NEW FEATURES o The new functions rTraitCont and rTraitDisc simulate continuous and discrete traits under a wide range of evolutionary models. + o The new function delta.plot does a delta plot following Holland et + al. (2002, Mol. Biol. Evol. 12:2051). + o add.scale.bar() has a new option 'ask' to draw interactively. diff --git a/DESCRIPTION b/DESCRIPTION index 3d42a84..1eec426 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape -Version: 2.4-2 -Date: 2010-01-11 +Version: 2.5 +Date: 2010-01-14 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/R/delta.plot.R b/R/delta.plot.R new file mode 100644 index 0000000..290c8fb --- /dev/null +++ b/R/delta.plot.R @@ -0,0 +1,39 @@ +## delta.plot.R (2010-01-12) + +## Delta Plots + +## Copyright 2010 Emmanuel Paradis + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +delta.plot <- function(X, k = 20, plot = TRUE, which = 1:2) +{ + if (is.matrix(X)) X <- as.dist(X) + n <- attr(X, "Size") + if (n < 4) stop("need at least 4 observations") + + ## add a category for the cases delta = 1 + ans <- .C("delta_plot", as.double(X), as.integer(n), + as.integer(k), integer(k + 1), double(n), + NAOK = TRUE, DUP = FALSE, PACKAGE = "ape") + counts <- ans[[4]] + ## add the counts of delta=1 to the last category: + counts[k] <- counts[k] + counts[k + 1] + counts <- counts[-(k + 1)] + + delta.bar <- ans[[5]]/choose(n - 1, 3) + + if (plot) { + if (length(which) == 2) layout(matrix(1:2, 1, 2)) + if (1 %in% which) { + barplot(counts, space = 0, xlab = expression(delta[q])) + a <- axTicks(1) + axis(1, at = a, labels = a/k) + } + if (2 %in% which) + plot(delta.bar, type = "h", ylab = expression(bar(delta))) + } + + invisible(list(counts = counts, delta.bar = delta.bar)) +} diff --git a/man/delta.plot.Rd b/man/delta.plot.Rd new file mode 100644 index 0000000..0a24cfe --- /dev/null +++ b/man/delta.plot.Rd @@ -0,0 +1,53 @@ +\name{delta.plot} +\alias{delta.plot} +\title{Delta Plots} +\usage{ +delta.plot(X, k = 20, plot = TRUE, which = 1:2) +} +\arguments{ + \item{X}{a distance matrix, may be an object of class ``dist''.} + \item{k}{an integer giving the number of intervals in the plot.} + \item{plot}{a logival specifying whether to draw the + \eqn{\delta}{delta} plot (the default).} + \item{which}{a numeric vector indicating which plots are done; 1: the + histogram of the \eqn{\delta_q}{delta_q} values, 2: the plot of the + individual \eqn{\bar{\delta}}{delta.bar} values. By default, both + plots are done.} +} +\description{ + This function makes a \eqn{\delta}{delta} plot following Holland et + al. (2002). +} +\details{ + See Holland et al. (2002) for details and interpretation. + + The computing time of this function is proportional to the fourth + power of the number of observations (\eqn{O(n^4)}), so calculations + may be very long with only a slight increase in sample size. +} +\value{ + This function returns invisibly a named list with two componente: + + \itemize{ + \item{counts}{the counts for the histogram of + \eqn{\delta_q}{delta_q} values} + \item{delta.bar}{the mean \eqn{\delta}{delta} value for each + observation} + } +} +\references{ + Holland, B. R., Huber, K. T., Dress, A. and Moulton, V. (2002) Delta + plots: a tool for analyzing phylogenetic distance data. + \emph{Molecular Biology and Evolution}, \bold{12}, 2051--2059. +} +\author{Emmanuel Paradis} +\seealso{ + \code{\link{dist.dna}} +} +\examples{ +data(woodmouse) +d <- dist.dna(woodmouse) +delta.plot(d) +delta.plot(d, 40, which = 1) +} +\keyword{hplot} diff --git a/man/mst.Rd b/man/mst.Rd index e2c0f60..06c9dfa 100644 --- a/man/mst.Rd +++ b/man/mst.Rd @@ -61,9 +61,8 @@ mst(X) \code{\link[stats]{dist}}, \code{\link[graphics]{plot}} } \examples{ -library(stats) -n <- 20 -X <- matrix(runif(n * 10), n, 10) +require(stats) +X <- matrix(runif(200), 20, 10) d <- dist(X) PC <- prcomp(X) M <- mst(d) diff --git a/man/write.tree.Rd b/man/write.tree.Rd index 70e1967..a90f4c5 100644 --- a/man/write.tree.Rd +++ b/man/write.tree.Rd @@ -6,7 +6,7 @@ write.tree(phy, file = "", append = FALSE, digits = 10, tree.names = FALSE) } \arguments{ - \item{phy}{an object of class \code{"phylo"}.} + \item{phy}{an object of class \code{"phylo"} or \code{"multiPhylo"}.} \item{file}{a file name specified by either a variable of mode character, or a double-quoted string; if \code{file = ""} (the default) then the tree is written on the standard output connection (i.e. the console).} diff --git a/src/delta_plot.c b/src/delta_plot.c new file mode 100644 index 0000000..66fadac --- /dev/null +++ b/src/delta_plot.c @@ -0,0 +1,70 @@ +/* delta_plot.c 2010-01-12 */ + +/* Copyright 2010 Emmanuel Paradis + +/* This file is part of the R-package `ape'. */ +/* See the file ../COPYING for licensing issues. */ + +#include + +void delta_plot(double *D, int *size, int *nbins, int *counts, double *deltabar) +{ + int x, y, u, v; /* same notation than in Holland et al. 2002 */ + int n = *size, nb = *nbins; + double dxy, dxu, dxv, dyu, dyv, duv, A, B, C, delta; + int xy, xu, xv, yu, yv, uv, i, k; + + for (x = 0; x < n - 3; x++) { + xy = x*n - x*(x + 1)/2; /* do NOT factorize */ + for (y = x + 1; y < n - 2; y++, xy++) { + yu = y*n - y*(y + 1)/2; /* do NOT factorize */ + dxy = D[xy]; + xu = xy + 1; + for (u = y + 1; u < n - 1; u++, xu++, yu++) { + uv = u*n - u*(u + 1)/2; /* do NOT factorize */ + dxu = D[xu]; + dyu = D[yu]; + xv = xu + 1; + yv = yu + 1; + for (v = u + 1; v < n; v++, xv++, yv++, uv++) { + dxv = D[xv]; + dyv = D[yv]; + duv = D[uv]; + /* if (dxy <= dxu && dxy <= dxv && dxy <= dyu && dxy <= dyv && dxy <= duv) k=1; */ +/* else if (duv <= dxy && duv <= dxu && duv <= dxv && duv <= dyu && duv <= dyv) k=1; */ +/* else if (dxu <= dxy && dxu <= dxv && dxu <= dyu && dxu <= dyv && dxu <= duv) k=2; */ +/* else if (dyv <= dxy && dyv <= dxu && dyv <= dxv && dyv <= dyu && dyv <= duv) k=2; */ +/* else if (dxv <= dxy && dxv <= dxu && dxv <= dyu && dxv <= dyv && dxv <= duv) k=3; */ +/* else if (dyu <= dxy && dyu <= dxu && dyu <= dxv && dyu <= dyv && dyu <= duv) k=3; */ + //Rprintf("%d\t%d\t%d\t%d\t%d\t%d\t\n", xy, xu, xv, yu, yv, uv); + //Rprintf("dxy = %f\tdxu = %f\tdyu = %f\n", dxy, dxu, dyu); + //Rprintf("D[xv] = %f\tD[yv] = %f\tD[uv] = %f\n", D[xv], D[yv], D[uv]); + /* switch (k) { */ +/* case 1 : A = dxv + dyu; B = dxu + dyv; C = dxy + duv; break; */ +/* case 2 : A = dxv + dyu; B = dxy + duv; C = dxu + dyv; break; */ +/* case 3 : A = dxu + dyv; B = dxy + duv; C = dxv + dyu; break; */ +/* } */ + //Rprintf("A = %f\tB = %f\tC = %f\n", A, B, C); + A = dxv + dyu; B = dxu + dyv; C = dxy + duv; + //Rprintf("A = %f\tB = %f\tC = %f\n", A, B, C); + if (A == B && B == C) delta = 0; else while (1) { + if (C <= B && B <= A) {delta = (A - B)/(A - C); break;} + if (B <= C && C <= A) {delta = (A - C)/(A - B); break;} + if (A <= C && C <= B) {delta = (B - C)/(B - A); break;} + if (C <= A && A <= B) {delta = (B - A)/(B - C); break;} + if (A <= B && B <= C) {delta = (C - B)/(C - A); break;} + if (B <= A && A <= C) {delta = (C - A)/(C - B); break;} + } + /* if (delta == 1) i = nb - 1; else */ + i = delta * nb; + //Rprintf("delta = %f\ti = %d\n", delta, i); + counts[i] += 1; + deltabar[x] += delta; + deltabar[y] += delta; + deltabar[u] += delta; + deltabar[v] += delta; + } + } + } + } +} -- 2.39.2