- CHANGES IN APE VERSION 2.4-2
+ CHANGES IN APE VERSION 2.5
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.
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 <Emmanuel.Paradis@ird.fr>
--- /dev/null
+## 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))
+}
--- /dev/null
+\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}
\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)
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).}
--- /dev/null
+/* 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 <R.h>
+
+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;
+ }
+ }
+ }
+ }
+}
projects / ape.git / commitdiff